Monte Carlo simulations of the critical properties of the restricted primitive model

Recent Monte Carlo simulations of the critical point of the restricted primitive model for ionic solutions are reported. Only the continuum version of the model is considered. A finite size scaling analysis based in the Bruce-Wilding procedure gives critical exponents in agreement with those of the three-dimensional Ising universality class. An anomaly in the scaling of the specific heat with system size is pointed out.


I. INTRODUCTION
The primitive model (PM) for electrolytes, molten salts, colloïds, etc is a mixture of M species of charged hard spheres living either on a lattice or within a continuous volume of real space. In this paper we shall focus only on the off-lattice version of the model. The simplest version of the PM consists of a binary mixture (i.e. M = 2) of positive and negative charged hard spheres ±q all with the same diameter σ. Under this form the model which is thought to be the prototype of many ionic fluids has been christened the restricted primitive model (RPM). A thermodynamic state of the RPM is entirely specified by a reduced density ρ * = Nσ 3 /V (N number of ions, V volume) and a reduced temperature T * = kT σ/q 2 (k Boltzmann's constant).
The RPM undergoes a liquid-vapor transition which has been studied extensively these last past years by means of Monte Carlo (MC) simulations and various theoretical approaches. In particular the behavior of the system at its critical point (CP) has been the subject of a huge amount of numerical and theoretical studies. The question is obviously of great importance since it is reasonable to assume that real electrolytes -or at last a large class of them-and the RPM belong to the same universality class which dictates a similar critical behavior.
It is perhaps the right place to note that the main feature of ionic solutions is that the pair potential between two ions i and j at a distance r ij = |r i − r j | which reads as v(r ij ) = q i q j r ij for r ij > σ = +∞ for r ij < σ (1.1) is a long range interaction. This fact would suggest classical (mean field) behavior, whereas the well-known screening of the interactions pleads in favor of an Ising-like criticality typical of systems with short range interactions.
On the experimental side it seems well established now that for many real electrolytes apparent mean field behavior applies with sharp crossover (much sharper than in nonionic fluids) to Ising criticality close to the critical temperature 1 . France and, at this point of the story, I must agree with him retrospectively. However many advances have been done since. Before giving an account of these new achievements some comments are in order.
• (i) All the MC studies confirm the existence of a liquid vapor transition for the RPM.
It seems to take place at unusually low densities and temperatures. Caillol and Weis give further support for such a low critical temperature 14 . Moreover it turns out that the coexistence curve is very dissymmetric 9,11,12 .
• (ii) The MC simulation of ionic systems is a numerical challenge due to the long range of Coulomb potential. In order to deal with this, some caution is needed. Thus, in the case of MC simulations performed in a cubic box with periodic boundary conditions (PBC), one must use Ewald potentials in order to obtain the correct physics 15,16,17 .
The point is that the Ewald potential is the solution of Poisson equation in a cubicoperiodical geometry 17  • (iii) None of the above mentioned studies took correctly into account finite size effects which are of an overwhelming importance near a CP. These effects affect the behavior of finite systems as soon as the correlation length of the critical density fluctuations is of the same order of magnitude than the size of the simulation box. In the simulations 9,11,12 some "apparent" critical temperature T * c has been measured which could be very different from its infinite volume limit T * c (∞).
In order to extract from MC simulations the critical behavior of the RPM in the thermodynamic limit (i.e. the critical exponents) and also the infinite volume limit of T * c and ρ * c it is necessary to perform an analysis of the MC data in the framework of the finite size scaling (fss) theory which is part of the renormalization group (RG) theory 18,19 . In this approach one needs to work in the Grand Canonical (GC) ensemble rather than in the Gibbs ensemble which is ill adapted for a fss analysis. Subsequent MC simulations on the RPM were thus all performed in this ensemble. Panagiotopoulos and coworkers turned their attention to the lattice version of the RPM whereas the Orsay group continued to work on its off-lattice version.

A. Scaling fields and operators
Starting with the seminal work of Bruce and Wilding (BW) 20,21,22 simulation results for the critical behavior of fluids have customarily been analyzed along the lines of the so-called revised scaling theory of Rehr and Mermin 23 . In this approach one first defines scaling fields and operators aimed at restoring the particle-hole symmetry and therefore to map the the fluid onto a magnetic sytem with Ising-like symmetry.
The two relevant scaling fields h (the strong ordering field) and τ (the weak thermal field) are assumed to be linear combinations of deviations from their critical values of the chemical potential µ and the inverse temperature β = 1/T (reduced values are assumed henceforward). One thus has where r and s are the field mixing parameters which define the mapping. Of course relations where Ξ is the GC partition function of the RPM, ρ the total number density, and u the internal energy per unit volume. Brackets < . . . > denote GC averages. M is the order parameter (magnetization) of the magnetic system associated with the fluid and E its mag- The revised scaling of Rehr and Mermin implies the analyticity of the coexistence chemical potential µ(T ) at T * c . Although this is the case for some peculiar lattice gas models with "hidden" symmetries there is no reason that in general, for fluid systems µ(T ) should lack a singularity as recognized already by Rehr and Mermin 23 and emphasized more recently by Fisher and co-workers 40,41,42 .

B. The scaling hypothesis
A central role in the subsequent fss analysis is played by the GC joint distribution P L (M, E) ∝ P L (ρ, u) for the scaling operators M and E. Following BW 20,21,22 we will assume that, in the immediate vicinity of the CP, P L (M, E) obeys to the following scaling law: in terms of the usual critical exponents: • β exponent of the ordering field, i.e. < δM >∼ |τ | β for T * < T * c at h = 0 • ν exponent of the correlation length, i.e. ξ ∼ |τ | −ν • θ Wegner's correction-to-scaling exponent (first irrelevant exponent).
The scaling hypothesis (3.3) was established on a solid RG basis for Ising-like systems 24 and received substantial supports from MC studies 25 . We stress once again that the coexistence curve is determined in this approach by the condition h = 0 and that, at coexistence, the order parameter distribution P L (M) should be an even function of M. In practice this symmetry requirement can be satisfied by tuning the two parameters (µ, s) at a given β.
We now concentrate our attention on the scaling behavior of the histogram P L (M).

C. The matching procedure
Integrating both sides of eq. (3.3) over E one finds that, along the coexistence line h = 0 one has where, in the r.h.s. the dependence of the universal function P upon h has been discarded for clarity. Let us define now x = a −1 M L d−y h δM, then, assuming τ ∼ 0 and L ∼ ∞ a Taylor expansion of eq. (3.5) yields where the various P * entering the r.h.s. are universal functions. Note that, for L = ∞ the normalized ordering field distribution P L (M) collapses onto an universal function P * (x) at τ = 0. For L finite but large P L (M) collapses approximately onto P * (x) at some apparent τ L ∝ L −yτ +y i . Since for h = 0 one has τ ∝ β −β c then the matching of the histogram P L (M) onto the universal function P * (x) should occur at some apparent temperature T * c (L) scaling with system size as where T * c (∞) denotes the infinite volume limit of the critical temperature.

D. Technical details
To assess the critical behavior and the critical parameters of the system, we need, in a first step, to locate the coexistence curve h = 0. At a given temperature β close to β c the ordering distribution function P L (M) depends solely on the chemical potential µ and the mixing parameter s. At coexistence, the value of (µ,s) can be obtained unambiguously by were realized by using the estimate of P * is (x) made by Hilfer and Wilding 32 for the 3D Ising model. Two new -and better-estimates of P * is (x) obtained by Tsypin and Blöte 33 for the 3D Ising model and the spin-1 Blume-Capel model were considered in ref. 31 . The discussion is postponed to next section.

A. General discussion
It turns out that the field mixing parameter s of the RPM is practically independent of the temperature and of the size L of the system. Its magnitude, s ∼ −1.46 29,30,31 , is much higher than for neutral fluids (typically s ∼ 0.02 for square well or Lennard-Jones fluids 34 ) which explains the large dissymmetry of the liquid-vapor coexistence curve of the RPM.
The collapse of the ordering operator distribution P L (M) onto the universal ordering distribution P * is (x) given by the Blume-Capel model 33 is depicted in Fig. 1 for four different values of the volume ranging from V /σ 3 = 5000 to V /σ 3 = 40000, i.e. up to a linear size L/σ = 34. At volume V /σ 3 = 5000 a mismatch is observed at the lowest values of M due to an inadequate sampling of of the low density configurations at small volume. The overall good agreement leads us to conclude that the universality class of the RPM is that of the 3D Ising model.
The reduced apparent critical temperature T * c (L) versus the size L of the system (in reduced units) has been plotted in Fig. 2. Depending on the choice made for the universal ordering distribution P * is (x) one obtains two sets of values of T * c (L) from which T * c (∞) can be obtained by making use of eq. (3.7). One obtains T * c (∞) = 0.04917 ± 0.00002 using P * is (x) derived from the Blume-Capel model and T * c (∞) = 0.04916 ± 0.00002 using P * is (x) obtained for the 3D-Ising model. The approximate P * is (x) of Hilfer and Wilding yields slightly different results. Note that in all cases we have used the Ising values ν = 0.630 35 and θ = 0.53 36 of the critical exponents.
The previous analysis merely establishes the compatibility of the MC data with an Isinglike criticality. One can try to go beyond by considering the scaling behavior of the Binder As a consequence of the scaling hypothesis (3.6) one can show that, at coexistence (h = 0), Q B (L) should scale with system size as where the last term takes into account contributions from irrelevant fields and q 1 , q 2 , q 3 , and b 1 are non universal constants. If the contribution of irrelevant fields could be neglected then the curves Q B (L) would intersect at the fixed point Q c . As apparent in Fig. 3 this is clearly not the case and corrections to scaling must be taken into account.
Recall that for the 3D-Ising model the fixed point value is Q c = 0.623 37 and that the exponent of the correlation length has the value ν = 0.630 35 . We have attempted to fit all our MC data with eq. (4.2). If all the parameters in the RHS of eq. (4.2) are kept free such an ambitious fit turns out to be impossible. Various other less satisfactory strategies can be considered however.
The variations of Q B (L) as a function of β for the different volumes is shown in Fig. 3.
Although there is considerable spread in the intersection points due to correction-to-scaling contributions, the corresponding values of Q c are close to the Ising value Q c = 0.623 and permit to rule out mean field behavior (i.e. Q c = 0.457 38 ).
Further support for Ising criticality is provided by the behavior of < δM 2 > at T * c (L). According to the scaling hypothesis (3.6) it should scales as L 2β/ν with system size. From the slope of the curve displayed in Fig. 4 one obtains β/ν = 0.52 in accord with the 3D Ising value (0.517) and in clear contrast with the classical value 1.
In summary, our fss analysis leads to an estimate of the critical exponents ν and β/ν and the Binder cumulant Q c based on the sole knowledge of the critical temperature and the renormalization exponent θ. Within the numerical uncertainties these values are compatible with Ising-like criticality. Our conclusion is that the RPM, as ordinary neutral fluids, belongs to the universality class of the Ising model.
A complete discussion of our MC data is out of the scope of the present paper and can be found in ref. 31 . For completeness I give below the values obtained for the critical temperature, chemical potentials and densities (the critical pressure is largely unknown): • T ⋆ c = 0.04917 ± 0.00002 • ρ ⋆ c = 0.080 ± 0.005 • µ ⋆ c = −13.600 ± 0.005

B. The specific heat
The revised scaling theory of Rehr and Mermin which is the framework of our fss analysis is however not the most general scaling theory which can be proposed for a fluid system lacking the "particle-hole" symmetry. Its main weakness, as recognized already by Rehr and Mermin 23 , Yang and Yang 39 , and more recently by Fisher and co-workers 40,41,42 , is that it assumes the analyticity of the chemical potential at coexistence µ(T ) at the critical point.
The more general scaling assumption should yield singularities for both µ(T ) and p(T ) as T * → T * c . Let us examine the consequences of these singularities on the behavior of the specific heat capacity at constant volume C V . In the two phase region it can be rewritten as 39 where C p (not to be confused with the specific heat capacity at constant pressure) and C µ (not to be confused with the specific heat capacity at constant chemical potential) denote the two contributions to C V . I stress that, in eq. (4.3) p(T ) and µ(T ) denote the pressure and the chemical potential at coexistence. The formula can be used for any density ρ g (T ) < ρ < ρ l (T ) within the two phase region (ρ g (T ) and ρ l (T ) being the densities of the gas and the liquid at coexistence respectively). In the revised scaling theory only C p diverges as |T * − T * c | −α whereas one expects a divergence of both C p and C µ (both as |T * − T * c | −α ). In Fig. 5 I display the curves C µ (T ) and C V (T ) along the locus