Monte Carlo simulation of phase equilibria in Ising fluids and their mixtures

The mean field theory for the pure Ising fluid was recently extended to binary mixtures of an Ising and a van der Waals fluid. Depending on the relative interaction strengths, their three dimensional phase diagrams show lines of tricritical consolute and plait points, lines of critical end points and magnetic consolute point lines. Our current efforts are to compare these mean field results with different Monte Carlo simulation techniques, investigating both first order (liquid-vapor and demixing) and second order (paramagneticferromagnetic) phase transitions. We show the resulting ρ, T phase diagrams of the pure Ising fluid for different magnetic interaction strengths R and constant pressure cross-sections of the x, T, p phase diagrams of Ising mixtures for different relative interaction strengths. The methods we have used include Gibbs Ensemble MC, Multihistogram Reweighting, Hyper-parallel Tempering, the cumulant intersection method and the newly developed Density of States MC technique.


Introduction
The investigation of continuum fluid models with coupled translational and spin degrees of freedom is of current theoretical interest [1][2][3].The importance of such models lies in their display of a rich variety of transitions between solid, liquid, and gas, as well as magnetically ordered and disordered phases, which may occur in real systems.The liquid-gas and magnetic para-ferro phase transitions in spin fluids were previously studied using the mean field theory [4][5][6], the method of integral equations [7][8][9][10][11], and Monte Carlo (MC) simulation techniques [2,7,10,[12][13][14][15].
In [3] we presented the mean field theory for a binary mixture of an Ising fluid and a van der Waals fluid and studied the topology of the resulting three-dimensional phase diagrams.Our interest in this system is founded, on the one hand, on its extremely rich phase behavior, including, for example, lines of tricritical points and consolute point lines in the ferromagnetic phase, and, on the other hand, in its connection to ternary mixtures which is due to the fact that the Ising fluid can be mapped onto a symmetric binary mixture [16].Moreover, the Ising system exhibits features that are also found in other more realistic models, like, for example, the restricted primitive model (RPM) of ionic fluids [17].
The outline of this paper is as follows: In section 2 we summarize the mean field theory of Ising mixtures from [3].Section 3 deals with MC simulations for the pure fluid system, the applied methods and the results, and section 4 contains the corresponding treatment for binary mixtures.

Mean Field theory
We consider the general case of a binary mixture in the molar volume V consisting of a van der Waals fluid (fluid 1) and an Ising fluid (fluid 2).The mole fraction of the second component shall be denoted as x and its magnetization per particle as m.The case of a pure Ising fluid follows for x = 1.We describe our system with a van der Waals-like equation of state, where the attraction parameter a is defined according to the quadratic mixing rule as and the size parameter b is assumed constant.In equation ( 2), a 11 and a 22 denote the nonmagnetic interactions between particles of the same kind, a 12 the nonmagnetic interaction between unlike particles and a m the magnetic interaction in the Ising fluid.For the magnetization, the equation of state at zero magnetic field reads The corresponding molar Helmholtz free energy of the system described by ( 1) and (3) with respect to the reference state of an ideal unmixed gas with molar volume V 0 = b is given by where is the entropy part of the free energy of the Ising fluid component.The mixture can be described by the three reduced parameters ζ, Λ and R m , and we introduce dimensionless variables scaled by the critical parameters of component 2, As was shown in [4], the phase diagram of the pure van der Waals Ising fluid can have three different topologies: • type I: For R m > 0.304 the second order magnetic transition line T λ ends in a tricritical point below which a first order phase transition in the liquid appears.This is the scenario also valid in the "ideal" case when → ∞ [18].
• type II: For 0.211 < R m < 0.304 the second order magnetic transition line T λ ends in a tricritical point below which a first order phase transition in the liquid appears which ends in a triple point on the gas liquid first order phase transition.This first order line ends in a critical point.Thus, a tricritical and a critical point coexist.
• type III: For R m < 0.211 the second order magnetic transition line T λ ends in a critical end point on the gas liquid first order phase transition.This first order line ends in a critical point.
In the mixture, the conditions for the coexistence of two phases α and β at a temperature T 0 and pressure p 0 are given by By solving these equations numerically, one can find the first order phase transition surfaces in x, T, p-space.Lines of second order phase transitions (plait point and consolute point lines) are calculated from the equations where G m = A m −∂A m /∂V is the Gibbs free energy.In the case of a tricritical line the first of these conditions taken in the limit of m → 0 yields an analytic expression for the tricritical temperature as a function of the mole fraction x.It turns out that depending on the three mixture parameters, the tricritical line has either the character of a consolute line and goes to infinite pressure, or it resembles a plait point line and eventually becomes unstable in a tricritical endpoint.In the three dimensional parameter space, the boundary between these two types of behavior is described by the equation

Monte Carlo simulations I -Pure Ising fluid
To investigate the phase diagram of a pure magnetic fluid we have modelled the pair potential between two particles i and j as a sum of a nonmagnetic and a magnetic interaction part.For the first contribution we used a conventional Lennard-Jones potential, while for the second we used a Yukawa potential multiplied by the two Ising spins connected with the particles: with Here, σ is the size of the particles, r = | r i − r j | is the inter-particle distance and s i is the spin of particle i, which is either +1 or −1.The potential was truncated at a cutoff-distance r c equal to half the box length.This truncation was accounted for by adding a tail correction to the potential energy of each particle.The parameter R was used to vary the relative strength of the magnetic interaction.Simulations were performed for systems with R = 0 (pure Lennard-Jones fluid), R = 0.05, R = 0.075, R = 0.1, R = 0.115 and R = 0.1333.The phase diagrams of these fluids were investigated via three different simulation methods: Gibbs Ensemble MC, Grand Canonical simulations with Multihistogram reweighting and Density of States MC.
In the GEMC [19,20] simulations, the total number of particles N was usually 500 or 1000.The simulations were performed in cycles consisting of N trial particle displacements, N spin flip attempts, one volume change attempt and a number of particle swap attempts which was chosen to yield about 1-3% acceptance and depended strongly on the temperature.Typically, the number of the performed MC cycles was about 5 × 10 5 .
Additionally, simulations in the Grand Canonical Ensemble have been performed.These were conducted with a system size of L = 8σ or L = 10σ.Typical runs consisted of 2 × 10 9 MC steps, where one step was either an attempted deletion or insertion of a random particle.Particle displacements and spin flip attempts were not implemented explicitly.The resulting histograms were reweighted using the algorithm from [21,22].In order to speed up the simulation process, the simulation runs at different state points for the same system were performed in parallel and configurations of two such runs could be interchanged according to the Hyper-Parallel Tempering scheme described in [23].
Since at low temperatures and high densities of the liquid phase, both the GEMC and the GC simulations suffer from a poor acceptance rate of particle insertion attempts, we also applied another Monte Carlo method which was recently introduced for lattice [24,25] and off-lattice systems [26,27] and which is especially well suited for low temperatures, namely the Density of States (DOS) method.This approach uses the so-called uniform ensemble in which all configurations are equally probable.The considered MC moves consist of random particle displacements, particle insertions and deletions and in our case also spin flip attempts.The discretization of energy values was made with a bin size equal to the reduced unit of energy, ε.For each system and particle number, the energy range corresponding to the temperature interval we were interested in was determined beforehand from short simulation runs in the canonical ensemble at the highest and the lowest considered temperature.The paramagnetic-ferromagnetic phase transition line (Curie line) was investigated using the Binder crossing technique [28].During simulations in the canonical ensemble, the two-dimensional energy-magnetization histograms were obtained and extrapolated to other temperatures via the multi-histogram reweighting technique.Thus, the fourth-order cumulant ratio of the magnetization distribution P (m) could be calculated as a function of temperature.This was done for three different system sizes, N = 150, N = 300 and N = 500, and the intersection point of the resulting curves was determined as an estimate for the critical temperature of the infinite system.We found that the slope of the Curie line is in fact proportional to the parameter R, as predicted by mean field theory.From the systems with R = 0.05 and R = 0.075, where we calculated two points of the Curie line, we obtained a proportionality constant c = 27.1.For the other parameters we assumed the Curie line to be given by T r = cRρ * and determined one transition point in order to check if it really falls onto this line.
The phase diagrams we found in the simulations (figure 2 (a)) correspond well to the ones predicted by mean field theory (figure 2 (b)).For small values of R, 0.05 and 0.075, we find a critical end point where the Curie line intersects the phase coexistence curve.Above the critical end point temperature, these phase diagrams are almost indistinguishable from the case of the pure Lennard-Jones fluid (R = 0) which is also included in figure 2 (a).For larger values of R, e. g.R = 0.133, the shape of the coexistence curve changes and there is no critical point any more.Instead, a tricritical point seems to occur where the Curie line meets the binodal.For R values that lie in between these two regimes, mean field theory predicts the existence of both a critical and a tricritical point.For the parameter we considered, R = 0.1, our simulations in fact showed a triple point at T * = 1.32 with para-para as well as para-ferro transitions directly above the triple point temperature.However, a direct comparison between phase diagrams with a certain value of R in MC simulations and R m in mean field theory is not possible since there is no direct correspondence between the two ratios.We demonstrate in figures 2 (a) and (b) that the same topologies are found in the global phase diagram.
At H = 0 the Ising fluid may be mapped onto a symmetric binary mixture and comparison with other methods is possible.Our results show the same scenario as MC simulations in an Ising liquid with both interactions, the non-magnetic and the magnetic one, of square well [16] or Lennard-Jones type [29].It turned out that the hierarchical reference theory (HRT) [30] did not agree with the scenario found in MC simulations.In the language of the magnetic fluid, for small values of R no critical end point [29], but a phase diagram of type II was found [30].
The complete phase diagram of a magnetic fluid includes the magnetic field H.Recently the magnetic field dependence of the critical point location in an ideal Ising fluid has been considered by MC simulations [31] supporting the existence of a tricritical point in the Ising fluid at H = 0.The H-dependence at small fields is in agreement with the H 2/5 law for the wings from MF theory.In the limit H → ∞, a monotonic decrease of the wing line to the critical point of the Yukawa fluid has been found, being also in agreement with mean field theory.Moreover, it turned out that the MF phase diagram topology of type II splits into two parts depending on whether the critical point or the tricritical point is connected to the critical point at saturation (H → ∞).This leads to the appearance of a van Laar point in the global phase diagram and at this point to a phase diagram containing an unsymmetrical tricritical point at finite magnetic field [32].This is in agreement with the scenario found in mean field theory for symmetrical mixtures [33] and the HRT [30].

Monte Carlo simulations II -Mixture with an Ising fluid
For binary mixtures, we used the same interaction potential as for the pure fluid, However, the spin interaction part only enters if both particles belong to the magnetic component.Also, the Lennard-Jones interaction strength ε now depends on the species a and b of the two particles, and the Yukawa potential is accordingly given by We can then define parameters ζ and Λ in analogy to (5) with ε 11 , ε 12 and ε 22 .Again, we have added a tail correction term to the energy to compensate the truncation of the potential.
In our simulations we looked at the two systems with parameter values we expected to lead to similar phase diagram topologies as the two types in our mean field theory, namely ζ = 0.5, Λ = −0.05,R = 0.05 and ζ = 0.5, Λ = −0.25,R = 0.1333.These values of R were chosen as representatives of type I and III Ising fluids according to the results of the simulations for the pure Ising fluid.
The phase diagrams of these systems were investigated with GEMC simulations (N = 1000 or 500) at constant pressure in the case of first order phase transitions, as well as simulations in the constant N pT ensemble (N = 150, 300 and 500) and the cumulant crossing technique for the second order para-ferro phase transitions.Moreover, the first order transition was also studied within an isobaric Semigrand Canonical ensemble (N = 300) at constant N, ∆µ, p, T using Multihistogram reweighting.The MC moves performed in this ensemble consist of particle displacements, spin flips, particle identity changes and volume rearrangements.Thus, both liquid-vapor and demixing transitions could be studied.Figures 3 and 4 show the results for the first system at p * = pσ 3 /ε 22 = 0.08 and p * = 0.15, figure 5 for the second one at p * = 0.24.It turns out that the observed phase behaviors indeed resemble those of two mean field systems belonging to either of the two topology classes (at least in the temperature ranges accessible to our simulations).In the figures, isobaric sections of the mean field x, T, p-diagrams are presented together with the corresponding Monte Carlo results in each case.Regarding the values of the ratios R and R m the same remark as in the case of pure Ising fluids has to be made.We can only say that similar topological features are recovered by MC simulations.
At p * = 0.08, the first system exhibits a critical end point where the Curie line intersects the phase coexistence curve, separating first order para-para phase transitions occurring at temperatures above the critical end point temperature from first order para gas -ferro liquid transitions below this critical temperature.At higher pressure, p * = 0.15, the Curie line does not connect to the para-para liquidvapor two-phase region.Instead, it obviously passes into a para-ferro first order demixing transition via a tricritical point.

Figures 1 (
a) and (b) show the complete three-dimensional phase diagrams of systems corresponding to the two topologies.

Figure 2 .
Figure 2. (a) Results of GEMC simulations (diamonds), multihistogram reweighting (squares) and DOS MC simulations (lines) for Ising fluids with different interaction ratios R. The paramagnetic-ferromagnetic phase transition points (circles) were found via the Binder crossing technique, the dashed lines are linear extrapolations.T * = k B T /ε, ρ * = ρσ 3 .(b) Mean field phase diagrams of pure Ising fluids with different interaction ratios R m = 0, 0.1333, 0.2, 0.26, 0.35, 0.5.M, tricritical point; C, critical point; solid lines, first order phase transitions; dashed lines, second order (magnetic) phase transitions.The line of all tricritical points for R m > 0.211 is also shown.