# Accurate quantum-corrected cubic equations of state for helium, neon, hydrogen, deuterium and their mixtures

Ailo Aasen, Morten Hammer, Silvia Lasala, Jean-Noël Jaubert, Øivind Wilhelmsen

► To cite this version:

Ailo Aasen, Morten Hammer, Silvia Lasala, Jean-Noël Jaubert, Øivind Wilhelmsen. Accurate quantum-corrected cubic equations of state for helium, neon, hydrogen, deuterium and their mixtures. Fluid Phase Equilibria, 2020, 524, pp.112790. 10.1016/j.fluid.2020.112790. hal-02937126

HAL Id: hal-02937126

https://hal.univ-lorraine.fr/hal-02937126

Submitted on 12 Sep 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Accurate quantum-corrected cubic equations of state for helium, neon, hydrogen, deuterium and their mixtures

<figure><img src="image_3.png" /></figure>

Ailo Aasen<sup>a</sup>, Morten Hammer<sup>a</sup>, Silvia Lasala<sup>b</sup>, Jean-Noël Jaubert<sup>b</sup>, Øivind Wilhelmsen<sup>a,c,*</sup>

## ARTICLE INFO

Article history:
Received 1 June 2020
Revised 8 August 2020
Accepted 19 August 2020
Available online 20 August 2020

Keywords:
Cubic equation of state
Quantum corrections
Hydrogen
Helium
Deuterium
Neon

## ABSTRACT

Cubic equations of state have thus far yielded poor predictions of the thermodynamic properties of quantum fluids such as hydrogen, helium and deuterium at low temperatures. Furthermore, the shape of the optimal $\alpha$ functions of helium and hydrogen have been shown to not decay monotonically as for other fluids. In this work, we derive temperature-dependent quantum corrections for the covolume parameter of cubic equations of state by mapping them onto the excluded volumes predicted by quantum-corrected Mie potentials. Subsequent regression of the Twu $\alpha$ function recovers a near classical behavior with a monotonic decay for most of the temperature range. The quantum corrections result in a significantly better accuracy, especially for caloric properties. While the average deviation of the isochoric heat capacity of liquid hydrogen at saturation exceeds 80% with the present state-of-the-art, the average deviation is 4% with quantum corrections. Average deviations for the saturation pressure are well below 1% for all four fluids. Using Péneloux volume shifts gives average errors in saturation densities that are below 2% for helium and about 1% for hydrogen, deuterium and neon. Parameters are presented for two cubic equations of state: Peng–Robinson and Soave–Redlich–Kwong. The quantum-corrected cubic equations of state are also able to reproduce the vapor-liquid equilibrium of binary mixtures of quantum fluids, and they are the first cubic equations of state that are able to accurately model the vapor-liquid equilibrium of the helium–neon mixture. Similar to the quantum-corrected Mie potentials that were used to develop the covolume corrections, an interaction parameter for the covolume is needed to represent the helium–hydrogen mixture to a high accuracy. The quantum-corrected cubic equation of state paves the way for technological applications of quantum fluids that require models with both high accuracy and computational speed.

© 2020 Published by Elsevier B.V.

## 1. Introduction

Since the first cubic equation of state (EoS) was introduced by van der Waals [1] more than 100 years ago, cubic EoS have remained the preferred choice when computational speed in combination with a reasonable accuracy is needed. They can all be formulated in terms of the pressure

$$P(T, v) = \frac{RT}{v-b} - \frac{a}{(v-r_1b)(v-r_2b)}, \quad (1)$$

where $T$ and $v$ are the temperature and molar volume, $R$ is the gas constant, and $r_1$ and $r_2$ are constants that define the particular EoS. The covolume, $b$, and the attractive parameter, $a$, are functions of temperature and mixture composition.

Cubic EoS are usually preferred for computationally demanding process simulations, optimization studies, and in computational fluid dynamics [2,3]. Two of the most frequently used cubic EoS are Soave–Redlich–Kwong (SRK) and Peng–Robinson (PR), named after the authors of the papers where these EoS were introduced [4,5]. A variety of modifications have contributed to improve the accuracy and widespread use of cubic EoS, for non-polar fluids [6], polar fluids [7], electrolytes [8] and many other applications [9,10]. One type of fluids that so far has eluded the successful modeling by cubic EoS is fluids that exhibit vapor-liquid equilibrium (VLE) at very low temperatures, such as hydrogen and helium. Since hydrogen has the potential to become a widespread environmentally

<footnote><sup>a</sup>SINTEF Energy Research, Trondheim NO-7465, Norway</footnote>

<footnote><sup>b</sup>Université de Lorraine, Ecole Nationale Supérieure des Industries Chimiques, Laboratoire Ré et Génie des Procédés (UMR CNRS 7274), 1 rue Grandville, Nancy 54000, France</footnote>

<footnote><sup>c</sup>Norwegian University of Science and Technology, Department of Energy and Process Engineering, Kolbjørn Hejes vei 1B, Trondheim NO-7491, Norway</footnote>

<footnote>* Corresponding author at: Norwegian University of Science and Technology, Department of Energy and Process Engineering, Kolbjørn Hejes vei 1B, NO-7491 Trondheim, Norway.</footnote>

<footnote>E-mail address: oivind.wilhelmsen@ntnu.no (Ø. Wilhelmsen).</footnote>
friendly energy carrier in the future, the interest in such fluids is increasing [11]. For instance, liquefaction is being investigated as an alternative for large-scale distribution of hydrogen across long distances [12]. Mixtures consisting of helium, neon and hydrogen have been touted as promising refrigerants with the potential to significantly improve the energy efficiency of the hydrogen liquefaction process [12,13]. Unfortunately, even with advanced mixing rules [14], cubic EoS have been unable to reproduce the VLE of helium-neon [12]. The reason is that these fluids are strongly influenced by quantum-mechanical effects, as quantified by the de Broglie thermal wavelength $\lambda_B$ [15,16]:

$$ \lambda_B = \frac{h}{\sqrt{2\pi mk_B T}} \qquad (2) $$

where h is Planck's constant, m is the particle mass, and $k_B$ is Boltzmann's constant. When $\lambda_B$ becomes comparable to the typical intermolecular separation, quantum effects can influence the fluid properties significantly.

Recently, the statistical associating fluid theory of quantum-corrected Mie potentials (SAFT-VRQ Mie) [17,18] was presented. This EoS models quantum effects by adding the so-called Feynman-Hibbs-corrections to a Mie potential, where the Mie potential represents the classical behavior of the fluid [19]. Excellent agreement with experiments was obtained for the VLE of helium-neon, except in the critical region due to the incomplete development of thermodynamic perturbation theory for mixtures [20]. However, molecular simulations with the underlying interaction potentials of the EoS were in excellent agreement with the experiments, also in the critical region. The interaction potentials revealed that a crucial effect to account for is “quantum swelling”, namely that the effective size of the particles increases at low temperatures due to the wave-like nature of particles that becomes more apparent as $\lambda_B$ increases. In cubic EoS, the size of the particles is represented by the covolume b in Eq. (1), which is usually assumed to be temperature-independent and therefore does not account for quantum swelling. In this work, we shall derive quantum corrections for the covolume parameter of cubic EoS based on the quantum swelling of the Feynman-Hibbs-corrected Mie potentials that were successful in previous work.

The shortcomings of cubic EoS for describing fluids that exhibit quantum effects has for long been known. In his celebrated work that formed the basis for the SRK EoS [4], Soave pointed out that less accurate results were obtained for mixtures with hydrogen. Some efforts have been directed towards improving cubic EoS for such applications. Thomas et al. [21] evaluated the applicability of several EoS for modeling pure helium, demonstrating that PR in general gives inaccurate predictions at low temperatures, except for pressures below 5 bar. Kaviani et al. [22] made a cubic EoS for pure helium by fitting five-parameter correlations for both the $\alpha$ function and the covolume, which were fitted to thermodynamic properties below 15 K and up to 16 bar. Although high accuracy was obtained at low temperatures, their EoS is not applicable beyond the range where it was regressed.

Tkaczuk et al. [23] recently developed a Helmholtz energy mixture model for the helium-neon mixture. Their model relies on multiparameter pure fluid EoS combined with multiparameter mixture models. Tkaczuk et al. use more than 40 binary interaction parameters to regress the properties of the binary helium-neon mixture model. Their representation of fluid properties is therefore expected to be much more computationally demanding than the description provided in this work.

The attractive parameter in Eq. (1), a, is usually assumed to depend on temperature, and is written as the product of its value at the critical point, $a_c$ and the $\alpha$ function, as: $a(T) = a_c\alpha(T)$. Which functional form to use for the $\alpha$ function has been widely debated. Recently, Le Guennec et al. [24] introduced consistency criteria for

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Pair interaction potential for the Lennard-Jones (LJ) potential (black, solid line) and the LJ potential with second order Feynman-Hibbs (FH) quantum corrections (red, dashed). The temperature is given by $k_B T/\epsilon = 1$, and the particle mass is that of helium. The FH correction increases the effective diameter $\sigma_{eff}$, and decreases the well depth $\epsilon_{eff}$, compared to the classical LJ parameters, $\sigma$ and $\epsilon$. The inset shows that the effective range of the potential increases slightly.</figcaption></figure>

the $\alpha$ functions of cubic EoS. They showed that by enforcing these criteria, the predictions in the supercritical domain can be significantly improved, with a negligible loss of accuracy in the subcritical domain [25]. One of the consistency criteria was that the $\alpha$ function should exhibit a monotonic decay as a function of temperature. Following these recommendations, cubic EoS have been reparametrized for thousands of components using consistent alpha functions fitted to saturation data, including for helium, neon, hydrogen and deuterium [6,26]. Le Guennec et al. [24] showed that the $\alpha$ function of hydrogen and helium did not exhibit a monotonic decay with temperature, and did thus not follow the consistency criteria derived for other fluids. The reason for this behavior has remained unexplained, although it has been linked to their accentric factors being negative [6]. We will show that a near classical behavior with a monotonic decay can be obtained for the $\alpha$ function by incorporating temperature-dependent quantum corrections for the covolume parameter that explicitly accounts for quantum swelling.

The physical validity of a temperature-dependent covolume has been questioned by Kalikhman et al. [27], who demonstrated that any repulsive term involving temperature-dependent covolumes results in a negative infinite value for the isochoric heat capacity at infinite pressure. We will show that the quantum corrections derived in this work will introduce negative heat capacities in agreement with the derivations by Kalikhman, but only at pressures well outside regions of industrial applications of the substance. We also verify that thermodynamic stability criteria are fulfilled and that pressure-volume isotherms do not cross close to the critical point for the covolume corrections derived in this work. Parameters for the quantum fluids will be provided for the most frequently used cubic EoS, SRK and PR, but can easily be extended to other cubic EoS. The quantum corrections yield a vast improvement in the accuracy of cubic EoS, both for single-component fluids and mixtures. The VLE of all examined fluid mixtures are reproduced to a high accuracy without the need for advanced mixing rules [28].

## 2. Theory

We showed in previous work that Feynman-Hibbs-corrected Mie potentials are capable of capturing the thermodynamic properties of fluids with quantum effects to a relatively high accuracy [17,18]. As shown in Fig. 1, the quantum corrections influence a classical interaction potential mainly through two effects:

1. They increase the distance where the classical interaction potential is zero, $\sigma$, to a larger, temperature-dependent value, $\sigma_{eff}(T)$.
2. They reduce the well-depth, $\epsilon$ of the classical interaction potential to a smaller, temperature-dependent value $\epsilon_{\text{eff}}(T)$.

Numerous studies have shown that cubic EoS are capable of representing the thermodynamic properties of classical fluids to a reasonable accuracy. The main hypotheses of this work are (1) that the quantum swelling can be incorporated into the cubic EoS through a temperature-dependent covolume correction, and (2) that the reduction in the well-depth can be captured by a suitable modification of the $\alpha$ function. In the following we will derive functional forms suitable to describe these quantum corrections based on the Feynman-Hibbs-corrected Mie potentials described in Refs. [17,18].

## 2.1. The covolume of cubic EoS

Van der Waals derived his famous equation of state under the assumption that the excluded volume per particle, known as the covolume, is four times the particle volume:

$$b' = \frac{1}{2} \left( \frac{4\pi\sigma^3}{3} \right), \qquad (3)$$

where $\sigma$ is the particle diameter. The covolume, $b$, used in the equation of state can be found by multiplying $b'$ with Avogadro's number,

$$b = N_A b'. \qquad (4)$$

From the reduced form of the van der Waals equation of state, the covolume can be related to the critical properties of the fluid by imposing that $(dP/dV)_{T_c} = 0$ and $(d^2P/dV^2)_{T_c} = 0$, where $T_c$ is the critical temperature. This gives the following expressions:

$$b_{vdW} = 0.125 \frac{T_{\text{crit}} R}{P_{\text{crit}}}, \qquad (5)$$

$$b_{\text{SRK}} = 0.08664 \frac{T_{\text{crit}} R}{P_{\text{crit}}}, \qquad (6)$$

$$b_{\text{PR}} = 0.07780 \frac{T_{\text{crit}} R}{P_{\text{crit}}}, \qquad (7)$$

for the covolumes of the van der Waals, the SRK and the PR equations of state respectively.

## 2.2. Relating the covolume of cubic EoS to interaction potentials of classical fluids

The intermolecular potential between two nearly spherical molecules can be approximated by a Mie potential:

$$u_{\text{Mie}}(r) = C \epsilon \left[ \left(\frac{\sigma}{r}\right)^{\lambda_r} - \left(\frac{\sigma}{r}\right)^{\lambda_a} \right], \qquad (8)$$

where $\epsilon$ is the well depth, $\sigma$ is the finite distance at which the potential is zero, $\lambda_a$ and $\lambda_r$ are the attractive and repulsive exponents, and

$$C = \left( \frac{\lambda_r}{\lambda_r - \lambda_a} \right) \left( \frac{\lambda_r}{\lambda_a} \right)^{\frac{\lambda_a}{\lambda_r - \lambda_a}}. \qquad (9)$$

For impenetrable spheres as represented by the hard sphere potential, the excluded volume is given uniquely as the volume of the spheres. The definition of "excluded volume" corresponding to the repulsive part of the Mie potential is not unique, as the finite opposing force of the potential at radii smaller than $\sigma$ allows particles to be at intermolecular distances smaller than $\sigma$. However, as the repulsive part of the interaction potential is usually very steep, the excluded volume can be approximated by the positive part of

<figure><img src="image_5.png" /><figcaption>Fig. 2. The covolume $b$ of SRK and PR divided by the excluded volume of the Mie potential models $b_{\text{Mie}}$ in Ref. [30], plotted as a function of molecular weight for methane, argon, krypton and xenon.</figcaption></figure>

the potential. For Mie potentials, the particle diameter is then simply $\sigma$ in Eq. (8), where the interaction potential goes from positive to negative. Fig. 2 shows that the excluded volume of the Mie potential, given by $b_{\text{Mie}} = \pi\sigma^3/6$ is proportional to the covolume parameter $b$ of cubic EoS for classical fluids for both SRK and PR. This observation, and the fact that the proportionality factor is close to 0.5, has been made before [29].

## 2.3. Interaction potentials and quantum corrections

Assuming that the Mie potential can be used to describe classical fluids, Feynman-Hibbs corrections can be added to the classical potential to represent quantum effects,

$$u_{\text{MieFH}}(r)/(C\epsilon) = \left(\frac{\sigma}{r}\right)^{\lambda_r} - \left(\frac{\sigma}{r}\right)^{\lambda_a} + \frac{D}{T} \left( Q_1(\lambda_r) \frac{\sigma^{\lambda_r}}{r^{\lambda_r+2}} - Q_1(\lambda_a) \frac{\sigma^{\lambda_a}}{r^{\lambda_a+2}} \right) + \frac{D^2}{T^2} \left( Q_2(\lambda_r) \frac{\sigma^{\lambda_r}}{r^{\lambda_r+4}} - Q_2(\lambda_a) \frac{\sigma^{\lambda_a}}{r^{\lambda_a+4}} \right). \qquad (10)$$

where the term that has prefactor $D$ is the first-order, and the term with prefactor $D^2$ is the second-order Feynman-Hibbs quantum correction, and

$$Q_1(\lambda) = \lambda(\lambda - 1), \qquad (11)$$

$$Q_2(\lambda) = \frac{1}{2}(\lambda + 2)(\lambda + 1)\lambda(\lambda - 1), \qquad (12)$$

$$D = \frac{\hbar^2}{12mk_B}. \qquad (13)$$

For Mie potentials with Feynman-Hibbs quantum corrections, the repulsive region increases in size at lower temperatures (cf. Fig. 1). This can be quantified by the effective particle diameter, $\sigma_{\text{eff}}$, given by the separation where the quantum-corrected potential equals zero. The effective diameter was found by solving

$$u_{\text{MieFH}}(r; T) = 0 \qquad (14)$$

for $r$ using a Newton-Raphson solver with analytical derivatives and the value of $\sigma$ as initial guess.

## 2.4. The quantum-corrected covolume parameter

We will relate the quantum covolume $b_q(T)$ to the classical covolume $b$ by

$$b_q(T) = b\beta(T), \qquad \beta(T) = \left( \frac{\sigma_{\text{eff}}(T)}{\sigma_{\text{eff}}(T_c)} \right)^3. \qquad (15)$$
Eq. (15) is based on the assumption that the covolume parameter b of a classical cubic EoS is correct at the critical point, as it is determined from the measured critical temperature, T_c, and pressure, P_c. The quantum correction, β(T) captures the quantum swelling of the particle volume with respect to its volume at the critical point. Since β(T) is given from the quantum-corrected Mie potential, it should thus be independent of the choice of cubic EoS. This will be studied in further detail in Section 4.

Clearly $\sigma_{\text{eff}}(T^{-1} = 0) = \sigma$; in other words, the classical value is recovered at high temperatures. Differentiating $\sigma_{\text{eff}}$ as defined by Eq. (14) with respect to $T^{-1}$ gives that the following slope should be satisfied in the limit $T^{-1} \to 0$:

$$ \left. \frac{\partial \sigma_{\text{eff}}}{\partial T^{-1}} \right|_{T^{-1} \to 0} = \frac{D(Q_1(\lambda_r) - Q_1(\lambda_a))}{\sigma(\lambda_r - \lambda_a)}. \quad (16) $$

Eq. (16) holds for both the first and the second order quantum-corrected Mie potentials. In the limit $T^{-1} \to \infty$ however, the first order quantum-corrected Mie potential reaches the maximum value

$$ \sigma_{\text{eff}}^{\text{max}} = \sigma \left( \frac{Q_1(\lambda_r)}{Q_1(\lambda_a)} \right)^{1/(\lambda_r - \lambda_a)}. \quad (17) $$

The following expression can satisfy all these constraints:

$$ \sigma_{\text{eff}} = \sigma + \frac{D}{T + T_{\text{adj}}} \frac{(Q_1(\lambda_r) - Q_1(\lambda_a))}{\sigma(\lambda_r - \lambda_a)}, \quad (18) $$

where

$$ T_{\text{adj}} = c_{FH} \frac{D(Q_1(\lambda_r) - Q_1(\lambda_a))}{\sigma(\sigma_{\text{eff}}^{\text{max}} - 1)(\lambda_r - \lambda_a)}. \quad (19) $$

The reason for the prefactor $c_{FH} = 1.4$ for first order-corrections and $c_{FH} = 0.5$ for the second order corrections, is to compensate for the inability of the linear extrapolation in $T^{-1}$ at higher temperatures. The above expression has the right derivative in the high temperature limit, but underpredicts the maximum effective diameter slightly. However, since the limit $1/T \to \infty$ can never be reached in practice due to solid formation, this is of little practical relevance.

Eq. (18) is plotted together with the increase in the effective diameter from the first and the second order corrected Mie potentials with parameters for hydrogen, deuterium, neon and helium in Fig. 3, where the curves are from Eq. (18) and the circles are from the quantum-corrected Mie potentials. The figure displays excellent agreement between the derived expression and the exact value for $\sigma_{\text{eff}}/\sigma$ obtained by solving Eq. (14), in particular at supercritical temperatures.

The derived $\beta$ function can be written as

$$ \beta(T) = \left(1 + \frac{A}{T + B}\right)^3 / \left(1 + \frac{A}{T_c + B}\right)^3, \quad (20) $$

where A and B are component-specific parameters. Since cubic EoS are not based on Mie fluids, and since the Feynman-Hibbs corrections become inaccurate at very low temperatures, we will also evaluate the empirical correlation given by Eq. (20) where A and B are fitted to experimental data.

## 2.5. The $\alpha$ function and quantum effects

In addition to the hypothesis that quantum swelling can be captured by Eq. (15), we also hypothesize that the reduction of the well-depth given by the quantum corrections (see Fig. 1) can be modeled by properly modifying the $\alpha$ function of cubic EoS. Since the quantum-corrected covolume is larger at subcritical temperatures, the attraction and thus the $\alpha$ function must be reduced, e.g. to recover the same saturation pressures. This is because the two

<caption>**Table 1**</caption>

Triple point temperature $T_t$, critical temperature $T_c$ and critical pressure $P_c$ for the four components considered. For helium $T_t$ represents the superfluid transition.

<table><thead><tr><th></th><th>H<sub>2</sub> [35]</th><th><sup>4</sup>He [36]</th><th>Ne [37]</th><th>D<sub>2</sub> [38]</th></tr></thead><tbody><tr><td>T<sub>t</sub> (K)</td><td>13.957</td><td>2.17</td><td>24.556</td><td>18.724</td></tr><tr><td>T<sub>c</sub> (K)</td><td>33.19</td><td>5.1953</td><td>44.492</td><td>38.34</td></tr><tr><td>P<sub>c</sub> (bar)</td><td>12.964</td><td>2.276</td><td>26.79</td><td>16.796</td></tr></tbody></table>

terms of cubic EoS (cf. Eq. (1)) are regressed together to approximate the experimental data. For supercritical temperatures, the co-volume is smaller than the classical covolume, and thus the quantum $\alpha$ function should be smaller than the classical one. These effects will change the “bell-shaped” $\alpha$ function found previously for hydrogen and helium for the classical PR. However, it is not guaranteed that the change will make the $\alpha$ function consistent in the classical sense [24], and not even monotonically decreasing. We will explore this further by comparing to experimental data in Section 4.1.

Throughout this work, the cubic $\alpha$ function is correlated using the three-parameter Twu function [31]:

$$ \alpha(T) = T_r^{N(M-1)} \exp[L(1 - T_r^{MN})]. \quad (21) $$

Eq. (21) is designed such that $\alpha(T_c) = 1$ for all (L, M, N). Having three parameters, we assume that this functional form is flexible enough to capture the quantum effects on the attractive term. In Section 4.1 we will investigate whether the optimal $\alpha$ function becomes a monotonic function when incorporating a quantum-corrected covolume.

## 2.6. Mixing and combining rules

For $a$ we apply the usual quadratic mixing rules:

$$ a = \sum_{i} \sum_{j} x_{i}x_{j}a_{ij}, \quad a_{ij} = \sqrt{a_{ii}a_{jj}(1 - k_{ij})}, \quad (22) $$

where $x_i$ and $x_j$ are the mole fractions of components $i$ and $j$, and $k_{ij}$ is a binary interaction parameter.

The mixture covolume is calculated as

$$ b = \sum_{i} \sum_{j} x_{i}x_{j}b_{ij}, \quad b_{ij} = \frac{b_{ii} + b_{jj}}{2}(1 - l_{ij}). \quad (23) $$

A common choice is to set $l_{ij} = 0$, which converts the mixture co-volume into a linear combination of the single-component covolumes.

In addition to the $\alpha$ and $\beta$ functions, for each component we also fitted a (temperature-independent) Péneloux volume shift $c$ [32]. Some properties such as enthalpies are affected by such a volume shift [33], and care should be taken to also shift these properties correctly. The impact of volume shifts on mixtures is more complicated [34], but one can show that a linear mixing rule of Péneloux parameters will not impact the Pxy phase envelopes. Such a linear mixing rule is adopted in this work:

$$ c = \sum_{i} x_{i}c_{ii}. \quad (24) $$

# 3. Parameter estimation methodology

In this work, we have employed the common method of first fitting single-component parameters, and in a subsequent step fitting the binary interaction parameters. Table 1 lists the triple point temperature, critical temperature, and critical pressure for the considered components.
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** The relative increase in effective diameters for helium, hydrogen, neon and deuterium as computed exactly from the Mie-FH potentials (circles) and the analytic approximation Eq. (18) (curves).</figcaption></figure>

<table><caption>**Table 2**
Weights in the objective function for the fitting pure fluid parameters. Square brackets indicate an interval within which the weight varied among the different components.</caption><thead><tr><th colspan="6">Saturation</th></tr></thead><tbody><tr><th>Property</th><th>p<sup>sat</sup></th><th>ρ<sup>sat</sup></th><th>c<sub>i</sub><sup>sat</sup></th><th>c<sub>p</sub><sup>sat</sup></th><th>Δh<sup>sat</sup></th></tr><tr><th>Weight</th><td>1.0</td><td>0.5</td><td>[0,0.5]</td><td>[0,0.5]</td><td>0.5</td></tr></tbody><thead><tr><th colspan="6">Single-phase</th></tr></thead><tbody><tr><th>Property</th><th>ρ<sup>sup</sup></th><th>c<sub>i</sub><sup>sup</sup></th><th>c<sub>p</sub><sup>sup</sup></th><th>h<sup>sup</sup></th><th>w<sup>sup</sup></th></tr><tr><th>Weight</th><td>0.5</td><td>0.5</td><td>0.5</td><td>0.5</td><td>[0,0.1]</td></tr></tbody></table>

## 3.1. Objective functions

### 3.1.1. Pure fluids

The parameters (L, M, N) for the Twu $\alpha$ function (Eq. (21)), the Péneloux volume shift $c$, and optionally the parameters (A, B) in the covolume beta function (Eq. (20)), were fitted to pseudo-experimental data. These data were generated using the highly accurate reference equations of state for helium [36], neon [37], hydrogen [35], and deuterium [38], which were accessed through CoolProp [39].

The objective function is given by a weighted sum of relative deviations:

$$O^{\text{pure}}(L, M, N; A, B; c) = \sum_{i \in \{\text{states}\}} \sum_{\xi \in \{\text{properties}\}} w_{\xi} \frac{|\xi_{i}^{\text{expt}} - \xi_{i}^{\text{calc}}|}{|\xi_{i}^{\text{expt}}|}, \quad (25)$$

where $i$ indexes the states where the deviations are calculated, $\xi$ represents one of the ten thermodynamic properties considered, and $w_{\xi}$ is the weight given to the property $\xi$. Table 2 lists the properties and weights used in the objective function. We used a simplex algorithm [40] for the minimization.

For the saturation properties, we chose 20 states equispaced in temperature. The upper temperature limit was 4.8 K for helium, 30.0 K for hydrogen, 34.5 K for deuterium and 41.0 K for neon. The

lower temperature limits for neon and deuterium coincide with their triple point temperatures. Currently there are no technical applications of hydrogen at temperatures below the boiling point, 20.3 K; to avoid sacrificing accuracy at higher temperatures, the lower temperature for hydrogen was chosen to be the same as for deuterium. The lower limit of helium was chosen as 3 K, which is above the superfluid transition at 2.2 K but still well below its boiling point temperature of 4.2 K.

The single-phase properties were evaluated over a rectangular grid in temperature-pressure space, with five supercritical temperatures and 20 pressures from 1 bar and up. For neon, hydrogen and deuterium the temperatures were 50, 100, 150, 200 and 300 K, and the maximum pressure was 500 bar. For helium the temperatures were 20, 30, 40, 100 and 150 K and the maximum pressure was 300 bar.

Residual properties of pure fluids, such as the residual enthalpy $h^R$ and the residual isochoric heat capacity $c_i^R$, are sometimes included in the objective function when regressing EoS parameters. We strongly advise against this practice: in our fits we have encountered mean average percentage errors (MAPEs) of exceeding 200% for $h^R$ in the supercritical domain, while the MAPE in $h^sup$ is less than 2%. This is because these properties may become zero, as shown in Figs. 5 and 6.

### 3.1.2. Binary mixtures

Using the optimal parameters for the regression for pure fluids, binary interaction parameters ($k_{ij}$, $l_{ij}$) were fitted by minimizing the total absolute deviation of measured (superscript expt) and calculated (superscript calc) mole fractions, corresponding to the objective function

$$O^{\text{binary}}(k_{ij}, l_{ij}) = \sum_{i} |x_{i}^{\text{expt}} - x_{i}^{\text{calc}}| + \sum_{i} |y_{i}^{\text{expt}} - y_{i}^{\text{calc}}|. \quad (26)$$

Here $x_i$ and $y_i$ are the mole fractions in the liquid and the vapor. The calculations were performed at the same temperature
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** The optimal $\alpha$ and beta functions for four different cases of PR (defined in Section 4), all fitted to the same objective function. The lower temperature limits for each component are given by $T_r$ of Table 1.</figcaption></figure>

and pressure as the measurements. Only regular VLE states were included in the objective function, whereas liquid-liquid equilibrium (LLE) and vapor-liquid-liquid equilibrium (VLLE) measurements were excluded.

## 4. Results and discussion

To evaluate on a neutral basis the influence of the quantum corrections for the covolume parameter, we consider five cases:

**Classic:** All calculations are performed with the present state-of-the-art for classical cubic EoS for quantum fluids, the tc-PR EoS [6].

**Classic-fit:** We use the classical PR EoS but fit the parameters of the $\alpha$ function and the volume shift. This case allows us to evaluate whether the improvement in comparison to the tc-PR EoS comes from re-fitting against a different objective function.

**Case FH1:** The quantum correction derived in Section 2 based on first-order Feynman-Hibbs corrections (FH1) with parameters from Ref. [18] are used. New parameters are regressed for the $\alpha$ function and the volume shift.

**Case FH2:** The quantum correction derived in Section 2 based on the second-order Feynman-Hibbs corrections (FH2) with

parameters from Ref. [18] are used. New parameters are regressed for the $\alpha$ function and the volume shift.

**Case Empirical:** We regress all parameters: (A, B) in Eq. (20), the $\alpha$ function, and the volume shift.

For the case Classic, we use the same values for the critical pressure and temperature as Ref. [6]. For all other cases, we use the values in Table 1.

The cases Classic and Classic-fit are benchmarks for gauging the improvement offered by quantum corrections. For Case FH1 and Case FH2, the quantum corrections for the covolume parameter are predicted. These cases have the same number of fitting parameters as classical cubic EoS. Case Empirical allows assessing whether anything can be gained by also fitting the quantum corrections to the covolume.

Table 3 lists the optimal pure-component parameters (L, M, N) for the Twu $\alpha$ function (Eq. (21)), the parameters (A, B) for the covolume correction (Eq. (20)), and the volume-shift parameter. We also obtained the optimal binary interaction parameters ($k_{ij}$, $l_{ij}$) presented in Table 4.

### 4.1. Single-component systems

Results for the pure components hydrogen, helium, neon and deuterium are summarized for all the cases listed above in Table 5.
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Thermodynamic properties of normal hydrogen (H₂): saturation diagrams (a and b) and supercritical isobars (c–f). The solid curves are calculated from the quantum PR EoS, and the dashed curves are calculated with the reference multiparameter EoS from Ref. [35]. The isobars are plotted between 20 and 300 K.</figcaption></figure>

Evidently the quantum-corrected cubic EoS are vastly superior to the classical cubic EoS, even for the cases with no new fitting parameters (Cases FH1 and FH2). The lowest value of the objective function was in all cases achieved for Case Empirical, since it has more adjustable parameters. Interestingly, for hydrogen, deuterium and neon, this does not appreciably decrease the MAPEs in comparison to Case FH1, and seems to reflect the particular choice of objective function instead of constituting an improved representation of quantum fluids. For these fluids, we recommend the co-volume corrections based on FH1 corrections, since these have a sound theoretical basis. The exception is helium, where Case Empirical is clearly superior and is the recommended model. This is expected, as we showed in previous work that both the FH1 and FH2 corrections were unable to capture the properties of helium below 20 K [18].

The temperature-dependence of the covolume corrections $\beta$ and the regressed $\alpha$ functions are displayed in Fig. 4. The figure shows that the shape of the $\alpha$ function for the classic PR (Case Classic-fit) deviates from the classical monotonically decaying behavior for helium and hydrogen. However, the $\alpha$ functions for all cases that employ quantum corrections for the covolume parameter recover a classical behavior, except for helium at the lowest temperatures. Using a quantum-corrected covolume for helium, it is in fact possible to achieve a good representation with a monotonically decreasing $\alpha$ function, as shown in the row named "Decreasing $\alpha$" in Table 5. Although this improves the accuracy of a few properties, most properties are less accurately predicted. In particular, the MAPE in the saturation pressure increases from 0.67% to 1.18%. In a previous work [18], we showed that for helium, which has the strongest quantum effects, the Feynman-Hibbs corrections are only accurate down to about 20 K. We therefore hypothesize that the functional form of Eq. (20) must be further modified in order to describe helium to an even higher accuracy at the lowest temperatures. Such a modification however, falls beyond the scope of the present work.

The optimal parameters (L, M, N) for the $\alpha$ function of hydrogen and deuterium in Table 3 are notably different from the other components. In particular, the value of L is one to two orders of magnitude larger. To examine this more closely, we have attempted to refit the parameters with fixed values of L, which revealed that the objective function can be reduced further by increasing L. However, the achievable reduction of the objective function is insignificant in comparison to the objective function corresponding to the parameters in Table 3. We hypothesize that this stems from a suboptimal form of the Twu $\alpha$ function, but it seems unlikely that a more flexible form would yield a significantly better fit.

The thermodynamic properties of deuterium are known less accurately than for the other components considered in this work, and the MAPEs in Table 5 should be interpreted with caution. For example, the uncertainties in the vapor pressures are 2% and the uncertainties in the saturated liquid densities are 3% [38]. For deuterium, speeds of sound were not included in the objective function as there is large scatter in the experimental data [38].

We have also regressed parameters for the SRK EoS, and the corresponding parameters and MAPEs are provided in the supplementary information. The main conclusion from this is that, also for SRK, the quantum corrections yield greatly improved agreement with experimental data with no new fitting parameters. This strongly indicates that the hypothesis stated in Section 2.4 is valid; namely that the quantum correction of the covolume parameter is independent of the functional form of the cubic EoS. We find that the performance of SRK is worse than for PR, especially for the liquid-phase densities.

Figs. 5 and 6 illustrate the high accuracy of the quantum-corrected cubic EoS for helium and hydrogen, which are the two components exhibiting the strongest quantum effects.

## 4.2. Thermodynamic consistency

Temperature-dependent covolume corrections are usually avoided in the development of cubic EoS since they can lead
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Thermodynamic properties of helium (⁴He): saturation diagrams (a and b) and supercritical isobars (c–f). The solid curves are calculated from the quantum PR EoS, and the dashed curves are calculated with the reference multiparameter EoS from Ref. [36]. The isobars are plotted between 20 and 150 K.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Contour plots of the molar isochoric heat capacity $c_v$ (Figs. (a)–(d)), and the molar isobaric heat capacity $c_p$ (Figs. (e)–(h)) in the temperature–density plane. The saturation curve is given by the black curve, and negative values are illustrated by gray regions. The lower temperature is the triple point temperature, except for ⁴He where it is 2.5 K. Isobars for 100 bar, 250 bar and 500 bar are also indicated.</figcaption></figure>

to unphysical results. In particular, they lead to diverging and negative isochoric heat capacities at sufficiently high pressures [27]. In Fig. 7, we show that this is also the case for the quantum corrected cubic EoS. However, this only occurs at pressures that are well beyond the domain of industrial interest, and into the solid-formation regime for some of the components. We have also included contour plots of the isobaric heat capacity $c_p$. A negative $c_p$ indicates thermodynamic instability of the pure fluid [12,41],

and such unstable regions should not occur in the regions of known fluid stability; Fig. 7 shows that no spurious instabilities are predicted.

Interestingly, Fig. 7f and h illustrate spurious regions of stability for hydrogen and deuterium. Since the isobaric heat capacity becomes positive within the spinodal region, it follows that the equation of state predicts the existence of a pseudostable phase within the spinodal region for adiabatic and isobaric systems (PS
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Pxy diagrams for the binary mixtures of the fluids considered in this work. The crosses are experimental measurements. The full lines are VLE, the dashed lines are LLE, and the dotted lines indicate VLLE, all calculated with the optimal parameters given in Tables 3 and 4. The colors signify different temperatures.</figcaption></figure>

ensemble) [12,41]. Such a pseudostable phase is in fact also predicted by the classical SRK and PR equations of state, and is therefore not an artefact of quantum corrections. Aursand et al have also found a region of positive $C_p$ within the spinodal region for methane, even for the simpler van der Waals cubic equation of state [12].

Another consistency criterion is that pressure-density isotherms should not cross for pure components close to the critical point. Although this is not a rigorous consistency criterion, such a crossing of near-critical isotherms has never been experimentally observed. For the fluids considered in this work, we have indeed checked that this criterion is satisfied.

### 4.3. The performance for mixtures

Fig. 8 compares calculations using the quantum-corrected PR with available experiments [42–53] for the six binary mixtures: hydrogen-neon, deuterium-neon, hydrogen-deuterium, helium-neon, helium-deuterium and helium-hydrogen. The figure shows that quantum-corrected cubic EoS is in excellent agreement with available experimental results, even for pressures up to 200 bar.

For the hydrogen-helium mixture, we found that it was necessary to include an interaction parameter, $l_{ij}$, that accounts for non-ideal mixing of the covolume parameter. This was expected, as a similar parameter is needed for the cross-interaction parameter of
<caption>**Table 3**</caption>

Optimal parameters (L, M, N) for the Twu α function, (A, B) for the beta function, and c for the Péneloux volume shift. The parameters for the recommended quantum-corrected PR is shaded gray.

<table><caption>**Table 4**</caption><thead><tr><th>Model</th><th>L (-)</th><th>M (-)</th><th>N (-)</th><th>A (K)</th><th>B (K)</th><th>c (cm<sup>3</sup>/mol)</th></tr></thead><tbody><tr><td colspan="7"><strong>Normal hydrogen (H<sub>2</sub>)</strong></td></tr><tr><td>Classic-fit</td><td>2.8994</td><td>-0.61791</td><td>-0.42846</td><td>0</td><td>0</td><td>-4.1101</td></tr><tr><td>Classic</td><td>1.5147</td><td>-3.7959</td><td>-0.1377</td><td>0</td><td>0</td><td>-5.3386</td></tr><tr><td>FH1</td><td>156.21</td><td>-0.0062072</td><td>5.047</td><td>3.0696</td><td>12.682</td><td>-3.8139</td></tr><tr><td>FH2</td><td>347.52</td><td>-0.0027936</td><td>8.2946</td><td>5.8821</td><td>14.791</td><td>-2.9125</td></tr><tr><td>Empirical</td><td>158.54</td><td>-0.0061196</td><td>5.2105</td><td>3.477</td><td>15</td><td>-3.8140</td></tr><tr><td colspan="7"><strong>Helium (<sup>4</sup>He)</strong></td></tr><tr><td>Classic-fit</td><td>-0.046019</td><td>1.2618</td><td>0.69755</td><td>0</td><td>0</td><td>-3.4875</td></tr><tr><td>Classic</td><td>0.0063</td><td>1.2175</td><td>1.0909</td><td>0</td><td>0</td><td>-4.8915</td></tr><tr><td>FH1</td><td>0.18976</td><td>1.3964</td><td>0.58143</td><td>1.8774</td><td>7.7564</td><td>-2.9291</td></tr><tr><td>FH2</td><td>1.1393</td><td>93.272</td><td>0.0044747</td><td>2.7979</td><td>5.2677</td><td>-3.9406</td></tr><tr><td>Empirical</td><td>0.48558</td><td>1.7173</td><td>0.30271</td><td>1.4912</td><td>3.2634</td><td>-3.1791</td></tr><tr><td colspan="7"><strong>Neon (Ne)</strong></td></tr><tr><td>Classic-fit</td><td>0.40805</td><td>0.98441</td><td>0.78674</td><td>0</td><td>0</td><td>-2.6039</td></tr><tr><td>Classic</td><td>0.1887</td><td>0.947</td><td>1.4698</td><td>0</td><td>0</td><td>-2.3573</td></tr><tr><td>FH1</td><td>0.40453</td><td>0.95861</td><td>0.8396</td><td>0.4673</td><td>2.4634</td><td>-2.4665</td></tr><tr><td>FH2</td><td>0.38356</td><td>0.94695</td><td>0.87127</td><td>0.4679</td><td>0.88094</td><td>-2.4556</td></tr><tr><td>Empirical</td><td>0.3981</td><td>0.96535</td><td>0.82696</td><td>0.22069</td><td>-0.65243</td><td>-2.5676</td></tr><tr><td colspan="7"><strong>Normal deuterium (D<sub>2</sub>)</strong></td></tr><tr><td>Classic-fit</td><td>0.3089</td><td>1.0716</td><td>0.6551</td><td>0</td><td>0</td><td>-4.4250</td></tr><tr><td>Classic</td><td>0.1486</td><td>0.9968</td><td>1.0587</td><td>0</td><td>0</td><td>-4.5555</td></tr><tr><td>FH1</td><td>55.007</td><td>-0.016981</td><td>3.1621</td><td>1.6501</td><td>7.309</td><td>-3.8718</td></tr><tr><td>FH2</td><td>63.647</td><td>-0.014525</td><td>3.283</td><td>1.9086</td><td>3.4071</td><td>-3.6319</td></tr><tr><td>Empirical</td><td>52.586</td><td>-0.017779</td><td>3.2179</td><td>2.2117</td><td>12.768</td><td>-3.8717</td></tr></tbody></table>

Optimal interaction parameters ($k_{ij}$, $l_{ij}$) for the quantum PR EoS.

<table><thead><tr><th></th><th>D<sub>2</sub></th><th>H<sub>2</sub></th><th><sup>4</sup>He</th></tr></thead><tbody><tr><td>H<sub>2</sub></td><td>(0,0)</td><td>—</td><td>—</td></tr><tr><td><sup>4</sup>He</td><td>(0.45,0)</td><td>(0.17, -0.16)</td><td>—</td></tr><tr><td>Ne</td><td>(0.18,0)</td><td>(0.18,0)</td><td>(-0.17, 0)</td></tr></tbody></table>

hydrogen-helium as described by quantum-corrected Mie potentials [17]. The semiclassical force field for the helium-hydrogen interaction uses a non-additivity for the interaction diameter of 1.05 [17], i.e. $\sigma_{12} = 1.05(\sigma_1 + \sigma_2)/2$. This translates into a non-additivity for the covolume of $1.05^3 \approx 1.16$, which corresponds to $l_{ij} = -0.16$. Using $l_{ij} = -0.16$ and refitting $k_{ij}$ yielded a significantly improved fit for the helium-hydrogen mixture. We found that little can be gained by further tuning the $l_{ij}$ parameter. For the helium-neon mixture, we verified that including $l_{ij}$ in the fitting process yields a negligible improvement. This is once again in agreement with the

<caption>**Table 5**</caption>

Mean average percentage error (MAPE) relative to the reference equations of state [35–38] for each component. The properties are saturated liquid density $\rho^{\text{sat}}$, saturation pressure $p^{\text{sat}}$, saturated liquid isobaric heat capacity $c_p^{\text{sat}}$, supercritical density $\rho^{\text{sup}}$, supercritical isobaric heat capacity $c_p^{\text{sup}}$, and supercritical speed of sound $w^{\text{sup}}$. Results for the model “Classic” were generated with parameters from Ref. [6]. The MAPEs for the recommended quantum-corrected PR is shaded gray.

<table><thead><tr><th>Model</th><th>p<sup>sat</sup></th><th>ρ<sup>sat</sup></th><th>c<sub>v</sub><sup>sat</sup></th><th>c<sub>p</sub><sup>sat</sup></th><th>Δh<sup>sat</sup></th><th>ρ<sup>sup</sup></th><th>c<sub>v</sub><sup>sup</sup></th><th>c<sub>p</sub><sup>sup</sup></th><th>h<sup>sup</sup></th><th>w<sup>sup</sup></th></tr></thead><tbody><tr><td colspan="11"><strong>Normal hydrogen (H<sub>2</sub>)</strong></td></tr><tr><td>Classic-fit</td><td>0.52</td><td>3.59</td><td>40.72</td><td>14.93</td><td>1.26</td><td>3.18</td><td>4.12</td><td>1.39</td><td>0.52</td><td>5.14</td></tr><tr><td>Classic</td><td>2.02</td><td>1.98</td><td>83.88</td><td>19.12</td><td>2.12</td><td>6.64</td><td>4.75</td><td>3.37</td><td>2.34</td><td>6.76</td></tr><tr><td>FH1</td><td>0.33</td><td>1.10</td><td>4.11</td><td>11.16</td><td>0.93</td><td>0.71</td><td>1.04</td><td>1.05</td><td>0.60</td><td>3.29</td></tr><tr><td>FH2</td><td>0.60</td><td>1.00</td><td>6.58</td><td>14.77</td><td>1.99</td><td>3.43</td><td>1.53</td><td>1.79</td><td>1.01</td><td>6.22</td></tr><tr><td>Empirical</td><td>0.36</td><td>1.11</td><td>3.94</td><td>10.90</td><td>0.90</td><td>0.78</td><td>0.97</td><td>1.07</td><td>0.61</td><td>3.38</td></tr><tr><td colspan="11"><strong>Helium (<sup>4</sup>He)</strong></td></tr><tr><td>Classic-fit</td><td>1.94</td><td>7.61</td><td>2.21</td><td>65.61</td><td>5.18</td><td>6.38</td><td>4.52</td><td>3.39</td><td>3.78</td><td>15.08</td></tr><tr><td>Classic</td><td>2.52</td><td>4.30</td><td>18.06</td><td>51.47</td><td>3.79</td><td>10.50</td><td>8.77</td><td>2.76</td><td>4.74</td><td>23.24</td></tr><tr><td>FH1</td><td>0.81</td><td>4.67</td><td>4.79</td><td>33.59</td><td>2.87</td><td>1.55</td><td>2.58</td><td>1.16</td><td>0.49</td><td>4.62</td></tr><tr><td>FH2</td><td>0.29</td><td>4.88</td><td>5.15</td><td>10.46</td><td>2.84</td><td>1.12</td><td>1.64</td><td>1.50</td><td>0.39</td><td>1.46</td></tr><tr><td>Empirical</td><td>0.67</td><td>1.70</td><td>2.17</td><td>12.26</td><td>1.76</td><td>0.45</td><td>1.64</td><td>0.74</td><td>0.48</td><td>2.57</td></tr><tr><td>Decreasing α</td><td>1.18</td><td>2.18</td><td>1.74</td><td>12.19</td><td>1.45</td><td>0.51</td><td>1.77</td><td>0.70</td><td>0.54</td><td>2.59</td></tr><tr><td colspan="11"><strong>Neon (Ne)</strong></td></tr><tr><td>Classic-fit</td><td>0.57</td><td>1.51</td><td>1.69</td><td>8.66</td><td>0.96</td><td>0.38</td><td>1.66</td><td>0.53</td><td>0.24</td><td>1.53</td></tr><tr><td>Classic</td><td>0.86</td><td>1.78</td><td>13.18</td><td>6.13</td><td>0.82</td><td>1.07</td><td>5.66</td><td>2.22</td><td>0.98</td><td>1.74</td></tr><tr><td>FH1</td><td>0.25</td><td>1.18</td><td>1.99</td><td>8.16</td><td>0.59</td><td>0.57</td><td>2.25</td><td>0.65</td><td>0.23</td><td>2.01</td></tr><tr><td>FH2</td><td>0.26</td><td>1.17</td><td>2.17</td><td>8.46</td><td>0.61</td><td>0.59</td><td>2.30</td><td>0.67</td><td>0.24</td><td>2.04</td></tr><tr><td>Empirical</td><td>0.29</td><td>1.29</td><td>1.74</td><td>7.83</td><td>0.65</td><td>0.40</td><td>1.92</td><td>0.56</td><td>0.21</td><td>1.71</td></tr><tr><td colspan="11"><strong>Normal deuterium (D<sub>2</sub>)</strong></td></tr><tr><td>Classic-fit</td><td>2.04</td><td>1.89</td><td>6.27</td><td>20.14</td><td>2.31</td><td>2.07</td><td>3.09</td><td>0.62</td><td>0.67</td><td>15.79</td></tr><tr><td>Classic</td><td>1.89</td><td>1.43</td><td>7.67</td><td>18.64</td><td>1.88</td><td>2.82</td><td>3.65</td><td>0.88</td><td>0.78</td><td>17.10</td></tr><tr><td>FH1</td><td>0.61</td><td>0.83</td><td>6.55</td><td>14.23</td><td>0.90</td><td>0.60</td><td>0.90</td><td>0.84</td><td>0.44</td><td>10.47</td></tr><tr><td>FH2</td><td>0.78</td><td>1.17</td><td>11.45</td><td>20.52</td><td>1.22</td><td>1.07</td><td>1.01</td><td>0.99</td><td>0.35</td><td>9.07</td></tr><tr><td>Empirical</td><td>0.60</td><td>0.81</td><td>5.66</td><td>12.80</td><td>0.85</td><td>0.70</td><td>0.82</td><td>0.86</td><td>0.44</td><td>10.05</td></tr></tbody></table>
semiclassical force field, where non-additivity effects were found to be negligible [17].

The helium-neon mixture has been highlighted as a particu-
larly difficult example where most available EoS, including cubic
EoS with advanced mixing rules have fallen short [12]. Tkaczuk
et al. [23] recently developed a Helmholtz energy mixture model
for the helium-neon mixture that has been included in Fig. 8 for
comparison. Their model relies on multiparameter pure fluid EoS
combined with multiparameter mixture models. Tkaczuk et al. use
more than 40 binary interaction parameters to regress the proper-
ties of the binary helium-neon mixture model. The quantum cor-
rected PR EoS requires only a constant binary interaction parame-
ter to achieve a good agreement with experimental results for the
VLE. Moreover, care should be taken in this comparison, as the ex-
perimental data for helium-neon have significant uncertainty, e.g.
they do not extrapolate to the correct saturation pressure for pure
neon [17].

The excellent agreement displayed in Fig. 8 shows that the quantum-corrected cubic EoS are among the most accurate EoS available for calculating the properties of mixtures of fluids that exhibit strong quantum effects. A moderate overestimation of the critical pressure for the helium–hydrogen mixture at 29.0 K (blue curve) is observed, which is surprising considering that the critical pressure agrees well at 31.5 K (brown curve). We found that this disagreement can be explained by the reported temperature from experiments having a 0.2 K uncertainty in the temperature measurement. The calculated phase envelope at 29.2 K yields excellent agreement with the experimental results for the critical pressure, which suggests that uncertainty in the reported temperature can explain the observed deviation.

5. Conclusion

Cubic EoS are often the preferred choice when computational
speed in combination with a reasonable accuracy is required. Un-
til now, they have yielded inaccurate predictions for fluids that
exhibit strong quantum effects, such as helium and hydrogen at
low temperatures. Even with advanced mixing rules, cubic EoS
are unable to reproduce the vapor-liquid equilibrium (VLE) of the
helium-neon mixture to a high accuracy. Moreover, the α func-
tions of helium and hydrogen published in the literature are not
monotonically decreasing like those of other fluids.

In this work, we have derived quantum corrections for the co-
volume parameter of cubic EoS. The derivation of their functional
form and temperature dependence was based on Feynman–Hibbs-
corrected Mie potentials. The corrections were shown to be inde-
pendent of the choice of cubic EoS. Regression of the Twu α func-
tion after incorporating the quantum correction for the covolume
parameter recovered a monotonically decaying behavior for hydro-
gen, and also for helium except at the lowest temperatures. New
α functions and additional volume-shifts were developed for both
the Peng–Robinson and the Soave–Redlich–Kwong cubic EoS.

The quantum corrections were shown to significantly improve
the accuracy of the predictions compared to the present state-of-
the-art (tc-PR), even with no new fitting parameters for the cases
when the covolume parameters were theoretical predictions. The
final values for the quantum corrections were those that most ac-
curately reproduced the experimental data. The quantum-corrected
PR improved the average error in the isochoric heat capacity at sat-
uration from above 80% to 4.1% for hydrogen, and the average er-
ror in the predictions of the supercritical speed of sound of helium
from 23% to 2.6%. Most other thermodynamic properties, both sub-
and supercritical, were also predicted much more accurately after
introducing the quantum correction. The additional introduction of
a constant volume correction reduced the error in density predic-

tions to below 2% for helium and around 1% for hydrogen, neon
and deuterium.

The VLE of all investigated mixtures were reproduced to a high accuracy with simple mixing rules and temperature-independent interaction parameters. In particular, the quantum-corrected cubic EoS is the first cubic EoS presented in open literature that is able to reproduce the experimental VLE for helium-neon. The hydrogen-helium mixture required an additional interaction parameter for the covolume, similar to the quantum-corrected force fields for the very same mixture.

The quantum-corrected cubic EoS offers superior accuracy
while retaining the unparalleled computational speed and numeri-
cal simplicity of cubic EoS. It therefore opens the door for a wide
range of technological developments involving fluids that exhibit
strong quantum effects. One attractive possibility is the develop-
ment of the next generation of hydrogen liquefaction processes.
Helium-hydrogen-neon mixtures have been touted as a promising
"quantum refrigerant" with the potential to significantly enhance
the efficiency of the liquefaction process.

The description of saturated heat capacities for quantum fluids
still has room for improvement, but is likely sufficiently accurate
for most industrial applications, such as hydrogen liquefaction pro-
cess design.

Declaration of Competing Interest

The authors declare that they have no known competing finan-
cial interests or personal relationships that could have appeared to
influence the work reported in this paper.

Acknowledgments

Øivind Wilhelmsen was supported by the ENERSENSE initiative.

**Appendix A. A comment on the implementation**

Temperature-dependent covolume parameters are not new in the literature on cubic EoS, and occur for example in the popular Wong-Sandler mixing rule [54]. They are conveniently implemented using the approach by Michelsen and Mollerup [55] based on the chain rule: derivatives of the Helmholtz energy are computed first with respect to the parameters such as *b* and then multiplied with derivatives of state variables such as temperature. In this way, the temperature dependence of the covolume can be handled similarly to its composition dependence.

**Supplementary material**

Supplementary material associated with this article can be
found, in the online version, at 10.1016/j.fluid.2020.112790

**References**

[1] J.D. Van Der Waals, J.S. Rowlinson, On the Continuity of the Gaseous and Liquid States, Courier Corporation, 2004.

[2] G. Skaugen, M. Hammer, P.E. Wahl, Ø. Wilhelmsen, Constrained non-linear optimisation of a process for liquefaction of natural gas including a geometrical and thermo-hydraulic model of a compact heat exchanger, Comput. Chem. Eng. 73 (2015) 102.

[3] P. Mathias, J. Boston, S. Watanasiri, Effective utilization of equations of state for thermodynamic properties in process simulation, AIChE J. 30 (2) (1984) 182-186, doi:10.1002/aic.690300203.

[4] G. Soave, Equilibrium constants from a modified Redlich-Kwong equation of state, Chem. Eng. Sci. 27 (6) (1972) 1197-1203, doi:10.1016/0009-2509(72) 80096-4.

[5] D.Y. Peng, D.B. Robinson, A new two-constant equation of state, Ind. Eng. Chem. Fund. 15 (1976) 59, doi:10.1021/i160057a011.

[6] Y. Le Guennec, R. Privat, J.-N. Jaubert, Development of the translated-consistent tc-PR and tc-RK cubic equations of state for a safe and accurate prediction of volumetric, energetic and saturation properties of pure compounds in the sub- and super-critical domains, Fluid Phase Equilib. 429 (2016) 301-312, doi:10.1016/j.fluid.2016.09.003.
[7] G.M. Kontogeorgis, M.L. Michelsen, G.K. Folas, S. Derawi, N. von Solms, E.H. Stenby, Ten years with the CPA (Cubic-Plus-Association) equation of state. Part 1. Pure compounds and self-associating systems, Ind. Eng. Chem. Res. 45 (14) (2006) 4855–4868, doi:10.1021/ie051305v.

[8] B. Maribo-Mogensen, Development of an Electrolyte CPA Equation of state for Applications in the Petroleum and Chemical Industries, Technical University of Denmark, 2014 Ph.D. thesis.

[9] J.S. Lopez-Echeverry, S. Reif-Acherman, E. Araujo-Lopez, Peng-Robinson equation of state: 40 years through cubics, Fluid Phase Equilib. 447 (2017) 39–71, doi:10.1016/j.fluid.2017.05.007.

[10] G.M. Kontogeorgis, G. Folas, Thermodynamic Models for Industrial Applications, Wiley, 2010, doi:10.1002/9780470747537.

[11] Fuel Cells and Hydrogen, The European Green Deal: Hydrogen is a priority area for a clean and circular economy, 2019, (www.fch.europa.eu/news/european-green-deal-hydrogen-priority-area-clean-and-circular-economy). Accessed: 21.12.19.

[12] Ø. Wilhelmsen, D. Berstad, A. Aasen, P. Nekså, G. Skaugen, Reducing the exergy destruction in the cryogenic heat exchangers of hydrogen liquefaction processes, Int. J. Hydrogen Energy 43 (10) (2018) 5033–5047, doi:10.1016/j.ijhydene.2018.01.094.

[13] H. Quack, Conceptual design of a high efficiency large capacity hydrogen liquefier, in: AIP Conf.Proc., 613, 2002, pp. 255–263, doi:10.1063/1.1472029.

[14] A. Aasen, M. Hammer, G. Skaugen, J.P. Jakobsen, Ø. Wilhelmsen, Thermodynamic models to accurately describe the PVTxy-behavior of water/carbon dioxide mixtures, Fluid Phase Equilib. 442 (2017) 125–139, doi:10.1016/j.fluid.2017.02.006.

[15] J.-P. Hansen, J.-J. Weis, Quantum corrections to the coexistence curve of neon near the triple point, Phys. Rev. 188 (1969) 314–318, doi:10.1103/PhysRev.188.314.

[16] J.-P. Hansen, I.R. McDonald, Theory of Simple Liquids: With Applications to Soft Matter, Academic Press, 2013.

[17] A. Aasen, M. Hammer, E.A. Müller, Ø. Wilhelmsen, Equation of state and force fields for Feynman–Hibbs-corrected Mie fluids. II. Application to mixtures of helium, neon, hydrogen, and deuterium, J. Chem. Phys. 152 (2020) 074507, doi:10.1063/1.5136079.

[18] A. Aasen, M. Hammer, A. Ervik, E.A. Müller, Ø. Wilhelmsen, Equation of state and force fields for Feynman–Hibbs-corrected Mie fluids. I. Application to pure helium, neon, hydrogen, and deuterium, J. Chem. Phys. 151 (6) (2019) 064508, doi:10.1063/1.5111364.

[19] R.P. Feynman, A.R. Hibbs, D.F. Styer, Quantum Mechanics and Path Integrals, Emended Edition, McGraw-Hill, New York, 2005.

[20] M. Hammer, A. Aasen, Å. Ervik, Ø. Wilhelmsen, Choice of reference, the influence of non-additivity and challenges in thermodynamic perturbation theory for mixtures, J. Chem. Phys. 152 (2020) 134106, doi:10.1063/1.5142771.

[21] R.J. Thomas, R. Dutta, P. Ghosh, K. Chowdhury, Applicability of equations of state for modeling helium systems, Cryogenics 52 (7–9) (2012) 375–381, doi:10.1016/j.cryogenics.2012.03.002.

[22] S. Kaviani, F. Feyzi, B. Khosravi, Modification of the Peng–Robinson equation of state for helium in low temperature regions, Phys. Chem. Liq. 54 (4) (2016) 542–551, doi:10.1080/00319104.2015.1126836.

[23] J. Tkaczuk, I.H. Bell, E.W. Lemmon, N. Luchier, F. Millet, Equations of state for the thermodynamic properties of binary mixtures for helium-4, neon, and argon, J. Phys. Chem. Ref. Data 49 (2) (2020) 023101, doi:10.1063/1.5142275.

[24] Y. Le Guennec, S. Lasala, R. Privat, J.-N. Jaubert, A consistency test for α-functions of cubic equations of state, Fluid Phase Equilib. 427 (2016) 513–538, doi:10.1016/j.fluid.2016.07.026.

[25] Y. Le Guennec, R. Privat, S. Lasala, J.-N. Jaubert, On the imperative need to use a consistent α-function for the prediction of pure-compound supercritical properties with a cubic equation of state, Fluid Phase Equilib. 445 (2017) 45–53.

[26] I.H. Bell, M. Satyro, E.W. Lemmon, Consistent twu parameters for more than 2500 pure fluids from critically evaluated experimental data, J. Chem. Eng. Data 63 (7) (2018) 2402–2409, doi:10.1021/acs.jced.7b00967.

[27] V. Kalikhman, D. Kost, I. Polishuk, About the physical validity of attaching the repulsive terms of analytical EOS models by temperature dependencies, Fluid Phase Equilib. 293 (2) (2010) 164–167, doi:10.1016/j.fluid.2010.03.003.

[28] M.-J. Huron, J. Vidal, New mixing rules in simple equations of state for representing vapour-liquid equilibria of strongly non-ideal mixtures, Fluid Phase Equilib. 3 (4) (1979) 255–271, doi:10.1016/0378-3812(79)80001-1.

[29] G.M. Kontogeorgis, R. Privat, J.-N. Jaubert, Taking another look at the van der Waals equation of state–almost 150 years later, J. Chem. Eng. Data 64 (11) (2019) 4619–4637, doi:10.1021/acs.jced.9b00264.

[30] T. Lafitte, A. Apostolakou, C. Avendaño, A. Galindo, C.S. Adjiman, E.A. Müller, G. Jackson, Accurate statistical associating fluid theory for chain molecules formed from mie segments, J. Chem. Phys. 139 (15) (2013) 154504, doi:10.1063/1.4819786.

[31] C.H. Twu, D. Bluck, J.R. Cunningham, J.E. Coon, A cubic equation of state with a new alpha function and a new mixing rule, Fluid Phase Equilib. 69 (1991) 33–50, doi:10.1016/0378-3812(91)90024-2.

[32] A. Péneloux, E. Rauzy, R. Fréze, A consistent correction for Redlich-Kwong-Soave volumes, Fluid Phase Equilib. 8 (1) (1982) 7–23, doi:10.1016/0378-3812(82)80002-2.

[33] J.-N. Jaubert, R. Privat, Y. Le Guennec, L. Coniglio, Note on the properties altered by application of a Péneloux-type volume translation to an equation of state, Fluid Phase Equilib. 419 (2016) 88–95, doi:10.1016/j.fluid.2016.03.012.

[34] R. Privat, J.-N. Jaubert, Y. Le Guennec, Incorporation of a volume translation in an equation of state for fluid mixtures: which combining rule? Which effect on properties of mixing? Fluid Phase Equilib. 427 (2016) 414–420, doi:10.1016/j.fluid.2016.07.035.

[35] J.W. Leachman, R.T. Jacobsen, S.G. Penoncello, E.W. Lemmon, Fundamental equations of state for parahydrogen, normal hydrogen, and orthohydrogen, J. Phys. Chem. Ref. Data 38 (2009) 721–748, doi:10.1063/1.3160306.

[36] D.O. Ortiz-Vega, A new wide range equation of state for helium-4, Texas A&M University, 2013 Ph.D. thesis.

[37] M. Thol, R. Beckmüller, R. Weiss, A. Harvey, E. Lemmon, R. Jacobsen, R. Span, Thermodynamic properties for neon for temperatures from the triple point to 700K at pressures to 700MPa, J. Phys. Chem. Ref. Data (2020). (to appear)

[38] I.A. Richardson, J.W. Leachman, E.W. Lemmon, Fundamental equation of state for deuterium, J. Phys. Chem. Ref. Data 43 (1) (2014) 013103, doi:10.1063/1.4864752.

[39] I.H. Bell, J. Wronski, S. Quoilin, V. Lemort, Pure and pseudo-pure fluid thermophysical property evaluation and the open-source thermophysical property library coolprop, Ind. Eng. Chem. Res. 53 (6) (2014) 2498–2508, doi:10.1021/ie4033999.

[40] J.A. Nelder, R. Mead, A simplex method for function minimization, Comput. J. 7 (4) (1965) 308–313, doi:10.1093/comjnl/7.4.308.

[41] P. Debenedetti, Metastable Liquids: Concepts and Principles, Princeton University Press, Princeton, 1996.

[42] W.B. Streett, C.H. Jones, Liquid phase separation and liquid-vapor equilibrium in the system neon–hydrogen, J. Chem. Phys. 42 (11) (1965) 3989–3994, doi:10.1063/1.1695872.

[43] C.K. Heck, P. Barrick, Liquid-vapor phase equilibria of the neon-normal hydrogen system, in: Adv. Cryog. Eng., Springer, 1966, pp. 349–355, doi:10.1007/978-1-4757-0522-5_38.

[44] M. Knorn, Vapour–liquid equilibria of the neon–helium system, Cryogenics 7 (1967) 177, doi:10.1016/S0011-2275(67)80063-8.

[45] C. Heck, P. Barrick, Liquid-vapor equilibria of the neon–helium system, in: Adv. Cryog. Eng., Springer, 1967, pp. 714–718, doi:10.1007/978-1-4757-0489-1_75.

[46] M. Simon, On the phase separation in the liquid system Ne-nD₂, Phys. Lett. 2 (1962) 234, doi:10.1016/0031-9163(62)90241-X.

[47] W. Streett, Liquid phase separation and liquid-vapor equilibrium in the system neon–deuterium, in: Proc. 2nd Int. Cryog. Eng. Conf., Iliffe Science and Technology Publ. Ltd., Guildford, England, 1968, pp. 260–263.

[48] M. Hiza, Liquid-vapor equilibria in binary systems containing He-4 or He-3 with nH₂ or nD₂, Fluid Phase Equilib. 6 (3–4) (1981) 203–227, doi:10.1016/0378-3812(81)85005-4.

[49] R. Newman, L. Jackson, The P,T,x relationships of H₂+HD and H₂+D₂ mixtures between 18 and 28 K, Trans. Faraday Soc. 54 (1958) 481, doi:10.1039/TF9585401481.

[50] W.B. Streett, R.E. Sonntag, G.J. Van Wylen, Liquid-vapor equilibrium in the system normal hydrogen–helium, J. Chem. Phys. 40 (5) (1964) 1390–1395, doi:10.1063/1.1725322.

[51] R.E. Sonntag, G.J. Van Wylen, R.W. Crain, Liquid-vapor equilibrium in the system equilibrium hydrogen–helium, J. Chem. Phys. 41 (8) (1964) 2399–2402, doi:10.1063/1.1726277.

[52] C.M. Sneed, R.E. Sonntag, G.J. Van Wylen, Helium–hydrogen liquid-vapor equilibrium to 100 atm, J. Chem. Phys. 49 (5) (1968) 2410–2414, doi:10.1063/1.1670414.

[53] W. Streett, Phase equilibria in molecular hydrogen-helium mixtures at high pressures, Astrophys. J. 186 (1973) 1107–1125, doi:10.1086/152575.

[54] D.S.H. Wong, S.I. Sandler, A theoretically correct mixing rule for cubic equations of state, AIChE J. 38 (5) (1992) 671–680, doi:10.1002/aic.690380505.

[55] M.L. Michelsen, J.M. Mollerup, Thermodynamic Models: Fundamentals & Computational Aspects, second ed., Tie-Line Publications, Holte, Denmark, 2007.