Kinetic dielectric decrement revisited: phenomenology of finite ion concentrations

Abstract

With the help of a recently developed non-equilibrium approach, we investigate the ionic strength dependence of the Hubbard–Onsager dielectric decrement. We compute the depolarization of water molecules caused by the motion of ions in sodium chloride solutions from the dilute regime (0.035 M) up close to the saturation concentration (4.24 M), and find that the kinetic decrement displays a strong non-monotonic behavior, in contrast to the prediction of available models. We introduce a phenomenological modification of the Hubbard–Onsager continuum theory, which takes into account the screening due to the ionic cloud at the mean-field level and, which is able to describe the kinetic decrement at high concentrations including the presence of a pronounced minimum.

---

More than thirty years ago, in what was one of the last articles written by Onsager, he and Hubbard made a captivating prediction that has eluded direct observation until now. They stated that in a saline solution, due to the motion of ions, polar solvent molecules should show a tendency to orient against any external, static electric field, in apparent contradiction with electrostatics.^1–3 According to the continuum model presented by Hubbard and Onsager, the rotational current induced in the solvent by ionic currents should generate a net solvent depolarization that survives in the zero frequency limit. As a consequence, a decrement of the static permittivity of the solution should be observed, even though the effect is purely dynamic, and as such cannot be explained in terms of molecular configurations only. Systems containing free charges can be thought as being infinitely polarizable, since the positive and negative ions in a pair can be brought at an infinite distance by a static electric field. Strictly speaking, the static dielectric permittivity of an electrolyte solution cannot be measured; instead, the real part of the dielectric spectrum at low (but finite) frequencies is used as a measure of the permittivity. To date, however, no direct experimental proof of the kinetic decrement exists, because its detection is complicated by the presence of dielectric saturation, from which it cannot be easily separated.^4–8

A quantitative picture of the kinetic contribution to the dielectric decrement is therefore key to the investigations of ion solvation properties, which rely on a correct estimate of the static contribution of the decrement.^9–11

The continuum theory of the kinetic decrement predicts that the static permittivity ε₀ of a solvent should change upon addition of salt by an amount (in CGS units) due to a subtle interplay between motion of ions and rotational orientation of the solvent molecules,

Δε_HO = −4πpστ(ε₀ − ε_∞)/ε₀, (1)

Here, σ denotes the conductivity of the solution and τ is the time constant of the Debye relaxation process that characterizes the dielectric susceptibility of the solvent, and ε_∞ is the infinite frequency dielectric constant.^1,2 The factor p can take values between 2/3 and 1, depending on the type of boundary condition at the surface of the ion (full slip and no slip, respectively). Strictly speaking, the continuum theory is valid only in the infinite dilution limit, and for large ionic radii.³

Despite these limitations, the formula for the decrement bears an enthralling elegance, and explains qualitatively the dependence of the dielectric permittivity of electrolyte solutions on their conductivity, even well within the concentrated solution regime.¹² However, the kinetic decrement is not the only effect that is expected to lower the dielectric permittivity of electrolyte solutions. The strong local electric field in the vicinity of the ions tends to polarize solvent molecules more than any external electric field in the linear regime. Such a high field saturates the dielectric response of solvent molecules next to ions, effectively reducing the dielectric permittivity of the solution. This effect depends on the salt concentration c and, implicitly, on the conductivity σ. For this reason it becomes hard, if not impossible, to separate the kinetic contribution from dielectric saturation experimentally.^4,5,7,8 This situation prevents not only a direct observation of the kinetic decrement, but also a precise evaluation of the effect of saturation⁶ and, consequently, a correct estimate of the number of solvent molecules coordinated by the ions.

Here, we use a non-equilibrium molecular dynamics approach to compute the kinetic decrement over an unprecedented wide range of concentrations, which is not accessible with conventional, equilibrium approaches.¹³ Moreover, we present a simple phenomenological theory that gives a quantitative account of the features of the kinetic decrement at higher concentrations.

The kinetic decrement can be seen as the sum of two complementary contributions: the first is the depolarization due to the moving ions, which exerts a torque on the solvent molecules; the second, more subtle effect, is a change in the imaginary part of the ion conductivity, a lag in the response of ions induced by the rotation of solvent molecules which are oriented along the external electric field. The two contributions must have exactly the same value, as a consequence of Onsager reciprocal relations. A straightforward way to see this is through the Green–Kubo expression for the kinetic decrement. The change in solvent polarization current J_p due to the ionic one J_i leads to a contribution to the conductivity spectrum , where β is the inverse thermal energy, L³ the simulation box volume, and 〈·〉 is a suitable ensemble average. This change in conductivity reflects a change in permittivity, since ε(ω) − 1 = i4πσ(ω)/ω,¹⁴ and results in the first contribution to the kinetic decrement Δε_pi = lim_ω→04πiΔσ_pi(ω)/ω. Owing to the symmetry of the current cross-correlation function, the second contribution, which originates from the action of the rotating solvent molecules on the ions, is Δε_ip = Δε_pi. The total kinetic decrement is therefore twice the first contribution, Δε = 2Δε_pi.

However elegant, the Green–Kubo expression is not very much suited for the computation of the kinetic decrement, because the signal-to-noise ratio at extreme dilutions would be too small for any practical purposes. A much more efficient way to compute Δε_pi consists instead in applying an external fictitious field E_f, that couples to the ions only, and in calculating the resulting polarization of the solvent P, so that Δε_pi = 4πP/(L³E_f).¹⁵ This is evidently the out-of-equilibrium counterpart of the Green–Kubo formula for Δε_pi, because J_i is the current that couples to the external field E_f, and J_p the one generating the polarization P.

We applied this non-equilibrium calculation to an aqueous solution of sodium chloride at 11 different salt concentrations. In our simulations we model water molecules using the three-sites SPC/E potential¹⁶ and sodium and chloride ions using the thermodynamics consistent Kirkwood–Buff potential.¹⁷ The salt concentration c varies from 0.035 to 4.24 M, keeping the water content fixed at 1621 molecules per simulation box and changing the number of salt pairs from 1 to 140. We kept constant temperature (300 K) and pressure (1 atm) using the Nosé–Hoover^18,19 and Parrinello–Rahman²⁰ algorithms with relaxation times of 5 ps. Electrostatic interactions were computed using the smooth particle mesh Ewald method²¹ with tin–foil boundary conditions, a 4th order spline interpolation on a grid with spacing not larger than 0.12 nm and a relative interaction strength of 10⁻⁵ at 0.9 nm. We switched the short range part of the electrostatic interaction and of the Lennard-Jones potential smoothly to zero between 0.9 and 1.2 nm using a fourth-degree polynomial. Simulations were performed using an in-house modified version of Gromacs²² for the online calculation of the currents associated to the different species, and an integration time step of 1 fs was used.

To check that we are performing calculations in the linear regime, we applied our non-equilibrium approach for the system of 100 ion pairs using applied external fictitious fields of magnitude varying from 0.1 to 2 V nm⁻¹, sampling the resulting depolarization of water −4πP_x/L³ during 500 ps long simulation runs (Fig. 1) as a function of the applied field. The linear behavior holds even for rather high fields, and this allows us to collect meaningful statistics also for very dilute solutions with relatively short simulation runs, making this non-equilibrium approach the key to calculating the kinetic decrement over an unprecedented broad range of concentrations.

Fig. 1 Calculated water depolarization as a function of the applied fictitious field (points) and the linear fit (line).

Using an external fictitious field of 0.5 V nm⁻¹ we proceeded to calculate the kinetic decrement for different salt concentrations in 1 ns long runs, obtaining the values of dielectric decrement presented in Fig. 2. The kinetic decrement shows a marked non-monotonic dependence on the concentration, with a clear minimum right before c = 2 M. To test the Hubbard–Onsager formula, eqn (1), we also calculated the molar conductivity Λ of the solution at the different concentrations, as reported in Fig. 3. A best fit to the Kohlrausch law, , allows us to extrapolate the molar conductivity to infinite dilution, and estimate the limiting molar conductivity Λ₀ separately for the sodium and chloride ions. Therefore, we can write the Hubbard–Onsager decrement for the mixture of sodium and chloride in the form Δε_HO = −4πpτ(Λ^Na₀ + Λ^Cl₀)c(ε₀ − ε_∞)/ε₀. The Hubbard–Onsager decrement so calculated (Fig. 2, light line) is not compatible with the simulation data above a concentration of 0.2 M. The presence of a pronounced minimum and the subsequent increase of Δε cannot be explained even qualitatively with the continuum Hubbard–Onsager theory.

Fig. 2 Kinetic contribution to the static dielectric permittivity. Squares: simulation results; light line: the Hubbard–Onsager theory, Δε_HO(c); dark line: αΔε_DS(c), with α = 2; dashed line: αΔε_DD(c), with α = 1.64.

Fig. 3 Molar conductivity of sodium and chloride ions, as a function of the concentration. Squares: sodium; circles: chloride; solid lines: best fit to Kohlrausch law.

Physical intuition suggests that the local field of the ions, which determines the torque on water molecules, should be screened by the presence of oppositely charged ions in its vicinity. To formalize this, we introduce a mean-field correction to the Hubbard–Onsager theory, along the lines of the Debye–Hückel theory. The crucial step in the Hubbard–Onsager theory is the calculation of the rotational current J_R induced in the polar medium by an ion travelling with speed u. The coupling between electrostatics and Navier–Stokes equations allows us to express the rotational current as a functional of the local field generated by the ion, E₀, as J_R = ∫(χ/2ε₀)E₀ × (∇ × v)dr, where v is the velocity field of the solvent surrounding the ion and χ its dielectric susceptibility.¹ For large ionic radii R, the velocity field can be approximated (in the low salt concentration limit) by Stokes' solution, v(r) = (3R/4r³)[r²u + (r·u)r]. If, instead of using the Coulomb field, we use the Debye–Hückel one, E₀ = (q/ε₀r³)exp(−κr)(1 + κr)r, the rotational current can be evaluated analytically as J_R = (2π/ε₀)uχq exp(−κR). Here is the inverse Debye screening length.

Since the ratio between the ion speed u and the driving electric field E_x is u/E_x = σ/q, it is possible to express the dielectric decrement (which we denote here as Δε_DS, and where the suffix hints at the use of the Debye–Hückel field and of the Stokes' flow solution) in terms of the rotational current

Here, the imaginary part of a quantity is denoted by double-primes. The susceptibility of the dipolar medium is assumed to be characterized by a single Debye relaxation, so that 4πχ(ω) = (ε₀ − ε_∞)(1 + iωτ_D), from which one derives

The solution resembles the classical Hubbard–Onsager one, but features an additional factor exp(−κR) which depends on the (effective) ion size. This difference is an important one, because it shows that even at the mean-field level there is an additional length scale, the screening length κ⁻¹, that governs the non-monotonic behaviour of the kinetic decrement. The screening, however, affects not only the electric field generated by the ion, but also the velocity of the fluid around it, so that at finite concentrations the Debye–Hückel-corrected velocity field is expressed by²³v(r) = (3R/4r³)exp(−κr)(4/3)r²u + 2κ⁻²[1 + κr + κ²r²/3 − exp(κr)][u − 3(r·u)r/r²]. This flow can be integrated in a similar way to the Stokes' one, yielding the decrement

where the subscript underlines that the Debye–Hückel approximation has been used both for the electric field and for the flow solution.

In Fig. 2 we compare the simulation results with the mean-field results eqn (3) and (4), summed over the contributions of the two ionic species, and multiplied by a phenomenological scaling factor α, that takes into account the effect of ionic correlations arising at high salt concentration. The relaxation time τ has been computed from a fit of the Debye process χ(ω) to the spectrum of the pure solvent, and the solution conductivity has been calculated from the limiting molar ones as for the Hubbard–Onsager case. The relaxation of water is not described by a single Debye process, and often a Cole–Cole relaxation or two Debye processes²⁴ are used to fit experimental data. However, the dominant water contribution at lower frequencies comes from the main Debye relaxation, which is the one we are using here to calculate the decrement.

A very good qualitative agreement is already achieved by using the Stokes' flow based solution, eqn (3), and α = 2 (dark curve in Fig. 2), taking as a value for the effective ionic radius R the size of the first hydration shell of the ion. R is defined as the sum of the position of the first minimum in the ion–water radial distribution function and of the Lennard-Jones diameter of a water molecule. Also taking into account the screening of the flow due to the counterion clouds introduces a qualitative difference in the formula for the decrement, eqn (4), which is characterized by the presence of an additional term proportional to κ. In fact the smaller scaling factor α = 1.64 is sufficient to obtain an equivalently good fit, if at the same time the value of both effective ionic radii is increased by 1.5 Å. Both models are thus able to give a qualitatively correct picture of the concentration dependence of the kinetic dielectric decrement at finite concentrations, underlining the importance of the screening induced by the counterion clouds in determining the effective water–ion interaction in the hydrodynamic regime.

As a final remark, we note that due to the presence of the scaling factor α, both curves do not converge, for κR ≪ 1, to the solution of Hubbard and Onsager, the latter being a better approximation at low concentrations. Nevertheless, the simulation data show that the applicability range of the Hubbard–Onsager theory is limited to concentrations smaller than approximately 0.2 M, a condition which has been often not fulfilled when searching for experimental evidence of the kinetic decrement.¹² Our simulation results thus resolve the doubts which were cast on the attribution of the measured dielectric decrement⁵ in favor of the hypothesis of a static effect arising from dielectric saturation,⁷ which is definitely the largest contribution to the dielectric decrement.

M. S. acknowledges support from the European Community's the Seventh Framework Programme (FP7-PEOPLE-2012-IEF) funded under grant No. 331932 SIDIS. S. S. K. acknowledges support from RFBR grant mol-a-ved 12-02-33106, the Ministry of Science and Education of RF 2.609.2011, the URFU Stimulating Program, and, the Austrian Science Fund (FWF): START-Projekt Y 627-N27. The authors thank Christian Schröder and Othmar Steinhauser for useful discussions. A. A. would like to thank the German Research Foundation (DFG) for financial support through the Cluster of Excellence in Simulation Technology (EXC 310/1) at the University of Stuttgart.

References

1. J. B. Hubbard and L. J. Onsager, J. Chem. Phys., 1977, 67, 4850.
2. J. B. Hubbard, J. Chem. Phys., 1978, 68, 1649.
3. J. B. Hubbard, P. Colonomos and P. G. Wolynes, J. Chem. Phys., 1979, 71, 2652–2661.
4. P. Winsor IV and R. H. Cole, J. Chem. Phys., 1982, 86, 2491–2494.
5. Y.-Z. Wei, P. Chiang and S. Sridhar, J. Chem. Phys., 1992, 96, 4569.
6. K. Nörtemann, J. Hilland and U. Kaatze, J. Phys. Chem. A, 1997, 101, 6864–6869.
7. U. Kaatze, J. Solution Chem., 1997, 26, 1049–1112.
8. U. Kaatze, J. Mol. Liq., 2011, 162, 105–112.
9. K. Tielrooij, N. Garcia-Araez, M. Bonn and H. Bakker, Science, 2010, 328, 1006–1009.
10. R. Buchner, G. T. Hefter and P. M. May, J. Phys. Chem. A, 1999, 103, 1–9.
11. N. Nandi, K. Bhattacharyya and B. Bagchi, Chem. Rev., 2000, 100, 2013–2046.
12. J. Hubbard, L. Onsager, W. Van Beek and M. Mandel, Proc. Natl. Acad. Sci. U. S. A., 1977, 74, 401–404.
13. A. Chandra, J. Chem. Phys., 2000, 113, 903.
14. J. P. Hansen and I. R. McDonald, Theory of Simple Liquids, Academic Press, London, 1986.
15. M. Sega, S. Kantorovich, C. Holm and A. Arnold, J. Chem. Phys., 2014, 140, 211101.
16. H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, J. Phys. Chem., 1987, 91, 6269–6271.
17. S. Weerasinghe and P. E. Smith, J. Chem. Phys., 2003, 119, 11342–11349.
18. S. Nosé, Mol. Phys., 1984, 52, 255–268.
19. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695–1697.
20. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
21. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. Pedersen, J. Chem. Phys., 1995, 103, 8577.
22. B. Hess, C. Kutzner, D. Van Der Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435–447.
23. D. Long and A. Ajdari, Eur. Phys. J. E: Soft Matter Biol. Phys., 2001, 4, 29–32.
24. J. Barthel, K. Bachhuber, R. Buchner and H. Hetzenauer, Chem. Phys. Lett., 1990, 165, 369–373.

---