Role of Viscosity in Deviations from the Nernst-Einstein Relation.

Deviations from the Nernst-Einstein relation are commonly attributed to ion-ion correlation and ion pairing. Despite the fact that these deviations can be quantified by either experimental measurements or molecular dynamics simulations, there is no rule of thumb to tell the extent of deviations. Here, we show that deviations from the Nernst-Einstein relation are proportional to the inverse viscosity by exploring the finite-size effect on transport properties under periodic boundary conditions. This conclusion is in accord with the established experimental results of ionic liquids.

Introduction

“Physical chemistry of ionically conducting solutions” is one of the cornerstones for energy storage applications in supercapacitors and lithium-ion batteries.[1] Historically, most electrolytes were regarded as incompletely dissociated, and the dissociation constant is related to the factor α, which can be expressed as the conductivity ratio Λc/Λ0 according to Arrhenius, where Λ0 is the value at the infinite dilution.[2]

Following the idea of using transport properties to quantify the extent of ion dissociation (“ionicity”), Angell and co-workers proposed the use of the classical Walden rule for the purpose of classification.[3] In the Walden plot of log Λ versus log η, the product of the conductivity Λ and the viscosity η of KCl solution measured at 0.1 m was set as the reference point. Downward deviations from the KCl line are usually regarded as the formation of charge-neutral ion pairs.

The concept of ionicity was put forward further by Watanabe and co-workers.[4−6] They suggested the use of the molar conductivity ratio Λimp/ΛNMR measured by the impedance spectroscopy (imp) and the pulse-field gradient NMR to quantify the self-dissociativity of ionic liquids (ILs).

Despite its conceptual simplicity, the nature of ionicity is by no means simple. Apart from the static picture of charge-neutral ion pairs, other factors may alter the interpretation of the experimentally measured ratio Λimp/ΛNMR. For instance, ΛNMR was obtained via the Nernst–Einstein relation and the charge transfer effect can lead to the deviation from the formal charge of ions.[7,8] Exceptional case can also be found in which the ionicity goes up with the concentration counterintuitively where the chelate effects become important.[9]

A more general point related to ionicity is that ion pairing is a subset of ion–ion correlations.[10,11] In one of a series of classic studies on the dense ionized matter from Hansen and McDonald,[12] they commented that “It is also clear that deviations from the Nernst–Einstein relation are not necessarily the result of a permanent association of ions of opposite charge”. However, the remaining question is still what determines the deviation from the Nernst–Einstein relation if not ion pairing.

In this work, we used the finite-size effect in molecular dynamics (MD) simulation of transport properties to investigate the deviation from the Nernst–Einstein relation in the case where permanent ion pairing is excluded. It is found that while the Nernst–Einstein conductivity depends strongly on the system size, the Green–Kubo conductivity is system-size independent. We showed that these two types of conductivities crossover at certain simulation box size Lmin for both NaCl solutions and [BMIM][PF6] IL. Furthermore, this observation suggests that the deviation from the Nernst–Einstein relation, i.e., (ΛN–E – ΛG–K) is inversely proportional to the viscosity η resembling the classical Walden rule, with Lmin being a system-specific parameter. We verified this relation with published experimental data for a variety of ILs. These results indicate that viscosity is a dominating factor for the deviation from the Nernst–Einstein relation and provide a new avenue to gauge the extent of ion–ion correlations in electrolyte systems.

Theoretical Background of Ionic Conductivity

At low salt concentrations, the ionic conductivity of a 1:1 symmetric electrolyte can be described by the Nernst–Einstein (N–E) equation (eq ), in which the ionic conductivity is only linked to the self-diffusion of ions.[13]where β = 1/kbT is the inverse temperature, q is the formal charge of each ion, and ρ = N/Ω is the number density of the electrolyte (in the formula unit) with Ω as the system volume. σ+s and σ–s are contributions to the ionic conductivity from self-diffusion coefficients D+s and D–s of cations and anions, respectively.

The Nernst–Einstein relation becomes approximated at high salt concentrations, where the effect of ion–ion cross-correlation starts to show up.[11,12,14−19] In this case, the ionic conductivity can be formally defined by the Green–Kubo (G–K) formula[20]where J is the total current density and P is the itinerant polarization in ionic solutions[21] or the Berry phase polarization in solids.[22,23] ⟨···⟩ indicates the ensemble average.

The difference between σG–K and σN–E can be decomposed into contributions from so-called distinct diffusion coefficients of cations D++d, anions D––d, and cation–anion D+–d.[15,18] The name “distinct” emphasizes the nature of cross-correlation between different ions either in the same species or in different species. Subsequently, this allows rewriting the Green–Kubo conductivity aswhere σ++d, σ––d, and σ+–d are the distinct ionic conductivities from the corresponding distinct diffusion coefficients.

Deviations from the Nernst–Einstein relation, i.e., (σN–E – σG–K), can be quantified by either experiments or MD simulations. In experiments, they can be obtained as the difference between the pulsed-field gradient NMR and the impedance spectroscopy measurements for the same system under the same conditions.[5,10,24,25] In MD simulations, one can compute σG–K using either eqs or 4 and σN–E with self-diffusion coefficients obtained from either velocity autocorrelation functions eq or mean squared displacement eq .[20]where t is the time, N is the number of cations or anions in solution, α ∈ {+, −}, v is the velocity of the ith cation or anion, and r is the corresponding position.

It is worth noting that MD results of σN–E obtained by computing self-diffusion coefficients D+s and D–s contain a significant finite-size error because of the hydrodynamic self-interaction in periodic systems.[26,27] To obtain the corrected self-diffusion coefficients Dαs(L → ∞) that are system-size independent, the following formula can be applied.where Dαs(L) is the self-diffusion coefficient obtained from eq with the box length L, ξ is about 2.837297 for cubic simulation boxes,[27] and η is the shear viscosity.

Model Systems and Molecular Dynamics Simulations

Following the spirit of using the ideal potassium chloride (KCl) line for the Walden plot,[3] here we took sodium chloride (NaCl) electrolyte solution as a prototype system for aqueous electrolytes. Water molecules were described by the simple point charge/extended (SPC/E) model,[28] and Na+/Cl– ions were modeled as point charge plus Lennard-Jones potential using the parameters from Joung and Cheatham,[29] which is suitable for highly concentrated solution.[30−32] The stoichiometry of three different simulation boxes (large, medium, and small) is listed in Section A of the Supporting Information. The molecular dynamics simulations were performed with the LAMMPS code.[33] The sizes of cubic simulation boxes were determined by experimental densities.[34] The long-range electrostatics was computed using the particle–particle particle–mesh (PPPM) solver.[35] Short-range cutoffs for the van der Waals and Coulomb interactions in direct space are 9.8 Å. For computing ionic conductivities, NVT (constant number of particles, constant volume, and constant temperature) simulations ran for 20 ns with a timestep of 2 fs and trajectories were collected every 0.5 ps. The Bussi–Donadio–Parrinello thermostat[36] was used to maintain the given temperature of 20 °C.

Because the [BMIM][PF6] system does not show permanent ion pairing,[37] we picked up this model system in our investigation of ILs here. The interaction model and parameters derived from the OPLS-based force field[38] for ILs (OPLS-2009IL)[39] were used for the [BMIM][PF6] system. A charge scaling factor of 0.8e was applied to account for the electronic polarization effects,[40] which was shown to improve the prediction of self-diffusion coefficients.[41] The stoichiometry of three different simulation boxes (large, medium, and small) is listed in Section B of the Supporting Information. Short-range cutoffs are 13 Å for the [BMIM][PF6] system. NPT (constant number of particles, constant pressure, and constant temperature) simulations ran for 100 ns with a timestep of 1 fs, and trajectories were collected every 0.5 ps. The Bussi–Donadio–Parrinello thermostat and the Parrinello–Rahman barostat[42,43] were used to maintain the selected temperatures constant and the pressure at 1.0 atm.

Results and Discussion

System-Size Dependence of the Green–Kubo Conductivity

As shown in Figure , we found that the ionic conductivities computed using the Green–Kubo formula show no system-size dependence. Such a characteristic is similar to that of the viscosity η, which is also a system-size-independent quantity.[27]

Figure 1

Ionic conductivities calculated from the Green–Kubo formula for both (a) NaCl solutions at 20 °C and different concentrations and (b) [BMIM][PF6] ILs at different temperatures.

Despite that there is no obvious reason why this should be the case, note that both the supercell polarization P used for computing the Green–Kubo conductivity and the pressure tensor p used for computing the viscosity are collective properties of the whole system rather than the average of individual particle’s properties.

With periodic boundary conditions, point charge density and point force density are modified by the compensating background as q(δ(r – r) – 1/Ω) and F(δ(r – r) – 1/Ω), respectively. This gives the supercell polarization and the virial part of the pressure as P = ∑q(δ(r – r) – 1/Ω) * r and p = ∑F(δ(r – r) – 1/Ω) ⊗ r. Considering the mathematical similarity between these expressions, it may not be a total surprise that the resulting Green–Kubo conductivity and the viscosity from the linear response theory have the same system-size dependence.

Another angle of looking into this problem may be through the connection between the ionic conductivity σG–K and the Maxwell–Stefan diffusion coefficient D+–M–S. The Maxwell–Stefan diffusion coefficient D+–M–S describes the mutual diffusion between cations and anions, which is independent of the reference frame. In binary systems, it is linked to the Green–Kubo conductivity asRecently, it has been proposed that the system-size dependence of the Maxwell–Stefan diffusion coefficient DM–S in molecular binary mixtures follows the expression[44]Apart from the familiar expression given in eq , the new ingredient is the inclusion of the thermodynamic factor Γ as a correction. When Γ is significantly larger than 1, which happens when the two species like to associate with each other, DM–S becomes effectively system-size independent. What we observed in the case of IL [BMIM][PF6] may fall into this category, where cations and anions attract each other naturally. However, eq simply does not hold for the case of binary electrolyte solution (cations, anions, and solvent molecules). This makes eq not applicable to NaCl solutions. Moreover, the thermodynamic factor calculated from the experimental mean activity coefficient is not much larger than 1 over the whole concentration range of NaCl solutions (see Section C in the Supporting Information), which further indicates that eq may not be suitable to explain Figure .

System-Size Dependence of the Distinct Conductivities

As we found that σG–K is system-size independent (Figure ) and we knew that D+s and D–s (therefore σ+s and σ–s) are system-size dependent (eq ), these together imply that some if not all of the distinct conductivities in eq should also be system-size dependent.

To verify this, we calculated the distinct conductivities for both NaCl solution and [BMIM][PF6] IL with different box sizes, and the results are shown in Figure . It is found that the distinct conductivity σ++d of cations (or σ––d of anions) has a very similar and strong system-size dependence to that of the corresponding σ+s (or σ–s) coming from the self-diffusion of ions. These system-size dependencies are more apparent in the case of NaCl solutions than the case of [BMIM][PF6]. This is likely due to the fact that the viscosities of NaCl solutions are much smaller than those of [BMIM][PF6] ILs, following the relation in eq . In contrast, the cation–anion distinct conductivity σ+–d shows little or no system-size dependence.

Figure 2

System-size dependence of different contributions to Green–Kubo ionic conductivities for NaCl solutions at 20 °C and different concentrations: (a) cation contributions, (b) anion contributions, and (c) cation–anion distinct diffusion contribution; the same for [BMIM][PF6] ILs at different temperatures: (d) cation contributions, (e) anion contributions, and (f) cation–anion distinct diffusion contribution.

Why do these distinct conductivities have different system-size dependencies? A simple argument would be that it is for the sake of symmetry. Since there are five terms in eq , σ++d (or σ––d) should be paired up with σ+s (or σ–s) and this leaves σ+–d on its own. In fact, this is not just an intuition. By connecting Onsager’s phenomenological transport equations with the linear response theory,[45] Schönert showed that the distinct diffusion coefficient and the self-diffusion coefficient have the following general relation for 1:1 electrolytes.where α ∈ {+, −}, β ∈ {+, −}, and Ωαβ are the barycentric-fixed Onsager coefficients; NA is the Avogadro constant; and δαβ is the Kronecker delta function. A barycentric-fixed reference frame means the velocity of the center of mass of the system is set to zero, which is the most suitable reference frame for MD simulations.

For ILs, these Ωαβ coefficients are not independent but follow the expression below because of the conservation of momentum.[45]where M+ and M2 are the molecular weights of cations and anions, respectively.

Putting eqs , 12, and 13 together, one can arrive at the following expressionThis means σ+–d has the same system-size dependence as σG–K in the case of ILs. Since σG–K is system-size independent, therefore, σ+–d is also system-size independent. This theoretical prediction is exactly what is shown in Figure f. Subsequently, eq also indicates that Ω+–, Ω++, and Ω–– are all system-size independent quantities.

Nevertheless, one needs to be aware that there is no such simple relation as eq for a solution made of simple salt and solvent, e.g., NaCl solutions. Therefore, similar behavior of σ+–d shown for NaCl solutions in Figure c remains as simulation observations.

Crossover Box Length between the Nernst–Einstein Conductivity and the Green–Kubo Conductivity

It is known that the self-diffusion coefficients have strong system-size dependence (eq ), and therefore, one would expect that the Nernst–Einstein ionic conductivity has the same tendency. Indeed, it is the case for both NaCl solutions and [BMIM][PF6] ILs as shown in Figure .

Figure 3

Comparison between the system-size dependent Nernst–Einstein ionic conductivity and the Green–Kubo ionic conductivity: (a) NaCl solutions at 20 °C and different concentrations and (b) [BMIM][PF6] ILs at different temperatures.

What is interesting is that for a small enough simulation box, there exists a crossover box length between the Nernst–Einstein ionic conductivity σN–E and the Green–Kubo ionic conductivity σG–K. This is clearly seen in both cases of NaCl solutions and [BMIM][PF6] ILs. The corresponding crossover box length for NaCl solutions is 12.3 Å, and it is 18.0 Å for [BMIM][PF6] ILs.

Of course, one would immediately argue that the actual crossover box length depends on the force field used even for the same type of systems. However, this is not the question that we will dwell on in this work. Instead, the question that matters here is: Can this crossover between σN–E and σG–K be always achieved? Supposing that all cations and anions in the system are paired up permanently, then σG–K will be absolutely zero, while σN–E is not. In other words, if the system has permanent ion pairing, then the crossover between σN–E and σG–K will never happen. Therefore, we restrict our following discussions to the cases where there is no permanent ion pairing.

Implication of the Crossover Box Length for Ion Transport in ILs

The system-size dependence as discussed in previous sections is usually considered as a finite-size error that needs to be corrected. However, here we turn the tables and use it as a tool instead to investigate the role of viscosity in deviation from the Nernst–Einstein relation.

The observation of the system-size independence of σG–K and the crossover box length Lmin implies thatCombining eqs and 9, one can getInserting eq into eq , one arrives at the following expressionSince Λ = σ/c with c as the molar concentration and c = ρ/NA, the above equation in terms of the molar conductivity Λ can be expressed asEquation suggests that deviations from the Nernst–Einstein relation have a linear relation with respect to 1/η. In other words, it states that the system with a high viscosity will have a small deviation from the Nernst–Einstein relation or vice versa. Furthermore, it has not escaped our notice that eq is reminiscent of the well-known Walden rule Λη = k.[46]

Note that Lmin depends on the specific system. For the prototype systems NaCl and [BMIM][PF6] used here, Lmin values are about 12.3 and 18.0 Å, respectively. We notice that Lmin/2 for NaCl solutions is about 6.1 Å, which is close to the Kirkwood correlation length in bulk liquid water.[47]

Then, it is exciting to know what the boundaries set by these prototype systems would look like according to eq when compared to experiments. For this purpose, we took a few of seminal experimental studies on ILs, which promoted the idea of ionicity[5,24,25] and made the following mapping: ΛN–E(L → ∞) ↔ ΛNMR and ΛG–K ↔ ΛImp. This leads to the results shown in Figure .

Figure 4

(a) Experimental deviation from the Nernst–Einstein relation for 13 types of ILs extracted from ref (5), together with the theoretical boundaries (gray area) set by prototype systems NaCl and [BMIM][PF6] using eq . (b) Experimental deviation from the Nernst–Einstein relation for 6 types of ILs at various temperatures using fitting coefficients of VFT equations from refs (24, 25), together with the theoretical boundaries (gray area) set by prototype systems NaCl and [BMIM][PF6] using eq .

Figure a contains 13 different types of ILs measured at 30 °C (see Table 1 in ref (5) and the list of names in Section D of the Supporting Information), and in Figure b, the temperature dependence of molar conductivities and viscosities of 6 types of ILs were measured experimentally and fitted to Vogel–Fulcher–Tammann (VFT) equations.[24,25] It is interesting to see that most experimental data and lines fall into the boundaries set by NaCl and [BMIM][PF6] and follow eq . Since we now know that the value of Lmin depends on the specific system under the investigation, this agreement is somehow fortuitous.

Furthermore, we notice that the corresponding Lmin/2 follows the order [MMIM][TFSI] > [EMIM][TFSI] > [BMIM][TFSI] > [C6MIM][TFSI] > [C8MIM][TFSI], which is in reverse to the alkyl chain length. This suggests that Lmin/2 should be regarded as an effective ion–ion correlation length that goes down as the size of the cation becomes larger, following the attenuation of electrostatic interactions. Alternatively, this trend may originate from the reduction of the dielectric constant of the corresponding ILs with an increase of the alkyl chain. Verifying these implications should be the topics of future studies.

Before closing this section, it is necessary to make a connection to the quantities related to the ionicity in ILs. For example, the deviation in ionicity Δ.[10]According to the Stokes–Einstein relation, ΛN–E can be expressed as follows.where r+ and r– are the hydrodynamic radii for cations and anions, respectively. r̅ is the mean hydrodynamic radius.

Combining eqs , 19, and 20, one arrives at a succinct expression of Δ.As shown in eq , Δ does not explicitly depend on the temperature and the pressure. Therefore, one would expect Δ to be a constant for one specific system. This agrees with the experimental pieces of evidence.[11,48]

When Δ = 1, this implies that ΛG–K = 0 and Lmin is about 3 times of the mean hydrodynamic radius of ions (see the text around eq regarding the constant ξ). Note that in this limit, eq is no longer applicable.[49] This in turn suggests that Δ = 1 limit will never be met when permanent ion pairing is not considered. On the other hand, when Δ = 0, this means that ΛN–E simply equals to ΛG–K, regardless of the system size. For the infinite dilute solution, this means the correlation length will diverge and Lmin → ∞.

Conclusions

The system-size dependence of the Nernst–Einstein conductivity σN–E and the Green–Kubo conductivity σG–K in NaCl solutions and [BMIM][PF6] IL was investigated using MD simulations. It is found that σN–E is strongly system-size dependent as expected, while σG–K does not depend on the system size.

By analyzing the contributions from the distinct diffusion coefficient we further showed that σ+,–d between cations and anions have the same system-size dependence as σG–K, which is exact for the case of ILs and effective for electrolyte solutions.

Due to different system-size dependencies of the Nernst–Einstein conductivity and the Green–Kubo conductivity, there exists a crossover box length where these two types of conductivities become effectively the same. This leads to an expression for the deviation from the Nernst–Einstein relation (ΛN–E – ΛG–K), showing that a low viscosity leads to a strong deviation and a high viscosity leads to a weak deviation (for systems without permanent cation–anion associations), following eq .

This new expression was verified against published experimental data of different types of ILs and the system-specific crossover box length Lmin may provide a new avenue to gauge the ion–ion correlation in the electrolyte system. Future studies should focus on extending the current formulation to the cases that contain permanent ion pairing and investigating the relationship between the hydrodynamic radius of ions, the Lmin, and nano-scale confinement.