An analytical longitudinal dielectric function of primitive electrolyte solutions and its application in predicting thermodynamic properties

Accepted Manuscript Title: An analytical longitudinal dielectric function of primitiveelectrolyte solutions and its application in predicting thermodynamicproperties Author: Tiejun Xiao PII: DOI: Reference:

S0013-4686(15)30060-8 http://dx.doi.org/doi:10.1016/j.electacta.2015.06.145 EA 25270

To appear in:

Electrochimica Acta

Received date: Accepted date:

20-6-2015 30-6-2015

Please cite this article as: Tiejun Xiao, An analytical longitudinal dielectric function of primitiveelectrolyte solutions and its application in predicting thermodynamicproperties, Electrochimica Acta (2015), http://dx.doi.org/10.1016/j.electacta.2015.06.145 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.

Tiejun Xiao

a,b,∗

Guizhou Provincial Key Laboratory of Computational Nano-Material Science,Guizhou Normal College, Guiyang 550018, P. R. China b Guizhou Synergetic Innovation Center of Scientific Big Data for Advanced Manufacturing Technology, Guizhou Normal College, Guiyang, 550018, P. R. China

us

cr

a

ip t

An analytical longitudinal dielectric function of primitive electrolyte solutions and its application in predicting thermodynamic properties

an

Abstract

pt

ed

M

In this paper, the longitudinal dielectric function ϵl (k) of primitive electrolyte solutions is discussed. Starting from a modified mean spherical approximation, an analytical dielectric function in terms of two parameters is established. These two parameters can be related to the first two decay parameters k1,2 of the dielectric response modes of the bulk system, and can be determined using constraints of k1,2 from statistical theories. Furthermore, a combination of this dielectric function and the molecular Debye-H¨ uckel theory[J. Chem. Phys. 135(2011)104104] leads to a self-consistent mean filed description of electrolyte solutions. Our theory reveals a relationship between the microscopic structure parameters of electrolyte solutions and the macroscopic thermodynamic properties, which is applied to concentrated electrolyte solutions.

Ac ce

Keywords: electrostatics, dielectric function, electrolyte solution, mean spherical approximation, Debye-H¨ uckel theory 1. Introduction

The longitudinal dielectric function ϵl (k) of an isotropic solution is the wave-number-dependent longitudinal component of the wave-vector-dependent static dielectric tensor ϵαβ (k), which describes the charge response of a bulk ∗

[email protected]

Preprint submitted to Electrochimica Acta

June 19, 2015

Page 1 of 19

Ac ce

pt

ed

M

an

us

cr

ip t

system to an external electric potential [1, 2]. This dielectric function has important applications in solvation processes [3, 4, 5], in electron transfer processes [6, 7, 8, 9], in physics of electrostatic absorbtion [10, 11, 12, 13], and in ionic criticality [14, 15]. ϵl (k) can be evaluated experimentally or theoretically, however, and in most cases only numerical values of ϵl (k) are available. So a suitable model of the dielectric function would not only provide a reasonable fitting function of ϵl (k), but would also be a good starting point to study Coulomb-interaction-related phenomena. For polar fluids, the simplest model of ϵl (k) is to take ϵl (k) = ϵs with ϵs the bulk dielectric constant, where the k -dependence is neglected as in the Born model of solvation [16, 17, 18]. Later studies reveals the existence of dramatic k -dependence of ϵl (k), that is, the dielectric function becomes negative in certain region of k [19, 20], so more sophisticated models are proposed to capture such features. For example, a phenomenological Landau Hamiltonian which contains two independent polarization modes, is used to derive an analytical expression of ϵl (k) for polar fluids [21, 22]. This theory can well explain the main feature of dielectric function from experiments or from molecular dynamics simulations, and has been successfully applied to polar fluids such as water, acetonitrile and dimethyl [23]. It is also possible to build ϵl (k) with other number of polarization modes, and the interested readers can refer to Ref [24] for a comprehensive review. For electrolyte solutions, an expression of the dielectric function ϵs /ϵl (k) = 2 kD uckel(DH) theory [25], which describes 2 +k 2 can be found from the Debye-H¨ kD the dielectric response by a linearized Poisson-Boltzmann equation as ▽2 ϕ(r) = −k r 2 kD ϕ(r), with ϕ(r) = ϵqs e rD a Yukawa electric potential and kD the inverse Debye length that related to the ionic strength [26]. However, such a dielectric function is only valid for dilute solution, and is not suitable for concentrated electrolytes. Another systematic way to derive analytical expressions for ϵl (k) is to use analytical radial distribution function gij (r) from the integral equation theory, such as the mean spherical approximation(MSA) and its extensions [2, 27, 28]. In our previous studies on electron transfer processes in electrolyte solutions, a half empirical analytical functional has been used to fit the numerical ϵl (k) calculated from Hypernetted chain(HNC) theory or molecular dynamics simulations [8, 9]. In this study we demonstrate that the empirical function used in our previous studies can be derived from rigorous microscopic statistical theory of liquids. The correlation function hij (r) = gij (r)−1 and the direct correlation 2

Page 2 of 19

pt

ed

M

an

us

cr

ip t

function cij (r) are two important concepts in the theory of simple liquids, and are related to each other via the Ornstein-Zernike(OZ) equation [2]. Furthermore, the OZ equation combined with a suitable closure of cij (r) can be used to find the correlation function hij (r). Our strategy is to use a modified mean spherical approximation as the closure of cij (r), where the asymmetric part of cij (r) inside the core domain is a linear function of r, from which gij (r) and hence the longitudinal dielectric function ϵl (k) can be evaluated analytically. This ϵl (k) has two parameters which can be related to the first two decay parameters k1,2 of the dielectric response modes of the bulk solutions. Using analytical constraints on k1,2 from various statistical theories, actually we find a family of analytical dielectric function. Furthermore, this dielectric function is combined with our molecular Debye-H¨ uckel(MDH) theory [25] to give a mean field description of the electrolyte solutions, where the dielectric response is characterized by a linear combination of DH-like response modes, and another analytical expression of dielectric function is derived. Our study establishes a relationship between the microscopic structure parameters such as the first two decay parameters k1,2 and the thermodynamic properties of an electrolyte solution. Our strategy is applied to primitive model of electrolyte solutions, and the validity is tested for the thermodynamic properties. This article is organized as following : in section 2 an analytical expression for ϵl (k) of an electrolyte solution is built from a modified mean spherical approximation, in section 3 the MDH theory of electrolyte solutions on the basis of ϵl (k) is introduced. Application to concentrated electrolyte solutions is shown in section 4. Finally a brief conclusion is drawn in section 5.

Ac ce

2. The static longitudinal dielectric function ϵl (k) from a modified mean spherical approximation In this section, we will show how to derive an analytical ϵl (k) of electrolyte solutions based on a modified mean spherical approximation(MMSA). Consider a restricted primitive model of electrolyte solutions. Ions are immersed in a continuum with bulk dielectric constant ϵs , the reduced temperature is β = 1/kB T , e is the element charge and Zj = qj /e is the charge number of the j -th ion. The particle number density of the i -th species is ∑ ni , so n = i ni is the total number density, and xi = ni /n is the molar fraction of i -th type particle. The interaction potential between two ions is uij (r) = usij (r) + qϵisqrj , with usij (r) a hard sphere potential such that usij (r) = ∞, r ≤ σij and usij (r) = 0, r > σij . The cations and anions have the 3

Page 3 of 19

us

cr

ip t

same diameter and opposite charges such that σ++ = σ−− = σ+− = σ−+ = σ, q+ = −q− = q = Ze, where + and − denotes the cation and anion species, respectively. According to statistical theory of simple liquids, the microscopic structure of a liquid can be described by the correlation function hij (r) = gij (r) − 1 or the direct correlation function cij (r), and it is known that hij (r) and cij (r) are related to each other via the Ornstein-Zernike(OZ) equation [2]. We will use the MMSA closure combined with the OZ equation to find hij (r) of an electrolyte solution, and the details are as following. Due to the symmetry of the system, one can split the correlation function hij (r) and the direct correlation function cij (r) into a symmetric part and an asymmetric part [27], that is

an

h++ (r) = h−− (r) = hs (r) + ha (r), h+− (r) = h−+ (r) = hs (r) − ha (r), c++ (r) = c−− (r) = cs (r) + ca (r), c+− (r) = c−+ (r) = cs (r) − ca (r),

(1) (2) (3) (4)

ed

M

∫ ∫ Hereafter we denote fij (k) = dreik·r fij (r) = 0∞ fij (r)4π Sin(kr) r2 dr as the kr three dimensional Fourier transform of fij (r) with fij (r) = hij (r) and cij (r). The symmetric part cs (r) can be evaluated from the Percus-Yevick theory of hard sphere fluids as [2]

pt

r ηt1 r 3 cs (r) = −t1 − 6ηt2 ( ) − ( ) ,r ≤ σ σ 2 σ cs (r) = 0, r > σ

2

(5) (6)

2

Ac ce

(2+η) π 3 with t1 = (1−2η) , t2 = − 4(1−η) is the packing fraction of 4 , and η = 6 nσ (1−η)4 the electrolytes. This analytical cs (r) can be used to evaluate the symmetric part of the correlation function hs (r) based on the OZ equation [1 + nhs (k)][1 − ncs (k)] = 1. A modified mean spherical approximation is proposed to evaluate the asymmetric part ca (r) , such that

ca (r) = a + br, r ≤ σ βZ 2 e2 ca (r) = − , r > σ. ϵs r

(7) (8)

4

Page 4 of 19

The Fourier transformation of ca (r) reads ∫

Sin(kr) Cos(kσ)(a1 k 2 − a2 ) + kSin(kσ)a3 + a2 ca (r) = , kr k4n 0 (9) 2 2 with a1 = 4πaσn − 4πbnσ − kD , a2 = −8πbn and a3 = 8πbσn − 4πan. From now on we show how to evaluate ϵl (k) with ca (k). The static charge structure factor Szz (k) of the bulk system can be evaluated as [2] dr4πr2

Szz (k) =

∑

Zi2 xi + n

i

∑

ip t

∞

cr

ca (k) =

Zi Zj xi xj hij (k),

i,j

(10)

an

us

∫ with hij (k) = dreik·r hij (r) the Fourier transform of hij (r). One can evaluate ha (k) using the OZ equation [1 + nha (k)][1 − nca (k)] = 1, and then Szz (k) reads Z2 . (11) Szz (k) = Z 2 (1 + nha (k)) = 1 − nca (k)

M

According to the linear response theory in ionic fluids [2], the bulk dielectric function ϵl (k) is related to the static charge structure factor as

ed

ϵs 4πβne2 =1− Szz (k), ϵl (k) εs k 2 and then a response function defined as χ(k) ≡ 1 −

pt

χ(k) =

4πβne2 Szz (k). ϵ0 k 2

(12) ϵs ϵl (k)

can be evaluated as (13)

Ac ce

Combine Eq.(9), Eq.(11) and Eq.(13), the response function χ(k) is found to be 2 2 kD k χ(k) = 4 . (14) 2 k − (a1 k − a2 )Cos(kσ) − a3 kSin(kσ) − a2 Eq.(14) represents a family of dielectric function with various control parameters, and is our key result. As one can see, this equation can be taken as a justification for the empirical fitting function used in our previous studies [8, 9]. It is well known that the mean electrostatic potential ϕ(r) of an ion immersed in an electrolyte solution satisfies an asymptotic behavior as ϕ(r) ∼ ∑ e−kl r l αl r , with kl a decay parameter of the l -th Yukawa potential [29, 30, 5

Page 5 of 19

31, 32]. The decay parameters kl can be determined from the pole structure of the dielectric function or equivalently χ(k = ikl )−1 = 0,which now reads kl4 + (a1 kl2 + a2 )Cosh(kl σ) + a3 kl Sinh(kl σ) − a2 = 0.

(15)

us

cr

ip t

According to Eq.(15), k1,2 can be determined uniquely for a given combination of a and b or vice versa. So conceptually, any constraints for k1,2 can be used to determine a and b in the modified MSA theory, and leads to an analytical ϵl (k). It is noted that k1,2 could become complex numbers, which leads to oscillatory profiles of electric potentials [9]. To illustrate how the constraints on k1,2 from various statistical theories can be used to evaluate a and b in the modified MSA theory, we will list four statistical theories. The first constraints on kl is an equation derived by Kirkwood and Poirier(KP) [29], which reads

an

2 kl2 = kD Cosh(kl σ).

(16)

ed

M

2 This equation equals to set a = 0 and b = 0, so that a1 = −kD , a2 = 0 and a3 = 0. KP theory neglects the direct correlation function cD (r) inside the core domain, and one can expect it would be valid only for dilute solutions. The second one comes from a modified DH theory by Kjellander(KJ) [33], which reads 2 (17) ekl σ kD = (1 + kl σ)kl2 .

pt

In this case no analytical form of a and b is available and hence numerical calculations are required. The third one comes from the modified PoissonBoltmann theory by Outhwaite(OU) [31], which reads 2 kl Cosh(kl σ) + kD Sinh(kl σ) = kl3 (1 + kD σ)/kD .

(18) k2

k2

Ac ce

kD D and b = 0, so that a1 = − 1+kDD σ , This equation equals to set a = − 4πn 1+kD σ k3

a2 = 0 and a3 = − 1+kDD σ . As one can see, Eq.(16) , Eq.(17) or Eq.(18) combined with modified MSA leads to a and b which are not reported in the original theories. The forth constraints comes from the conventional MSA theory [27], which reads

with Γ =

MSA sets −8Γ4 and

kl (kl − 2Γ) = 2Γ2 (ekl σ − 1),

√ 1+2kD σ−1 an energy parameter of the solvent. 2σ 2 2 kD kD 2Γ Γ and b = 4πn ( 1+Γσ )2 , such that a = − 4πn 1+Γσ 3 a3 = −8Γ .

(19) It is known that a1 = −4Γ2 , a2 =

6

Page 6 of 19

In the next section, we will show how to use the analytical longitudinal dielectric function ϵl (k) to evaluate thermodynamic properties of the electrolyte solutions based on our MDH theory.

ip t

3. A Molecular Debye-H¨ uckel theory of electrolyte solutions based on ϵl (k)

l

an

us

cr

In this section, we briefly review the main results of our MDH theory. Our MDH theory is a reformulation of the dressed ion theory [10, 33], and can be build with the dielectric function ϵl (k) as input. We consider a restricted primitive model of electrolyte solution and denote ϕj (r) as the mean electric potential of a single ion j with charge qj and diameter σ. The electric potential ϕj (r) of the ion is expanded by a set of DH-like modes [25], such that ∑ ϕj (r) = Cl ϕjl (r), (20)

M

where each DH mode ϕjl (r) has its own decay parameter or effective Debye parameter kl and is defined as 4π qj δ(r), r ≤ σ ϵs ∇2 ϕjl (r) = kl2 ϕjl (r), r > σ.

(21) (22)

ed

∇2 ϕjl (r) = −

Ac ce

pt

Eq.(20), Eq.(21) and Eq.(22) are the central results of our MDH theory. Under this context, the solute is described by a hard body with certain molecular shape and charge distribution, while the solvent is characterized by a continuum with a set of DH-like response modes. {kl } and {Cl } are the key parameters in our MDH theory, and can be determined from the bulk dielectric function ϵl (k). kl (l = 1, 2, 3, ..., L) can be determined by the root structure of ϵl (k), such that ϵl (k = ikl ) = 0.

(23)

It is known that ϵl (k) undergoes a small k expansion as ∞ ∑ ϵs = Dn k 2n , ϵl (k) n=1

(24)

and then Cl satisfies a set of linear equations as [25] TC = Q,

(25)

7

Page 7 of 19

where T = {Tij } is a matrix with element Tij = √

4πβ

∑

2 kD

∑2j−1

ki2j−2

(ki σ)l /l! , (1+ki σ)

l=0

kD =

ni qi2

is the conventional inverse Debye length, and Q = {Qj } is a 2 vector with element Q1 = kD , Q2 = 1 and Qj = (−1)j−1 Dj−2 , (j ≥ 3). Eq.(25) can be solved as C = T−1 Q. In our previous study, it is found that the inclusion of several DH modes (typically L = 2 or 4) are capable to predict the thermodynamics of electrolyte solutions or molten salts [25, 8, 9], and so one can always safely truncate the elements in Eq.(25) to a finite number. With the sets of {kl } and {Cl }, the mean electric potential ϕj (r) of the j-th ion can be evaluated explicitly as [25]

us

cr

ip t

i ϵs

qj + ψj , r≤σ ϵs r ∑ qj Cl e−kl (r−σs ) ϕj (r) = ,r > σ r l ϵs (1 + kl σs )

an

ϕj (r) =

(26) (27)

where ψj is an induced potential in the ion center as qj ∑ Cl kl . ϵ s l 1 + kl σ

M

ψj = −

(28) ∑

pt

ed

Note that our MDH theory leads to analytical form for i qi ni hij (r) [25], from which the static structure factor Szz (k) and hence the response function χ(k) can be rewritten as function of in terms of Cl and kl . After a straightforward calculation, it is found that 2 ∑ Cl kl2 Cos(kσ) + kl Sin(kσ)/k kD [1 − ]. k2 k 2 + kl2 l 1 + kl σ

(29)

Ac ce

χ(k) =

This analytical expression for dielectric function is another key result of this study, where the contribution from each DH-like mode is singled out. Similar physical picture can be found in the dielectric function of polar fluids as noted by Medvedev, where the dielectric response cames from a linear combination of polarization modes [23]. ex The electrostatic energy per particle βUN can be evaluated with the induced potential, which reads [25] βU ex ∑ xj qj ψj k 2 ∑ Cl kl = =− D . N 2 8πn l 1 + kl σ j

(30)

8

Page 8 of 19

In the weak electrostatic coupling limit where the first Debye mode is dominant, the above equation reduces to the familiar DH result as [25] 3 kD βU ex =− . N 8πn(1 + kD σ) βU ex N

ip t

The electrostatic energy per particle from the dielectric function as [25]

(31)

can also be evaluated directly

cr

βU ex 1 ∫∞ 2 = 2 [χ(k)k 2 − kD ]dk. N 4π n 0

(32)

an

us

One can verify that Eq.(29) together with Eq.(32) lead to the same electrostatic energy as in Eq.(30). The hard sphere interaction also has a contribution to the thermodynamic properties such as the free energy and the pressure. The excess Gibbs free energy per particle can be evaluated as [25] (33)

M

β∆Gs βU ex β∆Gex = + , N N N s

ed

with β∆G a contribution from the hard sphere interactions, which can be N evaluated using the Canahan-Starling formulas as long as the electrostatic coupling is not too strong [34]. The reduced pressure can be evaluated from the virial equation as (34)

pt

1 βU ex βP 2πnσ 3 ∑ xi xj gij (σ) + =1+ , n 3 3 N ij

Ac ce

where gij (σ) is the contact value of the radial distribution function, and could be determined from statistical theories such as exponential approximation [35]. The excess Helmholtz free energy can be evaluated from the thermodynamic relations as β∆Gex βP β∆Gs 2πnσ 3 ∑ β∆Aex 2 βU ex = −( −1) = − . (35) xi xj gij (σ)+ N N n N 3 3 N ij

Now all excess thermodynamic properties contributed from Coulomb interaction are evaluated explicitly. One can see that a family of dielectric function under MMSA is established, and hence a self-consistent theory of electrolyte solutions can be built. 9

Page 9 of 19

cr

ip t

Firstly, k1,2 is found from an analytical constrains such as Eq.(16), Eq.(17), Eq.(18) or Eq.(19), or any numerical or experimental method. Secondly, k1,2 is used in Eq.(15), from which a and b can be determined analytically or numerically, and hence ϵl (k) is available according to Eq.(14). Thirdly, {kl } and {Cl } can be solved according to Eq.(25) with ϵl (k) as input. Those parameters combined with Eq.(20) and Eq.(21) lead to a mean field description of the electrolytes, and can be used to predict the thermodynamic properties of electrolyte solutions. 4. Application to electrolyte solutions

Ac ce

pt

ed

M

an

us

In this section, our theory is applied to electrolyte solutions. Specifically we compare several statistical theories that leads to k1,2 , and test their ability to predict thermodynamic properties of electrolyte solutions. We consider a restricted primitive model of electrolyte√solution with Z = 2 1, σ = 1, ϵs = 1, β = 1, where the Debye parameter kD = 4πβnZ is used as ϵs the control parameter. When particle number density n varies from 0.01 to 0.8, kD varies from 0.35 to 3.17. As a comparison, for a model NaCl aqueous solution with a relative dielectric constant ϵr = 78.5, temperature T = 300K and ion diameter σ = 4.2˚ A at the concentration of 1 mol/l and 5 mol/l, the reduced Debye parameter kD σ is about 1.4 and 3.1, respectively. As the first test of the validity of our modified MSA theory, we consider the bulk dielectric function. The HNC theory is known to be very accurate in a wide range of the parameter space [36], and is used as a benchmark in this study. The correlation function hij (r) of the solvent species from HNC calculation is used to evaluate Szz (k) and the response function χ(k). The s numerical value of χ(k) = 1 − ϵlϵ(k) from HNC theory is fitted to Eq.(14) with a and b as adjustable parameters. For particle number density n = 0.1, 0.4, 0.7, the Debye parameters are kD = 1.12, 2.24, 2.97, respectively. As shown in Fig.1(a), Fig.1(b) and Fig.1(c), our proposed dielectric function model works well compared with the HNC theory, as long as suitable a and b or equivalently the first two decay parameters k1,2 are used. We also note that χ(k) from MMSA theory becomes less accurate for stronger electrostatic coupling, e.g., noticeable difference between MMSA and HNC theory appears for a large number density n = 0.7. This may be related to the fact cij (r) from MMSA is a linear function of r for r ≤ σ, while cij (r) from HNC is a nonlinear function of r for r ≤ σ, and the nonlinearity becomes stronger for larger electrostatic coupling. 10

Page 10 of 19

1.0

HNC MMSA

0.8

(k )

n = 0.1

0.4

0.2

0.0 2

4

k

6

8

cr

0

ip t

0.6

us

(a) n = 0.1

1.0

HNC

MMSA

an

0.8

(k )

n = 0.4

0.6

0.2

0.0 2

4

k

6

8

ed

0

M

0.4

(b) n = 0.4

1.2

pt

HNC MMSA

Ac ce

(k )

1.0

n = 0.7

0.8

0.6

0.4

0.2

0.0 0

2

4

k

6

8

(c) n = 0.7

s Figure 1: The response function χ(k) = 1 − ϵlϵ(k) of a restricted primitive model of electrolyte solution under n = 0.1, 0.4, 0.7 from HNC theory(filled square) and our MMSA theory(hollow star). The lines are guides to the eye.

11

Page 11 of 19

Ac ce

pt

ed

M

an

us

cr

ip t

As the second test, we focus on the dependence of k1,2 on kD . Theoretical prediction for k1,2 from Eq.(16) of KP theory, Eq.(17) of KJ theory and Eq.(18) of OU theory are plotted in Fig.2(a) and Fig.2(b). Eq.(19) of MSA theory has been tested in our previous study [25], which is found to be have a very similar k1,2 structure compared with Eq.(18) and will not show here. k1,2 from HNC calculation is also shown as a reference, where the response function is fitted to the function in Eq.(14), and then k1,2 is determined from Eq.(15). As one can see, these theories has similar decay parameter structures, that is, k1,2 are real numbers for small kD and become complex conjugates when kD is beyond some critical value kc . kc from KP, KJ, OU and MSA theory are about 1.03, 1.35, 1.24 and 1.23 respectively, which could be compared with a value kc = 1.29 from HNC. KP theory has neglect the ionic correlation, and leads to k1,2 which is unsatisfactory for large kD . ex As the third test, we consider the electrostatic energy βUN of the electrolyte solution. Using k1,2 from various theories, a and b in the modified MSA can be determined using Eq.(15), and then the dielectric function ϵl (k) can be evaluated from Eq.(14). With ϵl (k), one can solve Eq.(23) and Eq.(25) to find the set of {kl } and {Cl } with l = 1, 2, ..., L. Hereafter we test three cases with L = 1, 2, 4, where singe mode, two modes and four modes are considered in the MDH theory. In general, one can expect that a theory with kl structure more close to the HNC theory would be more satisfactory. ex For the single DH mode case, the electrostatic energy βUN for 0.1 < kD < 1.4 is shown in Fig.3(a). Note that k1 undergoes a transition from a real number to a complex value as kD increases, while the single mode theory requires k1 > 0, one can expect that the single mode theory is valid only in a limited region of kD < kc . When kD is small, constraints mentioned in this study all lead to k1 ≃ kD , and hence predict electrostatic energies which are almost indistinguishable from HNC and DH theory. ex For two DH modes case, βUN for 0.3 < kD < 3.2 is shown in Fig.3(b). The two-modes theory does not require k1,2 to be real value, and hence is applicable to a much wider range of kD than for the single mode case. Due to the fact that k1,2 from KP is unsatisfactory for large kD , KP leads to a result which is only comparable with DH theory, and hence is not useful for practical applications. KJ and OU all lead to electrostatic energies which are much more accurate than the conventional DH theory, where the maximum difference compared to HNC theory is 7 and 10 percent, respectively, while DH theory has an error up to 25 percent. ex For four DH mode case, βUN for 0.5 < kD < 3.2 is shown in Fig.3(c). Note 12

Page 12 of 19

Ac ce

pt

ed

M

an

us

cr

ip t

that the evaluation of thermodynamic properties of the solvent relies on the accuracy of four kl , it is expected that KP theory has the worst performance among the four theories mentioned in this study. As one can see, KP theory even leads to electrostatic energies which are much worse than DH theory. The electrostatic energies predicted from KJ and OU theory have a maximum error of 16 percent and 10 percent when compared with HNC, and is still much better than DH theory. MSA leads to electrostatic energies which are very close to that from OU theory (where the difference is smaller than 2 percent for 0.5 < kD < 3.2), and is not shown. In our previous study on electrolyte solutions, χ(k) from MSA or HNC is used as input of MDH theory, and it is found that MDH theory gets more accurate as the number L of DH modes increases [25]. Herein we find something different when χ(k) comes from other theories. When χ(k) from KP theory is used, MDH theory with smaller L is found to be better, i.e., MDH theory with L = 1 works best among L = 1, 2, 4. When χ(k) from KJ theory is used, we find that MDH theory with L = 2 has be best performance. When χ(k) from OU theory is used, it is noted that MDH with larger L is more accurate, and MDH with L = 4 leads to the most accurate results. Note that the thermodynamic properties are explicit functions of decay parameters according to MDH theory, one may understand the phenomenon observed in this study, that is, when {kl }(l = 1, 2, ..., L) is accurate, one expect that including more DH modes will improve the results, however, when {kl } differs quite bit from that of HNC theory, including more DH modes in the MDH theory may even leads to less accurate thermodynamic properties. According to the above discussions, our MMSA together with suitable k1,2 does lead to a realistic model of the dielectric function ϵl (k), which is further used to establish a self-consistent mean field theory of electrolyte solutions. Note that the analytical constraints of kl mentioned in this study depends only on kD , while HNC theory shows that kl depends on both kD and the reduced inverse temperature β [36], one can expect that a more sophisticated theoretical constraint of k1,2 will lead to a better dielectric description of electrolyte solutions. 5. Concluding Remarks In conclusion, an analytical expression for the dielectric function ϵl (k) of an electrolyte solution is derived using the MMSA theory, where the two parameters a and b in the dielectric function can be related to the first two 13

Page 13 of 19

ip t

9

KJ

7

Re(

1

Re(

2

)

cr

k k

8

)

6

HNC 4

OU

3

KP

2 1 0 0.5

1.0

1.5

kD

2.0

2.5

3.0

3.5

an

0.0

us

5

3

M

(a) Re(k1,2 ) from various theories

Im(k ) 1

KJ

OU

2

ed

HNC

KP

1

0

0.5

pt

0.0

1.0

1.5

kD

2.0

2.5

3.0

3.5

Ac ce

(b) Im(k1 ) from various theories

Figure 2: The first two decay parameters k1,2 from various theories. Re(kl ) and Im(kl ) denote the real part and the imaginary part of kl , respectively. The lines are guides to the eye.

14

Page 14 of 19

HNC DH

/N

-0.1

ex

KP

U

KJ OU

-0.3

-0.4 0.3

0.6

k

0.9

1.2

D

us

(a) MDH theory with one DH mode

cr

ip t

-0.2

-0.1

HNC

an KP KJ

U

ex

/N

DH

-0.2

OU

-0.3

M

-0.4

-0.5

1.0

1.5

k

2.0

2.5

3.0

D

ed

0.5

pt

(b) MDH theory with two DH modes

HNC DH KP

-0.2

KJ

ex

/N

-0.1

Ac ce

U

OU

-0.3

-0.4

-0.5 0.5

1.0

1.5

k

2.0

2.5

3.0

D

(c) MDH theory with four DH modes

Figure 3: Electrostatic energies from MDH theory with various DH modes. The lines are guides to the eye. 15

Page 15 of 19

cr

ip t

decay parameters k1,2 of the dielectric response modes of the bulk solution, and can be determined according to constraints of k1,2 from various statistical theories. ϵl (k) is further used as input of our MDH theory, which leads to a mean field theory of electrolytes, where the thermodynamic properties of the bulk solution can also be evaluated analytically using parameters of dielectric response modes. Our theory establishes a relationship between the parameters k1,2 and the bulk thermodynamic properties, and is applied to concentrated electrolyte solutions. ACKNOWLEDGMENTS

M

an

us

We acknowledge the financial support from the National Natural Science Foundation of China (Grant No. 21403041), the NSF of department of education of Gizhou province(No. QJTD[2013]16), the NSF of Guizhou province (No. QKJ[2014]2139), and Construction Project for Guizhou Provincial Key Laboratories (Z[2013]4009). This work is also supported by a startup package from Guizhou Normal College. References

[1] J. D. Jackson, Classical Electrodynamics, John Wiley & Sons, Inc, 1998.

ed

[2] J. Hansen, I. McDonald, Theory of Simple Liquids, 2nd, 1986.

pt

[3] B. Mennucci, R. Cammi, W. Interscience, Continuum solvation models in chemical physics: from theory to applications, Wiley Online Library, 2007.

Ac ce

[4] S. Bose, R. Adhikary, C. A. Barnes, D. B. Fulton, M. S. Hargrove, X. Song, J. W. Petrich, Comparison of the dielectric response obtained from fluorescence upconversion measurements and molecular dynamics simulations for coumarin 153- apomyoglobin complexes and structural analysis of the complexes by nmr and fluorescence methods, The Journal of Physical Chemistry A 115 (16) (2010) 3630–3641. [5] M. Maroncelli, X.-X. Zhang, M. Liang, D. Roy, N. P. Ernsting, Measurements of the complete solvation response of coumarin 153 in ionic liquids and the accuracy of simple dielectric continuum predictions, Faraday discussions 154 (2012) 409–424. 16

Page 16 of 19

[6] R. Marcus, Chemical and electrochemical electron-transfer theory, Annual Review of Physical Chemistry 15 (1) (1964) 155–196.

ip t

[7] A. Kornyshev, M. Tosi, J. Ulstrup, Electron and ion transfer in condensed media, World Sci.

cr

[8] T. Xiao, X. Song, Reorganization energy of electron transfer processes in ionic fluids: A molecular debye-h¨ uckel approach, The Journal of chemical physics 138 (11) (2013) 114105.

us

[9] T. Xiao, X. Song, A molecular debye-h¨ uckel approach to the reorganization energy of electron transfer reactions in an electric cell, The Journal of chemical physics 141 (13) (2014) 134104.

an

[10] R. Kjellander, D. J. Mitchell, Dressed-ion theory for electrolyte solutions: A debye–h¨ uckel-like reformulation of the exact theory for the primitive model, The Journal of chemical physics 101 (1) (1994) 603– 626.

M

[11] R. Kjellander, D. J. Mitchell, Dressed ion theory for electric double layer structure and interactions; an exact analysis, Molecular Physics 91 (2) (1997) 173–188.

ed

[12] J. Ulander, H. Greberg, R. Kjellander, Primary and secondary effective charges for electrical double layer systems with asymmetric electrolytes, The Journal of Chemical Physics 115 (15) (2001) 7144–7160.

Ac ce

pt

[13] B. Forsberg, J. Ulander, R. Kjellander, Dressed ion theory of sizeasymmetric electrolytes: Effective ionic charges and the decay length of screened coulomb potential and pair correlations, The Journal of chemical physics 122 (6) (2005) 064502. [14] M. E. Fisher, B. P. Lee, Ginzburg criterion for coulombic criticality, Physical review letters 77 (17) (1996) 3561. [15] J.-N. Aqua, M. E. Fisher, Ionic criticality: an exactly soluble model, Physical review letters 92 (13) (2004) 135702. [16] D. Bashford, D. A. Case, Generalized born models of macromolecular solvation effects, Annual Review of Physical Chemistry 51 (1) (2000) 129–152. 17

Page 17 of 19

[17] A. Onufriev, D. Bashford, D. A. Case, Exploring protein native states and large-scale conformational changes with a modified generalized born model, Proteins: Structure, Function, and Bioinformatics 55 (2) (2004) 383–394.

cr

ip t

[18] T. T. Duignan, D. F. Parsons, B. W. Ninham, A continuum model of solvation energies including electrostatic, dispersion, and cavity contributions, The Journal of Physical Chemistry B 117 (32) (2013) 9421– 9429.

us

[19] F. O. Raineri, H. Resat, H. L. Friedman., Static longitudinal dielectric function of model molecular fluids., The Journal of chemical physics 96 (4) (1992) 3068–3084.

an

[20] P. A. Bopp, A. A. Kornyshev, G. Sutmann, Static nonlocal dielectric function of liquid water, Physical review letters 76 (8) (1996) 1280.

M

[21] A. A. Kornyshev, G. Sutmann, The shape of the nonlocal dielectric function of polar liquids and the implications for thermodynamic properties of electrolytes: A comparative study, The Journal of chemical physics 104 (4) (1996) 1524–1544.

ed

[22] A. Kornyshev, S. Leikin, G. Sutmann, overscreening in a polar liquid as a result of coupling between polarization and density fluctuations, Electrochimica acta 42 (5) (1997) 849–865.

Ac ce

pt

[23] I. G. Medvedev, The analytical expression for the static nonlocal dielectric function of a polar liquid with due account of the overscreening effect, Electrochimica acta 49 (2) (2004) 207–217. [24] J. P. Bardhan, Gradient models in molecular biophysics: progress, challenges, opportunities, Journal of the mechanical behavior of materials 22 (5-6) (2013) 169–184. [25] T. Xiao, X. Song, A molecular debye-h¨ uckel theory and its applications to electrolyte solutions, The Journal of chemical physics 135 (10) (2011) 104104. [26] D. A. McQuarrie, Statistical Mechanics, Harper and Row, New York, Evanston, San Francisco, London, 1976.

18

Page 18 of 19

[27] L. Blum, Primitive electrolytes in the mean spherical approximation, in: Theoretical Chemistry: Advances and Perspectives, Vol. 5, Academic Press: New York, 1980, pp. 1–66.

ip t

[28] L. M. Varela, M. Garcıa, V. Mosquera, Exact mean-field theory of ionic solutions: non-debye screening, Physics reports 382 (1) (2003) 1–111.

cr

[29] J. G. Kirkwood, J. C. Poirier, The statistical mechanical basis of the debye–h¨ uekel theory of strong electrolytes, The Journal of Physical Chemistry 58 (8) (1954) 591–596.

us

[30] C. W. Outhwaite, Extension of the debye–h¨ uckel theory of electrolyte solutions, The Journal of Chemical Physics 50 (6) (1969) 2277–2288.

an

[31] C. Outhwaite, The linear extension of the debye-h¨ uckel theory of electrolyte solutions, Chemical Physics Letters 5 (2) (1970) 77–79.

M

[32] C. W. Outhwaite, M. Molero, L. B. Bhuiyan, Primitive model electrolytes in the modified poisson–boltzmann theory, Journal of the Chemical Society, Faraday Transactions 89 (9) (1993) 1315–1320.

ed

[33] R. Kjellander, Modified debye-h¨ uckel approximation with effective charges: an application of dressed ion theory for electrolyte solutions, The Journal of Physical Chemistry 99 (25) (1995) 10392–10407.

pt

[34] N. F. Carnahan, K. E. Starling, Equation of state for nonattracting rigid spheres, The Journal of Chemical Physics 51 (2) (1969) 635–636.

Ac ce

[35] H. C. Andersen, D. Chandler, Optimized cluster expansions for classical fluids. i. general theory and variational formulation of the mean spherical model and hard sphere percus-yevick equations, The Journal of Chemical Physics 57 (5) (1972) 1918–1929. [36] E. Guti´errez-Valladares, M. Luksic, B. Mill´an-Malo, B. Hribar-Lee, V. Vlachy, Primitive model electrolytes. a comparison of the hnc approximation for the activity coefficient with monte carlo data, Condensed Matter Physics 14 (3) (2011) 1–15.

19

Page 19 of 19