Incremental viscosity by non-equilibrium molecular dynamics and the Eyring model

Theoretical treatments of the shear viscosity, η, of liquids started with empirical expressions for its temperature and pressure dependence. The effect of shear is to distort the molecular arrangement of the molecules which induces an opposing shear stress that is non-linear in the shear rate (SR) when the shear rate has large values. This gives rise to shear thinning, which has also been represented by semi-empirical analytic expressions. The molecular level factors that govern the shear rate dependent viscosity are still poorly understood,^1,2 and while the distortion of the local structure of a liquid under shear is now well characterized for small molecules,^3,4 theoretical prediction of the shear rate dependence of the viscosity⁵ has proved to be more problematic. Phenomenological models of shear thinning are still prevalent in part for this reason.^6,7

In 1936, Eyring derived an expression for the shear rate dependence of η based on what today might be called a “trap” model,⁸ which is based on assumed processes that occur at the molecular level. This Eyring model (EM) considers a representative molecule to be in a potential well whose escape barrier height is lower in the flow direction owing to the anisotropic structural distortion of the surrounding liquid caused by the shear. The nonlinearity of the stress vs. shear rate in the EM leads to the prediction of shear thinning. The temperature and pressure dependence of the Newtonian viscosity part can be included largely as a separate exercise in the expression for the prefactor. The appealing feature of this model is that the assumed origins of the shear thinning are intuitively reasonable and the parameters could in principal be related to molecular properties, therefore making it potentially molecule specific to a certain extent. The EM for shear thinning has been used widely in fitting to experimental viscosity data, such as in high pressure lubrication,⁹ and even for solids undergoing plastic flow.^10,11 We note in passing that more recently, other trap models have been derived which are more relevant to slowly evolving or ageing “quenched” systems such as glasses and gels.¹² 

Nonequilibrium statistical mechanics has developed over recent decades to provide useful analytic expressions relating to the shearing of liquids, some aspects of which can be implemented in molecular dynamics simulation. In the low shear rate limit, the Green-Kubo equation,¹³ or the equivalent Helfand-Einstein formula,¹⁴ can be used to determine the viscosity. Another significant theoretical development describes the stress build-up from an equilibrium state after the application of a shear rate of arbitrary magnitude, which is cast in the form of a transient time correlation function (TTCF).^15–17 In 1993, Evans, Cohen, and Morriss proposed a formula expressing the probability of states for a system driven arbitrarily far out of equilibrium, which is called the Fluctuation Theorem (FT).^18–20 In fact there are many slight variants of the FT formula which reflect different conditions applied to the system.^15,16,20,21 Non-Equilibrium Molecular Dynamics (NEMD) simulation is a complementary tool which can be used to compute these functions for model fluids and thereby provide definite predictions for specific chemical systems.

Constitutive relations for shear viscoelasticity in principle provide a complete description of the time dependent rheology of a liquid and have been widely investigated.^22–24 This work focuses on the steady-state stress relaxation function, C(t), or SRF which requires a different procedure to compute by NEMD, and has not been investigated far from equilibrium by this technique, as far as we are aware. The steady state shear SRF currently falls outside the scope of the above mentioned non-equilibrium statistical mechanical advances. The SRF is the transient viscoelastic response of a sheared liquid to a suddenly applied additional infinitesimal shear rate and gives directly what might be called an incremental viscosity, which could be considered to be a better measure of the system at a given shear rate than the usual definition of the shear-rate dependent viscosity (i.e., the shear stress divided by the shear rate). This is discussed in Sec. II. It is shown how these quantities can be obtained by NEMD and also within the framework of the EM. Inter alia analytic expressions are derived for the SRF and TTCF of the Eyring model. In the case of the SRF, insights are gained into why it has proved to be so successful in representing the shear thinning behavior of a diverse range of chemical systems. The extent to which the Eyring model is consistent with the FT is also explored here. It is shown that one of the limitations of the viscoelastic Eyring model derived here is that it does not allow for a time delay in structural evolution after the start of shearing, which could give shear stress overshoot (as is found in a number of experimental systems). A simple modification of the derived equations which has the same consequence as a “fictive” quantity to implement this time-delay is proposed which remedies this problem in a semi-empirical way.

To date, the fields of NEMD and associated non-equilibrium statistical mechanics of fluids and rheological constitutive equations have been advanced largely in parallel and independently. The present study is a preliminary attempt to analyze the links and potential synergy between the two approaches, with the ultimate objective of formulating constitutive equations that are consistent and up to date with the advances and conclusions derived from NEMD and related non-equilibrium statistical mechanics developments.

In this section, the concepts of an incremental viscosity, and related time dependent quantities for an arbitrarily large shear rate, are introduced. Their calculation and investigation using NEMD and within the theoretical framework of the Eyring model are also discussed, particularly in regard to the constraints required by the Fluctuation Theorem (FT) in the latter case.

The viscosity of a fluid in the zero shear rate limit, called the Newtonian viscosity, η₀, is defined by

where σ is the time averaged shear stress at shear rate $γ̇$ carried out at steady state (the sign convention is such that σ > 0 if $γ̇>0$⁠). In experiment, NEMD, and many theories, the viscosity at a finite shear rate is arbitrarily defined to be $η(γ̇)=σ/γ̇$⁠, which is inconsistent with the definition of Eq. (1) as it is not based on the incremental change in steady state shear stress in response to a small change in the shear rate. A potentially more useful alternative definition of the viscosity at finite shear rates which is consistent with Eq. (1) at any shear rate is what might be called the incremental viscosity, η_i, defined by

which is a more sensitive indicator of the physical state of the liquid at $γ̇$⁠. Figure 1 gives a schematic diagram illustrating the difference between η and η_i. As shown in Eq. (2), the shear rate dependent viscosity, $η(γ̇)$⁠, is an accumulated average of the incremental viscosities exhibited by the system up to that value of the shear rate and is therefore not a specific measure of the state of the system at that final shear rate, $γ̇$⁠, unlike the incremental viscosity, $ηi(γ̇)$⁠. The incremental viscosity can be obtained from $η(γ̇)$ through

As (usually) $dη/dγ̇<0$⁠, the incremental viscosity is less than the viscosity (i.e., η) at a given shear rate, a trend which is accentuated with increasing shear rate, as may be seen directly from the final identity in Eq. (3). The formula is a restatement of the formal definition of η_i, as $σ=γ̇η$⁠, which indicates that the incremental viscosity quantity is definable at any shear rate and is strongly dependent on $dη/dγ̇$⁠.

Some liquids exhibit an increase in shear viscosity after an initial shear thinning stage, an effect known as shear thickening. The viscosity, η, has a minimum at this transition from shear thinning to shear thickening.²⁵ Equation (3) indicates that when $dη/dγ̇$ changes sign from negative to positive, the incremental viscosity goes from being below to being above the value of η. The incremental viscosity definition therefore indicates clearly the shear thinning to thickening transition. The incremental viscosity cannot be zero, except in a hypothetical infinite shear rate limit, where η would also tend to zero.

So far the analysis has been for systems at steady state. In Sec. II A, the associated time dependent or viscoelastic properties associated with the incremental viscosity definition are considered.

A phenomenological way of representing the linear stress response to a time dependent shear rate is to use the Boltzmann superposition principle,^26–28 in terms of a time-dependent elastic modulus, C, or “shear stress relaxation function” (SRF), which has units of stress,

Equation (4) represents the case where the shear stress is on average zero for t < 0, i.e., it is an equilibrium system. In this work, the situation where there is a background or reference shear rate (which can be zero), $γ̇R$ already established at t = 0, is considered. Assume that the system is at steady state at $γ̇=γ̇R$ with a time average stress of $σ¯R$⁠. At time t = 0, the shear rate is changed to $γ̇(t)=γ̇R+γ̇kf(t)$⁠, which can be considered to be a shear rate “kick,” and where the unitless function, f(t), takes into account the time dependence of the kick shear rate magnitude, $γ̇k$⁠,

as the $γ̇R$ shear rate part of the shear stress is already at steady state, by assuming $γ̇R$ was initiated sufficiently far in the past, t < 0, for the transient start-up response to have finished before t = 0. Therefore,

The last step in Eq. (6), which involves replacing $γ̇$ by $γ̇R$ in C, is made because $γ̇k$ can always be made arbitrarily small so it does not significantly affect the structural and dynamical state of the system and hence the time dependence of C. Therefore we have

where δ(t′) is the Dirac delta function. The function, $ηi(t,γ̇R)→ηi(γ̇R)$ in the t → ∞ limit, as lim_t→0C(t) = 0. The last equation in Eq. (7) is the corresponding formula in discrete time steps.

For f(t′) = 1, we also have

The last step derives a widely used relationship in rheology,²⁹ which is obtained by making the substitution x = t − t′, so dx = −dt′. The important point is that far from equilibrium, the appropriate stress relaxation function to be used in the superposition formula is the one which gives η_i and not η, which underscores the importance of the incremental viscosity and related functions. To emphasize this point, the shear stress relaxation function, C, used in the kick procedure is denoted by C_i. In the $γ̇R→0$ limit, C_i(t) tends to the shear stress autocorrelation function of an equilibrium liquid.¹⁶ 

Equations (6) and (8) are key to the methodology used here. The applied shear rate $γ̇$ is the sum of that of the reference system (which can be zero) and a very small increment of the shear rate. The shear rate increment can be made sufficiently small not to affect noticeably the relaxation characteristics of the system compared to its reference state. The shear stress relaxation function so determined is then that of the reference system, which is what we want. The shear-kick (SK) technique is to probe the viscoelastic state of the reference system. Equation (8) is a useful identity which emphasises a related point. It is shown that this same relaxation function can be obtained when the shear kick is applied continuously during the perturbed trajectory. Both of these approaches are used in the NEMD calculations described in Sec. II B.

Equations (4)–(8) can be written in a form suitable for numerical integration by NEMD by replacing the integrations by summations. For example, in Eq. (4),

where Δt is the time step, t = MΔt is the time duration of the kicked trajectory, and $γ̇j$ is the shear rate at time jΔt. The formula in Eq. (9) can be used for numerical evaluation of computer generated C. Also after the last line in Eq. (7),

The NEMD procedure used here to calculate the incremental viscosity and its corresponding stress relaxation function, C_i(t), for a model simple liquid, and by EM, are presented below.

As shown above, the incremental viscosity can be calculated by applying an infinitesimally small shear rate “kick,” $γ̇k$⁠, at an arbitrary time declared to be t = 0 to an equilibrium (i.e., unsheared) fluid, or one that has reached steady state at a constant shear rate, $γ̇R$⁠. The viscoelastic stress response, σ_k(t), associated with the shear kick (SK) can be written in a form involving the associated shear stress relaxation function, C_i(t), of the fluid in the $γ̇$ state,

where $γ̇$ can be 0. The incremental viscosity is $ηi(t→∞,γ̇)$⁠.

Assuming a monatomic fluid in which the molecules interact with a radially symmetric pair interaction, ϕ, the virial expression for the shear stress of the system is

where N is the number of molecules in the simulation cell of volume, V. For particle i, the mass is m_i (all the masses are the same here) and r_α,ij is the α component of the vector $r̲ij=r̲i−r̲j$ between molecules i and j. Also, $vα,i$ is the α Cartesian component of the peculiar velocity of molecule i (i.e., the excess above any advected component).

A series of trajectories was initiated from equilibrium states or states at steady state shear (⁠$γ̇$⁠) which were well separated in time. The shear rate kick was either for the duration of a single time step or applied at every time step of the generated new trajectory. The function $Ci(t,γ̇)$ was obtained directly from the kick lasting just one time step, or alternatively by taking the time derivative of the stress response to a shear rate kick added every time step, and using the right-hand side equation in Eq. (11). The many-step and one step applications of γ_k correspond to the second and third lines of Eq. (7). This procedure can be used because (a) in the $γ̇k→0$ limit, $Ci(t,γ̇)$ will in practice be independent of $γ̇k$ and (b) the reference shear rate, $γ̇R$⁠, is constant during the kick procedure action. We refer to this small kick variant of NEMD as Kick NEMD or KNEMD, which in slightly different formulations has been applied in the past to model by NEMD shear flow starting from an equilibrium state in Refs. 30 and 31, and by Evans to determine the time dependent structure factor.³² 

The shear rate of the reference trajectory was $γ̇=dvx/dy$⁠, where x is the flow direction and y is the shear gradient direction. Two trajectories starting from the same reference state point were carried out repeatedly (typically 2000 times). Let $q̲3N$ represent the 3N position degrees of freedom and $v̲3N$ the corresponding quantity for the peculiar molecular velocities (the particle masses are the same). The origin of the coordinate system was at the center of the origin cubic MD cell. Let the microstate set of defining parameters be ${q̲3N,v̲3N,xdis,γ̇R}$⁠, where x_dis is the Lees-Edwards periodic boundary misalignment distance in the flow direction between adjoining cell boundaries and $γ̇R$ is the steady state shear rate of the reference system. The initial reference state is labeled Γ₊. The reference trajectory and the kicked trajectory with shear rate $γ̇k$ were carried on from this starting state point. Another reference starting state, Γ₋, was obtained from Γ₊ by implementing y-reflection mapping (e.g., see Ref. 16, p. 192), which here is ${qy→−qy,vy→−vy,xdis→−xdis,γ̇R→−γ̇R}$⁠, to which an additional very small shear rate, $γ̇k$⁠, was added to implement the kicked trajectory. The shear rates for either a time step or all time steps in the kicked trajectory were $γ̇R+γ̇k$ and $−γ̇R+γ̇k$ for the “+” and “−” starting states, respectively. The stress response needed in Eq. (11) is

where k stands for a kicked trajectory and 0 stands for one of the two conjugate paired reference trajectories. The stress contributions from the two conjugate pair reference states cancel out and the sum of the two was found to add to zero within machine error for all t, a deliberate feature of the procedure which made the number of significant figures defining the kick stress response greater than if conjugate pairing had not been implemented. The magnitude of the kick shear rate should be small enough not to significantly change the fluid’s structure and physical properties. The kick induced stress response function, $Ci(t,γ̇)$⁠, characterizes or probes the viscoelastic response of the reference fluid and therefore the kick must not perturb this state to any important extent. The stress response must be much smaller than the root mean square stress about the mean stress of the reference fluid.

NEMD simulations of the Lennard-Jones (LJ) fluid were carried out subject to imposed shear flow using velocity rescaling thermostatted SLLOD equations of motion^33,34 to obtain η_i and C_i for that system. The quantities quoted in this work are expressed in the usual LJ pair potential units of energy and distance and the mass, m, of the molecule. For compactness of notation, $γ̇R$ is replaced by $γ̇$ for the rest of this work as the difference between $γ̇R$ and the total shear rate during a kick trajectory can be made arbitrarily small.

Figure 2 shows the shear stress relaxation function, C_i(t), associated with η_i for a LJ liquid at the widely used state point, ρ = 0.8442 and T = 0.722, which is near the triple point.³⁵ Two types of shear history were used to compute $Ci(t,γ̇)$⁠, (a) by applying a shear rate “pulse,” $γ̇k=10−3$ for 1 time step at t = 0 and (b) for $γ̇k$ applied at every time step for t > 0 and using the second equation in Eq. (11) to obtain the SRF. As may be seen in the figure, the two methods give statistically indistinguishable curves for $Ci(t,γ̇)$ at each value of $γ̇$⁠. Figure 2 also has in it the time correlation function, C_i(t) = (V β)⟨σ(0)σ(t)⟩ (where V is the volume of the system, β = 1/k_BT, k_B is the Boltzmann’s constant, and T is the temperature), for the equilibrium (⁠$i.e.,γ̇=0$⁠) reference state,³⁶ which is seen to agree very well with that generated by the two kick routes when $γ̇=0$⁠. The trend is for the C_i(t) to decay more rapidly at long times with increasing the reference fluid shear rate, which is consistent with the observation that the viscosity decreases with increasing shear rate. The value of C_i(0) increases to a small extent with increasing $γ̇$⁠. Note also that for $γ̇=1.5$ and 2.0, C_i(t) is negative at intermediate times, which is not inconsistent with the fact that the stress autocorrelation functions at these shear rates are oscillatory. The figure also shows (green curve) a least squares fit of the correlation function to the Rounded Stretched Exponential (RSE) function defined in Eq. (3) of Ref. 37 and given in the caption of Fig. 2. This analytic form fits well the C_i(t) for $γ̇=0$⁠. In Ref. 37, the relaxation function is normalised to 1 at t = 0, whereas here the t = 0 value, C_i(0), is G_∞ instead.

Figure 3 compares the values of η and η_i as a function of $γ̇$ for the same system as used for Fig. 2. The figure shows that η_i decays with increasing shear rate more rapidly than η, as expected for a shear thinning liquid. Some of the η data in the figure are taken from Ref. 35 for system sizes of mainly 2048 LJ particles. The function,³⁸ 

where p₀, p₁, p₂, and p₃ are arbitrary variables, can be fitted to the NEMD η data quite well, as shown in the figure. The fit parameters are given in the figure caption. The viscosities, $ηi(γ̇)$ calculated from Eq. (3) and $η(γ̇)$ from Eq. (14) are also presented in the figure. This analytic IV is seen to be in good agreement with the KNEMD simulation data and within its statistical uncertainty. The agreement seen between the incremental viscosity obtained by (a) just applying the kick shear rate, $γ̇k$⁠, for the first time step of the segment (labeled “K1” in the figure) and then using Eq. (11) and (b) adding $γ̇k$ for every time step of the kick (“KA” in the figure) gives strong support for the applicability of Boltzmann’s superposition assumption even when $γ̇$ is arbitrarily large.

In Sec. II C, exact formulas are derived for the incremental viscosity and related stress relaxation functions for the Eyring model.

The Eyring model (EM) for liquid shear thinning is extended below to encompass the intrinsic viscosity, associated stress relaxation function, and the transient time correlation function (TTCF).

The steady state shear thinning curve of a liquid can be represented by the generic analytic form

where σ is the shear stress at a given shear rate, σ₀ is a characteristic shear stress above which the liquid starts to shear thin to a noticeable extent, and F(⋯) is a generic function. The term, $σ0$⁠, is a somewhat imprecisely defined quantity and is often obtained by fitting experimental or NEMD simulation shear thinning data to the analytic form given in Eq. (15) and treating σ₀ as a fit parameter. Clearly F(σ/σ₀) → σ/σ₀ in the low shear rate limit as η → η₀ in this limit. A specific example of Eq. (15) is the Eyring formula,^8,9

where F(σ/σ₀) = sinh(σ/σ₀). In this case, σ₀ = k_BT/2$v0$⁠, where $v0$ is a characteristic volume traced out by the flowing “unit” (again not well defined on a molecular scale) during a cage of the molecule reorganisation event or “jump.” By rearrangement of Eq. (16),

where τ₀ = η₀/σ₀ is a characteristic relaxation time and $σ=η0γ̇$ has been used. Note this is not the Maxwell relaxation time, τ_M = η₀/G_∞, where G_∞ is the infinite frequency shear rigidity modulus. The Newtonian viscosity acts as a coefficient or reference state whose temperature dependence is often assumed to follow an activation energy or Arrhenius dependence although there are a number of alternative assumed flow mechanisms proposed in the literature which also give rise to the same inverse temperature dependence but do not invoke the concept of activated dynamics. Notably there is the so-called “shoving model,”^39–43 in which stress relaxation is taken to be caused by elastic deformation and ensuing structural rearrangement of the cage of molecules around an arbitrary molecule. The Arrhenius activation energy in the shoving model is replaced by the product of G_p, the plateau modulus (which is slightly lower than the true infinite frequency modulus, G_∞), and Ω_c, a characteristic volume of order of the molecular volume or smaller.

What is not so well known is that Eyring’s shear thinning formula can also be applied to predict time dependent shear rate viscoelastic effects such as during shear-start-up,

when it is incorporated in Maxwell’s equation including time which leads to viscoelasticity,⁴⁴ 

which expresses the total strain rate as the sum of elastic and viscous parts, the first two terms on the left of Eq. (19), respectively. Note that the infinite frequency shear rigidity modulus, G_∞, is assumed to be independent of $γ̇$ in the Eyring model, which is usually not a bad approximation.³⁵ 

Substitution of Eq. (16) in Eq. (19) gives the fundamental equation used in our extension of the EM,

The analytic solution of σ(t) from Eq. (20) is presented in Eq. (A8) derived in  Appendix. Figure 4 shows what might be called a “time dependent instantaneous viscosity,” $η(t)=σ(t)/γ̇$⁠, normalized by the Newtonian value, η₀, for several shear rates, starting from an equilibrium (i.e., σ = 0) state. Figure 5 is as is shown in Fig. 4 except that the starting stress at t = 0 corresponds to $γ̇=1.5$ and the shear rate for t > 0 is 1.5 + x where x is the number given in the figure. This results in some of the curves representing the effect of a decrease in shear rate at t = 0. The figure demonstrates that Eq. (A8) is capable of predicting the stress response for a step-down as well as a step-up in shear rate, and is one of the main results of this work.

The Transient Time Correlation Function (TTCF) for shear stress, $CTr(t,γ̇)$⁠, is^16,47,48

where ⟨⋯⟩ is a time correlation function, which is the product of the stress taken from an equilibrium (i.e., unsheared) system at t = 0 and one at a time t, σ(t), which has evolved at constant shear rate, $γ̇$⁠, during that period (note the shear stress is positive for, $γ̇>0$⁠).

The stress build up after the commencement of shear is more complicated if expressed in terms of an alternative stress relaxation formulation formalism used in the rheological literature, in which the SRF is usually expressed as a weighted sum of exponentials with different relaxation times. The current relaxation times are a function of the present physical state of the system as it evolves from an equilibrium to a sheared steady state. If we assume that the current shear stress relaxation function at time t′ is G_∞ϕ_s(τ_s(t′)), where τ_s is a characteristic shear stress relaxation time, then the time dependent stress from shear start-up at an arbitrary shear rate would be⁴⁵ 

which is more complicated than the TTCF approach given in Eq. (21). The two treatments in Eqs. (21) and (22) are both correct, but the two formulations are in terms of different time dependent functions.

Figure 6 shows $CTr(t,γ̇)$⁠, normalized by G_∞, for several $γ̇$ values using Eq. (A10) for the Eyring model. The function decays more rapidly with time with increasing shear rate. Substitution of Eq. (A9) in the definition of η_i in Eq. (2) gives for the Eyring model [i.e., where Eq. (17) applies],

which is, in fact, a special case of the Carreau shear thinning formula (i.e., with the exponent of the denominator set to 1/2),⁴⁹ sometimes used to represent the viscosity, η, itself. An advantage of the incremental viscosity in the present case is that the Carreau model (with α = 1/2) has a point of inflection at a shear rate of 0.71/τ₀ whereas the Eyring model has a point of inflection at 1.035/τ₀. This should facilitate a more accurate extrapolation of the shear rate dependent viscosity to zero shear rate and hence to the Newtonian viscosity (note that η and η_i have the same value in the $γ̇→0$ limit). The relationship between η and η_i for other simple analytic forms for η is discussed further in the supplementary material.

Figure 7 presents a comparison of η and η_i as a function of the shear rate for the Eyring model using Eqs. (17) and (23), respectively. The incremental viscosity, η_i, is seen to decay more rapidly with shear rate than η, as would be expected from the defining equations in Eq. (3).

The incremental viscosity can be calculated if an infinitesimally small shear rate “kick,” $γ̇k$⁠, is added at t = 0 to an equilibrium or steady state sheared system at shear rate $γ̇$⁠. Now, σ_k(t) is the difference in the shear stress between the kicked and reference trajectories and therefore σ_k(0) = 0. The viscoelastic response associated with the shear-kick can be written in a form involving a shear stress relaxation function appropriate to the $γ̇$ state,

Figure 8 shows the Eyring model incremental viscosity stress relaxation function, C_i, given in Eq. (A11) for several shear rates, revealing that it decays more rapidly with time with increasing (reference state) shear rate. The slower relaxation processes are more affected by the shearing mechanism, which is the same qualitative trend as found in past NEMD studies.⁵⁰ Comparison between Fig. 6 and Fig. 8 reveals that the $Ci(t,γ̇)$ are more separated than the $Ctr(t,γ̇)$ for the same $γ̇$⁠, notably at short times.

See the supplementary material for discussion of other model analytic forms for η found in the literature and their derived η_i.

Section III deals with the extent to which the NEMD data and Eyring model comply with the Fluctuation Theorem.

The average stress, σ_t, over a single time period from 0 to t is

Let us assume that the set of σ_t follows a shifted gaussian probability distribution,

where $σ¯$ is the simulation average stress taken over many segments and σ_D, is the standard deviation (SD) of the quantity, $[σt−σ¯]$⁠, over that period (⁠$>>$t). Then,

For an equilibrium system,^14,51 the Einstein-Helfand formula for the viscosity (⁠$i.e.,σ¯/γ̇$⁠) is, for times t, greater than the Maxwell relaxation time, with $σ¯=0$⁠,

where ⟨⋯⟩ indicates an average over a sequence of segments of length t. We assume here that it is still reasonably accurate for finite but not too large $σ¯$⁠. Equation (28) substituted in Eq. (27) (assuming the shifted gaussian stress PDF is applicable at a finite shear rate) gives the Fluctuation Theorem (FT) expression,

Despite making a number of assumptions, the above short approximate derivation of the FT for shear does lead one to expect that a (shifted) gaussian distribution of time-averaged shear stresses should still be compatible to a reasonable approximation with the FT. Note however that the FT does not rely on the PDF being gaussian. The FT is not able to provide the analytic form of $η(γ̇)$ but any model system that gives rise to a particular rheological constitutive equation should comply with the FT. Although Eq. (27) is only formally exact in the equilibrium (⁠$i.e.,γ̇→0$⁠) limit, it is a reasonable approximation for not too high shear rates, especially for long segment averaging times.^52,53 Deviations from the gaussian form, however, have important consequences which, for example, limit the use of the Green-Kubo expression for the viscosity to small values of the shear rate.⁵² 

A spatially local stress can be defined, and the effect of the stress sampling volume on the PDF for finite shear rates is investigated here, to explore what the analytic form of this stress distribution actually is. Basing this definition on Eq. (12), the contribution from molecule i to the virial expression for the shear stress of the whole periodic system is

where the average number density is ρ = N/V. The virial expression for the stress of the whole system is $σ=∑i=1Nσi/N$ and for a subvolume Ω in the simulation cell, containing on average N_Ω molecules,

which is the average of σ_i for those molecules inside a smaller volume Ω of a chosen shape. This Atomic Virial (AV)⁵⁴ definition is just one of many possible prescriptions for the local stress. It has the advantage that it is physically reasonable, being based on those molecules located in the subvolume, and it goes over to the usual virial expression when Ω → V, the volume of the simulation cell. The AV method is suitable for spatially homogeneous systems but gives rise to unphysical oscillatory stress fluctuations near walls, in which case the exact Method of Planes (MOP) method⁵⁵ and Volume Averaging (VA)⁵⁶ methods should be used.

The logarithm of the probability distribution functions (PDF) of σ_Ω for equilibrium and sheared liquids is presented in Fig. 9 using a spherical subvolume, Ω (see Refs. 54 and 57 for further details). The average number of molecules in the subvolume is indicated in the figure. The $γ̇=0$ distributions increasingly depart from a gaussian form as the sampling volume or N_Ω decreases. For not too small N_Ω, the main effect of the shear rate is to shift the gaussian along the abscissa by the average shear stress, $σ¯$⁠, value. For smaller N_Ω, the distribution adopts a noticeable skewed peaked shape with prominent tails in the high stress wings. These trends are the same for a cubically shaped subvolume. Figure 10 shows the PDFs as a function of the shear rate for spherical and cubic subvolumes of volume, 4πR³/3 and 8R³, respectively, where R = 2.41. The gaussian shape is reasonably well maintained but the standard deviation increases with $γ̇$⁠. This partly vindicates an early phenomenological molecular theory of viscosity for simple liquids which made this assumption.⁵⁸ 

The shape of the PDF can be characterized using the second to fourth moments of the shear stress distribution about the mean, σ₂, σ₃, and σ₄.⁵⁹ The standard deviation (SD), $σ2$⁠, skew (SK), $=σ3/σ23/2$⁠, and kurtosis (KU), $X4=σ4/σ22−3$⁠, measure the breadth, degree of slant, and extent of the tails, respectively, of the PDF. Figure 11 shows that these quantities increase with decreasing Ω, which is consistent with the PDF shapes seen in Fig. 9. The SD, SK, and KU increase in magnitude with diminishing subvolume, to a good approximation as ∼Ω^−ν, where ν is an exponent. There is a gradual departure from a gaussian form as the size of Ω decreases. Searles and Evans showed that the instantaneous and time averaged shear stress for the whole of a sheared liquid is close to being gaussian.⁵² The trends in Fig. 11 are consistent with the higher moments being larger than those of a gaussian.^60–64 The effect of shear is also to skew the PDF to the right. The figure also shows that the data for sphere and cubes fall on the same lines, so the magnitude of the subvolume, Ω, rather than its shape (within reasonable limits) is the key governing factor for this local stress model.

The EM assumes that, on the one hand, there is only one shear stress and one shear rate in the system, while on the other hand, the flow process involves molecular “jumps” which would imply strong local fluctuations in both quantities. The two scenarios sit rather uncomfortably within the same model. This issue is considered in this subsection and partly resolved. A modification of the EM is derived in which the effects of local shear stress fluctuations are incorporated through the introduction of a gaussian shear stress PDF.

Let the time and spatial average value of the local shear stress, σ_t, be $σ¯$⁠. The EM can be generalized by assuming that there is an assembly of such local states at any given time characterized by σ_t, each obeying

where A_u = σ₀/η₀. The σ_t distribution is taken to be gaussian from Eq. (26). The system is deemed to explore a distribution of P(σ_t) and hopping rates, which is in fact the Ree-Eyring extension to EM but with the average flow a sum of stressed units, rather than the other way around.^65,66 Previous shear stress relaxation models^67–69 have assumed such processes, which are also consistent with stretched exponential stress relaxation.⁷⁰ There is then a distribution of shear rates at any one time, whose average value over time and space is

The latter step in Eq. (33) was performed by integrating the two exponential terms in the sinh function separately and using $∫−∞∞⁡exp(±ax)exp(−b(x−c)2)dx=π/bexp([a2/4b]±ac)$⁠.⁷¹ In the zero shear rate limit, from Eq. (33),

to give $Au=exp(−σD2/2σ02)σ0/η0$⁠, which when combined with Eq. (33) leads to

This is the same analytic form as the original Eyring relation of Eq. (32) with the single stress, σ_t, replaced by the system’s average value, $σ¯$⁠. Therefore, even though there is a stress PDF which is much broader than the (Eyring) delta function distribution, the analytic form of the final formula relating average shear stress and shear rate is the same. The fact that the Eyring analytic form is consistent with a realistic distribution of stresses on the molecular scale might in part explain its success in reproducing finite shear rate experimental data, as σ₀ is almost always treated as a parameter to be fitted to experimental or simulation data.

The above treatment does not indicate how large these flowing units are. Insights into this can be gained for the present LJ system by considering a spatially resolved shear stress correlation function,^72–74 

which measures correlations between the configurational stress experienced by two particles i and j a distance r apart. The kinetic component of the stress is not included in the above definition, but for liquid-like densities, and not too high T, it is relatively small compared to the interaction part of the stress. Unlike g(r), it is defined here so that it is not normalized by 4πr²Δr and as a quantity is more compatible with the number of particles in a shell of thickness Δr at r. Figure 12 shows P_xy(r) for a dense equilibrium fluid state (the simulation details are given in the figure caption). Note that if P_xy(r) = 0, there is no shear stress correlation between particles i and k at that separation. The figure shows that stress correlations reflected in the oscillations extend to at least the third coordination shell as found also in Refs. 72 and 73, which contains of order 100 particles. The number of particles within a sphere of radius r, or n(r), is also given in the figure. Such long range correlations in the local shear stress might have been anticipated from the work of Stassen and Steele who decomposed the stress time correlation function used in the Green-Kubo formula for the Newtonian viscosity into its two-, three-, and four-body components.⁷⁵ They were all shown to be important and to decay slowly at long times.⁷⁵ This behavior also explains why there is a strong system size dependence of the shear viscosity of high density fluid states, especially near the triple point.⁷³ Figure 12 suggests that for simple monatomic liquids, the number of molecules involved in the flow process is rather large (∼100), which is probably much larger than would be the case for large organic molecules. The relative success of the Eyring model for large molecules is possibly in part attributable to a more localized stress field than for simple monotomic liquids.

Simple molecular fluids (e.g., Lennard-Jones) are one type of system where the Eyring equation has had mixed success in predicting the shear thinning curve.^6,76 This is illustrated using the present data in Fig. 13, which shows the shear thinning curve for a LJ system close to the triple point (note the log-log scale). The Eyring model does not match the simulation data for large shear rates, where a power law is better. This is possibly because the shear rates are relatively high (∼10¹² s⁻¹ for liquids composed of small molecules such as argon) for LJ molecules and long range structural changes induced by the shearing action are more readily achieved due to their spherical shape, which is expressed in its most extreme form at very high shear rates by the appearance of the so-called “string” phase.⁷⁷ The possibility of long-range reordering due to shear was not envisaged in the construction of the Eyring model, only local small shear deformations of the structure. In addition, for a LJ liquid or real liquid composed of small molecules, a significant proportion of the area under the shear stress relaxation function (and hence viscosity) at short times is due to kinetic or inertial stress decay, whose time expansion is in even powers of time [unlike the Eyring model stress relaxation decay which is based on exponentials; see $σ(t,γ̇)$ in Eq. (A8) and derived functions].

Another limitation of the Eyring model is that it assumes that for large shear rates, the stress continues to increase with shear rate, whereas in some studies (e.g., for glasses⁷⁸), a limiting stress value has been found.

Another limitation of the Maxwell-Eyring model of viscoelasticity is that it does not predict shear stress overshoot (where the shear stress exceeds at short time its long time value) at high shear rates and soon after the start-up of shearing, an effect observed in polymer solutions, melts,⁸⁰ metallic glasses,⁷⁹ and even in simple liquids.⁵⁰ This absence is because the instantaneous relaxation time soon after the application of a step in the shear rate is calculated from the steady state shear formula [see Eq. (20)], which is questionable as the system may not have had time to come to a structural and rheological steady state at that shear rate during this transient period. A real system takes time to evolve its structure (and hence stress) to that of the final steady state. Therefore at short times, the EM underestimates the structural relaxation time. The real system’s characteristic structural and property values would lag behind those predicted using the Maxwell-Eyring viscoelastic model, possibly resulting in a shear stress overshoot.

In the field of glass formation, a so-called “fictive” temperature, T_f, has been introduced to quantify such time lag effects. The fictive temperature is that of an equilibrium liquid which has the same structure and thermodynamic properties as a supercooled state. The fictive temperature is higher than or equal to the actual temperature of the system, i.e., T_f ≥ T during rapid cooling.⁸¹ This takes account of the lag in structural evolution as the system is forced out of equilibrium.

A similar construction can be adopted for the shear stress evolution as introduced in Ref. 45 to represent the physical state of a lubricant fluid element as it passes through the contact zone in high pressure elastohydrodynamic lubrication. The stress relaxation function in that study was represented by a stretched exponential in time, which can be written as a weighted sum of exponentials with a distribution of relaxation times. In the present simpler analytic model, there is only one effective relaxation process representing the system at any given time after shear-start-up [see Eq. (19)]. This time lag effect can be introduced empirically in the Eyring model by the inclusion of a time-delay function, $Fi(t,γ̇)$⁠, which has the limits of 0 and 1 at t = 0 and t = ∞, respectively. From Eq. (19),

The performance of one form chosen for $Fi(t,γ̇)=1−exp(−ϵγ̇t)$⁠, where ϵ is a positive constant, is seen in Fig. 14. The figure shows the case where shear start-up from an equilibrium liquid is initiated with $γ̇=1$ and ϵ = 4 and using the same Eyring parameters as for the previous figures. The functional form of $Fi(t,γ̇)$ and value of ϵ are physically reasonable as the peak in the overshoot stress profile is typically in the strain (⁠$i.e.,∫0tγ̇(x)dx$⁠) range between 0.1 and 1.0. If F_i = 1, we have the original Eyring model, and for the adopted time-delay functional form, this case might be called the fictive Eyring model. Although there is no specific fictive variable in Eq. (37), its effects are encapsulated in the additional semi-empirical term, F_i. Figure 14 shows that this model exhibits a peak in the stress at short times, unlike the original Eyring model itself (i.e., where $Fi(t,γ̇)=1$⁠) also shown in the figure, which always increases monotonically with time after shear start-up. An alternative approach to generate a time delay effect would be to relate the instantaneous relaxation time to an order parameter which evolves more slowly than, in this case, the Eyring shear stress. For temperature changes in the glassy state, a fictive temperature has been used to implement this time delay.^82,83 In Ref. 83 this was interpreted in terms of a so-called material time which is an effective time which is moving slower than the actual time. For a pressure jump, a fictive free volume has been used.²² In the present case, this order parameter would be evolved parallel in time with the shear stress itself. Further experimental and NEMD studies would be needed to identify a suitable candidate for this property.

Rheometers measure shear stress against shear strain or strain rate, $i.e.,σ(γ̇)$⁠, and the invocation of a “shear viscosity” may not be a necessary or desirable exercise as it is certainly not uniquely defined apart from in the Newtonian regime where the shear rates are very small on a molecular time scale. The usual definition of the non-Newtonian viscosity is the ratio of the shear stress divided by the shear rate (⁠$i.e.,η(γ̇)=σ(γ̇)/γ̇$⁠). We propose an alternative definition of a shear rate dependent viscosity, the derivative of shear stress against strain rate, which we call the “incremental viscosity” (IV) or $ηi(γ̇)$⁠. The IV is a better representation of the physical state of the system at a particular shear rate than $η(γ̇)$⁠. Also η_i converges to the Newtonian viscosity as the shear rate tends to zero, and because the IV is compatible with Boltzmann’s superposition principle, it is the natural extension of the Newtonian viscosity definition to systems where the shear rate is arbitrarily large. The IV is calculated here by NEMD using a time dependent shear rate or “shear-kick” scheme and analytically for the Eyring model where it was found to be a special case of the Carreau formula which has often been used to represent liquid $η(γ̇)$⁠. The associated stress relaxation function for the shear kick is also computed and derived within the NEMD and Eyring models.

The Eyring model formula is shown to apply also when a gaussian distribution of stressed regions in space is assumed, which is approximately compatible with the fluctuation theorem. This may in part explain the Eyring model’s success in reproducing experimental data for a wide range of liquid classes. A simple semi-empirical modification of the Eyring model for the time dependent stress after shear start-up, to include the effects of time delay in the system’s response, is also shown to lead to shear stress overshoot behavior which is often seen in glassy and polymer solution systems.^80,81

The concept of the incremental viscosity unifies the rheological description of both Newtonian and non-Newtonian regimes, as the viscosity in both cases is defined as the derivative of the stress with respect to shear rate. The Newtonian viscosity has to be defined that way as any shear rate dependence of the ratio of stress to strain rate (i.e., not a linear relationship) means that one is not in the Newtonian shear rate regime. In practice, the incremental viscosity could be obtained experimentally by the direct application of a small shear rate pulse or an oscillatory strain rate to an already sheared system, as in the NEMD calculations reported here. Alternatively, the incremental viscosity could be computed by fitting experimental or simulation shear stress vs. shear rate data to a function and then analytically differentiating this with respect to shear rate.

We would like to add that the usual definition for the viscosity (⁠$i.e.,σ/γ̇$⁠) has advantages, for example, in being directly measurable by experiment or NEMD simulation, and the stress of the system is readily obtained from it. The incremental viscosity has complementary advantages, however, in probing more accurately the physical state of the system at a given shear rate and in being compatible with viscoelastic constitutive equations based on Boltzmann’s superposition principle.

See supplementary material for further discussion of the relationship between η and η_i for other simple analytic forms for η found in the literature.

D.D. would like to acknowledge the support received from the EPSRC under the Established Career Fellowship Grant No. EP/N025954/1. All data can be made available by emailing the authors of the paper or [email protected].

The analytic solution of Eq. (20) for σ(t) is given here.

Let t = 0 be the time an equilibrium liquid is subjected to a constant shear rate. An analytic solution of Eq. (20) was derived in Ref. 45 via a Riccatti equation to give

where

An alternative more general solution for σ(t) resulting from the following strain rate history,

is derived here. This represents the transformation of the system with time from a steady stressed state with shear rate, $γ̇1$⁠, at t = 0⁻ to eventually a new steady state at $γ̇2$⁠, of arbitrary relative magnitude and sign. This is derived using the substitution, Y(t) = exp(σ(t)/σ₀) in Eq. (20), which after rearrangement gives

where A₁ = −G_∞/2η₀, $B1=γ̇G∞/σ0$⁠, and C₁ = G_∞/2η₀ are time independent parameters. Equation (A4) can be integrated by separation of variables to give⁴⁶ 

where

and G is the constant of integration. Therefore, if H = exp(SG),

where Y(0) = exp[σ(0)/σ₀] is the value of Y at t = 0, the starting point being a steady state system with stress, σ(0) = σ₁. From Eqs. (A5) and (A7),

which is the time dependent stress during the transition between the two steady states 1 and 2. Equation (A8) has the long time limit,