Double universality of the transition in the supercritical state

Using molecular dynamics simulations (see the “Methods” section for detail), we have simulated several fluids with different structure and chemical bonding in order to ascertain the effect in a wide range of systems. We simulate molecular (N₂, CO₂), metallic (Pb), hydrogen-bonded network fluid (H₂O) and noble Ar. Supercritical CO₂ and H₂O are particularly important from the industrial point of view due to their deployment in extracting, cleaning, dissolving, environmental and green energy applications Akiya and Savage (2002); Savage (1999); Huelsman and Savage (2013); Kiran et al. (2000). 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_{V}$ on the dynamical length $\lambda_{d}=c\tau$, where $\tau$ is liquid relaxation time J. Frenkel (1955) (for details of calculation and interpretation of $\tau$, see the Methods section) and $c$ is the transverse speed of sound. The specific heat, $c_{V}$ (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 $\lambda_{d}$ is that it sets the upper range of wavelengths of transverse phonons in the liquidlike regime of supercritical dynamics below the Frenkel line (FL) Cockrell et al. (2021a). Details of this mechanism are given in the Methods section. In the gaslike dynamics above the FL, $\lambda_{d}$ corresponds to the particle mean free path and sets the wavelength of the remaining longitudinal mode. This way, $\lambda_{d}$ governs the phase space available to phonons in the system. Since the energy of these phonons contributes to liquid $c_{V}$ Proctor and Maynard-Casely (2020); Proctor (2020); Wallace (1998); Chen (2022); Trachenko and Brazhkin (2016), we predict a unique universal relationship between $c_{V}$ and $\lambda_{d}$ in the supercritical state.

We show the calculated plots of $c_{V}$ on $\lambda_{d}$ in nitrogen, carbon dioxide, lead and water in Figs. 1-4. We set $k_{\rm B}=1$ everywhere in the paper. The variation of $c_{V}$ and $\lambda_{d}$ shown in these Figures corresponds to a very wide range of pressure and temperature. To illustrate this, we also plot several representative paths simulated on the pressure and temperature phase diagram using the argon data Cockrell et al. (2021b) in Fig. 5 and the corresponding plot $c_{V}$ vs $\lambda_{d}$ in Fig. 6.

The dependence of $c_{V}$ on $\lambda_{d}$ 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_{v}$ with respect to $\lambda_{d}$.

The origin of the inversion point is as follows. The dynamical length always has a minimum as a function of temperature when crossing from liquidlike to gaslike regimes of particle dynamics. Recall that this crossover is related to the dynamical crossover at the Frenkel line (FL) Cockrell et al. (2021a). In the liquidlike regime below the FL, particle dynamics combines oscillatory motion around quasi-equilibrium positions and flow-enabling diffusive jumps between these positions J. Frenkel (1955). In this regime, $\tau$ and $\lambda_{d}=c\tau$ decrease with temperature. In the gaslike regime above the FL, the oscillatory component of particle motion is lost, leaving the diffusive jumps only Cockrell et al. (2021a). In this regime, $\lambda_{d}$ becomes the particle mean free path which increases with temperature (see the Methods section for more details). The inversion point is therefore related to the transition between liquidlike and gaslike particle dynamics.

The values of $\lambda_{d}$ and $c_{V}$ at the inversion point are physically significant. The value of $\lambda_{d}=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 $\lambda_{d}$ matches this lengthscale, the fluid stops supporting all transverse phonons simply because the modes with shorter wavelength are non-existent. 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_{V}$ of about 2 in monatomic argon and lead is important too. $c_{V}=2$ corresponds to the loss of the contributions from the two transverse phonon branches, with only the kinetic part $\left(c_{V}=\frac{3}{2}\right)$ and the potential part of the longitudinal mode $\left(c_{V}=\frac{1}{2}\right)$ remaining. Since this loss corresponds to the disappearance of the oscillatory component of particle motion, $c_{V}=2$ is taken as a thermodynamic criterion of the FL Trachenko and Brazhkin (2016); Cockrell et al. (2021a). Phonon anharmonicity can change this result by a relatively small amount Trachenko and Brazhkin (2016), and the disappearance of transverse modes corresponds to $c_{V}=2$ approximately. The inversion point in nitrogen and carbon dioxide corresponds to $c_{V}=2.8-2.9$ (in molecular systems, $c_{V}$ is heat capacity per molecule) due to the additional rotational term contributing $1$ to $c_{V}$. Subtracting $1$ from calculated $c_{V}$, we arrive at $c_{V}=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_{V}(\lambda_{d})$ 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 Figure 4. This is not unexpected, given that water possesses many anomalies which continue to inspire enquiry and research Poole et al. (1992); Gallo et al. (2016). Water’s supercritical state is little understood despite extensive exploitation in industrial and environmental applications Akiya and Savage (2002); Savage (1999); Huelsman and Savage (2013); Kiran et al. (2000). 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 4 to 6, with the associated contribution to entropy and specific heat Eisenberg and Kauzmann (2005). 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_{V}$ at the inversion point varies in the range of about 5-6 per molecule. This higher $c_{V}$ 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 corresponds to $\lambda_{d}$ close to 1 Å as in previous fluids.

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_{V}$ as a function of parameters other than the dynamical length $\lambda_{d}$. In this study we performed simulations along isobars, isotherms, and isochores. The dynamical parameter $\tau$, the relaxation time introduced in the main article, provides a way to compare $c_{V}$ vs $\tau$ along different paths. In Fig. 7 we observe substantial path dependence of $c_{V}(\tau)$ 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_{V}$ or $\tau$. Third, there is no fixed inversion point. This is to be contrasted to the main sequence curves seen in Figures 1-4 and 6, wherein all phase diagram paths share a cross-system universal fixed inversion point at the minimal value of about $\lambda_{d}$ = 1 Å and $c_{V}=2$ ($c_{V}=2$ in monatomic fluids or appropriately modified $c_{V}$ in molecular fluids).

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_{V}$ vs $\lambda_{d}=c\tau$ as seen in Figures 1-4 and 6.

Similarly to $\lambda_{d}$, $\tau$ is a dynamical parameter. However the stark path dependence of $c_{V}(\tau)$ emphasises that the c-transition is not a consequence of simply reducing $c_{V}$ to dynamics, but that the introduction of the special new dynamical parameter $\lambda_{d}=c\tau$ is necessary to achieve a fixed inversion point, data collapse and observe double universality discussed in the next section.

We observe that although there is a moderate path dependence of $c_{V}(\lambda_{d})$ for lead in Figure 3, this path dependence of $c_{V}(\tau)$ in Figure 7 is much more profound. Hence in cases where the $c_{V}(\lambda_{d})$ plot does not achieve the full data collapse, it brings the paths significantly closer together.

We now come to the main finding of this work related to double universality of the c-transition. The first universality is that for each system, the c-transition plot has an inversion point which is fixed and corresponds to about $\lambda_{d}=1$ Å and $c_{V}=2$ ($c_{V}=2$ for monatomic systems or appropriately modified $c_{V}$ in molecular systems) for all paths on the phase diagram, including isobars, isochores and isotherms, spanning orders of magnitude of temperature and pressure. This inversion point provides an unambiguous, theory-independent and path-independent, transition between liquidlike and gaslike states in the sense discussed earlier. The second universality is that this behavior is generic on the supercritical phase diagram and is the same for all fluids simulated.

We now analyse our four systems (which excludes water) on the same set of axes. In order 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_{V}$ as mentioned earlier. This inter-system plot is presented in Fig. 8a.

The four fluids exhibit qualitatively similar main sequence curves: the “c”-shape is present in all curves, and the divergent liquidlike branches converge into almost the same gaslike branch. This plot exhibits what we are calling “double universality”: the function $c_{V}(\lambda_{d})$ across not only different phase diagram paths but also across different fluids converges at the universal inversion point of $c_{V}\approx 2$, $\lambda_{d}\approx 1$ Å. This inversion point therefore constitutes a system independent, path independent, and unambiguous model-free separation between liquidlike and gaslike states in the supercritical state.

To draw the analogy with ordinary phase transitions, we recall the behaviour 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 Guggenheim (1945). In this plot, the specific microscopic details of different systems are often irrelevant to the qualitative behaviour near a phase transition, and the transition falls into a universality class determined by system symmetries and dimensionality Sethna (2006).

The plots in Fig. 8a depict the relationship between a thermodynamic quantity, $c_{V}$, and a dynamical quantity, $\lambda_{d}$. 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 ($\lambda_{d}=1$ Å, $c_{V}=2$).

We also 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.

Finally, 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 Cockrell et al. (2021a). 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 Eckert et al. (1996); Sarbu et al. (2000); Kiran et al. (2000). Our current results therefore give a universal way to locate the inversion point where the performance of a supercritical fluid is optimised, improving the supercritical technologies.

We use DL_POLY molecular dynamics simulations package Todorov et al. (2006). For argon and nitrogen, we use the Lennard-Jones potential fitted to their properties. For nitrogen, we use a rigid two-site Lennard-Jones potential Powles and Gubbins (1976). The potential for carbon dioxide is a rigid-body non-polarisable potential based on a quantum chemistry calculation, with the partial charges derived using the distributed multipole analysis method Gao et al. (2017). 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, gas-liquid coexistence lines and extended to high-pressure and high-temperature conditions Gao et al. (2017). The potential used for water was TIP4P/2005 potential, which is optimised for high pressure and temperature conditions Abascal and Vega (2005). A careful analysis Vega et al. (2009); Vega and Abascal (2011) 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 temperature behaviour, 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 (EAM) potential Belashchenko (2017), which has been used to calculate the properties of molten lead at temperatures up to 25000 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 Xu et al. (2005) of their respective phase diagrams. Equilibration was performed in the NPT ensemble with the Langevin thermostat in order to generate the mean densities along the isobars and isotherms. For argon, system sizes between 500 and 108000 atoms were used with no discrepancy in calculated quantities, consistent with the earlier ascertained insensitivity of viscosity to system size Yeh and Gerhard (2004). 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 $10^{5}$. 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 were run for 1 ns in the NVE ensemble during which all properties were calculated. We calculated $c_{V}$ in the NVE ensemble as Allen and Tildesley (1991):

$$\langle K^{2}\rangle-\langle K\rangle^{2}=\frac{f}{2}NT^{2}\left(1-\frac{f}{2c_{V}}\right)$$

(1)

with $K$ the kinetic energy, and $f$ 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 Green-Kubo theory Zwanzig and Mountain (1965); Balucani and Zoppi (1994):

$$G_{\infty}=\frac{V}{T}\langle\sigma^{xy}(0)^{2}\rangle$$

(2)

$$\eta=\frac{V}{T}\int_{0}^{\infty}\mathrm{d}t\ \langle\sigma^{xy}(t)\sigma^{xy}(0)\rangle$$

(3)

The dynamical length $\lambda_{d}$ was calculated as $\lambda_{d}=c\tau$, where $\tau=\frac{\eta}{G_{\infty}}$, $c^{2}=\frac{G_{\infty}}{\rho}$ and $\rho$ is density.

In this section, we explain the physical origin of the inter-relationship between the specific heat and the dynamical length in the supercritical state. The specific heat, $c_{V}$, 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 Maxwell (1867); J. Frenkel (1955). A liquid has a combined response to shear stress:

$$\derivative{s}{t}=\frac{\sigma}{\eta}+\frac{1}{G_{\infty}}\derivative{\sigma}{t}$$

(4)

where $s$ is the shear strain, $\sigma$ is the shear stress, $\eta$ is the shear viscosity and $G_{\infty}$ is the high-frequency shear modulus. When the external perturbation stops, the internal stress relaxes according to:

$$\sigma(t)=\sigma_{0}\exp\left(-\frac{t}{\tau}\right)$$

(5)

having introduced the Maxwell relaxation time $\tau$:

$$\tau=\frac{\eta}{G_{\infty}}$$

(6)

Frenkel related this time to the average time between molecular rearrangements. This relationship is backed up by experiments and modelling Jakobsen et al. (2012); Iwashita et al. (2013) and has become an accepted view Dyre (2006).

Using Eq. (4), the Navier-Stokes equation can be generalised to include the elastic response of the liquid, yielding Trachenko and Brazhkin (2016):

$$c^{2}\partialderivative[2]{v}{x}=\partialderivative[2]{v}{t}+\frac{1}{\tau}\partialderivative{v}{t}$$

(7)

Seeking the solution of Eq. (7) in the form $v=v_{0}\exp\left(i\left(\omega t-kx\right)\right)$ gives

$$\omega=-\frac{i}{2\tau}\pm\sqrt{c^{2}k^{2}-\frac{1}{4\tau^{2}}}$$

(8)

For $k\leq\frac{1}{2c\tau}$, $\omega$ has no real solutions. For larger $k$, the plane waves decay according to the decay time $\tau$. We therefore define $k_{g}$ as

$$k_{g}=\frac{1}{2c\tau}$$

(9)

Here, we work in terms of the “dynamical length” featuring in Eq. (9), $\lambda_{d}$:

$$\lambda_{d}=c\tau$$

(10)

$\lambda_{d}$ sets the propagation range, or mean free path of transverse phonons, in the liquidlike regime below the FL Cockrell et al. (2021a). This is seen from Eq. (8) which gives the decay factor $\exp\left(-\frac{t}{2\tau}\right)$. Since $\tau$ 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\tau$ 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 gaslike regime of particle dynamics above the FL Cockrell et al. (2021a), $\lambda_{d}$ corresponds to the mean free path of particle motion, $l_{\mathrm{FP}}$ J. Frenkel (1955). Indeed, the shear modulus in a fluid with no interactions is $G_{\infty}=nT$ Zwanzig and Mountain (1965), where $n$ is the concentration. Meanwhile, the gaslike viscosity is Blundell and Blundell (2010) $\eta=\frac{1}{3}\rho v_{\mathrm{th}}l_{\mathrm{FP}}$, where $v_{\mathrm{th}}$ is the thermal velocity and $l_{\mathrm{FP}}$ is the particle mean free path. Noting that $\frac{1}{2}\rho v_{\mathrm{th}}^{2}=\frac{3}{2}nT$, we find $\tau=\frac{\eta}{G_{\infty}}=\frac{l_{\mathrm{FP}}}{v_{\mathrm{th}}}$. The dynamical length in the gaslike state is $\lambda_{d}=c\tau=v_{\mathrm{th}}\tau$ (since the speed of sound in the liquidlike state below the FL, $c$, approximately becomes thermal velocity of particles in the gaslike state above the FL, $v_{\mathrm{th}}$), and $\lambda_{d}=l_{\mathrm{FP}}$.

We can now see the role played by the dynamical length in the liquidlike and gaslike regimes of the supercritical state. $k_{g}$ in Eq. (9) increases with temperature because $\tau$ decreases, and the dynamical length in (10) becomes shorter. This reduces the phase space available for transverse phonons Trachenko and Brazhkin (2016). When $\tau$ approaches its shortest value comparable to the Debye vibration period $\tau_{\rm D}$, $k_{g}$ approaches the Brillouin zone boundary because $c\tau_{\rm D}=a$, where $a$ is the interatomic separation. At this point, all transverse modes disappear, corresponding to $c_{V}=2$ (this value is equal to the kinetic term $\frac{3}{2}$ and the potential energy of the remaining longitudinal mode $\frac{1}{2}$ Trachenko and Brazhkin (2016)). On further temperature increase, the system crosses over to the gaslike regime where the phonon phase space continues to reduce, albeit now for longitudinal phonons. In particular, the longitudinal phonons with wavelengths shorter than $l_{\mathrm{FP}}$ disappear because $l_{\mathrm{FP}}$ sets the shortest wavelength in the system. The associated potential energy of the longitudinal phonons reduces, eventually resulting in $c_{v}=\frac{3}{2}$ as in the ideal gas Trachenko and Brazhkin (2016). Since the phonon energy contributes to liquid $c_{V}$ Proctor and Maynard-Casely (2020); Proctor (2020); Wallace (1998); Chen (2022); Trachenko and Brazhkin (2016), we see that $\lambda_{d}$ directly affects $c_{V}$ because it governs the phonon states in the system.

We are grateful to V. V. Brazhkin and J. Proctor for discussions.