Negative Intrinsic Viscosity in Graphene Nanoparticle Suspensions Induced by Hydrodynamic Slip

TL;DR

Graphene nanoparticle suspensions can have lower viscosity than pure water due to hydrodynamic slip, which could help design better lubricants and inks.

Contribution

Molecular simulations confirm negative intrinsic viscosity in graphene suspensions due to hydrodynamic slip at the liquid–solid interface.

Findings

Intrinsic viscosity becomes negative for graphene sheets with aspect ratios above ≈5.5.

Hydrodynamic slip suppresses particle rotation and reduces viscous dissipation.

Viscosity initially decreases with concentration in the dilute regime before rising due to aggregation.

Abstract

The viscosity of nanoparticle suspensions is always expected
to
increase with particle concentration. However, a growing body of experiments
on suspensions of atomically thin nanomaterials such as graphene contradicts
this expectation. Some experiments indicate effective suspension viscosities
below that of pure solvent at high shear rates and low solid concentrations,
i.e., the intrinsic viscosity is negative. Using molecular dynamics
simulations, we investigate the shear viscosity of few-nanometer graphene
sheets in water at high Péclet numbers (Pe ≥ 100), for
aspect ratios from 4.5 to 12.0. These simulations robustly confirm
that the intrinsic viscosity decreases with increasing aspect ratio
and becomes negative beyond a threshold ≈5.5, providing a molecular-level
confirmation of this behavior in a realistic graphene–water
system. Comparison with continuum boundary integral modeling shows
quantitative agreement in the dilute regime, confirming the effect
is hydrodynamic in origin. We demonstrate that this anomalous behavior
originates from hydrodynamic slip at the liquid–solid interface,
which suppresses particle rotation and promotes stable alignment with
the flow direction, thereby reducing viscous dissipation relative
to dissipation in pure solvent. This slip mechanism holds for both
fully 3D disc-like and quasi-2D particle geometries explored in the
molecular simulations. As the concentration of graphene particles
increases in the dilute regime, the viscosity initially decreases,
falling below that of pure water. At higher concentrations, however,
particle aggregation becomes significant, leading to a rise in viscosity
after a minimum is reached. Our work has important implications for
the design of lubricants, inks, and nanocomposites with tunable viscosity.

Figures

Molecular dynamics simulation setup. The two panels illustrate the two nanoparticle geometries used: (a) several quasi-2D graphene sheets; (b) fully-3D disc-like graphene sheet. Graphene sheets are shown in pink. Selected water molecules near the graphene particle are rendered explicitly (blue and white spheres for O and H, respectively), with the remaining fluid shown as a translucent isosurface. The confining FeO walls (shown in orange/gray) translate along x̂ to apply a simple shear flow. The simulations are periodic along x̂ and ẑ. (c) Schematic definition of the particle orientation angles θ and φ, where the unit director vector corresponding to the particle’s symmetry axes is given by d̂ = (−sinθsinφ, sinθcosφ, cosθ).

Single-particle simulation results for slip (blue) and no-slip (red) graphene particles. (a) α vs a/b for fully-3D (open symbols) and quasi-2D (filled symbols) particles. Triangles indicate data from ref . Lines are 2D BI predictions. Shaded blue represents α < 0. (b) Conditional intrinsic viscosity α̃(φ) for quasi-2D particles with a/b = 10. The error bars represent the standard error. (c) Logarithm of the orientation probability distribution log p(φ) for the same systems. To illustrate how α̃(φ) reflects the hydrodynamic resistance generated by a particle, in (c) we include insets showing MD snapshots of particle configurations and the surrounding fluid velocity fields.

Orientation dynamics and hydrodynamic torque for a fully-3D disc-like particle. (a, b) Time evolution of (a) azimuthal angle φ and (b) polar angle θ for slip (blue) and no-slip (red) particles, for two initial orientations. (c, d) x-component of the hydrodynamic torque for a disc held at fixed orientation, plotted as a function of: (c) φ, for two values of θ; and (d) θ, for a fixed value of φ.

(a) Relative shear viscosity vs solid fraction for fully-3D and quasi-2D particles with a/b = 10. The crosses represent BI simulation results. (b) Reduced shear viscosity values corresponding to the data in (a). Dashed lines are least-squares fits to the data.

Tables

1. Coefficients α, β and minimum viscosity concentration c min (where applicable) for graphene particles with a/b = 10

particle type	α	β	c min = – α/(2β)
quasi-2D No-Slip	3.4	41.1	N/A
quasi-2D Slip	–1.4	5.4	0.13
fully-3D Slip	–1.2	21.1	0.03

Full text

Introduction

Two-dimensional (2D) materials such as graphene have attracted significant scientific interest over the past two decades due to their exceptional electronic, thermal, mechanical, and optical properties, resulting in new research directions and enabling the development of technologies that are now entering commercial applications. ?−? ? ? ? ? ? Many applications, ranging from lubricants and nanocomposites to conductive inks and coatings, involve the dispersion of 2D material nanoparticles in liquids. Predicting the flow behavior of such suspensions requires quantifying their rheological properties. ?,? Among these, the effective steady shear viscosity is the most important rheological property.

Recent experiments on suspensions of sheet-like materials of nanometric thickness in both simple and complex fluids have revealed rheological behaviors that cannot be explained by classical suspension hydrodynamics.? There is increasing evidence that the addition of graphene nanosheets to polymer nanocomposites can reduce the shear viscosity. ?−? ? Relative viscosities (defined as the ratio of effective viscosity to suspending fluid viscosity) well below unity have been reported for suspensions of graphene nanoparticles in lubricating oils in the regime of low particle concentrations. ?,? A similar trend has been observed for few-layer graphene dispersed in the ionic liquid [HMIM]­BF_4_. ?,? Plate-like nanoparticles of α-zirconium phosphate and yttrium oxide, with geometric aspect ratios around 10, have also been reported to reduce the shear viscosity of mineral oil in the dilute regime, with viscosity decreasing as the particle concentration increases. ?,?

These findings are unexpected. Classical suspension rheology predicts that shear viscosity always increases with solid concentration, irrespective of particle shape or Brownian effects.? Physically, this increase is due to the fact that suspended particles disturb the flow because they impose a boundary condition on the velocity field. The flow gradients associated with this disturbance translate to an enhanced viscous dissipation compared to the pure solvent. This classical result has been confirmed by extensive theoretical, experimental, and simulation studies since at least the pioneering work of Einstein, who calculated an expression for the viscosity of dilute suspensions of rigid spheres.? In contrast, the molecular dynamics (MD) simulations of this paper unambigously demonstrate that graphene nanoparticles, due to their combination of pronounced hydrodynamic slip and atomic-scale thickness, display a suspension viscosity lower than that of the pure fluid. To our knowledge, this phenomenon is unique to ultrathin, large-slip particles and has not been observed in conventional colloidal systems.

The anomalous viscosity reduction we report is not only of fundamental interest but also holds promise for technological applications. The development of rheological modifiers that can reduce viscosity, or at least mitigate its increase, could revolutionize the field of lubrication. Approximately 23% of the global energy consumption (over 100 EJ) is attributed to tribological contacts, with 20% of this energy lost to friction.? Even modest reductions in the viscosity of nanoparticle-based lubricants could, therefore, yield significant energy savings.

The effective shear viscosity η of a dilute suspension of small, rigid particles can be expanded in powers of the particle volume fraction c as

$$[eqn]$$

where η_0_ is the viscosity of the pure solvent, α is the intrinsic viscosity, and β quantifies interparticle interactions.? For c ≪ 1, the leading-order term dominates and η ≈ η_0_(1 + αc). A negative intrinsic viscosity (α < 0) thus implies that adding particles reduces the viscosity in the dilute regime. Recent boundary integral simulations on model plate-shaped particles in unbounded two-dimensional shear flow have shown that α can become negative, provided that (i) the particles are sufficiently thin and (ii) they exhibit hydrodynamic slip with a slip length exceeding a threshold value comparable to the particle thickness. ?,? While these continuum simulations offer valuable insights, they are based on idealized particle geometries and do not account for molecular-level effects such as water structuring at the liquid–solid interface or adhesive interactions between particles. Moreover, the continuum approach is only an approximation for graphene, whose thickness is comparable to the size of solvent molecules, challenging the validity of continuum hydrodynamics at this scale. Thus, it remains unclear whether negative intrinsic viscosity can arise in realistic molecular models of graphene. For example, previous molecular dynamics simulations of hexabenzocoronene, a disc-like nanographene with a geometric aspect ratio ≈3.2, did not display a negative intrinsic viscosity.?

In the current work, we systematically investigate the conditions required for observing a negative intrinsic viscosity. We employ MD simulations to compute the shear viscosity of suspensions of single-layer graphene nanoparticles with varying aspect ratios in water at high Péclet numbers. Nonfunctionalized, nonoxidized graphene is chosen as a model 2D material due to its atomic-scale thickness (≈0.5 nm for a monolayer in water)? and its characteristically large hydrodynamic slip length, which typically ranges from 10 to 100 nm in water, small-chain alcohols, and certain ionic liquids. ?−? ? A substantial slip length at the solid–liquid interface is a precondition for the occurrence of negative intrinsic viscosity, as demonstrated in recent theoretical studies. ?,? Water is selected as the solvent for our MD simulations owing to the availability of reliable force fields and the presence of high-quality MD data for carbon-based materials in aqueous environments. ?−? ? ? Beyond quantifying the intrinsic viscosity, we also investigate the effect of interparticle interactions on the second-order coefficient β in eq. Our analysis reveals that β is highly sensitive to particle geometry and highlights significant discrepancies between MD and continuum predictions, particularly at higher concentrations where molecular-scale effects and particle aggregation become important.

Results and Discussion

MD simulations of plate-like nanoparticles in water were carried out using LAMMPS? (Figure). The particles were suspended in a shear flow ** u ** ^∞^ = γ̇y x̂ produced by two translating walls. Here, y denotes the coordinate along the gradient direction (aligned with ŷ) relative to the center of the simulation domain, the flow is directed along x̂, and γ̇ is the applied shear rate. Periodic boundary conditions were applied along the x̂ and ẑ directions. To generate the shear flow, a constant translational velocity was imposed on two solid walls. The walls, composed of iron­(II) oxide (FeO), were characterized by a high friction coefficient, ensuring a negligible slip length (λ ≈ 0 nm) in water. The nanoparticles were modeled as pristine graphene monolayers and described with the OPLS all-atom force field, which includes explicit bond, angle, and dihedral interactions. ?,? Previous molecular dynamics simulations estimated the bending rigidity of graphene in water to be on the order of 1 eV. ?,? Based on our selected force fields, the hydrodynamic slip length at the graphene–water interface is estimated to be λ ≈ 60 nm, consistent with previous molecular dynamics studies. ?,? For comparison, no-slip nanoparticles were simulated by artificially increasing the graphene–water interaction energy, as described in the Methods section. This adjustment effectively eliminated hydrodynamic slip at the particle surface, providing a contrasting behavior to that of pristine graphene. The system was maintained at a constant temperature of 300 K using a Nosé-Hoover thermostat and the fluid was maintained at a pressure of 1 atm by applying a vertical force on the top wall. Further simulation details, including equilibration and production run protocols, are provided in the Methods section.

Molecular dynamics simulation setup. The two panels illustrate the two nanoparticle geometries used: (a) several quasi-2D graphene sheets; (b) fully-3D disc-like graphene sheet. Graphene sheets are shown in pink. Selected water molecules near the graphene particle are rendered explicitly (blue and white spheres for O and H, respectively), with the remaining fluid shown as a translucent isosurface. The confining FeO walls (shown in orange/gray) translate along x̂ to apply a simple shear flow. The simulations are periodic along x̂ and ẑ. (c) Schematic definition of the particle orientation angles θ and φ, where the unit director vector corresponding to the particle’s symmetry axes is given by d̂ = (−sinθsinφ, sinθcosφ, cosθ).

With this simulation setup, we investigated two nanoparticle geometries: (i) fully-3D disc-like graphene nanosheets (Figureb), which are free to rotate and translate in all spatial directions, and (ii) quasi-2D rectangular nanosheets that extend across the entire simulation box along the ẑ direction (Figurea), effectively confining their motion and rotation to the x̂ – ŷ plane due to periodic boundary conditions. In both geometries, the graphene sheets retain flexibility and can undergo deformation in response to the applied shear flow.

The effective suspension viscosity, η, was determined by measuring the shear stress on the confining walls. Simulations were first performed with a single particle, corresponding to a dilute solid fraction of approximately c ≃ 0.01. At this low concentration, the reduced viscosity, (η/η_0_ – 1)/c provides a good approximation to the intrinsic viscosity α. We investigated a range of particle aspect ratios (4.5 < a/b < 12) for a fixed shear rate of γ̇ = 5 × 10^10^ s^–1^, where 2a and 2b denote the effective length (or disc diameter) and thickness of the particle, respectively. We defined these two dimensions taking into account the effective size of the constituent atoms.? To maintain a constant c while varying particle length, the channel height H was adjusted while keeping the x–y box dimensions fixed, and using H/a ≈ 2.2 (this value should be sufficient to ensure that wall-particle hydrodynamic interactions are comparatively weak). ?,? The viscosity η_0_ was obtained from solvent-only simulations carried out in a computation domain of identical geometry as the particle-laden simulations. The corresponding Péclet number, Pe = γ̇/D r, ranged from 80 to 2800, with the lowest Péclet number corresponding to the shortest particle. The rotational diffusion coefficient, D r, was estimated for quasi-2D graphene particles using the expression D r ≈ k B T/(2πηa ^3^) developed for an infinitely thin 2D rigid plate.? For 3D particles we used D r = 3k B T/(32ηa ^3^).? Here, k B is the Boltzmann constant, T is the temperature, and a is the particle’s half-length. For $[eqn]$ , the rotational dynamics is controlled by hydrodynamics, and the influence of Brownian noise on the orientational distribution function is minimal.?

Figurea shows the intrinsic viscosity as a function of the particle aspect ratio for both slip and no-slip particles. For slip particles, α decreases monotonically with increasing aspect ratio, becoming negative at a/b ≃ 5.5. In contrast, no-slip particles consistently exhibit positive values of α, with a much weaker dependence on a/b compared to slip particles. Across the range of a/b ∈ [4, 10] considered, α remains approximately constant at α ≈ 4, with a minor reduction observed at the highest aspect ratio examined (a/b = 10). This reduction coincides with enhanced particle alignment in the flow direction and less frequent tumbling at higher aspect ratios, as evidenced by the molecular dynamics simulations presented in the Supporting Information.

Single-particle simulation results for slip (blue) and no-slip (red) graphene particles. (a) α vs a/b for fully-3D (open symbols) and quasi-2D (filled symbols) particles. Triangles indicate data from ref . Lines are 2D BI predictions. Shaded blue represents α < 0. (b) Conditional intrinsic viscosity α̃(φ) for quasi-2D particles with a/b = 10. The error bars represent the standard error. (c) Logarithm of the orientation probability distribution log p(φ) for the same systems. To illustrate how α̃(φ) reflects the hydrodynamic resistance generated by a particle, in (c) we include insets showing MD snapshots of particle configurations and the surrounding fluid velocity fields.

Overall, Figurea underscores the dramatic reduction in α due to hydrodynamic slip across all aspect ratios, with the effect becoming more pronounced at larger a/b. Notably, while previous boundary element simulations by Allison? for ellipsoidal particles at low Péclet numbers also reported a significant reduction in α due to slip, the intrinsic viscosity remained positive even for ‘perfect slip’ ellipsoids with λ → ∞. For example, for a/b = 10, Allison observed a decrease in α from 7.9 to 1.7 for the no-slip and perfect slip cases, respectively. In contrast, our results for high Péclet numbers reveal a change in sign of the intrinsic viscosity, with α dropping from 3.8 ± 0.2 to −1.6 ± 0.2 for quasi-2D particles and from 5.4 ± 0.3 to −1.0 ± 0.2 for 3D disks.

To further investigate the mechanism responsible for the negative intrinsic viscosity, we examine the orientation statistics of the particles and the intrinsic viscosity conditioned on the particle orientation. In simple shear flow, the xy-component of the time-averaged macroscopic shear stress, ⟨Σ_ xy _⟩, comprises two contributions: the solvent stress, η_0_γ̇, and the particle-induced stress, ⟨Σ_p,xy _⟩.? In the dilute limit, the particle contribution can be expressed as ⟨Σ_p,xy _⟩ = n⟨S _ xy _⟩, where ⟨S _ xy _⟩ is the xy-component of the time-averaged stresslet tensor for an isolated particle, and n is the particle number density.? At large Pe, the stresslet is determined primarily by hydrodynamic stresses and represents the deviatoric component of the hydrodynamic traction on the particle’s surface.? In Stokes flow, the instantaneous traction at time t depends solely on the particle’s configuration at that time. For simplicity, we restrict our attention to a quasi-2D particle whose surface normal (d̂) is constrained to lie within the flow-gradient (x̂–ŷ) plane. In this configuration, the particle’s orientation is fully characterized by the angle φ, defined as the counterclockwise angle between d̂ and ŷ. The time-averaged shear stresslet is expressed as an average of the stresslet S _ xy _(φ), conditioned on the particle adopting orientation φ and weighted by the corresponding probability distribution p(φ):

$$[eqn]$$

The orientational space φ ∈ [−π/2, π/2] represents the range of possible angles explored by the particle. Similarly, the intrinsic viscosity, α, is obtained by averaging the orientation-dependent intrinsic viscosity, α̃(φ):

$$[eqn]$$

This expression highlights the critical role of particle orientation in determining the shear viscosity of the suspension.

Using MD simulations, we calculated the conditional intrinsic viscosity α̃(φ) and the orientation probability distribution p(φ) for quasi-2D particles. The results, shown in Figureb,c, correspond to particles with a/b = 10, for γ̇ = 5 × 10^10^ s^–1^. The conditional intrinsic viscosity was determined by measuring the wall shear stress at instants when the particle adopted a given orientation angle φ, while the orientational distribution p(φ) was computed by histogramming the frequency of φ over the simulation trajectory.

No-slip particles exhibit an approximately cosinusoidal variation in α̃(φ), with a minimum at φ = 0 (Figureb). Unlike no-slip particles, which sample the full range of orientations due to continuous tumbling, slip particles perform fluctuations near a fixed angle, resulting in a peaked profile of p(φ) (Figurec). The value of φ about which the fluctuations occur is known to be governed by the particle aspect ratio and by the hydrodynamic slip length.? Because α̃(φ) < 0 within the narrow range for which p is nonzero, the time-averaged intrinsic viscosity α is also negative.

The insets of Figurec, which show representative snapshots from the MD simulations, describe typical particle configurations and the corresponding fluid velocity fields. These images clearly demonstrate that slip particles induce significantly less disturbance in the surrounding flow compared to no-slip particles. The deformation of slip graphene particles is small in comparison to that of no-slip graphene. This observation is explained by the higher shear stress for buckling of slip sheets with respect to no-slip sheets. ?,?

We compared our MD simulations with non-Brownian, single-particle boundary integral (BI) calculations for quasi-2D rigid particles. The details of BI calculations are given in the Methods section. Note that the use of non-Brownian BI calculations is justified by the weak dependence of α on Pe for Pe ≳ 100 regardless of the slip length; the effect of Pe is discussed further in the Supporting Information.

The MD results are in excellent agreement with BI predictions in the slip case (λ = 60 nm) across all aspect ratios (full lines in Figurea). For no-slip particles (λ = 0 nm), MD simulations predict larger values of α than the BI calculations. This difference could be explained by the large interaction energy between carbon and oxygen atoms used to enforce the no-slip boundary condition in our MD calculations. A high solid–liquid interaction energy results in water molecules being strongly adsorbed to the graphene platelet, effectively ‘sticking’ to it and leading to an increased effective thickness. This can be observed from the radial distribution function of oxygen atoms around carbon atoms presented in the Supporting Information. It is known that strongly structured fluid molecules around a particle result in higher viscosity.?

We also computed α̃(φ) using BI simulations by evaluating the instantaneous single-particle stresslet, given by α̃(φ) = S _ xy _(φ)/(η_0_γ̇A p), where A p is the particle’s cross-sectional area (solid lines in Figureb). The close quantitative agreement between the MD and BI predictions shows that both atomistic and continuum simulations describe the same orientation-dependent contributions to the viscosity, reinforcing the robustness of our approach.

Our MD simulations reveal minimal difference in shear viscosity between fully-3D and quasi-2D particle geometries. This agreement arises because of the propensity of a fully-3D slip graphene to align on average with its normal in the flow plane. We describe the instantaneous orientation of each fully-3D graphene particle using the unit vector d̂, which is aligned with the particle’s axis of symmetry. For brevity, we refer to the particle’s orientation simply as the orientation of d̂. To compute d̂, we apply principal component analysis (PCA) to the centered atomic coordinates of the disc and define d̂ as the eigenvector corresponding to the smallest eigenvalue of the covariance matrix, following Hoppe et al.? To uniquely define the director, we choose the sign of d̂ such that d̂ ·ŷ > 0. In a spherical coordinate system, the components of d̂ are

$$[eqn]$$

where θ is the polar angle with respect to the vorticity direction ẑ, and φ is the azimuthal angle in the flow-gradient (x̂ – ŷ) plane (see Figurec). Our MD simulations show that, regardless of their initial orientation, slip discs rapidly (within the first few hundred picoseconds) relax to an orientation for which θ ≃ π/2 and φ ≪ 1 (Figurea,b). This stable alignment of particle with the flow-gradient plane is consistent with previous MD simulations of nanographene reported by Gravelle et al.? The observation that a single 3D graphene disk eventually aligns such that its symmetry axis lies in the flow-gradient plane, performing a motion that is statistically two-dimensional, helps understanding why the quasi-2D particle simulations (Figurea) and the fully-3D disk simulations (Figureb) give quantitatively comparable results.

Orientation dynamics and hydrodynamic torque for a fully-3D disc-like particle. (a, b) Time evolution of (a) azimuthal angle φ and (b) polar angle θ for slip (blue) and no-slip (red) particles, for two initial orientations. (c, d) x-component of the hydrodynamic torque for a disc held at fixed orientation, plotted as a function of: (c) φ, for two values of θ; and (d) θ, for a fixed value of φ.

In contrast, the no-slip discs exhibit persistent tumbling and periodic oscillations in both φ and θ, consistent with Jeffery’s classical theory (Figurea,b). Occasionally, these particles transition between Jeffery orbits, gradually shifting from tumbling toward orbits closer to θ = ±π/2. This transition is facilitated by weak Brownian rotational diffusion, which enables the particle to explore a range of Jeffery orbits. Such behavior is well-described by the theoretical framework of Leal and Hinch,? which extends Jeffery’s theory to account for the influence of weak Brownian noise.

To further understand the tendency of graphene particle to reorient toward a stable orientation, we computed the hydrodynamic torque component along the flow direction, τ_X_, exerted by the fluid on a particle held at fixed orientation (shear rate γ̇ = 5 × 10^10^ s^–1^). Figurec compares τ_X_ as a function of the azimuthal angle φ for two polar angles. For θ = π/2, the torque remains nearly zero for the entire range of φ, consistent with a stable in-plane orientation. In contrast, for a more tilted configuration (θ = 3π/4), the torque exhibits a pronounced dependence on φ. Specifically, τ_X_ is negative for φ > 0, indicating a restoring torque toward the flow-gradient plane. While τ_X_ is positive for φ < 0, indicating a torque that drives the particle further from the flow-gradient plane. Figured illustrates the dependence of τ_X_ on θ for a fixed azimuthal angle φ = π/12, which corresponds to the steady-state orientation typically observed for slip discs in our simulations (see Figurea). The torque is positive for θ < π/2, negative for θ > π/2, and vanishes at θ = π/2, indicating a restoring hydrodynamic moment.

The origin of this torque can be understood by considering the distribution of hydrodynamic stresses on the particle surface. For slip particles, tangential stresses on the broad faces are strongly suppressed, so the dominant contribution to the torque comes from the edges.? The restoring torque appears only when the particle is simultaneously tilted out of the flow-gradient plane and rotated azimuthally. In such configurations, the shear forces acting on the opposite edges produce a net moment, which scales as sin­(2θ)­sin­(2φ), as predicted by slip-ellipsoid theory (see Supporting Information for details).

To explore the influence of interparticle interactions on the suspension viscosity, we simulated suspensions with solid fractions in the range c ∈ [0.02, 0.22] by varying the number of particles while keeping the wall separation constant. For both fully-3D and quasi-2D particles, the effective dimensions were 2a ≈ 5 nm (length) and 2b ≈ 0.5 nm (thickness). Figurea shows the relative viscosity η/η_0_ as a function of the solid fraction c. The simulations were performed at a shear rate of γ̇ = 7 × 10^9^ s^–1^, corresponding to a Péclet number on the order of 10^2^. Our findings indicate that, even in the semidilute concentration regime, suspensions of slip particles exhibit a relative viscosity below unity (blue symbols in Figurea). However, the dependence of η/η_0_ on the solid fraction c becomes nonmonotonic: the viscosity initially decreases from 1, in agreement with the dilute-limit prediction, but then increases again beyond a critical concentration.

(a) Relative shear viscosity vs solid fraction for fully-3D and quasi-2D particles with a/b = 10. The crosses represent BI simulation results. (b) Reduced shear viscosity values corresponding to the data in (a). Dashed lines are least-squares fits to the data.

In Figureb, we plot the reduced viscosity (η/η_0_ – 1)/c as a function of c. For both slip and no-slip particles, the reduced viscosity was found to increase linearly with c. We calculated the linear and quadratic coefficients α and β in eq by fitting the MD data for the reduced viscosity to a line (dashed lines in Figurea,b). The obtained values of α and β are listed in Table. In our simulations, we used a channel height to particle length ratio of H/(2a) ≈ 2.2. Additional tests with H/(2a) = 1.5 and H/(2a) = 3.0 showed that increasing the channel height leads to a higher reduced viscosity. Nevertheless, for the tested solid fraction c ≈ 0.03, the reduced viscosity for slip discs remains negative even in the largest simulation domains considered.

1: Coefficients α, β and minimum viscosity concentration c min (where applicable) for graphene particles with a/b = 10

For quasi-2D particles, β is significantly larger for no-slip particles compared to slip particles. This difference stems from different particle–particle interaction behaviors: no-slip particles exhibit collective orbital motion, leading to frequent collisions, whereas slip particles primarily slide past one another with minimal interactions (see Supporting Information). Interestingly, the values of α for quasi-2D and fully-3D particles are quantitatively similar, but this is not the case for values of β. Fully-3D slip discs exhibit significantly larger β values compared to quasi-2D slip particles. This difference results from the enhanced tendency of fully-3D discs to aggregate, because of their increased rotational freedom and stronger edge-mediated interactions.

In addition to single-particle BI calculations, we compare our MD data with multiparticle BI data presented in Kamal and Botto? (crosses in Figurea). Both BI and MD show that the reduced viscosity increases with solid fraction at a slower rate for slip particles compared to no-slip particles. However, the multiparticle BI calculations predict smaller values of the interparticle interaction term, β, and do not capture the minimum in effective viscosity observed in the MD simulations at intermediate concentrations. These differences are due to the numerical treatment of near-contact particle–particle interactions in the BI model, where a short-range repulsive force was used to avoid overlap between the particles. Consequently, the results of the multiparticle BI calculations are most applicable to graphene nanosheets dispersed in good solvents (or in the presence of small-molecule dispersants), where adhesive interactions between the sheet-like nanoparticles are minimal.

The negative value of α for slip particles leads to a minimum in the suspension viscosity at c min ≈ −α/(2β). This concentration is smaller for fully-3D particles, consistent with their larger β values. These findings suggest that, for achieving substantial reductions in relative viscosity at finite solid concentrations using slip particles, it is essential to minimize β by suppressing adhesive interparticle interactions.

Conclusions

We have systematically investigated the shear viscosity of suspensions of nanometer-scale graphene sheets using MD simulations, complemented by continuum BI simulations of the Stokes equations. Our results demonstrate that, for sufficiently large Péclet numbers, graphene particles (sheets) with aspect ratios above approximately 5.5 reduce the effective suspension viscosity below that of the pure solvent. This provides a direct, molecular-level confirmation of prior continuum predictions and supports recent experimental observations of anomalous viscosity reduction in graphene-based suspensions. In the dilute regime, where particle interactions are negligible, MD and BI simulations are in excellent agreement, confirming that the phenomenon is driven by hydrodynamic slip and particle alignment rather than simulation artifacts or experimental uncertainties. The absence of negative intrinsic viscosity in a prior MD study? is attributed to the use of nanoparticles with insufficiently large aspect ratio. At higher concentrations, MD simulations yield a larger interaction coefficient (β) than BI simulations, primarily due to the aggregation tendency of graphene nanoparticles in water, a feature not captured in the continuum BI model.

In dilute suspensions, particles with both quasi-2D and fully-3D geometries yield comparable intrinsic viscosity values, validating the quasi-2D approximation commonly used in previous studies. ?,?,?,? In both geometries, slip particles stabilize with their surface nearly aligned with the flow-vorticity plane, whereas no-slip particles undergo continuous tumbling as described by Jeffery’s theory. The degree of alignment of functionalized graphene is known to decrease as Pe is reduced, ?,? so a negative α is only possible at large Pe, where Brownian effects are comparatively small. The threshold Pe for negative α also depends on particle aspect ratio, since the influence of Brownian fluctuations increases with particle length.?

The negative intrinsic viscosity arises from the interplay between the hydrodynamic slip length λ and the particle geometry, captured by two dimensionless parameters: the slip length to thickness ratio (λ/b) and the aspect ratio (a/b). When λ/b is sufficiently large and the particle is thin (a/b ≳ 5.5), increased slip induces a decrease in the steady alignment angle. Thin, well-aligned slip particles minimally disturb the surrounding fluid. Velocity gradients along their slender surfaces are substantially reduced compared to those in the undisturbed flow. Consequently, for slip particles (λ/b ≫ 1), increasing slenderness (a/b) results in a monotonic increase in the magnitude of negative intrinsic viscosity. This mechanism is analogous to the one observed in high-capillary-number, low-Reynolds-number bubbles, which also exhibit negative intrinsic viscosity? (a difference is that bubbles, unlike graphene, have effectively an infinite slip length).

Beyond shear viscosity, slip-induced alignment also affects the suspension’s normal stress response. The microstructure of the slip particles is anisotropic, but this anisotropy is not intrinsic to the two-phase medium (a change in the shearing direction would simply reorient the alignment direction), so the effective viscosity remains isotropic. However, unlike no-slip particles, the slip particles generate a significant difference in normal stresses between the flow and flow-gradient directions. In the Pe → ∞ limit, this normal stress difference is positive and increases with λ/b or a/b. An incomplete analysis of this effect for two-dimensional particles was reported previously,? but a more complete study with three-dimensional MD or continuum simulations would be valuable and could be carried out using the present framework.

For quasi-2D particles with slip, the interaction coefficient β is reduced by approximately an order of magnitude compared to no-slip particles. This reduction arises because the stable alignment of slip particles allows them to slide past one another while maintaining nearly parallel orientations. As a result, interparticle interactions remain weak even when the center-to-center separation is comparable to the particle length. In our simulations, particle–particle interactions are more pronounced for fully-3D disc-like particles than for their quasi-2D counterparts due to the increased rotational freedom of the 3D particles. Previous theoretical studies have reported the influence of slip length on the hydrodynamic interaction coefficient β in idealized systems. For example, Luo and Pozrikidis? computed β for pairwise interactions of slip spheres, while Kamal and Botto? considered quasi-2D circular cylinders with infinite depth. In both studies, slip has been found to reduce the effective interaction strength; however, those analyses neglect aggregation and anisotropic orientation effects, which are instead present in the current work.

The close agreement between MD and BI simulations indicates that our dilute-limit, high-Péclet-number findings are broadly applicable to other solvents and 2D materials beyond graphene. In particular, any well-dispersed system of plate-like particles where the key continuum parameters, λ/b (slip length to thickness ratio) and a/b (the aspect ratio), are sufficiently large will display the same slip-induced alignment mechanism and negative intrinsic viscosity reported here. Examples include hydrophobic layered materials such as h-BN, which exhibits large slip when immersed in a variety of polar and nonpolar liquids, and MoS_2_, which reduces boundary friction in various nonpolar oils. ?−? ? Moreover, engineered van der Waals heterostructures such as MoS_2_/graphene and MXene-graphene hybrids can exhibit lubrication performance superior to that of their individual constituents, owing to interfacial charge transfer that alters surface electronic structure and potentially reduces fluid–solid friction. ?−? ? Graphene can also exhibit much larger slip lengths in certain ionic liquids and nonaqueous solvents compared to water (see Table 1 in ref ?. for a survey of slip lengths across various graphene-liquid systems). However, in the semidilute and dense regimes, the surface properties of the particles, and thus their tendency to adhere to each other, an effect which is particularly important for sheet-like particles, become significant.

A key challenge in realizing negative intrinsic viscosity at technologically relevant concentrations lies in the strong tendency of graphene to aggregate, especially in water. Notably, the viscosity-reduction effect is most often reported for experimental systems based on lubricating oils and ionic liquids, which are known to give good dispersion stability. ?−? ? ? Future work should therefore investigate how molecular dispersants or alternative solvents influence both the dispersion stability and the suspension viscosity of 2D materials. Another important factor is flexibility: while the nanometer-scale graphene flakes studied here are relatively rigid, larger micrometer-scale sheets may experience significant bending, wrinkling, or self-folding under flow, ?−? ? ? ? potentially leading to a trend toward a positive intrinsic viscosity. Since slip reduces compressional viscous forces, the buckling transition becomes less significant for slip sheets than for no-slip sheets under practical shear rates.? These considerations highlight the complex interplay between flexibility, slip, and hydrodynamics in large graphene sheetsa topic that remains to be fully understood. Our results offer a framework to enhance our understanding of the hydrodynamic behavior of atomically thin materials.

Methods

Molecular Dynamics Simulation Details

The TIP4P/2005 model was used for water,? and the OPLS-AA force field was used for the graphene nanosheets. ?,? By default, carbon–water, water-wall and carbon-wall interaction parameters were calculated using the Lorentz–Berthelot mixing rules. In that case, the hydrodynamics slip length at the graphene–water interface was determined to be λ = 60 ± 11 nm following the method described in Herrero et al.? To model no-slip boundary conditions at the graphene–water interface, the Lennard-Jones interaction energy between carbon atoms in graphene and oxygen atoms in water was adjusted to suppress interfacial slip, yielding a slip length of λ ≲ 0.1 nm (see Supporting Information for details). An instantaneous vertical force was applied to the top wall, while the bottom wall was constrained in the vertical direction, ensuring a pressure of 1 atm in the solution. A Nosé-Hoover temperature thermostat ?,? was applied to the degree of freedom normal to the direction of the flow to maintain a constant temperature of T = 300 K in the solution. The system was equilibrated without moving the walls for 200 ps. The walls were then moved at a constant velocity along the x̂ direction for a second equilibration phase of 100 ps, after which production runs were conducted. During these runs, the force exerted on the walls, along with the positions and velocities of all atoms in the system, were recorded.

Continuum Boundary Integral Formulation

We employ a boundary integral method to solve the incompressible Stokes equations for a non-Brownian, isolated, rigid, quasi-2D particle freely suspended in the imposed flow field ** u ** ^ ∞ ^. Since the quasi-2D particle is effectively two-dimensional, the method requires the specification of a continuous line boundary $[eqn]$ onto which the integral equations are discretized. Molecular dynamics simulations of the flow field around a single-layer quasi-2D graphene in water suggest that $[eqn]$ can be well-approximated as a rectangle with semicircular edges. ?,? The Navier slip boundary condition is applied on $[eqn]$ , where the slip velocity is given by

$$[eqn]$$

where ** n ** is the outward-pointing normal vector, ** f ** is the hydrodynamic traction distribution over $[eqn]$ , λ is the slip length, and η is the dynamic viscosity of the suspending fluid.

In the two-dimensional boundary integral method, the incompressible Stokes equation is recast as an integral over $[eqn]$ . For a point $[eqn]$ , the boundary integral equation is ?,?

$$[eqn]$$

where ** G ** and ** K ** are tensors associated with the two-dimensional Stokeslet and stresslet, respectively? and dL(** x ) is a boundary element. We have assumed that the center of the particle is positioned at the origin, so that the particle rotates with angular velocity Ωŷ**.

Equation is solved numerically, along with the condition of zero net torque, to find ** f ** and Ω for a specific orientation angle φ. The numerical scheme is described and validated in ref ?.

In what follows, we will show that the intrinsic viscosity coefficient α can be evaluated in terms of ** f **. For a system of two-dimensional particles, α can be expressed as

$$[eqn]$$

where A, B are dimensionless coefficients and the angled brackets ⟨ ⟩ represent an average over the steady-state orientation distribution function p(φ).

The coefficients A, B depend on the particle shape and slip length. Using matrix transformations shown in refs ? and ?, the coefficients A and B can be decomposed as

$$[eqn]$$

where A p is the cross-sectional area of the particle. These coefficients are expressed in terms of the hydrodynamic stresslet tensor S _ ij _. This tensor is evaluated in terms of ** f ** as?

$$[eqn]$$

In this equation, S _ ij _ is evaluated in the particle frame (ŝ, t̂), where the unit vectors ŝ and t̂ are parallel to the platelet’s major and minor axes, respectively.

For Pe → ∞, ⟨1 – cos 4φ⟩ in eq can be evaluated analytically in terms of k e as?

$$[eqn]$$

Here, $[eqn]$ is the square root of the ratio between the torques exerted on a particle held fixed parallel (T(0)) and perpendicular (T(π/2)) to the flow. For no-slip quasi-2D particles, k e is always real? and follows a power-law relationship with the geometric aspect ratio.? Increasing λ, reduces k e so that k e = 0 at a critical slip length λ_c_ ∼ b. ?,? The cause for this reduction is that slip reduces the tangential traction over the slender region of the particle surface when the particle is held fixed in the direction of flow, which decreases the value of T(0). As λ/a →∞, the contribution to T(0) originating from the slender portion of the particle vanishes, and the resulting torque comes from the contribution from the edges. The contribution from the edges has an opposite sign to T(π/2), resulting in $[eqn]$ in this limit.?

Supplementary Material