Double universality of the transition in the supercritical state.

Universality aids consistent understanding of physical properties and states of matter where a theory predicts how a property of a phase (solid, liquid, and gas) changes with temperature or pressure. Here, we show that the matter above the critical point has a remarkable double universality. The first universality is the transition between the liquid-like and gas-like states seen in the crossover of the specific heat on the dynamical length with a fixed inversion point. The second universality is the operation of this effect in many supercritical fluids, including N2, CO2, Pb, H2O, and Ar. Despite different structure and chemical bonding, the transition has the same fixed inversion point deep in the supercritical state. This advances our understanding of the supercritical state previously considered to be a featureless area on the phase diagram and a theoretical guide for improved deployment of supercritical fluids in green and environmental applications.

INTRODUCTION

Our view of the phase diagram of ordinary matter is dominated by the three states of solid, liquid, and gas, and the first-order phase transition lines between them, which branch out from the triple point. Of these phase transitions, two have coexistence lines that are finite in length, including the solid-gas sublimation line and the liquid-gas boiling line terminating at the critical point. The matter above the critical point, the supercritical matter, was not thought of as a distinct state of matter and instead seen as a homogeneous state intermediate to liquids and gases and lacking transitions. In particular, distinction between liquid-like and gas-like states within this region was thought to be impossible (, ). Critical anomalies, such as the heat capacity maxima, do not persist far beyond the critical point and furthermore depend on the path taken on the phase diagram (, ). Understanding both the supercritical and liquid states involves several fundamental problems related to dynamical disorder and strong intermolecular interactions (, –). However, such understanding is believed to enhance the deployment of supercritical fluids in important green and environmental applications (, –).

The Frenkel line (FL) separates two qualitatively dynamical regimes of particle motion: combined oscillatory and diffusive motion below the line and purely diffusive above the line (). Practically, the line is calculated from either the dynamical criterion based on the minima of velocity autocorrelation function or the thermodynamic criterion based on the disappearance of transverse modes. This separation of the supercritical state into two different states involves a physical model. It is interesting to ask whether this separation can also be done in a way that is model free. A related question is whether the separation involves universality across all supercritical systems in terms of suitably identified physical parameters.

Here, we show that a deeply supercritical state has a clearly identifiable transition between liquid-like and gas-like states seen in the dependence of the specific heat c on the dynamical length λ, which is doubly universal. The first universality is a fixed path-independent inversion point of the c(λ) crossover, seen as the change of the sign of the derivative of c with respect to λ. The second universality is that the location of the inversion point is similar in all simulated fluids, including supercritical N2, CO2, Pb, Ar, and, to some extent, H2O. Supercritical water has an anomaly, displaying not only similarities but also differences to the other systems in which the universal transition is identified. The inversion point therefore constitutes a system-independent, path-independent, and an unambiguous separation between two physically distinct supercritical states.

RESULTS

Specific heat and dynamical length

Using molecular dynamics (MD) simulations (see Methods for detail), we have simulated several fluids with different structures and chemical bonding to ascertain the effect in a wide range of systems. We simulate molecular (N2 and CO2), metallic (Pb), hydrogen-bonded network fluid (H2O), and noble Ar. Supercritical CO2 and H2O are particularly important from the industrial point of view due to their deployment in extracting, cleaning, and dissolving environmental and green energy applications (, –). We simulate these fluids along several isobars, isotherms, and isochores in the deep supercritical state.

We zero in on the dependence of the specific heat c on the dynamical length λ = cτ, where τ is liquid relaxation time () (for details of calculation and interpretation of τ, see Methods), and c is the transverse speed of sound. The specific heat, c (heat capacity per atom), is an obvious important choice of a thermodynamic quantity because it reflects the degrees of freedom in the system. The role of the dynamical length λ is that it sets the upper range of wavelengths of transverse phonons in the liquid-like regime of supercritical dynamics below the FL (). Details of this mechanism are given in Methods. In the gas-like dynamics above the FL, λ corresponds to the particle mean free path and sets the wavelength of the remaining longitudinal mode. This way, λ governs the phase space available to phonons in the system. Since the energy of these phonons contributes to liquid c (–), we predict a unique universal relationship between c and λ in the supercritical state.

We show the calculated plots of c on λ in nitrogen, carbon dioxide, lead, and water in Figs. 1 to 4. We set kB = 1 everywhere in the paper. The variation of c and λ shown in these figures corresponds to a very wide range of pressure and temperature in the supercritical part of the phase diagram. To illustrate this, we also plot several representative paths simulated on the pressure and temperature phase diagram using the data for argon in Fig. 5 () and the corresponding c versus λ in Fig. 6.

Fig. 1.

Specific heat as a function of the dynamical length in N2.

Specific heat c in the units of kB as a function of the dynamical length λ in supercritical nitrogen across seven phase diagram paths spanning the supercritical state up to 240 times the critical temperature and 3700 times the critical pressure, showing the collapse onto the main sequence with an inversion point at c ≈ 2.9 and λ = 1 Å. Here and elsewhere, kB 1.

Fig. 4.

Specific heat as a function of the dynamical length in H2O.

Specific c as a function of dynamical length λ in supercritical water along eight phase diagram paths spanning the supercritical state up to 15 times the critical temperature and 500 times the critical pressure, showing significant differences in behavior between different phase diagram paths including the location of the inversion point in c and λ.

Fig. 5.

Simulated paths on the phase diagram.

Paths on the phase diagram and simulated pressure and temperature points for argon. Also labeled are the state points where the transition at the inversion point takes place. The triple point and the critical point, together with the boiling and melting lines, are shown.

Fig. 6.

Specific heat as a function of the dynamical length in Ar.

Specific c as a function of the dynamical length, λ cτ across nine paths spanning the supercritical state of argon up to 300 times the critical temperature and 8000 times the critical pressure. All these paths collapse onto the main sequence curve, thereby undergoing a unified dynamic-thermodynamic transition at the path-independent inversion point c ≈ 1.9 and λ = 1 Å.

Specific heat as a function of the dynamical length in CO2.

Specific heat c as a function of the dynamical length λ in supercritical carbon dioxide across six phase diagram paths spanning the supercritical state up to 33 times the critical temperature and 550 times the critical pressure, showing the collapse onto the main sequence with an inversion point at c ≈ 2.9 and λ = 1 Å.

The dependence of c on λ across all simulated paths nearly collapses onto a group of “c”-shaped curves, which we refer to as the main sequence. The main sequence is c-shaped and has an inversion point corresponding to the change of the sign of the derivative of c with respect to λ.

The origin of the inversion point is as follows. The dynamical length always has a minimum as a function of temperature when crossing from liquid-like to gas-like regimes of particle dynamics. Recall that this crossover is related to the dynamical crossover at the FL (). In the liquid-like regime below the FL, particle dynamics combines oscillatory motion around quasi-equilibrium positions and flow-enabling diffusive jumps between these positions (). In this regime, τ and λ = cτ decrease with temperature. In the gas-like regime above the FL, the oscillatory component of particle motion is lost, leaving the diffusive jumps only (). In this regime, λ becomes the particle mean-free path, which increases with temperature (see Methods for more details). The inversion point is therefore related to the transition between liquid-like and gas-like particle dynamics.

The values of λ and c at the inversion point are physically meaningful. The value of λ = 1 Å corresponds to the ultraviolet cutoff, approximately equal to the shortest length scale in the condensed matter system: the interatomic separation set by the length of the chemical bond. When λ matches this length scale, the fluid stops supporting all transverse phonons simply because the modes with shorter wavelength are nonexistent. Concomitantly, particle dynamics can be viewed as the motion with the particle mean-free path approximately equal to the interatomic separation.

The value of c of about 2 in monatomic argon and lead is important, too. c = 2 corresponds to the loss of the contributions from the two transverse phonon branches, with only the kinetic part and the potential part of the longitudinal mode remaining. Since this loss corresponds to the disappearance of the oscillatory component of particle motion, c = 2 is taken as a thermodynamic criterion of the FL (, ). Phonon anharmonicity can change this result by a relatively small amount (), and the disappearance of transverse modes corresponds to c = 2 approximately. The inversion point in nitrogen and carbon dioxide corresponds to c = 2.8 − 2.9 (in molecular systems, c is heat capacity per molecule) due to the additional rotational term contributing 1 to c. Subtracting 1 from calculated c, we arrive at c = 1.8 − 1.9 as in monatomic fluids.

We note that the c plot is not limited in pressure and temperature as long as the system remains chemically unaltered (the same proviso as for the melting line), extending to the entire supercritical state of matter.

All phase diagram paths of different types in argon, nitrogen, and carbon dioxide collapse onto the main sequence curve. As discussed in the next section in more detail, c(λ) along different phase diagram paths follows the same c shape of the main sequence in lead, but moderate path dependence remains far from the inversion point. This could be related to the electronic contribution (not accounted for in the theory based on phonons) represented by the many-body empirical potential in classical MD simulations.

Water, however, shows a different behavior in Fig. 4. This is not unexpected, given that water has many anomalies that continue to inspire inquiry and research (, ). Water’s supercritical state is little understood despite extensive exploitation in industrial and environmental applications (, –). The specific heat of liquid water at the melting point at atmospheric pressure is almost twice as high as that of ice and is related to large “configurational” contribution to the liquid heat capacity. This contribution is related to water-specific hydrogen-bonded network undergoing the coordination change from four to six, with the associated contribution to entropy and specific heat (). This effect precludes the description of water’s heat capacity using phonons only as discussed earlier. Although different paths still result in the c-shaped curves, we see significant path dependence in Fig. 4. c at the inversion point varies in the range of about five to six per molecule. This higher c can be understood as a result of the additional configurational term in water mentioned earlier as well as the rotational term. Nevertheless, the inversion points correspond to λ close to 1 Å as in previous fluids.

Path dependence

The universality of the inversion point and c-shaped main sequence curves observed in the previous section is best taken in the context of the path dependence of c as a function of parameters other than the dynamical length λ. In this study, we performed simulations along isobars, isotherms, and isochores. The dynamical parameter τ, the relaxation time introduced in the main article, provides a way to compare c versus τ along different paths. In Fig. 7, we observe substantial path dependence of c(τ), which manifests in several different ways. The first is that different paths, particularly isochores, have different shapes from one another. The second is that these curves do not coincide at the values of c or τ. Third, there is no fixed inversion point. This is to be contrasted to the main sequence curves seen in Figs. 1 to 4 and 6, wherein all phase diagram paths share a cross-system universal fixed inversion point at the minimal value of about λ = 1 Å and c = 2 (c = 2 in monatomic fluids or appropriately modified c in molecular fluids).

Fig. 7.

Specific heat as a function of relaxation time.

c of simulated Ar, N2, CO2, and Pb as a function of relaxation time τ across different diagram paths.

In summary, we see very different plots depending on which path on the phase diagram is chosen: There is no fixed inversion point, and all curves are far away from each other. This variation is removed once we plot c versus λ = cτ as seen in Figs. 1 to 4 and 6.

Similar to λ, τ is a dynamical parameter. However, the stark path dependence of c(τ) emphasizes that the c transition is not a consequence of simply reducing c to dynamics, but that the introduction of the special new dynamical parameter λ = cτ is necessary to achieve a fixed inversion point and data collapse and to observe double universality discussed in the next section.

We observe that although there is a moderate path dependence of c(λ) for lead in Fig. 3, this path dependence of c(τ) in Fig. 7 is much more profound. Hence, in cases where the c(λ) plot does not achieve the full data collapse, it brings the paths significantly closer together.

Fig. 3.

Specific heat as a function of the dynamical length in Pb.

Specific heat c as a function of dynamical length λ in supercritical lead along eight phase diagram paths spanning the supercritical state up to 20 times the critical pressure and 3400 times the critical temperature, showing the collapse onto the main sequence but with path dependence remaining. The different paths exhibit the same qualitative behavior and share an inversion point of c ≈ 2.0 and λ ≈ 1 Å.

Double universality

We now come to the main finding of this work related to double universality of the c transition. Figures 1 to 6 show the first universality: For each system, the c-transition plot has an inversion point that is fixed and corresponds to about λ = 1 Å and c = 2 (c = 2 for monatomic systems or appropriately modified c in molecular systems) for all paths on the phase diagram, including isobars, isochores, and isotherms, and spanning orders of magnitude of temperature and pressure. This inversion point provides an unambiguous, theory-independent, and path-independent transition between liquid-like and gas-like states in the sense discussed earlier. The second universality seen in Figs. 1 to 6 is that this behavior is generic on the supercritical phase diagram and is the same for all fluids simulated.

To make the second universality more apparent, we analyze our four systems (which excludes water) on the same set of axes. To compare, we must remove the rotational degrees of freedom from the heat capacity of nitrogen and carbon dioxide, which amounts to subtracting 1 from c, as mentioned earlier. This inter-system plot is presented in Fig. 8A.

Fig. 8.

Specific heat as a function of the dynamical length in all systems and universality.

(A) c as a function of the dynamical length, λ, across different phase diagram paths and across four different fluids (Ar, CO2, N2, and Pb), showing the system-independent universal inversion point of c ≈ 2, λ = 1 Å (molecular CO2 and N2 have rotational degrees of freedom removed from its c). (B) The reduced density (ρ is the critical concentration) of coexisting liquids and gases at the boiling line at temperatures close to the critical temperature, T, in a variety of different systems, with all curves coinciding into the same shape [reproduced from ()].

The four fluids exhibit qualitatively similar main sequence curves: The c shape is present in all curves, and the divergent liquid-like branches converge into almost the same gas-like branch. This plot exhibits what we are calling “double universality”: the function c(λ) across not only different phase diagram paths but also different fluids converges at the universal inversion point approximately corresponding to c = 2, λ = 1 Å. The inversion point therefore constitutes a system-independent, path-independent, and unambiguous model-free separation between liquid-like and gas-like states in the supercritical state.

DISCUSSION

To draw the analogy with ordinary phase transitions, we recall the behavior of liquid and gas densities on the coexistence line as the critical point is approached, depicted in Fig. 8B. The experimental relationship between reduced density and reduced temperature of the coexisting liquids and gases is system independent for several small, noble, and molecular elements near the critical point (). In this plot, the specific microscopic details of different systems are often irrelevant to the qualitative behavior near a phase transition, and the transition falls into a universality class determined by system symmetries and dimensionality ().

The plots in Fig. 8A depict the relationship between a thermodynamic quantity, c, and a dynamical quantity, λ. It is in this sense that we consider the c transition to represent a dynamical-thermodynamic transition. The system independence of the main sequence for simple fluids is further suggestive of a universal transition operating in the supercritical state. The fixed point of this dynamical-thermodynamic transition approximately corresponds to (λ = 1 Å, c = 2).

We note that the c transition is not observed in proximity to the critical point. The critical anomalies, caused by diverging correlation lengths, present in this region disrupt the relationship between the dynamical length and the heat capacity in all systems studied here. As mentioned earlier, the inversion point is far above the critical.

We also note that the universal inversion point and the related dynamical transition at the FL corresponds to the solubility maxima (known as “ridges”) and optimal extracting and dissolving abilities of supercritical fluids (). This importantly addresses the widely held belief that improved and more efficient deployment of supercritical fluids will benefit from better theoretical understanding of the supercritical state (, , ). Our current results therefore give a universal way to locate the inversion point where the performance of a supercritical fluid is optimized, improving the supercritical technologies.

In summary, we have shown that the supercritical state has a remarkable double universality. First, the transition between the liquid-like and gas-like states is characterized by fixed inversion point and near path independence. Second, this effect universally applies to many supercritical fluids. This provides new understanding of the supercritical state of matter and a theoretical guide for improved deployment of supercritical fluids in green and environmental applications.

METHODS

Simulation details

We used DL_POLY MD simulations package (). For argon and nitrogen, we used the Lennard-Jones potential fitted to their properties. For nitrogen, we used a rigid two-site Lennard-Jones potential (). The potential for carbon dioxide is a rigid-body nonpolarizable potential based on a quantum chemistry calculation, with the partial charges derived using the distributed multipole analysis method (). The potential was derived and tuned using a large suite of energies from ab initio density functional theory calculations of different molecular clusters and validated against various sets of experimental data including phonon dispersion curves and PVT data. These data included solid, liquid, and gas states, as well as gas-liquid coexistence lines, and extended to high-pressure and high-temperature conditions (). The potential used for water was TIP4P/2005 potential, which is optimized for high-pressure and high-temperature conditions (). A careful analysis (, ) assigned this potential the highest score in terms of the extent to which the results agree with different experimental properties, including the equation of state, high-pressure and high-temperature behavior, and structure. The electrostatic interactions were evaluated using the smooth particle mesh Ewald method in the MD simulations of carbon dioxide and water. The potentials for water and carbon dioxide are rigid body potentials. Simulations of lead were performed using an embedded atom model potential (), which has been used to calculate the properties of molten lead at temperatures up to 25,000 K and 280 GPa, which include the range discussed here.

Systems are simulated along several isobars, isotherms, and isochores in the deep supercritical state, with all paths but named exceptions being far from the critical point and Widom line () of their respective phase diagrams. Equilibration was performed in the NPT ensemble with the Langevin thermostat to generate the mean densities along the isobars and isotherms. For argon, system sizes between 500 and 108,000 atoms were used with no discrepancy in calculated quantities, consistent with the earlier ascertained insensitivity of viscosity to system size (). System sizes of 512 molecules were used for water, nitrogen, and carbon dioxide simulations, and 5120 atoms for lead. The timestep used was 1 fs for water and carbon dioxide and 0.5 fs for lead, which conserved total energy under the Velocity-Verlet integrator in the NVE ensemble to one part in 105. Configurations at the target densities on all paths were then generated, which were then equilibrated with the NVT ensemble for 50 ps. Following this equilibration, we generated 20 independent initial conditions for each state point using seeded velocities, and each of these initial conditions was run for 1 ns in the NVE ensemble, during which all properties were calculated.We calculated c in the NVE ensemble as ()with K being the kinetic energy and f being the number of translational and rotational degrees of freedom available to the molecule in question.

The shear modulus at high frequency and shear viscosity were calculated using the molecular stress autocorrelation function, from the Green-Kubo theory (, )with σ being an off-diagonal component of the microscopic stress tensor. The integration of the long-time tails of autocorrelation functions was implemented using the Green-Kubo formulas (). The 20 independent initial conditions were used to average the autocorrelation function 〈σ(t)σ(0)〉 over these initial conditions. The end result for viscosity was insensitive to adding more initial conditions. The dynamical length λ was calculated as λ = cτ, where , , and ρ is density.

Theory: Specific heat and dynamical length

In this section, we explain the physical origin of the interrelationship between the specific heat and the dynamical length in the supercritical state. The specific heat, c, is an obvious important choice of a thermodynamic quantity because it reflects the degrees of freedom in the system. The dynamical length and its role are discussed below.

The choice of the dynamical parameter is informed by the Maxwell-Frenkel viscoelastic theory (, ). A liquid has a combined response to shear stresswhere s is the shear strain, σ is the shear stress, η is the shear viscosity, and G∞ is the high-frequency shear modulus. When the external perturbation stops, the internal stress relaxes according tohaving introduced the Maxwell relaxation time τ

Frenkel related this time to the average time between molecular rearrangements. This relationship is backed up by experiments and modeling (, ) and has become an accepted view ().

Using Eq. 4, the Navier-Stokes equation can be generalized to include the elastic response of the liquid, yielding ()where v is the transverse velocity field, c is the transverse speed of sound , and ρ is the density.

Seeking the solution of Eq. 7 in the form v = v0 exp (i(ωt − kx)) gives

For , ω has no real solutions. For larger k, the plane waves decay according to the decay time τ. We therefore define k aswhich sets the shortest wave vector for propagating transverse phonons and corresponds to the gap in the phonon momentum space ().

Here, we work in terms of the “dynamical length” featured in Eq. 9, λλ sets the propagation range, or mean-free path of transverse phonons, in the liquid-like regime below the FL (). This is seen from Eq. 8, which gives the decay factor . Since τ sets the time over which the shear stress decays in the liquid as discussed earlier or, in other words, the lifetime of transverse phonons, cτ is a measure of their mean-free path. This implies no phonons with wavelengths longer than the propagation range and is consistent with Eq. 9.

In the gas-like regime of particle dynamics above the FL (), λ corresponds to the mean-free path of particle motion, lFP (). The shear modulus in a fluid with no interactions is G∞ = nT (), where n is the concentration. Meanwhile, the gas-like viscosity is () , where vth is the thermal velocity and lFP is the particle mean-free path. Noting that , we find . The dynamical length in the gas-like state is λ = cτ = vthτ (since the speed of sound in the liquid-like state below the FL, c, approximately becomes thermal velocity of particles in the gas-like state above the FL, vth), and λ = lFP.

We can now see the role played by the dynamical length in the liquid-like and gas-like regimes of the supercritical state. k in Eq. 9 increases with temperature because τ decreases, and the dynamical length in Eq. 10 becomes shorter. This reduces the phase space available for transverse phonons (). When τ approaches its shortest value comparable to the Debye vibration period τD, k approaches the Brillouin zone boundary because cτD = a, where a is the interatomic separation. At this point, all transverse modes disappear, corresponding to c = 2 [this value is equal to the kinetic term and the potential energy of the remaining longitudinal mode ()]. On further temperature increase, the system crosses over to the gas-like regime where the phonon phase space continues to reduce, albeit now for longitudinal phonons. In particular, the longitudinal phonons with wavelengths shorter than lFP disappear because lFP sets the shortest wavelength in the system. The associated potential energy of the longitudinal phonons reduces, eventually resulting in as in the ideal gas (). Since the phonon energy contributes to liquid c (–), we see that λ directly affects c because it governs the phonon states in the system.