Viscosity bounds in liquids with different structure and bonding types

Recently, it was realised that liquid viscosity has a lower bound which is nearly constant for all liquids and is governed by fundamental physical constants. This was supported by experimental data in noble and molecular liquids. Here, we perform large-scale molecular dynamics simulations to ascertain this bound in two other important liquid types: the ionic molten salt system LiF and metallic Pb. We find that these ionic and metallic systems similarly have lower viscosity bounds corresponding to the minimum of kinematic viscosity of about 10$^{-7}$ $\frac{{\rm m}^2}{\rm s}$. We show that this agrees with experimental data in other systems with different structures and bonding types, including noble, molecular, metallic and covalent liquids. This expands the universality of viscosity bounds into the main system types known.


INTRODUCTION
Different liquids are considered for use as coolants in nuclear reactors including molten salts and metals [1][2][3][4].In these applications, viscosity η and thermal conductivity are two important properties characterising the performance of these liquids and governing their flow/diffusion and thermal transports properties.Understanding and predicting these properties over a range of temperatures and pressures is therefore important from the application point of view.This understanding is also of fundamental theoretical importance, in view that properties of liquids are strongly system-dependent and hence are considered to be not amenable to a general theory, in contrast to solids and gases [5].
Viscosity of low-temperature dense liquids is governed by the activation energy barrier U as: where U is set by the inter-molecular interactions and structure in the liquid, T is temperature and η 0 is the prefactor related to the high-temperature limit of η.Here and below, k B = 1.Eq. ( 1) applies in the low-temperature liquid-like regime where a molecule oscillates before undergoing a jump into a neighbouring quasi-equilibrium position [6,7].In this regime, η increases on cooling and, if crystallization is avoided, can vary by up to 15 orders magnitude in the viscous regime.At high temperature, the oscillatory component of particle motion is lost.This can occur in either the gas phase at pressures below the critical pressure or above the Frenkel line in the supercritical regime where particle dynamics become purely diffusive [8].In the gas-like regime of particle dynamics, η follows a different temperature dependence: where ρ is density, v is the average particle speed and L is the particle mean free path, and where η increases with temperature as η ∝ T . The increase of η at low temperature and its increase at high temperature implies that η has a minimum.It turns out that the value of the minimum of the kinematic viscosity ν = η ρ , ν m , can be approximately evaluated in terms of fundamental physical constants as [9,10]: where m is the mass of the molecule and m e is the electron mass.Deriving Eq. ( 3) involves two steps.First, it was shown that the minimum of ν depends only on two parameters characterising a condensed matter system: interatomic separation a and the largest, Debye, vibration frequency ω D .Second, characteristic values of a and ω D are set by the Bohr radius and Rydberg energy.This gives Eq. (3) [9,10].
Eq. ( 3) has explained the long-standing question considered by Purcell and Weisskopf, namely why the viscosity of all liquids never falls below a certain value comparable to the viscosity of water at room conditions?[11].The answer comes in two parts [12].First, viscosities stop decreasing because they have minima.Second, those minima are largely fixed by fundamental constants in Eq.
In addition to liquids, viscosity minima have been of interest in other areas of physics, including holographic models based on the correspondence between stronglyinteracting field and gravity theories [13].More generally, understanding the origin of bounds on system properties has enthralled physicists, including those interested in collective dynamics and systems where many interacting agents operate.Apart from the interest in the values and origins of the bounds themselves, there is another reason why these bounds are important: finding and understanding bounds on physical properties often means that we enhance our grasp of, or clarify the underlying physics of, the property in question [10,14,15].
In addition to viscosity, it was realised that Eq. (3) provides a lower bound for an unrelated property of liquids: thermal diffusivity α [16].As discussed earlier, this property is directly relevant to industrial properties where molten salts are used including the operation of coolants in nuclear reactors where the heat transfer processes are important.
The lower bound ν m in Eq. ( 3) is of the order of 10 −7 m 2 s (this value corresponds to m equal to the proton mass which sets the magnitude of atomic masses).This agrees with experimental viscosity minima for noble, molecular and network liquids to within a factor of 1-3 [9].It was also found to agree with the experimental high-temperature limiting value of viscosity of metallic alloys [15,17].This analysis was further extended in Ref. [18].
However, no estimations of viscosity minima in molten salts have been undertaken experimentally due to high melting points and hence very high temperatures required to reach the minima.It therefore remains unknown how the viscosities of molten salts conform to the presumably universal crossover between the liquidlike and gas-like regimes and to the theoretical minimum (3).Apart from theoretical importance, knowledge of this would be important from the application point of view: knowing the pressure and temperature conditions of the minima of kinematic viscosity and thermal diffusivity would enable predictions of the optimal state of operation of molten salts.
Here, we perform large-scale molecular dynamics (MD) simulations of the molten salt LiF as a case study (LiF is a common component in molten salt mixtures used in nuclear reactors [2]).We also simulate metallic Pb.We find that these ionic and metallic systems have lower viscosity bounds corresponding to the minimum of kinematic viscosity of about 10 −7 m 2 s .We show that this agrees with the experimental data in other systems with different structures and bonding types, including noble, molecular, metallic and covalent liquids.This expands the universality of viscosity bounds into the main types of systems known.

METHODS
We use the DL POLY MD package [19].For LiF, we used the empirical potential as in the previous work [20][21][22] with parameters derived earlier [23].For Pb, we used the potential from Ref. [24].We also simulated liquid Ar using the standard Lennard-Jones potential.
Unless otherwise stated we simulated for 1,000,000 time steps, with a fixed simulation time step of 0.001 ps.The system size used varied from 2000 atoms in LiF to 5120 in Pb.We have also simulated larger systems with 100,000 particles and found that the values of viscosity collected were consistent regardless of the system size.Similar behaviour of measured viscosities with system size has been found previously [25].We simulated a wide range of temperatures for each pressure.We first equilibrated the system at each pressure in the constantpressure ensemble for 20,000 steps and then performed production runs for 1,000,000 steps in the constant energy and volume ensemble where the data were collected for calculating properties.
The dynamic viscosity η was calculated using the Green-Kubo method [26,27] as where V is the volume of the system and P xy is the xy component of the stress tensor.
Obtaining accurate statistics via the Green-Kubo method is a well known computational issue [28], which we address by averaging statistics over 20 independent initial conditions.This has been sufficient in prior work involving viscosity calculations [25].
We also calculate the kinematic viscosity ν = η ρ governing the non-equilibrium flow and other properties such as, for example, fuel atomization quality [29].Density ρ was calculated at the same state points as η.

RESULTS AND DISCUSSION
We show the calculated η in Figure 1a for two different systems sizes.We observe a good agreement of η calculated in the two systems.We also observe the agreement with the earlier MD results in the low temperature range using the same potential [22].Little experimental data for LiF viscosity exists above 1000 K.We use the experimental data from Ref. [30] which contains the widest range of data and was found to be in agreement with other experiments in the low-temperature range [31,32].The MD results underestimate the experimental viscosity by a factor of about 1.3-5, however, the overall shape of η is similar in experimental and MD data.As discussed in Ref. [22], the simulated viscosity is lower than experimental viscosity due to the model not accurately matching the melting point, around 300 K lower than ex-perimental values (predicting melting points accurately is a more general problem in atomistic simulations due to several factors including short simulation times which can be insufficient to exceed slow kinetics of phase transformations, small system sizes compared to experimental ones and hence the absence of long-wavelength fluctuations with large amplitudes which are efficient in destabilising the solid phase and so on).Translating all temperatures upwards by this melting point discrepancy of 300 K while maintaining all values of viscosity results in a much better agreement, demonstrating that the empirical potential accurately model the dissipative dynamics of the liquid.We note that the underestimation of the melting point does not affect our results since we are interested in the values of viscosity minima.
We observe that the calculated η tends to an approximately constant value of about η = 2 × 10 −4 Pa•s at high temperature.This is close to viscosity minima in noble, molecular and network liquids [9].Simulating higher temperature results in the known instability of the Born-Mayer potential at short distances.This could be fixed by, for example, adding a short-range repulsive terms to the potential (e.g. in the form of Ziegler-Biersack-Littmark potentials [34]), however, this is not required for the purposes of the current work aimed at evaluating the limiting lower viscosity bounds.
Density as a function of temperature is shown in Fig. 1b for different pressures.We observe that the slope of density decrease becomes smaller at high pressure, corresponding to the decrease of the thermal expansion coefficient with pressure.This is consistent with the common behavior seen in solids [35].
The calculated kinematic viscosity ν = η ρ is shown in Figure 1c for different pressures.We make several observations.First, ν reaches a constant value at high temperature and shallow minima at low pressures.The minima become more apparent here, compared to in the dynamic viscosity in Figure 1a, due to the decrease of density at high temperature.Second, pressure increases the values of the minima and constant values at which ν tends to a constant at high temperature.This is similar to what is observed in noble, molecular and network liquids [9] and is related to the increase of the activation energy for molecular rearrangement with pressure and associated increase of viscosity.Third, the value of ν at their constant value or minima is in the range (2-5)•10 −7 m 2 s .This is of the same order, 10 −7 m 2 s as predicted by Eq. ( 3).This last point is important and implies that the viscosity of ionic systems and molten salts have the same behavior in terms of the value of their viscosity minima.We will make the comparison between calculated and theoretical values more quantitative later on, alongside the other liquids we study.
η and ν for the simulated metallic system, Pb, are shown in Fig. 2.These show very similar results to those of LiF.η increases with pressure and tends to about η = 2×10 −4 Pa•s at high temperature.ν has a minimum FIG.1: (a) Dynamic viscosity of LiF for small, 500 atoms, and large, 100,000 atoms, systems.We compare to earlier MD results [22] as well as experimental results in [30], which has been shown to have good agreement with more recent experiments.[33], [32]  s at low pressure as predicted by Eq. ( 3).To emphasise this similarity, we plot ν for Pb, LiF and Ar in Fig. 2c.
To compare our results for ionic LiF and metallic Pb with a wider data set, we add the experimental η and ν for noble Ar and molecular CH 4 and CO 2 as well as network H 2 O [36] to the plots in Fig. 3.It is appropriate to note here that, differently from noble and molecular liquids, the NIST database [36] does not include viscosities of ionic, molten salts and metallic systems.Part of the issue is the experimental difficulty related to high melting and boiling points of these systems.The value of our current simulations therefore includes provision of this data and serving as a guide for future high-temperature experiments.
The results collected in Fig. 3 are consistent with what is expected for the dependence of viscosity on temperature, as discussed in the Introduction: the decrease with temperature in the liquid-like regime, followed by its increase in the gas-like regime.We observe a great amount of variation in the values of viscosity collected for different systems, as well as the variability of the shape of viscosity curves.We also observe a consistency in the magnitude of the minimum values of viscosity, in agreement with Eq. ( 3) predicting close values of the lower viscosity bound in all liquids.This includes the ionic molten salt LiF and metallic Pb.
We now address the quantitative comparison between the predicted and calculated ν m .In Table I we show theoretical and calculated ν m for LiF, Pb and Ar.The ratio between calculated and theoretical values is in the range 1.7-4.6.This is similar to the range of ratios for the large set of noble, molecular and network liquids found earlier, where this ratio is in the range 0.5-3 [9].The difference between theoretical and observed ν m by a factor of 1-4 is related to a number of approximations involved in deriving Eq. ( 3).This includes approximating the interatomic separation by the Bohr radius and the bonding energy by the Rydberg energy.Recall that the main purpose of Eq. ( 3) is two-fold: first, it shows that viscosity minimum is largely the same for all liquids because it is set by fundamental physical constants.Second, Eq. ( 3) serves to evaluate a characteristic value, order of magnitude, of ν m .
We have discussed viscosity minima in ionic, metallic, noble, molecular and network liquids.This leaves the remaining important bonding type: covalent.Silica is a commonly discussed system in which covalency is strong: it is a system with mixed bonding where covalency and ionicity make similar contributions to the bonding type.Experimentally, viscosity of silica was measured up to about 3000 K, and yet no minimum was seen due to challenges involved in high-temperature experiments [38].Simulated silica, taken to yet higher temperature in excess of 6000 K, shows saturation to a constant value of about 10 −3 Pa•s [39].This corresponds to ν m of about 5•10 −7 m 2 s and falls in the range of ν m predicted theoretically by Eq. ( 3).
An interesting observation from Eq. ( 3) is that the minimal viscosity is a quantum property and approaches zero in the classical limit ℏ = 0.This might be perceived to be at odds with our thinking about liquids as mostly high-temperature classical systems.However, the nature and origin of interatomic forces and radii are quantum-mechanical, and Eq. ( 3) reminds us of this [12].This then brings the question of how our classical MD simulations reproduce the lower bound of liquid viscosity, the essentially quantum effect?We understand this if we recall that the parameters of empirical potentials in classical MD simulations are tuned to result in interatomic separations and energy values in real systems (either by fitting to experiments or quantum-mechanical simulations) where they are quantum-mechanical in origin.The classical simulations capture quantum effects, including the lower viscosity bound, through these potential parameters.

SUMMARY
We have explored the nature of viscosity minima in ionic molten salt liquid LiF, complemented by metallic Pb.We have found that these systems have lower viscosity bounds corresponding to the minimum of kinematic viscosity of about 10 −7 m 2 s .This agrees with the experimental data for other systems with different structures and bonding type, including noble, molecular, metallic and covalent liquids, and it expands the universality of viscosity bounds into the main types of systems known.
In future work, it may be interesting to develop more accurate potentials or employ quantum-mechanical molecular dynamics simulations to simulate molten salts.
We have previously found that Eq. ( 3) also give the minima of an unrelated, yet important property: thermal diffusivity [16].This was supported by the experimental data for noble and molecular systems.Our current findings therefore suggest that, similarly to the kinematic viscosity, the minima of thermal diffusivity are more universal and include other system types including the ionic molten salts.This is interesting to explore in future work.
We are grateful to V. V. Brazhkin for discussions, EP-SRC (grants No. EP/X011607/1 and EP/W029006/1) and Queen Mary University of London for support.Computing resources for the 10 5 LiF data were provided by STFC Scientific Computing Department's SCARF cluster.
For all other simulations conducted, this research utilised Queen Mary's Apocrita HPC facility, supported by QMUL Research-IT.http://doi.org/10.5281/zenodo.438045.The U.S. researchers acknowledge the financial support from the U.S. Department of Energy via Grant no.DE-NE0009288.
FIG.1:(a) Dynamic viscosity of LiF for small, 500 atoms, and large, 100,000 atoms, systems.We compare to earlier MD results[22] as well as experimental results in[30], which has been shown to have good agreement with more recent experiments.[33],[32].(b) Density at different pressures.(c) Kinematic viscosity of LiF simulated using 2000 atoms for a range of pressures.

TABLE I :
Comparison of predicted values calculated from the formula against simulated or experimental values of the kinematic viscosity minimum.