Shear Viscosity of Liquid Potassium and Cesium : a Simulation Study

The density and temperature dependences of the shear viscosity of liquid potassium and cesium are studied. The stress autocorrelation function is calculated from equilibrium molecular dynamics simulations. Using the GreenKubo formula, the shear viscosity is obtained. Interionic interactions are calculated by Fiolhais potential and are validated by comparison between simulation and experimental data along the liquid-gas coexistence curve for K and Cs. For both metals, three isochors and one isotherm are investigated. The recently proposed relation in [Phys. Rev. B 93, 214203 (2016)] is tested in the cases of K and Cs and it appears that this function reproduces qualitatively and quantitatively well the behavior of each element.


Introduction
In our recent article [1], the density and temperature dependences of the shear viscosity of liquid sodium were studied.Along isochoric lines, the viscosity presents a minimum, while it monotonically increases along isotherms.An expression was proposed for the viscosity as a function of temperature and density which reproduced our data for liquid sodium at any density in the range [1000; 2000 kg.m −3 ] and any temperature in the range [700; 7000 K].
Sodium is one of the alkali metals and it was reported that several properties such as the diffusion coefficient have an universal behavior among these metals.In this spirit, potassium and cesium are now studied.These metals have a high chemical reactivity especially at high temperatures.As a consequence, there are little or no experimental measurements of transport properties in the literature, even near the triple point.Adding to this difficulty, an increase of the density makes measurements very impractical to do.In this context, prediction of these properties using numerical simulation appears to be an interesting alternative to the experimental determination.
The purpose of this work is threefold: (1) investigate the behavior of the density and temperature dependences of the viscosity for K and Cs, (2) compare them with the results of sodium and (3) test the relation proposed for sodium on its counterparts.

Simulation method
We performed equilibrium molecular dynamics simulations, carried out in the NVE ensemble with N=2048 particles.The velocity Verlet algorithm was used to calculate positions and velocities of each particle.Fiolhais potential [2] has been chosen for our simulations, given that it has proven its capability to describe the alkali metals in the liquid state [3].In addition, the previous study on sodium showed its validity even in the high density region.The effective interionic potential in the metal can be expressed, in atomic units, by where Z is the ionic valency.This potential is density dependent, but does not depend on the temperature.It displays a minimum, followed by Friedel oscillations due to the electron gas screening.The first term in the right side of Eq. 1 corresponds to the direct Coulomb repulsive interaction between two ions while the second one describes electron-ion and electronelectron interactions.The total expression behaves asymptotically as r −3 , so that there is no need to take into account long range interactions and that the potential can be truncated (at the node just below half the box size in our case).Complete analytical expressions, as well as references to earlier studies, can be found in Ref. [1].Viscosity values are obtained by integration of the stress autocorrelation function (SACF) thanks to the Green-Kubo relation: where V is the volume of the system, T , its temperature and k B , the Boltzmann constant.δσ αβ (t) denotes the difference between the instantaneous value of the stress tensor αβ component, σ αβ (t), and its statistical average value <σ αβ > (with α, β = x, y, z) [1].The thermalisation stage lasted at least 60 000 steps during which the velocities were rescaled to the desired temperature every 50 steps.In order to reach the required accuracy when computing the viscosity, production steps are about 6 000 000 for K and Cs, that is a duration of 1 and 2 ns, respectively.During the production stage, we observed that the drift of the temperature was less than 1%.As a consequence, it was not necessary to perform simulation within NVT ensemble in order to control the temperature.The SACF fluctuates strongly and three techniques were used to improve its accuracy.(1) We considered many time origins (at least 590 000 in Eq. 2).( 2) An average was done over three off-diagonal elements of the stress tensor σ xy , σ yz , σ xz , and three directions obtained by 45 degrees rotations of the axes, namely 1 2 (σ xx − σ yy ), 1 2 (σ yy − σ zz ), and 1 2 (σ xx − σ zz ).(3) 8 independent microscopic states corresponding to the same macroscopic state point were considered.Thus, we reach an error of less than 5%.For further details about this point, see Ref. [1] where the method for viscosity calculation is exactly the same and more detailed.

Phase diagram
At ambient conditions, alkali metals are considered as simple because they only have one valence electron and crystallize in the body-centered cubic structure.Under pressure, they transform first into face-centered cubic.On further compression, a wide variety of more or less complex phases has been observed and reviewed in recent years.Degtyareva [4] summarized all the sequences of structural transformations on pressure increase for each  [6].From the melting point, the positive gradient of the melting temperature plateaus at (5.8 GPa; 550 K).
The melting curve decreases up to the bcc-fcc-liquid triple point found at (13.6 GPa; 466 K).For pressures higher than this point, the melting temperature remains constant up to 15.6 GPa.Further pressure increase induces a negative slope up to a clear minimum (19 GPa; 390 K).The experimental critical and melting points of K and Cs are summarized in the Table 1.For all alkali metals, simulations are prevented in the low density region by the metal -non metal transition at about 2ρ c .In order to remain in the metallic region, densities lower than 2.5ρ c are not considered in this study.
In the framework of molecular dynamics simulation, it is more convenient to use temperature T and density ρ.However, published experimental and ab-initio studies do not specify the density values as a function of T and P. Knowledge of the density values is needed because above a given pressure, the physical properties of liquid potassium and cesium are expected to change, namely above 19 GPa [6] and 3.9 GPa [5], respectively.The potential we used to describe the atomic interactions is a priori not able to account for these changes.To our knowledge, the only study which indicates approximately the density at high pressures was made by molecular dynamics simulation on liquid cesium using embedded atom model [8].At 493 K and 3330 kg.m −3 , the authors found a pressure corresponding to 2.98 GPa.This justifies our study of density dependence up to 3600 kg.m −3 .Concerning potassium, density specifications do not exist in the literature to the best of our knowledge.

Comparison with literature
In Figure 1, we compare the calculated shear viscosity with the experimental data of K (panel a) and Cs (panel d) along the liquid-gas coexistence.We can see that there exists a slight shift of the same order between both sets of curves.As outlined in our study of Na, the difference can be explained by the possible inaccuracies in the description of the interactions.However, it is important to note that measurements are done indirectly and that the uncertainties are very large.Moreover, other parameters such as the presence of impurities in the material can influence the viscosity value.Measures date back to the sixties or earlier and to our knowledge, there are no updated data available in the literature.So, these experimental data should be confirmed.Nevertheless, experimental and numerical studies show similar results  and the same trend is observed along this line.This constitutes the required condition for discussing the density and temperature dependences of the viscosity in the following.

Temperature dependence
Panels c and f of Figure 1 show the temperature dependence of the viscosity of K and Cs respectively.For both, the viscosity is studied along three isochors: at the density of the melting point ρ m , at a lower density corresponding to 0.75ρ m and at a higher density equal to 1.5ρ m .For densities equal or higher than ρ m , we obtain the same qualitative behavior of the viscosity.As temperature increases, viscosity decreases quickly at low temperatures, reaches a minimum, increases slowly in the high-temperature region.At low densities, the viscosity behavior of these metals can be modeled by a simple linear law.For both density ranges, the same typical behaviors were observed in the case of Na.Each isochoric line was fitted by the empiric relation proposed in the study of Na [1]: As we can see, a very good agreement between the data and the fit over the whole temperature interval is found.Thus, in addition to the results of Na, this expression is also able to reproduce qualitatively and quantitatively the temperature dependence of the viscosity for K and Cs over the whole temperature range studied.

Density dependence
In panels b and e of Figure 1, shear viscosity is represented as a function of density for K and Cs.The density dependence is studied along an isotherm corresponding approximately to 0.64T c (at 1430 K for potassium and 1230 K for cesium).For both elements, an increase of the density implies a higher viscosity.At high density, the behavior of the viscosity is linear while in the low density region, viscosity values are not accessible due to the metal-non metal transition.Nevertheless, at intermediate densities, the behavior does not follow a linear law anymore.According to Chapman-Enskog theory [17], the viscosity must remain positive as ρ tends to zero; this would not be recovered by extrapolating our high density results.In our study of Na, we obtained a relation from Eq. 3 which accounts for the density dependence as where a, b, c, and d were obtained from the density dependence of A, B, and C [1].In Figure 1, the isotherms are directly fitted using Eq. 4. At high density, the above mentioned linear behavior is recovered.For K and Cs, a good agreement between our simulations and Eq. 4 is observed.As in the case of Na, we are able to model the density dependence of the viscosity.Consequently, we believe that this relation can be applied to the different alkali metals.

Conclusion
We have presented a systematic investigation of shear viscosity for K and Cs, obtained by equilibrium MD simulations.We confirmed that the relation proposed previously (Eq.3) for Na is able to reproduce quantitatively the behavior of the viscosity for K and Cs over the wide temperature and density ranges studied.For both elements, the behavior of the temperature and density dependences is in qualitative agreement with those of Na.Investigations about the universality of viscosity behavior among the alkali metals are underway.

Figure 1 .
Figure 1.Shear viscosity along the liquid-gas coexistence line (panel a for K, panel d for Cs); along an isotherm line (panel b for K, panel e for Cs); along isochoric lines (panel c for K, panel f for Cs).Experimental data correspond to Refs [9] -[16].

Table 1 .
[5]erimental critical (T c , ρ c ) and melting (T f , ρ f ) points values[7]of K and Cs.In the work of Falconi et al.[5], the phase diagram of cesium is presented in the (T ,P) plane.The Cs melting curve shows drastic changes in slope with pressure.Starting from ambient pressure, a positive slope is observed then two maxima at (2.25 GPa; 470 K) and (3.05 GPa; 471 K).Further increase of the pressure induces a strongly negative slope of the melting curve until (4.2 GPa; 361 K) before again becoming strongly positive.study ends at the highest pressure investigated by experimental measurements, equivalent to 9.8 GPa.The phase diagram of K is drawn up to 22 GPa by Narygina et al.

Table 2 .
Parameters values of Eq. 3 for K and Cs.