Communication: Simple liquids’ high-density viscosity

The factors that affect the steady-state shear viscosity, η, of a fluid have been the subject of much investigation, yet no first-principles expression for the state-point dependence of η exists that allows for a straightforward calculation of η from the system’s Hamiltonian. This means that one cannot easily calculate the viscosity’s dependence on temperature T, particle number density ρ ≡ N/V, and pressure p for dense fluids. The literature is centered around semi-empirical expressions based on parameters such as the free volume that are not uniquely defined for real atoms and molecules, and viscosity data have typically been fitted to empirical polynomial expressions not obviously linked to the underlying physics.^1–5 The density and temperature ranges where they are applicable are also limited.

A significant early development in providing an accurate expression for η is Enskog’s expression for hard spheres,^6,7

Here σ is the sphere radius and B₂ = 2πσ³/3 is the second virial coefficient. The viscosity prefactor $η0HS$ is given by

The mass of the sphere is m and k_B is the Boltzmann constant. The value of the radial distribution function at contact of the hard spheres is given by g(σ) = (Z − 1)/(B₂ρ) in which Z = p/ρk_BT is the so-called compressibility factor.⁸ At liquid-like densities, the last term on the right-hand side of Eq. (1) dominates and is, apart from a numerical factor, the same as derived by Longuet-Higgins and Pople⁹ based on simple kinetic arguments.

The dimensionless viscosity η* in Eq. (1) involves the effective molecule radius σ. This concept appears also in the expression for the fluidity (1/η) of Dymond¹⁰ for high-density fluids in which the fluidity is taken to be linear in the volume per particle, $v$⁠. A linear dependence on $v$ was also considered by Hildebrand¹¹ and by van Loef,^12,13 who showed that it fits well to the viscosity of dense fluids composed of a wide variety of molecules. The molecular diameter is typically expressed in terms of the close-packed volume per particle of the crystal. This prescription has been found to reproduce well also the fluidity of model inverse-power fluids.¹⁴ 

An alternative definition of the characteristic length scale is ρ^−1/3. This has been used in scaling onto a single curve the viscosity trends for high-density fluids¹⁵ and hydrocarbon fluids.^16,17 Andrade in the 1930s proposed an early application of this type of viscosity scaling, for molten metals,^18–20 which lead to the “macroscopically” reduced (dimensionless) viscosity,

Even though σ and ρ^−1/3 are numerically close for dense liquids, this substitution is a significant conceptual leap. It makes the transition from molecular-potential focused equations for transport coefficients, which trace their origin to the kinetic theory of gases and the hard-sphere model, to ones based on macroscopic variables independent of molecular detail and characteristics. This treatment is an important step towards including transport coefficients within the category of thermodynamic state functions.²¹ 

The idea that viscosity and thermodynamic properties are subject to the same treatment is implicit in a number of previous studies.^21–23 Hoover et al.²⁴ showed theoretically that the density and temperature properties of model fluids formed by the inverse power-law pair potential, $v$(r) ∝ r⁻ⁿ, collapse onto a single curve where ρ^n/3/T is constant. Alba-Simionesco, Roland, and their respective co-workers more recently showed that the experimental reduced dynamical properties of many glass-forming molecular liquids can be collapsed in the same way, treating the exponent n as a fitting parameter.^25–27 The recent body of work on isomorph scaling has confirmed this assumption and provided a theoretical basis.^28–30 With isomorphs defined as the lines of constant excess entropy in the thermodynamic phase diagram,^28,31 the reduced viscosity $η̃$ is constant along an isomorph because the reduced-unit dynamics is. This fact is not reliant on a particular form of the potential as long as the system has strong virial potential-energy correlations at the state points in question.²⁸ In particular, it is not restricted to inverse-power potentials or even to pair potentials, and several molecular systems have been shown to conform to the isomorph theory in computer simulations^32,33 and experiments.^34–36 

Not all systems have isomorphs, though. While most metals and van der Waals bonded liquids are believed to have isomorphs in the condensed-phase part of their phase diagram, systems with significant directional bonding like covalently or hydrogen-bonded systems are not expected to have isomorphs.^32,37,38 The Lennard-Jones (LJ) liquid belongs to the former class of the so-called Roskilde (R)-simple systems,^31,39–43 which includes most systems with more or less spherical interaction symmetry, the property that traditionally defines a “simple” system.⁸ We use below the isomorph theoretical framework as a guide to arrive at an expression for how simple liquids’ viscosity varies throughout the high-density part of the thermodynamic phase diagram. A novel viscosity expression is suggested with reference to the LJ system and shown to describe well experimental data for methane and argon as well as simulation data for an effective pair-potential model of liquid sodium.

R-simple systems are characterized by a potential-energy function U(R) obeying the following scaling condition³¹ (in which R is the 3N coordinates of the N particles and λ > 0 is a scaling factor):

Equation (4) is strictly obeyed only for an Euler-homogeneous potential plus a constant, but it is a good approximation for several realistic atomic and molecular models.^31,44 For such systems, the scaling property is rarely obvious from the mathematical expression for U(R) for which reason the term “hidden scale invariance” is sometimes used.^{28,33,45–49}

The R-simple region of a given system is identified in computer simulations as the state points characterized by strong virial potential-energy correlations in constant-volume canonical-ensemble fluctuations.^28,50,51 For the LJ liquid, our simulations have shown that the Pearson correlation coefficient obeys R > 0.94 at all state points with density and temperature higher than those of the liquid at the triple point. Whenever Eq. (4) applies, the isomorphs are to a good approximation invariance curves for both structure and dynamics^28,31 if quantities have been made dimensionless using the “macroscopic” unit system with distance measured in units of ρ^−1/3, energy in units of k_BT, and time in units of $(m/kBT)1/2ρ−1/3$⁠.²⁸ 

In the simplest version of the isomorph theory,^28,52 a function of density h(ρ) exists for any R-simple system such that its isomorphs are given^47,52 by

For systems in three dimensions with pair potentials of the form $v$(r) = ∑_n$vn$r⁻ⁿ, the function h(ρ) has one term for each inverse power-law term, h(ρ) = ∑_nC_nρ^n/3.⁵² For the LJ pair potential $vLJ(r)=4ε(r/σ)−12−(r/σ)−6$⁠, one has if γ₀ is the so-called density-scaling exponent determined from a simulation at a single state point with density ρ₀⁴⁷ 

Because $η̃$ is isomorph invariant, Eq. (5) implies that $η̃$ is a function of h(ρ)/T. If the temperature variation along a reference isomorph as a function of density is denoted by T₀(ρ), Eq. (5) shows that h(ρ) ∝ T₀(ρ). Thus one can write $η̃(ρ,T)=F(T0(ρ)/T)$⁠.

The freezing line, which marks the limit between liquid and solid-liquid coexistence regions in the density-temperature phase diagram, defines a temperature function T_F(ρ) in which $ρ$ is the liquid density. The freezing line is an approximate liquid-state isomorph,^28,53,54 implying that the reference isomorph function T₀(ρ) may be taken to be T_F(ρ),³⁰ i.e.,

Rosenfeld applied this temperature scaling for the single-particle diffusion constant and the viscosity of the Yukawa system,⁵⁵ which is now known to be R-simple.⁴⁹ The function F in Eq. (7) is a priori system dependent, but the data presented below interestingly suggest that the same functional form applies for all simple liquids.

Simulations were carried out using Roskilde University Molecular Dynamics running on graphics processing units⁵⁸ (details are given in the supplementary material). At each state point in Fig. 1(a), a SLLOD simulation was run to evaluate the viscosity. Figure 1(b) shows viscosity data for the LJ liquid at the state points given in Fig. 1(a). Upon heating at constant density, the viscosity initially decreases and then increases. The existence of viscosity minima along isochores is well known.⁵⁹ 

In order to validate Eq. (7) we proceed to identify an expression for the LJ system’s freezing temperature as a function of density. A good overall representation of the freezing line is³⁰ (with density given in LJ units)

An expression of this form was used by Rosenfeld in 1976⁶¹ and more recently by Khrapak and others.^62,63 Since the freezing line is an approximate isomorph, for the LJ system the functional form Eq. (8) follows from Eq. (6).³⁰ 

Figure 2(a) plots the reduced viscosity of the LJ system as a function of T_F(ρ)/T using Eq. (8) for T_F(ρ). There is excellent collapse. Close to T_F we confirm the linear dependence of $lnη̃$ on T_F(ρ)/T predicted by Kaptay’s liquid-metal constant activation energy generalization of the Andrade equation⁶⁰ (dotted straight line),

Systematic deviations from this are seen moving away from the freezing line, however.

Figure 2(b) investigates the functional form of F in Eq. (7) by plotting the logarithm of the reduced viscosity in a log-log plot with $η̃0$ as a free parameter. Care must be exercised when subjecting data to two logarithms, but keeping this in mind the figure shows that for $η̃0=0.41$ most data collapse onto a line with slope −1/2. This means that the LJ system’s viscosity is well represented by the following expression:

To validate this directly, Fig. 2(c) plots the reduced viscosity data as a function of T/T_F(ρ) with Eq. (10) as the dashed green line. The full curves of Fig. 1(b) are the corresponding non-reduced viscosity predictions. Using instead of 1/2 the exponent 0.4⁵⁵ gives a worse fit; see the supplementary material.

Equation (10) involves just two dimensionless parameters, $η̃0$ and B. The reduced viscosity at freezing $η̃F$ determines B via

The two free parameters of Eq. (10) may be determined from one viscosity measurement at freezing and one measurement away from the freezing line. For the LJ system, we find $η̃F=5.2$ and $η̃0=0.41$⁠, corresponding to B = 2.54.

Figure 3 reproduces data for methane and argon in which Figs. 3(a) and 3(d) are density-temperature phase diagrams showing the state points at which the viscosity was measured. Figures 3(b) and 3(e) compare the reduced viscosity data with best fits to Eq. (10). In both cases, Eq. (10) reproduces the data within 10% (marked by the dashed lines). Finally, Figs. 3(c) and 3(f) compare the non-reduced viscosity data to the prediction along a few isotherms.

Recently, Meyer et al. simulated liquid sodium modeled by a density-dependent effective pair potential derived from the Fiolhais model of electron-ion interaction with self-consistent screening using a local field correction.⁷⁰ As is clear from Fig. 4, the reduced viscosity data of the sodium model conform well to Eq. (10).

Based on a hard-sphere argument, Andrade in 1934 derived an expression for the viscosity of liquids at freezing, which implied a constant reduced viscosity along the freezing line identified by the Lindemann melting criterion.⁷¹ In 2005, Kaptay generalized this expression to η ∝ m^1/2T^1/2V^−2/3 exp(CT_F/T), which is Eq. (9) in non-reduced units.⁶⁰ This expression was shown to work well for 18 metals at ambient pressure.⁶⁰ Equation (10) is qualitatively different by having $(TF/T)1/2$ instead of T_F/T in the exponential. The new viscosity equation represents data well in the region of the phase diagram studied by Kaptay and others, as well as at much higher temperatures.

In summary, Eq. (10) suggests a quasiuniversal viscosity expression for simple liquids,^22,44,72 which reduces the number of independent experiments/simulations needed to map out the viscosity as a function of state point in the phase diagram. Just two measurements of viscosity can now be used to predict its value throughout a fairly large liquid region of the phase diagram. This has potential applications in fields where the pressure can be very high, such as in some lubrication applications and magma geology. Obviously the chemical species there are more complicated than the small molecules considered in the present study. Nevertheless, the scaling treatment verified here may apply for more complex chemical systems and throw light on some important outstanding fundamental questions such as what is the viscosity of the Earth’s outer core and how it varies with depth.^73,74

The argument leading to Eq. (10) was based in part on the isomorph theory which is, however, also consistent with Kaptay’s expression Eq. (9). An important task for the future is to justify the temperature square root needed to fit data properly, which as far as we are aware is not yet derivable from theory. Another important task is to check Eq. (10) for more real liquids and computer-simulated simple liquids.

See the supplementary material for details on the computer simulations and comparisons of our findings to the viscosity expressions of van Loef, Longuet, and Rosenfeld, as well as showing that it is possible to use a reference isomorph that is not the melting line.

We thank Professor J.-F. Wax for making the numerical data of Ref. 70 available and are grateful to Vadim Brazhkin and Kenneth Harris for helpful correspondence. This work was supported by the VILLUM Foundation’s Grant Nos. VKR-023455 and 16515 (Matter).