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

# Application of thermodynamics at different scales to describe the behaviour of fast reacting binary mixtures in vapour-liquid equilibrium

Silvia Lasala<sup>a, *</sup>, Konstantin Samukov<sup>a</sup>, H. Mert Polat<sup>b</sup>, Véronique Lachet<sup>c</sup>, Olivier Herbinet<sup>a</sup>, Romain Privat<sup>a</sup>, Jean-Noël Jaubert<sup>a</sup>, Othonas A. Moultos<sup>b</sup>, Kevin De Ras<sup>d</sup>, Thijs J. H. Vlugt<sup>b, *</sup>

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

## ARTICLE INFO

**Keywords:**

Reactive mixtures

Vapour-liquid equilibrium

Monte Carlo simulations

Quantum Mechanics simulations

Equations of state

## ABSTRACT

The use of reactive working fluids in thermodynamic cycles is currently being considered as an alternative to inert working fluids, because of the preliminarily attested higher energy-efficiency potential. The current needs to simulate their use in thermodynamic cycles, which may operate in liquid, vapour or vapour-liquid state, are an accurate real-fluid equation of state and ideal gas thermochemical properties of each molecule constituting the mixture, to calculate the equilibrium constant. To this end, the appeal to a multi-scale theoretical methodology is paramount and its definition represents the objective of the present work. This methodology is applied and validated on the system $N_2O_4 \rightleftharpoons 2NO_2$. Firstly, the equations solved for simultaneous two-phase and reaction equilibrium are presented. Secondly, ideal gas thermochemical properties of $N_2O_4$ and $NO_2$ are computed at atomic scale by quantum mechanics simulations. Then, to apply the selected cubic equation of state, pure-component properties of the species forming the reactive mixture (critical point coordinates and acentric factor) are required as input. However, these properties are not measurable, since $NO_2$ and $N_2O_4$ do not exist in nature as pure components. To get around this difficulty, the methodology relies on molecular Monte Carlo simulations of the pure $N_2O_4$ and $NO_2$, as well as on the reactive $N_2O_4 \rightleftharpoons 2NO_2$, enabling the determination of those missing pure-component properties and thus the calculation, on a macroscopic scale, of the reactive mixture properties. Finally, the comparison of calculated mixture properties with available experimental data leads to validate the accuracy of the proposed methodology.

## 1. Introduction

There exist few reversible chemical reactions which evolve *very rapidly* towards chemical equilibrium, whether subjected to a modification of their intensive properties, such as temperature and pressure. According to the specific use of such fluids, the high reaction rate can be approximated as being infinite and thus the chemical equilibrium state considered as being instantaneously attained.

One of those few systems is formed by the molecules $N_2O_4$ and $NO_2$, whose proportion is dictated by the chemical equilibrium achieved by the reversible reaction $N_2O_4 \rightleftharpoons 2NO_2$, at the temperature and pressure of the system. At present, the $NO_2/N_2O_4$ gas mixture is used as a powerful hypergolic oxidizer of fuels like hydrazine, $N_2H_4$, in space and

military rockets [1–3]. Between the years 1950 and 1990, the possible exploitation of the chemical properties of such a reactive fluid has also been explored in power plants, proposing it as a working fluid of thermodynamic cycles [2–13]. Indeed, the spontaneous, fast and reversible evolution of the chemical reaction when the fluid undergoes a modification of its temperature and pressure, while crossing the unit operations (turbine, heat exchangers, pumps, compressors, etc.), implies the involvement of chemical energy in the thermodynamic transformations, with significant effects on the performance of the thermodynamic cycle and on the sizing of units, such as the turbine [10]. Many studies have been conducted on the effect of using such a working fluid (commonly called “dissociating gas” in the literature) in power cycles, or on the characterisation of its thermochemical and thermophysical properties.

<footnote><sup>a</sup> Université de Lorraine, CNRS, LRGP, F-54000 Nancy, France</footnote>

<footnote><sup>b</sup> Engineering Thermodynamics, Process & Energy Department, Faculty of Mechanical Engineering, Delft University of Technology, Leeghwaterstraat 39, Delft 2628CB, The Netherlands</footnote>

<footnote><sup>c</sup> IFP Energies nouvelles, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison, France</footnote>

<footnote><sup>d</sup> Ghent University, Laboratory for Chemical Technology (LCT), Technologiepark 125, B-9052 Ghent, Belgium</footnote>

<footnote>* Corresponding authors.</footnote>

<footnote>E-mail addresses: silvia.lasala@univ-lorraine.fr (S. Lasala), t.j.h.vlugt@tudelft.nl (T. J. H. Vlugt).</footnote>
The above-mentioned works have theoretically demonstrated that the thermodynamic efficiency of supercritical (Brayton) and subcritical (Rankine) cycles operating with reactive working fluids, mainly N₂O₄, is larger than the efficiency of cycles operating with specific inert fluids, such as water, ammonia, helium or CO₂. However, the comparison between, for example, the dissociating molecule N₂O₄ and the inert molecule H₂O is not straightforward, as multiple molecular aspects should be considered, not only related to the chemical reaction but also to molecular heterogeneity (N₂O₄/NO₂ vs. H₂O). To enable a fairer comparison between reactive and inert fluids, other more generalised studies have considered fictitious chemical reactions (Aₙ ⇌ (n/m) Aₘ) and investigated the effect of the reaction stoichiometry (n, m) and thermochemistry (enthalpy and entropy of reaction) on the performance of the Brayton power cycles and heat pumps, with respect to pure fluids (reactants or products) [14,15]. Two main conclusions were that the use of reactive working fluids may lead to the improvement of +30 % in the efficiency of power plants [14] and a more than doubled coefficient of performance of heat pumps based on a gaseous inversed Brayton cycle [15]. The reliability and advantages related to the use of N₂O₄ have been experimentally validated at the Nuclear Power Institute of the Academy of Science, in Byelorussia, in mock-ups of 100 – 1000 kWth [16]. After 1986, with the Chernobyl disaster, Byelorussia stopped activities related to nuclear research and thus those studies on dissociating fluids were discontinued too. Not a single experimental or full-scale power plant operates with N₂O₄, nor with other similarly proposed dissociating gases, such as Al₂Cl₆ ⇌ 2 AlCl₃ and Al₂Br₆ ⇌ 2 AlBr₃. A more complete review of past technological studies about reactive fluids and their applications is provided in Lasala et al. [14].

Theoretical studies on the use or characterization of reactive fluids potentially usable to enhance the performance of unit operations or thermodynamic cycles are nowadays again under development. Heat transfer properties of reactive N₂O₄ are investigated by computational fluid dynamics in various heat exchangers [17,18]. One of the main conclusions is that the presence of the N₂O₄/NO₂ chemical reaction in the heat exchange process enhances the heat transfer coefficient by 600 % with respect to non-reactive fluids. The problem underlying the limited spreading of using reactive working fluids in thermodynamic cycles is the scarce availability of similarly reversible fast chemical reactions, other than N₂O₄. Taking that into consideration, a research project named REACHER [19,20], funded by the European Research Council, aims to discover novel working fluids, characterise their thermodynamic and kinetic properties, to optimise the thermodynamic cycle architecture and, finally, to validate cycle's performance calculations on an experimental micro-power plant.

The high number of studies performed in the past has made many thermodynamic data available on the reactive N₂O₄/NO₂, Al₂Cl₆/AlCl₃ and Al₂Br₆/AlBr₃. Most of the measured data are collected in the Thermo Data Engine database (TDE 10.1) [21] and in the Dortmund Data Base (DDB2020) [22]. The thermodynamic modelling of N₂O₄ is more investigated than Al₂Cl₆ and Al₂Br₆. Calculations of thermodynamic properties of N₂O₄ have been performed in different works, with the use of semi-empirical equations of state [23–31] or Monte Carlo simulations [32]. Stai et al. [23] used a modification of the Wohl equation to describe volumetric properties of gaseous N₂O₄, but the specifications for the critical coordinates are not explicitly provided. Seshadri [24] used virial equation of state parametrized over experimental data points and calculated the fugacity coefficients with the model developed by Hougen and Watson [25]. McCarty et al. [26] performed calculations using the Peng-Robinson equation of state, and optimised the coordinates of the critical point of the molecules forming the reactive mixture over available volumetric experimental data and assigned, as a preliminary assumption, the same acentric factor to NO₂ and to N₂O₄. Also, McCarty [26] optimised the 32 parameters of a modified Benedict-Webb-Rubin equation of state to reproduce experimental P-ρ-T data of the equilibrium mixture, treating it as a pseudo-pure component. With the aim to perform calculations on a power cycle, Binotti et al. [27] and

Manzolini et al. [28] also treated the system as a pseudo-pure component and modelled it with the Peng-Robinson equation of state. Vapor-liquid equilibrium of N₂O₄/NO₂ was calculated by De Souza and Deiters [29], with both a hard sphere attractive model and a semiempirical equation of state optimized over experimental and Monte Carlo data. Then, Belkadi et al. [30] modelled VLE with soft-SAFT equation of state, without considering reactional equilibrium but treating the molecule NO₂ as a self-associating component.

Moreover, Monte Carlo simulations have been performed on similar reactive mixtures being in phase and reaction equilibrium. Monte Carlo codes including Gibbs and reaction Monte Carlo were introduced by Johnson et al. [33], who treated the equilibrium of the reactive (NO₂)/NO, and Smith and Triska [34]. The reactive system N₂O₄/NO₂ system has been previously characterised by MC reaction and Gibbs ensembles by Bourasseau et al. [32], who also provided force field parameters optimized over experimental data.

Recently, Lemmon et al. [31] produced a thorough review work on the thermodynamic modelling and experimental characterisation of this system, and also developed a highly accurate multiparameter equation of state being specific to the reactive N₂O₄. This equation of state also accounts for the occurrence of the higher temperature reaction NO₂ ⇌ 2NO + O₂, other than N₂O₄ ⇌ 2NO₂. In this model, authors used the pure-component critical temperature and density they determined for pure N₂O₄ and pure NO₂ by Monte Carlo Gibbs ensemble simulations, using the open source code Cassandra 1.2.5 [35] and force-field parameters from Bourasseau et al. [32]. Moreover, the parameters of the considered Helmholtz-based equation of state have been optimised on the basis of a high number of available experimental data that they reviewed in the same work [31].

As a matter of fact, at the basis of the development of equations of state there are experimental data. However, when designing and characterising the thermodynamics of novel fluids, experimental data are not available because not yet measured or not physically measurable. More specifically to this work, the design and characterisation of novel reversible chemical reactions to be used as working fluids in thermodynamic cycles, as an alternative to N₂O₄, requires to set up a predictive methodology.

In this context, this paper presents a methodology, that we have implemented to the mentioned scope, and applied to the well-known system N₂O₄ ⇌ 2 NO₂. It is worth noting that, in the presence of these two molecules, this is the only reaction taking place between 300 K and 410 K. However, between 410 K and 1150 K a second chemical reaction, 2NO₂ ⇌ 2NO + O₂, occurs. In the interval of about 410 K – 600 K these two reactions take place simultaneously. In this preliminary work, the second reaction is not considered because we limit our discussion and calculations to a temperature range where the occurrence of this reaction is negligible. In this work, the definition of a thermodynamic model to enable the characterization of these systems is based on the complementary use of quantum chemistry, force field-based Monte Carlo simulations and equations of state. The only semi-empirical inputs to the methodology here presented are the force fields of NO₂ and of N₂O₄, which are taken from Bourasseau et al. [32] and, as mentioned above, are optimised to reproduce experimental data.

The structure of the paper is presented in the following. Section 2.1 introduces the singular thermodynamic behaviour of bi-component reversible reactions, such as N₂O₄ ⇌ 2NO₂, being characterized by phase equilibrium features similar to the ones of pure components, such as a unique critical and triple point, and a single vapour-liquid equilibrium (VLE) curve in the pressure–temperature space [36,37]. Sections 2.2 and 2.3 present, respectively, the set of equations applied to solve for chemical equilibrium in the presence of two coexisting phases and one chemical reaction, and the applied predictive equation of state. Section 2.4 introduces the way the ideal gas thermochemical properties (standard enthalpy of formation at 298.15 K, standard molar entropy at 298.15 K, heat capacities) and critical point coordinates of pure NO₂ and of N₂O₄ (required as an input to the equation of state) have been
calculated by Quantum Mechanic (QM) simulations and force field-based Monte Carlo simulations. Section 2.5 presents the methodology applied to perform force field-based Monte Carlo (MC) simulations on the studied reactive binary systems at VLE conditions, with a focus on the reactive N₂O₄ ⇌ 2NO₂, to characterize their phase equilibrium and volumetric properties of the coexisting vapour and liquid phases. Finally, section 3 presents the results obtained from QM and MC, and the final exploitation of those results in the application of equations of state, with the aim to characterize the full thermodynamics of binary reactive mixtures. Our findings are finally summarized in Section 4.

## 2. Methodology

Let us consider a binary system (A-B) undergoing the following general reversible reaction at its chemical equilibrium state:

$$ \nu_A(A)_{\bar{g}/l} \rightleftharpoons \nu_B(B)_{\bar{g}/l} \quad (1) $$

The implementation of mass, energy and entropy balances in thermodynamic cycles operating with reactive working fluids requires (i) the availability of an equation of state for mixtures, (ii) algorithms enabling the resolution of chemical equilibrium, considering the presence of chemical reactions and, eventually, of coexisting vapour and liquid phases, (iii) ideal-gas standard enthalpy of formation at 298.15 K, ideal-gas standard molar entropy at 298.15 K and ideal-gas molar heat capacity as a function of temperature for each of the components present in the mixture. The latter properties are needed to calculate the equilibrium constant at any temperature, while the ideal gas heat capacities are also used to calculate the enthalpy and/or the entropy of a system containing the reactive fluid with an equation of state.

### 2.1. The number of degrees of freedom of the system

By definition, the number of degrees of freedom of a system, or its variance, is the number of intensive independent variables of its phases which must be fixed, in order to characterize all the intensive variables of the phases of the whole system. It can be demonstrated that the variance, v, of a mixture composed by c-molecular species, r-chemical reactions and φ-fluid phases, is given by [38]:

$$ v = 2 - \varphi + c - r \quad (2) $$

If we consider binary and reactive systems such as the one considered in this paper, N₂O₄ ⇌ 2NO₂, i.e. a system involving only one reaction, it turns out that the variance of the system is equal to:

$$ \begin{cases} v = 2 & \text{if } \varphi = 1 \\ v = 1 & \text{if } \varphi = 2 \end{cases}, \quad \text{with } c = 2 \text{ and } r = 1 \quad (3) $$

This result coincides with the variance of a system composed by one single molecular species (where c = 1 and r = 0). The latter explains why the pressure-temperature relation of this system in vapour liquid equilibrium conditions is qualitatively similar to the vapor pressure of a pure component (Fig. 1). In addition, other singular features will be shown in section 3.2.2.

### 2.2. Fluid phase and reaction chemical equilibrium

Let us consider the chemical reaction reported in eq. (1). The chemical equilibrium condition states that the Gibbs energy of reaction of the system, ΔR_G, equals zero:

$$ \Delta_R G \triangleq \sum_{i=A,B} \nu_i \bar{g}_i = 0 \quad (4) $$

That is,

$$ \nu_B \bar{g}_B - \nu_A \bar{g}_A = 0 \quad (5) $$

In which $\bar{g}_i$ is the partial molar Gibbs energy (chemical potential) of

<figure><img src="image_3.png" /><figcaption>Fig. 1. Vapour-liquid equilibrium pressure (P<sup>VLE</sup>) as a function of temperature (T) for the reacting binary mixture N₂O₄ ⇌ 2NO₂. Blue points are experimental data from Scheffer and Treub [39] and Reamer and Sage [40]; red line and its end critical point (circle) are from DIPPR [41]. As expected, such a curve is very similar to the vaporization curve of a pure component.</figcaption></figure>

species i and vᵢ the algebraic stoichiometric coefficient. If the system is single phase, eq. (5) applies to that phase (liquid or vapour) and eq. (5) becomes:

$$ \nu_B \bar{g}_B^{liq} - \nu_A \bar{g}_A^{liq} = 0 \quad (6) $$

Or

$$ \nu_B \bar{g}_B^{vap} - \nu_A \bar{g}_A^{vap} = 0 \quad (7) $$

In the case the system is in vapor-liquid equilibrium (VLE), the fluid-phase equilibrium condition between the vapour and the liquid phase dictates that the partial molar Gibbs energy of each species is the same in all phases [38]:

$$ \begin{cases} \bar{g}_B^{liq} = \bar{g}_B^{vap} \\ \bar{g}_A^{liq} = \bar{g}_A^{vap} \end{cases} \quad (8) $$

As a consequence, equations (6) and (7) are simultaneously satisfied which demonstrates that the chemical reaction takes place in the two phases, simultaneously; indeed, it results:

$$ \Delta_R G^{liq} = \Delta_R G^{vap} = 0 \quad (9) $$

In other words, for a two-phase system, the condition ΔR_G<sup>vap</sup> = 0 is sufficient to impose that ΔR_G<sup>liq</sup> = 0, and vice-versa.

The expression of the Gibbs energy of reaction, ΔR_G, is recalled in the Supplementary Materials (section S1) as a function of:

* ideal gas properties (ideal gas standard molar enthalpy of formation at 298.15 K, ideal gas standard molar entropy at 298.15 K, and ideal gas isobaric heat capacity of each species of the mixture) obtained at the atomic level using ab-initio calculations

and

* real fluid properties (the fugacity of each species) determined by an equation of state.

The system of equations that needs to be solved to determine the chemical equilibrium of the binary reactive system in VLE can be written as follows:
$$ 
\left\{ 
\begin{aligned}
\bar{g}_{A}^{liq} &= \bar{g}_{A}^{vap} \\
\bar{g}_{B}^{liq} &= \bar{g}_{B}^{vap} \\
\Delta_{R}G &= 0 \text{ i.e. } K - \mathcal{P} = 0
\end{aligned}
\right. 
\qquad (10) $$

where, $\mathcal{P} \triangleq \prod_{i=A,B} (\hat{f}_i^\dagger)^{v_i}$, $\hat{f}_i$ is the fugacity of component i, $P^\circ$ is the standard pressure, $P^\circ = 1bar$, and K is the reaction equilibrium constant, related as follows to the standard Gibbs energy of reaction, $\Delta_R G^\circ$,

$$ \ln K(T) \triangleq -\frac{\Delta_R G^\circ(T)}{RT} \qquad (11) $$

In our case, the previous system of equations will be solved by an equation of state so that the uniformity of the chemical potential equation ($\bar{g}_i^{liq} = \bar{g}_i^{vap}$) is simply replaced by the isofugacity equation ($\hat{f}_i^{liq} = \hat{f}_i^{vap}$). However, in doing so, 2 new variables appear: the molar volumes of the liquid ($v^{liq}$) and vapor ($v^{vap}$) phases since the fugacity $\hat{f}_i$, derived from an equation of state, depends on temperature, composition and molar volume. It is therefore necessary to add 2 mathematical equations to reflect the fact that the variables T, P, v and composition are linked in each phase by the EoS. Finally, we get:

$$ 
\left\{ 
\begin{aligned}
\hat{f}_{A}^{liq}(T, v^{liq}, \mathbf{x}) &= \hat{f}_{A}^{vap}(T, v^{vap}, \mathbf{y}) \\
\hat{f}_{B}^{liq}(T, v^{liq}, \mathbf{x}) &= \hat{f}_{B}^{vap}(T, v^{vap}, \mathbf{y}) \\
K(T) - \mathcal{P} &= 0 \\
P - P_{EoS}(T, v^{liq}, \mathbf{x}) &= 0 \\
P - P_{EoS}(T, v^{vap}, \mathbf{y}) &= 0
\end{aligned}
\right. 
\qquad (12) $$

In eq. (12), $\mathcal{P}(T, v, z)$ can be calculated for either the liquid or gas phase since $\mathcal{P}^{vap} = \mathcal{P}^{liq}$. $x$ and $y$ denote the composition of the liquid $x = \{x_A, x_B = 1 - x_A\}$ and vapor phase $y = \{y_A, y_B = 1 - y_A\}$ respectively. $P_{EoS}(T, v, z)$ is the pressure of the system, calculated as a function of temperature, T, molar volume, v and molar composition, z, either of the liquid ($z = x$) or of the vapour phase ($z = y$) with a specifically selected EoS. The EoS model applied in this work is presented in the section 2.3. To conclude this section, we recall that the reactive system in VLE is monovariant. It is therefore necessary to specify an intensive variable before solving the system of equations. As an example, if the temperature is specified, the 5 equations make it possible to calculate $P, v^{liq}, v^{vap}, x_A$ and $y_A$.

## 2.3. The proposed predictive equation of state

The cubic Peng-Robinson equation of state [42], given by eq. (13), is used in this work to estimate the fugacities of each component:

$$ P(T, v, z) = \frac{R \cdot T}{v - b_m} - \frac{a_m}{v(v + b_m) + b_m(v - b_m)} \qquad (13) $$

where $z$ is the vector of the molar composition, $T$ is the temperature, $R$ is the universal gas constant and $a_m$ and $b_m$ are the mixing energy and co-volume parameters. The latter are calculated in this work with advanced mixing rules [43]:

$$ 
\left\{ 
\begin{aligned}
b_m &= \sum_{i=1}^{NC} z_i b_i \\
a_m &= \sum_{i=1}^{NC} z_i \frac{a_i}{b_i} + \frac{a_{res}^{E,y}}{\Lambda_{EoS}} \\
b_m &= \sum_{i=1}^{NC} z_i \frac{a_i}{b_i} + \frac{a_{res}^{E,y}}{\Lambda_{EoS}}
\end{aligned}
\right. 
\qquad (14) $$

where, $a_{res}^{E,y}$ is the residual part of an excess Helmholtz energy model calculated from an activity coefficient ($\gamma$) model, $\Lambda_{EoS}$ is a numerical parameter being a function of the considered equation of state [43], the pure component energy and co-volume parameters, i.e. $a_i$ and $b_i$, respectively, are calculated using eq. (15).

$$ 
\left\{ 
\begin{aligned}
& R = 8.314472 J \cdot mol^{-1} \cdot K^{-1} \\
& X = \left[ 1 + \sqrt[3]{4 - 2\sqrt{2}} + \sqrt[3]{4 + 2\sqrt{2}} \right]^{-1} \approx 0.253076587 \\
& b_i = \Omega_b \frac{RT_{c,i}}{P_{c,i}} \text{ with: } \Omega_b = \frac{X}{X+3} \approx 0.0777960739 \\
& a_i(T) = a_{c,i} \alpha_i(T) \text{ with: } 
\left\{ 
\begin{aligned}
& a_{c,i} = \Omega_a \frac{R^2 T_{c,i}^2}{P_{c,i}} \text{ and } \Omega_a = \frac{8(5X+1)}{49 - 37X} \approx 0.457235529 \\
& \alpha_i(T) = \left[ 1 + m_i \left( 1 - \sqrt{\frac{T}{T_{c,i}}} \right) \right]^2 \\
& \text{if } \omega_i \le 0.491 \text{ then } m_i = 0.37464 + 1.54226\omega_i - 0.26992\omega_i^2 \\
& \text{if } \omega_i > 0.491 \text{ then } m_i = 0.379642 + 1.48503\omega_i - 0.164423\omega_i^2 + 0.016666\omega_i^3
\end{aligned}
\right. \\
& (15)
\end{aligned}
\right. $$

More precisely, a zero-residual excess Helmholtz energy part, $a_{res}^{E,y}$, is considered for this work. In other words, the NO₂/N₂O₄ binary system is assumed to be athermal because the intermolecular forces between non-polar N₂O₄ molecules are weak, as well as between N₂O₄ and NO₂. Clearly, the interaction between NO₂ radicals is strong at low temperature and gives rise to the association reaction, which is indeed accounted for as being a chemical reaction in the system. The use of a zero residual excess Helmholtz energy leads to a model free of semi-empirical binary interaction parameters to be adjusted (otherwise present in the residual excess Helmholtz energy model) and will thus work in a predictive way.

## 2.4. The determination of thermochemical and thermophysical properties of pure fluids (N₂O₄ and NO₂)

It is worth recalling that while the use of an equation of state is sufficient to perform VLE calculations, the ideal gas heat capacity, $c_{p,i}^{ig}(T)$, of each species is needed for the determination of derived thermodynamic properties, such as the enthalpy or the entropy of a given phase [44]. Moreover, as recalled in eq. (S7) of the Supplementary Materials, the calculation of the ideal-gas standard enthalpy and entropy of reaction requires the knowledge – for each species – of the standard enthalpy of formation at 298.15 K, $\Delta_f H_{i,298.15K}$, the molar entropy at 298.15 K, $S_{i,298.15K}^e$, and the heat capacity as a function of temperature.

The temperature-dependent ideal gas heat capacity, the standard enthalpy of formation and the standard molar entropy of the pure components can be obtained by QM calculations, as explained in section 2.4.1, where the case of N₂O₄ ⇌ 2NO₂ is considered as an example.

Also, as shown in the previous section 2.3, the application of the Peng-Robinson equation of state (eqs. (13)-(15)) requires the knowledge of the temperature and pressure of the critical point of the two species forming the mixture, $T_{c,i}$ and $P_{c,i}$, and their acentric factors, $\omega_i$. Since the species forming any fast and reversible chemical reaction, i.e. N₂O₄ and NO₂ in the N₂O₄ ⇌ 2NO₂ system, do not exist in their pure form, the experimental measurement of such properties is physically not feasible. To overcome this problem, force field-based Monte Carlo simulations have been performed in this work, with two available MC softwares, Brick-Continuous-Fractional-Component-Monte-Carlo (CFCMC) [45] and Monte Carlo GIBBS [46], to determine P-ρ-T properties of the pure NO₂ and N₂O₄ and, thus, to deduce their critical point properties and acentric factors, as explained in section 2.4.2.

### 2.4.1. Quantum mechanics simulations

The thermochemical properties of NO₂ and N₂O₄, i.e. the ideal gas standard enthalpy of formation at 298.15 K, standard molar entropy at 298.15 K, and heat capacity, have been calculated using software Gaussian 09 (revision B.01) [47] and the GPOP post-processing software suite from Miyoshi [48].
QM calculations were performed at the CBS-QB3 level of theory [49]. A conformational analysis was performed for N₂O₄ at the B3LYP/6-31G (d) level [50], to determine the most stable conformer, that is the conformer corresponding to the lowest electronic energy. The conformational analysis consists in varying the considered dihedral angle by a step of 10° and re-optimizing all degrees of freedom except the dihedral angle, at each step. Also, internal modes are treated as harmonic oscillators in Gaussian; the conformational analysis enables to understand the validity of the Gaussian's harmonic oscillator assumption for internal rotation of moieties around simple bonds, such as the N-N bond in the N₂O₄ molecule. In general, if the conformational analysis shows that the highest energy gap between adjacent maximum and minimum energies, represented as an example in Fig. 2, is lower than R·T (where R is the universal gas constant), then the harmonic oscillator model is not valid because internal rotation can occur. In such a case, the results have to be corrected accounting for the hindered rotor model [51]. In this work, this treatment is done with GPOP and requires as an input the Gaussian's rotational hindrance potential. The post-processing is also performed with GPOP and enables the estimation of the standard enthalpy of formation and of the standard molar entropy.

To determine the equilibrium constant of the reactive system, and thus perform Monte Carlo simulations for the reactive system as shown in section 2.5, it is necessary to calculate the partition function of the isolated molecules. Those can be obtained from available thermochemistry data (JANAF tables) or from quantum mechanical calculations (with Gaussian), from which it is possible to determine the atomisation energy and the rotational and vibrational temperatures used as an input to the calculation of the partition function of the isolated molecules.

### 2.4.2. Monte Carlo simulations with Brick-CFCMC software

To perform Monte Carlo simulations of N₂O₄ and NO₂, a classical force field is used in which interatomic interactions are described by both Lennard-Jones (LJ) and electrostatic interactions. NO₂ is treated as a rigid molecule. Concerning the N₂O₄ molecule, this molecule is rigid as well, except when the rotation of NO₂ moieties is considered and the rotational hindrance potential is then required as additional input for the MC simulations. The potential is deduced by optimising the parameters of a cosine series, on the rotation energy scan obtained with Gaussian.

Lennard-Jones interactions are truncated at 12 Å and analytic tail corrections are used [52]. Electrostatics are computed using the Ewald

<figure><img src="image_4.png" /><figcaption>Fig. 2. Example of correlated torsion potential energy obtained in Gaussian by varying the dihedral angle of the molecule.</figcaption></figure>

summation technique with a relative precision of 10⁻⁵. The force field parameters of NO₂ and N₂O₄ are taken from Bourasseau et al. [32] and are also listed in Table S1 of the Supplementary Materials. For atoms that are not identical, the Lorentz-Berthelot mixing rules [52,53] are used. The intermolecular coordinates of NO₂ and N₂O₄ are respectively listed in Table S2 and Table S3 of the Supplementary Materials. Three-dimensional systems with periodic boundary conditions are considered in all simulations.

Pure component Gibbs ensemble simulations using the Continuous Fractional technique [54–56] were performed using the Brick-CFCMC software [45,57]. In the constant volume-constant temperature (NVT) version of the Gibbs ensemble [52,53], there are two simulation boxes that can exchange volume, energy, and molecules, in such a way that the simulation boxes resemble the coexisting phases at equilibrium. Molecule transfers between both simulation boxes are facilitated using a single fractional molecule per component [58]. The interactions of such a fractional molecule are determined by an order parameter λ, with λ = 0 when the fractional molecule has no interactions with the surrounding molecules, and λ = 1 when the fractional molecule has full, unaltered, interactions with the surrounding molecules. Besides thermalization trial moves (molecule translations, molecule rotations, and volume changes), the simulation requires three additional trial moves [58]: (1) attempts to change the value of λ; (2) attempts to transfer a fractional molecule from one simulation box to the other one; (3) attempts to transform the fractional molecule into a whole molecule, and simultaneously, transform a “whole” molecule in the other simulation box into a fractional molecule, at the same value of λ. Biasing using the Wang-Landau algorithm [59] is needed to ensure that the probability distribution of λ is flat in both boxes, and that fractional molecules are equally likely to be in the simulation boxes. The chemical potential of component i in simulation box A follows directly from [45,57,60]:

$$ \mu_i^A = k_B T \left( - \ln \left( \frac{q_i}{\Lambda_i^3 \rho_0} \right) + \ln \left( \frac{\rho_i^A}{\rho_0} \right) - \ln \left( \frac{p(\lambda_i^A = 1)}{p(\lambda_i^A = 0)} \right) \right) \quad (16) $$

in which k_B is the Boltzmann constant, T is the absolute temperature, q_i is the intramolecular part of the isolated molecule partition function of molecule type i (excluding the translational part), Λ_i is the thermal wavelength of molecule type i, ρ_0 is an arbitrary reference density set to 1 molecule per cubic Angstrom, ρ_i^A is the number density of component i in simulation box A, and p(λ_i^A) is the Boltzmann probability distribution of the value of λ of component i in simulation box A. The term q_i/Λ_i^3 can be obtained either from QM simulation or thermochemical tables (e.g. the JANAF tables), see Refs. [45,57,61] for details. At phase equilibrium, chemical potentials of component i are equal in both simulation boxes. As the Brick-CFCMC software does not compute pressures, a series of simulations in the constant temperature-pressure (NPT) ensemble are carried out around the density of the gas phase box. The coexistence pressure is calculated by matching the gas phase density to the NPT simulations via interpolation. For more details on the simulation technique and the derivation of the accepted probabilities of the trial moves, the reader is referred to Refs. [45,57,58].

### 2.4.3. Monte Carlo simulations with GIBBS software

In this paper, other Monte Carlo results are presented for N₂O₄ and NO₂. These molecular simulations have been performed using NVT Gibbs ensemble Monte Carlo (MC) [62], with the “GIBBS Monte Carlo” code [46], to determine the saturation densities and pressures of pure N₂O₄ and pure NO₂. In this MC algorithm, the sampling of the configurational space is ensured by three Monte Carlo moves: rigid body translations or rotations of molecules, and volume changes. More specifically, differently with respect to the calculations performed in this work with Brick-CFCMC software, NO₂ and N₂O₄ molecules have been both represented as rigid molecules. Thus, in these simulations, N₂O₄ is not allowed to make torsional rotations, while in Brick the torsional potential is modelled as determined from Quantum Mechanics
simulations. However, the force field parameters in these simulations coincide with the ones from Bourasseau et al. [32], like in settings used in Brick. In all simulations performed with GIBBS code, a cutoff for LJ intermolecular interactions equal to the half of the simulation box length has been used to reduce the computing time. Beyond the cutoff, standard long-range corrections have been employed. For more details on the simulation technique or on the molecular models used to determine these latter MC results, the reader is invited to refer to some previous works [32,63–65].

### 2.4.4. The derivation of critical coordinates and acentric factors

The determination of the critical point properties by MC is not straightforward, since the Gibbs ensemble is not applicable in the vicinity of the critical region [52]. To localize the critical point, we have proceeded as described below.

Once P-ρ-T properties in saturation conditions are calculated by MC, the critical temperature $T_{c,i}$ of the pure compounds, NO₂ and N₂O₄, are obtained by extrapolation from the sub-critical simulation points, assuming the following scaling law:

$$ \rho_i^{\text{liq}}(T) - \rho_i^{\text{vap}}(T) = B_i(T_{c,i} - T)^{\beta} \quad (17) $$

where $\rho_i^{\text{liq}}(T)$ and $\rho_i^{\text{vap}}(T)$ are the densities of component i at temperature T in the liquid and the vapor phases, respectively, calculated with Monte Carlo simulations, while $\beta$ is the critical exponent, equal to 0.32 +/- 0.01 [66]. The critical temperature of each i-th component, $T_{c,i}$, and its parameter $B_i$, also component-dependent, have been regressed by least square method minimizing the following objective function, O.F.:

$$ O.F. : \min_{T_{c,i}, B_i} \left\{ \sqrt{\sum_{k=1}^{NP} \left( \rho_k^{\text{liq,MC}}(T_k) - \rho_k^{\text{vap,MC}}(T_k) - B_i(T_{c,i} - T_k)^\beta \right)^2} \right\} \quad (18) $$

where NP is the number of data points. The critical density, $\rho_{c,i}$, is then calculated with the application of the law of rectilinear diameter,

$$ \frac{\rho_i^{\text{liq}}(T) + \rho_i^{\text{vap}}(T)}{2} = \rho_{c,i} - A_i(T_{c,i} - T) \quad (19) $$

and the least square optimization,

$$ O.F. : \min_{\rho_{c,i}, A_i} \left\{ \sqrt{\sum_{k=1}^{NP} \left( \frac{\rho_k^{\text{liq,MC}}(T_k) + \rho_k^{\text{vap,MC}}(T_k)}{2} - \rho_{c,i} + A_i(T_{c,i} - T_k) \right)^2} \right\} \quad (20) $$

Moreover, from the calculated P-T saturation points of the two pure fluids, a correlation $P_i^{\text{sat}}(T)$ is determined

$$ \ln P_i^{\text{sat}}(T) = \alpha_i + \frac{\beta_i}{T} \quad (21) $$

where, parameters $\alpha_i$ and $\beta_i$ are determined by minimizing:

$$ O.F. : \min_{\alpha_i, \beta_i} \left\{ \sqrt{\sum_{k=1}^{NP} \left( \exp\left(\alpha_i + \frac{\beta_i}{T_k}\right) - P_i^{\text{sat}}(T_k) \right)^2} \right\} \quad (22) $$

and $P_{c,i}$ can be finally calculated from the optimized correlations, eq. (21), as a function of $T_{c,i}$.

Once $T_{c,i}$ and $P_{c,i}$ are obtained, the acentric factor $\omega_i$ is calculated from its definition:

$$ \omega_i = -\log_{10} P_{r,i}^{\text{sat}} - 1 \quad (23) $$

where $P_{r,i}^{\text{sat}}$ is the reduced vapor pressure ($P_i^{\text{sat}}/P_{c,i}$) at a reduced temperature ($T/T_{c,i}$) equal to 0.7.

### 2.5. The determination of thermodynamic properties of the reactive N₂O₄ ⇌ 2 NO₂

To study the VLE phase coexistence of the reactive system N₂O₄ ⇌ 2 NO₂, Monte Carlo simulations in the reactive Gibbs Ensemble [33,34,60,67,68] are performed with Brick-CFCMC [45]. In these simulations, the constant-volume Gibbs Ensemble setup of section 2.4.2 is coupled to a grand-canonical reservoir of molecules of NO₂ and N₂O₄. The imposed chemical potentials of NO₂ and N₂O₄ are set such that the system is at chemical equilibrium so $\mu_{N_2O_4} = 2 \cdot \mu_{NO_2}$ and the exchanges from the reservoir are such that the total number of atoms of each type in the system is conserved. A fractional group is added to the vapor box of the Gibbs Ensemble, which contains fractional molecules of either the reactants or the reaction products [60]. Note that this group of fractional molecules is completely independent from the fractional molecules required for molecule transfers between the simulation boxes (as used in the CFCMC version of the Gibbs Ensemble described in section 2.4.2). In addition to the trial moves of the CFCMC version of the Gibbs Ensemble (section 2.4.2), the following trial moves (which all operate in the vapor box) are used for the molecules in the fractional group: (1) attempts to change the order parameter $\lambda_R$ for the molecules in the fractional group. $\lambda_R = 0$ means that the molecules in the fractional group have no interaction with the surrounding molecules, and $\lambda_R = 1$ means that the molecules in the fractional group have full, unaltered interactions with the surrounding molecules; (2) attempts to transform the molecules in the fractional group from reactants into reaction products, or vice-versa. This requires the terms $q_i/\Lambda_i^3$ for NO₂ and N₂O₄; (3) attempts to transform the fractional molecules in the group into “whole” molecules, and simultaneously, transform “whole” molecules into a fractional group, in such a way that the fractional group changes from reactants to reaction products (or vice versa). In this trial move, the value of $\lambda_R$ is unchanged. Biasing is needed to ensure that it is equally likely that the fraction group contains a single fractional molecule of N₂O₄, or two fractional molecules of NO₂, and that each value of $\lambda_R$ is equally likely. At combined chemical equilibrium and phase equilibrium, we have $\mu_{N_2O_4} = 2 \cdot \mu_{NO_2}$ for both simulation boxes A and B, and $\mu_i^A = \mu_i^B$ for both NO₂ and N₂O₄. The coexistence pressure is obtained by performing a series of NPT simulations at the average composition of the gas phase box, and use interpolation to match the simulated gas phase density from the Gibbs ensemble simulation in the reaction ensemble. For more details about the simulation setup and the derivation of the acceptance rules, the reader is referred to Refs. [45,60]. The isolated molecule partition function for both NO₂ and N₂O₄ were obtained in this work both from QM simulations, by using Gaussians outputs reported in Table S4 of the Supplementary Materials, and the thermochemical results reported in JANAF tables [45]. The calculated partition functions are shown in Fig. 3, where $V_0 = 1 \AA^3$. Two correlations that correlate the JANAF values are proposed below,

$$ \ln \left( \frac{qV_0}{\Lambda^3} \right)_{NO_2} = \frac{110983.613}{T} + 16.403 \quad (24) $$

$$ \ln \left( \frac{qV_0}{\Lambda^3} \right)_{N_2O_4} = \frac{228591.177}{T} + 23.218 \quad (25) $$

## 3. Results

This section presents the results obtained from molecular simulations, as well as the calculations performed with the cubic equation of state and enabled by the calculated equation of state’s input critical temperature, pressure and acentric factor of the pure fluids NO₂ and N₂O₄.
<figure><img src="image_5.png" /><figcaption>Fig. 3. Natural logarithm of the partition functions, q, divided by the thermal wavelength, for NO₂ (fig. (a)) and N₂O₄ (b). Data obtained from the JANAF tables. Correlations are eqs. (24) and (25).</figcaption></figure>

## 3.1. Thermochemical and thermophysical properties of N₂O₄ and NO₂

### 3.1.1. Quantum Mechanics results

Fig. 4 shows the geometries of NO₂ and N₂O₄, optimised with Gaussian as discussed in section 2.4.1. It is worth noting that the length of N = O bonds and the N = O = N angles are very similar in both NO₂ and N₂O₄, which indicates that very limited rearrangements occur in moieties when N₂O₄ splits into two NO₂ moieties. The main specificity of N₂O₄ is the length of the N-N bond which is particularly long for such a bond: 1.797 Å compared to 1.487 Å for the N-N bond in hydrazine (value obtained at the CBS-QB3 level of theory with software Gaussian) and which can be attributed to an inductive effect of the vicinal oxygen atoms. This has consequently an impact on the energy of that particular bond as discussed further. Fig. 5 shows the result of the conformational analysis that was performed to verify that the optimized geometry obtained for N₂O₄ corresponds well to the conformer with the lowest energy. The optimal coefficients for the torsion potential N₂O₄ are shown in Table S5 of the Supplementary Materials. As expected, two distinct stable conformations exist for N₂O₄: a planar one with the lowest energy, and a skew one with an energy ca. 30 kJ/mol higher than that one the planar conformer. As it will be confirmed by MC simulation results considering the torsion of N₂O₄ (see section 3.1.2), the calculated potential barrier shows that the actually observed conformations are plane or almost plane. Some studies [69] focused on the origin of the greater stability of the planar conformer, but this is beyond the scope of the present work.

As enthalpies provided by the Gaussian software are relative to a reference that is specific to the method used, they need to be rescaled. For NO₂, we considered in the present work the fictive reaction (2) NO₂ + N₂

<figure><img src="image_6.png" /><figcaption>Fig. 5. Conformational analysis of N₂O₄. The points in the figure represent Gaussian results, while the line is the cosine series fitted to the QM results.</figcaption></figure>

≈ 4 NO leading to $\Delta_f H_{i,298.15\text{K}}$ equal to 34.20 kJ/mol, which compares well to the data available in the NIST database [70] (33.1 kJ/mol). For N₂O₄, the fictive reaction considered was N₂O₄ ⇌ 2 NO₂ leading to $\Delta_f H_{i,298.15}$ equal to 9.99 kJ/mol, to be compared with the value of 9.08 kJ/mol provided in the NIST database. This leads to an enthalpy of reaction of 56.2 kJ/mol for

<figure><img src="image_7.png" /><figcaption>Fig. 4. Optimal geometry of NO₂ and N₂O₄. Atomic coordinates are reported in Table S2 and S3 the Supporting Information.</figcaption></figure>
the reaction of dissociation of N₂O₄ into two NO₂ moieties (57.12 kJ/mol using the NIST data), confirming that the value of the N-N bond is usually low and in line with the bond length of 1.797 Å.

For the ideal-gas standard molar entropy, Ṡᵢ₂₉₸.₁₅K, assuming an external symmetry number of 2 for NO₂, the value derived with post-processing Gaussian results with GPOP, as briefly described in section 2.4.1, is 240.03 J/mol/K which compares very well to the data in the NIST database (240.04 J/mol/K). For N₂O₄, assuming an external symmetry number of 4, the calculated standard entropy is 297.65 J/mol/K, which is a bit lower than the data provided in the NIST database (304.38 J/mol/K) but still within the acceptable range. At such a low temperature, the difference in entropy for the –NO₂ rotor is only 0.54 J/mol/K between the harmonic oscillator approximation and the hindered rotor one, as a consequence of the fact that the highest energy gap between adjacent maximum and minimum energy (i.e., the torsional barrier) is higher than R·T = 2.47 kJ/mol at 298.15 K and thus the harmonic oscillator assumption is still valid (see section 2.4.1). The results obtained are summarised in Table 1.

Another result obtained from the postprocessing step is the ideal gas isobaric specific heat capacity, which has been correlated to the DIPPR correlation n. 107:

$$c_{P,i}^{ig}(T) = c_0 + c_1 \left( \frac{c_2/T}{\sinh(c_2/T)} \right)^2 + c_3 \left( \frac{c_4/T}{\cosh(c_4/T)} \right)^2 \quad (26)$$

The coefficients optimised in this work for NO₂ and N₂O₄ are reported in Table 2 and shown in Fig. S1 and S2 of the Supplementary Materials.

### 3.1.2. Monte Carlo results

The MC simulations performed with Brick-CFCMC enabled the computation of coexisting densities and saturation pressures of the two pure compounds, NO₂ and N₂O₄, at different temperatures. Results are presented in Table 3 and Table 4 as well as in Figs. 6–9. MC simulations of N₂O₄ have shown that the most probable conformer is the planar one, according to the potential barrier calculated in Fig. 5.

More precisely, Fig. 6 and Fig. 7 show the densities of the liquid and vapour phase of, respectively, NO₂ and N₂O₄, as a function of temperature. A comparison between the results of the two different sets of simulations show an excellent agreement between the two methods. Main quantitative differences between the results of Brick-CFCMC and GIBBS may derive from slightly different implementation of the Ewald summation, or from a different cutoff and from the additional consideration (in Brick-CFCMC calculations) of torsion potential of N₂O₄. Moreover, these figures show the predicted critical points (T<sub>c,bi</sub> ρ<sub>c,i</sub>) of both NO₂ and N₂O₄, optimized over the two sets of data. It can be observed that despite simulations performed with Brick-CFCMC did not converge close to the critical point, the prediction of the critical point obtained from this set of data is acceptably similar to the one obtained with GIBBS code which better converges close to the critical point. It is highlighted that the uncertainties of MC simulations from Brick-CFCMC and GIBBS are not shown in Fig. 6 and Fig. 7, for the sake of simplicity, but are reported in Table 3 and in Table 4.

The calculated saturation pressures for the pure components are shown in Fig. 8 and Fig. 9, for NO₂ and N₂O₄ respectively. As described in section 2.4.4, the correlation in eq. (21) has been optimized for both components, over the experimental data obtained by Brick-CFCMC, and this is also reported in these figures. Coefficients of the correlations are

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

Ideal gas standard properties at 298.15 K, obtained in this work by quantum mechanics.

<table><caption>**Table 2**</caption><thead><tr><th></th><th>NO<sub>2</sub></th><th>N<sub>2</sub>O<sub>4</sub></th></tr></thead><tbody><tr><td>Δ<sub>f</sub>H<sub>i,298.15</sub> (kJ/mol)</td><td>34.20</td><td>9.99</td></tr><tr><td>S<sub>i,298.15K</sub> (J/mol/K)</td><td>240.03</td><td>297.65</td></tr></tbody></table>

Coefficients for the calculation of the isobaric specific heat capacity of NO₂ and N₂O₄ by eq. (26).

<table><caption>**Table 3**</caption><thead><tr><th></th><th>c<sub>0</sub> [J/mol/K]</th><th>c<sub>1</sub> [K]</th><th>c<sub>2</sub> [K]</th><th>c<sub>3</sub> [K]</th><th>c<sub>4</sub> [K]</th></tr></thead><tbody><tr><td><strong>NO<sub>2</sub></strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>QM, this work</td><td>33.631</td><td>24.566</td><td>1.1684</td><td>10.4065</td><td>0.6063</td></tr><tr><td>DIPPR [41]</td><td>33.261</td><td>24.919</td><td>1.1122</td><td>9.2534</td><td>0.5592</td></tr><tr><td><strong>N<sub>2</sub>O<sub>4</sub></strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>QM, this work</td><td>67.574</td><td>56.723</td><td>0.72373</td><td>13.722</td><td>2.0934</td></tr><tr><td>DIPPR [41]</td><td>65.796</td><td>58.428</td><td>0.69602</td><td>11.887</td><td>2.0748</td></tr></tbody></table>

Results of MC simulations obtained for pure NO₂ and N₂O₄, obtained in this work by the use of Brick-CFCMC. “u” denotes the uncertainty of the calculated variables. These uncertainties correspond to the standard deviation between the results of MC simulations performed at each temperature to calculate saturation densities and pressures.

<table><caption>**Table 4**</caption><thead><tr><th>T<br>[K]</th><th>p<sup>sat</sup><br>[bar]</th><th>u(P)<br>[bar]</th><th>ρ<sub>liq</sub><br>[kg/m<sup>3</sup>]</th><th>ρ<sub>vap</sub><br>[kg/m<sup>3</sup>]</th><th>u(ρ<sub>liq</sub>)<br>[kg/m<sup>3</sup>]</th><th>u(ρ<sub>vap</sub>)<br>[kg/m<sup>3</sup>]</th></tr></thead><tbody><tr><td colspan="7"><strong>NO<sub>2</sub></strong></td></tr><tr><td>240</td><td>25.07</td><td>0.45</td><td>985.6</td><td>74.9</td><td>1.3</td><td>1.9</td></tr><tr><td>245</td><td>28.60</td><td>0.31</td><td>960.9</td><td>86.6</td><td>1.6</td><td>1.5</td></tr><tr><td>250</td><td>32.82</td><td>0.84</td><td>935.0</td><td>101.4</td><td>2.3</td><td>4.1</td></tr><tr><td>255</td><td>37.11</td><td>0.52</td><td>905.5</td><td>116.8</td><td>1.9</td><td>2.9</td></tr><tr><td>260</td><td>42.42</td><td>1.01</td><td>875.4</td><td>138.3</td><td>3.5</td><td>6.2</td></tr><tr><td>265</td><td>47.39</td><td>1.33</td><td>838.7</td><td>158.6</td><td>3.5</td><td>8.9</td></tr><tr><td>270</td><td>52.10</td><td>0.92</td><td>791.7</td><td>177.3</td><td>6.2</td><td>7.0</td></tr><tr><td colspan="7"><strong>N<sub>2</sub>O<sub>4</sub></strong></td></tr><tr><td>260</td><td>0.107</td><td>0.020</td><td>1537.6</td><td>0.454</td><td>1.9</td><td>0.085</td></tr><tr><td>270</td><td>0.175</td><td>0.036</td><td>1516.5</td><td>0.72</td><td>1.3</td><td>0.15</td></tr><tr><td>280</td><td>0.332</td><td>0.026</td><td>1494.1</td><td>1.31</td><td>1.1</td><td>0.10</td></tr><tr><td>290</td><td>0.457</td><td>0.038</td><td>1472.2</td><td>1.75</td><td>0.7</td><td>0.15</td></tr><tr><td>300</td><td>0.752</td><td>0.101</td><td>1449.0</td><td>2.80</td><td>0.5</td><td>0.38</td></tr><tr><td>310</td><td>1.046</td><td>0.093</td><td>1424.8</td><td>3.79</td><td>1.1</td><td>0.34</td></tr><tr><td>320</td><td>1.55</td><td>0.12</td><td>1401.8</td><td>5.47</td><td>1.1</td><td>0.43</td></tr><tr><td>330</td><td>2.23</td><td>0.20</td><td>1377.4</td><td>7.71</td><td>1.0</td><td>0.72</td></tr><tr><td>340</td><td>2.90</td><td>0.17</td><td>1352.3</td><td>9.85</td><td>0.3</td><td>0.63</td></tr><tr><td>350</td><td>4.11</td><td>0.14</td><td>1326.0</td><td>13.8</td><td>1.0</td><td>0.52</td></tr><tr><td>360</td><td>5.04</td><td>0.12</td><td>1299.3</td><td>16.6</td><td>0.7</td><td>0.45</td></tr><tr><td>370</td><td>6.67</td><td>0.13</td><td>1269.40</td><td>21.8</td><td>0.73</td><td>0.48</td></tr><tr><td>380</td><td>8.51</td><td>0.31</td><td>1240.68</td><td>27.6</td><td>2.40</td><td>1.13</td></tr><tr><td>390</td><td>10.63</td><td>0.35</td><td>1210.22</td><td>34.5</td><td>1.50</td><td>1.30</td></tr><tr><td>400</td><td>13.17</td><td>0.38</td><td>1176.49</td><td>42.4</td><td>1.66</td><td>1.48</td></tr><tr><td>410</td><td>16.55</td><td>0.18</td><td>1141.44</td><td>52.9</td><td>0.99</td><td>0.71</td></tr><tr><td>420</td><td>19.76</td><td>0.35</td><td>1098.91</td><td>63.2</td><td>1.39</td><td>1.47</td></tr><tr><td>430</td><td>22.97</td><td>0.40</td><td>1050.68</td><td>73.8</td><td>5.37</td><td>1.67</td></tr><tr><td>440</td><td>26.54</td><td>0.68</td><td>989.79</td><td>85.8</td><td>6.12</td><td>2.89</td></tr></tbody></table>

Results of MC simulations obtained for pure NO₂ and N₂O₄, obtained with GIBBS code. “u” denotes the uncertainty of the calculated variables. These uncertainties correspond to the standard deviation between the results of MC simulations performed at each temperature to calculate saturation densities and pressures.

<table><thead><tr><th>T<br>[K]</th><th>p<sup>sat</sup><br>[bar]</th><th>u(P)<br>[bar]</th><th>ρ<sub>liq</sub><br>[kg/m<sup>3</sup>]</th><th>ρ<sub>vap</sub><br>[kg/m<sup>3</sup>]</th><th>u(ρ<sub>liq</sub>)<br>[kg/m<sup>3</sup>]</th><th>u(ρ<sub>vap</sub>)<br>[kg/m<sup>3</sup>]</th></tr></thead><tbody><tr><td colspan="7"><strong>NO<sub>2</sub></strong></td></tr><tr><td>220</td><td>13.1</td><td>0.29</td><td>1087.5</td><td>38.7</td><td>1.02</td><td>2.77</td></tr><tr><td>230</td><td>17.7</td><td>0.4</td><td>1043.1</td><td>51.9</td><td>1.49</td><td>3.43</td></tr><tr><td>240</td><td>24.1</td><td>0.43</td><td>1001.8</td><td>71.8</td><td>1.83</td><td>3.98</td></tr><tr><td>250</td><td>31.2</td><td>0.53</td><td>947.9</td><td>94.4</td><td>2.54</td><td>4.86</td></tr><tr><td>260</td><td>40.3</td><td>0.68</td><td>889.3</td><td>127</td><td>3.42</td><td>6.13</td></tr><tr><td>270</td><td>51.0</td><td>0.84</td><td>819.4</td><td>174</td><td>4.83</td><td>6.84</td></tr><tr><td>275</td><td>56.5</td><td>0.73</td><td>765.8</td><td>200</td><td>5.91</td><td>9.41</td></tr><tr><td>280</td><td>63.8</td><td>0.66</td><td>708.9</td><td>242.7</td><td>7.31</td><td>18.53</td></tr><tr><td colspan="7"><strong>N<sub>2</sub>O<sub>4</sub></strong></td></tr><tr><td>300</td><td>0.7</td><td>0.04</td><td>1462.8</td><td>2.7</td><td>0.17</td><td>5.49</td></tr><tr><td>350</td><td>4.4</td><td>0.11</td><td>1346.0</td><td>15</td><td>0.45</td><td>3.7</td></tr><tr><td>400</td><td>13.6</td><td>0.46</td><td>1202.2</td><td>44.4</td><td>1.71</td><td>5.11</td></tr><tr><td>450</td><td>35.9</td><td>1.05</td><td>1018.8</td><td>129.4</td><td>6.12</td><td>8.57</td></tr><tr><td>475</td><td>52.8</td><td>1.39</td><td>848.2</td><td>219.8</td><td>11.62</td><td>30.35</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Liquid-vapor coexisting densities of pure NO₂ obtained by MC simulations in this work. Crosses denote the results obtained with Brick-CFCMC as described in section 2.4.2 and reported in Table 3, while points denote the results obtained with GIBBS reported in Table 4. Black symbols represent the critical point predicted as detailed in section 2.4.4, using the two sets of data.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Liquid-vapor coexisting densities of pure N₂O₄ obtained by MC simulations in this work. Crosses denote the results obtained with Brick-CFCMC as described in section 2.4.2 and reported in Table 3, while points denote the results obtained with GIBBS reported in Table 4. Black symbols represent the critical point predicted as detailed in section 2.4.4, using the two sets of data.</figcaption></figure>

shown in Table S6 of the Supplementary Materials. With those correlations, the critical pressure could then be determined (see section 2.4.4), and this is reported in these figures too.

The critical coordinates and acentric factors calculated as described in section 2.4.4 are reported in Table 5. It can be observed that the values determined in this work are very different with respect to the critical coordinates previously obtained by McCarty et al. [26], also reported in Table 5, who made the assumption of equal acentric factor between the two molecules.

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Vapor pressures of pure NO₂ obtained by MC simulations and the resulting correlations. The empty circles represent the critical point predicted as detailed in section 2.4.4, using the set of data obtained in this work.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Vapor pressures of pure N₂O₄ obtained by MC simulations and the resulting correlation. The empty circles represent the critical point predicted as detailed in section 2.4.4, using the set of data obtained in this work.</figcaption></figure>

**Table 5**

Critical coordinates and acentric factor of NO₂ and N₂O₄ obtained from the molecular simulation results.

<table><thead><tr><th></th><th>T<sub>c</sub> [K]</th><th>ρ<sub>c</sub> [kg/m<sup>3</sup>]</th><th>P<sub>c</sub> [bar]</th><th>ω [-]</th></tr></thead><tbody><tr><td>NO<sub>2</sub> / This work, from results of MC Brick</td><td>282.2</td><td>471.80</td><td>68.2</td><td>0.0565</td></tr><tr><td>N<sub>2</sub>O<sub>4</sub> / This work, from results of MC Brick</td><td>484.2</td><td>502.44</td><td>55.9</td><td>0.3212</td></tr><tr><td>NO<sub>2</sub> / This work, from results of MC GIBBS</td><td>285.0</td><td>471.20</td><td>69.9</td><td>0.0611</td></tr><tr><td>N<sub>2</sub>O<sub>4</sub> / This work, from results of MC GIBBS</td><td>488.6</td><td>524.79</td><td>66.7</td><td>0.3306</td></tr><tr><td>NO<sub>2</sub> / Results from McCarty et al. [26]</td><td>239.3</td><td>-</td><td>103.3</td><td>0.01413</td></tr><tr><td>N<sub>2</sub>O<sub>4</sub> / Results from McCarty et al. [26]</td><td>547.5</td><td>-</td><td>221.0</td><td>0.01413</td></tr></tbody></table>
## 3.2. Thermodynamic properties of N₂O₄ ⇌ 2 NO₂

### 3.2.1. Monte Carlo results

The thermodynamic properties of the reactive mixture determined by MC simulations with Brick-CFCMC, as detailed in section 2.5, are presented in Table 6 and in Figs. 10–12. Fig. 10 shows that the calculated VLE pressure of the reactive system is in very good agreement with available experimental data and with the correlation provided by DIPPR. The same consideration is valid for the other calculated properties: the composition of the vapour and liquid phase, in Fig. 11, which is compatible with the MC results obtained by Bourasseau et al. [32], and the T - ρᵢᵗᵢ - ρᵣᵗᵣ data points reported in Fig. 12, are also in good agreement with the experimental data.

It is worth noting that the P-{*x,y*} diagram reported in Fig. 11 has two unique features that derive from the mono-variance of this binary reactive system in VLE (see section 2.1): (i) The diagram is a *unique* and *non-isothermal* P-{*x,y*} diagram for this system; indeed, at each specified pressure (between the triple point pressure and the critical pressure of the reactive system) corresponds a different temperature and different mole fractions in each phase. (ii) The diagram is not closed at the lowest pressure (the triple point pressure of the reactive fluid), because the condition of pure N₂O₄ (x,yNO₂ = 0) is not compatible with the occurrence of the chemical reaction, whose extent of reaction is non-zero at those pressure levels.

### 3.2.2. Thermodynamic calculations

The determination of the ideal gas standard enthalpy of formation (at 298.15 K), standard molar entropy (at 298.15 K) and isobaric heat capacity (as a function of temperature) of NO₂ and of N₂O₄, by quantum mechanics (see sections 2.4.1 and 3.1.1), as well as their critical point coordinates and acentric factors by force-field based Monte Carlo simulations (see sections 2.4.2, 2.4.4, and 3.1.2) allows the further application of a predictive thermodynamic model to fully characterize the properties of this system: the equation of state presented in section 2.3.

The first computed results are shown in Fig. 13, where the global phase equilibrium diagram (GPED) [76] is represented for this system, and is calculated with two sets of critical coordinates: the solid lines represent the results of the equation of state with the critical coordinates obtained from MC simulations of Brick-CFCMC (see Table 5), also called in the following “EoS 1”, while the dashed lines are results of calculations performed with the critical coordinates obtained from MC data obtained from GIBBS, also called in the following “EoS 2”. Moreover, the GPED shows:

<figure><img src="image_12.png" /><figcaption>Fig. 10. VLE pressures of the coexisting vapour and liquid phases of the reactive mixture containing N₂O₄ and NO₂. Crosses are points calculated by MC with Brick-CFCMC. Data of Scheffer and Treub [39] and Reamer and Sage [40] are also shown for comparison.</figcaption></figure>

* the vaporisation curves of the pure NO₂, ending at their critical point (represented by the circle mark) – pink lines;

* the vapourisation curves of the pure N₂O₄, ending at their critical point – black lines;

* the locus of critical points of the binary non-reactive mixture formed by N₂O₄ and NO₂ – green lines;

* the VLE pressure of the reactive binary mixture, which is a line instead of a surface because of the mono-variance of this system (see section 2.1), ending at their critical point – blue lines; moreover, the red line and the light blue points on it represent, respectively, the DIPPR correlation and the experimental VLE pressures for the reactive system already shown in Fig. 1. Black crosses are the MC results obtained with Brick-CFCMC (see Table 6 and Fig. 10).

* the triple temperature of the system – the light blue dashed line.

Table 6

Results of MC simulations obtained for reactive N₂O₄ ⇌ 2 NO₂, obtained in this work by the use of Brick-CFCMC. “u” denotes the uncertainty of the calculated variables. These uncertainties correspond to the standard deviation between the results of MC simulations performed at each temperature to calculate saturation densities and pressures. Molar fraction uncertainties are determined from the calculated standard deviations of molecules number in the simulation boxes resulting from each simulation, by applying the propagation error law.

<table><thead><tr><th>T<br>[K]</th><th>p<sub>sat</sub><br>[bar]</th><th>u(P)<br>[bar]</th><th>ρ<sub>liq</sub><br>[kg/m<sup>3</sup>]</th><th>ρ<sub>vap</sub><br>[kg/m<sup>3</sup>]</th><th>u(ρ<sub>liq</sub>)<br>[kg/m<sup>3</sup>]</th><th>u(ρ<sub>vap</sub>)<br>[kg/m<sup>3</sup>]</th><th>x<sub>NO<sub>2</sub></sub><br>[-]</th><th>y<sub>NO<sub>2</sub></sub><br>[-]</th><th>u(x)<br>[-]</th><th>u(y)<br>[-]</th></tr></thead><tbody><tr><td>260</td><td>0.124</td><td>0.033</td><td>1537.59</td><td>0.47</td><td>1.01</td><td>0.13</td><td>0.000025</td><td>0.023</td><td>0.000037</td><td>0.012</td></tr><tr><td>270</td><td>0.230</td><td>0.041</td><td>1515.67</td><td>0.85</td><td>0.82</td><td>0.16</td><td>0.000088</td><td>0.045</td><td>0.000069</td><td>0.013</td></tr><tr><td>280</td><td>0.42</td><td>0.10</td><td>1493.65</td><td>1.43</td><td>0.68</td><td>0.35</td><td>0.00030</td><td>0.115</td><td>0.00012</td><td>0.014</td></tr><tr><td>290</td><td>0.65</td><td>0.10</td><td>1471.78</td><td>2.11</td><td>0.57</td><td>0.36</td><td>0.00082</td><td>0.205</td><td>0.00021</td><td>0.017</td></tr><tr><td>300</td><td>1.21</td><td>0.17</td><td>1448.61</td><td>3.77</td><td>1.07</td><td>0.59</td><td>0.00192</td><td>0.275</td><td>0.00024</td><td>0.018</td></tr><tr><td>310</td><td>1.76</td><td>0.10</td><td>1424.68</td><td>5.19</td><td>1.20</td><td>0.33</td><td>0.00358</td><td>0.3483</td><td>0.00010</td><td>0.0082</td></tr><tr><td>320</td><td>2.82</td><td>0.18</td><td>1400.51</td><td>7.92</td><td>1.06</td><td>0.59</td><td>0.00557</td><td>0.385</td><td>0.00018</td><td>0.013</td></tr><tr><td>330</td><td>4.12</td><td>0.16</td><td>1375.94</td><td>11.15</td><td>0.77</td><td>0.50</td><td>0.00820</td><td>0.4066</td><td>0.00036</td><td>0.0067</td></tr><tr><td>340</td><td>6.06</td><td>0.31</td><td>1348.55</td><td>15.52</td><td>0.99</td><td>0.96</td><td>0.01312</td><td>0.4577</td><td>0.00045</td><td>0.0093</td></tr><tr><td>350</td><td>8.75</td><td>0.33</td><td>1320.12</td><td>21.71</td><td>1.03</td><td>1.00</td><td>0.01960</td><td>0.4875</td><td>0.00070</td><td>0.0074</td></tr><tr><td>360</td><td>12.45</td><td>0.38</td><td>1290.56</td><td>30.23</td><td>1.14</td><td>1.15</td><td>0.02793</td><td>0.5060</td><td>0.00077</td><td>0.0032</td></tr><tr><td>370</td><td>17.81</td><td>0.45</td><td>1257.53</td><td>42.41</td><td>1.70</td><td>1.44</td><td>0.04080</td><td>0.5309</td><td>0.00170</td><td>0.0057</td></tr><tr><td>380</td><td>24.10</td><td>0.70</td><td>1219.54</td><td>56.36</td><td>1.00</td><td>2.27</td><td>0.05746</td><td>0.5573</td><td>0.00085</td><td>0.0076</td></tr><tr><td>390</td><td>32.63</td><td>0.55</td><td>1178.09</td><td>76.73</td><td>1.90</td><td>2.07</td><td>0.0784</td><td>0.5636</td><td>0.0014</td><td>0.0046</td></tr><tr><td>395</td><td>37.95</td><td>1.01</td><td>1152.09</td><td>89.27</td><td>2.20</td><td>3.68</td><td>0.0952</td><td>0.5768</td><td>0.0025</td><td>0.0080</td></tr><tr><td>400</td><td>44.34</td><td>0.24</td><td>1128.36</td><td>104.85</td><td>2.03</td><td>0.87</td><td>0.1096</td><td>0.5834</td><td>0.0014</td><td>0.0011</td></tr><tr><td>405</td><td>50.45</td><td>1.39</td><td>1096.56</td><td>121.37</td><td>3.27</td><td>6.79</td><td>0.1295</td><td>0.586</td><td>0.0038</td><td>0.012</td></tr><tr><td>410</td><td>58.34</td><td>1.26</td><td>1062.98</td><td>141.98</td><td>5.36</td><td>5.79</td><td>0.1511</td><td>0.5861</td><td>0.0051</td><td>0.0085</td></tr><tr><td>415</td><td>67.35</td><td>1.77</td><td>1029.43</td><td>170.26</td><td>3.90</td><td>9.75</td><td>0.1718</td><td>0.580</td><td>0.0030</td><td>0.012</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** P-{x,y} diagram of the reactive mixture containing N₂O₄ and NO₂. Crosses are the points determined with Brick-CFCMC, with the relative uncertainties; circles are points from MC simulations published in Bourasseau et al. [32].</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Vapour and liquid densities of the coexisting phases of the reactive mixture containing N₂O₄ and NO₂, as a function of temperature. Crosses represent results from Brick-CFCMC while the other are experimental data from the literature [40,71–75].</figcaption></figure>

Considering the fact that not a single mixture parameter has been optimised in the used equation of state for this work and that the model works in purely predictive form, the saturation pressure of the reactive mixture is appropriately represented.

Furthermore, although the quantitative representation of the Monte Carlo results on the P-{x,y} diagram in Fig. 14 is not well reproduced by the equation of state, the qualitative assessment of the phase behaviour of this system by the use of such a non-optimised equation of state is highly satisfactory. As shown in Fig. 15, the modelling of vapour and liquid-phase densities is also acceptable, considering that the representation of the liquid phase density by cubic equations of state -without a translation term [77,78]- is known to be inaccurate.

To conclude, the comparison of the results obtained by the two sets of critical coordinates and acentric factors leads to the observation that the impact of their accuracy on the P-ρ-T-{x,y} is non-negligible and that

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Global phase equilibrium diagram of the reactive mixture. The legend of this plot is detailed in the main text.</figcaption></figure>

their further optimisation over both pure component and mixture data is suggested. However, it is possible to observe that the impact on derived properties such as entropy (Fig. 16) and enthalpy (Fig. 17) is lower and that the energy and entropy balances will probably be less affected by the uncertainty on the critical coordinates of the pure components forming the mixture.

To conclude, Fig. 18 shows the unique and non-isothermal phase equilibrium diagram of the reactive NO₂ + N₂O₄ system (the black line – also shown in Fig. 14 with the continuous line) and 5 isothermal phase equilibrium diagrams of the binary mixture of NO₂ + N₂O₄ without reaction, calculated at 5 different temperatures (the coloured lines). This plot also shows isobaric lines corresponding to the saturation pressure of the reactive system at the 5 different temperatures at which the isotherms of the inert mixture have been traced. It can be seen that the phase diagram of the reactive system and each of the phase diagrams of the non-reactive system intersect in two points at equal pressure (the relative saturation pressure of the non-reactive isothermal system). Those two points are, respectively, the bubble and the dew point whose compositions are the solution of reactive system under vapour-liquid equilibrium at the considered temperature (and corresponding calculated pressure). For completeness, it is also specified that the equation of state used to generate all these diagrams is the model described in section 2.3, while the critical coordinates are the ones obtained in this work from results of MC Brick-CFCMC and reported in Table 5.

## 4. Conclusion

This work introduces a methodology enabling the assessment of the VLE and energetic properties of reactive fluids being characterised by one reversible reaction at chemical equilibrium, for example N₂O₄ ⇌ 2NO₂. These systems are known to be highly specific because, the reaction reaching instantaneously the chemical equilibrium composition dictated by the system's temperature and pressure, it is never possible to characterise experimentally the pure N₂O₄ and the pure NO₂ at pre-defined conditions. We have initially recalled in this paper some singular thermodynamic behaviours of these systems (such as the uniqueness of its single P-T curve in VLE condition and thus its unique critical point) and related that to the equal variance of such a binary and mono-reaction system with respect to a mono-component (inert) fluid.

To calculate thermodynamic properties, the Peng-Robinson equation of state associated to the predictive advanced mixing rules EoS + a<sub>res</sub><sup>Eγ</sup>, with a zero-a<sub>res</sub><sup>Eγ</sup> term characterising athermal solutions, is applied. The
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** P-{x,y} diagram of the reactive mixture containing N₂O₄ and NO₂. Fig. (a) represents the entire figure, Fig. (b) represents the ZOOM of Fig. (a). In the figures, crosses are the points determined with Brick-CFCMC, with the relative uncertainties; circles represent points published by Bourasseau et al. [32]; lines represent calculations with the equation of state associated to the two different sets of critical coordinates, solid lines represent EoS 1, while the dashed line represent EoS 2.</figcaption></figure>

required input critical temperature, pressure and acentric factors of the pure N₂O₄ and NO₂, being their experimental determination impossible, have been optimised in this work from force field-based Monte Carlo results of Gibbs and NPT ensemble simulations, using the force field-based simulation codes Brick-CFCMC and GIBBS. The comparison of the two sets of simulation results obtained with the two codes has proved their high compatibility. Although GIBBS allows performing simulations closer to the critical point than CFC methods, the impact on the prediction of the critical point remains limited.

Moreover, using reaction ensembles, Brick-CFCMC has also enabled the determination of reactive mixture properties, that we have then used to assess the capability of Monte Carlo calculations to reproduce experimental data. Considering the reported uncertainty of MC simulations, this comparison has resulted in a highly satisfactory reproduction of experimental data points. The calculation of energetic properties for this reactive system (enthalpies and entropies) requires the

<figure><img src="image_17.png" /><figcaption>**Fig. 15.** T-ρ<sup>vap</sup>-ρ<sup>liq</sup> diagram of the reactive mixture containing N₂O₄ and NO₂. Crosses are the points determined with Brick-CFCMC, with the relative uncertainties; lines represent calculations with the equation of state associated to the two different sets of critical coordinates. The other data points are the same already reported in Fig. 12.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 16.** T-s diagram of the reactive mixture containing N₂O₄ and NO₂, calculated with the thermodynamic model defined in this work.</figcaption></figure>

calculation of the thermochemical properties of the molecules forming the reactive mixture (ideal gas standard enthalpy of formation, molar standard entropy and heat capacity). Those properties have been determined from ab-initio quantum mechanics simulations.

Furthermore, thermodynamic calculations have been finally performed on the reactive system, with the selected predictive equation of state, and its input T_c, P_c and ω of NO₂ and of N₂O₄ optimised on MC results of the two respective systems, with the required ΔfH<sub>i,298.15K</sub>, S<sub>i,298.15K</sub> and the ideal-gas heat capacity temperature-dependent correlation obtained from Quantum Mechanics. The comparison with available P-ρ-T experimental data of the reactive system in VLE conditions has shown a very good agreement. However, a high sensibility of the accuracy of the equation of state to the critical coordinates used in this work, either obtained from MC simulations performed with Brick-CFCMC or with GIBBS, has been observed. That paves the way to a future investigation, intended to explore the possibility to optimise
<figure><img src="image_19.png" /><figcaption>**Fig. 17.** T-h diagram of the reactive mixture containing N₂O₄ and NO₂, calculated with the thermodynamic model defined in this work.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 18.** Phase equilibrium diagrams (solid lines) of: the **reactive** NO₂-N₂O₄ mixture (unique non-isothermal phase equilibrium diagram) - black line; the **inert** NO₂-N₂O₄ mixture - coloured lines at 270 K (red), 310 K (blue), 350 K (magenta), 390 K (green), 430 K (light blue). Isobars (dashed lines) corresponding to the VLE pressures calculated for the reactive system at the 5 temperatures above considered to build the phase equilibrium diagrams of the inert mixture.</figcaption></figure>

critical coordinates of NO₂ and N₂O₄ on MC of both pure components and reactive mixture data.

To conclude, the presented methodology has proven to be able to accurately determine the chemical equilibrium properties of N₂O₄ ⇌ 2NO₂, and is thus proposed to preliminarily predict the thermodynamic properties of other similar reactive systems.

### CRediT authorship contribution statement

**Silvia Lasala:** Writing – original draft, Visualization, Validation, Supervision, Software, Project administration, Methodology, Funding acquisition, Formal analysis, Conceptualization, Investigation, Data curation.
**Konstantin Samukov:** Visualization, Software.
**H. Mert Polat:** Software, Supervision.
**Véronique Lachet:** Writing – review & editing, Software, Investigation, Data curation.
**Olivier Herbinet:**

Writing – review & editing, Methodology, Investigation.
**Romain Privat:** Supervision, Software, Conceptualization, Formal analysis.
**Jean-Noël Jaubert:** Writing – review & editing, Supervision, Conceptualization, Formal analysis.
**Othonas A. Moultos:** Software.
**Kevin De Ras:** Writing – review & editing.
**Thijs J. H. Vlugt:** Writing – review & editing, Writing – original draft, Supervision, Software.

### Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Silvia Lasala reports financial support was provided by European Commission. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.

### Acknowledgements

This work has received funding from the European Research Council (ERC) under the European Union's Horizon Europe research and innovation program (grant agreement No. 101040994). We are grateful for the support by NWO Domain Science for the use of supercomputer facilities, with financial support from the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (Netherlands Organisation for Scientific Research, NWO). The authors acknowledge the use of computational resources of DelftBlue supercomputer, provided by Delft High Performance Computing Centre (https://www.tudelft.nl/dhpc).

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.cej.2024.148961.

### References

[1] E.T. Chang, N.A. Gokcen, Thermodynamic properties of gases in propellants and oxidizers. I. Solubilities of He, N₂, O₂, Ar, and N₂O₃ in liquid N₂O₄, J. Phys. Chem. 70 (1966) 2394-2399, https://doi.org/10.1021/j100879a050.

[2] G. Angelino, Performance of N₂O₄ gas cycles for solar power applications, Proceedings of the Institution of Mechanical Engineers 1847-1982 (Vols 1-196) 193 (1979) 313-320. doi: 10.1243/PIME_PROC.1979_193_033_02.

[3] H.-M. Huang, R. Govind, Optimisation of power plants using dissociating gases as working fluids, in: Proceedings of the 20th Intersociety Energy Conversion Engineering Conference, Warrendale, PA 15096, 1985: pp. 620-629.

[4] H.M. Huang, R. Govind, Use of dissociating gases in Brayton Cycle space power systems, Ind. Eng. Chem. Res. 27 (1988) 803-810, https://doi.org/10.1021/ie00077a015.

[5] M. Blander, L.G. Epel, A.P. Fraas, R.F. Newton, Aluminum chloride as a thermodynamic working fluid and heat transfer medium, Oak Ridge National Laboratory, Union Carbide Corporation, U.S. Atomic Energy Commission, 1959.

[6] A. Sorokin, Dissociating Nitrogen Dioxide (N₂O₄) as a working fluid in thermodynamic cycles, Nucl. Sci. Eng. 72 (1979) 330-346.

[7] K. Kesavan, The use of dissociating gases as the working fluid in thermodynamic power conversion cycles, Carnegie-Mellon University, Nuclear Science and Engineering Division, 1978.

[8] T.A. Jacobs, J.R. Lloyd, The influence of the equilibrium dissociation of a diatomic gas on brayton-cycle performance, J. Appl. Mech. (1963) 288-290.

[9] V.B. Nesterenko, V.P. Bubnov, E.N. Bunin, N.M. Shiryaeva, Thermodynamic cycles on chemically reacting fluids, Vesci Akad. Navuk BSSR, Ser. Fiz.-Energ. Navuk (1970) 74-76.

[10] A.K. Krasin, V.B. Nesterenko, Dissociating gases: a new class of coolants and working substances for large power plants, At. Energy Rev. 9 (1971) 177.

[11] R. Stochl, Potential performance improvement using a reacting gas (nitrogen tetroxide) as the working fluid in a closed Brayton cycle, DOE/NASA/1060-79/3 (1979), https://doi.org/10.2172/5746805.

[12] K. Kesavan, J.F. Osterle, Brayton Cycle Using Dissociating Nitrosyl Chloride (1981) 2204-2209.

[13] M.A. Bazhin, V.P. Bubnov, V.B. Nesterenko, N.M. Shiryaeva, Optimisation of parameters of power plants using dissociating working fluids, Nauka i Tehnika, Minsk, 1970.
[14] S. Lasala, R. Privat, O. Herbinet, P. Arpentinier, D. Bonalumi, J.-N. Jaubert, Thermo-chemical systems: unexploited high-potential energy converters, Energy Convers. Manag. 229 (2021) 113685, https://doi.org/10.1016/j.enconman.2021.113685.

[15] A. Barakat, S. Lasala, P. Arpentinier, J.-N. Jaubert, The original and impactful exploitation of chemical energy in heat pumps, Chemical Engineering Journal Advances 12 (2022) 100400, https://doi.org/10.1016/j.ceja.2022.100400.

[16] Dissociating gases as coolants and working substances in power plants (Proceedings of All-Union Conference), Nauka i Tehnika, Minsk, 1970.

[17] K. Zhang, Y. Shen, C. Duwig, Identification of heat transfer intensification mechanism by reversible N2O4 decomposition using direct numerical simulation, Int. J. Heat Mass Transf. 182 (2022) 121946, https://doi.org/10.1016/j.ijheatmasstransfer.2021.121946.

[18] K. Zhang, A. Laitinen, Y. Shen, V. Vuorinen, C. Duwig, Reactive cooling simulation of electronic components, Appl. Therm. Eng. 228 (2023) 120519, https://doi.org/10.1016/j.applthermaleng.2023.120519.

[19] S. Lasala, Reactive fluids for intensified thermal energy conversion, Project Repository Journal 13 (2022) 102–105, https://doi.org/10.54050/PRJ1318808.

[20] S. Lasala, ERC-REACHER, (2022). https://www.univ-lorraine.fr/erc-ercacher/.

[21] V. Diky, C. Muzny, A. Smolanitsky, A. Bazileva, R. Chirico, J. Magee, Y. Paulechka, A. Kazakov, S. Townsend, E. Lemmon, M. Frenkel, K. Kroenlein, ThermoData Engine (TDE) Version 10.1 (Pure Compounds, Binary Mixtures, Ternary Mixtures, and Chemical Reactions): NIST Standard Reference Database 103b, (2016).

[22] Dortmund Data Bank, (2023). http://www.ddbst.com.

[23] D.F. Stai, F. Bizjak, S.E. Stephanou, Thermodynamic properties of nitrogen tetroxide, J. Spacecrr. Rocket. 2 (1965) 742–745, https://doi.org/10.2514/3.28272.

[24] D.N. Seshadri, D.S. Viswanath, N.R. Kuloor, Thermodynamic properties of the system N2O4=2NO2=2NO + O2, AIChE J. 16 (1970) 420–425, https://doi.org/10.1002/aic.690160319.

[25] O.A. Hougen, K.M. Watson, Chemical Process Principles, John Wiley & Sons Inc, New York, 1952.

[26] R.D. McCarty, H.-U. Steurer, C.M. Daily, The thermodynamic properties of nitrogen tetroxide, Thermophysics Division, Center for Chemical Engineering, National Engineering Laboratory, National Bureau of Standards, Boulder, 1986.

[27] M. Binotti, C.M. Invernizzi, P. Iora, G. Manzolini, Dinitrogen tetroxide and carbon dioxide mixtures as working fluids in solar tower plants, Sol. Energy 181 (2019) 203–213, https://doi.org/10.1016/j.solener.2019.01.079.

[28] G. Manzolini, M. Binotti, D. Bonalumi, C. Invernizzi, P. Iora, CO₂ mixtures as innovative working fluid in power cycles applied to solar plants, techno-economic assessment, Solar Energy 181 (2019) 530–544, https://doi.org/10.1016/j.solener.2019.01.015.

[29] L.E.S. de Souza, U.K. Deiters, Modeling of the N₂O₄–NO₂ reacting system, Phys. Chem. Chem. Phys. 2 (2000) 5606–5613, https://doi.org/10.1039/b005464j.

[30] A. Belkadi, F. Llovell, V. Gerbaud, L.F. Vega, Modeling the vapor–liquid equilibrium and association of nitrogen dioxide/dinitrogen tetroxide and its mixtures with carbon dioxide, Fluid Phase Equilib. 266 (2008) 154–163, https://doi.org/10.1016/j.fluid.2008.01.026.

[31] E.W. Lemmon, A.H. Harvey, J.M. Young, Thermodynamic Model for the Reactive Mixture Comprising Dinitrogen Tetroxide (N2O4), Nitrogen Dioxide (NO2), Nitric Oxide (NO), and Oxygen (O2), (2022).

[32] E. Bourrasseau, V. Lachet, N. Desbiens, J.-B. Maillet, J.-M. Teuler, P. Ungerer, Thermodynamic behavior of the CO₂ + NO₂/N2O4 mixture: a Monte Carlo simulation study, J. Phys. Chem. B 112 (2008) 15783–15792, https://doi.org/10.1021/jp8068255.

[33] J.K. Johnson, A.Z. Panagiotopoulos, K.E. Gubbins, Reactive canonical Monte Carlo, Mol. Phys. 81 (1994) 717–733, https://doi.org/10.1080/00268979400100481.

[34] W.R. Smith, B. Triska, The reaction ensemble method for the computer simulation of chemical and phase equilibria. I. Theory and basic examples, J. Chem. Phys. 100 (1994) 3019–3027, https://doi.org/10.1063/1.466443.

[35] J.K. Shah, E. Marin-Rimoldi, R.G. Mullen, B.P. Keene, S. Khan, A.S. Paluch, N. Rai, L.L. Romanielo, T.W. Rosch, B. Yoo, E.J. Maginn, Cassandra: An open source Monte Carlo package for molecular simulation, J. Comput. Chem. 38 (2017) 1727–1739, https://doi.org/10.1002/jcc.24807.

[36] M.J. Molina, S.B. Rodriguez-Reartes, M.S. Zabaloy, A theoretical study on the simultaneous vapor–liquid and chemical equilibria in a highly restricted system, Fluid Phase Equilib. 557 (2022) 113439, https://doi.org/10.1016/j.fluid.2022.113439.

[37] M.J. Molina, S.B. Rodriguez-Reartes, M.S. Zabaloy, Computation and analysis of reactive isopleths involving fluid phases, Fluid Phase Equilib. 574 (2023) 113889, https://doi.org/10.1016/j.fluid.2023.113889.

[38] J.M. Smith, H.C.V. Ness, M.M. Abbott, Introduction to chemical engineering thermodynamics, 7 edizione, McGraw-Hill Education, Boston, 2004.

[39] F.E.C. Scheffer, J.P. Treub, Determinations of the vapour tension of nitrogen tetroxide, KNAW, Proceedings (1912) 166–178.

[40] H.H. Reamer, B.H. Sage, Volumetric behavior of nitrogen dioxide in the liquid phase, Ind. Eng. Chem. 44 (1952) 185–187, https://doi.org/10.1021/ie5050a502.

[41] DECHEMA, DIPPR 801 database, (2023). https://dechema.de/en/dippr801.html.

[42] D.-Y. Peng, D.B. Robinson, The characterization of the heptanes and heavier fractions for the GPA Peng-Robinson programs, Gas Processors Association, 1978.

[43] S. Lasala, P. Chiesa, R. Privat, J.-N. Jaubert, VLE properties of CO₂-based binary systems containing N₂, O₂ and Ar: experimental measurements and modelling results with advanced cubic equations of state, Fluid Phase Equilib. 428 (2016) 18–31, https://doi.org/10.1016/j.fluid.2016.05.015.

[44] J.W. Qian, R. Privat, J.N. Jaubert, P. Duchet-Suchaux, Enthalpy and heat capacity changes on mixing: fundamental aspects and prediction by means of the PPR78 cubic equation of state, Energy Fuel 27 (2013) 7150–7178, https://doi.org/10.1021/ef401605c.

[45] R. Hens, A. Rahbari, S. Caro-Ortiz, N. Dawass, M. Erdös, A. Poursaeidesfahani, H. S. Salehi, A.T. Celebi, M. Ramdin, O.A. Moultos, D. Dubbeldam, T.J.H. Vlugt, Brick-CFCMC: open source software for monte carlo simulations of phase and reaction equilibria using the continuous fractional component method, J. Chem. Inf. Model. 60 (2020) 2678–2682, https://doi.org/10.1021/acs.jcim.0c00334.

[46] P. Ungerer, B. Tavitan, A. Boutin, Applications of Molecular Simulation in the Oil and Gas Industry - Monte-carlo Methods, Editions TECHNIP, 2005.

[47] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, G. Scalmani, V. Barone, G.A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. Marenich, J. Bloino, B.G. Janesko, R. Gomperts, B. Mennucci, H.P. Hratchian, J. V. Ortiz, A.F. Izmaylov, J.L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throllsell, J.A. Montgomery Jr., J.E. Peralta, F. Ogliaro, M. Bearpark, J.J. Heyd, E. Brothers, K.N. Kudin, V.N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J.C. Burant, S.S. Iyengar, J. Tomasi, M. Cossi, J.M. Millam, M. Klene, C. Adamo, R. Cammi, J.W. Ochterski, R.L. Martin, K. Morokuma, O. Farkas, J.B. Foresman, D.J. Fox, Gaussian 09 (2016). https://gaussian.com/.

[48] A. Miyoshi, GPOP software, (2002). http://akrmys.com/gpop/.

[49] J.A. Montgomery, M.J. Frisch, J.W. Ochterski, G.A. Petersson, A complete basis set model chemistry. VI. Use of density functional geometries and frequencies, J. Chem. Phys. 110 (1999) 2822–2827, https://doi.org/10.1063/1.477924.

[50] A.D. Becke, Density-functional thermochimology. III. The role of exact exchange, J. Chem. Phys. 98 (1993) 5648–5652, https://doi.org/10.1063/1.464913.

[51] K.S. Pitzer, W.D. Gwinn, Energy levels and thermodynamic functions for molecules with internal rotation I. Rigid frame with attached tops, The Journal of Chemical Physics 10 (2004) 428–440, https://doi.org/10.1063/1.1723744.

[52] D. Frenckel, B. Smit, Understanding Molecular Simulation. From Algorithms to Applications., 3rd Edition, 2002.

[53] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, 2nd edition, Oxford University Press, 2017.

[54] W. Shi, E.J. Maginn, Continuous fractional component Monte Carlo: an adaptive biasing method for open system atomistic simulations, J Chem Theory Comput 3 (2007) 1451–1463, https://doi.org/10.1021/ct7000039.

[55] W. Shi, E.J. Maginn, Improvement in molecule exchange efficiency in Gibbs ensemble Monte Carlo: development and implementation of the continuous fractional component move, J. Comput. Chem. 29 (2008) 2520–2530, https://doi.org/10.1002/jcc.20977.

[56] A. Rahbari, R. Hens, M. Ramdin, O.A. Moultos, D. Dubbeldam, T.J.H. Vlugt, Recent advances in the continuous fractional component Monte Carlo methodology, Mol. Soln. 47 (2021) 804–823, https://doi.org/10.1080/08927022020182858.

[57] H.M. Polat, H.S. Salehi, R. Hens, D.O. Wasik, A. Rahbari, F. De Meyer, C. Houriez, C. Coquelet, S. Calero, D. Dubbeldam, O.A. Moultos, T.J.H. Vlugt, New features of the open source monte carlo software brick-CFCMC: thermodynamic integration and hybrid trial moves, J. Chem. Inf. Model. 61 (2021) 3752–3757, https://doi.org/10.1021/acs.jcim.1c00652.

[58] A. Poursaeidesfahani, A. Torres-Knoop, D. Dubbeldam, T.J.H. Vlugt, Direct free energy calculation in the continuous fractional component gibbs ensemble, J. Chem. Theory Comput. 12 (2016) 1481–1490, https://doi.org/10.1021/acs.jctc.5b01230.

[59] F. Wang, D.P. Landau, Efficient, multiple-range random walk algorithm to calculate the density of states, Phys. Rev. Lett. 86 (2001) 2050–2053, https://doi.org/10.1103/PhysRevLett.86.2050.

[60] A. Poursaeidesfahani, R. Hens, A. Rahbari, M. Ramdin, D. Dubbeldam, T.J. H. Vlugt, Efficient application of continuous fractional component Monte Carlo in the reaction ensemble, J. Chem. Theory Comput. 13 (2017) 4452–4466, https://doi.org/10.1021/acs.jctc.7b00092.

[61] H.M. Polat, F. de Meyer, C. Houriez, O.A. Moultos, T.J.H. Vlugt, Solving chemical absorption equilibria using free energy and quantum chemistry calculations: methodology, limitations, and new open-source software, J. Chem. Theory Comput. 19 (2023) 2616–2629, https://doi.org/10.1021/acs.jctc.3c00144.

[62] A.Z. Panagiotopoulos, Direct determination of phase coexistence properties of fluids by Monte Carlo simulation in a new ensemble, Mol. Phys. 61 (1987) 813–826, https://doi.org/10.1080/00268978700101491.

[63] V. Lachet, T. de Bruin, P. Ungerer, C. Coquelet, A. Valtz, V. Hasanov, F. Lockwood, D. Richon, Thermodynamic behavior of the CO₂+SO₂ mixture: experimental and Monte Carlo simulation studies, Energy Procedia 1 (2009) 1641–1647, https://doi.org/10.1016/j.egypro.2009.01.215.

[64] V. Lachet, B. Creton, T. de Bruin, E. Bourrasseau, N. Desbiens, Ø. Wilhelmsen, M. Hammer, Equilibrium and transport properties of CO₂+N2O and CO₂+NO mixtures: molecular simulation and equation of state modelling study, Fluid Phase Equilib. 322–323 (2012) 66–78, https://doi.org/10.1016/j.fluid.2012.03.011.

[65] B. Creton, C. Nieto-Draghi, T. de Bruin, V. Lachet, E. El Ahmar, A. Valtz, C. Coquelet, S. Lasala, R. Privat, J.-N. Jaubert, Thermodynamic study of binary systems containing sulphur dioxide and nitric oxide: Measurements and modelling, Fluid Phase Equilib. 461 (2018), https://doi.org/10.1016/j.fluid.2017.12.036.

[66] J.S. Rowlinson, F.L. Swinton, Liquids and liquid mixtures, 3rd edition, Butterworths, 1982.

[67] C. Heath Turner, J.K. Brennan, M. López, W.R. Smith, J. Karl Johnson, K.E. Gubbins, Simulation of chemical reaction equilibria by the reaction ensemble Monte Carlo
method: a review, Mol. Simul. 34 (2008) 119–146, https://doi.org/10.1080/08927020801986564.

[68] S.C. Glotzer, D. Stauffer, N. Jan, Monte Carlo simulations of phase separation in chemically reactive binary mixtures, Phys. Rev. Lett. 72 (1994) 4109–4112, https://doi.org/10.1103/PhysRevLett.72.4109.

[69] Y. Mo, W. Wu, Q. Zhang, Valence bond studies of N2O4, J. Mol. Struct. (Thoechem) 315 (1994) 173–178, https://doi.org/10.1016/0166-1280(94)03778-J.

[70] P.J. Linstrom, W.G. Mallard, NIST Chemistry WebBook (2023), https://doi.org/10.18434/T4D303.

[71] A. Mittasch, E. Kuss, H. Schlueter, Dichten und Dampfdrucke von wäßrigen Ammoniaklösungen und von flüssigem Stickstofftetroxyd für das Temperaturgebiet 0° bis 60°, Z. Anorg. Allg. Chem. 159 (1927) 1–36, https://doi.org/10.1002/zaac.19261590102.

[72] N.G. Polikhronidi, R.G. Batyrova, I.M. Abdulagatov, Isochoric heat capacity measurements of nitrogen tetroxide system at temperatures between 410 and 484 K and pressures up to 35 MPa, Fluid Phase Equilib. 175 (2000) 153–174, https://doi.org/10.1016/S0378-3812(00)00457-X.

[73] A.B. Verzhinskaya, V.I. Tsurbelev, P.M. Klepatskiy, Impact of nitrogen monooxide on shape of the coexistence curve of N₂O₄ and vicinity of the critical point liquid-vapor, in Application of mathematical methods for description and studying of the phisico-chemical equilibria (extended abstracts of reports), in: Novosibirsk, 1985.

[74] P. Gray, P. Rathbone, Dissociation of liquid dinitrogen tetroxide; Henry's law coefficients, heats and entropies of solution, and the thermodynamics of homolytic dissociation in the pure liquid, J. Chem. Soc. (1958) 3550, https://doi.org/10.1039/jr9580003550.

[75] C.C. Addison, B.C. Smith, Volume changes on mixing organic liquids with dinitrogen tetroxide: comparison with sulphur dioxide systems, J. Chem. Soc. (1958) 3664, https://doi.org/10.1039/jr9580003664.

[76] R. Privat, J.-N. Jaubert, Classification of global fluid-phase equilibrium behaviors in binary systems, Chem. Eng. Res. Des. 91 (2013) 1807–1839, https://doi.org/10.1016/j.cherd.2013.06.026.

[77] 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, https://doi.org/10.1016/j.fluid.2016.03.012.

[78] 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, https://doi.org/10.1016/j.fluid.2016.07.035.