Phase equilibria of symmetric Lennard-Jones mixtures and a look at the transport properties near the upper critical solution temperature

This study investigates phase equilibria and transport properties of five symmetric binary Lennard-Jones mixtures using molecular simulation and equation of state models. The mixtures are selected for their representation of different types of phase behavior and the research contributes to the development of simulation techniques, mixture theories and understanding of thermophysical mixture properties. A novel method is introduced for determining the critical end point (CEP) and critical azeotropic end point (CAEP) by molecular simulation. The van der Waals one-fluid theory is assessed for its performance in conjunction with Lennard-Jones equation of state models, while addressing different phase equilibrium types simultaneously. An empirical correlation is introduced to account for deviations between the equation of state and simulation that arise when using the same binary interaction parameter. This study also investigates the influence of the liquid-liquid critical point on thermophysical properties, which are found to exhibit no significant anomalies or singularities. System-size effects of diffusion coefficients are addressed by extrapolating simulation data to the thermodynamic limit and applying analytical finite-size corrections.


LLE and vapor phase simulations
LLE and vapor phase simulations were carried out using Monte Carlo (MC) simulations in the isobaric-isothermal (N pT ) ensemble containing 1000 to 4000 molecules, with an average of 8 · 10 5 production cycles. The chemical potential was sampled with Widom's insertion method, S6 using 3000 test particles per cycle.

VLE simulations
VLE simulations were conducted in two subsequent steps: (1) liquid run and (2) vapor run, following the grand equilibrium method by Vrabec and Hasse S7 . The liquid run was carried out with MD simulation in the N pT ensemble, followed by the vapor run using MC simulation in the pseudo-grand canonical (µV T ) ensemble. Both phases consisted of 1372 molecules and were sampled for 8 · 10 5 production steps.

Transport property simulations
Transport property simulations were performed in two consecutive steps: (1) MC simulation in the N pT ensemble to obtain the density and (2) MD simulation in the canonic (N V T ) ensemble to obtain the diffusion coefficients, shear viscosity and thermal conductivity, as well as the thermodynamic factor through Kirkwood-Buff integrals (KBI). S8 Simulations in the N V T ensemble were first equilibrated over 1.5 · 10 5 time steps, followed by a production run of 5 · 10 6 time steps. To account for finite size effects, up to four system sizes were investigated containing 2500, 5000, 10,000 and 20,000 molecules, respectively.

Helmholtz energy derivatives
The molar Helmholtz energy a = F/N can be used to determine all thermodynamic equilibrium properties, such as the isochoric heat capacity c v , speed of sound w or the Joule-Thomson coefficient µ JT . It also serves as the basis for the development of EOS. S9 The formalism for the calculation of thermodynamic properties from Helmholtz energy-based EOS is described in detail in Refs. S10, S11. The partial derivatives of the reduced Helmholtz energy α = a/(k B T ) can be expressed as where m and n denote the order of derivation with respect to τ = 1/T and ρ. The reduced Helmholtz energy derivatives are usually split into an ideal gas A o mn and a residual A r mn contribution. Direct sampling of the latter is implemented in ms2 and allows to extend this study to further properties near the critical solution point.
Reduced residual Helmholtz energy : A r 00 =

Thermodynamic factor
The thermodynamic factor is a key component for the description of both phase equilibria and mass transport properties. It can be used to identify the composition where phase stability ceases to exist (Γ = 0) and allows to determine the critical solution point, as demonstrated in this work. Moreover, it is also relevant for diffusive mass transport processes, since it S-3 T where g ij (r) is the radial distribution function (RDF) integrated over a spherical volume.
Since the infinite integration limit in Eq. (4) cannot be reached with molecular simulations because they are limited to finite radii R, a geometric weight function w(r, x) was proposed by Krüger et al. S12 where x = r/2R. Like the self-diffusion coefficient, the integrals G ij scale linearly with the inverse of the system size 1/R. Thus, the extrapolation to the thermodynamic limit (1/R → 0) is straightforward. Three types of RDF were used to calculate the KBI: the S-4 standard RDF, and two modifications -RDF corrected according to Ganguly and van der Vegt S13 (vdV) and a shifted version thereof (vdV + shf). S14 This amounts to three values for Γ per simulation run, plus three additional values for their extrapolated versions. The thermodynamic factor of a binary mixture was then calculated by 3 Transport properties Transport properties were obtained with equilibrium molecular dynamics (MD) simulations, mainly with the Einstein formalism. Three types of diffusion coefficients were considered: (1) self-diffusion D i , (2) Maxwell-Stefan (MS) Ð ij and (3) Fick diffusion D ij . Shear viscosity η, thermal conductivity λ and thermodynamic factor Γ were investigated, as the latter is an important auxiliary quantity to convert between Ð ij and D ij .

Self-diffusion coefficient
Self-diffusion is based on the Brownian motion of individual molecules in absence of driving forces. The self-diffusion coefficient of species i D i is given by the Einstein relation with the displacement where N i is the number of molecules of species i, r k,i (t) is the position vector of molecule k of species i and the angle brackets . . . denote the ensemble average. S-5

Maxwell-Stefan and Fick diffusion
MS and Fick diffusion both describe the collective motion of molecular systems. However, they differ in the underlying driving force that is assumed for the diffusive flux. Fick diffusion assumes the concentration gradient to be the cause of mass transport, whereas MS assumes the chemical potential gradient. From a physical standpoint, the MS approach is preferable.
Krishna S15 demonstrated its advantages, being unambiguous and more general. For a binary mixture, the MS diffusion coefficient can be obtained by on the basis of the Onsager coefficients where N i and N j are the number of molecules of component i and j, respectively, and N is the total number of molecules. The Onsager coefficients are given by the cross-correlation of the displacement (9) between species i and j.
Despite the advantages of the MS diffusion framework, Fick's law is still the most common approach for describing diffusive mass transport because of its simple structure and the experimental accessibility of the according coefficient. The conversion between the MS and Fick diffusion coefficients can be achieved through the thermodynamic factor Γ. For binary S-6

Shear viscosity and thermal conductivity
The shear viscosity η was obtained through the Einstein relation where V is the molar volume, α, β = x, y, z are Cartesian coordinates and J αβ p is an off- The thermal conductivity λ was obtained by means of the Green-Kubo formalism where J x q are the elements of the microscopic heat flow. S2

Equation for ξ EOS
The linear function modeling the temperature dependence of the binary interaction parameter needed for the EOS of Kolafa and Nezbeda S16 is given as Its parameters c 1 and c 2 were fitted for each mixture and temperature range individually and are reported in Table S1. S-7

Deviations
The temperature dependence is qualitatively the same for all three mixtures of type II: for temperatures below the T CEP , ξ EOS increases linearly with rising temperature. For temperatures above T CEP , ξ EOS shows only a minor temperature dependence. Figure S2 shows the relative deviation between ξ EOS and ξ as a function of the temperature reduced by T CEP for all three studied type II mixtures (B,C,D). Overall, mixtures of type II show an increasing deviation between ξ EOS and ξ. With decreasing ξ, the relative deviation between ξ EOS and ξ increases. Hence, the predictive capabilities of the EOS decline.  Figure S2: Relative deviation between the binary interaction parameter ξ EOS of the EOS by Kolafa and Nezbeda S16 and its true value ξ as a function of temperature T reduced by the temperature of the critical endpoint (CEP) T CEP . Results were obtained by fitting ξ EOS to VLE and LLE simulation data for mixtures B (purple), C (dark red) and D (orange). S-8

Simulation results
The simulations results for the studied mixtures are listed below. The nomenclature of the mixtures corresponds to 0.60 0.75 0.80 0.85 1.20 All thermodynamic properties were reduced to dimensionless numbers on the order of unity. For brevity, the asterisk is omitted in the following with the understanding that reduced quantities are reported throughout.         Table S3: Simulation results for the vapor-liquid-liquid equilibrium of mixtures A to D. Temperature T , pressure p and mole fraction x 1 are tabulated, where V, L1 and L2 denote the vapor and the two liquid phases.  Table S4: Simulation results for the azeotropic line of mixtures B to E. Temperature T , pressure p, density ρ and mole fraction x 1 are tabulated, where V and L denote the vapor and liquid phase, respectively.  Table S4: Simulation results for the azeotropic line of mixtures B to E. Temperature T , pressure p, density ρ and mole fraction x 1 are tabulated, where V and L denote the vapor and liquid phase, respectively.  Table S5: Simulation results for the transport properties near the critical solution point of the equimolar mixture C at p =0.0724. Temperature T , density ρ, MS Ð ∞ ij and self-diffusion coefficient D i as well as the shear viscosity η and thermal conductivity λ are tabulated. Results correspond to four system sizes N = 2500, 5000, 10,000 and 20,000.  Table S5: Simulation results for the transport properties near the critical solution point of the equimolar mixture C at p =0.0724. Temperature T , density ρ, MS Ð ∞ ij and self-diffusion coefficient D i as well as the shear viscosity η and thermal conductivity λ are tabulated. Results correspond to four system sizes N = 2500, 5000, 10,000 and 20,000.  Table S6: Simulation results for the UCST of mixtures B to D. Temperature T , pressure p and mole fraction x 1 are tabulated, where L1 and L2 denote the two liquid phases.