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

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 α 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 α 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.


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 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 λ B [15,16] : where h is Planck's constant, m is the particle mass, and k B is Boltzmann's constant. When λ B becomes comparable to the typical intermolecular separation, quantum effects can influence the fluid properties significantly.
Recently, the statistical associating fluid theory of quantumcorrected 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 heliumneon, 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 λ 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 effort s 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 α 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 α function, as: a (T ) = a c α(T ) . Which functional form to use for the α function has been widely debated.
Recently, Le Guennec et al. [24] introduced consistency criteria for 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 / = 1 , and the particle mass is that of helium. The FH correction increases the effective diameter σ eff , and decreases the well depth eff , compared to the classical LJ parameters, σ and . The inset shows that the effective range of the potential increases slightly. the α 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 α 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 α 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 acentric factors being negative [6] . We will show that a near classical behavior with a monotonic decay can be obtained for the α function by incorporating temperature-dependent quantum corrections for the covolume parameter that explicitly account 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] .

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, σ , to a larger, temperature-dependent value, σ eff (T ) .
2. They reduce the well-depth, of the classical interaction potential to a smaller, temperature-dependent value 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 α 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] .

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: where σ is the particle diameter. The covolume, b , used in the equation of state can be found by multiplying b with Avogadro's number, 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 where T c is the critical temperature. This gives the following expressions: for the covolumes of the van der Waals, the SRK and the PR equations of state respectively.

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: where is the well depth, σ is the finite distance at which the potential is zero, λ a and λ r are the attractive and repulsive exponents, and 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 σ allows particles to be at intermolecular distances smaller than σ . However, as the repulsive part of the interaction potential is usually very steep, the excluded volume can be approximated by the positive part of the potential. For Mie potentials, the particle diameter is then simply σ 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 Mie = πσ 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] .

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, 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 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, σ eff , given by the separation where the quantum-corrected potential equals zero. The effective diameter was found by solving for r using a Newton-Raphson solver with analytical derivatives and the value of σ as initial guess.

The quantum-corrected covolume parameter
We will relate the quantum covolume b q ( T ) to the classical co- 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 σ eff (T −1 = 0) = σ ; in other words, the classical value is recovered at high temperatures. Differentiating σ eff as defined by Eq. (14) with respect to T −1 gives that the following slope should be satisfied in the limit T −1 → 0 : Eq. (16) holds for both the first and the second order quantumcorrected Mie potentials. In the limit T −1 → ∞ however, the first order quantum-corrected Mie potential reaches the maximum value .
The following expression can satisfy all these constraints: where The reason for the prefactor c F H = 1 . 4 for first order-corrections and c F H = 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 → ∞ 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 σ eff /σ obtained by solving Eq. (14) , in particular at supercritial temperatures.
The derived β function can be written as where A are 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.

The α 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 α function of cubic EoS. Since the quantum-corrected covolume is larger at subcritical temperatures, the attraction and thus the α function must be reduced, e.g.
to recover the same saturation pressures. This is because the two  (1) ) are regressed together to approximate the experimental data. For supercritical temperatures, the covolume is smaller than the classical covolume, and thus the quantum α function should be smaller than the classical one. These effects will change the "bell-shaped" α function found previously for hydrogen and helium for the classical PR. However, it is not guaranteed that the change will make the α 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 α function is correlated using the three-parameter Twu function [31] : Eq. (21) is designed such that α(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 α function becomes a monotonic function when incorporating a quantumcorrected covolume.

Mixing and combining rules
For a we apply the usual quadratic mixing rules: 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 A common choice is to set l i j = 0 , which converts the mixture covolume into a linear combination of the single-component covolumes.
In addition to the α and β 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:

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.

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.  20) ), were fitted to pseudoexperimental 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:

Saturation
where i indexes the states where the deviations are calculated, ξ represents one of the ten thermodynamic properties considered, and w ξ is the weight given to the property ξ . 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, 20 0 and 30 0 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 R v , 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 .

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 Here x i and y i are the mole fractions in the liquid and the vapor. The calculations were performed at the same temperature  Table 1 . 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.

Results and discussion
To evaluate on a neutral basis the influence of the quantum corrections for the covolume parameter, we consider five cases: 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 α 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 .

Single-component systems
Results for the pure components hydrogen, helium, neon and deuterium are summarized for all the cases listed above in Table 5 . 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 covolume 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 β and the regressed α functions are displayed in Fig. 4 . The figure shows that the shape of the α function for the classic PR (Case Classic-fit) deviates from the classical monotonically decaying behavior for helium and hydrogen. However, the α 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 α function, as shown in the row named "Decreasing α" 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 α 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 α 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 quantumcorrected cubic EoS for helium and hydrogen, which are the two components exhibiting the strongest quantum effects.

Thermodynamic consistency
Temperature-dependent covolume corrections are usually avoided in the development of cubic EoS since they can lead   (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 4 He where it is 2.5 K. Isobars for 100 bar, 250 bar and 500 bar are also indicated.
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. 7 f 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 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.  [42][43][44][45][46][47][48][49][50][51][52][53] for the six binary mixtures: hydrogen-neon, deuterium-neon, hydrogen-deuterium, heliumneon, 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.

The performance for mixtures
For the hydrogen-helium mixture, we found that it was necessary to include an interaction parameter, l ij , that accounts for nonideal mixing of the covolume parameter. This was expected, as a similar parameter is needed for the cross-interaction parameter of Table 3 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-correcte d PR is shaded gray.  Table 4 Optimal interaction parameters ( k ij , l ij ) for the quantum PR EoS.
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. σ 12 = 1 . 05(σ 1 + σ 2 ) / 2 . This translates into a non-additivity for the covolume of 1.05 3 ≈ 1.16, which corresponds to l i j = −0 . 16 . Using l i j = −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 Table 5 Mean average percentage error (MAPE) relative to the reference equations of state [35][36][37][38] for each component. The properties are saturated liquid density ρ sat , saturation pressure p sat , saturated liquid isobaric heat capacity c sat p , supercritical density ρ sup , supercritical isobaric heat capacity c sup p , and supercritical speed of sound w sup . Results for the model "Classic" were generated with parameters from Ref. [6] . T he MAPEs for the recommended quantum-corrected PR is shaded gray. semiclassical force field, where non-additivity effects were found to be negligible [17] . The helium-neon mixture has been highlighted as a particularly 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 properties of the binary helium-neon mixture model. The quantum corrected PR EoS requires only a constant binary interaction parameter to achieve a good agreement with experimental results for the VLE. Moreover, care should be taken in this comparison, as the experimental 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.

Conclusion
Cubic EoS are often the preferred choice when computational speed in combination with a reasonable accuracy is required. Until 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 α functions 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 covolume parameter of cubic EoS. The derivation of their functional form and temperature dependence was based on Feynman-Hibbscorrected Mie potentials. The corrections were shown to be independent of the choice of cubic EoS. Regression of the Twu α function after incorporating the quantum correction for the covolume parameter recovered a monotonically decaying behavior for hydrogen, 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-ofthe-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 accurately reproduced the experimental data. The quantum-corrected PR improved the average error in the isochoric heat capacity at saturation from above 80% to 4.1% for hydrogen, and the average error in the predictions of the supercritical speed of sound of helium from 23% to 2.6%. Most other thermodynamic properties, both suband 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 hydrogenhelium 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 numerical 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 development 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 process design.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.