A corresponding-state framework for the structural transition of supercritical fluids across the Widom delta

TL;DR

This paper introduces a novel classification algorithm using radical Voronoi tessellation and percolation theory to define the Widom delta in supercritical fluids, revealing structural transitions and a corresponding state theorem.

Contribution

It presents a new classification method for supercritical fluid phases and establishes a three-parameter corresponding state theorem based on structural analysis.

Findings

The method detects the influence of attractive interactions on structural transitions.

The Widom delta overlaps with response function maxima ridges.

A three-parameter corresponding state theorem is derived.

Abstract

This work proposes a classification algorithm based on the radical Voronoi tessellation to define the Widom delta, supercritical gas-liquid coexistence region, of polyatomic molecules. In specific, we use a weighted mean-field classification method to classify a molecule into either gas-like or liquid-like. Classical percolation theory methods are adopted to understand the generality of the structural transition and to locate the Widom delta. A structural analysis on various supercritical fluids shows that the proposed method detects the influence of the attractive interaction on the structural transition of supercritical fluids. Moreover, we demonstrate that the supercritical gas-liquid coexistence region of water overlap with the ridges of the response function maxima. From the pressure-temperature relation, a three-parameter corresponding state theorem is derived, which states that the fraction of gas-like molecules of a substance is equal to that of another if their reduced pressure, reduced temperature and the critical compressibility factor are the same.

Tables

Table 1. Table 1: Critical constants estimated based on the flat top proposal.

Substances	$T_c$ [K]	$p_c$ [bar]	${\rho }_c$ [$kg/m^3$]	$z_c$
Ar	164.31	63.93	471.30	0.397
${\text{N}}_2$	129.97	38.61	291.91	0.343
${\text{O}}_2$	158.02	58.39	398.78	0.357
${\text{CH}}_4$	199.52	56.55	147.21	0.371
${\text{C}}_2{\text{H}}_6$	312.08	60.09	195.52	0.356
${\text{CH}}_2{\text{OCH}}_2$	470.68	97.16	358.11	0.305
${\text{CO}}_2$	311.99	86.50	452.35	0.324
${\text{NH}}_3$	414.60	132.42	215.09	0.304
${\text{H}}_2\text{O}$	654.10	178.90	284.23	0.208
${\text{CH}}_3\text{OH}$	537.07	89.14	247.32	0.259

Equations

$${\sigma }_{ij}=\frac{{\sigma }_{ii}+{\sigma }_{jj}}{2};{ϵ}_{ij}=\sqrt{{ϵ}_{ii}{ϵ}_{jj}}.$$

$${\rho }_i=\frac{1}{{\sum }_{j=1}^{n_a}{v}_i^j}$$

$$\overline{{\rho }_i}=\frac{1}{N_i+1}\sum _{j\in \mathrm{NN}(i)}\left(\frac{1}{N_j+1}\sum _{k\in \mathrm{NN}(j)}{\rho }_k\right)$$

$${\Pi }_{\mathrm{gas}}=\frac{N_{\mathrm{gas}}}{N_{\mathrm{gas}}+N_{\mathrm{liq}}};{\Pi }_{\mathrm{liq}}=1-{\Pi }_{\mathrm{gas}}$$

$$\{\begin{array}{rl}{\Pi }_{\mathrm{gas}}& =1.0(\mathrm{vapor})\\ {\Pi }_{\mathrm{gas}}& =0.0(\mathrm{liquid})\end{array}$$

$${\rho }_{pa}=\frac{\log (a)+2}{b};{\rho }_{pb}=\frac{\log (a)-2}{b}$$

$$CI=\frac{2}{N_{conf}}\sqrt{p_{inf}{N}_{conf}\left(1-p_{inf}\right)}$$

$$p_{inf}^i={\left[1+\exp (a{\Pi }_i+b)\right]}^c$$

$${\Pi }_i^{av}=\int {\Pi }_i\left(\frac{d{p}_{inf}^i}{d{\Pi }_i}\right)d{\Pi }_i$$

$${\Delta }_i^2=\int ({\Pi }_i-{\Pi }_i^{av}{)}^2\left(\frac{d{p}_{inf}^i}{d{\Pi }_i}\right)d{\Pi }_i$$

$${\Pi }_i^{av}=a{\Delta }_i+{\Pi }_i^{\infty }$$

$$A(\text{mbox}gas)\rightleftarrows A(\text{mbox}liquid)$$

$${\Pi }_{\mathrm{gas}}=1-\frac{1}{1+a\exp (-b{\rho }_r)}$$

$${\Pi }_{\mathrm{gas}}^{pa}=0.09163\pm 0.0014;{\Pi }_{\mathrm{liq}}^{pb}=0.1489\pm 0.0088$$

$$p_r{v}_r=\frac{z{T}_r}{z_c};\frac{p_r}{T_r}=\frac{{\rho }_{r}z}{z_c}$$

$$\frac{p_{r,sgl}}{T_{r,sgl}}=\frac{0.958}{z_c}\text{mbox}asT\to \infty$$

$$p_{r,sgl}=(1.633z_c^{-1}+0.061)(T_{r,sgl}-1)+1$$

Full text

I Introduction

Supercritical fluid, a state of matter beyond the liquid-gas critical point, is regarded as one of the most promising alternatives to conventional liquid solvents used in chemical processes Nalawade, Picchioni, and Janssen (2006); Perrut (2000). One of the useful physical characteristics for the industrial application is its great tunability: small changes in the thermodynamic conditions result in the continuous yet substantial adjustment of thermophysical properties and thermodynamic response functions Vega (2017); Yoon et al. (2017a). It is currently believed that this tunability is deeply related to the inhomogeneous nature of supercritical fluids, refuting the conventional point of view which sees it as a monotonic and homogeneous state of matter.

Thus, a considerable amount of theoretical and experimental works were devoted to understanding the anomalous behavior of supercritical fluids. In scattering and spectroscopic measurements, a large density fluctuation was observed, which imply the existence of the structural inhomogeneities Nishikawa, Tanaka, and Amemiya (1996); Pipich and Schwahn (2018); Tucker (1999). Anomalous behaviors of the response functions and transport properties in the vicinity of the critical point also supported the existence of two different states in the supercritical region Simeoni et al. (2010); Brazhkin et al. (2011); Banuti (2015), where the continuous crossover between these states is reminiscent of subcritical vapor-liquid transition. In order to capture this “supercritical gas-liquid transition”, or the continuous structural transition and the resultant thermodynamic anomalies beyond the critical point, the following thermodynamic methods were suggested.

First, thermodynamic response functions and the correlation length were used to locate the thermodynamic states where the supercritical gas-liquid transition occurs. Simeoni et al. proposed that the Widom line, which was originally defined as a loci of the maximum correlation length, distinguishes gas-like and liquid-like behaviors of supercritical fluids. Simeoni et al. (2010). Banuti et al. suggested that the line of heat capacity maxima, one of the so called pseudo-boiling lines, follows the three-parameter corresponding state principle Banuti, Raju, and Ihme (2017). This approach, which assumes the existence of a single line that demarcates the supercritical fluid region, has also been discussed to understand the anomalous behavior of pure supercritical fluids Imre et al. (2012, 2015); Gallo, Corradini, and Rovere (2014) and their mixtures (supercritical liquid-gas boundary) Raju et al. (2017), and extended to understand the liquid-liquid transition Dey, Roy, and Sarkar (2013); Fuentevilla and Anisimov (2006); Abascal and Vega (2010). However, the response-function based separation of the supercritical region has shown its limits. First, the magnitude of the correlation length and response function maxima rapidly decays and disappears when far from the critical temperature Brazhkin et al. (2011). Second, Widom line and the pseudo-boiling lines only overlap near the critical point; they diverge from each other departing from the critical pointBrazhkin et al. (2011); May and Mausbach (2012); Brazhkin et al. (2014); Fomin et al. (2015). Third, the location of the response function maxima shows a path-dependence Fomin, Tsiok, and Ryzhov (2011); Fomin et al. (2015); Schienbein and Marx (2018). The lines of the response function maxima on isobars and isotherms do not overlap with each other. From these remarks, Schienbein and Marx Schienbein and Marx (2018) recently addressed that no unique separating “line” can be deduced from the response functions.

Independent from these response-function based methods, the structural characteristics of the supercritical fluid have also been studied, which can be dated back to the works of Bernal Bernal (1959, 1964). In the spirit of the Bernal’s work, Finney Finney and Woodcock (2014), Woodcock, and Heyes Heyes and Woodcock (2013) located the percolation transition loci based on the rigidity calculation of the monatomic fluid. The pair correlation function integrals were also used to understand the nature of the inhomogeneous structure of near-critical fluids Tucker (1999); Song, Biswas, and Maroncelli (2000); Skarmoutsos, Dellis, and Samios (2009); Ploetz, Pallewela, and Smith (2017). Several molecular-level classification schemes have also been proposed to define the gas-like and the liquid-like structures in supercritical fluids. Traditionally, the fixed-distance cutoff methods Stillinger Jr (1963); Idrissi et al. (2013); Heyes and Melrose (1988); Škvor and Nezbeda (2009) were used to analyze the inhomogeneous structure of supercritical fluids based on the cluster analysis, but the fixed cutoff methods were susceptible to the selection of the cutoff radius Škvor and Nezbeda (2009). Hence, a series of recent works have used the Voronoi tessellation as a tool to define the local structure of supercritical fluids Idrissi et al. (2010, 2011); Škvor and Nezbeda (2011). For instance, Ovcharov et al. proposed a method for phase identification (MPI) based on the Voronoi tessellation to define liquid, gas, and surface atoms in the two-dimensional system Ovcharov et al. (2017).

Our recent works have focused on the application of the Voronoi tessellation to analyze the structural transition of monatomic supercritical fluids. Implementing machine learning Ha et al. (2018) and probabilistic classification Yoon et al. (2018) strategies, we have successfully identified distinct ‘liquid-like’ and ‘gas-like’ structures coexisting in supercritical fluid, and have located regions where the two microstates coexist, which we termed the ‘Widom delta’. This approach is based on the theory of fluid polyamorphism Anisimov et al. (2018) or the heterogeneous-homogeneous fluctuation Nilsson and Pettersson (2015), where the macroscopic anomalies of fluid phase is understood from the interplay of the two microscopic structures.

This geometric definition of the supercritical gas-liquid coexistence region, called the Widom delta, has shown the following advantages compared to the conventional approaches. First, the boundaries of the Widom delta agree with those obtained from the rigidity calculation Heyes and Woodcock (2013). They are also similar to the coexistence lines detected from the MPI algorithm Ovcharov et al. (2017). Second, neither path-dependence nor termination of the coexistence region was observed. Thus, a contradiction which comes from the path-dependence of the response functions does not occur, and the coexistence region does not vanish even in the deeply supercritical conditions far from the critical point Yoon et al. (2018). Third, the algorithm depends only on the critical density as a classification criterion, which is an intrinsic property of a fluid. Since selecting a proper cutoff distance comes with a certain degree of arbitrariness, a classification scheme without the cutoff distance is meritorious. Lastly, defining a coexistence “region” instead of the “ridges of the response function maxima” for the demarcation of the supercritical fluid phase is adequate to describe the continuous structural variation of supercritical fluids Hestand et al. (2019); no discontinuous (first-order) phase transition occurs above the critical point.

This work aims to extend the classification method, originally demonstrated in the simple monatomic Lennard-Jones fluid system Yoon et al. (2018), to examine the structural transition of polyatomic supercritical fluids and to discover their general structural characteristics. We first compute the fraction of gas-like molecules ($\Pi_{gas}$) by classifying molecules into gas-like and liquid-like applying the geometric procedure based on the radical Voronoi (Laguerre) tessellation. The characteristics of the gas-like and liquid-like structures are examined based on the percolation theory. The classical percolation theory enables us to shed light on the generality of the structural evolution based on the fraction of gas-like molecules. Specifically, we locate the supercritical gas-liquid coexistence region including the supercritical gas-liquid boundary111The term supercritical gas-liquid boundary, which we introduced in previous works, or any other terminologies used by previous works (supercritical liquid-gas boundary Raju et al. (2017), pseudo-boiling line Banuti, Raju, and Ihme (2017) and hypercritical lineBernal (1964)) do not exactly match with the subcritical gas-liquid transition which involves both orientational (topological) and spatial variation of the system. Rather, these terms in addition to ‘gas-like’ or ’liquid-like’ regard the large variation of spatial order observed across the critical density as the remark of the pseudo-boiling., which is defined as the thermodynamic states where $\Pi_{gas}$ becomes 0.5 following the philosophy of the fluid polyamorphism Anisimov et al. (2018), and the percolation transition lines from the finite-size scaling analysis.

We further demonstrate the generality of the supercritical gas-liquid coexistence region on the phase diagrams. The density-temperature diagram shows that the supercritical gas-liquid coexistence region is quasi-universal; the coexistence regions of non-polar and polar substances almost overlap with each other except strongly hydrogen-bonded substances. A comparison between the supercritical gas-liquid coexistence region of water and its lines of the response function maxima reveals that the supercritical gas-liquid coexistence region defined in this work agrees well with most of these demarcation lines. On the pressure-temperature diagrams, we show that the deltoid shape of the supercritical gas-liquid coexistence region is preserved for all substances. From simple relationships based on the compressibility factors, we further demonstrate that supercritical gas-liquid transition of various substances follows the extended corresponding state principle.

II Computational details

II.1 Molecular Dynamics (MD) simulations

Ten substances are selected to investigate the general characteristics of the Widom delta. The interaction parameters for methane, ethane, carbon dioxide, nitrogen, oxygen, and ethylene oxide are modeled using the Transferable Potentials for Phase Equilibria (TraPPE) forcefields Martin and Siepmann (1998); Potoff and Siepmann (2001); Zhang and Siepmann (2006); Ketko et al. (2008). For water and methanol, the potential parameters are adopted from the TIP4P/2005 Abascal and Vega (2005) and OPLS/2016 Gonzalez-Salgado and Vega (2016) forcefields, respectively. The interaction parameters for ammonia are obtained from work by Eckl et al. Eckl, Vrabec, and Hasse (2008). The critical points of these substances are estimated based on the flat top proposal. According to the flat top proposal Heyes and Woodcock (2013), the rigidity and its first derivative of a system become zero at the critical point ($d{\rho}/dp=0$ and $d{\rho^{2}}/dp^{2}=0$).

To calculate the rigidity and its first derivative, we perform NVT simulations Plimpton (1995) in the vicinity of the critical points reported in the earlier works. The number of molecules is 1,000. The timestep is 1 fs. The interatomic potentials are truncated at the cutoff radius of $16.0$ Å and the analytic tail correction term is added. For polyatomic species, the Lorentz-Berthelot mixing rules are used to model the interaction between different atoms:

[TABLE]

The systems are equilibrated for 100,000 steps, and the system pressures are collected every step during 1,000,000 steps. By fitting the cubic equations, the inflection points at different isotherms were obtained. The critical temperature is then calculated as the temperature where $d\rho/dp$ at the inflection point becomes zero. Then, the critical density and pressure are computed. Table 1 shows the critical points of the substances calculated in this manner. After the critical point of each substance is estimated, we perform NVT simulations at $T_{r}=T/T_{c}=1.00-3.70$ and $\rho_{r}=\rho/\rho_{c}=0.2-2.4$ with the same cutoff radius and the number of molecules. The systems are equilibrated for 100,000 steps, and the system configurations are collected every 10,000 steps during the production run (1,000,000 steps). The pressure data are sampled every step. We also simulate supercritical argon at $T_{r}=3.70$ and $\rho_{r}=0.75-1.17$ to examine the system size effect and to carry out the finite-size scaling analysis. The numbers of molecules are $N=1,000,2,000,4,000,8,000,$ and $16,000$. The simulations for the finite-size scaling analysis include 100,000 steps of equilibration run followed by 500,000 steps of data production. The trajectories are collected every 1,000 steps. Besides, we also perform the subcritical NVT simulations of the Lennard-Jones particles ($\sigma=1.0$ and $\epsilon=1.0$) to revalidate the classification strategy. These results can be found in the supplementary material.

II.2 Classification strategy

A probabilistic classification algorithm on the LJ monatomic fluid Yoon et al. (2018) is extended to the polyatomic substances as follows. Since the atomic radii of polyatomic molecules are different, we adopt the radical Voronoi tessellation (Laguerre tessellation) to partition the space into per-atom cells and define the local density of a molecule. In the radical Voronoi tessellation, a configuration of $N$ polydisperse spheres is partitioned into $N$ cells based on the power distance between each other Rycroft (2009). The power distance of a point with respect to a sphere is defined as the squared length of the tangential line from the point to the sphere, hence is dependent on the radius of the sphere (Fig. 1b). The Laguerre cell of a sphere S is then defined as the set of points that are ’closer’ to the S than any other spheres, in terms of the power distance metric. Note that the radical Voronoi tessellation is identical to conventional Voronoi tessellation when the spheres in the system are monodisperse (Fig. 1a).

Thus, a system with $N$ molecules which has $n_{a}$ atoms each has $n_{a}N$ Laguerre cells obtained from radical Voronoi tessellation. We adopt the atomic radii as the LJ distance parameter, $\sigma_{jj}$ for atom type $j$ except water and methanol. Since the atomic radii of hydrogen atoms are given as zero in OPLS/2016 and TIP4P/2005 models, their configurations are partitioned into simple Voronoi cells using the oxygen atoms as the seeds. The local volume of $i^{\rm{th}}$ molecule is defined as the sum of Laguerre cells of its constituting atoms, and its local density is defined as the reciprocal of its local volume:

[TABLE]

The neighbors of the $i^{\rm{th}}$ molecule is defined as the molecules of which atomic Laguerre cells share faces with those of the $i^{\rm{th}}$ molecule, and the second neighbors are defined as the neighbors of neighbors.

After the local density data of the molecules are obtained, the following weighted mean-field strategy is used to determine whether a molecule is gas-like or liquid-like. This weighted mean-field strategy exploits the chemical distance to give a weight on the central particle. When the weighted mean-field local density of the $i^{\rm{th}}$ molecule is supposed to be obtained, it is given as:

[TABLE]

where $N_{m}$ is the number of neighbors of the $m^{th}$ molecule, and $\rm{NN}(\it{m})$ is the set of neighbors of the $m^{\rm{th}}$ molecule. For instance, when counting cells in the two-dimensional systems shown in Fig. 1, the local density of the central atom (dark blue cell) which has ten nearest neighbors is counted eleven times. The local densities of the nearest neighbors are counted four times. The second nearest neighbors are counted either once or twice. If the second nearest neighbor shares two faces with two first nearest neighbors, it is counted twice; otherwise, it is counted once. By averaging the local densities with these weights, one can obtain the weighted mean-field local density of the central atom. The same procedure is applied for the three-diemensional systems. Finally, if this averaged local density is higher than the threshold density, the $i^{\rm{th}}$ molecule is classified as liquid-like. Otherwise, it is gas-like.

After the weighted mean-field classification is carried out, the gas-like fraction of a system is calculated as:

[TABLE]

where $N_{\rm{gas}}$ and $N_{\rm{liq}}$ mean the numbers of gas-like and liquid-like molecules, respectively. This classification result depends on the choice of the threshold density. When various threshold densities are applied to subcritical vapor and liquid phases, we can find that the following result is obtained when the critical density is selected as the threshold density.

[TABLE]

When a threshold density lower than the critical density is used to classify molecules, a considerable amount of vapor molecules are classified as liquid-like ones, and vice versa (see the Supplementary material for the results). Taking into account that the saturated vapor/liquid systems should be clearly distinct from each other, we chose the critical density as the classification density for all substances in this work.

II.3 Percolation analysis

In our last work, we estimated the percolation transition densities $\rho_{pb}$ and $\rho_{pa}$ of the monatomic LJ fluid from the gas-like fraction curve. They were defined as the densities where $\Pi_{\rm{gas}}$ starts to change steeply.

[TABLE]

where $a$ and $b$ are the parameters obtained by fitting a sigmoid equation to the gas-like fraction data. Here, $\rho_{pa}$ was called the percolation transition density of the available volume and $\rho_{pb}$ was the percolation transition density of the bonded clusters (liquid droplets) Heyes and Woodcock (2013). This work estimates $\rho_{pa}$ and $\rho_{pb}$ in a more elaborate manner based on the classical percolation theory. Percolation theory states that a spanning (infinite) cluster connected across the periodic boundary emerge as the occupation probability of a particle increases Stauffer and Aharony (2014). Here, the term “cluster” is not limited to denote the high-density (liquid-like) particles. Instead, it is mathematically defined as a group of neighboring atoms whose identity are same to each other following the classical percolation theory Stauffer and Aharony (2014). We refer to the cluster of gas-like and liquid-like molecules as ‘bubble‘ and ‘droplet’, respectively.

The occupation probability at which an infinite cluster appears is called the percolation threshold. To examine the percolation behavior, we apply the clustering algorithm proposed by Stoll Stoll (1998) to the Voronoi network. In this algorithm, two molecules are regarded to be connected if their classification results are identical and they are the nearest neighbors to each other. Hence, two molecules distant from each other are included in the same cluster if there exists a string of molecules of the same type linking them. The algorithm first counts the independent clusters without regard for the periodic boundary conditions. Second, it tests whether a cluster is infinite or not. The algorithm deems a cluster as an infinite one if one of its molecules at the periphery of the simulation box is connected to another one at the opposite side. After the infinity test, the independent clusters obtained in the first step are connected if they are linked across the periodic boundary conditions.

After the infinity test, we obtain the percolation probability ($p_{inf}$), the probability to find an infinite cluster (gas-like bubble or liquid-like droplet) by dividing the number of configurations which contain at least one infinite cluster by the number of total configurations at the given condition. Since each configuration is independent from each other, the confidence interval (CI) of the percolation probability can be defined from the binomial distribution.

[TABLE]

where $N_{conf}$ is the number of configurations (snapshots). A finite-size scaling analysis is then carried out to determine $\Pi_{thr}^{\infty}$ in the following manner Stauffer and Aharony (2014). We first fit the following sigmoidal equation to $p_{inf}$ data.

[TABLE]

where $i$ is the type of the molecule ($i=\mbox{gas},\mbox{liq}$). Then, the average percolation threshold ($\Pi_{i}^{av}$) and the transition region width ($\Delta$) are computed as:

[TABLE]

The percolation threshold for an infinite system ($\Pi_{i}^{\infty}$) is then obtained as an intercept of the following linear equation.

[TABLE]

where $a$ and $\Pi_{i}^{\infty}$ are fitting parameters.

III Results and Discussion

III.1 Molecular-level classification

Fig. 2 shows the fraction of gas-like molecules at $T_{r}=1.05$ and $T_{r}=3.50$. As observed in our earlier work on the monatomic LJ fluid, the dependence of gas-like fraction on the bulk density could be expressed by an inverse sigmoid function based on the fluid polyamorphism formulation Anisimov et al. (2018); Yoon et al. (2018); Ha et al. (2018). In the theory of fluid polyamorphism, the structural transition between two different amorphous structures is regarded as a reaction.

[TABLE]

The equilibrium constant $K_{eq}$ is then calculated as $K_{eq}=\Pi_{liquid}/\Pi_{\rm{gas}}=\exp(-\beta{\Delta}G^{\ddagger})$ where $\beta$ is the thermodynamic beta ($1/k_{B}T$) and ${\Delta}G^{\ddagger}$ is the Gibbs energy difference between the amorphous structures. By assuming that ${\Delta}G^{\ddagger}$ is proportional to the bulk density, the following equation can be derived.

[TABLE]

At $T_{r}=1.05$, the inverse sigmoid functions of polar substances do not agree with those of non-polar/weakly polar substances (Fig. 2a). $\Pi_{\rm{gas}}$ of polar substances start to decrease at lower densities than those of non-polar ones. When the temperature increases, the inconsistency of $\Pi_{\rm{gas}}$ curves of different substances diminishes (Fig. 2b). This temperature dependence substantiates that the supercritical gas-liquid transition across the Widom delta is largely influenced by the competition between the attractive force and the kinetic energy. Near the critical point, the magnitude of the kinetic energy is comparable to the attraction due to its neighbors. Hence, the magnitude of the intermolecular attraction matters. As the temperature increases, a significant proportion of molecules possess kinetic energy that can overcome the attractive forces among the molecules. As a result, the supercritical gas-liquid transition density increases, and the discrepancy of the gas-like fraction curves among the substances decreases. This competitive relationship between the attractive energy and the kinetic energy can also be understood by examining the compressibility factor ($z{\equiv}p/\rho{k_{B}}T$). When the attractive interaction dominates, $z$ becomes smaller than unity. As shown in Table 1 and the thermodynamic properties given in the supplementary information, the compressibility factors of polar substances at the critical point are lower than the non-polar molecules. As the density increases, they decrease. On the contrary, at higher temperatures far from the critical point ($T_{r}>2$), $z$ becomes larger than one for all materials and increases as the density increases. These results indicate that the influence of the attractive interaction on the structural transition diminishes.

It should also be noted that the $\Pi_{\rm{gas}}$ curves of different species at $T{\sim}T_{c}$ only overlap with each other at high densities. This result indicates that the Weeks-Chandler-Andersen (WCA) liquid perturbation theory Weeks, Chandler, and Andersen (1971) cannot be used to analyze the structural characteristics of low-density and near-critical fluid Ben-Amotz and Stell (2004). As Toxværd stated Toxværd (2015), the equilibrium structure of the gas and the near-critical fluid is strongly affected by the presence of the attractive interaction.

III.2 Structural analysis

We further characterize the structural characteristics based on the local density distributions and the percolation theory. Fig. 3 shows the local density distributions of carbon dioxide and water at $T_{r}=1.05$. The local density distributions of liquid-like molecules are more symmetric than those of gas-like ones. Compared to the local density distributions of carbon dioxide, those of water at the same reduced conditions (Fig. 3d-f) show stronger local density augmentation. Moreover, the distributions of water are broader, which means that more high-density molecules exist as a result of the strong intermolecular interactions ($z_{\rm{CO_{2}}}>z_{\rm{H_{2}O}}$). This result indicates that the classification method proposed in this work also detects the effect of the strong polar forces on the local density augmentation as the integration of the pair correlation function did Yoon et al. (2017b); Song, Biswas, and Maroncelli (2000); Skarmoutsos, Dellis, and Samios (2009).

Fig. 4a-c shows the dependence of $p_{inf}$ and ${\langle}n_{c}\rangle$ on the bulk density at $T_{r}=1.05$. Here, ${\langle}n_{c}\rangle$ is the number of independent gas-like bubbles or liquid-like droplets. At $T_{r}=1.05$, $p_{inf}^{liq}$ of the substances except strongly polar molecules collapse to a single line. $p_{inf}^{liq}$ of water and methanol starts to increase at lower densities than the other fluids. The $p_{inf}^{gas}$ curves show an opposite dependence on the bulk density (see the supplementary material). At high temperatures ($T_{r}\sim 3.50$), $p_{inf}^{liq}$ and $p_{inf}^{gas}$ of all substances collapse to single lines (Fig. 4d). This temperature dependence agrees with the behavior of $\Pi_{\rm{gas}}$ shown in Fig. 2. The steepness of $p_{inf}$ curves became high like the WCA fluids Yoon et al. (2018).

The narrower and higher peaks of ${\langle}n_{c}\rangle$ at high temperatures (Fig. 4e and f) reflect that the thermal agitation prevents the formation of large bubbles and/or droplets far from the critical point. In near-critical fluids, the inhomogeneous fluid structure fluctuates slowly Goodyear and Tucker (1999). As the temperature increases, the thermal fluctuation becomes so vigorous that the inhomogeneous fluid structure cannot remain stable. As a result, the mean bubble (droplet) size decreases and the maximum number of bubbles (droplets) increases. It can also be observed that the maximum of $\langle{n_{c}^{gas}}\rangle$ is always lower than that of $\langle{n_{c}^{liq}}\rangle$. This comes from the fact that the mean local density of liquid-like molecules is higher than that of gas-like ones. Due to the large volume, an infinite gas-like bubbles can be easily formed compared to liquid-like droplets.

The similar dependence of $p_{inf}$ and $\Pi_{\rm{gas}}$ on the temperature indicates that $\Pi_{\rm{gas}}$ can work as an order parameter to describe the supercritical gas-liquid coexistence region; $p_{inf}^{liq}$ and $p_{inf}^{gas}$ curves of all substances at all temperatures collapse to single lines when $\Pi_{\rm{gas}}$ is used as an order parameter (see the supplementary material). Therefore, we conduct the finite-size scaling analysis to calculate the percolation threshold ($\Pi_{i}^{\infty}$).

Fig. 5a demonstrates that no significant finite-size effect exists for $\Pi_{\rm{gas}}$. This result again supports our idea that the fraction of gas-like molecules can be used as a robust parameter to define the supercritical gas-liquid coexistence region. Fig. 5b shows the dependence of $p_{inf}^{gas}$ on $\Pi_{\rm{gas}}$. Unlike $\Pi_{\rm{gas}}$, $p_{inf}^{i}$ largely depends on the system size. As the number of molecules in a system increases, the $p_{inf}^{gas}$ curve becomes like a step function. $p_{inf}^{liq}$ shows a similar behavior as a function of $\Pi_{\rm{liq}}$ (Fig. 5c). By applying Eqn. (10) and Eqn. (11), we can obtain $\Pi_{\rm{gas}}^{\infty}$ and $\Pi_{\rm{liq}}^{\infty}$ as

[TABLE]

The percolation thresholds are comparable to the bond percolation and site percolation thresholds of a pruned Voronoi network Jerauld, Scriven, and Davis (1984). The percolation threshold of liquid-like droplets is larger than that of gas-like bubbles. Again, this result indicates that the percolation of gas-like bubbles can occur at lower concentration of gas-like molecules since the gas-like molecules have large local Voronoi volumes which enable them to be easily connected with each other throughout the entire system.

III.3 Location of the supercritical gas-liquid coexistence region

Since we obtain the percolation thresholds of gas-like bubbles and liquid-like droplets for an infinite system, we can calculate the percolation transition densities $\rho_{pa}$ and $\rho_{pb}$ by solving Eqn. (13) for $\rho$. In addition, we can calculate the supercritical gas-liquid boundary as a set of densities where the gas-like fraction becomes 0.5 from Eqn. (4) ($\rho_{sgl}=\log(a)/b$). Hence, the Widom delta, the supercritical gas-liquid coexistence region which is enclosed with two percolation transition loci and bisected by the supercritical gas-liquid boundary, is completely defined from the gas-like fraction data.

Fig. 6a shows the supercritical gas-liquid coexistence region on the $\rho_{r}$-$T_{r}$ diagram. As expected from the $\Pi_{\rm{gas}}$ curves, the crossover densities ($\rho_{sgl}$, $\rho_{pa}$ and $\rho_{pb}$) of different substances mostly agree with each other except for strongly polar molecules including water, ethylene oxide, ammonia, and methanol near the critical temperature. In the case of methanol and water, $\rho_{pb}$ is lower than that of non-polar substances, whereas $\rho_{pa}$ lines agree with the others. This discrepancy would reflect the asymmetry of the influence of the attractive force on the structure of supercritical fluid. In contrast to dilute gas, the structural characteristics of dense supercritical fluids are mainly determined by the repulsive forces Weeks, Chandler, and Andersen (1971); Toxværd (2015). Since the steepness of the repulsive part of the LJ potential is not relevant to the magnitude of the attractive interaction, $\rho_{pa}$ does not show a significant dependence on the type of the substances. In other words, $\rho_{pa}$ determines the boundary where the structure of supercritical fluid can be interpreted without consideration of the attractive force. On the other hand, $\rho_{pb}$ shows a substance-dependence because the influence of the attractive force becomes significant in the low-density region.

Of special interest is the case of water with hydrogen bonding network. Fig. 6b compares the supercritical gas-liquid coexistence region of water with the ridges of the response function maxima Schienbein and Marx (2018); Strong, Shi, and Skinner (2018), and the hydrogen bonding percolation lines Strong, Shi, and Skinner (2018). The ridge of the compressibility maxima obtained from the E3B3 model Strong, Shi, and Skinner (2018) is included in the supercritical gas-liquid coexistence region. Moreover, the ridges of isothermal compressibility ($\kappa_{T,max}$), isobaric heat capacity ($c_{p,max}$), and thermal expansion coefficient ($\alpha_{p,max}$) maxima on isotherms determined from the IAPWS95 equation of state Wagner and Pruß (2002) are within supercritical gas-liquid coexistence region in spite of the incompleteness of the TIP4P/2005 forcefield Shvab and Sadus (2015). The line of the isochoric heat capacity maxima ($c_{v,max}$), which is known to show an oscillatory behavior due to the limit of the equation of state and terminate at a certain thermodynamic state, deviates from the Widom delta at higher temperatures. Note that hydrogen-bond percolation lines show different behaviors depending on the Hydrogen-bond criteriaStrong, Shi, and Skinner (2018). One of the percolation lines, which only considers the distance $R$ between oxygen atoms of water molecules, is within the transition region, but the others which consider either the intermolecular OH distance or the orientation of water molecules are away from the supercritical gas-liquid coexistence region. This result indicates that the response function maxima and the supercritical gas-liquid coexistence region are mostly relevant to spatial heterogeneity; note that the orientational anomalies discovered in the water systems were observed far below the critical density Hestand et al. (2019); Skarmoutsos, Guardia, and Samios (2017); Ma and Ma (2011).

Next, we explore the pressure-temperature relation. As shown in our earlier work on the Lennard-Jones particle Yoon et al. (2018); Ha et al. (2018), they form a deltoid region regardless of the substances (see the numerical data in the Supplementary material). This linear relation is surprising since the classification algorithm does not consider any thermodynamic variables to classify a molecule; it only exploits the geometrical characteristics of the Voronoi cells and their relationships with neighbors. Fig. 7a shows the supercritical gas-liquid boundaries of various substances. The transition pressures at which $\Pi_{gas}$ becomes 0.5 linearly depend on the temperatures, but their slopes are different from each other. Non-polar molecules have lower slopes than polar or hydrogen-bonded ones. This substance dependence of the transition pressures can be understood based on the following relation. The pressure-temperature relation is formulated as:

[TABLE]

where $z_{c}$ is the critical compressibility factor. By analyzing the limiting behavior of the pressure-temperature relation at high temperature ($T\rightarrow\infty$), we can derive a simple expression from Eqn. (15) to understand the pressure-temperature relation of supercritical gas-liquid boundary. As shown in Fig. 6, the reduced crossover densities ($\rho_{r,sgl}$) converge to 0.958 at high temperature regardless of the substances. Provided that the linear relation between the pressure and the temperature holds at the temperature limit, we can obtain the following equation.

[TABLE]

Fig. 7b shows that this simple relation can be regarded as a good approximation. The slope of the supercritical gas-liquid boundary is well expressed as a linear function of the compressibility factor ($dp_{r,sgl}/dT_{r,sgl}=1.633z_{c}^{-1}+0.061$). Therefore, the supercritical gas-liquid boundary can be expressed as:

[TABLE]

This linear relation also holds for the percolation transition pressures. Hence, this result shows that the fractions of gas-like molecules (gas-likeness) of two substances are equal to each other at the same reduced pressure and temperature if their compressibility factors are equal. Given that the compressibility factor is a function of the acentric factor ($\omega$) Kulinskii (2013), this result also implies that the supercritical gas-liquid transition follows the three-parameter corresponding state principle, which states that all fluids with the same value of $\omega$ (or $z_{c}$) will have the same fraction of gas-like molecules at the same $T_{r}$ and $p_{r}$. That is, supercritical gas-liquid coexistence regions of two different substances on the $p_{r}-T_{r}$ diagram agree with each other if their critical compressibility factors (or acentric factors) are the same.

IV Conclusions

In conclusion, the classification algorithm based on the radical Voronoi tessellation in conjunction with the percolation theory reveals the generality of the Widom delta among simple molecular fluids. The fraction of gas-like molecules does not show significant substance-dependence except for strongly polar or hydrogen bonding molecules near the critical temperature. At $T{\gg}T_{c}$, the percolation transition lines and the supercritical gas-liquid boundary of all substances including the hydrogen bonded ones collapse to single lines. When the percolation probabilities are expressed as functions of $\Pi_{\rm{gas}}$, they collapse into a single line. A finite-size scaling can successfully locate the percolation thresholds of gas-like clusters (bubbles) and liquid-like clusters (droplets). By comparing the supercritical gas-liquid coexistence region of water and the lines of response function maxima, we demonstrate that the Widom delta include most of these lines. This result indicates that the notion of the supercritical gas-liquid coexistence region is suitable to describe the anomalous behavior of supercritical fluids since it does not violate the continuity picture proposed by van der Waals as well as connects the structural characteristics with the thermodynamic anomalies. The linear relationship between the reduced pressure and the reduced temperature observed in the Widom Delta of the monatomic Lennard-Jones system still holds for those of all substances studied in this work. The supercritical gas-liquid boundary and the percolation transition lines follow the three-parameter corresponding state theorem. Noting that the near-critical heat capacity maxima can be framed into the similarity law Banuti, Raju, and Ihme (2017), the generality of the supercritical gas-liquid coexistence region implies that a general two-state model Russo, Akahane, and Tanaka (2018) would be possibly constructed to link the relation between the response functions and the structural characteristics. Hence, the algorithm designed in this work would be helpful to build a structure-property relation to understand the thermophysical behavior of supercritical fluids.

Supplementary material

In supplementary material, we provide the numerical data that can help understand and reproduce the results in the main article. It includes the algorithm validation result; pressure data of near-critical fluids used for the estimation of the critical point; the fractions of gas-like molecules; the probabilities of finding an infinite droplet (bubble); the percolation transition densities (pressures); the number of independent gas-like bubbles and liquid-like droplets in $N=1,000$ systems; and the finite-size scaling data.

Acknowledgements

M.Y.H. and W.B.L. acknowledge the support of National Research Foundation of Korea Grants funded by Korean Government (NRF-2015R1A5A1036133, NRF-2017H1A2A1044355).