Vibrational Model of Heat Conduction in a Fluid of Hard Spheres

Abstract

Application of a vibrational model of heat transfer to a fluid made of hard spheres is discussed. The model was originally proposed to describe heat conduction in fluids with soft pairwise interactionsHere, it is shown that only minor modifications are required to apply the model in the opposite limit of hard sphere interactions. Good agreement with recent results from molecular dynamics simulation is documented in the moderately dense regime. Near the freezing point, however, the model overestimates the thermal conductivity coefficient (by ≃$50\%$). The new approach is compared with other simple models for the thermal conductivity coefficients such as Bridgman’s expression and the Enskog formula. The value of the coefficient in the Bridgman’s expression, appropriate for the hard sphere fluid, is determined. A new expression for the dependence of the reduced thermal conductivity coefficient on the reduced excess entropy is proposed. The obtained results can be useful for rough estimates of the thermal conductivity coefficient of simple fluids with steep interactions when more accurate experimental results are not available.

1. Introduction

Our understanding of dynamical and transport properties of liquids remains incomplete and fragmented even in the case of the simplest atomic systems. Despite considerable progress achieved over the years [1,2,3,4,5,6], we still often need to apply various phenomenological approaches, semi-quantitative models, and scaling relations. Among the renowned examples are the Stokes-Einstein relation between the self-diffusion and the shear viscosity coefficients [4,7,8,9,10,11,12,13,14,15,16,17], the excess entropy scaling of reduced transport coefficients [18,19,20,21,22], and several variants of the freezing temperature scaling [23,24,25,26,27]. Recently, a freezing density scaling of reduced transport coefficients of the Lennard-Jones and related fluids has been discussed as a new useful addition to the existing approaches [28,29,30].

Heat transfer in fluids (throughout this paper the term “fluid” is applied to liquids and supercritical fluids) is an important topic of contemporary research with diverse interdisciplinary applications [31,32,33,34,35,36]. A vibrational model of heat transfer has been recently applied to fluids with soft pairwise interactions [37]. In this model, the thermal conductivity coefficient is proportional to the effective frequency of atomic vibrations around the local temporary equilibrium positions. The effective frequency is related to the properties of the liquid collective excitation spectrum. The model has been applied to quantify heat transfer in a strongly coupled one-component plasma fluid, dense Lennard-Jones liquid, and plasma-related screened Coulomb (Yukawa) fluid. An impressive agreement with the available results from numerical simulations has been documented [37,38,39]. The purpose of this work is to consider the applicability of the vibrational model to a new physical system—a fluid made of hard spheres (HS fluid). The vibrational picture of atomic dynamics is clearly not very appropriate for extremely anharmonic HS system. Nevertheless, collective excitations are supported in the dense HS fluid and averaging over collective modes can formally be performed. The question is whether the results have any relation to reality. It will be shown below that a formal expression based on the vibrational model of heat transfer is in agreement with the results from molecular dynamics simulations in a moderately dense regime, but overestimates the numerical results near the freezing point. The accuracy of some other simplistic approaches will be also examined. In particular, the coefficient in the Bridgman’s expression, appropriate for the HS fluid, will be identified. The excess entropy scaling of the thermal conductivity coefficient in the HS fluid will be briefly discussed and a new expression will be proposed. The results can be useful for simple rough theoretical estimates when more accurate data are unavailable.

2. Hard Sphere Fluid

The HS interaction potential is defined as:

$$\varphi \left(r\right)=\left\{\begin{array}{cc}\infty ,\hfill & r<\sigma \hfill \\ 0,\hfill & r\ge \sigma ,\hfill \end{array}\right.$$

(1)

where $\sigma$ is the sphere diameter and r is the distance between the centres of two spheres. The potential ensures that spheres cannot overlap. The HS system is a very important simple model for the properties of condensed matter in its various states [40,41,42,43,44,45,46,47].

The thermodynamic and transport properties of HS systems depend on a single reduced density parameter, ${\rho }^{*}=\rho {\sigma }^3$, or the packing fraction, $\eta =\pi \rho {\sigma }^3/6$, where $\rho$ is the density. Transport properties of the HS fluid have been extensively studied (see e.g., Ref. [42] for a review). Here we use recent molecular dynamics (MD) simulation results for the thermal conductivity coefficient reported by Pieprzyk et al. [48]. The use of large simulation systems and long simulation times allowed accurate prediction of the thermal conductivity coefficient in the thermodynamic limit.

Throughout the paper we use macroscopically reduced units for the thermal conductivity coefficient:

$${\lambda }_{\mathrm{R}}=\lambda \frac{{\rho }^{-2/3}}{v_{\mathrm{T}}},$$

(2)

where $\lambda$ is the dimensional thermal conductivity coefficients, $v_{\mathrm{T}}=\sqrt{T/m}$ is the thermal velocity, T is the temperature in energy units ($\equiv k_{\mathrm{B}}T$), and m is the atomic mass. This normalization is essential in Rosenfeld’s excess entropy scaling approach [18,19], therefore the subscript $\mathrm{R}$ is often used. This normalization does not contain the Boltzmann constant $k_{\mathrm{B}}$, because the temperature is measured in energy units. Another normalization that is commonly employed is:

$${\lambda }^{*}=\lambda \frac{{\sigma }^2}{v_{\mathrm{T}}}.$$

(3)

These two normalizations are related via ${\lambda }^{*}={\lambda }_{\mathrm{R}}{\left({\rho }^{*}\right)}^2$.

3. Vibrational Model

The vibrational paradigm of atomic dynamics in dense liquids has been discussed by many authors over the years, see e.g., Refs. [1,7,49,50,51]. The main assumptions involved are as follows [16]: Atoms exhibit solid-like oscillations about temporary equilibrium positions corresponding to a local minimum on the system’s potential energy surface [1,51]. These positions do not form a regular lattice like in crystalline solids, but correspond to a liquid-like structure [50]. They are also not fixed and change (diffuse or drift) with time, which allows liquids to flow. However, these rearrangements occur on much longer time scales than those characterising solid-like oscillations. Effectively, one can assume that the local configurations of atoms are preserved for some time until a fluctuation in the kinetic energy allows the positions of some of these atoms to be rearranged towards a new local minimum in the multidimensional potential energy surface. This picture allows to make important approximations about the properties of atomic motion and the mechanisms of momentum and energy transport in the liquid state. For instance, the Stokes-Einstein (SE) relation without the hydrodynamic radius naturally emerges within this vibrational paradigm under a few additional assumptions [7,13,15,16].

For energy transfer, separation of time scales corresponding to fast solid-like atomic oscillations and their slow drift plays a dominant role. Namely, it is reasonable to assume that a vibrating atom transports energy from its hotter to its cooler neighbors with a characteristic energy exchange rate equal to its average vibrational frequency, $\nu =\langle \omega \rangle /2\pi$. Each atom controls the energy transfer through an area of order ${\Delta }^2$, where $\Delta ={\rho }^{-1/3}$ is the average interatomic separation. A liquid can be approximated by a quasi-layered structure with quasi-layers that are perpendicular to the temperature gradient (applied along x-axis) and are separated by the distance $\Delta$ (see e.g., Figure 1 from Ref. [37]). The energy difference between two neighbouring layers is $(dU/dx)\Delta$. The energy flux density can be approximated as:

$$j\simeq -\frac{\nu }{{\Delta }^2}\left(\frac{dU}{dx}\right)\Delta =-\frac{\langle \omega \rangle }{2\pi \Delta }{c}_{\mathrm{p}}\frac{dT}{dx},$$

(4)

where $c_{\mathrm{p}}={(dU/dT)}_{\mathrm{p}}$ is the specific heat at constant pressure and the minus sign implies that the heat flow is down the temperature gradient. On the other hand, Fourier’s law of the heat transfer reads:

$$j=-\lambda \frac{dT}{dx}.$$

(5)

By comparing Equations (4) and (5) we immediately obtain:

$$\lambda =c_{\mathrm{p}}\frac{\langle \omega \rangle }{2\pi \Delta }.$$

(6)

This is essentially the expression derived in Ref. [37], except that the specific heat at constant volume, $c_{\mathrm{v}}$, appeared there. The present choice seems in general more appropriate because pressure should be constant in equilibrium. The difference is insignificant for soft spheres, because dense fluids can be considered as essentially incompressible in a wide portion of their phase diagram not too far from the freezing curve and thus $c_{\mathrm{p}}\simeq c_{\mathrm{v}}$ holds. For example, Equation (6) with $c_{\mathrm{p}}\to c_{\mathrm{v}}$ works very well for soft plasma-related Coulomb and Yukawa interactions and the Lennard-Jones fluid [37,38,39]. In such cases it is more appropriate to use $c_{\mathrm{v}}$ for practical estimates (since $c_{\mathrm{v}}$ is normally easier to evaluate). This is not the case, however, in the HS limit. For hard spheres $c_{\mathrm{v}}\equiv 1.5$ and $c_{\mathrm{p}}$ and $c_{\mathrm{v}}$ can differ very considerably.

Before we discuss heat conduction further, it should be emphasized that the vibrational picture of atomic dynamics is not particularly suitable for the HS fluid. In contrast to softer interactions, atomic motion does not exhibit a pronounced oscillatory character [52,53,54]. The HS fluid is extremely anharmonic and dynamics is dominated by short hard-core-like collision events [55]. Nevertheless, similar to other simple fluids, the dense HS fluid does support the longitudinal and transverse collective excitations [56] (with a forbidden long-wavelength region for the transverse mode, the so-called “k-gap” [56,57,58,59,60]). Therefore, averaging over collective modes can be performed and in this sense $\langle \omega \rangle$ remains a meaningful quantity. The Stokes-Einstein relation holds in the HS fluid, just as it does in soft sphere fluids, where it emerges naturally as a consequence of the vibrational picture of atomic dynamics [7,16]. Vibrational model allows to estimate the excess entropy of various fluids with pairwise repulsive interactions near the fluid–solid phase transition, including the HS limit [61]. Under these circumstances it is not very unreasonable to compare the predictions of Equation (6) with the available results from MD simulations in the HS limit.

Since the actual frequency distribution is often not known in liquids, and can vary from one type of liquid to another, some simplifying assumptions are necessary to evaluate $\langle \omega \rangle$ in Equation (6). In the simplest Einstein approximation all atoms vibrate with the same (Einstein) frequency ${\Omega }_{\mathrm{E}}$. We arrive at the expression derived by Horrocks and McLaughlin [62]. But this approximation is obviously not applicable to the HS fluid. The conventional expression for the Einstein frequency,

$${\Omega }_{\mathrm{E}}=\frac{\rho }{3m}\int d\mathbf{r}{\nabla }^2\varphi \left(r\right)g\left(r\right),$$

(7)

where $g\left(r\right)$ is the radial distribution function (RDF), is simply undefined in the HS limit.

One can use an acoustic spectrum, $\omega =k{c}_s$ with an appropriate maximum cutoff wavenumber $k_{\max }$. Here $c_s$ is the conventional adiabatic sound velocity. Then, the standard averaging procedure would result in:

$$\lambda \sim \frac{c_s}{{\Delta }^2}.$$

(8)

This can be considered as an analog of Bridgman’s equation [63]. The numerical coefficient remains undefined, because of the qualitative character of the arguments involved (the values between 2 and 3 can be found in the literature, see e.g., Ref. [35]). We shall see that for the HS fluid the numerical coefficient is in fact close to unity.

As a more physically sound approximation, assume that a dense liquid supports one longitudinal and two transverse collective modes. Debye-like averaging in k-space can be performed using acoustic asymptotes ${\omega }_l\left(k\right)=c_{l}k$ and ${\omega }_t\left(k\right)=c_{t}k$ in the long-wavelength domain, where $c_l$ and $c_t$ are the longitudinal and transverse sound velocities, respectively. The thermal conductivity coefficient becomes [37]:

$$\lambda \simeq \frac{1}{4}{\left(\frac{3}{4\pi }\right)}^{1/3}{c}_{\mathrm{p}}\frac{c_l+2c_t}{{\Delta }^2}.$$

(9)

If we substitute $c_{\mathrm{p}}\simeq c_{\mathrm{v}}\simeq 3$ near the freezing point (according to Dulong-Petit law), we get a formula, which is similar to that of the minimal thermal conductivity model proposed by Cahill and Pohl for amorphous solids [64,65,66]. It turns out that a similar formula obtained using the vibrational model works rather well in the liquid regime. It is the expression (9) that will be compared with the recent MD data on the thermal conductivity coefficient of the HS fluid. Note that there are no free parameters in the discussed approach.

4. Methods

To evaluate the thermal conductivity coefficient from Equation (9) we need to know the longitudinal and transverse sound velocities ($c_l$ and $c_t$) and the specific heat at constant pressure ($c_{\mathrm{p}}$). The later is obtained from the equation of state of the HS fluid. Here a simple but sufficiently accurate equation of Carnahan and Starling (CS) is used [67]. The pressure is written as:

$$P(\rho ,T)=\rho TZ\left(\eta \right),Z\left(\eta \right)=\frac{1+\eta +{\eta }^2-{\eta }^3}{{(1-\eta )}^3},$$

(10)

where $Z\left(\eta \right)$ is the CS compressibility.

The various sound velocities are related to the elastic moduli. For example, the adiabatic sound velocity $c_s$ is related to the adiabatic bulk modulus $K_s$ via $c_s=\sqrt{K_s/m\rho }$. The adiabatic bulk modulus $K_s$ is determined from the CS compressibility factor Z [68,69]:

$$K_s=\rho T\left[Z\left(\eta \right)+\eta dZ\left(\eta \right)/d\eta +\frac{2}{3}{Z}^2\left(\eta \right)\right].$$

(11)

Similarly, the longitudinal and transverse sound velocities are obtained from the longitudinal and transverse instantaneous (infinite frequency) elastic moduli, $c_l=\sqrt{M_{\infty }/m\rho }$ and $c_t=\sqrt{G_{\infty }/m\rho }$. Moreover, the longitudinal elastic modulus is related to the bulk and shear moduli via:

$$M_{\infty }=K_{\infty }+\frac{4}{3}{G}_{\infty }.$$

(12)

The behavior of elastic moduli when approaching the HS interaction limit has recently been discussed in detail [60]. For the instantaneous shear and bulk moduli, the expressions derived by Miller [70] are employed, with minor modifications discussed in Ref. [71]. The shear modulus is [70,71]:

$$G_{\infty }=\rho T\left[1-\frac{8}{5}\eta g^{\prime }\left(1\right)\right],$$

(13)

where $g^{\prime }\left(1\right)$ denotes the reduced derivative at contact, $g^{\prime }\left(1\right)={lim}_{ϵ\to 0}{\left[dg\left(x\right)/dx\right]}_{x=1+ϵ}$ with $x=r/\sigma$. For $g^{\prime }\left(1\right)$ we use the approximation proposed in Ref. [72]:

$$g^{\prime }\left(1\right)=-\frac{9\eta (1+\eta )}{2{(1-\eta )}^4}.$$

(14)

The bulk modulus is [70,71]:

$$K_{\infty }=2P-\frac{8}{3}\rho T+\frac{2}{3}\frac{P^2}{\rho T}+\frac{5}{3}{G}_{\infty }.$$

(15)

The longitudinal modulus is then obtained from Equation (12). This allows to evaluate all necessary quantities and compare the results of various approaches with those from MD simulations. Results are reported in the next Section.

5. Results

Figure 1 shows the comparison between the MD data on the reduced thermal conductivity coefficient and the discussed theoretical approximations. Symbols correspond to the MD data [48]. The quantitative picture is conventional and resembles that in other simple fluids, see e.g., Refs. [30,73]. In the low-density domain ${\lambda }_{\mathrm{R}}$ decreases with density. In the first approximation, the thermal conductivity coefficient of a dilute HS gas is ${\lambda }_0=(75/64{\sigma }^2){(T/\pi m)}^{1/2}$ [74], which leads to ${\lambda }_{\mathrm{R}}\simeq 0.661/{\left({\rho }^{*}\right)}^{2/3}$ [28]. The minimum of ${\lambda }_{\mathrm{R}}$ is reached at ${\rho }^{*}\simeq 0.25$. The amplitude of this minimum is known to be quasi-universal for many model and real fluids. Numerically, a rough estimte ${\lambda }_{\mathrm{R}}^{\min }\sim 3$ holds for various liquids, except for plasma-related fluids with very soft Coulomb-like interactions [30,73]. This minimum corresponds to the crossover between two different mechanisms of energy transfer: pairwise collisions in the dilute gas regime and collective vibrations in the dense fluid regime. Note that a similar minimum appears in the reduced viscosity coefficient [30,73]. It is the location of the reduced viscosity minimum that previously served as one of the definitions of the gas-liquid dynamical crossover [75,76]. After the minimum, ${\lambda }_{\mathrm{R}}$ increases monotonously with density. The typical value of the thermal conductivity coefficient at the freezing point of various simple liquids is ${\lambda }_{\mathrm{R}}^{\mathrm{fr}}\sim 10$ [30,73]. We observe in Figure 1 that the thermal conductivity coefficient of the HS fluid at the freezing point is around $\simeq 30\%$ higher. A similar observation was previously made in Ref. [28].

The dashed curve in Figure 1 corresponds to the vibrational model of heat transfer as formulated above. The agreement is reasonably good at moderate densities (after minimum in ${\lambda }_{\mathrm{R}}$ is reached), but becomes worse at ${\rho }^{*}\gtrsim 0.7$. As the freezing density (${\rho }_{\mathrm{fr}}^{*}\simeq 0.939$ according to Ref. [77]) is approached the vibrational model overestimates ${\lambda }_{\mathrm{R}}$ by $\simeq 50\%$. This reflects the fact that the vibrational picture is not very appropriate for the dense HS fluid. It is to some extent surprising that the expression based on the averaging over collective modes remains a relatively good approximation in the moderately dense regime.

The dash-dotted curve corresponds to the Bridgman’s expression. Plotted is the value ${\lambda }_{\mathrm{R}}=c_s/v_{\mathrm{T}}$, where $c_s$ is evaluated using Equation (11). We see that the numerical factor equal to unity is appropriate for the dense HS fluid. Only near the freezing point it leads to some underestimation of ${\lambda }_{\mathrm{R}}$.

The solid curve corresponds to the Enskog theory, which is an extension of the gas kinetic theory for the transport coefficients to finite densities. It accounts for the excluded volume effects, but assumes that the successive collisions between the atoms are uncorrelated. The Enskog expression for the thermal conductivity coefficient of the HS fluid reads [48]:

$$\lambda =1.02513{\lambda }_0\rho B_2\left(\frac{1}{Z_{\mathrm{ex}}}+1.2+0.7574Z_{\mathrm{ex}}\right),$$

(16)

where ${\lambda }_0$ is the first approximation for $\lambda$ in the infinitely dilute gaseous limit (see above), the coefficient $1.02513$ is the Sonine polynomial correction up to 4th term, $B_2=2\pi {\sigma }^3/3$ is the second virial coefficient, and $Z_{\mathrm{ex}}=Z-1$ is the excess compressibility factor. It is observed that overall the thermal conductivity coefficient of the HS fluid is rather well described by the Enskog theory. Deviations from MD data of Ref. [48] do not exceed ∼5% in the entire fluid regime. This is in contrast to other transport properties: Agreement with Enskog theory for the self-diffusion and viscosity coefficients only takes place at sufficiently low densities [77].

6. Excess Entropy Scaling

In view of the importance of the excess entropy scaling for transport coefficients in fluids, the reduced thermal conductivity ${\lambda }_{\mathrm{R}}$ is replotted in Figure 2 as a function of the excess entropy. Excess entropy is defined as the actual entropy minus the entropy of an ideal gas at the same temperature and density, divided by the number of particles and Boltzmann constant $k_{\mathrm{B}}$ to make it dimensionless. It is negative in fluids, because they are more ordered than the ideal gas. The excess entropy can be easily calculated using the CS compressibility:

$$s_{\mathrm{ex}}={\int }_0^{\eta }\frac{1-Z\left(\eta \right)}{\eta }d\eta =\frac{\eta (3\eta -4)}{{(1-\eta )}^2}.$$

(17)

In Figure 2 symbols are the MD numerical results from Ref. [48]. The dashed curve corresponds to a “quasi-universal” exponential scaling proposed by Rosenfeld, ${\lambda }_{\mathrm{R}}\sim 1.5exp(-0.5s_{\mathrm{ex}})$ [19]. The advantage of this scaling is that it allows to estimate reasonably well (within about $\simeq 30\%$) the thermal conductivity coefficient in various systems with quite disparate pair interactions [19]. On the other hand, different potentials can be fitted better by somewhat different functional forms (see e.g., Ref. [22]). For the HS fluid we propose a simple cubic polynomial fit of the form

$${\lambda }_{\mathrm{R}}=2.3845+0.0272s_{\mathrm{ex}}+0.6550s_{\mathrm{ex}}^2-0.0307s_{\mathrm{ex}}^3.$$

(18)

This fit does a very good job for $s_{\mathrm{ex}}\lesssim -1$, as shown by the solid curve in Figure 2.

7. Conclusions

The system composed of hard spheres is a very important model in statistical mechanics and condensed matter physics. It allows to understand many generic mechanisms behind quasi-universal structural and dynamical properties of real materials. In this paper we discussed the properties of thermal conduction in the HS fluid. It has been demonstrated that the expression obtained within the vibrational model of atomic transport provides a reasonable description of the thermal conductivity coefficient at moderate densities, but leads to a significant overestimation (up to ≃$50\%$) near the freezing point. Overall, this is not particularly surprising, because the extremely anharmonic HS system does not fit into the vibrational paradigm and the vibrational model was not originally intended for the HS fluid. On the other hand, it should be noted that the Stokes-Einsten relation, emerging naturally within the vibrational paradigm of atomic transport, is satisfied to a good accuracy in the dense HS fluid. The vibrational model allows to estimate the excess entropy for fluids with repulsive interactions, including the HS limit. Thus, the vibrational approach seems not completely irrelevant and this represents an important lesson from this study. We have also demonstrated that the reduced thermal conductivity coefficient correlates well with the reduced sound velocity, and the Bridgman’s expression applies. The coefficient of proportionality is close to unity, much smaller than usually quoted. The Enskog expression is known to work rather well in the entire HS fluid density range, but it is not directly applicable to other system as the vibrational model does. Yet, the vibrational model is clearly much better suited for soft pairwise interactions.