Introduction

Statistical mechanics is the art of predicting the behaviour of a system with a large number of degrees of freedom, given the laws governing its microscopic behaviour. The statistical description of liquids, in comparison with the solid and gas phases, is incomplete. The problem of formulating a rigorous mathematical description of liquids has always been regarded as much more difficult than that of the kinetic theory of gases or solid-state theory, stimulating the ongoing research1,2,3,4,5,6,7,8,9. Owing to the simultaneous presence of strong interactions and large atomic displacements, common models and approximations used for gases and solids do not apply to liquids. For this reason, liquids do not generally fall into any simple classification and have been mostly treated as general many-body systems as a result.

In recent years, a significant effort has been devoted to investigation of various properties of supercritical fluids10,11,12,13,14. This has been an exciting field with a long history since 1822 when Baron Charles Cagniard de la Tour discovered supercritical fluids while conducting experiments involving the discontinuities of the sound in a sealed cannon barrel filled with various fluids at high temperature15. More recently, supercritical fluids have started to be deployed in several important applications, ranging from the extraction of floral fragrance from flowers to applications in food science, such as creating decaffeinated coffee, functional food ingredients, pharmaceuticals, cosmetics, polymers, powders, bio- and functional materials, nano-systems, natural products, biotechnology, fossil and bio-fuels, microelectronics, energy and environment13,14,16. Much of the excitement and interest of the past decade is because of the enormous progress made in increasing the power of relevant experimental tools17,18,19,20. The development of new experimental methods and improvement of existing ones continues to have an important role in this field22,23,24,25,26, with recent research focusing on dynamic properties of fluids27,28,29,30,31,32.

High density and high thermal motion are two main properties responsible for efficient cleaning, dissolving and extracting abilities of supercritical fluids in the above industrial applications. From the point of view of practical applications, supercritical fluids have got the best of both worlds: high density comparable to ordinary liquids and solids, and high thermal motion and diffusivity approaching that of gases. Notably, it is this very combination that presents a formidable problem to the theory: high density and strong interactions mean that theories and approximations used for dilute gases do not apply33. Enskog’s33 and related early kinetic approaches to gases were followed by more extensive developments, yet they do not adequately describe dense systems with strong interactions and many-body correlations, such as supercritical fluids. One general issue with extending gas-like approaches to fluids was noted earlier: in a system with strong interactions, the system energy strongly depends on the type of interactions, and is therefore system-specific, ruling out the possibility to develop a theory that is universally applicable to many fluids, in contrast to gases and solids34.

In addition to theoretical challenges, the lack of fundamental understanding is seen as an obstacle towards wider deployment of supercritical fluids in industrial applications, primarily because of the absence of guidance regarding pressure and temperature at which the desired properties are optimized, as well as the possibility to use new systems13.

In this paper, we focus on the thermodynamic properties of the supercritical state. On the basis of molecular dynamics simulations, we find that specific heat shows a crossover between two different dynamic regimes of the low-temperature rigid liquid and high-temperature non-rigid gas-like fluid. The crossover challenges the currently held belief that no difference can be made between a gas and a liquid above the critical point, and that the supercritical state is homogeneous in terms of physical properties35. We subsequently formulate a theory of system thermodynamics and heat capacity above the crossover. In this theory, energy and heat capacity are governed by the minimal length of the longitudinal mode in the system only, and do not depend on system-specific structure and interactions. We further study the predicted relationship between supercritical exponents of heat capacity and viscosity. A good agreement is demonstrated between calculated and experimental data for noble and molecular supercritical fluids with no free-fitting parameters.

Results

Dynamic crossover of the specific heat

We start with molecular dynamics simulations of a model liquid. Our primary aim here is to show that specific heat, cV, shows a crossover in the supercritical region of the phase diagram. This result is unexpected in view of currently perceived homogeneity of supercritical state in terms of physical properties.

Using molecular dynamics simulations (see Methods), we have simulated the binary Lennard-Jones (LJ) fluid. We have simulated the system with 64,000 atoms using constant-volume (nve) ensemble in the wide temperature range (see Fig. 1) well extending into the supercritical region. Indeed, the temperature range in Fig. 1 is between about 2Tc and 70Tc, where Tc is the critical temperature of Ar, Tc≈150 K; the simulated density, 2,072 kg m−3, corresponds to approximately four times the critical density of Ar. From the energy of the system E at each temperature, we calculate constant-volume specific heat, cV, as (kB=1).

Figure 1: Heat capacity of binary LJ fluid.
figure 1

Calculated cV showing the crossover and continuous dynamical transition around cV≈2, (kB=1). The crossover takes place between different dynamical regimes of the rigid liquid and non-rigid supercritical fluid.

Full size image

We observe that cV decreases steeply from the solid-state value of about 3kB at low temperature to ~2kB around 2,000 K. The steep decrease is followed by crossing over to a considerably weaker temperature dependence. This crossover is a new effect not reported in previous molecular dynamics (MD) simulations. We further observe that the crossover takes place around cV≈2kB. This value of cV=2kB is non-coincidental and corresponds to the crossover taking place across the Frenkel line36,37 as discussed below.

Crossing the Frenkel line corresponds to the qualitative change of atomic dynamics in a liquid. In liquids, atomic motion has two components: a solid-like, quasi-harmonic vibrational motion about equilibrium locations and diffusive gas-like jumps between neighbouring equilibrium positions. As the temperature increases, a particle spends less time vibrating and more time diffusing. Eventually, the solid-like oscillating component of motion disappears; all that remains is the gas-like ballistic motion. That disappearance, a qualitative change in particle dynamics, corresponds to crossing the Frenkel line, the transition of the substance from the liquid dynamics to the gas dynamics. This transition takes place when liquid relaxation time τ (τ is liquid relaxation time, the average time between consecutive atomic jumps at one point in space38) approaches its minimal value, τD, the Debye vibration period. As recently discussed36,37, crossing the Frenkel line is accompanied by qualitative changes of most important properties of the system, including diffusion, viscosity, thermal conductivity and dispersion curves36,37. One of these properties is the crossover of the speed of sound and appearance of the high-frequency ‘fast’ sound20,39,40,41. At the microscopic level, the Frenkel line can be identified by the loss of oscillations in the velocity autocorrelation function37, a point studied in detail elsewhere.

Figure 1 implies that in addition to the dynamic properties discussed above, the thermodynamics of the system changes at the Frenkel line too. This is an important general insight regarding the behaviour of the supercritical state.

The initial steep decrease of cV from about 3kB to 2kB can be quantitatively explained by the progressive loss of two transverse waves with frequency ω>1/τ (refs 42, 43, 44). Physically, this picture is based on Frenkel’s prediction that on timescale shorter than τ, liquid is a solid, and therefore supports rigidity and solid-like transverse waves at short times or at frequency larger than ω>1/τ (ref. 45). When τ approaches its maximal value, τD, the liquid can not sustain transverse waves at any frequency46. Consequently, the potential energy of the system is due to the longitudinal mode only, giving the total energy of 2NT and specific heat of 2kB42,44. Hence, the decrease of cV from about 3kB to 2kB corresponds to the region of a ‘rigid’ liquid, where short-time solid-like rigidity and high-frequency transverse waves exist. On the other hand, the liquid is unable to sustain transverse waves at any available frequency above the Frenkel line. Instead, the liquid enters a new dynamic ‘non-rigid’ gas-like regime, where oscillatory component of particles is lost and the motion becomes purely collisional as in a gas.

We therefore need to develop a new theory capable of describing thermodynamics of supercritical matter above the Frenkel line, where the system enters the new dynamic regime and where cV falls below 2kB and approaches the ideal-gas value of 3/2kB at high temperature (see Fig. 1).

Thermodynamic theory of supercritical state

We now focus on the theory of the non-rigid gas-like liquid above the Frenkel line, and add a new proposal regarding how the system energy can be evaluated. As temperature rises in the ballistic gas-like regime and kinetic energy increases, the mean free path l, the average distance between particle collisions, increases. At the Frenkel line where the ballistic regime starts, l is comparable to interatomic separation. In the limit of high temperature where the particle’s kinetic energy is much larger than potential energy, l tends to infinity as in the non-interacting ideal gas. Our proposal is that l determines the shortest wavelength of the longitudinal mode that exists in the system, λ, because below this length the motion is purely ballistic and therefore can not be oscillatory, λ=l. On the other hand, the longitudinal modes with larger wavelength are supported and they represent the excitations existing in the supercritical system.

We note that the existence of long-wavelength longitudinal waves in a gas (sound) is well known. What is new here is that we propose that the contribution of the longitudinal waves to the energy of the gas-like supercritical system starts from very short wavelengths comparable to interatomic separation a. In this sense, we are extending the solid-state concepts (for example, short-wavelength solid-like phonons with Debye density of states, see below) to the new area of gas-like supercritical state, where these ideas have not been hitherto contemplated. Indeed, it is well established experimentally that dynamics in subcritical liquids shows solid-like character, in that liquids can sustain high-frequency propagating modes down to wavelengths on the atomic scale, with solid-like dispersion relations46. Importantly, recent experimental evidence shows that the same applies to supercritical fluids20,47.

Here and elsewhere, our discussion of liquid vibrational states includes an important point. Namely, a disordered system, such as glass or liquid, supports non-decaying collective excitations obtainable from the secular equation involving the force matrix constructed from the amorphous glass or liquid structure. Harmonic plane waves naturally decay in liquids as in any non-crystalline systems, yet importantly they are clearly seen in fluids experimentally as quasi-linear solid-like dispersion relations even in low-viscous liquids26,46,48, leading to the quadratic density of states g(ω)ω2. A detailed discussion of this point is forthcoming.

In the proposed theory, the energy of the non-rigid supercritical fluid per particle includes the contribution from the kinetic energy, , and the potential energy of the longitudinal phonons with wavelengths larger than l. Using the equipartition theorem , where El is the energy of the longitudinal phonons, we write

where the upper integration limit ω0 is given by the shortest wavelength in the system, λ: , c is the speed of sound, ε(ω, T) is the mean energy of harmonic oscillator, , or ε(T)=kBT is the classical case and Eanh is the anharmonic contribution to the phonon energy.

The second term in equation (1) can be calculated using the Debye density of states, , giving , or , where a is the interatomic separation. The use of the quadratic density of states is supported by the experimental evidence showing solid-like quasi-linear dispersion relationships in supercritical fluids20,47 similar to the subcritical liquids. Eanh can be evaluated in the Grüneisen approximation from the softening of phonon frequencies with temperature, with the result that the energy is modified as , where α is the coefficient of thermal expansion43,44. Then, the energy of non-rigid gas-like supercritical fluid becomes:

We observe that when λ≈a at the Frenkel line, equation (2) gives E=2kBT and the specific heat of 2 (here we neglect the small term αT). This corresponds to the crossover of cV in Fig. 1 to the gas-like regime as discussed above. When λ increases on temperature increase in the gas-like regime, equation (2) predicts that cV tends to 3/2, the ideal-gas value as expected.

Comparison with experimental data

As discussed above, λ is given by the particle mean free path in the non-rigid gas-like regime, and can therefore be calculated from , where η is viscosity, ρ is density and ū is average velocity. Therefore, equation (2) provides an important relationship between the energy of the non-rigid supercritical liquid and its gas-like viscosity. We now check this relationship experimentally by comparing the specific heat predicted by equation (2) and the experimental cV.

We have used the National Institute of Standards and Technology (NIST Chemistry WebBook, http://webbook.nist.gov/chemistry/fluid) database. We aimed to check our theoretical predictions and selected the isochoric data of several supercritical noble and molecular liquids in a wide range of temperature and density. We note that λ in equation (1) changes in response to both temperature and density: λ increases with temperature and decreases with density. Practically, the range of isochoric data in the NIST database is fairly narrow in terms of density compared with the range of temperature. We therefore analyse the temperature behaviour of cv and η along several isochores. Our choice of liquids is dictated by the availability of isochoric data in the supercritical region. For each density, we fit experimental viscosity to η=A0+A1TA2, calculate λ from and subsequently use λ in equation (2) to calculate cV.

We note that Debye model is not a good approximation in molecular liquids, where the frequency of intra-molecular vibrations considerably exceeds the rest of frequencies in the system (3340 K in N2 and 3572 K in CO). However, the intra-molecular modes are not excited in the temperature range of experimental cV (see Figs 2 and 3). Therefore, the contribution of intra-molecular motion to cV is purely rotational, crot. On the other hand, the rotational motion is excited in the considered temperature range, and is therefore classical, giving crot=R for linear molecules in N2 and CO. Consequently, cV for molecular liquids shown in Fig. 3 correspond to heat capacities per molecule, with crot subtracted from the experimental data.

Figure 2: Heat capacity of noble gas liquids.
figure 2

Experimental and calculated cV (kB=1) for noble non-rigid supercritical fluids. Experimental cv and η are taken from the NIST database at different densities as shown in the Figure (a–h). Values of α used in the calculation are 3.5 × 10−3 K−1 (Ne), 3.3 × 10−3 K−1 (Ar), 1.8 × 10−3 K−1 (Kr) and 8.2 × 10−4 K−1 (Xe). The uncertainty of both experimental heat capacities and viscosities is about 2–5%. Insets show viscosity fits.

Full size image
Figure 3: Heat capacity of molecular liquids.
figure 3

Experimental and calculated cV (kB=1) for molecular non-rigid supercritical fluids (a, b). Experimental cv and η are taken from the NIST database. Values of α used in the calculation are 6.5 × 10−3 K−1 (CO) and 8.5 × 10−3 K−1 (N2). The uncertainty of both experimental heat capacities and viscosities is about 2–5%. Insets show viscosity fits.

Full size image

We also note that experimental isochoric cV is affected by λ-like critical anomalies (see Figs 2, 3, 4) because the isochoric NIST data do not extend far above the critical point. Here we do not consider critical effects related to phase transitions, and therefore fit the data at temperatures that are high enough to be affected by the λ-anomaly at the phase transition. In Figs 2, 3, 4 we observe good agreement between experimental and predicted cV, in view of (a) 2–5% uncertainty of experimental cV and η (see NIST Chemistry WebBook), (b) approximations introduced by the Debye model and (c) increased curvature of cV at low temperature because of proximity of λ-anomalies that are not taken into account by the theory. Notably, the agreement is achieved without using free-fitting parameters because ρ, α, ū and a are fixed by system properties. Values of these parameters used in equation (1) are in good agreement with their experimental values.

Figure 4: Heat capacity of liquid Ne.
figure 4

Experimental and calculated cV (kB=1) for non-rigid supercritical fluid Ne. Experimental cv and η are taken from the NIST database at density 700 kg m−3. Value of α used in the calculation is 3.5 × 10−3 K−1.

Full size image

Another way to compare our theory and experimental data is to study supercritical exponents of specific heat and viscosity. Indeed, if η in the non-rigid gas-like supercritical region can be approximated as a power law, ηTγ, then equation (2) makes two predictions. First, E and cV should also follow power laws. Second, equation (2) provides a specific relationship between the scaling exponents of η on one hand and E and cV on the other hand, the relationship that we check below.

Apart from comparing theory and experiments, studying the supercritical scaling exponents is interesting in the wider context of scaling behaviour of physical properties. In the area of phase transitions, the scaling behaviour idea has been a crucial element in the subject of liquids and other systems from the nineteenth century onwards49. Critical points occur in a great variety of systems50,51,52,53, yet there is a considerable degree of similarity in the way in which systems approach the critical point54. Here the calculation of critical exponents serves as one of the main test of the theories.

Once analysed within our theory below, experimental temperature dependence of η, ū, E and cV (see NIST Chemistry WebBook) can be fairly well approximated by the power law above the Frenkel line. Therefore, we use the following power-law relationships:

From (kB=1), λ=const × Tγ−0.5. Using λ in equation (2) gives for the energy (neglecting the small term αT):

and specific heat:

Then, from equations (3), (7) and (8), we find the following relationships between the power-law scaling exponents:

We note that equations (7) and (8), together with the general requirement for cV to be positive and the experimental requirement for cV to decrease with temperature, imply that γ should be in the range . This is the case for the experimental data we examined, as will be shown below. We also note that corresponds to the non-interacting ideal gas and ηT0.5.

We now check equation (10) experimentally. Using the NIST database and equation (6), we calculate ξ from the slope of log(cV=−1.5) vs. , where TF corresponds to the temperature on the Frenkel line and 1.5 is the asymptotic gas value of specific heat. Using equation (10), we calculate the predicted theoretical value γth and compare it with the experimental γexp obtained from fitting the experimental viscosity to ηTγexp.

In Table 1, we show ξ, γexp and γth, averaged over several data sets taken along the isochores at several different densities. We observe the overall good agreement between γexp and γth. We further observe that the power exponent of specific heat, ξ, is close to 0.2 for different systems. In our theory, the similarity of temperature behaviour of cV is due to temperature dependence of λ. This point is discussed in the next section in more detail.

Table 1 Comparison between experiment and theory.
Full size table

Discussion

Our first important observation in this work is that contrary to the current belief, the thermodynamic properties of the supercritical state are not homogeneous. Instead, the specific heat shows a crossover related to the change of particle dynamics, which we attributed to the recently introduced Frenkel line.

We have subsequently focused on thermodynamic properties of supercritical fluids above the Frenkel line. Here, we faced the problem of strong interactions, the long-persisting challenge in condensed matter physics. Indeed, strong interactions imply that approximations used for dilute gases do not apply to real dense liquids33. If we consider realistic strong interactions (assuming that interactions are known and can be represented in analytical form, the assumption that is valid for a relatively small number of simple systems only) and structural correlations that often include those beyond two-body correlations, we quickly find that the problem becomes intractable. Further, strong interactions, coupled with their specificity in different systems, have been suggested to preclude the calculation of energy and heat capacity in general form from the outset34.

In this paper, we addressed the problem in a different way, by substituting all potentially complicated effects of interactions and structural correlations by one physical quantity, the minimal wavelength of the longitudinal mode in the system λ. This has enabled us to rationalize the experimental behaviour of cV as well as to provide the relationship between different physical properties and experimental outcomes (for example, relationship between cV and η). Notably, our approach unveils similarity of thermodynamics of supercritical state in the following sense. First, cV does not explicitly depend on system details such as structure and interactions, but on λ only. Fluids may have very different structure and interactions, yet our theory predicts the similarity of their thermodynamic behaviour as long as λ behaves similarly in those systems. Second, and more specifically, our approach predicts that supercritical scaling of thermodynamic properties, such as heat capacity, is governed by viscosity scaling. Consequently, similar temperature scaling of viscosity gives similar temperature scaling of thermodynamic properties. We note here that we have mostly dealt with systems with fairly simple interatomic interactions, whereas the found similarity of thermodynamic behaviour may not hold in systems with the hierarchy of interactions and non-trivial structural transformations such as water55.

Methods

Molecular dynamics

We have used DL-POLY molecular dynamics simulation code56 to run the system with 64,000 atoms in the constant-volume (nve) ensemble at different temperatures. We have used 3,000 processors of the high-throughput cluster to simulate 500 temperature points in the temperature range of about 0–12,000 K shown in Fig. 1. We have used the model LJ potential57 to simulate model liquids. Each structure was equilibrated for 25 ps, and the system energy and other properties were averaged for the following 25 ps of the production run. The timestep of 0.001 ps was used.

Additional information

How to cite this article: Bolmatov, D. et al. Thermodynamic behaviour of supercritical matter. Nat. Commun. 4:2331 doi: 10.1038/ncomms3331 (2013).