Comment on ‘Pseudo hard-sphere viscosities from equilibrium molecular dynamics’

In a recent article, Nicasio-Collazo et al (2023 J. Phys.: Condens. Matter 35 425401) explore the viscosity of the pseudo-hard-sphere (PHS) model. In this comment, we highlight some discrepancies with expected behavior, and compare their results to new simulations of the same model as well as to true hard spheres. In contrast to the results of Nicasio-Collazo et al, our results follow the relation between shear, bulk, and longitudinal viscosity expected for isotropic fluids. Moreover, we observe clear differences in behavior between PHS and true hard sphere, and encourage future hard-sphere studies to focus on the true hard sphere model whenever possible.

In a recent article [1], Nicasio-Collazo et al explore the viscosity of the pseudo-hard-sphere (PHS) model.The authors observe shear viscosities that are significantly higher than literature data for true hard spheres.Additionally, they find that the expected relation Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. between the longitudinal viscosity η l , shear viscosity η and bulk viscosity η b is significantly violated for the investigated PHS system.This is surprising, as equation ( 1) should hold for isotropic fluids [2], which includes both PHS and true hard spheres.
Here, we check their results against both true hard spheres and PHS and find significantly different results.In particular, we simulate true hard spheres using event-driven molecular dynamics simulations using the approach of Alder et al [3].The code used is an adaptation of the one in [4], modified to periodically output the stress tensor integrated over a short time interval (using equation (9) in [3]).Simulations were run for 7.5 × 10 4 τ after equilibration, with τ = √ βmσ 2 the time unit.Here, β = 1/k B T the inverse thermal energy, k B is Boltzmann's constant, T is the temperature, m is the mass of a single particle, and σ is the particle diameter.For the PHS system, we used the LAMMPS simulation package [8].Note that the PHS pair interaction is given by [9] U PHS (r) = { 50 ( 50 49 where r is the distance between two particles, ε the interaction strength, and σ the effective diameter.As was done in [1], we set the interaction strength such that k B T/ϵ = 1.4737.After equilibration, the thermostat was switched off, velocities were reinitialized, and the simulation was continued in the microcanonical ensemble, running for 10 4 τ with an integration step size of 2 × 10 −4 τ .As a check, we have confirmed that the simulations accurately reproduce the deviations from the Carnahan-Starling equation of state reported in [9].For both methods, each simulation contained N = 15 3 = 3375 particles, and results were averaged over ten independent simulations. In figure 1, we plot the measured viscosities for both PHS and true hard spheres as a function of the packing fraction ϕ, normalized using the Enskog expression (using the expressions in [5]).The top panel shows the shear viscosity η.The hard-sphere data (in red) is compared with the data from [5,6] for a similar system size, finding excellent agreement.The corresponding viscosity for the PHS model (in blue) is slightly but measurably lower, and agrees well with the data by Pousaneh et al [7] for ϕ ≳ 0.2.As already noted in [1], the data from [7] appears to underestimate the viscosity at lower packing fractions.Surprisingly, the results of Nicasio-Collazo et al agree more closely with true hard spheres than with PHS, in contrast to what is shown in [1].This discrepancy likely stems from an error in unit conversion in [1]: Nicasio-Collazo et al state that they use τ ϵ = √ mσ 2 /ϵ as their unit of time, in contrast to the unit τ = √ βmσ 2 normally used for studies of true hard spheres.However, in the figures in [1], the literature values for hard spheres do not appear to have been converted to this new set of units, and as a result the reported hard-sphere results appear too low by a factor √ k B T/ϵ ≃ 1.21.Figure 1(b) shows the bulk viscosity η b .For true hard spheres, we again find good agreement with the work of Sigurgeirsson and Heyes [5] if we compare to their system size of N = 2048 particles.For N = 4000 particles, their results deviate significantly from the smaller system sizes, as well as from our results.This discrepancy does not appear to be simply a finite-size effect: repeating our hard-sphere simulations at these densities using N = 8000 yields essentially the same viscosities.For PHS, we find consistently lower bulk viscosities than for true hard spheres.This conflicts with the measurements of Nicasio-Collazo et al, which instead showed higher bulk viscosities for the PHS system at high densities.Viscosities of the hard-sphere (red) and pseudo-hard-sphere fluids (blue), as a function of the packing fraction ϕ, compared to the data in [1] (black).Shown are the shear viscosity (a), bulk viscosity (b), and longitudinal viscosity (c), each normalized by their respective Enskog prediction.Error bars are based on the standard error of ten independent runs.Data are compared to literature values by Sigurgeirsson and Heyes [5] and Pieprzyk et al [6] for hard spheres, and Pousaneh and de Wijn [7] for PHS.
In figure 1(c), we plot the longitudinal viscosity, calculated both directly and via equation (1).For both true hard spheres and PHS we obtain indistinguishable results using the two methods, as expected for an isotropic fluid.In contrast, [1] found significant differences between the results from the two methods, with discrepancies approaching a factor of two at some densities.Note that Nicasio-Collazo et al [1] propose fitting the longitudinal viscosity with a function that vanishes at zero density.However, since the bulk viscosity vanishes in the zero-density limit [5], equation (1) shows that the zerodensity limit of η L is given by 4η 0 /3 (with η 0 the ideal-gas shear viscosity).
In short, an independent recalculation of the viscosities of the model studied in [1] yields qualitatively different results.Most importantly, the recalculated data satisfies equation ( 1), as expected for an isotropic fluid.Additionally, it may be good to highlight the fact that all quantities measured in [1] are also readily accessible in a true hard-sphere model system.This includes the kinetic, potential, and cross-components of the viscosities (see e.g.[5]).Moreover, in particular the bulk viscosity of PHS appears to deviate significantly from that of true hard spheres, showing that in this respect, PHS are not a great model for hard spheres.Hence, for studying the behavior of hard-sphere-like systems, the high efficiency and simplicity of true hard spheres implies that they should be the preferred choice whenever possible.
Data files containing the measured autocorrelation functions for the PHS system, as well as tables of the measured viscosities, are provided as supplemental data.

Figure 1 .
Figure 1.Viscosities of the hard-sphere (red) and pseudo-hard-sphere fluids (blue), as a function of the packing fraction ϕ, compared to the data in [1] (black).Shown are the shear viscosity (a), bulk viscosity (b), and longitudinal viscosity (c), each normalized by their respective Enskog prediction.Error bars are based on the standard error of ten independent runs.Data are compared to literature values by Sigurgeirsson and Heyes[5] and Pieprzyk et al[6] for hard spheres, and Pousaneh and de Wijn[7] for PHS.