In the late 1970s, Hubbard and Onsager predicted that adding salt to a polar solution would result in a reduced dielectric permittivity that arises from the unexpected tendency of solvent dipoles to align opposite to the applied field. Here we develop a novel non-equilibrium molecular dynamics simulation approach to determine this decrement accurately. Using a thermodynamic consistent all-atom force field we show that for an aqueous solution containing sodium chloride around 4.8 mol/l, this effect accounts for 12% of the total dielectric permittivity. The dielectric decrement can be strikingly different if a less accurate force field for the ions is used. Using the widespread GROMOS parameters, we observe in fact an *increment* of the dielectric permittivity rather than a decrement, caused by ion pairing and introduced by a too low dispersion force.

It is well known that addition of salt to a polar liquid changes its dielectric properties, in particular by lowering its static dielectric permittivity.^{1} Much less is known about the actual physical mechanisms which contribute to the decrement. Hubbard and Onsager^{2} first proposed that the motion of ions, dragged by an external electric field, should induce a depolarization of the solvent molecules, by orienting their dipole moments opposite to the electric field. This effect would be genuinely kinetic in nature, not having any static equivalent.^{3–5} However, the existence of this effect has never been proven experimentally, because it is masked by the decrement caused by dielectric saturation. In other words, solvent molecules in the vicinity of an ion are strongly polarized by its local field,^{6–8} and do not take part in the usual dielectric relaxation process, therefore contributing to the reduction of the dielectric permittivity. Molecular dynamics simulations represent an ideal tool to investigate this effect, which was first addressed in simulations by Chandra,^{9} but so far a coherent framework for the calculation of the Hubbard-Onsager decrement has not been proposed suited for equilibrium and out-of-equilibrium simulations. We use the out-of-equilibrium approach to estimate with high accuracy the static limit of the kinetic decrement, showing that there is indeed a sizable kinetic contribution, even though the Hubbard-Onsager continuum theory overestimates it. In addition, we scrutinize directly the impact of long-living ion pairs on the dielectric spectrum of the solution. Finally, we show that the ionic current correlation results in a dielectric increment. This kinetic contribution is often disregarded, but can counterbalance the Hubbard-Onsager kinetic decrement, even in the absence of ion pairing.

Several approaches have been proposed to date to compute the dielectric spectrum ε(ω) in presence of free ions (e.g., Refs. 10,11). The common denominator of their theoretical background is the separate coupling of the dipole vector of neutral moieties and of the current of charged ones to the electrostatic scalar and vector potentials, respectively. As only the curl-free part of the electric field survives in classical simulations, all these approaches are equivalent^{12} to calculating the dielectric spectrum from that of the conductivity σ(ω) as ε(ω) = 1 + 4π*i*σ(ω)/ω, where the quantity 4π*i*σ(ω)/ω is sometimes called the generalized dielectric function and denoted by Σ(ω). The conductivity spectrum, in turn, can be computed in the linear response regime from the time correlation current^{13} σ(ω) = β⟨**j**_{T}(*t*)**j**_{T}(0)⟩_{ω}/3*V* of the total electric current **j**_{T}, generated by both partial charges on solvent molecules and by ionic charges. The symbols ⟨…⟩ and ⟨…⟩_{ω} denote the canonical average and its Fourier-Laplace transform,

The total current can be split arbitrarily into different contributions. For aqueous solutions of atomic ions, however, it is natural to decompose it into solvent contributions, **j**_{w}, and ionic contributions, **j**_{i}. The two currents stem from motion of the partial charges of water molecules and from the ionic charges, respectively. This splitting has a particularly important meaning in the context of dynamical contributions. In the linear response relation Δ**j**_{a}(ω)∝⟨**j**_{a}(*t*)**j**_{b}(0)⟩_{ω} = *C*_{ab}(ω) (where subscript *a* and *b* are completely general), the current **j**_{b} is the derivative of the dipole moment that couples to the perturbing field in the Hamiltonian. Hence, for example, in an aqueous solution of ions, the correlation *C*_{ww} represents the response of water in a fictitious case where only water molecules – and not the ions – are coupled to the electric field. Similar considerations apply to the other possible permutations *C*_{wi}, *C*_{iw}, and *C*_{ii}. Evidently, the zero-frequency limit of Δε_{iw}(ω) = 4πβ*iC*_{wi}(ω)/(3ω*V*) is precisely the first dynamic mechanism identified by Hubbard and Onsager, namely, the depolarization of water molecules due to the motion of ions. The symmetry between *C*_{wi} and *C*_{iw} is nothing else but an instance of Onsager reciprocal relations, and shows directly that the second dynamical mechanism identified by Hubbard and Onsager (namely, a reduction in ion conductivity due to the rotation of water molecules) leads to the same contribution to the static permittivity. The total kinetic decrement is therefore Δε_{HO} = 2Δε_{wi}(0). The static dielectric permittivity can be written as ε = 1 + Δε_{wT} + Δε_{iT} = 1 + Δε_{ww} + Δε_{ii} + 2Δε_{wi}, where here and in the following the absence of a frequency dependency implies that the limit to zero frequency has to be taken.

We have investigated the dielectric spectrum of aqueous solutions of sodium chloride, performing extensive MD simulations. The systems were composed of 1621 water molecules (SPC/E potential^{14}) and 160 ion pairs, corresponding to a salt concentration of about 4.7–4.9 mol/l, depending on the average volume yielded by the use of different potentials for the ions. Distinct simulations were performed using the GROMOS^{15} and the Kirkwood–Buff force field (KBFF)^{16} force fields (see Table I for the parameters). The ensemble of choice is the isothermal/isobaric one at 298 K and 1 atm, realized using the Nosè–Hoover^{17,18} and Parrinello–Rahman coupling schemes.^{19} The electrostatic interaction has been computed using the Particle Mesh Ewald method in its smooth variant^{20} with tin-foil boundary conditions, 4th order spline interpolation on a grid with 0.12 nm spacing, and a relative strength of the interaction at 0.9 nm of 10^{−5}. Both the Lennard-Jones and short-range part of the electrostatic potentials were switched smoothly to zero using a fourth-degree polynomial in the range from 0.9 to 1.2 nm. Several 8 ns long runs with different initial conditions were performed, up to a cumulative time of about 4 μs for each system. Water and ionic currents were saved to disk every time step (1 fs) for off-line analysis. Current correlations were calculated on the complete 8 ns datasets using a fast Fourier transform based approach.^{21–23} The correlations computed on each dataset were eventually averaged, and spectra were calculated using delays of up to 2 ns of the averaged correlations, as at longer delay times the noise in the correlation degraded the quality of the spectrum noticeably. No tapering or padding procedure has been used in order to avoid any bias in the spectral features, especially in those at low frequencies.

Model . | Atom . | ε (kJ/mol) . | σ_{LJ} (nm)
. | Charge (e)
. |
---|---|---|---|---|

KBFF | Na | 0.320^{a} | 0.245 | 1 |

Cl | 0.470 | 0.440 | −1 | |

GROMOS | Na | 0.062 | 0.258 | 1 |

Cl | 0.446 | 0.445 | −1 | |

SPC/E | O | 0.650 | 0.316 | −0.8476 |

H | 0.0 | 0.0 | 0.4238 |

Model . | Atom . | ε (kJ/mol) . | σ_{LJ} (nm)
. | Charge (e)
. |
---|---|---|---|---|

KBFF | Na | 0.320^{a} | 0.245 | 1 |

Cl | 0.470 | 0.440 | −1 | |

GROMOS | Na | 0.062 | 0.258 | 1 |

Cl | 0.446 | 0.445 | −1 | |

SPC/E | O | 0.650 | 0.316 | −0.8476 |

H | 0.0 | 0.0 | 0.4238 |

^{a}

Geometric mixing rules for both σ_{LJ} and ε are employed.

^{b}

The interaction energy between oxygen and sodium atoms has to be scaled by a factor 0.75 in KBFF.

At first, we investigated the solution using the GROMOS force field for Na and Cl ions,^{15} obtaining the spectra reported in Fig. 1. The spectrum has been split into the contributions coming from water *C*_{wT} and from the ions *C*_{iT}, where the index *T* again denotes the total current. In addition, the positive part of the cross correlations *C*_{wi} has also been reported. A clear, main relaxation peak of water (left panel of Fig. 1) appears at ω ≃ 78 GHz (τ = 14.6 ps), but the ionic contribution shows the appearance of an important, unexpected relaxation process at ω ≃ 14 GHz (τ = 70 ps, d.c. contribution subtracted). A simultaneous fit of the real and imaginary part of the spectra to the phenomenological Cole-Cole functional form augmented by a conductivity term,

led to the best fit parameters (amplitude Δε, relaxation time τ exponent α and conductivity σ) reported in Table II (cf. the experimental results in Ref. 24, supplementary material sample 12a for concentration *c* = 5.0105 mol/l and *T* = 298.15 K: Δε = 43.48, τ = 6.29, α = 0.801, ε_{∞} = 5.65).

. | GROMOS . | KBFF . | KBFF+pairs . |
---|---|---|---|

Δε_{wT} | 42.1 ± 0.3 | 27.4 ± 0.1 | 39.9 ± 0.3 |

Δε_{iT} | 18.0 ± 2.0 | … | 36.0 ± 4.0 |

Δε_{wi} | 6.0^{a} | −1.7 ± 0.1^{a} | 5.9^{b} |

Δε_{HO, cont} | −6.2 | −8.0 | −8.6 |

τ_{wT}/ps | 14.6 ± 0.2 | 10.6 ± 0.1 | 17.1 ± 0.3 |

τ_{iT}/ps | 70.0 ± 1 | … | 110 ± 10 |

α_{wT} | 0.820 ± 0.005 | 0.830 ± 0.003 | 0.800 ± 0.005 |

α_{iT} | 0.91 ± 0.01 | … | 0.94 ± 0.01 |

σ/(S/m) | 5.1 ± 0.1 | 10.2 ± 0.1 | 5.4 ± 0.1 |

. | GROMOS . | KBFF . | KBFF+pairs . |
---|---|---|---|

Δε_{wT} | 42.1 ± 0.3 | 27.4 ± 0.1 | 39.9 ± 0.3 |

Δε_{iT} | 18.0 ± 2.0 | … | 36.0 ± 4.0 |

Δε_{wi} | 6.0^{a} | −1.7 ± 0.1^{a} | 5.9^{b} |

Δε_{HO, cont} | −6.2 | −8.0 | −8.6 |

τ_{wT}/ps | 14.6 ± 0.2 | 10.6 ± 0.1 | 17.1 ± 0.3 |

τ_{iT}/ps | 70.0 ± 1 | … | 110 ± 10 |

α_{wT} | 0.820 ± 0.005 | 0.830 ± 0.003 | 0.800 ± 0.005 |

α_{iT} | 0.91 ± 0.01 | … | 0.94 ± 0.01 |

σ/(S/m) | 5.1 ± 0.1 | 10.2 ± 0.1 | 5.4 ± 0.1 |

^{a}

From NEMD simulation.

^{b}

Value of Σ(ω) at the smallest frequency.

The contribution coming from the cross-correlation *C*_{wi} is large and positive, in open contrast to the prediction of a kinetic *decrement* by the Hubbard-Onsager mechanism. In the continuum approximation, i.e., in the limit of infinite dilution and large ionic radius, this decrement can be estimated^{25} as Δε_{HO, cont} = −8πτ_{wT}σ(ε_{0} − ε_{∞})/3ε_{0}, where σ, ε_{0}, and ε_{∞} are the static conductivity, static permittivity, and infinite frequency dielectric permittivity of the solution, and τ_{wT} is the (Debye-like) relaxation time of the solvent. Note that the single Debye relaxation in the Hubbard-Onsager model is an approximation, as the solvent relaxation in this case would be better described by either a superposition of Debye processes or by a Cole-Cole function.^{24,26}

In the ionic contribution to the spectrum (right panel of Fig. 1), however, a feature which is compatible with a Cole-Cole process with a characteristic time of about 70 ps is present. This suggests that ion pairing is taking place on this relatively large timescale which is consistent with earlier reports that the effective attraction between solvated Na and Cl ions in the GROMOS force field is too strong.^{27} Sodium chloride, instead, is a strong electrolyte which is generally believed not to lead to ion pairing on such a timescale.^{28}

To test our hypothesis, we simulated the NaCl solution using the KBFF,^{16} which is more accurate than the GROMOS one in reproducing structural and energetic properties of solvated sodium chloride. The water contribution is qualitatively similar to the GROMOS case, but the ion contribution differs dramatically. The prominent ion relaxation process at about 14 GHz (70 ps), which characterized the GROMOS case, disappears completely, and the zero-frequency contribution becomes so small that it can hardly be estimated (see Fig. 2).

To check explicitly that the phenomenon underlying the relaxation observed in the GROMOS case is ion pairing, we performed another series of simulations using the KBFF parameters. This time, however, we created artificial ion pairs out of 50% of the ions in solution by introducing a stiff bond between the Na and Cl ions with a bond length corresponding to the first peak in the radial distribution function of the non-bonded case. The effect of persistent pairing on the dielectric relaxation of ions modeled using KBFF can be appreciated in Fig. 3. A relaxation process appears again at low frequencies (ω ≃ 9 GHz, τ = 110 ps), qualitatively very similar to the one observed using GROMOS force field, a clear indication that the process has to be attributed to ion pairing. This also shows, in particular, that it is the contact ion pair mechanism which is involved in the present case, as the first peak of the radial distribution function is located slightly closer than 0.3 nm, and therefore a water molecule cannot be accommodated in between the two ions. In this respect, KBFF proves to be superior to GROMOS not only in describing the structural properties, but also the dynamical ones, as it correctly reproduces the absence of ion pairs that are stable on time scales longer than 10–100 ps.

Nevertheless, the failure of the GROMOS force field in reproducing the spectrum of aqueous NaCl solutions has provided new important knowledge. On the one hand, this is the first direct observation of the relaxation process associated with ion pairing, with a characteristic time of about 100 ps, well separated from the water relaxation time (approximately 10 ps). On the other hand, these results show that kinetic contributions can lead not just to the Hubbard-Onsager decrement, but also to kinetic *increments*, which arise from the completely different mechanism of pairing. It is interesting to investigate the origin of the excessive pairing in the GROMOS force field. The small Lennard-Jones interaction energy of sodium (0.062 kJ/mol) suggests a weak attraction between ions. Paradoxically, this results in excessive pairing, as also the hard-core repulsion of the ion is small, allowing two oppositely charged particles to come so close that strong electrostatic attraction sets in, favoring the pairing.

The reliability of the KBFF results offers the opportunity to test quantitatively the Hubbard-Onsager continuum theory^{2,29} and the Hubbard-Colonomos-Wolynes molecular theory^{25} for the kinetic decrement. Obtaining a precise estimate of the kinetic decrement from the analysis of the spectrum *C*_{wi}(ω), however, can be computationally quite demanding, especially if the static contribution is small. In the present case, in fact, the fluctuations of the spectrum of the cross-correlation function at low frequencies did not allow us to perform a reliable extrapolation to ω = 0. The interpretation in terms of linear response of the correlation *C*_{wi} suggests, however, that it is possible to compute the decrement using a non-equilibrium MD (NEMD) approach. This consists in applying a constant (fictitious) electric field acting only on ions, and in measuring the induced polarization in water by the ion flow. This method realizes directly the thought experiment proposed by Hubbard and co-workers^{25} to explain the first dynamic mechanism for the decrement. This NEMD approach has obvious advantages over the linear-response approach: relatively high fields can be used to increase the signal-to-noise ratio, and no extrapolation procedure is needed. Indeed, a 500 ps long run with an electric field of about 0.25 V/nm is enough to determine the dynamic contribution Δε_{wi} = −1.7 with an accuracy of 5%. This term is only half of the contribution to the total dielectric decrement, which therefore amounts to −3.4. Two separate simulations with a different applied voltage were used to control the validity of linear response assumption.

The other kinetic contribution to the static permittivity is due to the ion-ion current correlation, but unlike the water-ion cross term this cannot be computed using a NEMD approach, since the ionic cloud is infinitely polarizable. As it is seen from the spectrum of *C*_{iT} in Fig. 2, the scatter of the data at low frequencies is rather large. This makes a precise extrapolation to zero frequency very difficult, but a rough estimate suggests Δε_{ii} ≃ 2−5. Such a contribution would in large part counterbalance the total kinetic decrement 2Δε_{wi}.

Coming back to the ion-water contribution, we now compare the simulation results with the available theoretical models. The continuum theory^{2,29} requires as input parameters the main dielectric relaxation time and the conductivity of the solution. Using the values computed by the MD simulation, the continuum model predicts for the KBFF case a considerably larger decrement Δε_{HO, mol} ≃ −8 (to be compared with 2Δ_{wi} = −3.4). The situation does not improve if the molecular theory developed by Hubbard, Colonomos, and Wolynes^{25} is used. For a system with many ions in solution the molecular theory estimate of the dielectric decrement is

where *μ*_{j} is the dipole moment of the *j*th water molecule, located at the relative position **r**_{ij} from the *i*th ion of radius *R*_{i} and specific conductance σ_{i} (i.e., the contribution to the total conductance stemming from the single ion current). Applying Eq. (2) to the KBFF case leads to Δε_{HO, cont} ≃ −19, which is in absolute value even larger than Δε_{HO, cont}. This is because Δε_{HO, mol} is very sensitive to the local structure of water around each ion, due to the 1/*r*_{ij} term in the statistical average.

In conclusion, we made a coherent derivation of the Green–Kubo formulas for the two mechanisms identified by Hubbard and Onsager as responsible for the kinetic dielectric decrement, and computed the different (ionic, dipolar, cross) contributions to the dielectric spectrum of sodium chloride aqueous solutions. Using a highly accurate non-equilibrium approach, we showed that the phenomenon of kinetic decrement can account for a considerable fraction of the static dielectric permittivity in salt solutions. However, both the continuum and molecular theories for the dielectric decrement fail describing the decrement at the present concentrations. The reason for this failure has to be sought in the inadequacy of the infinite dilution approximation. The high salt density screens the electrostatic interaction, therefore reducing the efficiency of the Hubbard-Onsager mechanism. This has been confirmed by a supplementary simulation at lower concentration (60 ion pairs and same number of water molecules), for which the kinetic contribution, Δε_{wi} = −2.3, is higher (in absolute value), despite the lower conductivity of the solution. In the limit of a continuum charged background, in fact, no Hubbard-Onsager mechanism can take place.

M.S. acknowledges support from the European Community's Seventh Framework Programme (FP7-PEOPLE-2012-IEF) funded under Grant No. 331932 SIDIS. S.S.K. acknowledges support from RFBR Grant Nos. mol-a 1202-31-374 and mol-a-ved 12-02-33106, from the Ministry of Science and Education of RF 2.609.2011 and, from Austrian Science Fund (FWF): START-Projekt Y 627-N27. A.A. and C.H. acknowledge support by the cluster of excellence SimTech of the University of Stuttgart. The authors thank Friedrich Kremer, Christian Schröder, and Othmar Steinhauser for useful discussions.