Theory and computer simulation of hard-core Yukawa mixtures: thermodynamical, structural and phase coexistence properties

In all our calculations, we used σ as unit length. We also define $\epsilon_{11} = \epsilon_{22} = \epsilon_{l}$, and $\epsilon_{12} = \epsilon_{u}$ as the depth of the like and unlike interaction, respectively, and introduce the ratio $\alpha = \epsilon_{u}/\epsilon_{l}$ to characterize the relative strength of interaction in the mixture.

Results for the excess internal energy per particle $U/N$ (with N total number of particles), the pressure P, and the chemical potential μ are reported in table 1 for $c_{1}=c_{2}=0.5$, $\alpha=0.75$ and different temperatures, densities and z values. By first looking at the MC data at fixed z, it appears that $U/N$ decreases (i.e. increases in absolute value) with the increasing density. For given z and $T^{*}$, such a decrease is nearly linear in ρ. This is not surprising, given that one has

$$\begin{align} \displaystyle \frac{U}{N}=\frac{\rho}{4}\int {\rm d}{\bf r}\, [v_{11}(r)g_{11}(r)+\alpha v_{11}(r)g_{12}(r)] \,. \label{eint} \nonumber \end{align} \tag{ 6 }$$

Table 1. Internal energy per particle $U/N$, pressure P and chemical potential μ for a symmetric equimolar HCYM with $\alpha=0.75$ at various z values and temperatures. Here $\epsilon = \epsilon_{11}$. Values in parentheses refer to error bars in the simulation data.

z	$T^{*}$	$\rho^{*}$	$U/(N \epsilon)$			$P \sigma^3/\epsilon$			$\mu/\epsilon$
			MC	HRT	MHNC	MC	HRT	MHNC	MC	HRT	MHNC
1.8	1.08	0.45	$-2.680(0.0019)$	$-2.707$	$-2.677$	0.146(0.0030)	0.140	0.133	$-3.429(0.0031)$	$-3.530$	$-3.281$
1.8	1.08	0.55	$-3.278(0.0042)$	$-3.303$	$-3.282$	0.249(0.0003)	0.233	0.218	$-3.469(0.0020)$	$-3.347$	$-3.385$
1.8	1.04	0.55	$-3.438(0.0032)$	$-3.316$	$-3.290$	0.103(0.0045)	0.154	0.140	$-3.851(0.0060)$	$-3.474$	$-3.500$
2.4	0.85	0.40	$-1.817(0.0020)$	$-1.862$	$-1.809$	0.139(0.0004)	0.141	0.135	$-2.402(0.0017)$	$-2.587$	$-2.590$
2.4	0.85	0.50	$-2.248(0.0023)$	$-2.278$	$-2.234$	0.191(0.0012)	0.181	0.177	$-2.508(0.0021)$	$-2.498$	$-2.499$
3.0	0.75	0.40	$-1.491(0.0011)$	$-1.542$	$-1.479$	0.170(0.0010)	0.168	0.166	$-1.873(0.0016)$	$-2.045$	$-2.035$
3.0	0.75	0.50	$-1.846(0.0033)$	$-1.892$	$-1.833$	0.234(0.0019)	0.222	0.223	$-1.929(0.0016)$	$-1.927$	$-1.909$

For the states reported in the table, the change in density is rather small and will therefore have a minor effect on the integral. Indeed, a semi-quantitative estimate of $U/N$ for these states can already be obtained by mean-field theory, in which correlations are disregarded altogether.

Inspection of entries at equal densities in the lower part of table 1 shows that $U/N$ increases (i.e. decreases in absolute value) on increasing z, despite the decrease in $T^*$. Such a behavior is also accounted for by (6), since as z increases, the off-core part of the interaction decreases in absolute value, and this effect dominates over the decrease in $T^*$, given the small amount of the latter. Again, mean-field theory provides a semi-quantitative estimate of the change of $U/N$.

As for theoretical predictions, it emerges from table 1 that the MHNC energies are almost quantitative in comparison to the MC data for the different set of parameters explored, namely in the ranges $0.4 \leqslant \rho^* \leqslant 0.55$, $0.75 \leqslant T^* \leqslant 1.08$, and $1.8 \leqslant z \leqslant 3$. The MHNC pressures also appear in overall quite good agreement with MC. The overall performances of the MHNC for energy and pressure appear thus comparable to those for asymmetric mixtures, as can be verified by looking at figure 2 of [13]. As for the MHNC chemical potentials, these appear moderately overestimated (less negative) at $z=1.8$ (slowly decaying potential) irrespective of the density and temperature, and underestimated (more negative) for higher z values (short-ranged potential) at relatively low temperatures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

HRT results for pressure and chemical potential are rather close to those of MHNC, and overall at least equally accurate. The internal energies, while still remaining accurate, show deviations from the MC data slightly larger than those of MHNC. As a side remark, we may observe that the HRT results for the internal energy reported in the table, which have been obtained by numerical differentiation of the excess Helmholtz free energy $\beta A$ with respect to β, agree very closely with the values determined via (6) according to the internal energy route, the deviations ranging from a fraction of a percent for $z=1.8$ to  ∼$1\%$ for $z=3$. Hence, consistency with the internal energy route is very satisfactory, even though it has not been explicitly imposed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We turn now to structural results by examining the RDFs of the system. These are shown in figures 1–7. We first focus on MC results: the like RDFs $g_{11}(r)=g_{22}(r)$ do not display relevant variations for the $T^* - \rho^* - z$ combinations and intervals considered: indeed, contact values of the RDFs are almost uniformly ≃$3.5$ and oscillations appear comparable in all the various panels displayed. A similar remark holds for unlike correlations described by $g_{12}(r)$, with contact values ≃$2.5$ for all the cases considered.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Turning to the theoretical $g_{ij}(r)$ it clearly appears that the MHNC patterns agree quantitatively, almost everywhere, with the simulation curves. Compared to MHNC, the HRT shows some deviations from the MC results, especially at contact, where it underestimates $g_{ij}(r)$. Elsewhere, the agreement with simulation is good. The underestimation of the contact value does not come as a surprise since, as recalled in section 2, the form of the correlations adopted to close the HRT equation is similar to that of the RPA, which is also known to underestimate the RDFs at contact, especially for short-range interactions.

Another feature implied by the RPA-like form of the correlations assumed in HRT is that they do not fulfill the exact condition $g_{ij}(r)=0$ for $r<\sigma$, which stems from the hard-core interaction in (1). Even though this has not been pursued here, it is worth observing that, for the HCY potential, the above condition could be implemented analytically by turning to the smooth cut-off formulation of HRT, which differs from the sharp cut-off version used here in the prescription by which long-wavelength Fourier components of the interaction are built into the system as Q evolves from $\infty$ to 0. This formulation has been applied so far to one-component HCY fluids [30] and screened Coulomb mixtures [31]. The latter are still an instance of symmetric, equimolar HCYM, but correspond to the rather special case $\epsilon_{l} = -\epsilon_{u} < 0$, so that the interaction is attractive only for particles of different species.

Phase diagrams are finally shown in figures 8–10. Before presenting them, let us recall that symmetric equimolar mixtures display both LV and LL phase equilibrium. The topology of the phase diagram depends on the relative location of the LV coexistence curve and the line of consolute points, usually referred to as λ-line, and has been thoroughly investigated by mean-field theory [6], numerical simulation [6, 32], and by the SCOZA that we mentioned in section 1 [7, 9–11], all of which predict the same qualitative scenario.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When α is relatively large, while still being smaller than 1, the λ-line meets the liquid side of the LV binodal at a point of first-order coexistence, thereby terminating into a critical endpoint (CEP). At temperatures lower than T_cep, the fluid coexisting with the vapor actually consists of an equal amount of two demixed fluids with concentrations symmetric with respect to $c_{1}=1/2$, say $\bar{c}_{1}$ and $1-\bar{c}_{1}$. The resulting phase diagram will be hereafter denoted as type I.

At small and positive α, instead, the phase diagram changes to type III. In this case, the λ-line meets the binodal at its critical point, which is then replaced by a tricritical point. At the tricritical point, one has the simultaneous coalescence of three phases, namely, the mixed vapor and the two demixed fluids.

Moreover, in a narrow interval of α intermediate between those corresponding to type III and type I, the phase diagram displays both a LV critical point and a tricritical point at higher density. In this situation, which we shall refer to as type II, the phase coalescing with the two demixed fluids at the tricritical point is the mixed liquid instead of the mixed vapor as in type III. The three cases described above correspond, respectively, to figures 1(a), (c) and (b) of [6].

The results discussed below show that, according to both theory and simulation, the phase diagram belongs to type I for all the interaction parameters considered here. We shall further comment about the occurrence of type I topology within HRT in section 3.2.

We first examine the case $z=1.8$, $\alpha=0.75$ of figure 8 and primarily focus our attention on simulation results therein shown.

Our MC results for LV coexistence can be compared with grand canonical MC (GCMC) results obtained by other authors [11] for the same z and α: as visible, the two lines are fairly coincident with each other over their liquid side portions, while on the vapor side we estimate relatively lower coexisting densities. We thus predict a somewhat wider LV binodal; however, the overall curvatures of the two phase portraits seem consistent with almost equal critical temperatures which, upon visual inspection, appear to fall in the range $1.05 < T^*_{\rm cr} < 1.065$. A similar agreement emerges also with the binodal curve (not shown here) obtained by other authors in [12] via grand canonical transition matrix Monte Carlo (GCTMMC) simulation. In relation to the above estimates of the LV coexistence locus, we can confirm that at $z=1.8$, $T^*=1.08$ and $T^*=1.04$ are (slightly) super-critical and sub-critical temperatures, respectively.

We then consider theoretical binodal curves: the MHNC compares satisfactorily well with our MC data and with those of [11], although some deviation from simulation tends to manifest in the very low density region. The HRT binodal is in good agreement with both MHNC and simulation for relatively high densities. Along the vapor branch, the agreement between HRT and our MC results is good at low temperature, while at higher temperature the HRT curve is shifted towards higher densities, most likely because its critical temperature is lower than that predicted by the present simulation and, a fortiori, by that of [11]. The SCOZA binodal, which we reproduce from [11] for a direct comparison, is in overall agreement with our theoretical results, especially those of MHNC, for which the agreement is nearly quantitative.

As far as LL coexistence is concerned, we report in figure 8 GCMC results of [11] for the $\lambda$-line and compare it with our theoretical derivation. We recall that the MHNC λ-line is based on the vanishing of $\Delta s$, as discussed in section 2 and illustrated in figure 11; for the HRT derivation see [8] and section 3.2 for details. It appears that the MHNC and the HRT λ-lines, in rather close agreement with each other, are also in good agreement with the λ-line of [11], by only moderately overestimating the consolution densities at each temperature. The SCOZA λ-line, which is displaced to lower densities with respect to the MC one, as already known from [11], turns out to be in qualitative agreement with our theories.

Zoom In Zoom Out Reset image size

The intersection of the λ-line with the binodal determines the CEP at which a crossover occurs from LV to LL coexistence. As visible in figure 8, CEP parameters can thereby be estimated from the intersection of the MC and MHNC binodals with the respective λ-lines. In particular, for the MHNC our prediction is $T^*_{\rm cep}=0.955$ $\rho^*_{\rm cep}=0.62$ in quite good agreement with the simulation values of [11], namely $T^*_{\rm cep}=0.965$ $\rho^*_{\rm cep}=0.617$ (see table 2). As far as the HRT performance is concerned, we found that, for the present $z=1.8$ case, the numerical algorithm failed to converge in the interval $0.90 < T^{*} < 0.95$, the latter temperature being still slightly higher than that at which the LV binodal intersects the λ-line. As a consequence, even though the phase behavior down to $T^{*}=0.95$ is fully consistent with the presence of a CEP, the latter had to be obtained via extrapolation. We estimate in this manner that, within two significant figures, a CEP should occur at $T^*_{\rm cep}=0.94$ and $\rho^*_{\rm cep}=0.63$ (see table 2).

Table 2. CEP parameters $\rho^*_{\rm cep}$, $T^*_{\rm cep}$ for various α and z: MC simulation data and theoretical predictions according to various bibliographic sources.

$\alpha$	z	$\rho^{*}_{\rm cep}$	$T^{*}_{\rm cep}$	Source
0.75	1.8	0.617	0.965	GCMC [11]
		—	0.97	GCTMMC [12]
		0.59	0.97	SCOZA [11]
		0.63a	0.94a	HRT this work
		0.62	0.96	MHNC this work
0.75	2.4	—	0.73	GCTMMC [12]
		0.67	0.72	HRT this work
		≃$0.58$b	≃$0.755$b	MHNC this work
0.75	3.0	—	0.61	GCTMMC [12]
		0.71	0.61	HRT this work
		≃$0.57$b	≃$0.63$b	MHNC this work

^aHRT estimates obtained through extrapolation (see text). ^bobtained through the extrapolated intersection of the MHNC $\Delta s=0$ locus with the GCTMMC binodal of [12].

We turn now to the $z=2.4$ case. Our HRT results for the binodal line, together with those obtained via GCTMMC [12], are shown in figure 9. It appears that the HRT prediction is in fairly good agreement with simulation up to rather high densities. As for the MHNC, we did not determine the related binodal line. We have experienced in fact similar problems as those encountered in the already cited work by two of us [5]: we performed in that case MHNC and Gibbs Ensemble MC investigations of HCYM characterized by $z=2.45$ and $\alpha=0.7$, parameters only slightly different form the present ones. Although we showed in [5] a LV coexistence line, the related calculations turned out to be affected by instabilities of the solution algorithm, a circumstance which prompted us to conjecture that the pattern we exhibited was actually corresponding to a metastable equilibrium. We have experienced here even greater difficulties of the same type, and for this reason we do not report MHNC results for LV coexistence.

A comparison of figures 8 and 9 allows us then to appreciate the effect of varying z on the binodal line. It appears that the MC LV coexistence lines shifts to lower temperatures as z increases. In fact, one goes from $T^*_{\rm cr} \simeq 1.05$ at $z=1.8$, to $T^*_{\rm cr}\simeq 0.77$ at $z=2.4$. This effect is described qualitatively (although not quantitatively) already at the mean-field level, because of the effect of z on the integrated intensity of the interaction. The curve also sensibly flattens with the increase of z. A similar shift of the LV binodal has been observed and discussed in GCTMMC calculations for z in the range $1.8$–$3.0$ (see [12]), with critical temperatures which compare favorably with those here visible. It is worth mentioning that a similar downward shift in temperature, and flattening of the binodal curve, on increasing z (that is, under the reduction of the attractive potential range) is well known to occur in the one-component HCY fluid case (see [29, 33] and [18] for a review).

We now consider the HRT λ-line for $z=2.4$ displayed in figure 9. With respect to the $z=1.8$ case, the HRT LL coexistence appears now shifted to significantly lower temperatures, similarly to what happens for the LV curve. It also emerges that, at variance from the $z=1.8$ case, the intersection point of the HRT λ and binodal lines can now clearly be determined without resorting to extrapolation: we find $\rho^*_{\rm cep}=0.67$ and $T^*_{\rm cep}=0.72$. The latter temperature value is in fairly close agreement with $T^*_{\rm cep}=0.73$ of [12], and also fairly similar to the SGCMC estimate for relatively close z and α parameters as reported in [5]. We also show in figure 9 a MHNC determination of the λ line according to the $\Delta s=0$ locus procedure. Since, as said before, a MHNC binodal is not available for the present z value, we resort to the intersection of the MHNC λ-line with either the GCTMMC or the HRT binodal: as appreciable from table 2, the crossing occurs at a temperature in good agreement with that of [12]. As for the CEP density, the MHNC $\rho^*_{\rm cep} \simeq 0.58$ is quite close to the SGCMC result of [5].

We have finally explored the case $z=3.0$. Results are shown in figure 10. As for $z=2.4$, the HRT again reproduces fairly well the GCTMMC binodal line from [12]; the latter further shifts down in temperature with respect to figure 9. The HRT λ-line also shifts to lower temperatures with the increase of z, with respect to that shown in figure 9. HRT gives $T^*_{\rm cep}=0.61$, $\rho^*_{\rm cep}=0.71$. The MHNC estimate, similar to what discussed for the $z=2.4$ case, is also attainable and displayed in table 2.

The simulation evidences hitherto cumulated, for different z values and from various sources, about CEP parameters prompt us to conclude that the increase of z, that is the decrease of the attractive Yukawa potential range, moderately affects the crossover density from LV to LL equilibrium, but significantly reduces the crossover temperature. The latter trend appears in full agreement with what found by other authors [12] through the GCTMMC investigation. As for the capability of liquid state theories to correctly predict such a key behavior, the HRT appears with certainty able to reproduce the changes of the CEP temperature for z in the range $1.8 \leqslant z \leqslant 3.0$. The MHNC, which yields rather accurate predictions for the CEP parameters at $z=1.8$, suffers of convergency problems at higher z values in relation to the LV binodal determination; nonetheless, its predictions for the λ line, when combined with knowledge of the LV binodal from other sources, lead to rather accurate estimates of the CEP temperature; As for the CEP density, both HRT and MHNC show a good agreement with simulation for $z=1.8$. At higher z values a sharp assessment is made difficult by the paucity of simulation data for $\rho_{\rm cep}$. A related issue worth being ascertained is the influence of z on the CEP density at fixed α: as shown in table 2, on increasing z, HRT predicts an increase of $\rho^{*}_{\rm cep}$, even though this effect is less marked than the corresponding decrease of $T^{*}_{\rm cep}$. Conversely, MHNC supplemented by the $\Delta s=0$ criterion predicts a decrease of $\rho^{*}_{\rm cep}$. At this stage more simulation estimates of the CEP density at different z values would be helpful in order to establish which trend does actually take place; in this respect, a check of the entropic criterion accuracy for the determination of the λ-line at high z values might also be in order.

A HRT study of the phase diagram of symmetric HCY mixtures for $z=1.8$ as a function of α has been carried out in [8]. As stated in section 2, the HRT PDE involves both density and concentration as independent variables, so that the phase diagram was investigated at all concentrations. In this way, it was possible to establish how the features of the three classes mentioned above were reflected in the overall phase portrait in the $(\rho^*, c_{1}, T^*)$ space. However, that study failed to give a clear evidence of the CEP scenario of type I expected for relatively large α, suggesting instead that type II might persist even for α approaching 1. Such a situation is at variance with the results reported here, according to which HRT agrees with both simulation and other theories in predicting type I in the regime of α where it is expected to appear. Therefore, the origin of this discrepancy must be explained.

To this end, we have considered the evolution of the HRT phase diagram on the $\rho^*$–$c_{1}$ plane as $T^*$ is decreased. Here we shall refrain from a detailed analysis of all the possibile scenarios as α is changed, for which we refer the reader to [8], and focus on one of the cases studied here, namely, $z=2.4$, $\alpha=0.75$ whose phase diagram at $c_{1}=1/2$ has been displayed in figure 9.

Figure 12 displays four isothermal sections of the phase diagram, each of which has been obtained by integrating the HRT PDE on a $200\times 200$ grid. Panel (a) features two disconnected coexistence regions. The region at lower densities spanning the whole concentration axis, hereafter denoted as region A, results from two LV domains, which originate from the coexistence regions of the pure species, and get larger and larger as $T^*$ is decreased. At temperatures lower than $T^*_{\rm cr}$, such as those of the figure, these domains merge at $c_{1}=1/2$, at which LV equilibrium is azeotropic, thereby generating the LV coexistence curve of figure 9.

Zoom In Zoom Out Reset image size

The λ-line is instead related to the domain at higher densities in panel (a), hereafter denoted as region B. At high $T^*$, that domain involves only one kind of phase equilibrium, namely, a demixing transition such that the coexisting phases have equal densities, and concentrations symmetric with respect to $c_{1}=1/2$. As one moves to lower densities, the concentrations of the coexisting phases get closer and closer, until at $c_{1}=1/2$ they merge at a demixing critical point, which is located at the low-density boundary of B, and gives the position of the λ-line.

The situation at low $T^*$ is displayed in panel (d). Here, A and B are no longer disconnected, but have merged into one single coexistence region. The vapor at $c_{1}=1/2$ is now in equilibrium with two fluids of equal densities and symmetric concentrations, so that the fluid at equimolar concentration coexisting with the vapor actually consists of an equal amount of two demixed fluids, as stated in section 3.1.

The scenario described so far is common to both type I and type II phase diagrams. The difference between the two classes is determined by how the two disconnected regions A and B of panel (a) merge into the single region of panel (d). In type I, the kind of phase equilibrium on B remains the same as that described above at all temperatures such that it has not yet met A. In particular, the low-density boundary of B at $c_{1}=1/2$ remains a critical mixing-demixing point, until at $T^*_{\rm cep}$ it meets A at the CEP. In type II, instead, this point first changes into a tricritical point at a temperature $T^*_{\rm tri}$ higher than that at which A and B touch, and then bifurcates into two critical points at $c_{1}\neq 1/2$ for $T^*<T^*_{\rm tri}$. These are the points which come first into touch with A as T is further lowered. Because contact initially occurs at non-equimolar concentration, there exists a narrow temperature interval in which the phase diagram shows an 'island' of homogeneous fluid enclosed in the coexistence region. This feature is characteristic of type II.

If we now consider panels (b) and (c) of figure 12, we find that they apparently conform to the scenario just described for type II, even though the homogeneous island is very narrow in density. This situation is similar to that represented in panels (c) and (d) of figure 13 of [8] for a HCY mixture with $z=1.8$, $\alpha=0.8$. As discussed there, it suggests that according to HRT, type II topology may survive at values of α significantly higher than those predicted by simulation, mean-field theory, or SCOZA, and in principle it could always prevail over type I even though, for α approaching 1, the size of the island would become so small, that telling one type from the other would become virtually impossible.

However, other features seem to conflict with such a behavior: first, the aforementioned appearance of two critical points below $T^*_{\rm tri}$ should be detectable by plotting the phase diagram in the P–$\Delta\mu$ plane, $\Delta\mu$ being the difference of the chemical potentials of the species. On this plane, domain B collapses into a line at $\Delta\mu=0$ which below $T^*_{\rm tri}$ should bifurcate into two symmetric branches at $\Delta\mu\neq 0$, but we failed to observe such a bifurcation for the cases displayed in figure 12. This situation was already reported in [8], where it was attributed to the numerical difficulties in detecting a non-vanishing $\Delta\mu$ when the densities of coexisting phases are very close. Moreover, at temperatures just below T_tri, the equimolar fluid at the low-density boundary of B should coexist with a demixed equimolar fluid of slightly higher density, but again we were not able to observe this. Once more, one could blame the very small difference between the densities at coexistence. Nevertheless, the resolution allowed by the grid employed in figure 12 should be high enough to resolve this difference, even though it might amount to just a few gridpoints.

Prompted by this conflicting evidence, we have reinvestigated the phase diagram of the same $z=2.4$, $\alpha=0.75$ mixture using an even larger $500\times 200$ grid in order to represent accurately the very flat shape of domains A and B near contact. The results in the interval $0.5\leqslant\rho^{*}\leqslant 0.8$ are displayed in figure 13. Clearly, the alleged intersections of the domains at $c_{1}\simeq 0.4$ and $c_{1}\simeq 0.6$ and the resulting island of mixed fluid have now disappeared, consistently with the scenario characteristic of type I. We are then led to conclude that HRT does predict a CEP for this mixture, and that the situation displayed in panels (b) and (c) of figure 12 does not correspond to a real contact between domains A and B. Rather, the two domains are actually still disconnected, but they run so close to each other, that their small distance cannot be resolved by the density grid used in the figure. Similar considerations hold for the aforementioned case $z=1.8$, $\alpha=0.8$, which in [8] was investigated by means of a $150\times 150$ grid. Hence, in the light of the present investigation, HRT conforms to simulation, mean-field theory, and SCOZA in predicting that, for large enough α, the phase diagram will display a CEP at equimolar concentration.

Zoom In Zoom Out Reset image size