Critical point calculation for binary mixtures of symmetric non-additive hard disks

We have calculated the values of critical packing fractions for the mixtures of symmetric non-additive hard disks. An interesting feature of the model is the fact that the internal energy is zero and the phase transitions are entropically driven. A cluster algorithm for Monte Carlo simulations in a semigrand ensemble was used. The finite size scaling analysis was employed to compute the critical packing fractions for infinite systems with high accuracy for a range of non-additivity parameters wider than in the previous studies.


Introduction
Two-dimensional fluid mixtures are quite common in soft-matter and in biological systems. Important examples are particles adsorbed at interfaces, on surfaces of pores in porous materials and biological membranes. When the adsorbed fluid forms a monolayer, it may be modelled as a two-dimensional system. The phenomenon of adsorption has been intensively studied, and one of the main contributors to the field is Stefan Sokolowski [1][2][3][4][5][6]. An interesting question that is not fully solved yet is the phase separation in a binary mixture on surfaces with different curvatures. This question is important not only for the adsorption on curved surfaces present in porous materials, but also for the properties of biological membranes surrounding organelle. In living organisms, the membranes are close to the critical point of the demixing transition [7,8]. Therefore, the phase behavior of multicomponent two dimensional fluids and, particularly, the critical behavior may be of biological importance. During the phase transition, a small change of thermodynamic parameters causes a large change of the composition. It may result in significant shape transformations of the biological membranes since their shape is linked with their composition [9]. Especially interesting are the membranes which form triply periodic bilayers [10,11], as well as porous materials whose internal surfaces are similar to minimal surfaces (the average curvature vanishes at every point, and the Gaussian curvature is negative).
The phase separation in binary fluid mixtures belongs to the Ising universality class, and the universal properties of the two-dimensional Ising model are well known from exact results [12]. The nonuniversal properties, however, should be determined for each particular system, and the only feasible method for a surface with arbitrary curvature is by computer simulations. Thus, it is important to develop a simulation procedure that is fast, efficient and accurate. Moreover, it is important to very accurately determine the critical parameters of a generic model on a flat surface, so that the results of approximate theories or the results of simulations on curved surfaces could be compared with reliable data. The lattice model is not appropriate for investigations of the properties of curved surfaces. Therefore, in this work we have chosen a generic continuous model of non-additive hard core mixtures. Some real phenomena, which can be modelled by a mixture of non-additive hard disks, are discussed, for example, in reference [13] and in references cited herein. The behavior of the hard core mixtures has been studied in bulk [14][15][16][17][18][19][20][21] and in restricted geometry [3,[22][23][24] in three dimensional systems. Much less attention was paid to the two dimensional systems [13,20,[24][25][26][27].
We study the mixture of symmetric non-additive hard disks with the interaction potential defined by: where r is the distance between the centers of two disks, indexes α ∈ {A, B} and γ ∈ {A, B} describe the species. The length scale is set by the A-component hard-disks diameter, σ A A = 1. For symmetric nonadditive mixtures, where ∆ is the non-additivity parameter. We study the mixtures of positive non-additivity with different values of the parameter ∆. This potential is an idealization of interactions in a mixture of identical colloid particles with surfaces covered with polymeric brushes of two types, A and B. The polymeric brushes of different type effectively repel each other, but the polymers of the same type can interpenetrate. For this reason, the separation between the like particles can be shorter than between different particles. Quite evidently, this is an athermal model, because all the allowed configurations are of the same energy. Nevertheless, such mixtures are capable of separating into two phases, i.e., one phase rich in component A and the other one rich in component B. The phase separation is not induced by competition between the internal energy and the entropy as in standard systems, but rather by competition between the entropy of mixing and the entropy associated with the area available for the particles. An increase of the packing fraction defined as where S is the surface area of the system and N A and N B are the numbers of the A and B particles, leads to a larger decrease of the available area in a homogeneous mixture than in a phase-separated system. This effect starts to dominate over the decrease of the entropy of mixing in a phase-separated system for η > η c . Thus, η c plays the role analogous to the critical temperature T c , and (η c − η)/η c plays a role analogous to (T − T c )/T c in standard mixtures with interacting particles. Simulation results [19] show that this model system belongs to the Ising universality class. The universal properties of the model are known from the exact solution of the Ising model in two dimensions, but the nonuniversal properties, such as η c , should be obtained by simulations.
At the absence of an external field, the symmetry of interactions implies, for the two coexisting phases I and II, the following relations: and where µ I α and x I α are the chemical potential and the composition of the component α in the i -th phase, respectively [19,28]. Here, the difference of the chemical potentials h = µ A − µ B plays the role of an external field.
We are interested in the critical properties of a mixture when the external field is zero. Along the symmetry line h = 0, the composition of the coexisting phases is symmetric, and therefore the critical point is at the concentration

The method
We model an open system in contact with a reservoir by using the Semigrand [19,29] Monte Carlo [30,31] simulation method. Using this method, the system is simulated under a constant total number of particles N , total volume V , temperature T , and the difference of the chemical potential of one species with respect to an arbitrarily chosen species ∆µ. Thus, the number of molecules of each species fluctuates, while the total density remains constant. The semigrand ensemble is superior to the grand ensemble in simulating dense fluids because the particle insertion moves are inefficient for dense fluids. For a symmetric binary hard disks mixture, the internal energy and the chemical potential difference are both zero.
The realization of the semigrand ensemble Monte Carlo simulation for symmetric non-additive hard disks requires two kinds of moves: translation and identity change. The identity change moves can be implemented in the following way. A molecule can be chosen randomly from all the molecules, and the identity change move is always accepted if there is no overlap between the particles after the identity change. Such a procedure works quite well for a small number of molecules, up to a few hundred. The simulations near the critical point might be, however, time consuming. Therefore, we use a cluster algorithm to perform the simulations for larger system sizes [16,24]. The idea of the cluster moves is as follows. The system is divided into a set of clusters. The molecules belong to a cluster if the distance from any molecule to any other molecule is smaller than σ A A + ∆. When the clusters are formed, the identity of all molecules in each cluster is randomly changed with the probability p = 0.5. We identify the clusters using the SLINK algorithm [32][33][34], which is fast and does not require a huge amount of computer memory. We have performed the calculations using either local MC moves or cluster moves. The results obtained by both methods were consistent, but the time of calculations was much shorter for the cluster algorithm.
For the translation moves, the maximum displacement is chosen to obtain 50% acceptance ratio. The identity-change and the translation moves are chosen randomly with the ratio of N translations per one cluster identity change move. We have performed calculations with a square or with a rectangular simulation box. The aspect ratio of the rectangular box is taken as 3 : 2 to allow for the arrangement of the molecules into a triangular lattice. Periodic boundary conditions are employed. We have not observed any dependence of the results of calculations on the shape of the simulation box when the simulations were performed for fluid mixtures. The averages are taken over 10 6 Monte Carlo cycles, where a cycle consists of N translation and one cluster identity change move.
In molecular simulations, the results of calculations depend on a system size. The dependence is more pronounced for the calculations near the critical point since the correlation length becomes larger and larger when the system is closer and closer to the critical point. When the correlation length is larger than the size of the computational box, the results of the calculations become biased. That is why it may be very difficult to perform exact calculations of the critical point parameters.
When the system is in the two phase region, we do not always have only one phase in the simulation box during the simulation in the semigrand Monte Carlo method. It is possible that in the simulation box we will have either the first or the second phase. With some frequency, the first phase disappears and the second phase appears and vice versa. The higher is the frequency of this change, the closer the system is to the critical point. That is why it is hard to calculate the concentration at the coexistence as an ensemble average. One may try to overcome this problem by taking the most probable value of the concentration instead of the ensemble average, but near the critical point, this procedure may be problematic especially for two dimensional fluids. The shape of the coexistence curve for the two dimensional fluid is much flatter than for the three dimensional fluids. The critical exponent for the two dimensional fluids is β = 1/8 while for the three dimensional fluids, it is β ≈ 0.3258. The distribution of the concentration is not sharply peaked and it is difficult to obtain the most probable value of the concentration with a sufficiently high accuracy.
Fortunately, it is possible to use the calculations performed for the systems of finite size to get the data on the infinite systems by using the concept of finite-size scaling [19,[35][36][37]. For each value of the non-additivity parameter ∆, the critical packing fraction of the infinite system, η c (∞), can be calculated from the set of apparent critical packing fractions in systems with N particles, η * c (N ), according to the relation [19,36,37]: where the critical exponent ν is equal to 1 for the 2D Ising universality class, and d = 2 for the two dimen- and x c is the critical concentration which is exactly 1/2 for symmetric mixtures. The distribution P N (m) is rescaled in such a way as to have a unit norm and a unit variance, and the rescaled distribution is denoted by P * N (y), where the rescaled order parameter is y = a −1 m N β/(d ν) m with a m denoting a proportionality constant. For the apparent critical packing fraction η * c (N ), P * N (y) has an universal shape P * (y) [19,37]. Since this model belongs to the Ising universality class, the universal function P * (y) is known [38,39].

Results
The order-parameter distribution P N (m) is calculated in the simulations as a histogram, with the number of bins equal to the number of particles N . The function P N (m) obtained in simulations for fixed N and various η is then rescaled in such a way as to have a unit norm and a unit variance. The rescaled function P * N (y) is compared with P * (y) for several values of η, and the best fit gives us η * c (N ). Figure 1 shows the order-parameter distribution function P * N (y) calculated in the simulations for the system with ∆ = 0.2 and N = 324 disks, compared with the distribution function for the 2D Ising model [39]. P * N (y) agrees very well with the universal distribution P * (y) for the value of η that we identify with the apparent critical volume fraction η * c (N ). In the same way, we obtain apparent critical packing fractions for a set of systems with a different size. In figure 2, we present the plot of apparent critical densities for different system sizes calculated for the non-additivity parameter ∆ = 0.2. The apparent critical packing fractions are estimated for a set of finite systems with N = {18, 20, 22, 24, 26, 28, 30}. The critical packing fraction as a function of N −1/2 for an infinite system was obtained by fitting the set of apparent critical packing fractions to a straight line and extrapolating the value of the critical packing fraction at infinity. The same procedure was used to determine the critical packing fractions for all the values of the non-additivity parameter ∆. In practice, the shape of the rescaled distribution function P * N (y) is compared with the universal distribution of the order parameter, P * (y), for the first estimation of the apparent critical packing fraction.
In order to determine the precise value of the apparent critical packing fraction, the fourth order cumulants, The n-th moment 〈m n 〉 can be easily calculated from the distribution of the order parameter P N (m): The apparent critical packing fractions η * c (N ) satisfy the equation U N (η * c ) = U * , and have been read off from the interpolated line for the value of U * = 0.61069 [41]. In figure 3, we present the cumulants, U N (m), as functions of the packing fraction for different system sizes N . The curves cross at the values of the cumulant higher than the universal value U * . Similar behavior was observed in references [20,24].
In figure 4 we present the results of our calculation of the critical packing fractions for a set of the nonadditivity parameter ∆, compared with the results already reported in the literature. The calculations reported in reference [26] significantly overestimate and in reference [27] significantly underestimate the results of other works [13,20,24,25]. This is most probably the result of inappropriate simulation methods employed in those calculations. It can be easily noticed that the results of our calculations are in very good agreement with the results of the calculations reported in reference [20] and reference [24]. In all these works, the authors use cluster algorithms in the simulations and the critical packing fractions are calculated for infinite systems where different kinds of finite size scaling analysis were used. show the calculations from this work, the triangles up -reference [26], the stars - [27], the triangles down - [13], the squares - [24], the diamonds - [20], the triangle right-hand - [25]. The dashed line is just to guide the eye.
In reference [25], the critical point packing fraction was obtained using the method of crossing the reduced second moment. This method is not as precise as the method of finite size scaling analysis. In reference [13], the critical packing fraction was estimated from the simulation of finite systems of the order of N = 2000 molecules. Thus, the values of the critical packing fractions should be lower than the values of the critical packing fractions obtained for infinite systems.
In our work, we employ the same cluster algorithm as in reference [24]. We use, however, a different numerical algorithm to identify the clusters, and different scaling procedure to obtain the critical packing fraction. In reference [20], a different cluster algorithm is used, where the clusters are built by superimposing two configurations and identifying the clusters by selecting the groups of overlapping molecules in the sense of hard core interactions. This method works very well for sufficiently large values of the non-additivity parameter ∆ but is not as good for small ∆. In reference [20], calculations were performed for ∆ ∈ {0.5, 1.0, 2.0, 4.0}. The cluster algorithm used in reference [24] and in our calculations works very well for any value of the non-additivity parameters. In reference [24], calculations were performed for ∆ ∈ {0.1, 0.2, 0.5, 1.0}. In this work, we have performed calculations for all the values of the non-additivity parameter ∆ for which the simulation results already existed, and we have extended the calculation to additional large and small non-additivities. We have performed the calculation for  figure 4, the results of our calculations are indistinguishable from the results of calculations in references [20] and [24]. An increase of the non-additivity parameter ∆ results in a decrease of the critical packing fraction as shown in figure 4. In the mixture of non-additive hard disks, we have a competition between the entropy of mixing and the excluded volume effects. When the non-additivity parameter is large, the excluded volume effects are strong and the

13002-6
demixing process takes place at a lower packing fraction. When the non-additivity parameter is small, the entropy of mixing dominates and the demixing takes place at a higher packing fraction.

Summary and conclusions
We have calculated the critical packing fractions with high accuracy for the non-additive hard disks mixtures for a wide range of the non-additivity parameter. We have used Monte Carlo simulations with cluster algorithm which resulted in rejection-free moves. We have used finite size scaling analysis based on the universal distribution of the order parameter to determine the critical packing fraction for infinite systems. The results of our calculations agree very well with the results of the previous calculations where the finite size scaling analysis was used. The proposed simulation method allows for accurate and unambiguous determination of the critical packing fraction for any value of the non-additivity parameter ∆. We hope that the results of our calculations will be the reference point for testing the results of approximate theories for hard disks systems. They should also allow one to compare the critical properties of the particles adsorbed at flat or at curved surfaces, and to determine in this way the effect of the curvature of the underlying surface on the critical properties of the adsorbed fluid mixture.