# Vapour-liquid phase behaviour of primitive models of ionic liquids confined
in disordered porous media

T. Hvozd  $|$ ${}^{\textsf{\footnotesize{\mbox{?}}}}$, T. Patsahan  $|$
${}^{\textsf{\footnotesize{\mbox{?},\penalty 1000\mbox{?}}}}$, Yu. Kalyuzhnyi
$|$ ${}^{\textsf{\footnotesize{\mbox{?},\penalty 1000\mbox{?}}}}$, O. Patsahan
$|$ ${}^{\textsf{\footnotesize{\mbox{?}}}}$, M. Holovko  $|$
${}^{\textsf{\footnotesize{\mbox{?}}}}$ Corresponding author:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Received April 17, 2024, in final form June 12, 2024)

###### Abstract

We develop a theory for the description of ionic liquids (ILs) confined in a
porous medium formed by a matrix of immobile randomly placed uncharged
particles. The IL is modelled as an electroneutral 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 theory
combines a generalization of the scaled particle theory, Wertheim’s
thermodynamic perturbation theory, and the associative mean-spherical
approximation and allows one to obtain analytical expressions for the pressure
and chemical potentials of the matrix–IL system. Using the theory, we
calculate the vapourliquid phase diagrams for two versions of the IL model,
i.e., when the cation is modelled as a dimer and as a chain, in a complete
association limit. The effects of the matrix confinement and of the non-
spherical shape of the cations on the vapour-liquid phase diagrams are
studied.

Key words: ionic liquids, disordered porous medium, chain-like cations,
vapour-liquid phase diagrams

###### Abstract

Ìè ðîçðîáëÿìî òåîðþ äëÿ îïèñó îííèõ ðäèí ó ïîðèñòîìó ñåðåäîâèù, ñôîðìîâàíîìó
ìàòðèöåþ íåðóõîìèõ õàîòèíî ðîçìùåíèõ íåçàðÿäæåíèõ àñòèíîê. îííà ðäèíà
ìîäåëþòüñÿ ÿê åëåêòðîíåéòðàëüíà ñóìø àíîíâ ñôåðèíî¿ ôîðìè òà ëàíöþãîâoïîäáíèõ
êàòîíâ, ïðåäñòàâëåíèõ òàíãåíöàëüíî çâ’ÿçàíèìè òâåðäèìè ñôåðàìè ç çàðÿäîì,
ðîçòàøîâàíèì íà îäíé ç êíöåâèõ ñôåð. Òåîðÿ ïîäíó â ñîá óçàãàëüíåííÿ òåîð¿
ìàñøòàáíî¿ àñòèíêè, òåðìîäèíàìíî¿ òåîð¿ çáóðåíü Âåðòãàéìà òà àñîöàòèâíîãî
ñåðåäíüî-ñôåðèíîãî íàáëèæåííÿ äîçâîëÿ îòðèìàòè àíàëòèí âèðàçè äëÿ òèñêó òà
õiìiíèõ ïîòåíöiàëiâ ñèñòåìè ‘‘ìàòðèöÿîííà ðäèíà’’. Âèêîðèñòîâóþè öþ òåîðþ, ìè
ðîçðàõóâàëè ôàçîâ äàãðàìè ‘‘ãàç-ðiäèíà’’ äëÿ äâîõ âåðñé ìîäåë îííî¿ ðäèíè, à
ñàìå, êîëè êàòîí ìîäåëþòüñÿ ÿê äèìåð ÿê ëàíöþã, â íàáëèæåíí ïîâíî¿ àñîöàö¿.
Äîñëäæåíî âïëèâ íåâïîðÿäêîâàíî¿ ìàòðèö òà íåñôåðèíî¿ ôîðìè êàòîíâ íà ôàçîâ
äàãðàìàõ ‘‘ãàç-ðiäèíà’’.

Ключов слова: îíí ðäèíè, íåâïîðÿäêîâàíå ïîðèñòå ñåðåäîâèùå, ëàíöþãîâîïîäáí
êàòîíè, ôàçîâ äàãðàìè ‘‘ãàç-ðiäèíà’’

##  1 Introduction

We dedicate the article to the memory of Ihor Yukhnovsky, who recently passed
away, and who, starting from the 1950s, made a significant contribution to the
development of the theory of electrolyte solutions. He founded a scientific
school on the theory of ionic systems in Lviv, of which the authors of this
article have the honor to be a part.

Ionic liquids (ILs) are a subclass of molten salts with a melting temperature
lower than 373 K. In general, ILs are composed of organic cations with either
inorganic or organic anions. There are many possibilities of making ILs due to
the variety of ions and the variations in side chains of the ions. ILs and
ionic solutions of complex molecular ions, due to their unique physical and
chemical properties, play an important role in fundamental research and have
an important applied value in modern technological processes [1, 2, 3, 4, 5].
A significant part of technological applications is based on the properties of
ionic liquids in porous materials. In particular, ILs in porous electrodes are
potentially important in innovative electrochemistry [6, 7]. The presence of a
porous medium can substantially affect the properties of ILs including their
phase behavior. The development of a theory for the study and prediction of
the phase behavior of ILs with a complex structure of ions under disordered
porous confinement still remains a relevant problem.

Phase behavior of the bulk ionic solutions with Coulomb-dominated interactions
has been a subject of active research for the last few decades. The simplest
and most frequently used model for such systems is a restricted primitive
model (RPM) which consists of equal numbers of equisized positively and
negatively charged hard spheres (HSs) immersed in a structureless dielectric
continuum. The RPM undergoes a vapour-liquid-like phase transition at low
temperature and at low density (see review [8] and the literature cited
therein). Further studies focused on the modified versions of the RPM that
takes into account charge and/or size asymmetry of ions.

Significant progress in a theoretical description of the vapour-liquid phase
behavior in the above-mentioned models was made within the approaches that
take into account ion association. Among them there is a generalized Debye-
Hückel theory proposed by Fisher and Levin [9, 10], the associative mean
spherical approximation (AMSA) [11, 12, 13], and the binding mean spherical
approximation (BIMSA) [14, 15, 16, 17, 18]. It is worth noting that both the
AMSA and the BIMSA are, in fact, identical. In addition, the method of
collective variables (CVs) appeared to be successful for the description of
these models [19, 20, 21]. The advantage of this approach is that one can
derive an analytical expression for the relevant chemical potential which
includes the effects of higher-order correlations between the ions. The
historical discussion of the vapour-liquid phase behaviour and criticality in
the RPM is available in [8, 22].

Despite a significant progress in understanding the phase behavior of
relatively simple models of electrolyte solutions, today’s challenges require
consideration of more complex ionic systems. Real ILs and ionic solutions are
characterized by a specific intermolecular interaction that is important at
both short and long distances. In particular, all these objects, in addition
to the presence of charges, have a complex molecular structure and geometry of
ions which are far from spherical shape. The location of charges is also
important, i.e., the charge is frequently not located at the center of mass of
molecular ions. These features can have a significant impact on the properties
of ILs. This also complicates the development of theoretical approaches
capable of properly predicting the properties of these systems.

Several simple IL models have been proposed over the past fifteen years [23,
24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. Molecular ions within the
framework of these models are represented either as hard spheres with off-
center point charges [24, 28, 29, 32], as dimers with point charges located on
one or both ends [23, 26, 27, 28, 34], or as hard spherocylinders with a point
charge located on one of their ends [25]. These models are mostly studied
using computer simulation methods. The vapour-liquid phase behavior of model
ILs with chain-like molecular ions was studied theoretically in [31, 34].
However, these studies were limited to the bulk case.

The most frequently used model for the study of fluids in disordered porous
media is the model proposed by Madden and Glandt [35]. In this model, the
porous medium is represented as a quenched disordered matrix of HSs while the
fluid is distributed inside the matrix. In this case, statistical-mechanical
averages used for calculations of thermodynamic properties become double
ensemble averages: the first average is taken over all degrees of freedom of
fluid particles keeping the quenched particles fixed, and the other average is
performed over all realizations of a matrix. Using this approach, the scale
particle theory (SPT) was developed and analytical expressions for the
thermodynamic functions of the HS fluid in a disordered HS matrix were
obtained [36, 37, 38, 39, 40]. This allowed one to use the model of a HS fluid
in a disordered HS matrix as a reference system (RS). Therefore, to study the
effect of the presence of a porous medium on the phase behavior in the RPM
fluid, both the CV method [41] and the AMSA approach [22] were used. It was
shown that the presence of a porous medium results in a shift of the vapour-
liquid phase diagram of the RPM towards lower densities and temperatures.
Later on SPT was generalized and applied to the study of the properties of the
multicomponent HS fluid confined in the HS disordered matrix [42]. This
version of the SPT was used to extend the theory proposed in [41, 22] for the
case of primitive models with ions of a spherical shape and different sizes
[43, 44], as well as the models involving oppositely charged spherical ions
differing in both size and valence [45].

In this work, we study the vapour-liquid phase diagrams of ILs in a disordered
porous medium. Following [34], we model an IL as a two-component mixture of HS
anions and flexible linear chain cations, represented by the tangentially
bonded HSs with the charge located on one of the terminal beads. The existing
theoretical approaches used for describing the chain ions provide a fairly
good description of the thermodynamics of these systems in the bulk case.
However, in the presence of a disordered matrix, this problem is much more
difficult and has remained unsolved until recently. We propose a theoretical
approach that allows one to obtain the thermodynamic functions of the system
of chain ions in a disordered porous matrix. The approach combines the
extension of the AMSA theory for chain-forming fluids [46, 47, 48, 49, 11] and
the extension of the SPT theory for HS fluids in a disordered HS matrix [50,
51, 52, 53, 37, 40]. Using this approach we calculate the vapour-liquid phase
diagrams of the IL models with molecular ions, represented by homo- and
hetero-nuclear charged chains confined in a disordered HS matrix of different
porosities.

The remainder of the paper is organised as follows. In section 2, we present
the models and develop the theoretical formalism. The results are presented
and discussed in section 3. We conclude in section 4.

##  2 Models and theory

###  2.1 Models

We present an IL as an electroneutral mixture of cations and anions with the
number densities $\rho_{c}$ and $\rho_{a}$ ($\rho_{i}=N_{i}/V$), respectively,
immersed in a structureless dielectric continuum. The anions are modelled as
charged HSs with diameter $\sigma_{1}=\sigma_{-}$ and charge $ez_{1}=ez_{-}$.
The cations are represented by $m-1$ tangentially bonded HS monomers with the
charge located on one of the terminal beads. The charged monomer of the cation
has a diameter $\sigma_{2}=\sigma_{+}$ and a charge $ez_{2}=ez_{+}$, while all
uncharged monomers have the same diameter $\sigma_{n}$. It is assumed that the
charges are in the center of the spheres.

The pair potential acting between the particles is represented by the sum of
site-site HS potentials $U_{ij}^{\text{hs}}(r)$:

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

and Coulomb potential $U_{ij}^{\text{C}}(r)$:

$U_{ij}^{\text{C}}(r)=\frac{e^{2}z_{i}z_{j}}{\varepsilon r},$ (2.4)

between charged HSs. Here, $\varepsilon$ is the dielectric constant of the
continuum. In general, we assume that

$\sigma_{+}=\sigma_{-}=\sigma\neq\sigma_{n}$

and $|z_{1}|=|z_{2}|=z$. The total number density of the system is
$\rho_{t}=2\rho$, where $\rho=\rho_{-}=\rho_{+}$. For
$\sigma_{+}=\sigma_{-}=\sigma=\sigma_{n}$, the cation model reduces to a homo-
nuclear chain. Two versions of the above described model of ILs are presented
in figure 1.

Figure 1: (Colour online) Two models of ILs distinguished by the model for the
cation: hetero-nuclear dimer (model A) and mono-nuclear chain (model B). In
both cases, the anion is modelled as a single charged hard sphere.

The IL model is confined in a disordered porous matrix formed by hard spheres
of diameter $\sigma_{0}$. The interaction potentials between an ion and a
matrix particle $U_{0i}^{\text{hs}}(r)$ ($i=+,-,n$) and between two matrix
particles $U_{00}^{\text{hs}}(r)$ are described by the HS potentials:

$\displaystyle
U_{0i}^{\text{hs}}(r)=\left\\{\begin{array}[]{ll}\infty,&r<\sigma_{0i}\\\
0,&r\geqslant\sigma_{0i}\end{array}\right.\,,\qquad
U_{00}^{\text{hs}}(r)=\left\\{\begin{array}[]{ll}\infty,&r<\sigma_{0}\\\
0,&r\geqslant\sigma_{0}\end{array}\right..$ (2.9)

In (2.9), $\sigma_{0i}=\frac{1}{2}\left(\sigma_{0}+\sigma_{i}\right)$ and, in
general, $\sigma\neq\sigma_{n}\neq\sigma_{0}$.

The matrix is characterised by different types of porosity, namely,
geometrical porosity $\phi_{0}$ and probe-particle porosity $\phi_{i}$ for the
$i$th species of the HS mixture. The probe-particle porosity $\phi_{i}$ is
defined by the excess value of the chemical potential of a fluid particle with
diameter $\sigma_{i}$ in the limit of infinite dilution and, hence, takes into
account the size of adsorbate species [42]. Explicit expressions for
$\phi_{0}$ and $\phi_{i}$ will be given below. In order to distinguish a pure
geometrical effect of confinement, we restrict the model only to a hard-core
interaction between the fluid and matrix particles with the latter being
uncharged.

We present Helmholtz free energy $f=F/V$ of the matrix-IL system [equations
(2.3)–(2.9)] in the form:

$\beta f=\frac{\beta F}{V}=\beta f^{(\text{id})}+\beta
f^{(\text{ref})}+\beta\Delta f,$

where $\beta f^{\text{id}}$ is an ideal-gas contribution, $f^{\text{ref}}$ is
the free energy of a reference system (RS), which is represented by a multi-
component HS fluid confined in the HS matrix. $\Delta f$ is the part connected
with the ionic subsystem, and $\beta=1/k_{\text{B}}T$. Below we present
expressions for these contributions to the free energy and the corresponding
contributions to the pressure and chemical potentials. Here, we provide a
brief description of each theoretical approach upon which our theory is based.
A detailed presentation of the generalized SPT theory and the extension of the
AMSA theory for chain-forming fluids can be found elsewhere [54, 40, 42, 11,
55, 56, 57, 58].

###  2.2 Reference system

To calculate thermodynamic properties of the RS we use a generalized version
of the SPT theory, in particular, the approach, proposed recently for a two-
component mixture of hard convex bodies confined in a disordered HS matrix
[59]. In order to characterize each particle of the system, we use three
geometrical parameters: the volume $V$ of a particle, its surface area $S$ and
the mean curvature $r$ taken with a factor ${1}/{4\piup}$ [54]. For a
multicomponent HS mixture and for HS matrix particles we have:

$V_{i}=\frac{1}{6}\piup\sigma_{i}^{3}\;,\;\;\;S_{i}=\piup\sigma_{i}^{2}\;,\;\;\;r_{i}=\sigma_{i}/2\;,\;\;i=0,\ldots,m.$
(2.10)

The geometrical porosity depends on a matrix structure and is related to the
volume of a void existing between the matrix particles, i.e.,

$\phi_{0}=1-\eta_{0},$ (2.11)

where $\eta_{0}=\piup\rho_{0}\sigma_{0}^{3}/6$ is the packing fraction of the
matrix particles. The second type of porosity is defined by the chemical
potential of a fluid in the limit of infinite dilution and it is called a
probe-particle porosity $\phi$. This porosity characterises the adsorption of
a fluid particle in an empty matrix. For a multicomponent HS fluid confined in
a hard-sphere matrix, we have the probe-particle porosity $\phi_{i}$ for each
species $i$. The probe-particle porosity $\phi$ depends only on the nature of
the $i$th species and has the form [42]:

$\displaystyle\phi_{i}$ $\displaystyle=$
$\displaystyle(1-\eta_{0})\exp\bigg{[}-3k_{i0}\left(1+k_{i0}\right)\frac{\eta_{0}}{1-\eta_{0}}-\frac{9}{2}k_{i0}^{2}\frac{\eta_{0}^{2}}{(1-\eta_{0})^{2}}$
$\displaystyle-$ $\displaystyle
k_{i0}^{3}\frac{\eta_{0}}{(1-\eta_{0})^{3}}\left(1+\eta_{0}+\eta_{0}^{2}\right)\bigg{]},$

where $k_{i0}=\sigma_{i}/\sigma_{0}$.

We present the pressure of the HS fluid mixture confined in a disordered HS
matrix in the form [59]:

$\displaystyle\beta P^{\rm{ref}}=\beta P^{\rm{SPT2b3}^{\star}}+\beta\Delta
P^{\rm{CS}},$ (2.12)

where $\beta P^{\rm{SPT2b3}^{\star}}$ is the contribution due to the SPT2b3*
approximation [40]. This approximation provides a correct description of the
thermodynamic properties of the RS in a wide range of a fluid packing fraction
$\eta=\piup\sum_{i=1}^{m}\rho_{i}\sigma_{i}^{3}/6$, starting from the smallest
ones and ending with the packing fractions close to the maximum values
possible in a given matrix. $\beta\Delta P^{\rm{CS}}$ is the Carnahan-Starling
correction for the case of a mixture in a porous medium [53]. Following [59],
we obtain:

$\displaystyle\frac{\beta P^{\rm{SPT2b3}^{\star}}}{\rho_{m}}$ $\displaystyle=$
$\displaystyle\frac{1}{1-\eta/\phi_{0}}+\frac{A}{2}\frac{\eta/\phi_{0}}{\left(1-\eta/\phi_{0}\right)^{2}}+\frac{2B}{3}\frac{\left(\eta/\phi_{0}\right)^{2}}{\left(1-\eta/\phi_{0}\right)^{3}}$
(2.13) $\displaystyle+$
$\displaystyle\frac{\phi_{0}-\phi^{\star}}{\phi^{\star}}\frac{\phi_{0}}{\eta}\left[\ln(1-\eta/\phi_{0})+\frac{\eta/\phi_{0}}{1-\eta/\phi_{0}}\right]$
$\displaystyle+$
$\displaystyle\frac{\phi^{\star}-\phi}{\eta}\left[\ln(1-\eta/\phi^{\star})+\frac{\eta/\phi^{\star}}{1-\eta/\phi^{\star}}\right],$
$\displaystyle\frac{\beta\Delta
P^{\rm{CS}}}{\rho_{m}}=-\frac{\left(\eta/\phi_{0}\right)^{3}}{\left(1-\eta/\phi_{0}\right)^{3}}\frac{q_{m}s_{m}^{2}}{9v_{m}^{2}},$
(2.14)

where
$\phi^{\star}=\frac{\phi_{0}\phi}{\phi_{0}-\phi}\ln\frac{\phi_{0}}{\phi}$ is
the maximum value of the packing fraction of the fluid in the disordered
matrix [39], $\phi$ is the total thermodynamic porosity of the matrix for a
given fluid mixture:

$\frac{1}{\phi}=\frac{1}{\eta}\sum_{i=1}^{m}\frac{\rho_{i}V_{i}}{\phi_{i}}.$
(2.15)

$A$ and $B$ are given by

$A=\sum_{i=1}^{m}x_{i}a_{i},\qquad B=\sum_{i=1}^{m}x_{i}b_{i},$ (2.16)

where the coefficients $a_{i}$ and $b_{i}$ define the porous medium structure
[42, 59], the expressions for them are too cumbersome to be presented here
(they can be found in [59]), $x_{i}=\rho_{i}/\rho_{m}$ denotes the fluid
composition, and $v_{m}$, $s_{m}$, $q_{m}$, and $\rho_{m}$ have the form [53]:

$\displaystyle
v_{m}=\sum_{i=1}^{m}x_{i}V_{i}\,,\,\,\,\,s_{m}=\sum_{i=1}^{m}x_{i}S_{i}\,,\,\,\,\,q_{m}=\sum_{i=1}^{m}x_{i}r_{i}^{2}\,,\,\,\,\rho_{m}=\sum_{i=1}^{m}\rho_{i},$
(2.17)

where $V_{i}$, $S_{i}$, and $r_{i}$ are given in (2.10), $\rho_{i}$ is the
number density of the $i$th species of the HS mixture.

Similarly, we present the partial chemical potential $\beta\mu_{i}$ as
follows:

$\displaystyle\beta\mu_{i}^{(\text{ref})}=\beta\mu_{i}^{\rm{SPT2b3^{*}}}+\beta\Delta\mu_{i}^{\rm{CS}},$
(2.18) $\displaystyle\beta\mu_{i}^{\rm{SPT2b3}^{\star}}$ $\displaystyle=$
$\displaystyle\beta\mu_{i}^{\rm{SPT2a}}+\frac{\eta(\phi_{0}-\phi^{\star})}{\phi_{0}\phi^{\star}\left(1-\eta/\phi_{0}\right)}+\frac{\eta\left(\phi^{\star}-\phi\right)}{\phi^{\star}\phi^{\star}\left(1-\eta/\phi^{\star}\right)}$
(2.19) $\displaystyle+$ $\displaystyle\left(\frac{\rho
V_{i}}{\eta}-1\right)\bigg{[}\frac{\phi_{0}-\phi}{\eta}\ln(1-\eta/\phi_{0})+\frac{\phi(\phi_{0}-\phi^{\star})}{\phi_{0}\phi^{\star}\left(1-\eta/\phi_{0}\right)}$
$\displaystyle+$
$\displaystyle\frac{\phi\left(\phi^{\star}-\phi\right)}{\phi^{\star}\phi^{\star}\left(1-\eta/\phi^{\star}\right)}\bigg{]}-\frac{\rho
V_{i}}{\eta}\left(\frac{\phi}{\phi_{i}}-1\right)\bigg{[}\frac{\phi}{\eta}\ln\left(1-\eta/\phi_{0}\right)$
$\displaystyle-$
$\displaystyle\frac{\phi(\phi_{0}-\phi^{\star})}{\phi_{0}\phi^{\star}\left(1-\eta/\phi_{0}\right)}-\frac{\phi\left(\phi^{\star}-\phi\right)}{\phi^{\star}\phi^{\star}\left(1-\eta/\phi^{\star}\right)}+1\bigg{]},$
$\displaystyle\beta\Delta\mu_{i}^{\rm{CS}}$ $\displaystyle=$
$\displaystyle-\frac{V_{i}}{v_{m}}\frac{\left(\eta/\phi_{0}\right)^{3}}{\left(1-\eta/\phi_{0}\right)^{3}}\frac{q_{m}s_{m}^{2}}{9v_{m}^{2}}+\frac{s_{m}}{9v_{m}^{3}}\big{[}\left(q_{i}s_{m}+2S_{i}q_{m}\right)v_{m}-2V_{i}q_{m}s_{m}\big{]}$
(2.20) $\displaystyle\times$
$\displaystyle\bigg{[}\ln(1-\eta/\phi_{0})+\frac{\eta/\phi_{0}}{1-\eta/\phi_{0}}-\frac{1}{2}\frac{\left(\eta/\phi_{0}\right)^{2}}{\left(1-\eta/\phi_{0}\right)^{2}}\bigg{]},$

where $\beta\mu_{i}^{\rm{SPT2a}}$ is given in [60]. The other notations are
the same as in (2.13)–(2.14).

###  2.3 Ionic subsystem

We consider an ionic subsystem consisting of an electroneutral mixture of
chain cations and monomeric anions in a structureless dielectric medium. The
model can be considered as a partial case of the $m$-component HS mixture with
two “sticky” spots (patches) of the type $A$ and $B$, randomly placed on the
surface of each particle forming the cation chain, simultaneously assuming an
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. Importantly, the size of each patch is small enough to ensure
that only one $A$-$B$ bond can be created. Thermodynamic properties for such a
generalized model were obtained using a modified Wertheim’s multidensity
Ornstein-Zernike (OZ) equation supplemented by the AMSA theory [11] formulated
for the case of chain-forming fluids [55, 57]. Here, we present only the final
expressions for the IL model under consideration.

The solution of the AMSA theory for our IL model can be reduced to the
solution of a nonlinear algebraic equation for Blum’s screening parameter
$\Gamma$ [55, 56, 57, 34]:

$\Gamma^{2}=\frac{\piup\beta
e^{2}}{\varepsilon}\rho\sum_{i=1}^{m}\mathbf{X}_{i}\mathbf{\alpha}\mathbf{X}_{i}^{T},$
(2.21)

where $\mathbf{X}_{i}=(X_{i}^{0},X_{i}^{A},X_{i}^{B})$ is a row matrix and
$\mathbf{\alpha}$ is a symmetrical $3\times 3$ matrix

$\displaystyle\mathbf{\alpha}=\begin{pmatrix}1&1&1\\\ 1&0&1\\\
1&1&0\end{pmatrix}.$

A general form of the matrix elements ${X}_{i}^{\alpha}$ with $i=1,2,\ldots,m$
and $\alpha=0,A,B$ was obtained in [34].

As will be shown below, the elements of the matrix $\mathbf{X}_{i}$ contain
the quantity $t$

$t=2\piup\sigma_{12}^{2}x^{2}K_{\text{as}}^{(0)}\exp\left(G_{00}(\sigma_{12}^{+})-\beta
U^{(C)}(\sigma_{12}^{+})\right)g_{12}^{00}(\sigma_{12}^{+})|_{z_{i}=0},$
(2.22)

where $g_{12}^{00}(\sigma_{12}^{+})|_{z_{i}=0}$ is the contact value of the
radial distribution function at zero charges on the anion and cation bead,
$G_{00}(\sigma_{12}^{+})=g_{12}^{00}(\sigma_{12}^{+})-g_{12}^{00}(\sigma_{12}^{+})|_{z_{i}=0}$,
$K_{\text{as}}^{(0)}$ is the association constant, $x$ is the fraction of free
anions (or cations) determined from the solution of the mass action law type
of equation

$2\rho t\sigma_{12}+x-1=0.$ (2.23)

In order to obtain (2.22), the exponential approximation [58] is used.

We consider two versions (model A and model B) of the IL model which differ by
the cation shape (see figure 1). In model A, the cations are modelled as a
hetero-nuclear dimers ($\sigma_{-}=\sigma_{+}=\sigma\neq\sigma_{n}$), where
indices $-,+,n$ indicate anion, cation, and neutral monomer, respectively) and
in model B the cations are modelled as homo-nuclear chains
($\sigma_{1}=\sigma_{2}=\ldots=\sigma_{m}=\sigma$). In both models, the anions
are modelled as single charged HSs. The expressions for ${X}_{i}^{\alpha}$ for
each of the models are presented below.

#### Model A: the cation modelled as hetero-nuclear dimer.

In this case, the matrix elements $\mathbf{X}_{i}$ have the form:
for $\mathbf{X}_{1}=\mathbf{X}_{-}$,

$\displaystyle X_{-}^{0}$ $\displaystyle=$
$\displaystyle-\Gamma_{\sigma_{-}}(z+\eta^{B}\sigma_{-}^{2}),\qquad
X_{-}^{A}=0,$ (2.24) $\displaystyle X_{-}^{B}$ $\displaystyle=$
$\displaystyle\Gamma_{\sigma_{-}}\Gamma_{\sigma_{+}}\sigma_{-}\rho
t\bigg{(}z-\eta^{B}\sigma_{+}^{2}-\eta^{B}\Gamma_{\sigma_{n}}\frac{\sigma_{+}\sigma_{n}^{2}}{2\sigma_{+n}}\bigg{)},$

for $\mathbf{X}_{2}=\mathbf{X}_{+}$,

$\displaystyle X_{+}^{0}$ $\displaystyle=$
$\displaystyle\Gamma_{\sigma_{+}}(z-\eta^{B}\sigma_{+}^{2}),$ (2.25)
$\displaystyle X_{+}^{A}$ $\displaystyle=$
$\displaystyle-\Gamma_{\sigma_{-}}\Gamma_{\sigma_{+}}\sigma_{+}\rho
t(z+\eta^{B}\sigma_{-}^{2}),$ $\displaystyle X_{+}^{B}$ $\displaystyle=$
$\displaystyle-\Gamma_{\sigma_{+}}\Gamma_{\sigma_{n}}\eta^{B}\frac{\sigma_{+}\sigma_{n}^{2}}{2\sigma_{+n}},$

for $\mathbf{X}_{3}=\mathbf{X}_{n}$,

$\displaystyle X_{n}^{0}$ $\displaystyle=$
$\displaystyle-\Gamma_{\sigma_{n}}\eta^{B}\sigma_{n}^{2},$ $\displaystyle
X_{n}^{A}$ $\displaystyle=$
$\displaystyle\frac{\Gamma_{\sigma_{+}}\Gamma_{\sigma_{n}}\sigma_{n}}{2\sigma_{+n}}\left[z-\eta^{B}\sigma_{+}^{2}-\rho
t\Gamma_{\sigma_{-}}\sigma_{+}(z+\eta^{B}\sigma_{-}^{2})\right],$
$\displaystyle X_{n}^{B}$ $\displaystyle=$ $\displaystyle 0.$

The expression for $\eta^{B}$ is as follows:

$\eta^{B}={{\piup\rho\over
2\Delta}\left\\{\sigma_{+}\Gamma_{\sigma_{+}}-\sigma_{-}\Gamma_{\sigma_{-}}+\left[{\sigma_{n}^{2}\Gamma_{\sigma_{n}}\over
2\sigma_{n+}}+\rho\Gamma_{\sigma_{-}}t\left(\sigma_{-}^{2}-\sigma_{+}^{2}-{\sigma_{+}\sigma_{n}^{2}\Gamma_{\sigma_{n}}\over
2\sigma_{n+}}\right)\right]\Gamma_{\sigma_{+}}\right\\}\over 1+{\piup\rho\over
2\Delta}\left\\{\sum_{i=-}^{n}\sigma_{i}^{3}\Gamma_{\sigma_{i}}+\sigma_{+}^{2}\left[{\sigma_{n}^{2}\Gamma_{\sigma_{n}}\over
2\sigma_{n+}}+\sigma_{-}^{2}\rho\Gamma_{\sigma_{-}}t\left(2+{\sigma_{n}^{2}\Gamma_{\sigma_{n}}\over\sigma_{+n}\sigma_{+}}\right)\right]\Gamma_{\sigma_{+}}\right\\}},$
(2.26)

where $\sigma_{+n}=(\sigma_{+}+\sigma_{n})/2$, $\Delta$ and
$\Gamma_{\sigma_{i}}$ have the form:

$\Delta=1-\piup\rho\sum_{i=1,\ldots,m}\sigma_{i}^{3}/6,$ (2.27)
$\Gamma_{\sigma_{i}}=\left(1+\sigma_{i}\Gamma\right)^{-1}.$ (2.28)

#### Model B: a chain cation with monomers of the same size.

In the case where $m$ monomers forming a chain cation are HSs of the same
diameter, the elements ${X}_{i}^{\alpha}$ of the matrix $\mathbf{X}_{i}$ can
be presented in the following compact form [34]:

$\displaystyle X_{i}^{0}$ $\displaystyle=$
$\displaystyle\left[z_{i}-\eta^{B}\sigma_{i}^{2}\right]\Gamma_{\sigma_{i}},$
(2.29) $\displaystyle X_{i}^{\alpha}$ $\displaystyle=$
$\displaystyle\sigma_{i}\left[\tau_{i}^{\alpha}(z)-\eta^{B}\tau_{i}^{\alpha}(\sigma^{2})\right],\;\;\;\;\;\;\alpha\neq
0\;\;(\alpha=A,B),$ $\displaystyle\eta^{B}$ $\displaystyle=$
$\displaystyle{{\piup\over
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\\}\over
1+{\piup\over
2\Delta}\rho\sum_{i=1}^{m}\sigma_{i}^{2}\left[\sigma_{i}\Gamma_{\sigma_{i}}+\tau^{A}_{i}(\sigma^{2})+\tau_{i}^{B}(\sigma^{2})\right]},$

where the following notations are introduced:

$\tau^{A}_{1}(y)=0,\quad\tau_{2}^{A}(y)=\rho\Gamma_{\sigma_{-}}\Gamma_{\sigma_{+}}y_{-}t,$
$\tau_{3}^{A}(y)={\Gamma_{\sigma_{n}}\Gamma_{\sigma_{+}}\over
2\sigma_{n+}}\left(y_{+}+\rho\sigma_{+}y_{-}\Gamma_{\sigma_{-}}t\right),$
$\tau_{i}^{A}(y)={\Gamma_{\sigma_{n}}^{2}\over
2\sigma_{n}}\left[\left(y_{+}+\rho\sigma_{+}y_{-}\Gamma_{\sigma_{-}}t\right){\sigma_{n}\Gamma_{\sigma_{+}}\over
2\sigma_{n+}}\left({\Gamma_{\sigma_{n}}\over
2}\right)^{i-4}+y_{n}\sum_{l=4}^{i}\left({\Gamma_{\sigma_{n}}\over
2}\right)^{i-l}\right],4\leqslant i\leqslant m,$
$\tau_{1}^{B}(y)=\rho\Gamma_{\sigma_{-}}\Gamma_{\sigma_{+}}t\left[y_{+}+y_{n}{\sigma_{+}\over\sigma_{+n}}\sum_{l=3}^{m}\left({\Gamma_{\sigma_{n}}\over
2}\right)^{l-2}\right],$
$\tau_{2}^{B}(y)={\Gamma_{\sigma_{+}}\Gamma_{\sigma_{n}}\over
2\sigma_{+n}}y_{n}\sum_{l=3}^{m}\left({\Gamma_{\sigma_{n}}\over
2}\right)^{l-3},$ $\tau_{i}^{B}(y)={\Gamma_{\sigma_{n}}^{2}\over
2\sigma_{n}}y_{n}\sum_{l=i+1}^{m}\left({\Gamma_{\sigma_{n}}\over
2}\right)^{l-i-1},\;\;\;\;\;3\leqslant i<m,\quad\tau_{m}^{B}(y)=0,$

and $y$ is taking on the values either $z$ or $\sigma^{2}$.
$\Gamma_{\sigma_{i}}$ and $\Delta$ entering the above equations are presented
in (2.27) and (2.28), respectively.

###  2.4 Structural properties

Equations (2.22)–(2.23) include contact values of the radial distribution
function $g_{+-}^{00}(\sigma_{+-})$ between oppositely charged ions of the
matrix–IL model, which includes charges and the model with $z_{i}=0$,
respectively [55, 57, 34]

$g_{+-}^{00}(\sigma_{+-})=g_{+-}^{00}(\sigma_{+-})|_{z_{i}=0}-\frac{\beta
e^{2}}{\varepsilon\sigma_{ij}}X_{+}^{0}X_{-}^{0},$ (2.30)

where the expressions for $X_{+}^{0}$ and $X_{-}^{0}$ are the same for both
models (see (2.24), (2.25), and (2.29)).

The contact value of the radial distribution function of the model at zero
charge within the framework of the AMSA theory coincides with the contact
value of the corresponding radial distribution function of the $m$-component
HS mixture. For our model, this will be the contact value of the radial
distribution function of a two-component HS mixture in the presence of
immobile matrix particles. It can be obtained by generalizing the result of
[61, 62] for the case of a two-component HS fluid in a disordered porous
medium. According to the SPT, the contact value of the radial distribution
function between a scaled particle of infinitesimally small size $\sigma_{s}$
and the $i$th species particle of HS mixture confined in a porous matrix can
be presented as follows [36, 42]:

$g_{is}^{\text{hs}}(\sigma_{is})=\frac{1}{p_{0}^{i}(\sigma_{s})-\sum_{b}\eta_{b}(1+\sigma_{s}/\sigma_{i})^{3}},$

where $\sigma_{\alpha s}=(\sigma_{\alpha}+\sigma_{s})/2$,
$p_{0}^{\alpha}(\sigma_{s})$ is the probability of finding a cavity created by
the scaled particle in the matrix in the absence of fluid particles [38]. For
a point scaled particle, $p_{0}^{\alpha}(\sigma_{s}=0)$ is the geometric
porosity of the matrix $\phi_{0}$ [42]. To obtain the contact value of the
radial distribution function between two fluid particles of type $i$ and $j$,
we follow [61, 59] and express $g_{is}^{\text{hs}}(\sigma_{is})$ as

$g_{is}^{\text{hs}}(\sigma_{is})=G_{is}^{(0)}+G_{is}^{(1)}\frac{\sigma_{s}}{\sigma_{s}+\sigma_{i}}+\frac{1}{2}G_{is}^{(2)}\left(\frac{\sigma_{s}}{\sigma_{s}+\sigma_{i}}\right)^{2},$

where $G_{is}^{(0)}$, $G_{is}^{(1)}$ and $G_{is}^{(2)}$ are found from the
continuity of $g_{is}^{\text{hs}}(\sigma_{is}^{+})$ and the first and second
derivatives with respect to $\sigma_{s}$ at $\sigma_{s}=0$.

As a result, we get for $g_{+-}^{(\text{hs})}(\sigma_{+-})$

$\displaystyle g_{+-}^{(\text{hs})}(\sigma_{+-})$ $\displaystyle=$
$\displaystyle\frac{1}{\phi_{0}-\eta}+\frac{3}{2}\frac{1}{(\phi_{0}-\eta)^{2}}\left(k_{10}\eta_{0}+\eta_{-}+\eta_{+}+\frac{1}{k_{1}}\eta_{n}\right)$
(2.31) $\displaystyle+$
$\displaystyle\frac{1}{2}\frac{1}{(\phi_{0}-\eta)^{3}}\left(k_{10}\eta_{0}+\eta_{-}+\eta_{+}+\frac{1}{k_{1}}\eta_{n}\right)^{2},$

where $\phi_{0}$ is the geometrical porosity of a HS matrix (2.11),
$\eta=\sum_{i=-,+,n}\eta_{i}$, $\eta_{i}=\piup\rho_{i}\sigma_{i}^{3}/6$ is the
packing fraction of the $i$th species of HS fluid mixture. For the model with
a hetero-nuclear cation, $\sigma_{-}=\sigma_{+}=\sigma\neq\sigma_{n}$ and we
have:

$\eta_{-}=\eta_{+},\quad k_{10}=\sigma_{0}/\sigma,\quad
k_{1}=\sigma_{n}/\sigma.$ (2.32)

Similarly, we obtain the expression for the contact value of the radial
distribution function $g_{+n}^{(\text{hs})}(\sigma_{+n})$

$\displaystyle g_{+n}^{(\text{hs})}(\sigma_{+n})$ $\displaystyle=$
$\displaystyle\frac{1}{\phi_{0}-\eta}+\frac{3}{2}\frac{1}{(\phi_{0}-\eta)^{2}}\left[k_{10}\eta_{0}+\frac{2}{1+k_{1}}\left(k_{1}(\eta_{-}+\eta_{+})+\eta_{n}\right)\right]$
(2.33) $\displaystyle+$
$\displaystyle\frac{1}{2}\frac{1}{(\phi_{0}-\eta)^{3}}\left[k_{10}\eta_{0}+\frac{2}{1+k_{1}}\left(k_{1}(\eta_{-}+\eta_{+})+\eta_{n}\right)\right]^{2},$

where the notations are the same as in (2.31).

The contact value of the radial distribution function
$g_{nn}^{(\text{hs})}(\sigma_{nn}^{+})$, has the form:

$\displaystyle g_{nn}^{(\text{hs})}(\sigma_{nn}^{+})$ $\displaystyle=$
$\displaystyle\frac{1}{\phi_{0}-\eta}+\frac{3}{2}\frac{1}{(\phi_{0}-\eta)^{2}}\left(k_{10}\eta_{0}+k_{1}(\eta_{-}+\eta_{+})+\eta_{n}\right)$
(2.34) $\displaystyle+$
$\displaystyle\frac{1}{2}\frac{1}{(\phi_{0}-\eta)^{3}}\left(k_{10}\eta_{0}+k_{1}(\eta_{-}+\eta_{+})+\eta_{n}\right)^{2}.$

For the model of chain cations with monomers of the same diameter $k_{1}=1$.
In this case, from (2.31)–(2.34) we have
$g_{+-}^{(\text{hs})}(\sigma)=g_{+n}^{(\text{hs})}(\sigma)=g_{-n}^{(\text{hs})}(\sigma)$.

###  2.5 Thermodynamic properties

We write the excess internal energy $\Delta E$ of our model at $z=1$ as [55,
57, 34]

$\beta\frac{\Delta E}{V}=\frac{\beta
e^{2}}{\varepsilon}\frac{\rho}{\sigma}\left[\left(\sum_{\alpha=0}^{B}X_{+}^{\alpha}-1\right)-\left(\sum_{\alpha=0}^{B}X_{-}^{\alpha}+1\right)\right].$
(2.35)

It is worth noting that equations (2.30) and (2.35) include $X_{i}^{\alpha}$,
which, in turn, depend on the parameter $\Gamma$.

The free energy of the ionic subsystem $\beta\Delta f$ is presented in the
form of the sum of two terms, namely, the contribution due to bonding or the
so-called mass action law (MAL) [15, 17] and the contribution due to
electrostatic interactions:

$\beta\Delta f=\beta\Delta f^{(\text{MAL})}+\beta\Delta f^{(\text{el})},$
(2.36)

where

$\beta\Delta f^{(\text{MAL})}=\rho\left(\ln{x}-{1\over 2}x+{1\over
2}\right)-\rho\sum_{i=2}^{m-1}\ln{\left[g_{i,i+1}^{(\text{hs})}(\sigma_{i,i+1})\right]}.$
(2.37)

Here, $x$ is the fraction of free anions (or cations).
$g_{i,i+1}^{(\text{hs})}(\sigma_{i,i+1})$ is the contact values of the radial
distribution functions between monomers in the chain cation [see
(2.33)–(2.34)]. For a model with a cation having the form of a hetero-nuclear
dimer, the second term in (2.37) is given by (2.33).

According to [17] and in the spirit of Wertheim’s multidensity thermodynamic
perturbation theory [46, 48] we assume different approximations for different
terms in the expression for Helmholtz free energy (2.36), i.e., we calculate
$\Delta f^{(\text{el})}$ using the parameter $\Gamma$ obtained in the complete
dissociation limit ($K^{(0)}_{\text{as}}=0$ and $t=0$). In this approach, the
effects connected with ion screening are taken into account exactly, while the
effects of ion association are neglected. As a result, $\beta\Delta
f^{(\text{el})}$ is presented as follows:

$\beta\Delta f^{(\text{el})}=\frac{\beta\Delta
E^{(0)}}{V}+\frac{(\Gamma^{(0)})^{3}}{3\piup},$ (2.38)

where, $\Gamma^{(0)}=\Gamma|_{K_{\text{as}}^{(0)}=0}$ and $\Delta
E^{(0)}=\Delta E|_{K_{\text{as}}^{(0)}=0}$. Accordingly, the equation for
$\Gamma^{(0)}$ is obtained from (2.21) if one put $t=0$. It should be noted
that $\Gamma^{(0)}$ contains a contribution related to the presence of the
neutral beads of the cation chain, and cannot be reduced to the screening
parameter $\Gamma$ in the MSA theory.

For $\beta\Delta E^{(0)}$ we have:

$\beta\frac{\Delta E^{(0)}}{V}=\frac{\beta
e^{2}}{\varepsilon}\frac{\rho}{\sigma}\left(\bar{X}_{+}^{0}+\bar{X}_{+}^{B}-\bar{X}_{-}^{0}-2\right),$
(2.39)

where $\bar{X}_{i}^{\alpha}={X}_{i}^{\alpha}|_{\Gamma=\Gamma^{(0)}}$.

The pressure and chemical potentials of the ionic subsystem can be obtained
using standard thermodynamic relations, namely:

$P=-\frac{\partial F}{\partial V},\qquad\rho(\beta\mu_{-}+\beta\mu_{+})=\beta
f+\beta P,$ (2.40)

where $F=fV$.

#### Thermodynamic functions in the complete association limit.

Consider the complete association limit [17, 12, 13] when all anions and
cations are dimerized. In this limit $K_{\text{as}}^{(0)}\to\infty$ and we get
[see (2.22)–(2.23)]

$t=\frac{1}{2\rho\sigma}.$ (2.41)

Taking into account this assumption, we obtain [34]:

$\displaystyle\beta\Delta f^{(\text{MAL})}|_{K_{\text{as}}^{(0)}\to\infty}$
$\displaystyle=$
$\displaystyle-\rho\left\\{\ln\rho-1+\ln\left[g_{-+}^{(\text{hs})}(\sigma_{-+})\right]+\ln\left[g_{+n}^{(\text{hs})}(\sigma_{+n})\right]\right.$
(2.42) $\displaystyle+$
$\displaystyle\left.\sum_{i=4}^{m-1}\ln{\left[g_{i,i+1}^{(\text{hs})}(\sigma_{i,i+1})\right]}+\beta
U_{-+}^{(C)}(\sigma_{-+})+G_{00}^{(\infty)}(\sigma_{-+})\right\\},$

where

$G_{00}^{(\infty)}(\sigma_{-+})=G_{00}(\sigma_{-+})|_{K_{\text{as}}^{(0)}\to\infty}=\frac{1}{T^{*}}\frac{\left[1-(\eta^{B(\infty)}\sigma_{-+})^{2}\right]}{(1+\sigma_{-+}\Gamma^{(\infty)})^{2}},$
(2.43)
$\Gamma^{(\infty)}=\Gamma|_{K_{\text{as}}\to\infty},\qquad\eta^{B(\infty)}=\Gamma|_{K_{\text{as}}\to\infty}.$

Here, $\Gamma^{(\infty)}$ is obtained from (2.21) under the condition (2.41).
In the same way, we obtain $\eta^{B(\infty)}$.

As the result, the pressure of the original matrix-IL system in the complete
association limit can be presented as follows:

$\beta P=\rho+\beta P^{(\text{ref})}+\beta\Delta P^{(\text{MAL})}+\beta\Delta
P^{(\text{el})},$ (2.44)

where $\beta P^{(\text{ref})}$ is given by (2.12)–(2.17). Taking into account
(2.40) and (2.42), we obtain for $\beta\Delta P^{(\text{MAL})}$

$\displaystyle\beta\Delta P^{(\text{MAL})}$ $\displaystyle=$
$\displaystyle-\rho^{2}\left[\frac{\partial\ln[g_{-+}^{(\text{hs})}(\sigma_{-+})]}{\partial\rho}+\frac{\partial\ln[g_{+n}^{(\text{hs})}(\sigma_{+n})]}{\partial\rho}\right.$
(2.45) $\displaystyle+$
$\displaystyle\left.\sum_{i=4}^{m-1}{\partial\ln{\left[g_{i,i+1}^{(\text{hs})}(\sigma_{i,i+1})\right]}\over\partial\rho}+\frac{\partial
G_{00}^{(\infty)}(\sigma_{-+})}{\partial\rho}\right],$

where $g_{-+}^{(\text{hs})}(\sigma_{-+})$,
$g_{+n}^{(\text{hs})}(\sigma_{+n})$,
$g_{i,i+1}^{(\text{hs})}(\sigma_{i,i+1})$, and
$G_{00}^{(\infty)}(\sigma_{-+})$ are given in (2.31)–(2.34), and (2.43),
respectively.

Electrostatic contribution $\beta\Delta P^{(\text{el})}$ is

$\frac{\beta\Delta
P^{(\text{el})}}{\rho}=-\frac{(\Gamma^{(0)})^{3}}{3\piup\rho}-\frac{2\beta
e^{2}}{\piup\varepsilon\rho}(\eta_{0}^{B})^{2}.$ (2.46)

Similar to (2.44), we represent the chemical potential
$\beta\mu=(\beta\mu_{-}+\beta\mu_{+})$ in the form:

$\beta\mu=\ln\rho+\beta\mu^{(\text{ref})}+\beta\Delta\mu^{(\text{MAL})}+\beta\Delta\mu^{(\text{el})},$
(2.47)

where
$\beta\mu^{(\text{ref})}=(\beta\mu_{-}^{\text{ref}}+\beta\mu_{+}^{\text{ref}})$
can be obtained from (2.18)–(2.20). For $\beta\Delta\mu^{(\text{MAL})}$, using
(2.40), we have

$\beta\Delta\mu^{(\text{MAL})}=\frac{\beta\Delta
f^{(\text{MAL})}}{\rho}+\frac{\beta\Delta P^{(\text{MAL})}}{\rho},$ (2.48)

where $\beta\Delta f^{(\text{MAL})}$ and $\beta\Delta P^{(\text{MAL})}$ are
given in (2.42) and (2.45), respectively. The contribution
$\beta\Delta\mu^{(\text{el})}$ has the form:

$\beta\Delta\mu^{(\text{el})}=\beta\frac{\Delta E^{(0)}}{V\rho}-\frac{2\beta
e^{2}}{\piup\varepsilon\rho}(\eta_{0}^{B})^{2}.$ (2.49)

##  3 Results and discussion

Using the theory proposed in the previous section, we have studied the vapour-
liquid phase behaviour of two models of ILs, model A and model B, confined in
a disordered matrix in the approximation of complete ion association [see
equations (2.44)–(2.49)].

The phase diagrams were obtained at subcritical temperatures using the
conditions of two-phase equilibrium

$\displaystyle\mu(\rho^{(\text{v})},T)$ $\displaystyle=$
$\displaystyle\mu(\rho^{(\text{l})},T),$ $\displaystyle
P(\rho^{(\text{v})},T)$ $\displaystyle=$ $\displaystyle
P(\rho^{(\text{l})},T),$

where the subscripts “v” and “l” refer to the vapour and the liquid phases,
respectively.

We introduce dimensionless units for the temperature, pressure, cation (anion)
number density, and volume fraction

$\displaystyle T^{*}=\frac{k_{\text{B}}T\varepsilon\sigma}{\rm{e}^{2}},\quad
P^{*}=\frac{P\varepsilon\sigma^{4}}{\rm{e}^{4}},\quad\rho^{*}=\rho\sigma^{3},\quad\eta=\frac{\piup}{6}\rho\sum_{\alpha=-,+,n}\sigma_{\alpha}^{3},$

where $\sigma=\sigma_{-}=\sigma_{+}$ is the diameter of the charged monomer
which is the same for anion and cation, and $\sigma_{n}$ is the diameter of
the uncharged monomer of the cation. The ionic model is characterised by the
parameter $k_{1}=\sigma_{n}/\sigma$ [equation (2.32)].

Figure 2: (Colour online) Vapour-liquid phase diagrams of model A with
$k_{1}=1$ (upper left-hand panel), $k_{1}=1.5$ (upper right-hand panel),
$k_{1}=2$ (lower left-hand panel), $k_{1}=2.5$ (lower right-hand panel),
confined in a porous matrix with different porosity. From the top to the
bottom: $\phi_{0}=1.0$ (red lines and squares), $\phi_{0}=0.95$ (black lines),
$\phi_{0}=0.9$ (blue lines). The diameter of matrix particles
$\sigma_{0}=1.5\sigma$. Here, the lines represent theoretical results, and the
symbols are the results of computer simulations [34]. Notations:
$T^{*}=k_{\text{B}}T\varepsilon\sigma/e^{2}$, $\rho^{*}=\rho\sigma^{3}$,
$k_{1}=\sigma_{n}/\sigma$, and $\phi_{0}=1-\eta_{0}$, where
$\eta_{0}=\piup\rho_{0}\sigma_{0}^{3}/6$.

###  3.1 Model A

First, we present the results for model A (see figure 1). Calculations are
done for three values of the matrix porosity $\phi_{0}$, i.e., $\phi_{0}=1$
(bulk case), $\phi_{0}=0.95$ and $\phi_{0}=0.9$. We fix the diameter of the
matrix particles $\sigma_{0}$ and study the effect of the porous medium on the
phase behavior of model A with different $k_{1}$. In figure 2, we present the
vapour-liquid phase diagrams of model A confined in a disordered matrix with
the diameter of matrix particles $\sigma_{0}=1.5$. Comparing the results
obtained for different $k_{1}$ in the absence of a porous matrix, the case
$\phi_{0}=1$ in figure 2 (red curves), we see that an increase of $k_{1}$
leads to a shift of the phase diagrams to lower temperatures and to lower
number densities which agrees with the results obtained in [34]. For the fixed
porosity $\phi_{0}<1$ (black and blue curves), an increase of $k_{1}$ also has
a similar effect: the coexistence curves shift to lower temperatures and
densities when $k_{1}$ increases, and it can be concluded that the presence of
the porous medium strengthens the effects previously observed for the bulk
case. On the other hand, for the fixed $k_{1}$, a decrease of matrix porosity
$\phi_{0}$ also leads to a shift of the phase diagram to lower temperatures
$T^{*}$ and to lower densities $\rho^{*}$, and the phase coexistence region
becomes narrower. For comparison, figure 2 (upper left-hand panel) shows the
results of computer simulations for $k_{1}=1$ in the absence of a porous
medium. It should be noted that the agreement of our theoretical results for
the IL model in the absence of a porous matrix with the results of computer
simulations is within the same order of accuracy as the theoretical
predictions made earlier for the RPM model (see [34]).

Figure 3: (Colour online) Vapour-liquid critical temperature $T_{c}^{*}$ of
model A confined in the porous medium as a function of the porosity
$\phi_{0}$. From the top to the bottom: $k_{1}=1$, $k_{1}=1.5$, $k_{1}=2$,
$k_{1}=2.5$. Notations: $T_{c}^{*}=k_{\text{B}}T_{c}\varepsilon\sigma/e^{2}$,
$k_{1}=\sigma_{n}/\sigma$ and $\phi_{0}=1-\eta_{0}$, where
$\eta_{0}=\piup\rho_{0}\sigma_{0}^{3}/6$.

Figure 4: (Colour online) Critical density $\rho_{c}^{*}$ (left-hand panel)
and critical volume fraction $\eta_{c}$ (right-hand panel) of model A confined
in the porous medium as a function of the porosity $\phi_{0}$ at different
values of $k_{1}$ (as indicated on the plots). Notations:
$\rho_{c}^{*}=\rho_{c}\sigma^{3}$,
$\eta_{c}=\piup\rho_{c}\sum_{\alpha=-,+,n}\sigma_{\alpha}^{3}/6$,
$k_{1}=\sigma_{n}/\sigma$, and $\phi_{0}=1-\eta_{0}$, where
$\eta_{0}=\piup\rho_{0}\sigma_{0}^{3}/6$.  Figure 5: (Colour online) Vapour-
liquid phase diagrams of model A with $k_{1}=2.5$, confined in a porous matrix
with a matrix porosity $\phi_{0}=1.0$ (red curves ), $\phi_{0}=0.95$ (black
curves), $\phi_{0}=0.9$ (blue curves). The diameter of matrix particles
$\sigma_{0}=2.5\sigma$. Notations:
$T^{*}=k_{\text{B}}T\varepsilon\sigma/e^{2}$, $\rho^{*}=\rho\sigma^{3}$,
$k_{1}=\sigma_{n}/\sigma$ and $\phi_{0}=1-\eta_{0}$, where
$\eta_{0}=\piup\rho_{0}\sigma_{0}^{3}/6$.

The dependence of the critical temperature $T_{c}^{*}$ on the matrix porosity
$\phi_{0}$ for $k_{1}=1$, $1.5$, $2$ and $2.5$ is shown in figure 3. It is
seen that the critical temperature $T_{c}^{*}$ decreases almost linearly with
a decreasing porosity for all considered values of $k_{1}$, although a small
deviation from linearity can be observed for $k_{1}=2.5$. Moreover, the
critical temperature $T_{c}^{*}$ decreases with an increase of $k_{1}$.

In figure 4, we show the dependence of the critical density $\rho_{c}^{*}$ and
the critical volume fraction $\eta_{c}$
($\eta_{c}=\piup\rho_{c}\sum_{\alpha=-,+,n}\sigma_{\alpha}^{3}/6$) on the
matrix porosity $\phi_{0}$. Again, the critical density $\rho_{c}^{*}$
decreases almost linearly with a decreasing matrix porosity for all values of
$k_{1}$ except for $k_{1}=2.5$, for which some deviation from a linear
behavior is observed: a decrease of $\rho_{c}^{*}$ becomes steeper for
$\phi_{0}\leqslant 0.925$. Similar to $T_{c}^{*}$, the critical density
$\rho_{c}^{*}$ decreases when $k_{1}$ increases. It should be noted that a
similar behavior of the critical density with a decreasing porosity was
obtained for a size-asymmetric primitive model of spherical ions if the
difference between the ion diameters increases [44]. Regarding the critical
volume fraction $\eta_{c}$, its behavior with a decrease of matrix porosity is
similar to the behavior of $\rho_{c}^{*}$, except that the deviation from the
linear dependence of $\eta_{c}$ is more pronounced for $\phi_{0}\leqslant
0.925$. On the other hand, the dependence of $\eta_{c}$ on $k_{1}$ is opposite
to that obtained for $\rho_{c}^{*}$, i.e., when $k_{1}$ increases $\eta_{c}$
also increases.

In figure 5, we show the vapour-liquid phase diagrams obtained for model A
with $k_{1}=2.5$ confined in a porous matrix with the matrix particles of the
diameter $\sigma_{0}=2.5\sigma$. Comparing these phase diagrams with the
corresponding phase diagrams in figure 2 (lower right-hand panel) for
$\sigma_{0}=1.5\sigma$, one can see that an increase of the diameter of the
matrix particles for the fixed matrix porosity $\phi_{0}$ leads to a shift of
the phase diagrams to higher temperatures $T^{*}$ and to higher densities
$\rho^{*}$, as well as to the broadening of the range of phase coexistence.

Figure 6: Critical temperature $T_{c}^{*}$ (left-hand panel) and the critical
density $\rho_{c}^{*}$ (right-hand panel) of model A with $k_{1}=2.5$ confined
in a porous matrix as a function of matrix porosity $\phi_{0}$ for two sizes
of matrix particles: $\sigma_{0}=1.5\sigma$ (solid curve) and
$\sigma_{0}=2.5\sigma$ (dashed curve). Notations:
$T_{c}^{*}=k_{\text{B}}T_{c}\varepsilon\sigma/e^{2}$,
$k_{1}=\sigma_{n}/\sigma$ and $\phi_{0}=1-\eta_{0}$, where
$\eta_{0}=\piup\rho_{0}\sigma_{0}^{3}/6$.

The dependence of the critical temperature $T_{c}^{*}$ and the critical number
density $\rho_{c}^{*}$ for model A with $k_{1}=2.5$ on the matrix porosity
$\phi_{0}$ for two diameters of matrix particles $\sigma_{0}=1.5\sigma$ (solid
curve) and $\sigma_{0}=2.5\sigma$ (dashed curve) is shown in figure 6. It is
seen that an increase of the diameter of matrix obstacles does not generally
change the behavior of $T_{c}^{*}$ and $\rho_{c}^{*}$: both $T_{c}^{*}$ and
$\rho_{c}^{*}$ decrease with a decrease of the matrix porosity, as it was
observed in the case of $\sigma_{0}=1.5\sigma$. However, the decrease of both
critical parameters is much slower for $\sigma_{0}=2.5\sigma$ and the
difference between $T_{c}^{*}$ ($\rho_{c}^{*}$) obtained for
$\sigma_{0}=1.5\sigma$ and $\sigma_{0}=2.5\sigma$ increases with a decreasing
porosity. Such a behaviour of the critical parameters with the change of the
size of matrix particles is consistent with the results obtained earlier for
the RPM model [41, 22] as well as with the results for a size-asymmetric
primitive models [44, 45].

Figure 7: (Colour online) Vapour-liquid phase diagrams of model B with cation
chain length: $m_{+}=3$ (upper panel), $m_{+}=5$ (lower left-hand panel),
$m_{+}=8$ (lower right-hand panel), confined in a porous matrix with porosity
$\phi_{0}=1.0$ (red lines and squares), $\phi_{0}=0.925$ (black lines),
$\phi_{0}=0.85$ (blue lines). Here, the lines represent theoretical results,
and the symbols are the results of computer simulation [34]. Notations:
$T^{*}=k_{\text{B}}T\varepsilon\sigma/e^{2}$ and $\rho^{*}=\rho\sigma^{3}$.

###  3.2 Model B

Now we focus on model B. In this model, the cations are presented as
mononuclear chains with $m_{+}$ beads (figure 1). The vapour-liquid phase
diagrams are calculated for four values of the length of the cation chain
$m_{+}=m-1=3,5,8$ and for three values of the matrix porosity
$\phi_{0}=1.0,0.925,0.85$. In our calculations, we fix the diameter of matrix
particles $\sigma_{0}=1.5$. The phase diagrams are shown in figure 7. In the
upper panel of these figures, in addition to the theoretical results, we also
show the simulation results for $m_{+}=3$ and $\phi_{0}=1$ (bulk case) [34].
In this case, our theoretical results for the vapour-liquid phase diagram are
in a good agreement with the results obtained using the Monte Carlo
simulations. Moreover, the agreement for the cation with $m_{+}=3$ is much
better than in the case of mononuclear dimer corresponding to $m_{+}=2$ (see
figure 2, upper left-hand panel). We can expect that in the presence of a
porous medium and longer cation chains, our results are satisfactory. As for
model A, the presence of a porous medium shifts the phase diagrams of model B
to lower temperatures and to lower number densities when the matrix porosity
decreases (figure 7). Both the critical temperature and the critical number
density decrease almost linearly when the matrix porosity decreases (figure
8). An increase of the cation chain length in the absence of a porous medium
($\phi_{0}=1.0$) has a similar effect. The presence of the porous matrix
significantly enhances this effect.

It should be noted that there is a small difference between the present
results for the bulk case and the results obtained in [34]. Comparing the
results, one can see that the critical temperatures for $\phi_{0}=1$ in figure
2 (upper left-hand panel) and figure 7 (upper panel) are about $3$% higher
than the corresponding critical temperatures obtained for the same models in
[34]. The difference appears due to different approximations used to calculate
the contact values, i.e., in this work the expressions for the contact values
of the radial distribution functions [equations (2.31)–(2.34)] include the
term corresponding to the Boublik-Mansoori-Carnahan-Starling-Leland
approximation while in [34] the Percus-Yevick contact values are used. Here,
contrary to [34], we use the same approximation for the thermodynamic
functions of the RS and for the contact values.

Figure 8: Dependence of the critical temperature (left-hand panel) and the
critical density (right-hand panel) on the matrix porosity $\phi_{0}$ for
model B with different lengths of the cation chains $m_{+}=2,3,5,8$ confined
in a porous medium. Notations are the same as in figure 7.

##  4 Conclusions

We have studied the vapour-liquid phase behaviour of ILs with complex
molecular structures of ions (cations) confined in a disordered matrix formed
by neutral HSs. The IL is modelled as an electroneutral mixture of HS anions
and flexible linear chain cations, formed by the tangentially bonded HSs with
the charge located on one of the terminal beads [34]. To this end, we have
developed the theoretical approach which combines a generalization of the
scaled particle theory (SPT), Wertheim’s thermodynamic perturbation theory,
and the associative mean-spherical approximation (AMSA). Using this theory, we
obtained analytical expressions for the Helmholtz free energy, pressure, and
chemical potentials of ions in the complete association limit. The
thermodynamic functions contain contributions from the RS, represented by a
multicomponent HS fluid of different sizes, confined in the HS matrix, and the
ionic subsystem, represented by spherical anions and chain cations. It should
be noted that the contribution of the ionic subsystem to the thermodynamics of
the full IL-matrix system depends on the contact values of the radial
distribution functions of the RS, particularly, on the contact values of the
radial distribution functions between the HSs belonging to the anion and the
cation and on the contact values of the radial distribution function between
the HSs belonging only to the cation.

We have calculated the vapour-liquid phase diagrams of two versions of the
above described IL model, i.e., models with cations represented by dimers with
charged and neutral beads of different size and by chains with neutral beads
of the same size, respectively. Our results show that a decrease in the
geometric porosity leads to a shift of the phase diagrams to lower
temperatures and to lower values of the number density of anions (cations), as
well as to lower values of volume fractions of IL for all considered models of
the cation, while the region of the coexistence of the two phases narrows.
This agrees with the results obtained for both simple (nonionic) liquids and
model ILs with spherical ions. On the other hand, the results for IL models
with molecular cations represented by dimers with the neutral bead of
different size and by chains with the neutral beads of the same size in the
bulk showed [34] that an increase of the diameter of the neutral bead or chain
length also leads to a shift of the vapour-liquid phase diagram to lower
temperatures and densities. Therefore, the presence of a disordered porous
medium enhances these effects on the vapour-liquid phase diagram of ILs.
Accordingly, the critical temperature and critical number density decrease
when the matrix porosity decreases and the diameter/length of the neutral part
of the cation increases. An increase in the size of the matrix particles at a
fixed matrix porosity leads to the opposite effect, i.e., the critical
temperature $T_{c}^{*}$ and the critical number density $\rho_{c}^{*}$
increase with an increase of the diameter of the matrix particles and,
accordingly, the phase diagram shifts to higher values of temperature and
density.

It should be noted that in this paper the thermodynamic functions of the
reference system and the contact values of the radial distribution functions
are treated within the framework of the same approximation. This leads, in
general, to a bit better agreement between the theoretical and simulation
results in the bulk case compared to [34]. We hope that this is also valid in
the presence of a disordered matrix.

## Acknowledgements

We thank the Ministry of Education and Science of Ukraine for its financial
support (Agreement No PH/16-2023). YVK acknowledge financial support through
the MSCA4Ukraine project, which is funded by the European Union.

## References

  * [1] Buzzeo M. C., Evans R. G., Compton R. G., ChemPhysChem, 2004, 5, No. 8, 1106–1120, doi:10.1002/cphc.200301017.
  * [2] Silvester D. S., Compton R. G., Z. Phys. Chem., 2006, 220, No. 10, 1247–1274, doi:10.1524/zpch.2006.220.10.1247.
  * [3] Suo L., Borodin O., Gao T., Olguin M., Ho J., Fan X., Luo C., Wang C., Xu K., Science, 2015, 350, No. 6263, 938–943, doi:10.1126/science.aab1595.
  * [4] Correia D. M., Fernandes L. C., Martins P. M., García-Astrain C., Costa C. M., Reguera J., Lanceros-Méndez S., Adv. Funct. Mater., 2020, 30, No. 24, 1909736, doi:10.1002/adfm.201909736.
  * [5] Egorova K. S., Gordeev E. G., Ananikov V. P., Chem. Rev., 2017, 117, No. 10, 7132–7189, doi:10.1021/acs.chemrev.6b00562.
  * [6] Singh M. P., Singh R. K., Chandra S., Prog. Mater. Sci., 2014, 64, 73–120, doi:10.1016/j.pmatsci.2014.03.001.
  * [7] Kondrat S., Feng G., Bresme F., Urbakh M., Kornyshev A. A., Chem. Rev., 2023, 123, No. 10, 6668–6715, doi:10.1021/acs.chemrev.2c00728.
  * [8] Patsahan O., Mryglod I., In: Order, Disorder and Criticality, Vol. 3, Holovatch Yu. (Ed.), World Scientific, Singapore, 47–92.
  * [9] Fisher M. E., Levin Y., Phys. Rev. Lett., 1993, 71, 3826–3829, doi:10.1103/PhysRevLett.71.3826.
  * [10] Levin Y., Fisher M. E., Physica A, 1996, 225, No. 2, 164–220, doi:10.1016/0378-4371(95)00336-3.
  * [11] Holovko M. F., Kalyuzhnyi Yu. V., Mol. Phys., 1991, 73, No. 5, 1145–1157, doi:10.1080/00268979100101831.
  * [12] Kalyuzhnyi Yu., Mol. Phys., 1998, 94, No. 4, 735–742, doi:10.1080/00268979809482366.
  * [13] Kalyuzhnyi Yu. V., Holovko M. F., Vlachy V., J. Stat. Phys., 2000, 100, No. 1/2, 243–265, doi:10.1023/a:1018699914319.
  * [14] Blum L., Bernard O., J. Stat. Phys., 1995, 79, No. 3–4, 569–583, doi:10.1007/bf02184871.
  * [15] Bernard O., Blum L., J. Chem. Phys., 1996, 104, No. 12, 4746–4754, doi:10.1063/1.471168.
  * [16] Jiang J. W., Blum L., Bernard O., Mol. Phys., 2001, 99, No. 20, 1765–1767, doi:10.1080/00268970110072764.
  * [17] Jiang J., Blum L., Bernard O., Prausnitz J. M., Sandler S. I., J. Chem. Phys., 2002, 116, No. 18, 7977–7982, doi:10.1063/1.1468638.
  * [18] Qin Y., Prausnitz J. M., J. Chem. Phys., 2004, 121, No. 7, 3181–3183, doi:10.1063/1.1770651.
  * [19] Caillol J.-M., Patsahan O., Mryglod I., Condens. Matter Phys., 2005, 8, No. 4, 665–684, doi:10.5488/cmp.8.4.665.
  * [20] Patsahan O., Mryglod I., Patsahan T., J. Phys.: Condens. Matter, 2006, 18, No. 45, 10223–10235, doi:10.1088/0953-8984/18/45/009.
  * [21] Patsahan O. V., Patsahan T. M., Phys. Rev. E, 2010, 81, No. 3, 031110, doi:10.1103/physreve.81.031110.
  * [22] Holovko M., Patsahan T., Patsahan O., J. Mol. Liq., 2017, 228, 215–223, doi:10.1016/j.molliq.2016.10.045.
  * [23] Malvaldi M., Chiappe C., J. Phys.: Condens. Matter, 2007, 20, No. 3, 035108, doi:10.1088/0953-8984/20/03/035108.
  * [24] Spohr H. V., Patey G. N., J. Chem. Phys., 2009, 130, No. 10, 104506, doi:10.1063/1.3078381.
  * [25] Martín-Betancourt M., Romero-Enrique J. M., Rull L. F., J. Phys. Chem. B, 2009, 113, No. 27, 9046–9049, doi:10.1021/jp903709k.
  * [26] Fedorov M., Georgi N., Kornyshev A., Electrochem. Commun., 2010, 12, No. 2, 296–299, doi:10.1016/j.elecom.2009.12.019.
  * [27] Wu J., Jiang T., Jiang D., Jin Z., Henderson D., Soft Matter, 2011, 7, No. 23, 11222, doi:10.1039/c1sm06089a.
  * [28] Ganzenmüller G., Camp P., Condens. Matter Phys., 2011, 14, No. 3, 33602, doi:10.5488/cmp.14.33602.
  * [29] Lindenberg E. K., Patey G. N., J. Chem. Phys., 2014, 140, No. 10, 104504, doi:10.1063/1.4867275.
  * [30] Lindenberg E. K., Patey G. N., J. Chem. Phys., 2015, 143, No. 2, 024508, doi:10.1063/1.4923344.
  * [31] Guzmán O., Ramos Lara J. E., del Río F., J. Phys. Chem. B, 2015, 119, No. 18, 5864–5872, doi:10.1021/jp511571h.
  * [32] Silvestre-Alcantara W., Bhuiyan L., Lamperski S., Kaja M., Henderson D., Condens. Matter Phys., 2016, 19, No. 1, 13603, doi:10.5488/cmp.19.13603.
  * [33] Lu H., Li B., Nordholm S., Woodward C. E., Forsman J., J. Chem. Phys., 2016, 145, No. 23, 234510, doi:10.1063/1.4972214.
  * [34] Kalyuzhnyi Yu. V., Reščič J., Holovko M., Cummings P. T., J. Mol. Liq., 2018, 270, 7–13, doi:10.1016/j.molliq.2018.01.109.
  * [35] Madden W. G., Glandt E. D., J. Stat. Phys., 1988, 51, No. 3-4, 537–558, doi:10.1007/bf01028471.
  * [36] Holovko M., Dong W., J. Phys. Chem. B, 2009, 113, No. 18, 6360–6365, doi:10.1021/jp809706n.
  * [37] Patsahan T., Holovko M., Dong W., J. Chem. Phys., 2011, 134, No. 7, 074503: 1–11, doi:10.1063/1.3532546.
  * [38] Holovko M., Patsahan T., Dong W., Pure Appl. Chem., 2013, 85, No. 1, 115–133, doi:10.1351/pac-con-12-05-06.
  * [39] Holovko M. F., Patsahan T., Dong W., Condens. Matter Phys., 2012, 15, No. 2, 23607: 1–13, doi:10.5488/cmp.15.23607.
  * [40] Holovko M., Patsahan T., Dong W., Condens. Matter Phys., 2017, 20, No. 3, 33602: 1–14, doi:10.5488/cmp.20.33602.
  * [41] Holovko M. F., Patsahan O., Patsahan T., J. Phys.: Condens. Matter, 2016, 28, No. 24, 414003, doi:10.1088/0953-8984/28/41/414003.
  * [42] Chen W., Zhao S. L., Holovko M., Chen X. S., Dong W., J. Phys. Chem. B, 2016, 120, No. 24, 5491–5504, doi:10.1021/acs.jpcb.6b02957.
  * [43] Holovko M. F., Patsahan T. M., Patsahan O. V., J. Mol. Liq., 2017, 235, 53–59, doi:10.1016/j.molliq.2016.11.030.
  * [44] Patsahan O. V., Patsahan T. M., Holovko M. F., Phys. Rev. E, 2018, 97, No. 2, 022109, doi:10.1103/physreve.97.022109.
  * [45] Patsahan O., Patsahan T., Holovko M., J. Mol. Liq., 2018, 270, 97–105, doi:10.1016/j.molliq.2017.12.033.
  * [46] Wertheim M. S., J. Stat. Phys., 1984, 35, No. 1-2, 19–34, doi:10.1007/bf01017362.
  * [47] Wertheim M. S., J. Stat. Phys., 1984, 35, No. 1-2, 35–47, doi:10.1007/bf01017363.
  * [48] Wertheim M. S., J. Stat. Phys., 1986, 42, No. 3-4, 459–476, doi:10.1007/bf01127721.
  * [49] Wertheim M. S., J. Stat. Phys., 1986, 42, No. 3-4, 477–492, doi:10.1007/bf01127722.
  * [50] Reiss H., Frisch H. L., Lebowitz J. L., J. Chem. Phys., 1959, 31, No. 2, 369–380, doi:10.1063/1.1730361.
  * [51] Reiss H., Frisch H. L., Helfand E., Lebowitz J. L., J. Chem. Phys., 1960, 32, No. 1, 119–124, doi:10.1063/1.1700883.
  * [52] Lebowitz J. L., Helfand E., Praestgaard E., J. Chem. Phys., 1965, 43, No. 3, 774–779, doi:10.1063/1.1696842.
  * [53] Boublík T., Mol. Phys., 1974, 27, No. 5, 1415–1427, doi:10.1080/00268977400101191.
  * [54] Holovko M., Shmotolokha V., Patsahan T., In: Physics of Liquid Matter: Modern Problems, Vol. 171, Bulavin L., Lebovka N. (Eds.), Springer International Publishing, Heidelberg, 3–30.
  * [55] Protsykevytch I., Kalyuzhnyi Yu. V., Holovko M. F., Blum L., J. Mol. Liq., 1997, 73-74, 1–20, doi:10.1016/s0167-7322(97)00053-6.
  * [56] Bernard O., Blum L., J. Chem. Phys., 2000, 112, No. 16, 7227–7237, doi:10.1063/1.481287.
  * [57] Kalyuzhnyi Yu. V., Cummings P. T., J. Chem. Phys., 2001, 115, No. 1, 540–551, doi:10.1063/1.1376426.
  * [58] Holovko M., In: Ionic Soft Matter: Modern Trends in Theory and Applications, Vol. 206, Henderson D., Holovko M., Trokhymchuk A. (Eds.), Springer Netherlands, Dordrech, 45–81.
  * [59] Hvozd M., Patsahan O., Patsahan T., Holovko M., J. Mol. Liq., 2022, 346, 117888, doi:10.1016/j.molliq.2021.117888.
  * [60] Hvozd M., Patsahan T., Holovko M., J. Phys. Chem. B, 2018, 122, No. 21, 5534–5546, doi:10.1021/acs.jpcb.7b11834.
  * [61] Kalyuzhnyi Yu. V., Holovko M., Patsahan T., Cummings P. T., J. Phys. Chem. Lett., 2014, 5, 4260–4264, doi:10.1021/jz502135f.
  * [62] Holovko M. F., Korvatska M. Y., Condens. Matter Phys., 2020, 23, No. 2, 23605, doi:10.5488/cmp.23.23605.

Ukrainian 0

Ôàçîâà ïîâåäíêà ãàç-ðäèíà ïðèìòèâíèõ ìîäåëåé îííèõ ðäèí ó íåâïîðÿäêîâàíîìó
ïîðèñòîìó ñåðåäîâèù Ò. Ãâîçäü${}^{\textsf{\footnotesize{\mbox{?}}}}$, Ò.
Ïàöàãàí${}^{\textsf{\footnotesize{\mbox{?},\penalty 1000\mbox{?}}}}$, Þ.
Êàëþæíèé${}^{\textsf{\footnotesize{\mbox{?},\penalty 1000\mbox{?}}}}$, Î,
Ïàöàãàí${}^{\textsf{\footnotesize{\mbox{?}}}}$, Ì.
Ãîëîâêî${}^{\textsf{\footnotesize{\mbox{?}}}}$

${}^{\textsf{\footnotesize 1}}$ Iíñòèòóò ôiçèêè êîíäåíñîâàíèõ ñèñòåì
Íàöiîíàëüíî¿ àêàäåìi¿ íàóê Óêðà¿íè 79011, ì. Ëüâiâ, âóë. Ñâíöiöüêîãî, 1,
Óêðà¿íà ${}^{\textsf{\footnotesize 2}}$ Iíñòèòóò ïðèêëàäíî¿ ìàòåìàòèêè òà
ôóíäàìåíòàëüíèõ íàóê, Íàöiîíàëüíèé óíiâåðñèòåò Ëüâiâñüêà ïîëiòåõíiêà, 79013,
Ëüâiâ, âóë. Ñ. Áàíäåðè, 12, Óêðà¿íà, ${}^{\textsf{\footnotesize 3}}$ Ôàêóëüòåò
õiìi¿ i õiìiíî¿ òåõíîëîãi¿, Óíiâåðñèòåò Ëþáëÿíè, âóë. Âåíà, 113, 1000 Ëþáëÿíà,
Ñëîâåíiÿ