Nonlocal nonlinear static dielectric response of polar liquids

Journal of Electroanalytical Chemistry 450 (1998) 143 – 156

Nonlocal nonlinear static dielectric response of polar liquids Alexei A. Kornyshev 1, Godehard Sutmann * Institut fu¨r Energie6erfahrenstechnik, Forschungszentrum Ju¨lich, D-52425 Ju¨lich, Germany Received 5 August 1997; received in revised form 6 October 1997

Abstract A linear, phenomenological theory of the nonlocal dielectric response is extended to the non-linear case. For small external electric fields, we give solutions within a perturbation approach. Results, obtained by molecular dynamics computer simulations for the non-linear response of water are in very good agreement with the theory predictions. The results are applied to the calculation of the electrostatic potential of ions in solution. For the electric field strengths, for which the perturbation theory is still applicable, there are only minor differences with respect to the linear case. A fitting procedure to the simulation data, which covers a wider range of electric fields shows, however, that at larger fields the non-linear effects may change the results considerably, relative to those obtained on the basis of the linear response approximation. We discuss a generalization of the nonlocal electrostatic approach to the calculation of solvation energies, which takes into account non-linear effects. © 1998 Elsevier Science S.A. All rights reserved. Keywords: Nonlocal dielectric response; Polar liquids; Solvation energies

1. Introduction The wave vector dependent dielectric response tensor, xab (k), determines the response of a polar liquid to an imbedded microscopic charge distribution at distances of the characteristic liquid structure. It appears in the thermodynamics of electrolytes, charge transfer kinetics, electrochemistry, adsorption and interaction of biological membranes and macromolecules. When the response is linear, only the longitudinal component, x(k) =kakbxab /k 2, appears in the description of electrostatic phenomena. Within the linear response approximation a variety of solvation determined phenomena were described, but its validity was never verified for realistic polar liquids. For the case of associated liquids, e.g. water or methanol, one will expect non-linear be* Corresponding author. Present address: Dipartimento di Fisica Universita’ di Trento, Povo, Trento, I-38050, Italy. E-mail: [email protected] 1 Also a member of The A.N. Frumkin Institute of Electrochemistry of the Academy of Sciences, 117071 Moscow, Russia. 0022-0728/98/$19.00 © 1998 Elsevier Science S.A. All rights reserved. PII S0022-0728(97)00622-0

havior at least, when the interaction energy between the molecular dipole and external field is so large that the hydrogen bonds are broken (mE0 \ 0.25 eV). When the response is linear it is possible to apply the linear response theory formalism to the calculation of x(k). In this case the dielectric function is related to the correlation function of equilibrium bound charge density fluctuations via the fluctuation dissipation theorem (FDT). In the last decade, a number of papers were devoted to the calculation of e(k) (which is related to the response function through x(k)=1− 1/e(k)), by using computer simulation methods or different statistical mechanical schemes2. In a recent work we calculated the linear response function for water by means of molecular dynamics computer simulation [11,12], which gave a more complete picture of x(k), in close proxim2

Molecular dynamics computer simulations to the calculation of e(k) for different types of polar liquids were used in Refs. [1– 3,11,12]. Statistical mechanical schemes were used in Refs. [4–9]. A combination of computer simulation and optimized cluster theory was developed in ref. [10].

144

A.A. Kornyshe6, G. Sutmann / Journal of Electroanalytical Chemistry 450 (1998) 143–156

ity to experimental data for this function which were extracted from neutron scattering measurements. In Ref. [13], a phenomenological theory was presented which was able to describe the main features of the linear response function, x(k), found from the simulation. In order to include non-linear effects in the theory of solvation, different field dependent dielectric functions were considered [14 – 24], which increased monotonically from the high frequency dielectric constant, e to the bulk dielectric constant e, with the variation of the ionic electric field with the distance from the ion. This semi-empirical approach does not take into account spatial correlations of polarization fluctuations and neglects the tensor character of xab. The description of the non-linear dielectric response with account for the spatial dispersion becomes in the general case a difficult task. There is no closed form fluctuation dissipation theorem, no superposition principle for different k modes (i.e. no decoupling of the modes) and there is no splitting of the response function into a longitudinal and a transversal components. In a homogeneous medium, the constitutive relation between polarization P(k) and external field E(0)(k) (which coincides with the induction in this case [25]) reads [26] (0) Pa (k) = x (0) ab (k)E b (k)

&

(0) (0) + dk%dk¦d(k −k%− k¦)x (1) abg (k%,k¦)E b (k%)E g (k¦) + ···

(1) where (···) denotes integrals of higher powers in E(0), which contain dielectric response tensors of increasing rank. In the linear response limit the rhs of Eq. (1) reduces to the longitudinal component of the first item. There is no general algorithm to evaluate the higher order rank tensors. However, to reveal the limits of the linear response assumption one may focus only on the diagonal components of the dielectric response tensor. Recently we performed a series of computer simulations [27], where a water system was exposed to an external, inhomogeneous electric field of the form E(0)(r) = k. E0 cos kr (k. =k/ k ). By varying the external field strength, we were able to study the crossover from the linear to the non-linear dielectric response as a function of the wavenumber k. It was found that the response at ˚ − 1 for the resonance-like maximum of x(k) (at k : 3 A water [11,12]) showed the first occurrence of non-linear ˚. behavior at field strengths of E0 :0.05 V/A In the present paper, we extend the phenomenological theory, proposed in Ref. [13] to the non-linear case. We use a perturbation approach in order to solve the resulting non-linear equation. For the case of water we compare our findings with computer simulation results. In order to explore the possible consequences of the non-linear response, we calculate the potential of hy-

drated ions. We also discuss results which are obtained with an approximation for the response function which reproduces the simulation results in the whole range of ˚ . Within applied electric fields, i.e. E0 [0.001,2.5] V/A this approximation, we calculate hydration energies of singly charged ions for which we extend the nonlocal electrostatic approach into the non-linear regime. When an ion is imbedded into the solvent, it also distorts the structure of the liquid through its finite size. In Appendix A, we present a formalism for a description of these kinds of effects and discuss their importance.

2. Mean field theory As well as in Ref. [13], we assume that the local polarization of a polar liquid may be described by three contributions P= P1 + P2 + P

(2)

Here, P stands for the high frequency, electronic polarizability. The correlations in space of the fluctuations of P are neglected, i.e. x = const. P1 and P2 describe two subsystems in the liquid, a slow one and a fast one, respectively. This description allows us to distinguish between different contributions, such as rotational diffusion and molecular reorientations (slow motions) and librations and intramolecular vibrations (fast motions). The fast polarization component P2 is assumed to be coupled to a scalar order parameter h, which represents local density fluctuations. Here, we assume a quadratic Hamiltonian for the slow component P1 and a fourth order expansion for the fast component, which describes therefore an anharmonic fluctuation pattern. The Landau free energy may then be written as

&

H= (h1 + h2 + h )dr

(3)

where a1 (4) h1 = [P21 + l 21(9P1)2]− P1D 2 a b g h2 = 2[P22 + l 22(9P2)2]+ [h 2 + l 2(9h)2]+ P29h −P2D 2 2 2 + sP42 h = − P D

(5) (6)

The parameters ai, b, g determine the weight of the contributions and li, l are characteristic correlation lengths. A nonzero s accounts for an anharmonicity of the polarization fluctuations. The minimum of the Landau free energy is found from the Euler-Lagrange equations, which read dH =b(h− l 292h)−g9P2 = 0 dh

(7)

A.A. Kornyshe6, G. Sutmann / Journal of Electroanalytical Chemistry 450 (1998) 143–156

dH =a1(P1 −l 2192P1) −D = 0 dP1

(8)

dH =a2(P2 −l 2292P2) +4sP32 +g9h − D =0 dP2

(9)

Via a Fourier transformation of Eq. (8), we may directly calculate the response function x1(k) = 4p P1(k) / D(k) which gives a Lorentzian with the correlation length l1 x1(k) =





1 1 1 − e  e 1+ l 21k 2

(10)

The prefactor, containing the macroscopic dielectric constant e and an intermediate dielectric constant, e , which separates the response of the slow and the fast polarization, results from a calculation of the factor a1 [13] In order to get a solution for P2 as a function of the induction field D, we differentiate Eq. (7) with respect to r and substitute 9h from Eq. (9), which leads to P2 − A 292P2 + B 492(92P2) + 1 = (D − l 292D) a2

4s 3 [P2 −l 292(P32)] a2

the polarization P a series solution in powers of the external field: (2) 3 (4) 5 P2(r) =P(0) 2 (rE0)+ P2 (rE 0)+ P2 (rE 0)+ ···

where we have introduced A = l 2 +l 2 −L 2, B = ll and L 2 =g 2/a2b. Eq. (11) has several important features: (i) the anharmonic term appears not as a single cubic term, but also as a second derivative; (ii) the induction field also appears with a second derivative; (iii) the anharmonic term, as well as the induction field term contain l, the characteristic length scale of the local density fluctuations. These effects are due to the coupling between h and P2, i.e. they vanish for a coupling parameter g=0. It is only due to the cubic term that Eq. (11) may not be Fourier transformed to get a closed form solution for x2(k). However, one may solve Eq. (11) within the framework of the perturbation theory. If we assume the external field, E(0)(r), to have the form E0(r) =k. E0 cos kr =D(r), where E0 is small, one may consider the term sP42 as a perturbation of the quadratic Hamiltonian and one may apply a perturbation theory to obtain a series solution in powers of E0. We present such a solution in the next section and evaluate the validity criterion for the perturbation approach.

3. Results

3.1. Perturbation solution For very small external fields one assumes the dielectric response to be linear, i.e. the forth order term in Eq. (5) will be negligible. This leads us to assume for

(12)

In an isotropic homogeneous medium, the polarization is an asymmetric function of the external field, i.e. P[E(0)]= − P[ −E(0)] and all even powers of E0 vanish in Eq. (12). Since the relation between P and D has, generally, nonlocal character, the dependence of the Pi s on E0 has a functional character. The expansion, Eq. (12), allows us to solve Eq. (11) iteratively. Here we give the results for the first three terms. The evaluation of higher terms is made in Section 3.2 in order to prove the convergence of the series. First order term ( E0): The polarization responds only to the ground wave with wave vector k, i.e. P(0) 2 = k. P01 cos kr P01 = x2(k)E0

(13)

Third order term ( E 30): A higher harmonic appears in the solution, i.e. P(2) 2 = k. P21 cos kr+k. P23 cos 3kr, where P21 = − 3x2(k)4E 30

(11)

145

P23 = − x2(k) x2(3k)E 3

(14) 3 0

(15)

Fifth order term (  E 50) contains three harmonics, P(4) 2 = k. P41 cos kr+k. PP43 cos 3kr+ k. P45 cos 5kr, where P41 = 3x2(k)6[9x2(k)+ x2(3k)]E 50

(16)

P43 = 3x2(k)5x2(3k)[3x2(k)+ 2x2(3k)]E 50

(17)

P45 = 3x2(k)5x2(3k)x2(5k)E 50

(18)

The solutions are expressed in terms of the linear response function, x2(k), which reads [13] x2(k)=





1 1 1+ l 2k 2 − e e  (1+l 22k 2)(1+ l 2k 2)− L 2k 2

(19)

From Eqs. (13)–(16), it is obvious that each higher, odd power of E0 adds higher harmonics to the ground wave of the polarization. The fact, that only odd partial waves appear in the solution is a direct consequence of the fourth order polarization term in the Hamiltonian. Only if one introduced an odd power of P in the Hamiltonian, would the solution contain both even and odd harmonics. However, the appearance of a cubic power of P would imply a non-vanishing equilibrium polarization, which does not exist in non-ferroelectric materials. If the linear response function has a large maximum at a certain wavenumber k, it is obvious from Eqs. (13)–(16), that the main contribution will come from those terms which contain the highest power of the ground wave, i.e. the Pi1s. Indeed, we find that for a parametrization which fits the response function of

A.A. Kornyshe6, G. Sutmann / Journal of Electroanalytical Chemistry 450 (1998) 143–156

146

water [11,12] (which has a large maximum at k: 3 ˚ − 1) the higher harmonic terms, Pij s ( j\ 1) are negliA gible. This tendency is weaker for a Lorentzian type function (L=0), but the higher harmonic terms still lead only to small corrections.

approximations with the calculated ci in Fig. 1, from where we find that the best fit is obtained with the convergent c (3) with parameters a˜ =14.3, b0 = 0.47. This i leads us to the conclusion that our perturbation approach to Eq. (11) is valid and the ci are convergent with cm = limi “ ci : 17.5.

3.2. Con6ergence of the perturbation series The validity criterion of the perturbation approach to Eq. (11) must define the maximum value of the external field E0 for which the series expansion of the polarization, Eq. (12), still converges, i.e. for which (i + 1) P(i) , Öi. After each iteration step, the structure 2 BP2 of the expressions for the Pij s become more complex: they consist of terms of the fundamental wave with wave number k and higher harmonic terms with the maximum wave number kmax =jk ( j refers to the index of the Pij s). Thus a complete analysis of the solution is a complicated task. However, since we have found that the higher harmonic terms are of minor importance, we may concentrate on the Pi1 terms. Moreover, we may focus here only on those items which contain the response function of the ground wave. With these simplifications in mind, we see that each iteration grows by a factor of fi = cisx(k)3E0(k)2. Here, ci is a number which varies in each iteration step. For a fixed value of s the task is to give an upper limit for the external field, E max so that fi B1 for all iteration steps i. In 0 order to get a criterion for all values of k, we may write + E max 0 =

'

1

c

max

s[x max (k)]3 2

a 32[l 2 −(L−l)2]3 c maxsl 6

(20) (21)

(k) is the maximum value of the linear where x max 2 response function x2(k), which may be expressed in terms of the parameters a2, l, L, l [13]. The factor c max is the maximum limit of the ci s. From Eq. (20) we see that the convergence criterion for the P-series is equivalent to the saturation of the ci upon successive iteration. If ci grows with the increasing number of iterations, limi “ E max =0 and the whole perturbation 0 approach would fail for any finite value of E0. To study the behavior of ci, we have calculated the first seven iteration steps of the perturbation series, i.e. solutions for the terms  E 15 0 . The set of ci s, shown in Table 1, demonstrate that dc =ci + 1 −ci decreases. However, this is not a sufficient criterion for convergence. Therefore, we tried to fit the ci to three different types of functions: a harmonic series (c (1) i ) =c1 + S n = 1i 6/n), a logarithmic function (c (2) =c ˜ log(b0 i )) i 1 +a and an exponential type function c (3) =c ˜ (1− i 1 +a (2) exp(b0 (i− 1))). The expressions for c (1) and c are i i divergent, while c (3) is convergent. We compare these i

4. Comparison with computer simulation In order to check the validity of the perturbation theory approach, we have performed a series of computer simulations for a system of 200 water molecules. The system geometry consists of a cubic box of side ˚ , to which periodic boundary condilength L= 18.15 A tions were applied. The mass density was r= 0.999 g cm − 3 and the simulations were performed for a temperature of 300 K. We used the flexible, BJH water model [28,29], which does not include electronic degrees of freedom. Due to the flexibility of the molecules, the time step of integration has to be small, Dt = 0.25 fs, compared to a simulation performed for rigid molecules. The total trajectory length differed between 100 and 500 ps, depending on the convergence of the mean value of polarization. In the computer simulation we followed the same strategy as we did in the perturbation approach. We imposed an external electric field on the water system, which had the form: E0(r) =k. E0 cos kr. In order to calculate the response function, we varied the magnitude of the field and its wavelength and detected the system polarization. Since this method requires a separate simulation run for each field strength and each wave number, this is a very CPU-time consuming task. We, therefore, performed simulations only for five selected k values, the characteristic points of the linear response spectrum, with the constraint kn = 2p n /L (n=(nx,ny,nz ) is an integer ˚ − 1 (the smallest number vector). These are: k1 = 0.34 A ˚ − 1, k value accessible in the simulation) k5 = 1.73 A Table 1 The parameters ci as a function of the iteration step i of the perturbation series Iteration step

ci

dc (1) i

dc (2) i

dc (3) i

1 2 3 4 5 6 7

3.00 9.00 12.00 13.75 14.89 15.69 16.26

0.00 0.00 0.00 0.25 0.61 1.01 1.44

0 0.02 0.54 1.29 2.09 2.87 3.64

0.00 −0.64 −0.29 0.06 0.23 0.25 0.19

Also shown are the dc (i j) =ci−c (i j), where the index j refers to the harmonic sum ( j =1), logarithmic function ( j=2) and the exponential function ( j =3) approximation. Note, that the divergent harmonic sum gives exact correspondence with the calculation up to i=3, from where it takes larger values than ci. This may be taken as a criterion for convergence.

A.A. Kornyshe6, G. Sutmann / Journal of Electroanalytical Chemistry 450 (1998) 143–156

rb (k) =

1

N

% qa eik(rj + draj )

147

(22)

aj

where N is the number of the molecules, qa is the partial charge of site a, rj is the center of mass position of molecule j and draj is the position of site a of molecule j. The bound charge density is directly related to the longitudinal part of the polarization k. P=rb (k)/ik, from where we may calculate the response function. This simple description of the polarization is strictly valid only in the linear case, since in the non-linear case the deconvolution of the longitudinal and transversal parts breaks down. Nevertheless, one may study the transition from linear to non-linear behavior and study the k dependence of the strength of the critical field, after which the response becomes non-linear. From our simulation data, we find that the bound charge density pattern as a function of the external field, may be well described by a Langevin function [27] Fig. 1. Comparison of the calculated values ci with the harmonic (2) series (c (1) with a˜ =1, b0 = 20) and the i ), the logarithmic type (c i (3) exponential type (c i with a˜ = 14.3, b0 = 0.47) approximation. From this graphical representation, convergence of the calculated ci, and, consequently, convergence of the expansion of P, Eq. (12), is evident.

Fig. 2. Comparison of the mean-field theory result for x(k, E0) and ˚ 3/eV. the computer simulation results for a parameter s= 5×10 − 4 A

˚ − 1 (the main maximum position), k14 = 4.84 k9 = 3.11 A −1 ˚ A (the satellite maximum position) and k20 = 6.92 ˚ − 1. For each k point we applied between 8 and 10 A different electric field strengths within the interval ˚ , corresponding to E0  [0.08; 0.005B E0 B5 V/A 80]kBT/m, where m= 2D is the average dipole moment of the water molecule. From the simulation data, we then calculated the Fourier components of the bound charge density distribution



rb (k)= Ak coth(BkE0)−



1 BkE0

(23)

where Ak and Bk are the k dependent fitting parameters. The non-linear behavior first starts at the resonance ˚ −1), position of the linear response function (k9 = 3.11 A ˚ . For fields as large as at field strengths of E0 : 0.05 V/A ˚ , the maximum is reduced to some 20% E0 = 1 V/A with respect to the linear response. On the other hand, the smallest k value remains responding linearly up to ˚ . For the case of water, our field strengths of E0 : 1 V/A analytical model thus predicts the correct behavior of the response. In Fig. 2, we compare the results of the mean field theory and the simulation data. We see here a nearly perfect agreement within the error bars of the simulation data [27]. Since the analytical solution is based on a perturbation approach, one cannot expect to reproduce the simulation data over the whole range of the external field strengths. For a parameter value of s= 5× 10 − 4 the theory reproduces the simulation data ˚ . For larger values of E0, up to field strengths of 0.1 V/A the series solution of P2 would not converge. It should be mentioned, however, that the phenomenological Hamiltonian, Eq. (3), was not designed to reproduce the high k fine structure of the response function, like ˚ − 1) in the case the satellite maximum of x(k) (k:5 A for water. Therefore, at high k our theory cannot give reliable results (for a discussion see also Ref. [13]) and we compared only the first three k points from the simulation data with our theoretical predictions. Until now we compared the simulation data with results coming from the perturbation theory approach. This is naturally limited to small values of the external field strength. In order to get a description of the response in the complete interval of external fields one may try to approximate the simulation data by Eqs. (10,19) with an appropriate choice of the parameters.

A.A. Kornyshe6, G. Sutmann / Journal of Electroanalytical Chemistry 450 (1998) 143–156

148

Fig. 3. Comparison between the computer simulation results for the ground wave response function, x(k, E0) and predictions of the closed form expression for the field dependent response function, Eq. (24). The reduction of the maximum in x(k, E0) with increasing field strength is parametrized via the field dependent coupling parameter L(E0) (shown in the inset) in Eq. (24).

Since the strongest non-linear effects were found in the ˚ − 1) one may concentrate on resonance k range (k \1 A Eq. (19) which describes the response here. It turns out that it is possible to cover the whole interval of external field strengths, by varying only the parameter L in Eq. (19), responsible for the coupling between polarization and density fluctuations. Small deviations are observed only for the highest field. Moreover, we find that an effective value of L depends linearly on the external field strength in the non-linear response regime, i.e. L(E0) = (1−u(E0 − E0,l ))L

The insert shows the field dependence of L. This description of the non-linear dielectric response has no theoretical justification. Nevertheless it might be useful to have an empirical, closed form expression for the deviation of the dielectric response from the linear behavior at hand3. We shall use it in the next section in order to calculate the potential and hydration Gibbs energies of ions in solution.

+u(E0 −E0,l )(L −0.035(E0 −E0,l )), where ˚ . In u(x) is a stepwise function and E0,l =0.05 V/A Fig. 3 the simulation results are plotted against the fitted field dependent response function

3 The parameters of the phenomenological theory cannot be directly expressed through the microscopic molecular parameters. They could, however, be fitted to the results of the molecular dynamics simulation (which is done is this paper) or to any kind of analytical microscopic theory. Since their results for the response functions are usually very sensitive to the model of the charge distribution inside the molecule [12,38], and the latter are usually quite crude, they cannot be considered as the truth in the last instance. Therefore, by slightly varying the parameters of the phenomenological theory, one may account for modifications of the model intramolecular charge form factors. This gives a certain flexibility to the formulae of the phenomenological theory.

x(k,E0)=





1 1 1 − e  e 1 + l 21k 2 1 1 1 + l 2k 2 + − e e (1 + l 22k 2)(1 +l 2k 2) − L(E0)2k 2 (24)





A.A. Kornyshe6, G. Sutmann / Journal of Electroanalytical Chemistry 450 (1998) 143–156

Born sphere model (BS)

5. Applications

sin(ka) ka

5.1. Potential of ions in solution

rBS (k)= q

The deviation from Coulomb’s law of the electrostatic potential of ions may be characterized by the screening function, S(r), which is defined as

smeared Born sphere model (SBS)

S(r) = F(r)er

=

(25)

where F(r) is the electrostatic potential of an ion centered at the coordinate origin and e is the macroscopic dielectric constant. For the evaluation of the electrostatic potential one may start with the expression F(r) = Fe (r) +4p

&

dr%

9P(r%) r −r%

(26)

where Fe is the external potential, defined through E0(r) = − 9Fe (r) and the second item of the rhs is the potential due to the medium polarization field P(r). The results of Section 3.1 indicate that the terms, containing higher harmonics to the ground wave of wave number k give only small contributions to the medium polarization. As an approximation, one may neglect these terms in Eqs. (14) – (18) and consider only those terms which are proportional to x(k). We then write the polarization as P = P0 +dP[E0], where P0 is the linear part and dP is the non-linear correction, depending on the external field strength. In order to calculate F(r) via Eq. (26), we perform a Fourier transformation of its rhs4 which gives F(r) = =

2 p

&

1 2p 2

&

149

dk ikr [r (k) + r (0) b (k) + drb (k)]e k2 e

dk[1−xeff(k,E0)]re (k)

sin(kr) kr

rSBS (k)

(29)

 

q sin(ka) 1 − a/h 1+ (h/a) (2−e ka 1+ h 2k 2 ) 2

+

h 2 2 cos(ka)− e − a/h a2 (1+h 2k 2)2

(30)

In the BS model it is assumed that the charge is uniformly distributed over a thin spherical layer with radius a (the radius of the ion). The SBS model assumes that the charge density has a maximum at a distance a from the center of the ion and that it decreases from there exponentially with a decay length h. For h=0, the SBS model reduces to the BS model. Since we have found that the non-linear effects are strongest at the resonance position of the linear response function, we expect for the electrostatic potential that a deviation from the linear case will strongly depend on the form factor re (k) of the ion. If re (k) has large contributions at the resonance position, the difference between the linear and non-linear case will be pronounced. For an illustration let us compare the Fourier representation of the BS and the SBS model in ˚ crosses Fig. 4. The BS form factor with radius a= 1 A zero at the maximum k position of x(k), while for ˚ the form factor has a considerable value at a= 1.5 A

(27) (28)

In Eq. (27) we have used the solution of the Poisson equation for the external charge density, (DFe (r)= 4pre (r)), and have defined a perturbation charge density, which results from the perturbation polarization (9dP[E0]= − drb (r)). In Eq. (28) we then introduced the effective response function, xeff(k, E0) = P(k)/E0(k), which contains non-linear corrections with respect to x(k). The quasi-linear equation, Eq. (28), gives the possibility to study the non-linear response of the solvent to a given ion charge distribution, characterized by its Fourier transform, re (k) In so doing we will compare the Born model of an ion with a smeared Born sphere model, studied for the linear case in Refs. [30,31]: 4

Since we have neglected all higher harmonic terms in the series solution of the perturbation approach, we consider only the response to the ground wave k, for which we may perform a Fourier transformation of the rhs of Eq. (26). This procedure may be called a quasi-linear approximation.

Fig. 4. The charge density form factors of the BS and the SBS ˚ and a = 1.5 A ˚ . The magnitude of models, compared for radii a = 1 A r(k) at the maximum position of the response function, x(k) varies with radius. This gives rise to different overlaps between r(k) and x(k) in the corresponding integrals.

150

A.A. Kornyshe6, G. Sutmann / Journal of Electroanalytical Chemistry 450 (1998) 143–156

Fig. 5. Screening functions calculated with the help of the linear response approximation and with the non-linear response function based on the analytical expression of the perturbation approach, ˚ and a parameter s = 5×10 − 4 applied for a field of E0 : 0.1 V/A ˚ 3/eV. Compared are results for BS and SBS models of radii a =1.5 A ˚ and a =1.5 A ˚ . For clarity, the results for the different radii are A ˚ ion, shifted with respect to each other. Left axis refers to the 1.5 A ˚ ion. right axis corresponds to the 1 A

this k. The SBS model shows the general tendency to smear out the form factors in real space and to diminish the contributions at high k in Fourier space. This leads to smaller contributions at the maximum position with respect to the BS model. When the ion is imbedded into the medium, the corrections to the linear response come not only from large electric fields but also due to the disturbance of the medium by its finite size. This effect is studied in Appendix A within the cut-out approximation for the BS- and SBS-models. Let us first discuss qualitatively the effect of the appearance of non-linearity in the response function. As a crude assumption we approximate xeff(k, E0) in Eq. (28) with the response function which results from perturbation theory for the maximum allowed external ˚ ). field strength (E0 :0.1 V/A The screening functions S(r), defined by Eq. (25), are compared in Fig. 5 for the BS and SBS models. A common feature of the screening functions are the

decaying oscillations around zero, the signature of the overscreening effect [30,32,33]. With increase of the ion radius their amplitudes are diminished. This is because the wave numbers which give prominent contributions to the integral in Eq. (28) become smaller, i.e. they move away from the resonance. There is also another, non-trivial consequence of the ion size dependence of the form factors. For instance, the differences in the screening functions, calculated on the basis of the linear and non-linear response, are small for an ion with ˚ . The only difference is that the ampliradius a=1 A tudes of oscillations are slightly decreased in the latter case. The effect of the non-linear response is stronger ˚ . The amplitudes of the for an ion with radius a=1.5 A oscillations are reduced more strongly and damped faster with respect to the linear response calculation. This is due to a stronger overlap between the ion form factor and the solvent response function (re (k) which gives a profound contribution around the resonance position of x(k)). Since the maximum value of xeff(k,E0) is reduced with respect to x(k), the effect appears to be ˚ ion than for the 1 A ˚ ion. For an stronger for the 1.5 A ˚ ion of 2 A radius, the effect is weaker than for the 1.5 ˚ ion because of the weaker overlap. A The effect of a smeared charge distribution is to reduce the oscillation amplitudes with respect to the BS model. This is because of the reduction of the high k components of the ion form factor, which reduces the overlap with the response function at the resonance k position. Since this effect is qualitatively the same for ions of different sizes, we observe no new qualitative features between linear and non-linear calculation for ˚ and 1.5 A ˚ ion. the smeared 1 A Thus the non-linear effects must not necessarily diminish monotonically with the increase of the ion size. We now discuss the effect of a variation of the electric field with distance from the center of an ion. ˚ [34], For realistic ions, e.g. Cl − with radius a= 1.64 A the field strength at the ionic surface in vacuum is ˚ , which must give rise to a strong non-linear :4.5 V/A response of the solvent in the first hydration shell. As an approximation, we insert the field dependent response function, Eq. (24) into Eq. (28), where now the parameter L is varied according to the radial field strength dependence of the ion5. A smeared charge 5 This is, of course, a crude assumption for the nonlocal response function. Indeed, in the non-linear case xab does not only depend on the relative distance between two coordinates, but also on their absolute positions, i.e. xab (r,r%)=xab (r −r%,(r+r%)/2)"xab ( r−r% ). A Fourier transformation of the second item leads to xab (k,(r+ r%)/2). The approximation therefore reduces the tensor to its longitudinal part, kakb /k 2xab (this was already assumed in the perturbation approach and in the computer simulation) and it assumes a linear relation between the distance r and the external field and does not distinguish between correlations of particles, the connecting vector of which (r− r%), has different orientations with respect to the distance vector (r +r%)/2.

A.A. Kornyshe6, G. Sutmann / Journal of Electroanalytical Chemistry 450 (1998) 143–156

distribution will exhibit a smaller field at the ionic surface than a BS model. This will be mirrored in a different response around the ion. As a consequence, the effective response function will depend on the charge distribution model, re (k) If we use Eq. (28) with these approximation we find a strong difference in the screening function with respect to the linear case (see Fig. 6). The overscreening effect is obviously removed. The oscillations of S(r) take only positive values, the amplitudes of which are reduced with respect to the linear case. This damping of the oscillations is due to the broadening of the k-space resonance of the effective response function, Eq. (24), with the decrease of L. The position of its maximum remains unchanged under the influence of an external field but the magnitude is drastically decreased for large fields. Consequently the spatial wavelength, observed in S(r) is not altered, but the spatial damping is strong, leading to an approximately exponential decrease of S(r) which is due to the Lorentzian x1(k) contribution to the response function. There are still substantial differences between linear and ˚ from the ion. non-linear results at a distance of 5 A This would imply that even fluctuations in the second hydration shell would be affected considerably by the ion. We stress, however, that this result can only point out a qualitative trend of the non-linear effects, since we have neglected all offdiagonal terms of the response tensor. Nevertheless, it shows that non-linear effects may considerably change the results obtained within the linear response approximation.

151

5.2. Gibbs energy of sol6ation We now consider the Gibbs energy of solvation which we call for brevity solvation energy. It is defined as the difference in the reversible work which has to be performed to create a charge distribution in vacuo and in the infinitely diluted solution. This work may be described by three main terms Ws = Wcav + Wdisp + DWel

(31)

Here Wcav is the work to create a cavity in the solvent. This energy may be approximated either by the enthalpy of evaporation of the solvent or the surface tension of the solvent multiplied by the cavity surface area. The latter approximation should work especially for large cavities. The term Wdisp stands for the energy of dispersion forces which contains the work to put the ion into the cavity without charging it. This may be approximated, e.g. by a Lenard-Jones type 6-12 potential. Finally DWel = W0 − Wmed is the difference in work which has to be performed to charge the ion in vacuo (W0) and in the dielectric medium (Wmed). W0 may be easily obtained from classical electrostatics. The second term may be calculated within a charging cycle, where the charge distribution, e.g. a BS or SBS model of an ion with radius a, is charged from 0 to q. Taking into account the spatial dispersion and the non-linear response of the solvent we may write the solvation energy as DWel =

&

1

dz z

0

&



dr 4pr 2r˜ (r)

0

2 p

&



dkxeff(k,zE0)r˜ (k)

0

sin kr kr

(32)

where r˜ (r) contains a correction to the charge distribution model due to the finite size of the ion (see Appendix A). For the case of a BS model, i.e. r(r)= d(r− a)/4pa 2, Eq. (32) may be written in a form of a modified Dogonadze-Kornyshev equation [35] DWel =

q2 p

&



dkxeff(k,E0)

0

&

sin2 ka k 2a 2

(33)

where q is the charge of the ion and x˜ eff(k,E0)= 2

1

dz zxeff(k,zE0)

(34)

0

Fig. 6. Screening functions for the BS and SBS models calculated on the basis of the linear response approximation and on the quasi-linear approximation, Eq. (28), which takes into account the distance dependence of the ionic field strength. We also show the effect of a smearing of the ionic charge-distribution.

and xeff(k,zE0) contains the electric field dependence of the response function in the course of the charging process. For an illustration, the field dependent response function, x˜ eff(k,E0), is shown in Fig. 7. Results, based on Eq. (32) for the SBS model and on Eq. (33) for the BS model, are shown in Fig. 8. They are compared to experimental values for singly charged ions [36] and to results, based on (i) the linear response

152

A.A. Kornyshe6, G. Sutmann / Journal of Electroanalytical Chemistry 450 (1998) 143–156

Fig. 7. The effective response function x˜ eff(k,E0), Eq. (34), which takes into account a variation of the field strength with distance from the ion. ˚ from the ion. x˜ eff(k,E0) reaches the linear response function at a distance r :15 A

function (L= const in Eq. (24)), (ii) a calculation with a Lorentzian response function (e  = e =1 and l1 = 3 ˚ in Eq. (24)), and (iii) a local electrostatic calculation, A using the Born formula [37]. As was already found in Refs. [30,38], the use of the linear response function, ˚ − 1, leads to a exhibiting a large maximum at k=3 A huge overestimation of the solvation energies, when applied to BS models. The resulting electrostatic part of the solvation energy for small ions is 3 times larger than the Born formula prescribes, which already predicts larger values, than found from experiment. In contrast, the Lorentzian type approximation underestimates the solvation energies by a factor of 3 for small ions. Since the Lorentzian interpolates the solvent dielectric function from the vacuum value at the ionic surface to the bulk dielectric constant e, it always predicts smaller solvation energies than the classical Born formula. Evaluating DWel with the use of the non-linear response function, x(k, E0), the results come close to the experimental values. Due to the strong electric field in the nearest surroundings of the ion, the resonance maximum in the response function is drastically reduced and consequently the overlap between ion form factor and response function in Fourier k space do not lead to the strong enhancement of DWel, as was found in the linear calculation. The non-linear calculation still overestimates the solvation energy slightly for larger ions. This could be due to the non-electrostatic parts Wcav and Wdisp, which will lead to a reduction of Ws. However, estimates based on noble gas solvation in water lead to corrections of only 0.3 eV which are too small to explain the difference [39].

The results discussed so far do not include the effect of a continuous charge distribution. Calculations, based on the SBS model show that the solvation energies may be reduced to values, comparable with experiment, ˚ for when assuming a smearing parameter h : 0.3 A large ions, which is a reasonable value for the delocalization of the ion excess charge. From Fig. 8 it is seen that DWel first increases when the ion charge is slightly smeared. We mention that with the introduction of the continuous charge distribution two opposite effects are balanced. On the one hand, smearing of the charge distribution leads to a reduction of the overlap between ion form factor and response function which would diminish DWel. On the other hand, since the electric field in the nearest surroundings of the ion is reduced the resonance maximum will be increased with respect to the BS model, which will lead to an increase of DWel. Thus the non-monotonic behavior of DWel with increasing h is not surprising. Note, also, that the values of h can vary from ion to ion, which may improve the agreement with experimental data. From the inset of Fig. 8, it follows that the increase of h with the radius of the ion will fully reproduce the experimental data.

6. Conclusions Are our results discouraging for the linear response theory if it is applied to electrostatic calculations of ionic fields in water? The Coulomb field of a singly ˚ radius is of the order of 1 V/A ˚ in charged ion of 2 A the first solvation shell. As we find from our computer

A.A. Kornyshe6, G. Sutmann / Journal of Electroanalytical Chemistry 450 (1998) 143–156

153

Fig. 8. Electrostatic parts of the hydration Gibbs energy for singly charged ions. The results from the non-linear calculation (solid line), based on Eq. (24), are compared to (i) experimental data ( anions, “ cations) [36], (ii) linear response approximation (L= const in Eq. (24)) (dashed line), ˚ (dotted line) and the dispersionless limit, i.e. Born and (iii) linear response calculation using a Lorentzian type response function with l = 3 A equation [37] (dash-dotted line). In the inset we show the effect of smearing the external charge distribution (SBS model) on the nonlinear result.

simulation for small k, non-linear effects start to become important from these field strengths. Thus an important criterion for the applicability of linear response theory to ion solvation will be the form factor of the ions, i.e. how the ion structure overlaps with the solvent structure and what local electric field strengths are produced by the ion. If the form factor has prominent contributions at the resonance k position, linear theory results for screening factors and interaction energies might become seriously inaccurate. As it was discussed in Ref. [31], the coupling between ion and solvent becomes weaker at the resonance k if the charge distribution of the ion is smeared out spatially, i.e. when it is not considered as a point charge or an infinitely thin charge layer as it is in the Born sphere model. Our calculations show that the hydration energies are too large when the ion is approximated by a Born sphere. Introducing a continuous charge distribution, both results based on a linear and a non-linear response function could be matched with

experimental hydration energies. Thus, not only the non-linear response of the solvent might be the crucial issue, but also the charge distribution of the solute. Having these considerations in mind, the application of linear response may still be a good approximation in some cases. In the present work we considered only the contribution of the diagonal elements of the response tensor. How much the off-diagonal elements will affect the results is still an open question, which should be explored in further simulations.

Acknowledgements We are grateful to S. Leikin and Ph.A. Bopp for useful and encouraging discussions. The financial support of Deutsche Forschungsgemeinschaft through grant Sti 74/3-2 is gratefully acknowledged.

A.A. Kornyshe6, G. Sutmann / Journal of Electroanalytical Chemistry 450 (1998) 143–156

154

Appendix A. The effect of a cut-out volume (‘The brass ball in a bathtub’) We derive here the linear response equations of the cut-out-approximation (COA), which corrects an immersed charge approximation (ICA). The latter assumes that the ionic charge distribution is immersed into the solvent not perturbing its structure. Within COA it is assumed that the solvent is excluded from the interior of a certain volume which is characteristic for the charge distribution geometry, e.g. the sphere of radius a for the BS model. Within COA the electrostatic part of the solvation energy of a singly charged ion, DWel =

1 8p

&

dr dr% xab (r,r%)Da (r)Db (r%)

(A1)

can be written as DWel =

1 8p

&

D0 (r) = D(r)u(r−a) as DWel = where D0 (k) =

1 1 8p (2p)3

&

&

dk x(k) D0 (k) 2

drD(r)e − ikr

(A3)

(A4)

(A5)

r \ a

Simple algebra shows that for a spherical cut-out volume with radius a D0 (k) = D(k) + dQD(Q)F( Q − k ,a)

&

(A6)

= D(k) + dD(k)

(A7)

The function F(x) is given by

dr dr% xab ( r − r% )Da (r)Db (r%)u(r− a)

u(r% −a)

Eq. (A2) may be rewritten in terms of a modified, effective induction

(A2)

F(x)=

1 [ j (x)+ j2(x)] 6p 2 0

(A8)

Fig. 9. Modified form factors of the SBS model due to the cut-out-approximation. The modified form factor is compared to the native r(k) of the SBS and BS models. The inset shows the real space analogue of the cut-out modified form factor.

A.A. Kornyshe6, G. Sutmann / Journal of Electroanalytical Chemistry 450 (1998) 143–156

where the jl (x) are spherical Bessel functions of order l [ j0 (x) =sin x/x, j2(x) = 3/x 2(sin x/x− cos x) −sin x/ x]. The correction term, dD(k), can be simplified further. Since it is only the longitudinal component of the induction which appears in electrostatics, we may write D(k) =k. D(k)

(A9)

where we have used the unit vector notation xˆ =x/ x . After integration over the azimuthal angle one finds

&

dD(k)= a 3 dqD(q)F(a q − k ,1)qˆk. = 4pa 3

&

(A10)



dqq 2D(q)f(a, q, k)

(A11)

0

f(a, q, k)=

1 2

&

(A12)

4p r(k) ik

(A13)

Following Eq. (A11) we may then establish the corresponding correction r˜ (k) = r(k) −dr(k), where dr(k)=4pka 3 =

2a p

!

&

&



dq q r(q)f(a, q, k)



(A18)

The resulting hydration energy may then be estimated via DWel :

q0 qi DW (SBS) + DW (BS) el el q q

(A19)

For SBS-models where the smearing parameter, h\0, is fixed to be the same within and outside the cut-out sphere with radius a, it is always q0 \ qi , i.e. the result will be closer to DW (SBS) , rather than to DW (BS) el el . Nevertheless, q0 = q0(h) and qi = qi (h), is sensible on the form of r(r), so that the corrected result will always depend on the special choice of the external charge density model.

References

dq r(q)

0

"

Thus, dD(k)= 4p dr(k)/(ik), where dr(k) is given by Eq. (A15) and the resulting expression for the electrostatic part of the solvation energy may be written as DWel =

&

q0 = q− qi

(A14)

(A15)

1 p

dr r 2 r(r),

0

q 2j 0(aq) cos ak−k 2j0(ak) cos aq −j0(aq)j0(ak) q2−k2

=

a

−1

It is often simpler to characterize the external charge density distribution by its form factor than by its induction field. These two quantities are related by the equation D(k) =

&

0

+1

dtF(a q 2 +k 2 −2qkt)t

energies with respect to the BS model. For the case of water, where the dielectric response function has a resonance like shape in k-space, the BS model leads to DWel which is a factor : 4 too large for small ions compared with experimental data. This tendency still holds for the nonlinear case (c.f. Fig. 8): for a fixed value of a, the SBS model leads to smaller DWel than the BS model. The correction term, dr, leads to DWel which is always between the results for the SBS (DW (SBS) ) and for the BS(DW (SBS) model. This becomes el el clear from the following considerations. The total charge quantity, q, may be subdivided into the regions within and outside the cut-out volume (c.f. Fig. 9), i.e. qi = 4p q

where we have introduced

155

1 1 (2p)3 2

&

dk k 2x(k)

( D(k) 2 + dD(k) 2 −2 D(k) dD(k) )

(A16)

dk x(k)(r(k)2 +dr(k)2 −2r(k)dr(k))

(A17)

where dD(k) is given by Eq. (A11) and dr(k) by Eq. (A15). In order to illustrate the procedure, outlined above, we show dr in Fig. 9 for the SBS-model of an ion, where it is assumed that the radius of the cut-out volume is a (c.f. Eq. (3)). In Ref. [30] it was shown that within the linear response approximation the SBS model reduces (depending on the smearing parameter h) the solvation

[1] M. Neumann, Mol. Phys. 57 (1986) 97. [2] M.S. Skaf, T. Fonseca, B.M. Ladanyi, J. Chem. Phys. 98 (1993) 8929. [3] D. Bertolini, A. Tani, Mol. Phys. 75 (1992) 1065. [4] A. Chandra, B. Bagchi, J. Chem. Phys. 90 (1989) 1832. [5] T. Fonseca, B.M. Ladanyi, J. Chem. Phys. 93 (1990) 8148. [6] P. Attard, D. Wei, G.N. Patey, Chem. Phys. Lett. 172 (1990) 69. [7] F.O. Raineri, Y. Zhou, H.L. Friedman, Chem. Phys. 152 (1991) 201. [8] F.O. Raineri, H. Resat, H.L. Friedman, J. Chem. Phys. 96 (1992) 3068. [9] M. Vossen, F. Forstmann, J. Chem. Phys. 101 (1994) 2379. [10] A.D. Trokhymchuk, M.F. Holovko, E. Spohr, K. Heinzinger, Mol. Phys. 77 (1992) 903. [11] P.A. Bopp, A.A. Kornyshev, G. Sutmann, Phys. Rev. Lett. 76 (1996) 1280. [12] P.A. Bopp, A.A. Kornyshev, G. Sutmann, J. Chem. Phys. (1998) in press. [13] A.A. Kornyshev, S. Leikin, G. Sutmann, Electrochim. Acta 42 (1997) 849. [14] D.C. Grahame, J. Chem. Phys. 18 (1950) 903. [15] F. Booth, J. Chem. Phys. 19 (1951) 391. [16] R.M. Noyes, J. Am. Chem. Soc. 84 (1962) 513. [17] E. Glueckauf, Trans. Faraday Soc. 60 (1964) 572. [18] K.J. Laidler, J.S. Muirhead-Gould, Trans. Faraday Soc. 63 (1967) 953. [19] J.S. Muirhead-Gould, K.J. Laidler, Trans. Faraday Soc. 63 (1967) 958.

156

A.A. Kornyshe6, G. Sutmann / Journal of Electroanalytical Chemistry 450 (1998) 143–156

[20] H. Block, S.M. Walker, Chem. Phys. Lett. 19 (1973) 363. [21] M. Abraham, J. Liszi, L. Meszaros, J. Chem. Phys. 70 (1979) 2491. [22] J. Liszi, L. Meszaros, I. Ruff, J. Chem. Phys. 74 (1981) 6896. [23] J. Liszi, I. Ruff, in: R.R. Dogonadze, E. Kalman, A.A. Kornyshev, J. Ulstrup (Eds.), The Chemical Physics of Solvation, part A, ch. 4, Elsevier, Amsterdam, 1985, p. 119. [24] E.W. Castner, G.R. Fleming, B. Bagchi, M. Maroncelli, J. Chem. Phys. 89 (1988) 3519. [25] Y.I. Kharkatz, A.A. Kornyshev, M.A. Vorotyntsev, JCS Faraday II 72 (1976) 361. [26] Y.A. Il’inskii, L.V. Keldysh, Electromagnetic Response of Material Media, Plenum Press, New York, 1994. [27] A.A. Kornyshev, G. Sutmann, Phys. Rev. Lett. 79 (1997) 3435. [28] P.A. Bopp, K. Jancso, K. Heinzinger, Chem. Phys. Lett. 98 (1983) 129. [29] G. Jancso, P.A. Bopp, K. Heinzinger, Chem. Phys. 85 (1984) 377.

.

[30] A.A. Kornyshev, G. Sutmann, J. Chem. Phys. 104 (1996) 1524. [31] A.A. Kornyshev, G. Sutmann, Electrochim. Acta 42 (1997) 2801. [32] A. Fasolino, M. Parrinello, M.P. Tosi, Phys. Lett. 66A (1978) 119. [33] M.P. Tosi, in: M.P. Tosi, A.A. Kornyshev (Eds.), Condensed Matter Physics Aspects of Electrochemistry, World Scientific, Singapore, 1991, p. 68. [34] B.S. Gourary, F.J. Adrian, Solid State Phys. 10 (1960) 127. [35] R.R. Dogonadze, A.A. Kornyshev, JCS Faraday II 70 (1974) 1121. [36] J.E.B. Randles, Trans. Faraday Soc. 52 (1956) 1573. [37] M. Born, Z. Phys. 1 (1920) 45. [38] A.A. Kornyshev, D.A. Kossakowski, M.A. Vorotyntsev, in: M.P. Tosi, A.A. Kornyshev (Eds.) Condensed Matter Physics Aspects of Electrochemistry, World Scientific, Singapore, 1991, p. 92. [39] M.H. Abraham, J. Liszi, JCS Faraday Trans. I 74 (1977) 1604.