Pseudo hard-sphere viscosities from equilibrium Molecular Dynamics

Transport coefficients like shear, bulk and longitudinal viscosities are sensitive to the intermolecular interaction potential and finite size effects when are numerically determined. For the hard-sphere (HS) fluid, such transport properties are determined almost exclusively with computer simulations. However, their systematic determination and analysis throughout shear stress correlation functions and the Green-Kubo formalism can not be done due to discontinuous nature of the interaction potential. Here, we use the pseudo hard-sphere (PHS) potential to determine pressure correlation functions as a function of volume fraction in order to compute mentioned viscosities. Simulation results are compared to available event-driven molecular dynamics of the HS fluid and also used to propose empirical corrections for the Chapman–Enskog zero density limit of shear viscosity. Moreover, we show that PHS potential is a reliable representation of the HS fluid and can be used to compute transport coefficients. The molecular simulation results of the present work are valuable for further exploration of HS-type fluids or extend the approach to compute transport properties of hard-colloid suspensions.

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. soft matter in general. Although HS model is an idealized representation of the interaction between atoms, molecules or colloids, its simple functional form describes very well structural [8,9], thermodynamic [10,11] and dynamical [12,13] behaviour of liquids. Furthermore, this simple interaction potential allows to determine physical properties for in-and out-equilibrium systems. Then, theoretical and simulation approaches to describe physical properties of HS fluids constitute an active research field [14][15][16][17] due to potential applications in a wide range of systems that could be characterised, as a first approach, with the HS fluid.
From the theoretical point of view, structural and thermodynamic properties of this kind of systems can be determined through the so-called perturbation theory [18][19][20], density functional theory [21,22] and integral equations [8,9,23]. Dynamic properties can be computed with the Chapman-Enskog (CE) theory [24][25][26] in the case of molecular liquids and with the mode coupling theory (MCT) [27][28][29] or the self consistent generalized Langevin equation (SCGLE) [30] for colloidal liquids. On the other hand, from computer simulations, structural and thermodynamic properties are easily determined with Monte Carlo techniques for both molecular and colloidal liquids. However, regarding to dynamical properties, in the molecular case the so-called event-driven molecular dynamics (ED-MD) [31,32] is commonly used, whereas, for colloidal systems the kinetic Monte Carlo [33] technique is usually employed.
Focusing on dynamical properties, it is worth mentioning that some drawbacks can be found in both theoretical and simulation approaches. For instance, in the molecular case, it is well known that the CE approach fails to reproduce transport coefficients at high fluid densities [34,35]. Furthermore, theoretical determination of colloidal dynamical properties strongly depends on the approaches used to determine the structure of the HS fluid, which is the input in both MCT and SCGLE theories. On the side of computer simulations, molecular results determined with ED-MD are commonly used to test and sometimes suggest corrections to theoretical approaches. Nevertheless, special care must be taken to handle periodic boundary conditions in the determination of transport coefficients [36,37]. In the case of colloidal suspensions, Brownian dynamics (BD) simulations require the knowledge of force between particles [38]. However, to the best of our knowledge, we have observed that only qualitative representations of HS interactions have been employed [39,40].
In addition to the technical issues mentioned above, systematic studies on dynamical properties and transport coefficients of HS are limited. In fact, the pioneering work by Alder, Gas and Wainwright [31], remains as the most influential work on this topic and has been revisited by other authors in order to improve results using more accurate simulations [41,42] and to propose corrections to the CE theory [13,[43][44][45]. However, the behaviour of potential, kinetic and cross (kinetic-potential) contributions of transport coefficients such as thermal conductivity or shear viscosity should be studied in more detail.
In a similar line, the acoustic spectroscopy technique is also used to determine longitudinal viscosity. However, it is theoretically predicted that both shear and bulk viscosities are related to longitudinal viscosity [64,65]. Therefore, is common to use such relationship to determine longitudinal viscosity [59,66,67].
In the present work, we use the GK formalism [68][69][70] to determine the shear, bulk and longitudinal viscosities as well as the potential, kinetic and cross contributions of the former viscosity as a function of the volume fraction. To this end, the corresponding time correlation functions of the pressure tensor [71] and its fluctuations are required [72], these quantities can be computed by means of simple molecular dynamics (MD) simulations. At this point it is worth recalling that in order to use the GK approach is fundamental to know the force between molecules, that in the case of HS fluid is unknown, since the corresponding interaction potential is discontinuous.
To survey the issue related to the force between hardspheres, along years, different approaches have been proposed to establish a continuous and differentiable function that represents the HS interaction potential [19,[73][74][75][76]. Although many of those approaches have been thought and designed to numerically reproduce structural or thermodynamic properties of the HS fluid [75], some authors have determined dynamic properties like diffusion [77] or shear viscosity [78] as a function of volume fraction. In the case of the last property, non-negligible differences can be found at high volume fractions when results are compared to ED-MD simulations [35,78]. Furthermore, those studies lack of a systematic analysis about correlation functions and shear viscosity contributions.
A continuous potential has been proposed using a physical criterion based on the extended law of corresponding states [79] and the Mie potential [73] Baéz et al mapped the HS potential to a continuous function [80]. This equivalent function is referred as pseudo hard-sphere (PHS) and has been used to successfully reproduce thermodynamic properties of the HS fluid in a wide range of volume fractions [80]. Here, we use the PHS pair potential to straightforwardly determine the force between HS particles in a standard MD algorithm. With this, we are able to compute pressure tensor contributions which allows to determine correlation functions of interest. Then, these data is employed in the GK formalism to compute viscosities of interest and their contributions as a function of the fluid volume fraction.
This work is organised into five sections as follows: In section 2 we discuss details about the interaction potential used to model HSs, also in this section, time correlation functions needed in the GK formalism to determine shear, bulk and longitudinal viscosities are presented. In addition, computer simulations conditions are discussed there. Results and discussion about shear stress autocorrelation and its contributions are detailed discussed in section 3. Then, in section 4 shear, bulk and longitudinal viscosities results of the HS fluid are summarised. Chapman-Enskog corrections and their predictions are discussed in detail within this section. Finally, in section 5, we provide some concluding remarks.

PHS potential
As it was previously mentioned, the HS pair potential can be mapped to a continuous function referred as PHS, which is given by [80] where with r = |r i − r j | being the separation distance between ith and jth particle of mass m and diameter σ, while, ϵ is the interaction strength between them and λ = 50. For a three dimensional, mono-disperse fluid, the dimensionless temperature T * ≡ kT/ϵ is set to T * = 1.4737, in order to accurate reproduce structural and thermodynamic properties of HS [80]. Here, k and T are the Boltzmann constant and absolute temperature, respectively.

Computation of time dependent correlation functions and viscosities
Transport coefficients in general and viscosity in particular can be computed towards MD simulations and the GK formalism. To this end, one can use two different approaches: one described by equilibrium conditions and the other by nonequilibrium conditions. For instance, in equilibrium MD simulations one need to compute integrals of time-dependent correlation functions [66] or analyse the long time-behaviour of the so-called Einstein relations [81,82]. Whereas in nonequilibrium MD simulations, the so called periodic perturbation method [83] is usually employed. Using a numerical calibrated continuous function that mimics the HS, Pousaneh and Wijn have determined shear viscosity as a function of density number [78].
Despite some technical issues, that need to be considered and carefully addressed in non-equilibrium dynamics, such as boundary conditions [36,84] and the extrapolation to the thermodynamic equilibrium [83], the major limitation of this approach is that only one transport coefficient can be computed during a simulation [85]. In contrast, in equilibrium molecular dynamics one can compute all thermodynamic and transport coefficients in a single simulation. Furthermore, since the interaction potential considered in this work is derived from physical arguments, we expect reliable and accurate results. This will be verified by comparing our predictions with HS results over a broad range of densities.
Thus, in this work we use equilibrium MD simulations to determine shear, bulk and longitudinal viscosities. The GK shear viscosity, η, is computed as where V is the volume, ⟨τ αβ (0) τ αβ (t ′ )⟩ is the shear stress autocorrelation function and angular brackets stands for an equilibrium ensemble average. The instantaneous stress tensor elements are defined as [86] where v i,α is the α component of i-th particle and r ij,α is the respective component of the separation distance between i and j particles, respectively. Besides, F ij,β is the β component of the force over i due to j particle. Its worth to stress out that the first and second terms of equation (4) are identified as the kinetic τ αβ, k and potential τ αβ, p contributions of the stress tensor, respectively. Hence, correlation function in equation (3) contains purely kinetic, potential and cross kinetic-potential contributions. Therefore, identifying and splitting shear stress autocorrelation function contributions, the shear viscosity can be expressed as the sum of those contributions. This allows us to write the shear stress as [85] η (ρ, where subscripts k, kp and p stands for kinetic, cross kineticpotential and potential contributions, respectively. In a similar way, bulk viscosity is defined as [87] where δP (t ′ ) is the instantaneous pressure fluctuation defined ∑ τ αα being the sum of the pressure tensor trace at time t' and ⟨P⟩ being the average pressure. Finally, the uniaxial volume viscosity or commonly named longitudinal viscosity is defined as where δτ αα terms are the fluctuations of the αα stress tensor trace elements, defined in a similar way that the instantaneous pressure fluctuation previously mentioned.
From an experimental point of view, longitudinal viscosity is considered as the counterpart of shear viscosity. In fact, from rotational invariance, it is theoretically predicted that longitudinal viscosity is related to shear and bulk viscosities as [64][65][66][67] giving us an indirect and alternative way to determine η L .

Computer simulation scheme
Properties of interest are determined using equilibrium MD simulations in the canonical ensemble, where the number of particles N, simulation domain volume V and system temperature T are fixed. All simulations consist of N = 10 976 spherical particles of mass m and diameter σ initially distributed in a FCC array with random initial velocities given by the Maxwell-Boltzmann distribution for the desired temperature T. The mass, diameter and the interaction energy parameter ϵ are used to define common dimensionless units of length, temperature and time as r * = r/σ, T * = kT/ϵ and t * = √ σ 2 m/ϵ, respectively. Using the density number, ρ * = σ 3 N/V, the size of the cubic simulation domain is defined. Periodic boundary conditions and minimum image convention were applied in all simulations. Additionally, results are analysed and discussed in terms of the volume fraction which is defined as φ = ρ * π /6. Results covers a range from 0.04 to 0.50 of φ values. In figure 1 we shown two snapshots of the pseudo hard-sphere fluid microscopic state at low and high volume fractions, namely (left) φ = 0.10 and (right) φ = 0.50, respectively.
The equations of motion were integrated with the standard velocity Verlet algorithm [88] employing a time step of ∆t * = 1 × 10 −3 t. In each simulation 1.2 × 10 7 integration time steps were carried out, from which, the first 2 × 10 6 time steps were used to reach thermal equilibrium, then, remaining steps were used to compute correlation functions of interest and gather statistics. Here, correlation functions were computed at each time step (∆t * = 1 × 10 −3 t), for a maximum correlation time of t * = 10. Finally, in order to enhance statistical results, ten

Shear stress correlation functions of pseudo hard-spheres
Using positions, velocities and forces of PHS particles provided by the MD simulations, the stress tensor is evaluated at every integration time step as it is indicated by equation (4). Then, shear stress, pressure fluctuation and trace fluctuation autocorrelation functions were computed. Additionally, kinetic, potential and cross kinetic-potential shear stress autocorrelation functions were also computed. Finally, these results are used in the GK formalism to determine the transport coefficients of interest.
Although shear viscosity is one of the most studied transport coefficients and despite that its computation is usually carried out within the GK approach, the behaviour of autocorrelation functions that defined it are usually not discussed. This could be, in general, due to its fast decay to zero or large fluctuations observed in small systems. Then, in figure 2 we summarize normalized results of shear autocorrelation function at different volume fractions ranging from low to high values. Albeit at low volume fractions (φ < 0.10) this functions decays in a very short time, one can observe a sudden decay at times t * < 10 −2 , then, function go to zero is a slow fashion. This is an indication of the low rate of interactions between particles in this regime, particularly, concerning to kinetic and kinetic-potential contributions of the shear stress.
It can be seen from equation (4), that the stress tensor is defined by a kinetic and potential contributions. Thus, in addition to calculating the complete shear stress autocorrelation function, it is also possible to separate and analyse the individual contributions from kinetic, potential and cross kinetic-potential components. Such contributions can be computed with the same relation than the one given by equation (5) and will allow to identify the influence of each contribution on shear viscosity as a function of the volume fraction. Therefore, kinetic, potential and kinetic-potential normalized shear stress autocorrelation functions are shown from top to bottom panels in figure 3, respectively. We stress out that this correlations are determined during the simulations in an independent way from the total correlation shown in figure 2. However, we corroborate that the last is recovered if we add up the three contributions.
Even though it is expected, the first remarkable fact is that both dynamic and potential components exhibits a nonzero correlation at short times. However, dynamic contribution decay depends on φ, while, decay of the potential contribution is almost independent of φ. This can be corroborated in the top and middle panels of figure 3, respectively. In contrast, a strong correlation in the kinetic-potential contribution is observed at times between 10 −2 and 10 0 . In fact, the correlations reach a maximum value near to t * < 10 −2 . After this, correlations goes to zero with a rate that depends on φ. The decay is slow for φ < 0.1 and very quick at φ > 0.4. It is worth mentioning that the observed magnitudes in the shear stress autocorrelation function depends on normalization, however, its observed behaviour does not depend on the normalization.
The determination and analysis of shear stress contributions allow us to identify the components that dominates the total shear stress over the time and φ. For instance, the slow decay observed at φ < 0.1 in figure 2, is mainly due to the kinetic contribution and in a second instance from kinetic-potential contributions (see the top and bottom panels of figure 3). Moreover, from this results, is expected that shear stress would have an important contribution coming form its kinetic part at low φ, while, at high densities, the potential contribution becomes more relevant and dominant. This will be more evident during the explicit discussion of shear stress and its contributions in the next section.

Shear viscosity
Results of shear viscosity as a function of φ for the PHS fluid are shown in figure 4, where, for comparison purposes, we have included results for the HS [35] and PHS [78] fluids. Although the magnitude of densities in these studies are not exactly the same as the ones considered in this work, we discovered a generally good agreement between the HS and PHS fluid results.
If we consider the results of Sigurgeirsson and Heyes [35] for N = 4000 as the reference, we can observe that the predictions by Pousaneh and Wijn [78] underestimate the target values. However, these differences seem to be smaller at high volume fractions, φ > 0.3. On the other hand, in the lowdensity regime, the results of Pousaneh and Wijn deviate from the zero-density limit. In contrast, our results capture mentioned limit with better agreement. However, we slightly overpredict the reference values. Nonetheless, these differences  [35] and numerically deducted PHS fluid [78], respectively. In addition, we show the zero density limit of this property (diamond) and the practical correction of the CE shear viscosity (dashed line), see equation (9). (bottom) Kinetic, η k , potential, ηp, and kinetic-potential, η kp , shear viscosities as a function of the volume fraction. become smaller as the density increases. Differences between approaches could be attributed to finite size effects, as discussed by Sigurgeirsson and Heyes [35]. In fact, they have demonstrated that the determination of viscosity can depend on the number of particles used in simulations.
On the other hand, according to the kinetic theory, it is well known that the shear viscosity in the limit of zero density is independent of density and is given by This limit value is included in figure 4 to illustrate the excellent agreement between our results and theoretical predictions.
Turning our attention to the explicit density dependence of this transport coefficient, we employ a similar approach than the one proposed by Krieger-Dougherty [89] to express the normalised shear viscosity in terms of its zero density limit.
We referrer to this approach as CE correction, and we express the shear viscosity as follows where both φ c and m are fitting parameters, while the upper script stands for the CE approach. Thus, by adjusting simulation results, we obtain φ c = 0.608 and m = 2.161, respectively. These values are in agreement with those previously reported by Sigurgeirsson and Heyes [35]. Moreover, it is worth stressing that HS results of these authors are limited to φ ⩾ 0.20, whereas we report values at lower densities. This allows us to obtain better predictions using equation (9); indeed, its results are shown by a dashed line in the top panel of figure 4. The shear viscosity, its contributions, as well as the bulk and longitudinal viscosities, are summarised in table 1 of the appendix. Now, utilizing the contributions of the shear stress autocorrelation, we determine the kinetic, potential and kineticpotential contributions as a function of φ. The results are depicted in the bottom panel of figure 4. The first noteworthy observation is that the kinetic and potential contributions differ by almost two orders of magnitude in both the low (φ < 0.10) and high (φ > 0.40) volume fraction regimes, respectively. In the former regime, the kinetic contribution dominates over the potential one, as expected from the kinetic theory. Conversely, in the latter regime, the potential term becomes the dominant and most significant contribution. This is attributed to the decreasing particle motility at high densities. Clearly, this is a direct consequence of the shear stress correlation behaviour previously discussed in figure 3.
Additionally, it is observed that the crossover between kinetic and potential contributions takes place near φ = 0.20. However, it is not until φ > 0.40 that shear viscosity is almost fully defined by its potential contribution. As for the kineticpotential viscosity, its magnitude lies between the kinetic and potential contributions throughout the entire range of studied densities, with only a small increase observed as the volume fraction rises. Despite the respective differences between both systems, these results support the assumption that the potential component is the main contributor to shear viscosity in nearhard-sphere colloids [39]. However, as we already pointed out, this in only valid at high volume fractions.

Bulk viscosity
Bulk viscosity is determined through equation (6), and its results as a function of φ are shown in figure 5. In the same figure, the results computed by Sigurgeirsson and Heyes [35] for HS are also presented. Is worth recalling that this transport coefficient only has collisional contributions, meaning that there are only potential contributions. In fact, its zero-density limit value is zero. Thus, even when we compute pressure fluctuations correlations at very low densities, the calculation of η B exhibits large fluctuations at φ < 0.10. This is due to both a very fast decay of correlation and possible finite size effects. We expect that such effects will be more pronounced for smaller systems than the one used here.
With above considerations in mind, and despite the mentioned density limit, we propose an equation similar to equation (9). By using the zero-density limit of the shear viscosity, we provide a practical correction of the CE approach, which can be expressed as where parameters φ 1 and m 1 are fitted to simulation results. Then, for a volume fraction range from φ = 0.10 to φ = 0.50, we obtain φ 1 = 0.571 and m 1 = 2.143, respectively. Results of equation (10) are displayed in figure 5 with a (red) dashed line.
As is shown in figure 5, our results are in good agreement with those for the HS fluid. However, in the highdensity regime, a slight deviation can be observed between the two simulation approaches. This deviation could be attributed to the different number particles used in each case. Moreover, in regard to the predictions of equation (10), is not surprising to find a lack of agreement. Firstly, we use the shear stress zero-density limit to express bulk viscosity, and secondly, equation (10) has only two adjustable parameters. Additionally, is worth to mention that the results of Sigurgeirsson and Heyes are based on simulations with only N = 500 particles, whereas we employ at least two orders of magnitude more particles. In fact, obtained values of φ 1 and m 1 strongly depends on N.Nevertheless, it is expected that finite size effects would be more significant in the high-density regime. However, for the time being, such a study is beyond the scope of the present work.

Longitudinal viscosity
Although longitudinal viscosity, η L , can be considered as the counterpart of the shear viscosity, it is by far the less studied viscosity. This is clearly due to its complex experimental determination [47,48]. For this reason, this transport coefficient is almost exclusively determined through computer simulations [90] or indirectly, as indicated by equation (8).
Here, the longitudinal viscosity is using equations (7) and (8). Results for this quantity as a function of φ are summarised in figure 6. In the same figure, ED-MD simulation data determined by Srivastava and colleagues [90] for the HS fluid are shown. Additionally, simulation results are used to propose an empirical correction to the Chapman-Enskog low-density limit η 0 , similar to the previous case for η L . Thus, the proposed CE longitudinal viscosity correction is given by where fitting parameters are φ 2 = 0.577 and m 2 = 2.389, respectively. Results of equation (11) are shown in figure 6 with a (red) dashed line. Although simulation results for η L are scare, in figure 6, we can observe a good agreement between our model and available results for thr HS fluid [90]. On one hand, it is observed that the predictions fo equation (8) overestimates both simulation results. On the other hand, the CE longitudinal viscosity describes the results very well, but only at moderate densities, φ > 0.20. However, at lower densities, significant deviations are seen. This behaviour is not surprising since we are using a very simple functional form to correlate simulation data and, most importantly, we are taking the zero-density limit value of shear viscosity as the reference point. Therefore, η CE L can only be considered as a qualitative description for this transport coefficient.
It is worth mentioning that Srivastava and co-workers also proposed an empirical equation to represents their simulation results [90]. However, their approach is based on the HS compressibility factor, namely, η L ≈ α exp (γZ), and a correction of shear viscosity at the zero-density limit. We verified that such an approach also represents our simulation results very well and can be used in a similar manner to express other viscosities of interest. In summary, equation (2) would be useful and accurate to determine not only structural and thermodynamic properties of the HS fluid but also for determining transport coefficients that depend sensitively on the number of particles.

Concluding remarks
By conducting equilibrium MD simulations and using the GK formalism, we have computed transport coefficients such as shear stress viscosity, as well as their contributions, bulk viscosity and longitudinal viscosity for the PHS fluid. Our results were compared to existing simulation predictions for the HS fluid, and we found excellent agreement. Furthermore, we were able to provide additionally values of the transport coefficients of interest across a wide range of densities.
In addition, based on the Krieger-Doughter shear viscosity equation for colloids, we use our simulation results to propose CE-like equations for shear, bulk and longitudinal viscosities. These equations represents density-dependent corrections of the zero-density limit of the shear viscosity. We have shown that such approaches works very well at moderate densities. However, their simple functional form is insufficient to accurately reproduce bulk and longitudinal viscosities in the low density regime. This emphasizes the need for further theoretical developments to accurately represent bulk and longitudinal viscosities in the low-density regime.
Our analysis also demonstrates that the PHS pair potential, derived from physical arguments, can accurately reproduce the dynamical behaviour of the HS fluid and its transport coefficients. Furthermore, due to the continuous functional form of the interaction potential, we were able to explicitly determine and analyse the contributions to shear stress. Through this analysis, we have confirmed that in the low-density regime, the shear viscosity kinetic contribution dominates this transport coefficient, while in the high-density regime, the potential contributions takes precedence. Moreover, we have identified the crossover between these to occur around φ = 0.20. These results supports the notion that high densities, shear viscosity is primarily determined by its potential contribution. Therefore, employing the approach described here, it will possible to systematically determine and study shear viscosity at high densities, such as near the freezing point of in the vicinity of a glass transition, in complex fluids such as polydisperse systems, rigid or flexible hard-sphere chains and hard-sphere colloidal systems. Work in this direction is in progress.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Acknowledgment
L A N-C and C A R-M acknowledges financial support provided by CONACYT through grant no. 709987 and 'Estancias Posdoctorales por México para la Formación y Consolidación de las y los Investigadores por México'; A T C acknowledge the postdoctoral-scholarship granted by CTIC-DGAPA-UNAM (2022). The authors thankfully acknowledge the computer resources, technical expertise and support provided by the Laboratorio Nacional de Supercómputo del Sureste de México, CONACYT member of the network of national laboratories.

Conflict of interest
The authors have no conflicts to disclose.

Appendix. Molecular dynamics viscosities of the pseudo hard-sphere fluid
In this appendix we summarise the molecular dynamics results of shear viscosity, its contributions, bulk viscosity and longitudinal viscosity of the pseudo hard-sphere fluid at different volume fractions.