What is liquid in random porous media: the Barker-Henderson perturbation theory

We apply the Barker-Henderson (BH) perturbation theory to the study of a Lennard-Jones fluid confined in a random porous matrix formed by hard sphere particles. In order to describe the reference system needed in this perturbation scheme, the extension of the scaled particle theory (SPT) is used. The recent progress in the development of SPT approach for a hard sphere fluid in a hard sphere matrix allows us to obtain very accurate results for thermodynamic properties in such a system. Hence, we combine the BH perturbation theory with the SPT approach to derive expressions for the chemical potential and the pressure of a confined fluid. Using the obtained expressions, the liquid-vapour phase diagrams of a LJ fluid in HS matrix are built from the phase equilibrium conditions. Therefore, the effect of matrix porosity and a size of matrix particles is considered. It is shown that a decrease of matrix porosity lowers both the critical temperature and the critical density, while the phase diagram becomes narrower. An increase of a size of matrix particles leads to an increase of the critical temperature. From the comparison it is observed that the results obtained from the theory are in agreement with computer simulations. The approach proposed in the present study can be extended to the case of anisotropic fluid particles in HS matrices.


Introduction
The original understanding of the nature of the liquid state of matter is connected with Van der Waals equation of state formulated nearly 150 years ago [1]. The Van der Waals picture focuses on different roles of the strong short-ranged repulsive and long-ranged attractive intermolecular interaction in forming the equilibrium properties of dense fluids. According to this picture, the harsh repulsive interactions fix the shape and size of molecules and essentially determine a high-density fluid structure. While the contribution of repulsive interaction is entropic, the contribution of attractive interaction is mainly energetic and can be treated as the perturbation. The first theory of liquids based on the Van der Waals idea was proposed by Barker and Henderson (BH) nearly fifty years ago [2][3][4]. Within this theory, the intermolecular potential is separated into the repulsive and attractive parts. The short-ranged repulsive contribution is described within the framework of the hard spheres model while the long-ranged attractive part is included within the Zwanzig high-temperature perturbation theory [5]. A few years later, Andersen, Chandler and Weeks (ACW) developed a somewhat different theory of liquid based on the Van der Waals approach [6][7][8][9]. In this theory, instead of intermolecular interactions, they separated intermolecular forces into repulsive and attractive parts. They also used the optimized cluster expansion (OCE) instead of the high-temperature perturbation theory for the treatment of attractive forces. However, the first term of perturbation related to the high temperature approximation (HTA) in both theories is identical except that the repulsive and attractive interactions are not exactly the same.
In this paper we make use of the Van der Waals ideas and expand BH theory for the description of liquids adsorbed in random porous media. To this end, we use the Madden-Gland model [10]. According to this model, a porous medium is presented as a quenched configuration of randomly distributed hard spheres forming a so-called matrix. A specific description of a fluid in such porous media is connected with double quenched-annealed averages: the annealed average is taken over all fluid configurations and the additional quenched average should be taken over all realizations of a matrix. One of the popular approaches to the solution of this problem is based on the replica method, which allows for the extension of many theoretical methods of liquid state physics to the case of a fluid confined in a random porous medium. For instance, using the replica Ornstein-Zernike (ROZ) integral equation theory [11], the statistical mechanics approach of liquid state was extended to the description of various fluids confined in random porous media [12,13], including the chemically reacting fluids [14,15]. However, unlike bulk fluids, no analytical result has been obtained from the ROZ integral equation approach even for such a simple model as a hard sphere (HS) fluid in a HS matrix. At the same time, the model of HS fluid has a peculiar importance, since similar to the case of a bulk fluid [16,17], it can be used as the reference system in the development of different perturbation schemes.
The first rather accurate analytical results for a HS fluid in a HS matrix were obtained quite recently in [18][19][20][21][22][23] by extending the scaled particle theory (SPT) [24,25] to the case of a HS fluid confined in random porous media. The SPT approach is based on the combination of the exact treatment of a point scaled particle in a HS fluid combined with the thermodynamical treatment of a macroscopic scaled particle. The exact result for a point scaled particle in a HS fluid in a random porous medium was obtained in [18]. However, the approach proposed in [18] referred to as SPT1 contains a subtle inconsistency appearing when the size of matrix particles is much larger than the size of fluid particles. Later on, this inconsistency was eliminated in a new approach referred to as SPT2 [20]. Starting from this formalism, a series of new approximations were developed [20][21][22][23]. Among these approximations we select only SPT2b and SPT2b1, which will be used in this paper for the description of the reference system in the BH perturbation theory.
The liquid-vapour phase diagrams play a central role in understanding the nature of liquids. In contrast to the bulk case, the results of investigations of liquid-vapour phase equilibrium of a simple fluid in random porous media are rather controversial. Computer simulations of a simple fluid confined in a HS matrix [26][27][28][29] demonstrate a possibility of the existence of two phase transitions. One is analogous to the bulk liquid-vapor transition, but with a narrower coexistence curve and lower values of the critical density and critical temperature. The second transition occurs at lower temperatures and at higher densities, and it is interpreted as a phenomenon related to the wetting effects in the fluid located in more confined regions of the matrix. On the other hand, in more thorough investigations [29], it was noticed that this second transition is extremely sensitive to a particular matrix configuration. In the case of HS matrix, it was shown that two phase transitions appear in some realizations of the matrix, while a single phase transition was observed in the others. In order to explain this observation, optimized cluster expansions were used in [30]. It was found that different approximations can lead to qualitatively different results. For example, the mean spherical approximation (MSA) gives only a single liquid-vapour transition, but the inclusion of the second and the third cluster coefficients result in two phase transitions. At the same time, it was shown in [31] that the association theory leads to one phase transition. Our recent investigation, in which we generalized the Van der Waals equation for simple fluids in random porous media, also demonstrates only one liquid-vapour transition [22].
In the present study we combine the BH theory with the previously developed SPT approach to describe a liquid-vapour phase behaviour of a Lennard-Jones fluid confined in a random matrix formed by hard spheres (HS matrix). For a comparison, the results obtained in different approximations are considered as well. As it was mentioned above for the reference system, the SPT approach is applied using the SPT2b approximation and its improved version SPT2b1. It should be noted that the SPT2b approximation is our first really successful result for confined fluids [20], although it may have essential problems at high fluid densities and/or for low matrix porosities. At the same time, the improved SPT2b1 approximation is based on the original SPT2b, but it is free of this shortcoming and gives an accurate description of the thermodynamics for a hard sphere fluid in a hard sphere matrix up to the close packing conditions [21][22][23]32]. Apart from the BH theory, we present the results obtained in the HTA approximation. Using the developed approach, the phase diagrams for a confined Lennard-Jones are built. Different ma-trix porosities and matrix particle sizes are considered. To check the accuracy of the approaches proposed in this paper, a comparison with the MSA results as well as with the results of Monte Carlo simulations is performed.

Reference system: HS fluid in HS matrix
We start our theoretical consideration with the description of the reference system. For this purpose we briefly recapitulate the main ideas of the SPT theory and present here the expressions for the chemical potential and pressure of a HS fluid in a HS matrix, which were obtained in our previous papers [20,21], and which are needed in the current study. The key point of the SPT theory consists in a derivation of the excess chemical potential of an additional scaled particle of a variable size inserted in a fluid. This excess chemical potential is equal to a work needed to create a cavity in a fluid which is free from any other particles. For a small scaled particle in a HS fluid in the presence of a porous medium, the expression for the excess chemical potential is equal to [20]: 1 is the fluid packing fraction, ρ 1 is the fluid density, σ 1 is the diameter of HS fluid particles. The term p 0 (λ s ) = exp(−βµ 0 s ) is defined by the excess chemical potential of the scaled particle confined in an empty matrix, µ 0 s . It has the meaning of probability to find a cavity created by the scaled particle in the matrix in the absence of fluid particles. We should note that here we use conventional notations [11][12][13][14][15][18][19][20][21][22], where the index "1" is used to denote a fluid component, the index "0" denotes matrix particles, while for the scaled particles the index "s" is used. For a large scaled particle, the excess chemical potential is presented by the thermodynamic expression for the work needed to create a macroscopic cavity inside the fluid, which at the same time is confined in a porous medium, and the corresponding expression can be presented as follows: where P is the pressure of fluid, V s is the volume of a scaled particle. The multiplier 1/p 0 (λ s ) appears due to an excluded volume occupied by matrix particles. In this context, it should be mentioned that the probability p 0 (λ) is directly related to two different types of porosity [20][21][22]. The first one corresponds to the case of λ s = 0 and provides the geometrical porosity φ 0 = p 0 (λ s = 0), (2.3) which depends only on the structure of a matrix and it is equal to the volume fraction of a void between the matrix particles. For a HS fluid in a HS matrix, it is equal to where η 0 = 1 6 πσ 3 0 ρ 0 , ρ 0 = N 0 V , N 0 is the number of matrix particles, σ 0 is the diameter of the matrix particles, V is the volume of the system.
Using the SPT theory [24,25] for the case of a HS fluid in a HS matrix, the following expression for φ can be derived:
According to the ansatz of SPT [18-22, 24, 25], w(λ s ) can be presented in the form of an expansion: Coefficients of this expansion can be found from the continuity of µ ex s and the corresponding derivatives ∂µ s /∂λ s and ∂ 2 µ s /∂λ 2 s at λ s = 0. After setting λ s = 1, the expression (2.2) yields the relation between the pressure P and the excess chemical potential µ ex 1 of a fluid: where the coefficients A and B determine the porous medium structure, and for a HS fluid in a HS matrix, they are as follows: (2.10) Using the Gibbs-Duhem equation, which relates the pressure of a fluid with its total chemical potential one derives the fluid compressibility as After dividing the expression (2.12) by ρ 1 and subsequently integrating it over ρ 1 , one obtains the chemical potential: It is worth noting that the second term − ln(φ) in (2.13) follows from the relation (2.5) and the corresponding substitution βµ 0 1 = − ln(φ). Similarly, integration of the right-hand side of expression (2.12) over ρ 1 leads to the pressure (2.14) The expressions (2.13) and (2.14) are considered as the expression derived within the framework of the SPT2 approach [20]. A simple analysis of (2.13) and (2.14) shows that they have two divergences at η 1 = φ and η 1 = φ 0 . Since φ < φ 0 , the divergence at η 1 = φ occurs at lower densities. However, from the geometrical point of view, this divergence should appear at higher densities near the maximum value of the fluid packing fraction available for a fluid in a given matrix. Different corrections improving the SPT2 approach were proposed in [20][21][22]. Here, we consider two of them, which provide rather accurate results in comparison with computer simulations. The first of them known as SPT2b can be derived if φ is replaced by φ 0 everywhere in (2.12) except for the first term. In this case, the chemical potential and the pressure of a confined fluid are as follows: The second approximation referred to as SPT2b1 can be derived from SPT2b by removing the divergence at η 1 = φ by an expansion of the logarithmic term in (2.15) As a consequence, one obtains the following expressions within the SPT2b1 approximation:

BH perturbation theory for simple fluid in random porous medium
The next step of theoretical treatment is connected with a consideration of an attractive part of interaction. We consider a simple fluid with an intermolecular interaction in the form v 11 (r ) = ∞, r < σ 1 , u 11 (r ), r > σ 1 , (2.20) where u 11 (r ) 0 is a pure attractive part of interaction.
In order to take into account the attractive part of interaction, in this subsection we generalize the BH perturbation theory for the case of a fluid in a random HS matrix. To this end, we use the replica trick [11] according to which a system of a fluid in a matrix of unmovable (frozen) particles can be replaced by an equilibrium mixture consisting of the movable (annealed) matrix particles and s identical copies (or replicas) of a fluid. The condition is also set that the fluid replicas from different copies do not interact with each other, but they interact with the matrix. Such a system can be described in a standard way using the liquid state theories, and the properties of a fluid can be obtained by considering the limit s → 0. Therefore, the Helmholtz free energy of a fluid in a matrix can be presented as [30]: where F (s) is the free energy of the (s + 1)-component equilibrium mixture.
Within the framework of the Zwanzig high-temperature perturbation theory [5], the first term of the free energy expansion corresponds to the high-temperature approximation: which is of the same form as in the optimized cluster expansions [30]. F 0 and g HS 11 (r ) are the free energy and the pair distribution function of a HS fluid in a HS matrix, respectively. The second correction term of Zwanzig expansion involves the three-and four-body distribution functions, for which it is difficult to find simple satisfactory approximations. Therefore, instead of this, we follow Barker and Henderson [3,4] and we write the free energy in the form: where we use the same semimacroscopic arguments as in [3,4] and neglect the cross-correlation terms between fluids from different replicas.
Differentiating the expressions (2.22) and (2.23) with respect to the fluid density, one derives the expression for the chemical potential of a fluid: where µ HS 1 is the HS contribution of the reference system, which can be given by equation (2.15) within the SPT2b approximation or by the equation (2.18) within the SPT2b1 approximation. The first term of (2.23) corresponds to the HTA approximation, and it has the following form: The contribution coming from the BH approximation as the second term of (2.23) is as follows: The expressions for the first and second derivatives of the isothermal compressibility with respect to the fluid density for the reference system,

13607-6
within the framework of the SPT2b or SPT2b1 approximations, correspondingly, and they are presented in Appendix. The functions I (ρ 1 ) and J (ρ 1 ) are the integrals from the first and second terms of (2.23) These integrals contain the pair distribution function g HS 11 (r ), which is unknown for this moment and will be considered separately in the next subsection. The pressure can be calculated by differentiating the expression (2.23) with respect to the volume of a system or from the general thermodynamical relation: (2.28) The general form of the pressure is as follows: where P HS is the HS contribution given by (2.16) within the framework of SPT2b approximation or by (2.19) within the framework of SPT2b1 approximation. The HTA term for the pressure is as follows: The contribution of BH term is as follows: (2.31)

The replica Ornstein-Zernike equations
The pair distribution function g HS 11 (r ) needed to calculate the integrals (2.27) can be obtained from a solution of the so-called replica Ornstein-Zernike (ROZ) equations, which were derived by Given and Stell [11] using the replica trick: h 00 = c 00 + ρ 0 c 00 ⊗ h 00 ,  As usual in the liquid state theory [16,17], the ROZ equations need additional closure relations. The Percus-Yevick (PY) approximation is used for a HS fluid in a HS matrix considered in this paper as the reference system. As it was shown in [33,34], this approximation gives results for the pair distribution functions which are in good agreement with computer simulations. In the PY approximation for a HS fluid in a HS matrix, the blocking direct correlation function is zero c b (r ) = 0, thus c c (r ) = c 11 (r ). Therefore, the closure conditions in our case are as follows: where σ 01 = 1 2 (σ 0 + σ 1 ).
The set of equations (2.32) in combination with the closure relations (2.34) is solved numerically using the hybrid Newton-Raphson procedure [35]. Some of the results for g HS 11 (r ) = 1 + h 11 (r ) are presented in figure 1. As one can see, the function g HS 11 (r ) shows a typical behaviour for a HS fluid. In the left-hand panel of figure 1 it is shown that for the fixed fluid density, an increase of the matrix density, i.e., lowering the matrix porosity, leads to a short-range order increase. The same effect is observed if the fluid density is increased, but the matrix density is kept constant (figure 1, right-hand panel). There is also a comparison of the ROZ results with the Monte-Carlo simulations for the low and high fluid densities of a HS fluid in a HS matrix. It is clearly seen that they fit very well, except the contact value g 11 (σ + 1 ) for the dense fluid, which is somewhat lower in the case of ROZ equations. Besides, there are two points of the contact value obtained from the SPT for a fluid in a matrix [32], which are a bit higher than the values of ROZ, thus closer to the corresponding simulation results. In this paper, for a comparison we also consider the mean spherical approximation (MSA) for the studied model of a fluid with the pair interaction between particles in the form (2.20). Similarly to the HTA and BH approximations, to calculate the chemical potential and the pressure within the MSA approximation, the correlation functions are needed. For this purpose, we can also use the ROZ equations in combination with the corresponding closure relations: Again, using the replica procedure [11] and the general expressions of Hoye and Stell [36] for thermodynamic properties of an equilibrium mixture of (s + 1) components in the limit s → 0, one obtains the chemical potential and pressure of a fluid within the framework of the MSA approximation: where "HS" denotes the reference system, i.e., the HS system. For the perturbation part within the MSA approximation, the expressions for the chemical potential is [30]: (2.37) In the same way the pressure of a fluid confined in a matrix can be found: where c HS 11 (r ) and c HS 10 (r ) are the fluid-fluid and fluid-matrix direct correlation functions for a HS fluid in a HS matrix, g 11 (σ + 1 ) and g 10 (σ + 10 ) are the contact values of pair distribution functions g 11 (r ) = 1 + h 11 (r ) and g 10 (r ) = 1 + h 10 (r ).

Some calculation details
In the present study we consider a fluid of spherical molecules confined in a HS matrix. According to (2.20), the fluid-fluid interaction is decomposed into a hard-sphere part and a Lennard-Jones (12-6) tail following Weeks, Chandler and Andersen [6], i.e., u 11 (r ) = − , σ 1 < r < 6 2σ 1 , 4 (σ 1 /r ) 12 − (σ 1 /r ) 6 , r > 6 2σ 1 . (2.39) Applying a hard core potential for the repulsive potential, we avoid difficulties connected with the temperature dependence of the fluid diameter, which requires a special treatment in the case of soft repulsion [3]. The Lennard-Jones attractive tail is truncated at r = 2.5σ 1 in order to compare our results with the corresponding computer simulations [28].
The functions I (ρ 1 ) and J (ρ 1 ) used in the generalized BH theory should be calculated from the integrals (2.27). The dependencies of I (ρ 1 ) and J (ρ 1 ) on the fluid density are illustrated in figure 2. As one

13607-9
can see, these dependencies are rather smooth, that is why they can be interpolated as polynomials. It is noticed that polynomials of the order of 9 provide a satisfactory fitting for the given functions. Since finally the functions I (ρ 1 ) and J (ρ 1 ) are polynomials, all first and second derivatives of these functions used in (2.25)-(2.26) and (2.30)-(2.31) are calculated analytically.
Having the chemical potential and pressure as functions of ρ 1 (or η 1 ) at different temperatures, one can calculate the coexistence curves of the liquid-vapour phase transition. For this purpose, we solve a set of two non-linear equations which follows from the conditions of thermodynamic equilibrium:

Results and discussions
We apply the theory presented in the previous section for the description of liquid-vapour phase coexistence of a simple fluid in a HS random porous medium. However, we first consider the bulk case, i.e., when ρ * 0 = 0. The fluid-fluid interaction v 11 (r ) is taken in the form (2.20) with the attractive potential u 11 (r ) (2.40). In figure 3, one can see the liquid-vapour phase diagrams in coordinates T −ρ 1 for a bulk fluid obtained within different approximations. The simulation results obtained in [28] using the method of grand-canonical Monte Carlo are shown for a comparison. As it is seen, the HTA approximation gives a good description only for low temperature and leads to the overestimation at higher temperatures. The BH approximation slightly improves the diagram, but still essentially overestimates the critical temperature. As expected, the MSA approximation provides the best description of the phase diagram among the considered ones, although the critical temperature is still higher than in the simulations.
Considering a fluid in a matrix within the model applied in our study, one should take into account that there is no attractive interaction between fluid and matrix particles. Thus, the only effect of the matrix is to confine the fluid in the void volume formed between matrix particles. Therefore, the most relevant parameters, which determine fluid properties, are the matrix density ρ * 0 or the corresponding matrix porosity φ 0 , and the size ratio of fluid and matrix particles τ = σ 1 /σ 0 . In figure 4, we present the liquid-vapour phase diagrams obtained for a fluid in matrices of different densities ρ * 0 = 0 (bulk), 0.046, 0.15 and 0.30. In this case, the fluid and matrix particles are of equal sizes, i.e., τ = 1. For this purpose, dashed lines correspond to the HTA approximation combined with SPT2b, solid lines -HTA combined with SPT2b1, symbols -GCMC results taken from [28]. Right-hand panel: dashed lines correspond to the BH theory combined with SPT2b, solid lines -BH theory combined with SPT2b1, dotted lines -MSA approximation, symbols -the GCMC results taken from [28].
we use the theoretical approaches considered in the previous section, i.e., the HTA and BH approximations in combination with the reference system obtained within the SPT2b and SPT2b1 approaches. The computer simulation results taken from [28] are presented in figure 4 for comparison. Similarly to the bulk case, one can see the MSA results for the considered systems [30]. It is observed that for a low matrix density ρ * 0 = 0.046, the coexistence curves calculated using the SPT2b and SPT2b1 approximations almost coincide. With an increase of the matrix density up to ρ * 0 = 0.15, the difference between the results obtained with SPT2b and SPT2b1 becomes more distinguishable. And finally, for the high matrix density ρ * 0 = 0.3, the diagrams differ essentially in these approximations. Moreover, the results obtained with the use of the SPT2b are rather anomalous, since they are far from any other approximation. This anomaly can be explained by the divergence contained in the expressions for the chemical potential and pressure of a fluid when η 1 → φ [21,22], and which become important at high matrix densities. To illustrate this problem, the dependencies of the chemical potential of a HS fluid in a HS matrix of densities ρ * 0 = 0.15 and ρ * 0 = 0.3, in comparison with the results of grand-canonical Monte Carlo are shown in figure 5. It is clearly seen that for the matrix density ρ * 0 = 0.3 the chemical potential of a fluid becomes wrong at the densities ρ * 1 > 0.2 and tends to infinity around 0.33. The latter value corresponds to the η 1 ≈ φ, where the divergence is expected. However, this is not the case for ρ * 0 = 0.15, where the SPT2b approximation is close to the result of SPT2b1 approximation. On the other hand, it is observed that the SPT2b1 approximation perfectly fits the simulation results for the both matrix densities up to the highest values of fluid densities. Therefore, STP2b1 should be considered as the best choice for the description of a HS fluid in HS matrix. Consequently, hereafter we restrict ourselves only to this approximation for the reference system.
All the approximations HTA, BH and MSA correctly reproduce the basic trends of the behaviour of liquid-vapour coexistence curves of a fluid in a matrix, i.e., a decrease of matrix porosity (or an increase of the matrix density) leads to a critical point shift toward lower fluid densities and lower temperatures, simultaneously the phase diagram becomes narrower ( figure 4). Furthermore, all of the considered approximations give only one critical point. A comparison of the diagrams obtained using the approximations presented in our paper with computer simulations ( figure 4) shows that the inclusion of the second term in the BH theory essentially improves the description of coexistence curves. The computer simulations data are in semiquantitative agreement with the theoretical prediction based on the MSA and BH approximation. It is worth mentioning that in contrast to MSA, which is mostly numerical approach, the BH approximation is more an analytical theory. repulsive part of LJ potential is a hard core potential. Thus, according to the form (2.20), we use v 11 (r ) =    ∞, r < σ 1 , 4 (σ 1 /r ) 12 − (σ 1 /r ) 6 , σ 1 < r < 2.5σ 1 , 0, r > 2.5σ 1 . 2) with computer simulations [29] (symbols) for the LJ potential (3.1) at τ = 1 and τ = 3/2, but for the fixed matrix porosity φ 0 = 0.95 (or the matrix packing fraction η 0 = 0.05). As one can see, the form of the phase diagrams obtained from the BH theory is rather close to the data of simulations, although for the both cases of τ they are notably narrower. This deviation is systematic and anticipated, since the simulations are performed for the conventional LJ potential, which has a soft core, while we use the reference system as a hard core fluid. According to (3.1), the repulsive potential is as follows: ϕ 11 (r ) = 4 (σ 1 /r ) 12 − (σ 1 /r ) 6 , r < σ 1 , 0, r > σ 1 .
It is acceptable to substitute this repulsive part of the LJ potential by a hard core, but instead of the diameter of HS particles σ 1 , one should take an effective diameter d 1 which is somewhat smaller than σ 1 and depends on the fluid temperature. To take into account a correct value of the effective diameter d 1 we use one of the successful and the most popular relations proposed by Barker and Henderson for a LJ fluid [3]: Simple calculations of d BH 1 (T ) depending on temperature show that it does not vary too heavily. For instance, for the temperature T * = 1.1, the effective diameter is d BH 1 = 0.9711 and for T * = 0.6, the diameter is d BH 1 = 0.9815. However, it still can have a strong effect on the thermodynamics of the system, hence on the curves of the liquid-vapour coexistence. Therefore, we need to substitute the diameter of HS particles σ 1 by d BH 1 (T ) in every place where it is needed in the expressions used for the reference system, i.e., the terms containing σ 1 or depending on it should be modified. The corrected expressions for the chemical potential and the pressure for a fluid in a HS matrix are used and calculated depending on the temperature. Again using the BH theory in combination with the SPT2b1 approximation, the phase diagrams of liquid-vapour transition are obtained. The results of this correction is presented as dashed lines in figure 7, and, as one can see, the theoretical curves coincide very well with the computer simulation data. It is worth noting that there is a quicker and more efficient way to improve the present phase diagrams. Our preliminary calculations show that coexistence curves change negligibly along the T -axis, and most deviations take place along the ρ 1 -axis. This is mainly related to the effect of an excluded volume, which depends on η 1 = πd 3 1 ρ 1 /6 and is overestimated in the case of the reference system with a hard sphere size d 1 = σ 1 . Therefore, to improve our results for the conventional LJ fluid, the packing fraction should be replaced by η BH 1 = π[d BH 1 (T )] 3 ρ 1 /6, and this is equivalent to the rescaling of the fluid density as: The correction made in such a way allows us to obtain the diagrams which are practically equivalent to those shown in figure 7.
Finally, in figure 8 we present the critical temperature T * cr and the critical density ρ * an increase of τ at the fixed porosity, and the change of the critical density is very small. For example, from the results presented in figure 7 it is found that for φ 0 = 0.95 and τ = 2/3, the critical temperature T * cr = 1.095 and the critical density is ρ * 1,cr = 0.273, while for τ = 1, the critical temperature is T * cr = 1.031 and the critical density is ρ * 1,cr = 0.275. If we estimate this in the limit τ → 0, one can see that the critical temperature T * cr shifts to T * (bulk) cr = 1.209 and the critical density to ρ * 1cr → ρ * (bulk)

Conclusions
In this paper, the Barker-Henderson (BH) perturbation theory is generalized for the a Lennard-Jones fluid confined in a random porous matrix. As the reference system, a hard sphere fluid in a hard sphere matrix is chosen. To describe the reference system, the extension of the scaled particle theory (SPT) is used, and two corresponding approximations are tested. It is shown that the SPT2b1 approximation, which was developed recently, makes it possible to achieve a very accurate description of the thermodynamics of confined hard sphere systems. Combining the SPT approach with the BH theory, the expressions for the chemical potential and the pressure of a simple fluid in a hard sphere matrix are derived. Based on the obtained expressions, the phase diagrams of liquid-vapour transition are calculated and compared with other theoretical approaches such as the high-temperature approximation and the meanspherical approximation. A comparison of our results with computer simulation data found in the literature is made as well. Different matrix porosities as well as the size ratios of fluid and matrix particles are considered in the present paper. The proposed extension of the BH theory for the case of a fluid in a matrix provides a good qualitative agreement with the computer simulations, and in some situations it provides a marvelous quantitative agreement, as it is observed in the case of low matrix porosities and at temperatures which are not very close to the critical point. The theory correctly reproduces the basic effects of porous media on the liquid-vapour phase coexistence of simple fluids, i.e., with a decrease of porosity, the critical point shifts toward lower fluid densities and lower temperatures. It is also observed that for a fixed matrix porosity, but for variable sizes of matrix particles, the critical temperature increases if the size of matrix particles becomes larger and moves to the value of critical temperature of a bulk fluid, while the critical density changes weakly, and in the limit τ → 0, it moves to the bulk critical value normalized by the porosity φ 0 .
The approach developed in this paper can be extended to the case of more complex fluid systems in a confinement. In future, we plan to generalize the BH theory for anisotropic fluids in random porous media. Our very recent investigation [37] connected with the extension of the Van der Waals theory to the case of anisotropic fluids in random porous matrices shows that due to anisotropic interactions, the orientational order causes a competition between isotropic and anisotropic interactions, and the effect of a matrix can essentially modify the liquid-vapour phase diagram.

A. Isothermal compressibility for a HS fluid in HS matrix
Here, we present the expression for the isothermal compressibility where A and B are given in equation (2.10).
The dependence of  figure 9. One can see that for low densities of ρ * 0 , the results in the both approximations are nearly identical, but for ρ * 0 = 0.30, the results within SPT2b show an odd behaviour for η 1 larger than 0.15.