The demixing phase instability in the ion-dipole mixture

By considering the variation of the grand potential functional Ω with respect to small density fluctuations δρα(1) we can determine the phase instability of a system from the correlation functions. Using the reference hypernetted chain (RHNC) approximation the correlation functions of an iondipole mixture have been calculated. With the obtained correlation functions the phase behaviour of ion-dipole mixtures has been investigated.


Introduction
A mixture of charged hard spheres and dipolar hard spheres is a simple model of an electrolyte.For such a model a series of statistical mechanical investigations have been done with integral equations as well as with simulations.In the mean spherical approximation (MSA) the analytic solution of the ion-dipole mixture can be obtained [1,2].Using the MSA results the phase transition in ion-dipole mixture has been discussed by Høye et al. [3] and Harvey [4].But it was known that MSA results are not sufficient for the ion-dipole mixture because of the strength of the interactions.The correlation functions of ion-dipole mixtures can also be calculated in the so-called linearized, quadratic [5] and reference hypernetted-chain [6] (LHNC, QHNC and RHNC) approximations with numerical iteration methods.Comparing with the Monte Carlo simulation results [7][8][9] only the RHNC approximation works very well [6,7], and here we will calculate the correlation functions in the RHNC approximation.
In the RHNC approximation, for moderate dipolar moments a convergent solution for ion-dipole mixtures cannot be obtained when the charge of the ions is still smaller than an elementary charge.This was already recognized by Dong et al. [10].Also the simulations encountered problems in this range of parameters [11].
By analyzing the phase stability of the system with respect to the small density fluctuations of ions and dipoles we have shown that this problem is related to the phase instability of the ion-dipole mixture [12].In this article we will give a more detailed investigation of the ion-dipole phase instability.
The article is organized as follows.In section 2 we define the ion-dipole mixture model in detail.Section 3 discusses the calculation of the correlation functions.In section 4 we describe the method for the phase instability analysis.We use the calculated correlation functions to discuss the phase instability of the ion-dipole mixture in section 5. Section 6 contains the conclusions.The appendix describes the treatment of the long range tails of the interactions and the correlations.

The model
Our model electrolyte is a mixture of hard sphere ions and hard sphere dipoles.The ions have the charge q i (q + = q, q − = −q) and number density ρ i .The dipoles have dipole moment µ and number density ρ d .Both ions and dipoles have the same hard sphere diameter σ.The system is taken to be electrically neutral: i q i ρ i = 0. (2.1) The pair potential between two particles in the system can be written as u ij (12) = u HS (12) + q i q j /r, (2.2) (2.4) The variables 1, 2 imply position and orientation (for the dipoles) with 2 ≡ ( r 2 , μ2 ), where μ2 is the unit vector in the direction of the dipole moment at r 2 with the Euler angles ω 2 ≡ (θ 2 , ϕ 2 ).r = | r 2 − r 1 | and r12 is the unit vector pointing from position 1 to position 2: r12 = ( r 2 − r 1 )/r, with its Euler angles denoted as ω r .u HS (12) is the hard sphere interaction defined by The function D (12) is the angle dependent part of dipole-dipole interaction and is defined as

The calculation of the correlation functions
The correlation functions of ion-dipole mixtures can be calculated with reliable results by the Ornstein-Zernike (O.Z.) equation and the reference hypernetted chain (RHNC) closure.This conclusion was reached by Caillol et al. [7].They have done Monte Carlo simulations for the ion-dipole mixture and have shown that the RHNC approximation works well.
The Ornstein-Zernike equation gives the relation between the total correlation functions h αβ (12) and the direct correlation functions c αβ (12) (32) . (3.1) Here ρ α (1) is the number of the particles α at the position r 1 pointing in the direction ω 1 per unit volume and unit space angle.The number density ρ α ( r 1 ) is obtained by integrating ρ α (1) over all orientations ω 1 .For the homogeneous isotropic system we write simply ρ α for the number density and ρ α (1) = ρ α /4π.For later angular expansions (section 4) it is convenient to extend these definitions also to the ion densities ρ + (1) and ρ − (1).We then have everywhere ρ ± (1) = ρ ± ( r 1 )/4π.Because the O.Z. equation contains two unknown functions, another relation between direct and total correlation functions is needed to calculate the correlation functions.This relation can be written as where b αβ (12) is the so-called bridge term.Usually the bridge term cannot be known exactly and an approximation must be made.If we take b αβ (12) = 0 the hypernetted chain (HNC) approximation is obtained.In the reference hypernetted chain approximation (RHNC) [13], we approximate b αβ in equation (3.2) by the bridge term of a reference system.Here we choose the hard sphere fluid as our reference system, which is obtained by switching off the Coulomb and dipole interactions.The bridge term of the reference system b R αβ (12) can be expressed as where h HS (12) and c HS (12) are the total and direct correlation function of hard spheres and can be calculated using Verlet and Weis's method [14].
Because of the angle dependence of the ion-dipole and dipole-dipole interactions the correlation functions of the ion-dipole and dipole-dipole relations are also angle dependent.To treat such correlation functions we will expand them by the socalled spherical invariants Φ l 1 l 2 l (ω 1 ω 2 ω r ), which are linear combinations of spherical harmonics [15][16][17].We follow the notation of Gray and Gubbins [16] where f = h or c and C(l 1 l 2 l r ; µνλ) is a Clebsch-Gordan coefficient.
In some situations it is convenient to use the " r-frame", which corresponds to choosing the polar axis along the intermolecular axis [16].In r-frame the correlation functions are expanded as ω ′ i are Euler angles of the dipoles relative to r 12 as polar axis.The r-frame coefficients f αβ (l 1 l 2 m, r) are related to the space-fixed frame coefficients f l 1 l 2 l αβ (r) via In the isotropic bulk problem, where the particle density ρ α (1) doesn't depend on position and orientation, the Ornstein-Zernike equation is used usually after Fourier transformation.After spherical harmonic expansion and Fourier transform the Ornstein-Zernike equation can be written in a matrix form where the matrices F (F = H, C) are defined by the k-frame coefficients [15][16][17] of correlation functions, fαβ (l β fαβ (l 1 l 2 m, k). (3.9) Here the indices i = (α, l 1 ), j = (β, l 2 ) are combinations of particle species index and angular index.The k-frame coefficients fαβ (l 1 l 2 m, k) are calculated in the following way from the r-frame coefficients where j l (kr) are the spherical Bessel functions.
In the closure equation (3.2) there are angular dependences in the exponent, therefore we will have difficulty when using the expansion by spherical harmonics.To overcome this difficulty Fries and Patey [6] proposed to differentiate the closure with respect to the particle distance r and differentiating it we get with functions c ′ αβ (12) and η ′ αβ (12) defined as Here we derive the equations in general.The RHNC approximation will be obtained when setting b αβ (12) = b R αβ (12) according to equation (3.3).After harmonic expansion the closure can be expressed with r-frame coefficients as where the factors F l 1 l 2 m l 3 l 4 m 1 ,l 5 l 6 m 2 are the products of Clebsch-Gordan coefficients Later Calliol [18] proposed to differentiate the closure with respect to the orientation of the particles.For example we can differentiate the closure with respect to the orientation of particle 2 by the angular momentum operator and obtain Using the r-frame coefficients the closure can be expressed as with the coefficients G l 1 l 2 m l 3 l 4 m 1 ,l 5 l 6 m 2 defined as The advantage of Caillol's method is that we get a set of equations directly for the expansion coefficients and not for their derivatives as in the procedure introduced by Fries and Patey.But equations (3.18) and (3.19) fail to calculate the coefficient belonging to l 1 = l 2 = 0.In our calculation here we use the combination of the two methods.For the l 1 = l 2 = 0 term we use Fries's method and for the other terms we use Caillol's method.

The method for analysis of the instability
The grand potential functional of an ion-dipole mixture can be written as where F id [ρ] is the Helmholtz free energy functional of the ideal gas  [19,20].
To study the stability of the equilibrium state we consider the fluctuation of the grand potential around the minimum caused by small density fluctuations δρ α (1) [17,21].We use the functional Taylor expansion to the second order The second functional derivatives of Ω are related to the direct correlation functions In order to treat also angular dependent correlation functions, as we need for dipolar interactions, we expand the particle density fluctuation δρ α (1) as well as the direct correlation functions c αβ (1, 2) in spherical harmonics The integral in equation ( 4.3) is a convolution, so we get a product of the corresponding Fourier components after the Fourier transform.Using the expansions (4.5) and (4.6) in equation (4.3) and finishing the angle integration over ω 1 and ω 2 we obtain Equation (4.7) is a quadratic form of the fluctuations vector δρ(k) with δρ µ (k) = . The index µ = (α, l 1 m) is a combination of particle species index and angular indices.The coefficients of the quadratic form yield the matrix M(k).The k-integral can be replaced by a sum over discrete k-values if we imagine periodic boundary conditions: In the case of angle dependent interactions, M(k) is diagonal with respect to the angular index m and the matrix has a block form.Therefore equation (4.8) can be decomposed into a sum over submatrix products for separate Since fluctuation δρ(m, 0) is proportional to V and C is independent of V , the homogeneous density fluctuations δρ 2 α 1 2 decay like V − 1 2 , as they should [23].Using the Ornstein-Zernike equation (3.8) the mean values of the fluctuation products can also be expressed by the total correlation function matrices (4.12) For complete insight into the phase instability of fluid mixtures we must diagonalize the matrix M(m, k) in equation (4.9).With the eigenvalues λ σ (m, k) and the normalized eigenvectors For stable phases all eigenvalues are positive and any fluctuation of a one-particle density will increase the grand potential Ω.If one eigenvalue λ σ (k) approaches zero, then there is a fluctuation δρ ′ σ (m, k) which doesn't increase Ω.In this case the system becomes unstable and the mean value of this δρ ′ σ (mk) becomes infinite.For a pure gas-liquid phase transition the most unstable fluctuation δρ σ (m, k) should be the fluctuation of the total density and for a purely demixing phase transition the most unstable fluctuation should be a concentration fluctuation.In general the most unstable fluctuation δρ ′ σ (m, k) (equation (4.14)) belonging to the smallest eigenvector λ σ , which eventually goes to zero,is a linear combination of all kinds of density fluctuations.The phase instability is characterized uniquely [17,21] by the eigenvector x σ (m, k) of the smallest eigenvalue.The essence of the derivation is that after expansion into spherical harmonics, all the expansion coefficients can be considered as independent fluctuating densities yielding a mixture with very many "components".
We now specialize more towards our case of the ion-dipole mixture.The dimension of the matrices M(m, k) depends on the truncation of the expansion in equations (4.5) and (4.6).If we trucate the expansion with l 1 , l 2 3, then M(0, k) is a 6 × 6 matrix.We introduce two new fluctuation variables which describe respectively the total ionic density fluctuation and the charge fluctuation.Using the new fluctuation variables the coefficient matrix of δΩ(0, k) becomes of block form of two 3 × 3 matrices.The total ionic density fluctuation δ ρion (k) couples only with δ ρ00 d (k) and δ ρ20 d (k).The charge fluctuation δ ρq (k) couples only with δ ρ10 d (k) and δ ρ30 d (k).For the coefficient matrix related to δ ρion (k), δ ρ00 d (k) and δ ρ20 d (k) we find the first vanishing eigenvalue at k = 0, which indicates an instability with respect to homogeneous density fluctuations.In this case, δ ρion (0) is only coupled with δ ρ00 d (0) by a 2 × 2 matrix because the matrix elements M ion,(d20) (0, k) and M (d00),(d20) (0, k) are proportional to k 2 when k → 0. This 2 ×2 matrix later contains the smallest and finally vanishing eigenvalue λ 1 .It is convenient to define the fluctuations of total density and ionic concentration where ρ = ρ + + ρ − + ρ d , c ion = c + + c − with c i = ρ i /ρ and δ ρd (0) = (4π) 1 2 δ ρ00 d (0).Using the new fluctuation variables the grand potential variations with respect to fluctuations δ ρion (0) and δ ρd (0) can be rewritten as where the coefficient matrix M(0) is symmetric and defined as Kirkwood and Buff have related the total correlation functions to thermodynamic functions.With the help of the O.Z. equation their relations can be written into the form (ρ α ρ β ) where αβ (0).Using the thermodynamic relations with the difference of the partial volumes ∆ = ρ[(ν After diagonalization of M(0) with the eigenvalues λ i and the related eigenvectors x i , the δΩ(m = 0, k = 0) can be written as a sum of pure squares where The smallest eigenvalue λ 1 decides the stability of the phase.The eigenvector x 1 determines the softest fluctuation mode δ ρ′ 1 (0) = x 1,1 δ ρ(0) + x 1,2 δc(0) which has the smallest restoring force (proportional to λ 1 ).Therefore x 1 characterizes properly the phase instability [17,21].The border of a stability region is indicated by λ 1 going to zero (from the positive side).In the case λ 1 = 0 there is also Det There are two cases where Det M(0) = 0.One is κ T −1 = 0.In this case M(0) is diagonal and the softest fluctuation mode is just the total density fluctuation and the phase instability is pure condensation.Another case is ∂ 2 G/∂c 2 ion | T,P,N = 0.Here the softest fluctuation mode δ ρ′ 2 ∆δ ρ(0) is generally a combination of total density and concentration fluctuation.How strong the total density will change in the softest fluctuation mode depends on the partial volume difference ∆.If ∆ = 0, the softest fluctuation mode is just a concentration fluctuation and the phase instability is pure demixing.When |∆| is large, the total density will fluctuate also.The sign of ∆ decides if the total density will increase or decrease.We have seen a case [21] where ∂ 2 G/∂c 2 = 0 and κ −1 T > 0 mark an instability with respect to essentially pure total density fluctuations with negligible concentration change.Therefore there exist instabilities, where ∆ is so large that ∂ 2 G/∂c 2 = 0 indicates a condensation instability and not a demixing.
In the following section we evaluate M(0) (equations (4.18)-(4.20))from the calculated correlation functions and find a demixing instability in our electrolyte model.

Results for the hard sphere ion-dipole electrolyte.
We now present our calculations for the model electrolyte.The total density is fixed at the high liquid value ρ * = ρ • σ 3 = 0.8 and the ion concentration is changed by exchange of ions versus dipoles.The dipolar interaction strength is chosen as µ * 2 = µ 2 /(σ 3 k B T ) = 2.5 because below µ * 2 = 2.25 the interesting missibility gap at very low ion charges and for neutral solutes does not appear [12].The angular expansions in equation (3.4) drastically increase the dimension of the correlation function calculation (the number of independent functions).We limit ourselves to l 1 , l 2 3, which leads to 23 unknown functions to be determined.l takes then values up to 6 due to l l 1 + l 2 .We proceed by fixing the ionic interaction q * 2 = q 2 /(σk B T ) and varying the ionic concentration c ion = (ρ + + ρ − )/ρ.For each system, we evaluate the correlation functions according to section 2, calculate the matrix M(k = 0) (equations (4.18)-(4.20))and determine the smallest of its eigenvalues λ 1 (equation (4.31)) and the associated eigenvector x 1 .We have checked occasionally that λ 1 (k = 0) < λ i (k) for all k > 0.
In figure 1 we plot the smallest eigenvalue λ 1 versus c ion for growing q * 2 .We find two regions where the solution procedure becomes unstable and where λ 1 strongly indicates an abrupt decrease to zero.We roughly extrapolate to zero and draw the "phase-diagram" in figure 2, where the lines mean spinodales separating stable or metastable regions from unstable systems.For two values of q * 2 we also give information about the eigenvector x 1 representing the most unstable fluctuation.We show in figure 3 the angle which this eigenvector forms in a δρ−δc-coordinate system with the positive δc axis [17,21].The angles are small and |α| < 10 • .This means that we find predominantly demixing instabilities.
The positive values of α for the small ion interaction q * 2 = 0.2 indicate that an increase of ion concentration is accompanied by a decrease of density for the most unstable fluctuation with the smallest restoring force.This shows that the attraction between dipoles is the driving force for the instability.We have proposed in [12] that favorable dipole structures are possible for certain concentrations of neutral solutes and that this is the reason for the missibility gap.It was later pointed out to us that indications of this mechanism had been seen earlier in MC-simulations for a Stockmeyer solvent (Lenard Jones plus dipole interactions) [24].
The negative values of α for the high ionic interaction q * 2 = 60 tell us that together with an increase of ionic concentration goes an increase in density, when the rise of the free energy should be minimal.Therefore the instability is driven by the attraction of the ions.We believe that this instability is related to the "condensation" of the purely Coulombic interacting fluid found in the restricted primitive model (RPM) (see discussion).
As explained in section 2, the fluctuations of ionic concentration and of net charge density are decoupled.When the concentration fluctuations become unstable, the charge fluctuations are still very strongly constrained, i.e. the related eigenvalues of the matrix M(k) are very large.Homogeneous charge fluctuations (k = 0) are not at all possible because the related eigenvalue is ∞.This can be seen from M q,q (k) = 1−ρ + c++ (k)+ρ − c+− (k) and equation (A.8) in the appendix.At large distances c ij (r) is proportional to the Coulomb potential and therefore its Fourier transform c ij (k) at small k goes like k −2 as shown in equation (A.8).Then M q,q (k) ∼ k −2 goes to infinity when k → 0 which suppresses homogeneous charge fluctuations completely.The total electrical neutrality of the system must be exactly satisfied.Figure 3.The angle α (in degrees) between the eigenvector x 1 and the direction of pure concentration fluctuations for q * 2 = 0.2 and q * 2 = 60.

Discussion and conclusions
We are aware of several objections to the analysis of phase equilibria and stability via correlation functions calculated from integral equations.Lovett has discussed that finding the correct solution from integral equations is a "marginal" event and that approximations to the kernel, which are always necessary, should lead into the realm of no solution [25].Belloni has demonstrated [26] that with pure HNC equations one always meets a point of bifurcation of two solutions and no solutions beyond, before the spinodal is reached where fluctuations or compressibilities become infinite or in our language the smallest eigenvalue of the matrix M (equation (4.8)) goes to zero.But on the other hand experience shows that usually one gets very resonable solutions from HNC or related equations with all kind of right physical properties; and more important there is no case known to us, where the indication of a phase transition by a sudden growth of fluctuations and related impossibility of finding solutions of the integral equations by iteration was not related to a real instability of the system.This can be seen for the cases discussed by Belloni [26].For the dipolar system we found the missibility gap with neutral solutes [12,17], the lower case in figure 1, which had been seen in simulations for dipoles and neutrals with additional Lenard-Jones interactions.The phase transition was also caused there by the dipole forces [24].
In recent years we also investigated the phase behaviour of the pure hard sphere dipole fluid.From the instability of fluctuations in the low density dipole gas we predicted the formation of aligned clusters or chains of dipoles [27].Later on, this phase was seen in simulations [28,29].At higher densities all the characteristics of the unstable fluctuations could be related [30] to phases seen in simulations at lower temperature [31][32][33] including the transition to a ferroelectric state.This transition was of first order and the coexistence lines well framed the instability line which again confirmed the interpretation of the no-solution-line as an indication of a spinodal.
Therefore there are good reasons to expect that in the near future there will be a simulation showing the phase separation for our model electrolyte.The phase instability of the ion-dipole mixture has also been found in a study using the mean spherical approximation (MSA) for the calculation of the correlation functions [3].This approximation underestimates the interactions and correlations and therefore finds the instability at much higher q * 2 values or lower temperatures.It is also true that the line separating a region of no solutions from one with physically reasonable solutions of the integral equations is not simply related to the real spinodal of the system (see the discussion in 25, also [34]).Therefore the critical point or the critical behaviour cannot be safely determined via integral equations.So we consider as the main result of our investigation, that the hard sphere ion dipole mixture is unstable with respect to demixing, when the ion interactions get strong enough.We expect that for a solvent with µ * 2 = 2.5 this instability is met before the ions have one full charge.Recently there was intensive research on ionic solutions, which indeed do show a demixing phase transition [35].With these electrolytes, one discusses the alternative of "hydrophobic" demixing or demixing due to condensation of the gas of ions towards a liquid.We have discussed in relation to figure 3 that the most unstable fluctuation indicates that the attractive forces between the ions lead to the instability, which points to the second alternative.The mechanism is then the same as in the pure Coulombic RPM of hard sphere ions.Belloni [26] shows the instability region as predicted by HNC at very low densities and q * 2 above 7.35 ("critical point" ρ * = 0.3 • 10 −3 , q * 2 = 7.35 or T * = q * −2 = 0.136).The ion dipole mixture in our RHNC treatment has the minimum of the instability curve (figure 2) near ρ * ion = 0.4 and q * 2 = 56.When the dipole solvent is considered as a medium with dielectric constant ǫ = 8, the effective interaction should be taken as q * 2 = q 2 /(ǫσk B T ) = 7 at the minimum.The simulations of the RPM yield the critical point at ρ * = 0.03, q * 2 = 17.2 (T * = 0.058) [36].
Compared to the RPM our phase transition is found at higher densities and at higher temperatures indicating that the solvation of the ions by the solvent dipoles changes the phase transition strongly.As already seen for neutral solvents, with the ions we can also expect that special dipole configurations which require certain concentrations of ions and dipoles will play a role in the determination of the phase diagram.

. 9 )
Now the index of vector and matrix is only a combination of particle species index and angular index l 1 .The matrices M(m, k) are Hermitian and can be expressed by the direct correlation function matrices M(m, k) = I − C(m, k) according to equation (3.9).They become real, if in equation (4.6) there exist only terms with even l-values.This is the case for pure dipolar systems or for mixtures of dipoles and neutral particles, but not for mixtures of dipoles and ions.The probability distribution P (δρ(m, k)) for fluctuations is GaussianP (δρ(m, k)) ∼ exp − δΩ(m, k) k B T = exp − 1 2V δρ(m, k)M(m, k)δρ † (m, k)(4.10) and the mean values of the fluctuation products are

1 Figure 1 .
Figure 1.The smallest eigenvalue λ 1 as function of ionic concentration for different values of q * 2 .

Figure 2 .
Figure 2. Regions of instablity and "spinodal" curves for the ion dipole mixture.
1/2the thermal de Broglie wavelength.F ex [ρ] is the excess Helmholtz free energy functional, which represents the contribution from the interactions between the particles.The sum over α in equation (4.1) is done for plus ion, minus ion and dipole.W α (1) = W α (1) − µ α with chemical potential µ α and external potential W (1)).The equilibrium density ρ α (1) is determined by the minimum of the grand potential functional for fixed T, V, µ α