Ionic fluids with association in dipoles and quadrupoles

Accepted Manuscript Ionic fluids with association in dipoles and quadrupoles Paulo Sérgio Kuhn

PII: DOI: Reference:

S0378-4371(18)30527-2 https://doi.org/10.1016/j.physa.2018.04.099 PHYSA 19531

To appear in:

Physica A

Received date : 29 August 2017 Revised date : 15 March 2018 Please cite this article as: P.S. Kuhn, Ionic fluids with association in dipoles and quadrupoles, Physica A (2018), https://doi.org/10.1016/j.physa.2018.04.099 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Highlights (for review)

HIGHLIGHTS • RPM fluid with quadrupole association is studied. • The Helmholtz free energy for quadrupole-ion interaction is constructed. • Phase diagrams are built and critical constants are estimated.

*Manuscript Click here to view linked References

Ionic Fluids with Association in Dipoles and Quadrupoles Paulo S´ergio Kuhn [email protected] Departamento de F´ısica Instituto de F´ısica e Matem´atica, UFPel Campus Cap˜ao do Le˜ao - Pr´edio 05 - 3o. andar CP 354 - 96001-970 - 55(53)3275-7343 Pelotas - RS - Brazil March 15, 2018 Abstract The Debye-H¨ uckel theory for an electrolyte solution is considered beyond dipole formation by including neutral clusters with four particles, through association between two positive and two negative ions. The Helmholtz free energy for the interaction between clusters and the ionic solution is evaluated according to Debye-H¨ uckel theory. A simple association constant for the quadrupole complex is used, considering the electrostatic interaction between the ions in a fixed configuration. Hard-core contributions with Carnahan-Starling and free-volume theories are also included. The complete free energy shows phase separation, with a critical point in the region of numerical simulations.

Keywords: 64.70.F- Liquid-vapor transitions 64.70.qd Thermodynamics and statistical mechanics 64.75.Gh Phase separation and segregation in model systems (hard spheres, Lennard-Jones, etc.)

1

Introduction

Ionic solutions are important systems in chemistry, physics and biology [1, 2]. The simplest theoretical description of these systems is the hard-sphere ionic 1

fluid, where an equal number of positive and negative charges are located at the centers of spheres with the same diameter σ, in a continuous medium of permitivitty ε. The model is named the restricted primitive model of electrolytes, or RPM fluid. A central step in the understanding of these systems is the Debye-H¨ uckel theory [3-6], which leads to a free energy that, in spite of its simplicity, shows the existence of a critical point and hence of a phase diagram, with or without hard-core contributions [7-10]. When dipole formation [11-17] is considered and the corresponding free energy is incorporated in the theory, the location of the critical point agrees better with the numerical simulation estimatives [18, 19]. The theory has been supplemented with hard-core contributions. In [20] we have chemical association into neutral dipolar pairs also considered. Another important theoretical approach to the problem of ionic fluids is by solving integral equations that lead to correlation functions [4-6]. The mean spherical approximation (MSA) [21-24] is a successful case in which the defining integral equations can be solved exactly for the RPM model. This is of extreme value, since exact solutions are rare in physics. In the limit of vanishing charge the theory reduces to the Percus-Yevick equation for hard spheres [5, 6]. An extensive analysis of the results of MSA for the RPM model is found in [19]. The reduced critical density is somewhat below simulation estimatives, while the critical temperature is located above them. This fact seems to be present in all MSA theories. This also occurs when hard-core terms are included [19]. The problem of ionic solutions has been studied through theory, experiments, and numerical simulations by many authors [25-58], hence we do not claim to do a full citation of all the works in the area. The best we can do is to read the works of M. E. Fisher and collaborators, some of them are listed in the references. In the last years, however, the location of the critical point estimated in numerical simulations has changed [59, 60]. The critical density particularly increased significantly from ∼ = 0.02 − 0.03 to ∼ = 0.06 − 0.08 (the total number density of particles). In this work we consider Debye-H¨ uckel theory including association of ions in neutral clusters of four ions, forming a quadrupole. The corresponding free energie for quadrupole-ion interaction is evaluated, and when the resulting free energy is supplemented with hard-core contributions, a critical point in the region of numerical simulation values is obtained. The inclusion of neutral clusters with four particles, corresponding to quadrupoles, is by no means a full description of the association phenomena that occur in the system. We may have of course complexes with three, four, five, etc. ions associated, hence complexes of various sizes can be formed, with many of them not neutral. However, the task of solving the resulting equilibrium equations, the mass law equations, is much more difficult, since 2

we have more equations to solve. Another possible approach would be to consider a cluster distribution of sizes, as for example in the polyelectrolyte case and other systems [61]. The problem is how to evaluate the free energy in this case for all kinds of clusters, which is an interesting question that may be investigated. We must also observe that we consider here only the interaction of clusters, dipoles and quadrupoles, with the ionic solution. Hence the interactions between dipoles among themselves, or quadrupoles, or between dipoles and quadrupoles, are not included in our approximation. This possibly may be done, at least in principle, but it is certainly not an easy task, therefore we have not tried to include these in the present approximation. The presence of clusters of several sizes has been observed in numerical simulations [62, 63]. The structure of this work is as follows. Section 2 presents succinctly the restricted primitive model for electrolyte solutions, and the corresponding Debye-H¨ uckel theory. Section 3 shows the extension of the theory by the inclusion of dipoles formed by association. Section 4 presents our evaluation of the Helmholtz free energy that accounts for the interaction between the ionic solution and the quadrupoles formed by association. This is done based on Debye-H¨ uckel ideas, evaluating an electrostatic interaction energy and then obtaining the free energy through the charging process. In this section we also state the approximation that we make for the association constant, or internal partition function, corresponding to a quadrupole. Section 5 presents the well-known free energies of Carnahan-Starling and free-volume approximation, which describe the hard-core interaction. In section 6 the total free energy of the system is stated and the equations that lead to the phase diagram are obtained. Finally, section 7 presents the results and the conclusions.

2

The RPM model and Debye-H¨ uckel theory

The Debye-H¨ uckel theory for the restricted primitive model of electrolytes (RPM) is very important in the physics and chemistry of charged solutions [3-19]. We make here a concise exposition, mainly to clarify the notation, and also because this is the starting point for the improvements of association, which is our main concern. The RPM model for electrolytes consists of Ns ionic pairs, hence there are 2Ns particles, Ns cationic ions and Ns anionic ions. The ions are hard spheres of diameter σ, and the medium is continuous with permitivitty ε. The charges are located at the centers of the spheres and the system is globally neutral. The Poisson equation gives the electrostatic potential around a given ion, 3

∇2 φ = 0 , r ≤ σ, ρq 2 ∇ φ=− , r > σ, ε

(1)

where ρq is the charge distribution and σ, besides being the diameter of an ion, is the radius of an exclusion sphere around any ion (fig. 1).

σ σ

Figure 1: Exclusion sphere of radius σ around a given ion. The Debye-H¨ uckel theory begins by approximating the charge distribution by a Boltzmann factor, X ρq = qj ρj exp{−βqj φ} , r > σ, (2) j

with β = 1/kT , k is the Boltzmann constant, and the sum is over the species, namely positive and negative ions. The Poisson equation is then named the Poisson-Boltzmann equation. The second approximation in Debye-H¨ uckel theory is the linearization of the charge distribution in r > σ, ρq ≈ qρ+ (1 − βqφ) − qρ− (1 + βqφ) = −βq 2 φρ1 , where q is the proton charge, and ρ1 = ρ+ + ρ− is the density of free ions. Without association we have ρ1 = 2ρs , but this is changed when dipoles and others clusters are formed. We have then the linearized Poisson-Boltzmann equation, ∇2 φ = 0 , ∇2 φ = κ 2 φ ,

r ≤ σ, r > σ,

(3)

for the electrostatic potential around an ion. The Debye screening length κ−1 is defined by 4

κ2 = 4πλB ρ1 ,

(4)

where λB = βq 2 /4πε is the Bjerrum length. The solution of equation (3) is qi κqi − , r < σ, 4πεr 4πε(1 + κσ) qi e−κ(r−σ) φ(r) = , r > σ, 4πεr(1 + κσ)

φ(r) =

(5)

with qi denoting the charge of the central ion. The second term in the potential for r < σ is the potential on the ion due to the ionic solution. This term represents the interaction between ions and it will be used in the charging process, in order to find the Helmholtz free energy. The result is [19] βFdh = −

V ω1 (κσ) , 4πσ 3

(6)

with ω1 (x) ≡ ln(1 + x) − x +

x2 , 2

(7)

and x = κσ.

3

Dipoles in solution

Dipole formation was considered by Levin and Fisher [18, 19], following Debye-H¨ uckel and Bjerrum theories. In this section we state the corresponding free energy they obtain for completeness, and also because the formation of other clusters may be studied in the same lines, as we make in the next sections. The electrostatic potential around a point dipole that satisfies the linearized Poisson-Boltzmann equation, ∇2 φ = 0 , ∇2 φ = κ 2 φ , is

5

r ≤ a2 , r > a2 ,

(8)



 A φ(r, θ) = + Br cos θ , r ≤ a2 , r2 φ(r, θ) = C k1 (κr) cos θ , r ≥ a2 ,

(9)

with the modified spherical Bessel function of third kind of order one [65], e−x (1 + x) . x2 The radius a2 of the exclusion sphere around a dipole (fig. 2) is evaluated through an angle average [18, 19], which gives a2 = 1.1619σ. k1 (x) =

a2 p

Figure 2: Exclusion sphere of radius a2 around a central point dipole. The constants in (9) are obtained in the usual way, by requiring continuity of the potential and its normal derivative at r = a2 , p , 4πε p y2 B=− , 4πεa32 3 + 3y + y 2 3p κ2 ey C= , 4πε 3 + 3y + y 2

A=

(10)

with p = qσ and y = κa2 . The second term in the potential for r < a2 , Brcos θ, is the potential on the dipole due to the ionic solution. This term represents the interaction between dipoles and ions, and it will be used in the charging process that leads to the Helmholtz free energy. The result is

6

βFdi = −N2

σ3 2 y ω2 (y) , T ∗ a32

(11)

with T ∗ = σ/λB and   y2 3 2 . ω2 (y) = 4 ln(1 + y + y /3) − y + y 6

(12)

We have also the association constant for dipoles, which is given by [1, 11, 18, 19]

ζ2 /σ

3

  2π e2 4π 1/T ∗ 1 1 = − 1+ e + 3 T ∗3 3 2T ∗ 2T ∗ 2 2π 1 − [Ei(2) − Ei(1/T ∗ )] , 3 T ∗3

(13)

for T ∗ < 0.5 according to Bjerrum’s choice. The corresponding free energy of association for N2 dipoles in solution is then F2 = −N2 kT ln(ζ2 /σ 3 ) .

4

(14)

Quadrupoles in solution

We consider a quadrupole formed by association as the complex in fig. 3, with s4 ≥ 1 a parameter to be determined. First, if the plane of the figure is the xy plane and the origin is at the symmetry center we have the quadrupole tensor [64]   0 −3qs24 σ 2 0 0 0 . Qjk =  −3qs24 σ 2 (15) 0 0 0

The electrostatic potential for an isolated quadrupole is φ0 =

1 1X lj lk Qjk , 4πεr3 2 jk

(16)

where the l′ s are direction cosines for vector r, and j, k = x, y, z. The result for our quadrupole is 3qs24 σ 2 sin2 θ sin2ϕ , 8πεr3 with the usual angles of spherical coordinates. φ0 = −

7

(17)

1 2

+

− d

d s4 σ

O d

d

−

+ 4

3 s4 σ

Figure 3: A quadrupole formed by association of two positive ions and two negative ions.

4.1

The association constant

The problem of evaluating the association constant for a given cluster is an old problem in the statistical mechanics of ionic solutions and other systems [19]. We make here the simple approximation of taking the electrostatic energy U4 of a given configuration of charges, namely that shown in fig. 3, U4 = √

X

Uij =

i
√ q2 (−4 + 2) , 4πεs4 σ

(18)

where we use d = s4 σ/ 2. Hence, βU4 =

√ λB (−4 + 2) . s4 σ

(19)

If we write the internal partition function of one cluster as ζ4 /σ 3 = exp{−βU4 }, the corresponding free energy for N4 quadrupoles is F4 = −N4 kT ln(ζ4 /σ 3 ) = N4 U4 , or √ λB λB (−4 + 2) = −N4 (2.58) . (20) s4 σ s4 σ This is our approximation for the free energy of association in quadrupoles. βF4 = N4

8

4.2

Helmholtz free energy for quadrupole-ion interaction

Here we evaluate the Helmholtz free energy for quadrupole-ion interaction according to Debye-H¨ uckel theory. The electrostatic potential for a point quadrupole in an ionic solution satisfy ∇2 φ = 0 , ∇2 φ = κ 2 φ ,

r ≤ a4 , r > a4 ,

(21)

where a4 is the radius of an exclusion sphere around a quadrupole. The diameter of a sphere with the same volume of a quadrupole is σ4 = 41/3 σ ∼ = 1.5874σ. Hence the exclusion sphere around a quadrupole has radius σ4 + σ σ = (1 + 41/3 ) ∼ = 1.2937σ . 2 2 The solution for the potential is a4 =

 B r < a4 , φ = Ar + 3 sin2 θ sin2ϕ , r φ = Ck2 (κr) sin2 θ sin2ϕ , r > a4 , 

2

(22)

where k2 is the modified spherical Bessel function of third kind of order two [65], e−x (3 + 3x + x2 ) . x3 The constants are evaluated through usual boundary conditions of continuity of potential and its normal derivative at the exclusion sphere, k2 (x) =

3qs24 σ 2 (1 + v)v 2 , 8πεa54 15 + 15v + 6v 2 + v 3 3qs24 σ 2 B=− , 8πε 5ev v 3 3qs24 σ 2 , C=− 8πεa34 15 + 15v + 6v 2 + v 3 A=

(23)

with v = κa4 . We may verify that both solutions goes to the potential φ0 of an isolated quadrupole as κ → 0, eq. (17). 9

The Helmholtz free energy for the quadrupole-ion interaction is evaluated by the Debye charging process of the electrostatic energy. The potential at the quadrupole due to the ionic solution, and that will be used in the charging process is ψ = A r2 sin2 θ sin2ϕ .

(24)

Note that this potential goes to zero if κ → 0. The potential energy of a quadrupole in an electric field is   ∂Ej 1X Qjk U =− , j, k = x, y, z. (25) 6 jk ∂k r=0

The derivatives are easily evaluated if we note that

ψ = A r2 sin2 θ sin2ϕ = 2 A x y .

(26)

The electric field corresponding to ψ is then E = −∇ψ = −2 A y ˆi − 2 A x ˆj ,

(27)

and the electrostatic energy of the quadrupole in this field is, using (25), 2 3q 2 s44 σ 4 (1 + v)v 2 U = Qxy A = − . 3 4πεa54 15 + 15v + 6v 2 + v 3

(28)

The charging process gives the free energy for quadrupole-ion interaction, if there are N4 quadrupoles in solution, as Z 1 3λB s44 σ 4 I(v) 2 A(λ)dλ = −N4 , (29) βFqi = N4 βQxy 3 a54 v2 0 with

v2 − 5v + 5 ln(v 3 + 6v 2 + 15v + 15) − 5 ln 15 . (30) 2 We have already considered quadrupoles in ionic solutions [66], but not in the context of association. I(v) =

5

The effect of the finite volume of ions

We have considered the dimension of particles in finding the electrostatic potentials around ions, dipoles and quadrupoles. This leads to exclusion radius around each type of particle. Here we consider the free energy corresponding 10

to this finite volume. Two forms will be used, the Carnahan-Starling free energy [67], and the free volume approximation as in van der Waals theory [4]. Since these are well-known theories, we just state the expressions. We have the Carnahan-Starling free energy, βFcs = with

X

Ni

 y (4 − 3y ) cs cs , (1 − ycs )2

(31)

π 3 X  σ ¯ ρi . (32) 6 The sums in (31) and in ycs extend over all species, ions, dipoles and quadrupoles. The mean diameter σ ¯ is defined as P Ni σi σ ¯= P . (33) Ni ycs =

The diameter of positive and negative ions is σ. For dipoles and quadrupoles we take the diameter of spheres of the same volume, hence σ2 = 21/3 σ and σ4 = 41/3 σ. In the free-volume approximation we have X    X βFhc = − Ni ln 1 − B i ρi , (34) where the sums extend over all species, and √ Bj /σj3 = 1/ 2 ,

√ 4/3 3 ,

1,

2π/3 ,

(35)

correspond, respectively, to fcc, bcc, and sc packings, and to the exact second virial coefficient.

6

Coexistence and phase diagrams

The system has Ns ionic pairs, hence a total of 2Ns ions, half positive and half negative. The whole system is neutral. At equilibrium we have association, and dipoles and quadrupoles are formed. We define in the usual way the degree of dissociation α, ρ+ = αρs ,

ρ− = αρs ,

(36)

where ρ+ is the number density of free positive ions, ρ− the density of free negative ions, and ρs = Ns /V is the density of ionic pairs. Note that the

11

0.065

0.06

T* 0.055

0.05

0.045 0

0.005

0.01

0.015

0.02 ρ s*

0.025

0.03

0.035

0.04

Figure 4: Phase diagrams for RPM fluid with monopoles only (see reference [19]): without any hard-core terms (solid line), with Carnahan-Starling free energy (dashed line); with free-volume free energy and Bj∗ = 2π/3 (dot√ dashed line) Bj∗ = 1 (dot-dot-dashed line), Bj∗ = 4/3 3 (dotted line), and √ Bj∗ = 1/ 2 (dot-dashed-dashed line). Note that ρs is half the total number density of particles. total density of particles is 2ρs . The fraction of particles forming dipoles is denoted by m2 , and the fraction of particles in quadrupoles by m4 , ρ2 2ρ4 , m4 = , (37) ρs ρs where ρ2 is the number density of dipoles, and ρ4 of quadrupoles. We have then the normalization condition α + m2 + m4 = 1. The densities satisfy m2 =

ρ+ = ρs − ρ2 − 2ρ4 ,

ρ− = ρ+ .

(38)

The inverse Debye screening length κ−1 is defined in the usual way by κ2 = 4πλB (ρ+ + ρ− ) .

(39)

The complete Helmholtz free energy for the system is then F = Fid + F2 + F4 + Fdh + Fdi + Fqi ,

(40)

supplemented with Fhc or Fcs , where Fid is the ideal gas free energy. The free energy is a function of the densities, 12

0.065

0.06

T* 0.055

0.05

0.045 0

1

2

3

κσ

4

5

Figure 5: Diagram of T ∗ × κσ for RPM fluid with monopoles only (see reference [19]); the curves are labeled as in figure 4. 0.055

0.055

0.05

0.05

0.045

0.045

T* ρ s*

0.04 0.035

0

0.05

0.1

0.15

0.2

α

0.04 0.035 0 0.25

0.055

0.055

0.05

0.05

0.045

0.045

0.1

0.2

0.3

0.4

T* m2

0.04 0.035 0.6

0.7

0.8

0.9

κσ

0.04 1

0.035

0

2

4

6

8

10

Figure 6: Phase diagrams for RPM fluid with monopoles and dipoles (see reference [19]); the curves are labeled as in figures 4 and 5. Note that ρs is half the total number density of particles.

F = F (ρ+ , ρ− , ρ2 , ρ4 ) .

(41)

Minimization of the free energy leads to the mass action law and to equations 13

0.065

0.055

(b)

(a) 0.05

0.06

T*

0.045 0.055 0.04 0.05

0

1

2

3

4

(c)

0.055

0

5

0.05

0.05

0.045

0.045

0.04

0.04 0

0.05

0.1

0.15

ρ s*

0.2

0.08

0.12

0.25

0.16

(d)

0.055

T*

0.035

0.04

0

0.1

0.2

ρs*

0.3

Figure 7: Phase diagrams for RPM fluid with monopoles, dipoles and quadrupoles. (a) without any hard-core terms; (b) with Carnahan-Starling free energy; (c) with free-volume free energy for Bj∗ = 2π/3 (solid line); for √ √ Bj∗ = 1 (dashed line); (d) for Bj∗ = 4/3 3 (solid line), and Bj∗ = 1/ 2 (dashed line). The values of s4 are given in tables 1 and 2. Note that ρs is half the total number density of particles. for ρ2 and ρ4 , µ2 = 2µ+ ,

µ4 = 4µ+ ,

(42)

where the µ’s are usual chemical potentials. Solving the above equations we have ρ2 and ρ4 , and then α, m2 and m4 . In order to find the phase diagram, we have to solve the above equations in each phase for a given temperature, together with equality of pressure and Gibbs free energy. This amounts to six equations, µ2v = 2µ+v , µ2l = 2µ+l , µ4v = 4µ+v , µ4l = 4µ+l , pv = pl , Gv = Gl ,

(43)

where v and l denotes vapour and liquid phases, respectively. The variables are ρs , α, m2 in each phase, performing six variables. Once the solution is found we have m4 , the others densities, and all variables we are interested 14

0.06

0.06

0.054

0.054

T* 0.048

0.048

α

0.042

m2

0.042 0

0.1

0.05

0

0.15

0.06

0.06

0.054

0.054

0.02

0.04

κσ

T* 0.048 0.042

0.06

0.048

m4

0.8

0.042 0.85

0.9

0.95

1

0

2

4

6

8

Figure 8: Values of α, m2 , m4 and κσ, corresponding to the diagrams in fig. 7: without any hard-core terms (solid line), with Carnahan-Starling free energy (dashed line); with free-volume free energy and Bj∗ = 2π/3 (dot-dashed line) √ √ Bj∗ = 1 (dot-dot-dashed line), Bj∗ = 4/3 3 (dotted line), and Bj∗ = 1/ 2 (dot-dashed-dashed line). in. Fortunately we have succeded in solving this system, using Octave in a Linux system. The results are discussed below. System − CS 2π/3 1√ 4/(3√ 3) 1/ 2 2π/3 (FL)

s4 1.02 1.05 1.05 1.04 1.04 1.04 −

Tc∗ 0.0619 0.0534 0.0511 0.0564 0.0577 0.0581 0.0522

ρ∗sc 0.0606 0.0396 0.0422 0.0424 0.0414 0.0421 0.0121

αc 0.0897 0.0346 0.0267 0.0455 0.0524 0.0572 0.145

m2c 0.0663 0.0496 0.0362 0.0559 0.0656 0.0687 0.855

m4c 0.844 0.916 0.937 0.898 0.882 0.874 −

Table 1: Critical values for the RPM fluid with dipoles and quadrupoles, for some values of s4 . The first line is without any hard-core terms. The several hard-core contributions are indicated as CS for Carnahan-Starling, and the free-volume terms are indicated by the respective values of Bj /σ 3 . The values obtained by Fisher and Levin [18, 19] with Bj /σ 3 = 2π/3 are in the last line. Note that ρs is half the total number density of particles.

15

System − CS 2π/3 1√ 4/(3√ 3) 1/ 2 2π/3 (FL)

s4 1.02 1.05 1.05 1.04 1.04 1.04 −

κc σ 1.486 0.803 0.744 0.927 0.972 1.021 0.922

p∗c 0.0162 0.0188 0.0212 0.0174 0.0163 0.0159 0.0068

Zc 0.134 0.238 0.251 0.205 0.196 0.190 0.282

ρ2c /ρ1c 0.369 0.717 0.678 0.614 0.626 0.601 2.940

ρ4c /ρ1c 2.352 6.617 8.774 4.937 4.208 3.820 −

Table 2: Other variables at the critical point. The systems are denoted as in table 1.

7

Results and Conclusions

The numerical simulations for the restricted primitive model locate the reduced critical temperature at Tc∗ ∼ = 0.049 and ρ∗s,c ∼ = 0.031 − 0.04 [59]. Note that reference [59] uses the total number density of particles, which in our notation is 2ρ∗s . We construct below the phase diagrams for the system, including different terms for the hard-core interaction. The results depend on the choice of the parameter s4 . The estimative of critical values obtained from the phase diagrams are collected in tables 1 and 2, including all systems. For comparison we reproduce the results for RPM fluid with monopoles only in figures 4 and 5, and with monopoles and dipoles in figure 6. This is done only for completeness, since the full analysis is made in [18] and [19]. The phase diagrams T ∗ -ρ∗s for the system with monopoles, dipoles and quadrupoles are shown in figure 7, and the corresponding diagrams for α, m2 , m4 and κσ are shown in figure 8. The reduced pressure is p∗ = βpσ 3 and the compressibility factor is Z = βp/(2ρs ). Figures 7 and 8 represent our main results. The values in tables 1 and 2 show a great relative number of quadrupole clusters, and the critical temperature and density are reasonable, according to simulation values, with our choice for the parameter s4 . Note that Fisher and Levin find a great number of particles associated in dipoles (last line of the tables). We remark that in our analysis the only kinds of clusters are dipoles and quadrupoles, and that the parameter s4 is used as a fitting parameter. Although in ionic solutions several different clusters are formed, as shown in simulations [62, 63], we expect that by considering only two kind of clusters, dipoles and quadrupoles, we have at least an approximation to the problem. We see that free ions are reduced in favor of association, and the critical density found here is increased relatively to the system without quadrupoles. Of course considering two kinds of clusters does not shows all the richness 16

of the complex formation process that numerical simulations have displayed. A more general theory would consider other clusters, including not neutral, but the mathematical difficulties involved obviously increase, since for one more kind of cluster we have two more equations in the equilibrium conditions as (43), in the task of building the phase diagram. The configuration of the system has one more equation of equilibrium to solve in (42). The full generalization includes the search for a size distribution, as stated in the introduction. This possibility is an interesting theoretical problem, although the inherent practical difficulties. A number of other refinements may be done, as the exact, or more precise, evaluation of association constant for quadrupoles. For dipoles there is a great number of works that account for precise values of the association constant [11-17]. For quadrupoles this is an interesting question that may be explored as a further improvement. The exclusion sphere around clusters is also a problem that may be addressed in a more precise way. For dipoles we have a precise formulation of a2 as an angular average [18, 19]. The problem for quadrupoles is evidently much more complex, hence we used the above approximations. The general features of the model, however, probably will be retained. In conclusion, a theory on the basis of Debye and H¨ uckel ideas for electrolyte solutions with association of ions in dipoles and quadrupoles, hardcore effects and an approximated association constant, results in the location of the critical parameters consistent with numerical simulation estimatives. The only parameter in the theory is the distance between associated ions in the quadrupole complex, or equivalently the parameter s4 , which is used as a fitting parameter. It is found that 1.02 ≤ s4 ≤ 1.05 is a good choice.

References [1] Y. Levin, Rep. Prog. Phys. 65 (2002) 1577. [2] H. L¨owen, Phys. Rep. 237(5) (1994) 249. [3] P. W. Debye, E. H¨ uckel, Phys. Zeits. 24(3), 185 (1923). [4] T. L. Hill, An Introduction to Statistical Thermodynamics, Dover, New York, 1986. [5] T. L. Hill, Statistical Mechanics, Dover, New York, 1956. [6] D. A. McQuarrie, Statistical Mechanics, Harper and Row, New York, 1976.

17

[7] V. McGahay, M. Tomozawa, J. Non-Cryst. Solids 109 (1989) 27. [8] V. McGahay, M. Tomozawa, J. Chem. Phys. 97 (1992) 2609. [9] W. Ebeling, Z. Phys. Chem. 247 (1971) 340. [10] M. E. Fisher, J. Stat. Phys. 75 (1994) 1. [11] N. Bjerrum, K. Dan. Vidensk. Selsk. Mat. Fys. Medd. 7, (1926) 1. [12] W. Ebeling, Z. Phys. Chem. 238 (1968) 400. [13] H. L. Friedman, B. Larsen, J. Chem. Phys. 70 (1979) 92. [14] W. Ebeling, M. Grigo, Ann. Phys. 37 (1980) 21. [15] D. Kremp, W. Beskrownij, J. Chem. Phys. 104 (1996) 2010. [16] H. Yokoyama, H. Yamatera, Bull. Chem. Soc. Japan 48(6) (1975) 1770. [17] K. S. Pitzer, J. Phys. Chem. 88 (1984) 2689. [18] M. E. Fisher, Y. Levin, Phys. Rev. Lett. 71 (1993) 3826. [19] Y. Levin, M. E. Fisher, Physica A 225 (1996) 164. [20] J. S. Høye, K. Olaussen, Physica A 104A, 447 (1980). [21] E. Waisman, J. L. Lebowitz, J. Chem. Phys. 56 (1972) 3086. [22] E. Waisman, J. L. Lebowitz, J. Chem. Phys. 56 (1972) 3094. [23] L. Blum, J. S. Høye, J. Phys. Chem. 81, (1977) 1311. [24] J. S. Høye, J. L. Lebowitz, G. Stell, J. Chem. Phys. 61 (1974) 3253 [25] J. P. Valleau, J. Chem. Phys. 95 (1991) 584. [26] A. Z. Panagiotopoulos, Fluid Phase Equil. 76 (1992) 97. [27] A. Z. Panagiotopoulos, Fluid Phase Equil. 92 (1994) 313. [28] J. M. Caillol, J. Chem. Phys. 100 (1994) 2161. [29] G. Orkoulas, A. Z. Panagiotopoulos, J. Chem. Phys. 101 (1994) 1452. [30] G. R. Stell, K. C. Wu, B. Larsen, Phys. Rev. Lett. 37 (1976) 1369. [31] M. J. Gillan, Molec. Phys. 49 (1983) 421. 18

[32] V. C. Weiss, W. Schr¨oer, J. Chem. Phys. 108(18) (1998) 7747-. [33] H. Weing¨artner, V. C. Weiss, W. Schr¨oer, J. Chem. Phys. 113 (2000) 762. [34] B. Guillot, Y. Guissani, Mol. Phys. 87(1) (1996) 37. [35] G. Stell, J. Stat. Phys. 78(1/2) (1995) 197. [36] G.-W. Wu, M. Lee, K.-Y. Chan, Chem. Phys. Lett. 307 (1999) 419. [37] K. S. Pitzer, J. Chem. Phys. 104 (1996) 6724. [38] M. Valisk´o, D. Boda, J. Chem. Phys. 140 (2014) 234508. [39] J. A. Given, G. Stell, J. Chem. Phys. 106 (1997) 1195. [40] F. Moucka et al., J. Phys. Chem. B 116 (2012) 5468. [41] W. Schr¨oer, J. Mol. Liq. 92 (2001) 67. [42] J. M. Romero-Enrique, L. F. Rull, A. Z. Panagiotopoulos, Phys. Rev. E 66 (2002) 041204. [43] A. Z. Panagiotopoulos, M. E. Fisher, Phys. Rev. Lett. 88(4) (2002) 045701. [44] J.-N. Aqua, S. Banerjee, M. E. Fisher, Phys. Rev. E 72 (2005) 041501. [45] M. E. Fisher, J. Phys.: Cond. Matt. 8 (1996) 9103. [46] H. L¨owen, Phys. Rev. Lett. 74(6) (1995) 1028. [47] J. W. Zwanikken, P. K. Jha, M. O. de la Cruz, J. Chem. Phys. 135 (2011) 064106. [48] J. S. Høye, G. Stell, J. Chem. Phys. 67 (1977) 1776. [49] H. Nadolny, H. Weingartner, J. Chem. Phys. 114 (2001) 5273. [50] K. S. Pitzer, Chem. Phys. Lett. 105(5) (1984) 484. [51] P. J. Camp, G. N. Patey, Phys. Rev. E 60(1) (1999) 1063. [52] W. Schr¨oer, H. Weing¨artner, Pure Appl. Chem. 76(1) (2004) 19. [53] K. S. Pitzer, Acc. Chem. Res. 23 (1990) 333. 19

[54] Y. Zhou, G. Stell, J. Chem. Phys. 96 (1992) 1504. [55] Y. Zhou, G. Stell, J. Chem. Phys. 96 (1992) 1507. [56] Y. Zhou, S. Yeh, G. Stell, J. Chem. Phys. 102 (1995) 5785. [57] Y. Zhou, G. Stell, J. Chem. Phys. 102 (1995) 5796. [58] L. B. Bhuiyan et al., J. Chem. Phys. 100(11) (1994) 8301. [59] A. Z. Panagiotopoulos, J. Phys.: Condens. Matter 17 (2005) S3205. [60] R. Fantoni, G. Pastore, Phys. Rev. E 87 (2013) 052303. [61] P. S. Kuhn, M. C. Barbosa, Y. Levin, Macromolecules 31 (1998) 8347. [62] F. Bresme et al., Phys. Rev. E 51(1) (1995) 289. [63] J.-M. Caillol, J.-J. Weis, J. Chem. Phys. 102 (1995) 7610. [64] R. K. Wangsness, Electromagnetic Fields, John Wiley & Sons, New York, 1986. [65] M. Abramowitz, I. E. Stegun, Handbook of Mathematical Functions, National Bureau of Standards, Washington, 1964. [66] P. S. Kuhn, Chin. J. Phys. 49(2) (2011) 647. [67] N. F. Carnahan, K. E. Starling, J. Chem. Phys. 51(2) (1969) 635.

20