Practical formula for the shear viscosity of Yukawa fluids

A simple practical formula for the shear viscosity coefficient of Yukawa fluids is presented. This formula allows estimation of the shear viscosity in a very extended range of temperatures, from the melting point to $\simeq 100$ times the melting temperature. It demonstrates reasonable agreement with the available results from molecular dynamics simulations. Some aspects of the temperature dependence of the shear viscosity and diffusion coefficients on approaching the fluid-solid phase transition are discussed.


I. INTRODUCTION
Studies of static and dynamical properties of Yukawa systems constitute an important interdisciplinary topic with applications to strongly coupled plasmas, complex (dusty) plasmas, and colloidal suspensions. Significant efforts have been made over the years to understand the transport properties of such systems. In particular, this includes diffusion, viscosity, and thermal conductivity. Extensive numerical simulations have been performed and a large amount of accurate data for the transport coefficients exist. What is often required in practical situations is a simple and accurate tool to estimate the transport coefficients in a broad range of parameters.
The main purpose of this article is to present a useful practical expression to estimate the shear viscosity coefficient of three-dimensional Yukawa fluids in a wide regime of coupling and screening. The proposed formula is shown to describe quite well (deviations are within ±10%) the available results from numerical simulations in a very broad temperature range, from the melting temperature to 100 times the melting temperature. It can be particularly useful in the context of complex (dusty) plasmas and related soft weakly dissipative systems.

II. BACKGROUND INFORMATION
The Yukawa systems are characterized by the repulsion between point-like charged particles immersed into neutralizing background, which provides screening. The corresponding interaction potential (also known as the Debye-Hückel or screened Coulomb potential) is φ(r) = (Q 2 /r) exp(−r/λ), (1) where Q is the particle charge and λ is the screening length. In the limit λ → ∞ (i.e. screening is absent), the pure Coulomb interaction potential is recovered, corresponding to the one-component plasma (OCP) limit. 1,2 Yukawa potential is widely used as a reasonable first approximation for actual interactions in three-dimensional isotropic complex plasmas and colloidal suspensions. [3][4][5][6][7][8][9] Yukawa systems are conventionally characterized by the two dimensionless parameters: 10 the coupling parameter Γ = Q 2 /aT , and the screening parameter κ = a/λ, where a = (4πn/3) −1/3 is the Wigner-Seitz radius and n is the particle number density. The screening parameter κ determines the softness of the interparticle interaction. It varies from the extremely soft and long-ranged Coulomb interaction at κ → 0 to the hard-sphere-like interaction limit at κ → ∞. In the context of complex plasmas and colloidal suspensions the relatively "soft" regime, κ ∼ O(1), is of particular interest. Thermodynamics and dynamics of 3D Yukawa fluids and crystals in this regime have been extensively studied in the literature. [10][11][12][13][14][15][16][17][18] Depending on the values of Γ and κ, Yukawa systems can form either a fluid or, at sufficiently high Γ, a solid phase. There is no gas-liquid phase transition, because attraction is absent. In the solid phase the body-centered-cubic (bcc) or face-centered-cubic (fcc) lattices can be thermodynamically stable (bcc lattice is thermodynamically favorable at weak screening). The values of the coupling parameter at which the fluid-solid phase transition occurs are usually denoted Γ m . The dependence Γ m (κ) obtained from molecular dynamics (MD) simulations is available. 10,11 A simple formula 19,20 Γ m (κ) 172 exp(ακ) is in a rather good agreement with the tabulated data in the regime κ 5. Here the constant α = (4π/3) of a glass phase on the phase diagram of Yukawa systems have been investigated. 21,22 This regime is not considered below, the attention is focused on the strongly coupled fluid phase.

A. Numerical data
Shear viscosity of single-component 3D Yukawa fluids has been extensively studied using MD simulations, 23-30 see in particular an overview in Ref. 31. Despite the apparent simplicity of the single-component Yukawa model, accurate determination of the shear viscosity by MD simulation is not trivial. This can be exemplified by the significant discrepancies in the results obtained over the years by different authors. 31 In the following, the two data sets tabulated in Refs. 28 and 31 are chosen as the most accurate results presently available.
The original data from Refs. 28 and 31 used in this work are plotted in Fig. 1. Here the conventional (for plasma physics) normalization is used: The reduced shear viscosity coefficient is where m is the particle mass and ω p = 4πQ 2 n/m is the plasma frequency. In Fig. 1, η * is plotted versus the coupling parameter; symbols are the tabulated MD data points 28,31 and the curves correspond to a practical interpolation formula derived later in this work.

B. Normalization
Historically, the use of the plasma frequency ω p in presenting reduced transport coefficients of OCP and Yukawa systems [like in Eq. (3)] originates from the pioneering works of Hansen and collaborators. 23,32 More general macroscopic reduction parameters, not limited to plasma physics context, have been suggested by Rosenfeld. 33,34 Namely, the mean interparticle separation ∆ = n −1/3 and the thermal velocity v T = T /m have been used as the units of length and velocity. In these units the reduced shear viscosity coefficient becomes Reduced coefficient of shear viscosity ηR of Yukawa fluids versus the reduced coupling parameter Γ/Γm. Symbols correspond to the data tabulated in Refs. 28 and 31 (the notation is the same as in Fig. 1). In (a) the entire range of coupling considered, 0.01 ≤ Γ/Γm ≤ 1, is shown. In addition, the interpolation formula (7) is plotted as a blue solid curve. The dashed curves mark ±10% deviation region. The red curve is the modified OCP fit of Eq. This form of the reduced viscosity coefficient is suggested by an elementary kinetic theory formula for viscosity, η ∼ nm v T , for a dense medium of particles with thermal velocities v T and a mean free path between collisions , which is of the order of the average interparticle distance. 34 The relation between the two reduced viscosity coefficients is η R = η * (4π/3) −5/6 √ 4πΓ. Rosenfeld used this normalization first to suggest a quasi-universal excess-entropy scaling of the properly reduced atomic transport coefficients of simple fluids. 33,35 Then he came to a conclusion that for sufficiently soft interaction potentials the excess-entropy scaling should result in a quasi-universal freezing-temperature scaling of the reduced transport coefficients (simply because the excess entropy itself scales quasi-universally with the temperature reduced by its value at freezing). 34,36 Various other arguments have been put forward in favor of the freezing-temperature scaling of the shear viscosity. 25,[37][38][39][40][41][42] Perhaps the most solid basis behind the freezing-temperature scaling is the isomorph theory developed recently. 43,44 With isomorphs defined as the lines of constant excess entropy in the thermodynamic phase diagram, the reduced viscosity is constant along an isomorph because the properly reduced dynamics is. 42 The freezing curve is an approximate liquid-state isomorph. Parallel curves, characterized by fixed ratios T m /T ≡ Γ/Γ m should be approximate isomorphs, too. Application of the isomorph theory to Yukawa systems has been recently discussed in detail. 45 Having all this in mind, the original data from Refs. 28 and 31 plotted in Fig. 1 have been re-scaled to produce the dependence η R on Γ/Γ m . The results are shown in Fig. 2 (a). Here the data corresponding to the fluid regime with 0.01 ≤ Γ/Γ m ≤ 1 are plotted (a few points corresponding to the very weakly coupled gaseous regime were not retained for further analysis). In contrast to Fig. 1, the data for η R versus Γ/Γ m appear to collapse onto a single universal curve. While some scattering of the data points is still present, no clear systematic dependence on the screening parameter κ is evident. The scattering is probably related to an uncertainty in evaluating shear viscosity from MD simulations. The quality of the collapse improves on approaching the fluid-solid (freezing) phase transition. Below we propose a practical expression of the form η R = F(x), where x = Γ/Γ m = T m /T , which provides a good estimate of the shear viscosity of Yukawa fluids in a wide region of the phase diagram, from T m to 100T m .
C. Practical formula Figure 2 (b) and (c) present closer look on the weaker (Γ/Γ m 0.06) and stronger (0.1 Γ/Γ m 1) coupling portions of the data, respectively. Two different functions F can fit these data sets separately. The first fit, is shown in Figure 2 (b). This is just an ad hoc approach for the intermediate coupling regime. There is no physical background behind this choice; the particular functional form has been chosen because of its simplicity. As a consequence, this fit does not (and is not expected to) reproduce the weak coupling asymptote 31 η R ∝ 1/(Γ 2 Λ) at Γ 1, where Λ is the Coulomb logarithm.
More attention should be given to the strongly coupled regime not too far from the fluid-solid phase transition, where some theoretical predictions are available. Here, although the functional form F(x) can be system-dependent, certain universality has been previously suggested. For example, Vaulina et al. 37,38 using the activation energy ideas suggested that the shear viscosity of Yukawa fluids scales as η R ∝ exp(bx) on approaching freezing. They proposed to use this dependence in the regime x 0.5, where they found b 2.9. Similar scaling was proposed by Kaptay 39 for the viscosity of pure liquid metals. He arrived at this scaling combining the Andrade's equation 46,47 with either the activation energy concept or the free volume arguments. Testing this scaling on the experimental data for 15 selected liquid metals the adjustable parameter b was obtained as b = 2.3 ± 0.2. In a recent study by Costigliola et al. 42 dealing with computer simulations of viscosity in the Lennard-Jones liquid and experimental data for argon and methane this scaling has been confirmed only in the close vicinity of the melting temperature. Systematic deviations have been observed when moving away from the melting point and a different scaling of the form F(x) ∝ exp(B √ x) has been put forward. Both expressions are consistent with the isomorph theory (freezing temperature scaling) 42 and we have a good opportunity to test which of the scalings performs better in the special case of Yukawa fluids. We therefore fitted the data corresponding to the strongly coupled regime 0.1 Γ/Γ m 1 using the two functional forms. The results are shown in Fig. 2 (c). The simple exponential scaling with b = 2.28 (dashed curve) does a good job, but only at Γ/Γ m 0.3. The Costigliola's formula with B = 3.64 (solid curve) allows to describe well the data in the entire range considered. Thus, this latter scaling is superior, at least for Yukawa fluids. The explicit expression for the reduced shear viscosity of strongly coupled Yukawa fluids is We can now combine the expressions for F 1 and F 2 in the following form, which then agrees with the individual expressions in their corresponding regimes of applicability and appropriately interpolates between them. The interpolation formula (7) is shown in Fig. 2(a) as the blue solid curve. The adjustable parameter γ is set to γ = 4. The viscosity is predicted correctly within 10% tolerance as indicated by the dashed curves (only 3 from about 90 data points are clearly outside the region marked by the dashed curves). Equation (7) can also be easily rewritten in the form η * (Γ). The corresponding curves are shown in Fig. 1.

D. Alternative formula
An alternative expression for the shear viscosity of Yukawa fluids can also be elaborated. Bastea proposed an accurate three-term fit to describe the behavior of the OCP fluid viscosity η * obtained in MD simulation in a wide range of coupling. 48 As we have observed in figure 2(a), the dependence of η R on x is practically insensitive to κ, at least for κ 3. Rewriting the original fit η * (Γ) into the form η R (x) we obtain the following generalization: This fit is shown by the red solid curve in Fig. 2(a). This fit can be particularly useful near the OCP limit at κ 1.

IV. SOME TENDENCIES IN THE STRONG COUPLING REGIME
At the melting point, both Eq. (6) and (8) yield the viscosity coefficient η R 4.8. Quasi-universality in the dependence η R (x) clearly implies quasi-universality of η R values at the melting point (x = 1). This can be rewritten as where the subscript "m" stands again for the melting point and C is approximately constant. This coincides with the scaling proposed by Andrade for liquid metals at the melting point. 46,47 However, the constant C is only approximately universal even for simple systems. For soft repulsive Yukawa systems considered here the value of the constant (C 4.8) is smaller than that for Lennard-Jones liquid (C 5.2) and argon (C 5.8), recently reported. 42 Liquid metals also demonstrate somewhat higher C, as illustrated in Tab. I, where the reduced shear viscosity coefficients of some liquid metals at the melting temperature have been evaluated using experimental data summarized by March and Tosi. 49 It was previously reported that the Yukawa viscosity model of liquid metals near melt predicts viscosities that are too low. 50 The Stokes-Einstein (SE) relation between the coefficients of viscosity and diffusion can be written as where D is the diffusion coefficient of a sphere of radius R immersed in a medium characterized by the shear viscosity η. Applying the SE relation to atomistic scales (although this does not always works satisfactory, see Ref. 51 and references therein), we arrive at where the characteristic interparticle separation a plays the role of the sphere radius R. Table I demonstrates that this conditions is satisfied to a reasonable accuracy by liquid metals at the melting temperature. Next, combine Eq. (11) with the de Gennes scaling of the self-diffusion coefficient in atomic liquids, 52 where Ω E is the Einstein frequency (this scaling has been recently verified for single-component Yukawa fluids 53 ). This yields which coincides (to within a numerical coefficient of order unity) with the expression obtained by Andrade using completely different arguments. 46,47 Having derived this expression he argued that "when a solid is melted it still retains in the liquid form sufficient of its crystalline character for the molecules to possess a frequency of vibration which is practically the same as that of a solid form at the melting point". 46 Quantitatively, this means that the Einstein frequency is not expected to change much across the fluid-solid phase transition. Indeed, for weakly screened Yukawa systems the Einstein frequency is only slightly higher in a fluid phase as compared to an ideal crystal, as has been recently shown theoretically 54,55 and documented experimentally (using a strongly coupled dusty plasma). 56 According to the Lindemann melting rule, Ω E ∝ T m /ma 2 at the melting point, and this immediately leads to the scaling of Eq. (9). This route provides an alternative derivation of the Andrade' scaling (9). Another consequence of the SE relation is that the reduced self-diffusion coefficient of Yukawa fluids, D R = Dn 1/3 /v T , is expected to scale as ∝ exp(−3.64 √ x) on approaching melting. This scaling is verified in Fig. 3 using the data tabulated in Refs. 32 and 57, which have been re-scaled to the present dimensionless form. 58 It is observed that the dependence D R exp(−3.64 √ x) describes the data quite well in the extended range of coupling. The value D R 0.03 at freezing is consistent with the values reported for several other simple model fluids at freezing (e.g. OCP, Hertzian, Gaussian-core, and inverse-power-law models). 53,59,60 For the strongly coupled Yukawa fluids the SE relation is of the form D R η R 0.13 and Dη(a/T ) 0.08. The latter value is somewhat smaller than those characterizing liquid metals at the melting temperatures (see Table I).

V. CONCLUDING REMARKS
A simple practical expression for the shear viscosity coefficient of 3D Yukawa fluids has been put forward. The proposed formula is applicable in a wide range of coupling and screening and demonstrates reasonable accuracy. The analyzed numerical results related to shear viscosity support the temperature scaling η ∝ √ T exp(B T m /T ) (at a constant density) on approaching the fluid-solid phase transition. 42 Combining this scaling with the Stokes-Einstein relation allows to reproduce quite well the behavior of the self-diffusion coefficient on approaching freezing. Note that this scaling operates in a very wide temperature (coupling) range. In the nearest vicinity of the melting point other scalings can potentially be more appropriate or accurate.
In the context of complex (dusty) plasmas many experiments and simulations have been dealing with twodimensional (2D) mono-layers of particles. [61][62][63][64][65][66] In this case Yukawa (Debye-Hückel potential) is also considered as a reasonable first approximation for in-plane interactions. It would be interesting, therefore, to elucidate whether the freezing temperature scaling is applicable (there are favorable indications 63 ) and what is the functional form of such scaling for 2D (Yukawa) systems. However, since the physics behind the transport coefficients in 3D and 2D is quite different, this topic merits a separate detailed consideration.