Primitive models of room temperature ionic liquids. Liquid-gas phase coexistence

TL;DR

This paper develops primitive models and a mean spherical approximation theory to describe the phase coexistence of room temperature ionic liquids, showing good agreement with simulations and qualitative reproduction of experimental phase behavior.

Contribution

It introduces new primitive models of RTILs and an AMSA-based theoretical framework to predict their phase diagrams, aligning well with simulation data.

Findings

The models qualitatively reproduce RTIL phase behavior.

Theoretical phase diagrams agree with computer simulations.

Increasing asymmetry lowers critical temperature and density.

Abstract

We propose several versions of primitive models of room temperature ionic liquids (RTILs) and develop a mean spherical approximation (MSA)-type theory for their description. RTIL is modeled as a two-component mixture of hard-sphere anions and flexible linear chain cations, represented by tangentially bonded hard spheres with the charge located on one of the terminal beads. The theoretical description of the model is carried out using the solution of the appropriately modified associative MSA (AMSA). Our solution reduces to solving one nonlinear algebraic equation for the Blum's screening parameter $\Gamma$, which in turn is used to express all thermodynamic properties of the models of interest. We calculate liquid-gas phase diagrams using theoretical and computer simulation methods for two versions of the model, represented by the dimer ($D$) and chain ($C$) models. Theoretical predictions for the phase diagrams appear to be in reasonably good agreement with computer simulation results. It is demonstrated that the models and theory are able to qualitatively reproduce experimentally observed phase behavior of RTILs, in particular the decrease of the critical temperature and critical density with increasing asymmetry of the model in its shape and position of the charge.

Figures

Equations

$$U_{ij}^{(hs)}(r)=\left\{\begin{array}[]{rl}\infty,&r<\sigma_{ij}=(\sigma_{i}+\sigma_{j})/2,\\ 0,&r>\sigma_{ij},\end{array}\right.$$

$$U_{ij}^{(C)}(r)=({\delta }_{i1}{\delta }_{j2}+{\delta }_{i2}{\delta }_{j1})\frac{e^2{z}_i{z}_j}{4\pi ϵ{ϵ}_{0}r},$$

$$\hat{\bf h}_{ij}(k)=\hat{\bf c}_{ij}(k)+\rho\sum_{l}\hat{\bf c}_{il}(k){\mbox{\boldmath${\alpha}$}}\hat{\bf h}_{lj}(k),$$

$${\hat{\bf h}}_{ij}(k))=\pmatrix{{\hat{h}}^{00}_{ij}(k)&{\hat{h}}^{0A}_{ij}(k)&{\hat{h}}^{0B}_{ij}(k)\cr{\hat{h}}^{A0}_{ij}(k)&{\hat{h}}^{AA}_{ij}(k)&{\hat{h}}^{AB}_{ij}(k)\cr{\hat{h}}^{B0}_{ij}(k)&{\hat{h}}^{BA}_{ij}(k)&{\hat{h}}^{BB}_{ij}(k)\cr},{\hat{\bf c}}_{ij}(k)=\pmatrix{{\hat{c}}^{00}_{ij}(k)&{\hat{c}}^{0A}_{ij}(k)&{\hat{c}}^{0B}_{ij}(k)\cr{\hat{c}}^{A0}_{ij}(k)&{\hat{c}}^{AA}_{ij}(k)&{\hat{c}}^{AB}_{ij}(k)\cr{\hat{c}}^{B0}_{ij}(k)&{\hat{c}}^{BA}_{ij}(k)&{\hat{c}}^{BB}_{ij}(k)\cr},$$

$$\hskip 11.38109pt\mbox{\boldmath${\alpha}$}=\pmatrix{1&1&1\cr 1&0&1\cr 1&1&0\cr},$$

$${\mathbf{c}}_{ij}(r)=-\mathbf{E}\beta U_{ij}^{(C)}(r)+\frac{{\mathbf{t}}_{ij}}{2\pi {\sigma }_{ij}}\delta (r-{\sigma }_{ij}),r\le {\sigma }_{ij}$$

$${\mathbf{h}}_{ij}(r)=-\mathbf{E},r<{\sigma }_{ij}$$

$$t_{ij}^{\alpha \beta }=\frac{1}{2\rho }\left[{\delta }_{\alpha A}{\delta }_{\beta B}\frac{{\delta }_{i,j+1}}{{\sigma }_{i,i-1}}+{\delta }_{\alpha B}{\delta }_{\beta A}\frac{{\delta }_{i,j-1}}{{\sigma }_{i,i+1}}\right],i,j\ge 2,$$

$$t_{ij}^{\alpha \beta }=\left({\delta }_{\alpha A}{\delta }_{\beta B}{\delta }_{i,j+1}+{\delta }_{\alpha B}{\delta }_{\beta A}{\delta }_{i,j-1}\right)t,i,j\le 2$$

$${t=2\pi {\sigma }_{12}^2{x}^2{K}_{ass}^{(0)}{e}^{G_{00}({\sigma }_{12}^{+})-\beta U^{(C)}({\sigma }_{12}^{+})}{g}_{12}^{00}({\sigma }_{12}^{+})|}_{z_i=0}.$$

$$h_{ij}(r)=\sum _{\alpha \beta }{h}_{ij}^{\alpha \beta }(r)$$

$${4\pi \rho {\sigma }_{12}^2{x}^2{K}_{ass}^{(0)}{e}^{G_{00}({\sigma }_{12}^{+})-\beta U^{(C)}({\sigma }_{12}^{+})}{g}_{12}^{00}({\sigma }_{12}^{+})|}_{z_i=0}+x-1=0$$

$$\Gamma^{2}={\beta e^{2}\over 4\epsilon\epsilon_{0}}\rho\sum_{i=1}^{m}{\bf X}_{i}\mbox{\boldmath${\alpha}$}{\bf X}_{i}^{T},$$

$$X_i^0=\left[z_i-{\eta }^B{\sigma }_i^2\right]{\Gamma }_{{\sigma }_i},$$

$$X_i^{\alpha }={\sigma }_i\left[{\tau }_i^{\alpha }(z)-{\eta }^B{\tau }_i^{\alpha }({\sigma }^2)\right],\alpha \ne 0(\alpha =A,B),$$

$${\eta }^B=\frac{\frac{\pi }{2\Delta }\rho {\sum }_{i=1}^m{\sigma }_i\left\{z_i{\Gamma }_{{\sigma }_i}+{\sigma }_i\left[{\tau }_i^A(z)+{\tau }_i^B(z)\right]\right\}}{1+\frac{\pi }{2\Delta }\rho {\sum }_{i=1}^m{\sigma }_i^2\left[{\sigma }_i{\Gamma }_{{\sigma }_i}+{\tau }_i^A({\sigma }^2)+{\tau }_i^B({\sigma }^2)\right]},$$

$${\tau }_1^A(y)=0,$$

$${\tau }_2^A(y)=\rho {\Gamma }_{{\sigma }_1}{\Gamma }_{{\sigma }_2}{y}_{1}t,$$

$${\tau }_i^A(y)=\frac{{\Gamma }_{{\sigma }_i}{\Gamma }_{{\sigma }_{i-1}}}{2{\sigma }_{i,i-1}}\{\sum _{l=3}^i{y}_{l-1}\left[\frac{\left(1-{\delta }_{li}\right)}{2^{i-l}}\prod _{k=l}^{i-1}\frac{{\sigma }_k{\Gamma }_{{\sigma }_{k-1}}}{{\sigma }_{k,k-1}}+{\delta }_{li}\right]$$

$$+\rho {\sigma }_2{y}_1{\Gamma }_{{\sigma }_1}t\left[\frac{\left(1-{\delta }_{3i}\right)}{2^{i-3}}\prod _{k=3}^{i-1}\frac{{\sigma }_k{\Gamma }_{{\sigma }_{k-1}}}{{\sigma }_{k,k-1}}+{\delta }_{3i}\right]\},3\le i\le m,$$

$${\tau }_1^B(y)=\rho {\Gamma }_{{\sigma }_1}{\Gamma }_{{\sigma }_2}t\sum _{l=2}^m{y}_l\left[\frac{\left(1-{\delta }_{l2}\right)}{2^{l-2}}\prod _{k=2}^{l-1}\frac{{\sigma }_k{\Gamma }_{{\sigma }_{k+1}}}{{\sigma }_{k,k+1}}+{\delta }_{l2}\right]$$

$${\tau }_i^B(y)=\frac{{\Gamma }_{{\sigma }_i}{\Gamma }_{{\sigma }_{i+1}}}{2{\sigma }_{i,i+1}}\sum _{l=i+1}^m{y}_l\left[\frac{\left(1-{\delta }_{l,i+1}\right)}{2^{l-i-1}}\prod _{k=i+1}^{l-1}\frac{{\sigma }_k{\Gamma }_{{\sigma }_{k+1}}}{{\sigma }_{k,k+1}}+{\delta }_{l,i+1}\right],2\le i<m,$$

$${\tau }_m^B(y)=0.$$

$${\tau }_1^A(y)=0,$$

$${\tau }_2^A(y)=\rho {\Gamma }_{{\sigma }_a}{\Gamma }_{{\sigma }_c}{y}_{a}t$$

$${\tau }_3^A(y)=\frac{{\Gamma }_{{\sigma }_n}{\Gamma }_{{\sigma }_c}}{2{\sigma }_{nc}}\left(y_c+\rho {\sigma }_c{y}_a{\Gamma }_{{\sigma }_a}t\right)$$

$${\tau }_i^A(y)=\frac{{\Gamma }_{{\sigma }_n}^2}{2{\sigma }_n}\left[\left(y_c+\rho {\sigma }_c{y}_a{\Gamma }_{{\sigma }_a}t\right)\frac{{\sigma }_n{\Gamma }_{{\sigma }_c}}{2{\sigma }_{nc}}{\left(\frac{{\Gamma }_{{\sigma }_n}}{2}\right)}^{i-4}+y_n\sum _{l=4}^i{\left(\frac{{\Gamma }_{{\sigma }_n}}{2}\right)}^{i-l}\right],4\le i\le m,$$

$${\tau }_1^B(y)=\rho {\Gamma }_{{\sigma }_a}{\Gamma }_{{\sigma }_c}t\left[y_c+y_n\frac{{\sigma }_c}{{\sigma }_{cn}}\sum _{l=3}^m{\left(\frac{{\Gamma }_{{\sigma }_n}}{2}\right)}^{l-2}\right],$$

$${\tau }_2^B(y)=\frac{{\Gamma }_{{\sigma }_c}{\Gamma }_{{\sigma }_n}}{2{\sigma }_{cn}}{y}_n\sum _{l=3}^m{\left(\frac{{\Gamma }_{{\sigma }_n}}{2}\right)}^{l-3},$$

$${\tau }_i^B(y)=\frac{{\Gamma }_{{\sigma }_n}^2}{2{\sigma }_n}{y}_n\sum _{l=i+1}^m{\left(\frac{{\Gamma }_{{\sigma }_n}}{2}\right)}^{l-i-1},3\le i<m,$$

Full text

1 Introduction

Much of success achieved in the statistical-mechanical theory of liquids is due to the availability of models that are simple enough to be described analytically and at the same time are able to reproduce the most important features of a targeted class of systems. In the case of electrolyte solutions, molten salts and liquid metals the models of this type are represented by the so-called ’primitive’ models (PMs). These are the models, which combine short-range hard-sphere repulsion and long-range Coulomb interaction. Sufficiently accurate and simple theoretical description of these models can be achieved in the framework of the mean spherical approximation (MSA). An important advantage of the MSA is that for most of the versions of PMs it can be solved analytically, yielding relatively simple expressions for the structure and thermodynamic properties of the system. Waisman and Lebowitz derived an analytic solution of the MSA for restricted PM (RPM) of electrolytes [1, 2] (equivalent two-component mixture of equal size charged hard spheres) using Laplace transform techniques. This solution was elaborated and extended to the general case of any number of components with arbitrary charges (subject to electroneutrality) and hard-sphere diameters by Blum [3]. He was able to reduce the problem to the solution of only one nonlinear algebraic equation for the famous Blum’s scaling (screening) parameter $\Gamma$, which appears to be the MSA analogue of Debye screening length. All MSA thermodynamic properties can be expressed in terms of this parameter [4].

The major goal of this study is to propose the extension of the PMs of electrolyte solutions for room temperature ionic liquids (RTILs) and develop their theoretical description. Although the vast majority of the previous studies have been focused on the description of a RTIL on the detailed atomistic level (see, e.g.[5]), over the last decade several simple models of RTIL have been proposed and investigated [6, 7, 8, 9, 10, 11, 12, 13]. Molecular ions (cations) in the framework of these models are represented either as a Lennard-Jones/hard spheres with off-center point charges [7, 9, 10, 11, 13], as a Lennard-Jones/hard-sphere dimers with point charges located on one or both sites [6, 9, 12] or as a hard spherocylinder with point charge located on one of its ends [8]. In these papers the structural and dynamic properties [6, 7, 9, 10, 11, 13], density and potential profiles in a planar electrical double layer [12], as well as liquid-solid [10, 11, 13] and gas-liquid [8, 9] phase behavior have been studied. The common feature of all these studies is that the properties of the models were investigated using computer simulation methods only, either Molecular Dynamics (MD) or Monte Carlo (MC) simulations.

We propose to model an RTIL as a two-component mixture of hard-sphere anions and flexible linear chain cations, represented by the tangentially bonded hard spheres with the charge located on one of the terminal beads. For a cation chain length of two, our model reduces to the model studied earlier [12]. More important is that in addition we propose the MSA type of the theory, which is able to provide an analytical description of the structural and thermodynamic properties of the model. The theory is based on the multidensity version of the MSA [14, 15, 16], the so-called associative MSA (AMSA), and represents its extension for chain-forming fluids [17, 18]. We are focused here on the liquid-gas phase behavior of the model, which appears to be one of the most important characteristics of the RTILs. Prediction of the phase diagram and critical parameters for the PM of electrolytes has long been a challenge for the theory and computer simulation since the late sixties and early seventies, when its existence for the RPM had been suggested both theoretically [19, 20, 21] and via computer simulation [22] (for more detailed historical review see Refs. [23, 24, 25]). Due to a number of computer simulation studies, the exact position of the critical point and the phase diagram of the RPM and PMs with different hard-sphere sizes and charge ratios have been recently identified [26, 27, 28]. The situation with theoretical estimates is less satisfactory, i.e. while the position of the critical point is predicted with reasonable accuracy the phase diagram is still too narrow in comparison with computer simulation phase diagrams [29, 30, 31, 32, 33, 34]. Perhaps the most accurate theoretical results for the RPM were obtained by Blum and co-workers [35]. They assumed that due to strong Coulomb interaction cations and anions form dimers and at low enough temperatures the system is completely dimerized. For their theoretical description, a combination of the usual MSA and AMSA was developed and applied. We note in passing that in this paper the authors refer to AMSA as to binding MSA, although both are identical. The theory yields fairly accurate description of the critical temperature and critical density when compared with computer simulation. The scheme suggested by Blum and co-workers was extended for size asymmetric PMs by Qin and Prausnitz [38]. Recently this scheme was used to describe phase behavior of the PMs confined in a disordered porous media by Holovko et al. [36, 37]. In our study we propose further extension of the Blum’s approach and apply it to primitive models of RTILs.

The remainder of the paper is organized as follows. In Section 2 we introduce the model and in Section 3 we discuss the AMSA theory. Solutions of the AMSA for the general version of the model and for two simplified versions are derived in Section 4 and in Section 5 we present expressions for the structure and thermodynamic properties of these models. In Section 6 we consider an extension of Blum’s approach to calculate the phase diagrams and in Section 7 we discuss details of computer simulation approach. Our results and discussion are presented in Section 8 and in Section 9 we collect our conclusions

2 The model

We are modeling an RTIL as a two-component mixture of hard-sphere anions with the number density $\rho_{a}$ and flexible linear chain cations with the number density $\rho_{c}$, represented by $m-1$ tangentially bonded hard spheres with the charge located on one of the terminal beads. The pair potential acting between the particles is represented by the sum of site-site hard-sphere potentials $U^{(hs)}_{ij}(r)$ :

[TABLE]

and Coulomb potential $U^{(C)}_{ij}$:

[TABLE]

which is valid between the charged hard-sphere sites. Here $\sigma_{i}$ is the size of the hard-sphere site $i$, $\delta_{ij}$ is the Kronecker delta, $\epsilon$ and $\epsilon_{0}$ are the dielectric constants of the continuum and vacuum, respectively, $ez_{i}$ is the charge of the site $i$ and the site species indices $i,j$ are taking the values $1,2,\ldots,m$, with $i=1,2$ denoting anion and charged bead of the cation, respectively (see figure 1). We assume also that $|z_{1}|=|z_{2}|=z$ and the total number density of the system is $\rho_{t}=2\rho$, where $\rho=\rho_{a}=\rho_{c}$.

3 Theory

The thermodynamic properties of the model are derived using an appropriately modified Wertheim’s multidensity Orstein-Zernike (OZ) equation supplemented by the associative mean spherical approximation (AMSA) [16] formulated for chain-forming fluids [17, 18]. Our model of an RTIL can be viewed as the complete association limit (CAL) of the $m$-component hard-sphere mixture with two sticky spots (patches) of the type $A$ and $B$, randomly placed on the surface of each particle forming the cation chain. The model at hand will be recovered assuming infinitely strong attraction between the patches of the type $B$ and $A$ located on the surface of the particles of the type $i$ and $i+1$, respectively. Note that the size of each of the patches is small enough to ensure that only one $A-B$ bond is formed. This feature of the model allows us to present the multidensity OZ equation and AMSA closure relation in the following form [17, 18]:

[TABLE]

where ${{\bf{\hat{h}}}}_{ij}(k)$, ${{\bf{\hat{c}}}}_{ij}(k)$ and ${\alpha}$ are the matrices,

[TABLE]

[TABLE]

with the elements ${\hat{h}}^{\alpha\beta}_{ij}(k)$, ${\hat{c}}^{\alpha\beta}_{ij}(k)$ ($\alpha,\beta=0,A,B$) being Fourier transforms of the partial correlation functions $h^{\alpha\beta}_{ij}(r)$, $c^{\alpha\beta}_{ij}(r)$, respectively. Note that $AA$ and $BB$ elements of the matrix ${\alpha}$ are equal 0, which is the consequence of one bond per patch restriction, mentioned above. Here

[TABLE]

[TABLE]

where $E^{\alpha\beta}=\delta_{\alpha 0}\delta_{\beta 0}$, $\beta=1/(k_{B}T)$, $k_{B}$ is Boltzmann constant,

[TABLE]

[TABLE]

and

[TABLE]

Note that in addition to the delta-function term in the closure relations for the correlation functions of the particles forming cation chain, eqs. (5) and (7), the delta-function term appears also in the correlation function for the anions and charged bead of the cations ${\bf c}_{12}(r)$, eqs. (5) and (9). This term is introduced to correct MSA closure for the effects of the ionic association [39]. Here $\left.g^{00}_{ij}(\sigma_{12}^{+})\right|_{z_{i}=0}$ is the contact value of the radial distribution function $g^{00}_{ij}(r)=h^{00}_{ij}(r)+1$ at zero charges on the anion and cation bead, $\left.G_{00}(\sigma_{12}^{+})=g^{00}_{12}(\sigma_{12}^{+})-g^{00}_{12}(\sigma_{12}^{+})\right|_{z_{i}=0}$, $K_{ass}^{(0)}$ is the association constant, $x$ is the fraction of free anions (or cations) and in eq. (9) we are using the exponential approximation [39]. The upper indices $\alpha,\beta$ in the partial correlation functions ${h}^{\alpha\beta}_{ij}(r),{c}^{\alpha\beta}_{ij}(r)$, which enter OZ equation (3), are taking the values $0,A$ and $B$ and denote bonding states of the corresponding particles [15, 17, 18]. The total partial correlation functions $h_{ij}^{\alpha\beta}(r)$ are related to the site-site total correlation functions $h_{ij}(r)$ by

[TABLE]

and the fraction of free anions(cations) $x$ follows from the solution of the mass action law type of equation

[TABLE]

The OZ equation (3), AMSA closure relations (5) and (6) and equation (11) for $x$ form a close set of equations to be solved.

4 Solution of the AMSA

The solution of the AMSA for the models similar to that discussed above have been derived earlier [17, 18] using Blum’s version [3, 40] of the Baxter factorization technique [41]. We shall therefore omit here any details and present only the final expressions, suitable for the model in question. We will consider the general version of the model with different sizes of all monomers in the system and two simplified versions, i.e. one with equal sizes of the neutral beads of the cation and the other with only one neutral cation bead

According to [17, 18, 42], the solution of the set of equations (3), (5), (6) and (11) can be reduced to the solution of one nonlinear algebraic equation for Blum’s screening parameter $\Gamma$

[TABLE]

where ${\bf X}_{i}=\left(X_{i}^{0},X_{i}^{A},X_{i}^{B}\right)$,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here $\beta^{*}=\beta e^{2}/(4\pi\epsilon\epsilon_{0})$, $\Gamma_{\sigma_{i}}=\left(1+\sigma_{i}\Gamma\right)^{-1}$ and $y$ is taking the values either $z$ or $\sigma^{2}$.

4.1 Model with equal sizes of the cation neutral beads (model $C$)

Substantial simplification of the above expressions for $\tau_{i}^{\alpha}(y)$ occurs for the model with equal hard-sphere sizes of all $m-2$ neutral beads of the cation, i.e. $\sigma_{i}=\sigma_{n}$ for $i\geq 3$. We have:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Note that in the above expressions for the sake of convenience we are using a slightly modified notation, i.e $\sigma_{1}=\sigma_{a}$, $\sigma_{2}=\sigma_{c}$, $y_{1}=y_{a}$, $y_{2}=y_{c}$ and $y_{i}=y_{n}\;(i>2)$.

4.2 Model with cation having one neutral bead (model $D$)

Further simplification is possible for the model with the cation represented by the dimer, i.e.

[TABLE]

[TABLE]

Using these expressions for $\tau_{i}^{\alpha}$ in expressions (13), (14) and (15) for $\eta^{B}$ and ${\bf X}_{i}^{\alpha}$, respectively, we have:

[TABLE]

where $q$ is taking the values $a,c,n$ and

[TABLE]

[TABLE]

[TABLE]

5 Structural and thermodynamic properties

The equation for the fraction of free anions $x$ (eq. (11)) includes the contact values of the radial distribution function $g_{12}^{00}(\sigma_{12}^{+})$ of the original version of the model and the version of the model with $z_{i}=0$, i.e $\left.g_{12}^{00}(\sigma_{12}^{+})\right|_{z_{i}=0}$. We have [17, 18]:

[TABLE]

In the framework of the present AMSA closure the contact value of the radial distribution function of the model at zero charges coincide with the Percus-Yevick contact value of the corresponding radial distribution function of the $m$-component mixture of hard spheres $g_{ij}^{(hs)}$, i.e.

[TABLE]

where $\eta=\pi\rho/6\sum_{l=1}^{m}\sigma_{l}^{3}$.

Following [17, 18] for the excess internal energy of the model $\Delta E$ we have

[TABLE]

The Helmholtz free energy of the model $A$ can be written as a sum of four terms, i.e. the ideal gas contribution, a contribution due to bonding or the so-called mass action law contribution [43, 35] $\Delta A^{(MAL)}$ and contributions due to the hard sphere and electrostatic interactions, $\Delta A_{hs}$ and $\Delta A_{el}$, respectively:

[TABLE]

where

[TABLE]

,

[TABLE]

For the hard-sphere contribution we are using the Boublik-Mansoori-Carnahan-Starling-Leland (BMCSL) expression [44, 45], and electrostatic contribution can be calculated numerically, using the coupling constant integration.

All the rest of thermodynamic properties can be obtained using the standard thermodynamic relations. In particular for the pressure $P$ and for the anion and cation chemical potentials $\mu_{a}$ and $\mu_{c}$ we have:

[TABLE]

6 Calculation of the phase diagram

The liquid-gas phase diagram was calculated extending the method proposed earlier [46, 47, 35]. We assume that at the temperatures close to the phase transition all anions and cations are dimerized, so that the system properties can be described using CAL [46], i.e. $K_{ass}^{(0)}\rightarrow\infty$ and $t=1/(2\rho\sigma_{ac})$. This assumption is based on the MC computer simulation observations, which suggest that in the coexisting phases the fraction of nonbonded ions is negligible and the phase diagram of the RPM can be reasonably well represented by the phase diagram of the corresponding fluid of dimers formed by the oppositely charged ions [48, 49, 50, 51]. In addition we follow Blum et al. [35] and in a spirit of Wertheim’s multidensity thermodynamic perturbation theory [14, 15] assume different approximations for different terms in the expression for Helmholtz free energy (23), i.e. we calculate $\Delta A^{(el)}$ using the $\Gamma$ parameter obtained in the complete dissociation limit (CDL) ($K^{(0)}_{ass}=0$ and $t=0$). According to Blum et al. [35] this combination of the CAL approximation for $\Delta A^{(MAL)}$ and CDL approximation for $\Delta A^{(el)}$ can be seen as an $ad\;hoc$ interpolation between the AMSA and the simple interpolation scheme of Stell and Zhou [52], which gives the most accurate prediction for the phase behavior of the PM of electrolytes [35, 38]. Note that in the CDL $\Gamma$ contains contribution due to the presence of the neutral beads of the cation chain and cannot be reduced to regular MSA $\Gamma$ parameter. Taking into account these two assumptions we have:

[TABLE]

[TABLE]

where $G_{00}^{(\infty)}(\sigma_{12}^{+})=G_{00}(\sigma_{12}^{+})|_{K_{ass}^{(0)}\rightarrow\infty}$, and

[TABLE]

where $\left.\Gamma^{(0)}=\Gamma\right|_{K^{(0)}_{ass}=0}$ and $E^{(0)}=E|_{K^{(0)}_{ass}=0}$. The corresponding expression for the pressure $P$ is:

[TABLE]

where

[TABLE]

For the hard-sphere contribution $\Delta P^{(hs)}$ we are using BMCSL expression [44, 45], and electrostatic contribution $\Delta P^{(el)}$ is:

[TABLE]

where $\eta^{B}_{0}=\eta^{B}|_{K_{ass}^{(0)}=0}$.

As usual the phase diagram was calculated from the equality of the pressure and chemical potentials in the coexisting phases.

7 Monte Carlo simulations

Monte Carlo computer simulations were performed in the canonical ensemble using Molsim software [53]. A model system had 100 particles of each type. Electrostatic interactions were calculated using the Ewald technique. During equilibration $4\cdot 10^{6}$ attempted moves per particle were performed and followed by $40-120\cdot 10^{6}$ attempts for production run. Cations were both displaced and rotated simultaneously during a trial move. The excess chemical potential was calculated via Widom’s test particle insertion method. An electroneutral pair (one anion and one cation) was randomly inserted after every 20 trial moves per particle. Inserted test cations were oriented randomly. At each temperature simulations were carried out at a total of 20 densities. Ideal, hard-sphere, and excess contribution were added together to obtain the total chemical potential. Coexisting densities were determined using the Maxwell construction in the $\mu-\rho$ plane. The Widom method is known to perform best at low to moderate densities. We were able to insert linear cations with $m_{c}=m-1=2$ and $m_{c}=3$ at all densities studied. However, a fraction of successful test particle insertions dropped significantly for larger cations ($m_{c}>3$ or $\sigma_{n}>\sigma_{a}$) at larger densities, preventing us from obtaining phase diagrams for these models.

8 Results and discussion

We have studied the phase behavior of two versions of the primitive models of RTIL proposed. The first version (model $D$) is represented by the model with cations modeled by dimers with the neutral bead of different sizes ($\sigma_{a}\leq\sigma_{n}\leq 3\sigma_{a}$) and the second one (model $C$) is a the model with cations represented by the flexible chains of $m_{c}$ tangentially bonded hard- sphere monomers ($1\leq m_{c}\leq 8$) of the same size. In all cases studied the sizes of the anions and charged beads of the cations were chosen to be equal, i.e. $\sigma_{a}=\sigma_{c}=\sigma$. In what follows the model parameters are expressed in terms of the dimensionless quantities: reduced density $\rho^{*}=\rho\sigma_{a}^{3}$ and reduced temperature $T^{*}=4\pi\epsilon\epsilon_{0}k_{B}T/e^{2}$.

In figure 2 we present our theoretical results for the model $D$ with $\sigma_{n}=\sigma,2\sigma,3\sigma$ and in figure 3 for model $C$ with $m_{c}=2,3,5,8$. These results are compared against computer simulation results for the models with $\sigma_{n}=\sigma$ (model $D$) and with $m_{c}=2,3$ (model $C$). In addition, we also show theoretical and computer simulation results for the phase diagram of electrolyte RPM ($\sigma_{n}=0$, $m_{c}=1$). In general theoretical predictions for the RTIL models and for RPM are of the same order of accuracy. The theory gives relatively accurate results for the critical density, predictions for the critical temperature are less accurate. For the RPM, the theory to be around 7% too high the critical temperature and with the increase of the neutral bead size $\sigma_{n}$ (figure 2) or cation chain length $m_{c}$ (figure 3) this disagreement gradually increases.

With the increase of the model asymmetry due to the increase of $\sigma_{n}$ or $m_{c}$ the phase envelope and critical point are moving towards lower temperatures and lower densities. This shift of the phase diagram is reflected in figures 4 and 5, where we show the critical density and critical temperature as a function of the neutral bead size $\sigma_{n}$ for the model $D$ and as a function of the cation chain length $m_{c}$ for the model $C$, respectively. Similar behavior was observed for the primitive model of RTIL with cations represented by spherocylinders [8], i.e. with the increase of spherocylinder length critical temperature and critical density decrease. According to Martin-Betancourt et al. [8] this effect is of entropic origin, since the presence of uncharged tails reduces the number of energetically favorable configurations of the ions. Increasing $\sigma_{n}$ and/or $m_{c}$ reduces the number of configurations in which when anion and charged bead of the cation are in contact. This feature of the model reduces its ability to form clusters and as a result both critical temperature and critical density decrease. Increase of the critical temperature, caused by the increase of dispersion attraction due to chain length increase is a common feature of nonionic fluids such as alkanes or alcohols [54]. This difference in the behavior of the critical temperature is a clear indication that Coulomb interaction is the major driving force of the phase transitions in RTIL and to a first approximation dispersion forces can be neglected (as it is done in primitive models of RTIL). In the same figures we present also the critical packing fraction $\eta_{cr}$ $vs$ $\sigma_{n}$ (figure 4) and $\eta_{cr}$ $vs$ $m_{c}$ (figure 5). For model $C$, $\eta_{cr}$ appears to be almost independent on the cation chain length (including RPM), i.e. computer simulation and theory give $\eta_{cr}\approx 0.04$ and $\eta_{cr}\approx 0.034$, respectively. The same value of the critical packing fraction $\eta_{cr}\approx 0.04$ was obtained for the primitive models of the RTIL studied earlier using computer simulation methods [8, 9]. For model $D$, the situation is different, i.e. here with the increase of the neutral bead size $\sigma_{n}$ for $\sigma_{n}/\sigma_{a}\leq 0.7$ $\eta_{cr}$ slightly decreases and for $0.7<\sigma_{n}/\sigma_{a}\leq 3$ it shows substantial increase.

Finally in figure 6 we compare our results for the critical temperature of the model $C$ against corresponding results obtained by extrapolating the experimental data for the surface tension and density of a homologous series of imidazolium-based ionic liquids, i.e. [Cnmim][BF4], [Cnmim][PF6] and [Cnmim][Ntf2] [55]. Since these ionic liquids are not stable at higher temperatures Rebelo et al. [55] have used the Guggenheim and Eötvos empirical relations [56, 57] to calculate the critical temperature. Our goal here is to verify the ability of the model and theory to give at least qualitatively correct description of the experimentally observed behavior. We have not made any attempts to optimize the model parameters and for the sake of simplicity follow Martin-Betancourt et al. [8] assuming the following values: $\sigma=4\;\AA$ and $L=1.3n_{C}(\AA)$, where $L$ is the chain length and $n_{C}$ is the number of carbons. Here $L$ is the distance between the centers of the terminal beads for completely stretched chain expressed in Angstroms. Note that for our model $L$ can take only values that are multiples of $\sigma$. The critical temperature for the intermediate values of $L$ was calculated via linear interpolation between the values obtained for the chain length $L$ being a multiple of $\sigma$. Although our results are about 60% higher than experimental results, the model and the theory proposed are able to give qualitatively correct behavior, i.e. with the increase of $n_{C}$ both the theoretical and computer simulation critical temperature decreases. For the longer chains this decrease becomes slightly less steep. Further improvement of our results can be achieved by recognizing that the dielectric permittivity $\epsilon$ in the expression for the Coulomb potential (2) can take different values. In the current study we assume that $\epsilon=1$. For imidazolium-based RTIL with ions lacking polar groups the dielectric permittivity is related to refractive index $n_{D}$, i.e. $\epsilon=n^{2}_{D}$. In the recent paper, Lu et al. [13] assume that $\epsilon=2$. Our theory appears to be in a very good agreement with experiment assuming that $\epsilon\approx 1.60-1.67$ (see figure 6).

9 Conclusions

In this paper we proposed several versions of primitive models for RTILs and developed a theory for their description using the analytical solution of the AMSA. Solution of the corresponding multidensity Ornstein-Zernike equation, supplemented by the AMSA, reduces to the solution of only one nonlinear algebraic equation for Blum’s screening parameter $\Gamma$. The theory is used to study the liquid-gas phase behavior of two versions of the model, i.e. models with cations represented by dimers with a neutral bead of different size and by chains with the neutral beads of the same size, respectively. We generated a set of computer simulation results for the liquid-gas phase diagrams of the models. Theoretical predictions for these phase diagrams appear to be in reasonably good agreement with computer simulation predictions. It is demonstrated that the models and theory are able to reproduce experimentally observed trends in the phase behavior of RTILs, in particular, the decrease of the critical temperature and critical density with the increase of the asymmetry of the model in its shape and position of the charge.

10 Acknowledgment

J.R. thanks the Slovenian Research Agency for its financial support through grant P1-0201 and M.H. acknowledges support from the European Union’s Horizon 2020 research and innovation pogramme under the Marie Sklodowska-Curie (grant No.734276) and the State Fund for Fundamental Research (project NF73/3/26-2017).