Activity coefficients of individual ions in LaCl$_{3}$ from the II+IW theory

TL;DR

This paper models the activity coefficients of ions in LaCl₃ by combining ion-ion and ion-water interactions without adjustable parameters, using Monte Carlo simulations and experimental data, achieving accurate predictions.

Contribution

It introduces a parameter-free model that combines Monte Carlo simulations and Born's solvation theory to predict ion activity coefficients in LaCl₃.

Findings

The model accurately predicts ion activity coefficients.

Large opposing ion-ion and ion-water contributions balance each other.

No adjustable parameters are needed, only experimental inputs.

Abstract

We investigate the individual activity coefficients of ions in LaCl$_{3}$ using our theory that is based on the competition of ion-ion (II) and ion-water (IW) interactions (Vincze et al., J. Chem. Phys. 133, 154507, 2010). The II term is computed from Grand Canonical Monte Carlo simulations on the basis of the implicit solvent model of electrolytes using hard sphere ions with Pauling radii. The IW term is computed on the basis of Born's treatment of solvation using experimental hydration free energies. This agreement is remarkable considering the facts that (i) the result is the balance of two terms that are large in absolute value (up to $20 kT$) but opposite in sign, and (ii) that our model does not contain any adjustable parameter. All the parameters used in the model are taken from experiments: concentration dependent dielectric constant, hydration free energies, Pauling radii.

Figures

Tables

Table 1. Table 1: Experimental parameters of ions studied in this work: the valence, z i subscript 𝑧 𝑖 z_{i} , the Pauling radius [ 42 ] , R i subscript 𝑅 𝑖 R_{i} , the hydration Gibbs free energy, Δ ​ G i s Δ superscript subscript 𝐺 𝑖 s \Delta G_{i}^{\mathrm{s}} (taken from Fawcett for Na + , Ca 2+ , and Cl - [ 4 ] , and from Marcus for La 3+ [ 66 ] ), and the Born radius, R i B superscript subscript 𝑅 𝑖 B R_{i}^{\mathrm{B}} (computed from Δ ​ G i s Δ superscript subscript 𝐺 𝑖 s \Delta G_{i}^{\mathrm{s}} on the basis of Eq.  10 with ϵ w = 78.37 subscript italic-ϵ w 78.37 \epsilon_{\mathrm{w}}=78.37 ).

Ion	$z_i$	$R_i/$Å	$\Delta G_i^s$/kJmol-1	$R_i^B/$Å
Na+	1	0.95	-424	1.62
Ca2+	2	0.99	-1608	1.71
La3+	3	1.05	-3145	1.96
Cl-	-1	1.81	-304	2.26

Equations

$${\mu }_i^{\mathrm{EX}}={\mu }_i^{\mathrm{II}}+{\mu }_i^{\mathrm{IW}},$$

$${\mu }_i^{\mathrm{EX}}=kT\ln {\gamma }_i$$

$${\mu }_i={\mu }_i^0+kT\ln c_i+{\mu }_i^{\mathrm{EX}},$$

$${\mathrm{C}}_{{\nu }_{+}}{\mathrm{A}}_{{\nu }_{-}}\rightleftharpoons {\nu }_{+}{\mathrm{C}}^{z_{+}}+{\nu }_{-}{\mathrm{A}}^{z_{-}},$$

$${\gamma }_{\pm }={\gamma }_{+}^{{\nu }_{+}/\nu }{\gamma }_{-}^{{\nu }_{-}/\nu },$$

$${\mu }_{\pm }^{\mathrm{EX}}=\frac{{\nu }_{+}}{\nu }{\mu }_{+}^{\mathrm{EX}}+\frac{{\nu }_{-}}{\nu }{\mu }_{-}^{\mathrm{EX}}.$$

$$u^{\mathrm{PM}}_{ij}(r)=\left\{\begin{array}[]{ll}\infty&\;\mbox{for}\;\;r<R_{i}+R_{j}\\ \dfrac{z_{i}z_{j}e^{2}}{4\pi\epsilon_{0}\epsilon(c)r}&\;\mbox{for}\;\;r\geq R_{i}+R_{j},\end{array}\right.$$

$${\mu }_i^{\mathrm{IW}}(c)=\Delta G_i^{\mathrm{s}}\frac{ϵ(c)-{ϵ}_{\mathrm{w}}}{ϵ(c)({ϵ}_{\mathrm{w}}-1)},$$

$${\mu }_i^{\mathrm{IW}}(c)=\frac{z_i^2{e}^2}{8\pi {ϵ}_0{R}_i^{\mathrm{B}}}\left(\frac{1}{ϵ(c)}-\frac{1}{{ϵ}_{\mathrm{w}}}\right).$$

$$\Delta G_i^{\mathrm{s}}=\frac{z_i^2{e}^2}{8\pi {ϵ}_0{R}_i^{\mathrm{B}}}\left(\frac{1}{{ϵ}_{\mathrm{w}}}-1\right),$$

Full text

This paper is dedicated to honor of the life time achievement of Professor Johann Fischer in the field of statistical mechanics of condensed phases.

1 Introduction

A well-known experimental fact in the world of electrolytes is that the activity coefficient (the excess chemical potential) depends on the electrolyte concentration non-monotonically; it decreases near infinite dilution according to the Debye-Hückel (DH) limiting law [1], goes through a minimum, and increases at high concentrations close to saturation. Many theoretical works tried to explain this behavior from various empirical modifications [2, 3, 4, 5, 6, 7, 8, 9] of the DH theory [1] to advanced statistical mechanical approaches based on the model of charged hard spheres immersed in an implicit-water dielectric solvent [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 5, 6, 20].

In most of these papers before 2010, the dielectric constant was either constant (fixed at its value at infinite dilution, e.g., the dielectric constant of water) or changing with concentration but the change in the interaction with the solvent ignored. A few notable exceptions are the papers of Abbas et al. [5] and Inchekel et al. [20]. In 2010, we proposed the II+IW theory [21, 22, 23, 24, 25] in which we split the excess chemical potential into two components

[TABLE]

where the II and IW terms describe the ion-ion and ion-water interactions, respectively (subscript $i$ refers to an ionic species). This equation implies the approximation that these two terms can be computed independently.

In this theory we took into account the experimental fact that the dielectric constant of the electrolyte changes with concentration due mainly to dielectric saturation (see Fig. 1) [3, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37]. If the properties of the surrounding solvent change, the interaction of the ion with this environment also changes. This is reflected in the IW term. This is the term that was ignored by many authors [12, 13, 19, 16, 17, 18]. We proposed, similarly to other authors [5, 20], that this term can be estimated with the Born-theory of solvation [38].

As far as the II term is concerned, several possibilities are available. It can be calculated with Grand Canonical Monte Carlo (GCMC) simulations [39, 40] as we did in our previous papers and in this work. The mean spherical approximation (MSA) can also be used as we demonstrated [21]. Recently, the concept of the concentration-dependent dielectric constant captured the imagination of Shilov and Lyashchenko [9] who followed Hückel [41] (who advocated this concept as early as 1925) and extended the DH theory for this case. In our comment [25] to the paper of Shilov and Lyashchenko, we have shown that our IW term is equivalent with their solvation term and shown that the extended DH theory gives better results if we use the Born radius in this solvation term. In another recent paper, the concentration dependent dielectric constant was built into the model in the Poisson-Fermi theory of Liu and Eisenberg [8].

The main message of the II+IW theory is that the non-monotonic behavior of the activity coefficient is the result of the competition of the II and IW terms. In earlier works, where the dielectric constant was unchanged, this behavior was reproduced by assigning unrealistically large radii to the hard sphere ions. The non-monotonic behavior in this case was the result of the balance of volume exclusion (hard sphere) and electrostatic terms.

It was said that the increased radius is a “solvated ionic radius” that takes solvation into account by including the hydration shell of tightly connected and oriented water molecules around the ion. We criticized the concept of the “solvated radius” in our previous papers [21, 22, 23, 24, 25] and pointed out that important configurations corresponding to cations and anions in contact are excluded from the statistical sample with this artificial concept. The enlarged solvated radius is primarily useful when we talk about the interactions of ions with water. That is why the Born-radius is used in the calculation of the IW term (this is equivalent with fitting the IW term to experimental hydration free energies). In the II term, however, we use the Pauling radii [42] of the “bare” ions.

Similarly to our previous papers [24, 25], we also present results for activity coefficients of individual ions in this work. We are aware of the turbulence about the measurability (and even the existence) of the single ion activities within the electrochemistry community. There is a recent debate in the literature involving Malatesta [43, 44, 45, 46, 47], Zarubin [48, 49, 50], and the group of Vera and Wilczek-Vera [51, 52, 53, 54, 55, 56, 57, 58, 59, 60]. We do not discuss this issue here; the reader is directed to our previous paper [24]. We just use the single-ion activity coefficients published by the Vera–Wilczek-Vera group [51] and Hurlen [61] for comparison with our theoretical results. When judging accuracy of theories, however, the reader should be aware that experiments for the determination of single-ion activities are uncertain and subjects of unknown errors hidden in the junction potential.

We emphasize that the II+IW theory in its present form does not contain any adjustable parameters. All the results shown in this work have been obtained using experimental hydration energies, experimental concentration dependent dielectric constant, and the Pauling radii. It is quite remarkable that the theory in its bare form can reproduce the non-monotonic behavior. It especially works well for 2:1 electrolytes, where the excess chemical potentials (and especially the II and IW terms) vary in a wide range measured in $kT$. We believe that it is more appropriate to build the competition of opposing forces into the theory through additive free energy terms. The fact that the sum of a large positive (IW) and a large negative (II) term estimates the EX term so well justifies using the concentration dependent dielectric constant as a link between these two terms.

In this paper, we make a step further and report results for a 3:1 electrolyte, LaCl3, for which we managed to find all the necessary experimental data. The range in which the balance of the IW and II terms produces the qualitatively correct behavior of the EX term is even wider, $\sim 40kT$ in the case of La3+.

Studying concentrated electrolytes containing multivalent ions has fundamental importance for many reasons. Concentrated electrolytes tend to be present in confined systems. The selectivity filter of ionic channels, for example, can contain crowded ions attracted by the charged groups of the protein. It is a common view that the electrolyte inside the pore has a reduced dielectric constant compared to the bulk solution [62]. Divalent ions play a special role in biology [63]. They are general messenger particles due to their low concentration and double charge. Energy storage in electrochemical double layer capacitors has a promising prospect when a superconcentrated electrolyte is present in narrow pores [64]. How ions behave under such circumstances is a relevant question especially that the confined systems are in equilibrium with bulk solutions. The correct description of bulk solutions is, therefore, necessary to describe the complete system properly.

2 The II+IW theory

The individual activity coefficient, $\gamma_{i}$, in an electrolyte solution describes the deviation from ideality through the excess chemical potential

[TABLE]

that is defined by

[TABLE]

where $\mu_{i}$ is the chemical potential of species $i$, $c_{i}$ is its concentration, $\mu_{i}^{0}$ is a reference chemical potential independent of the concentration, $k$ is Boltzmann’s constant, and $T$ is the temperature (298.15 K in this work). The reference point is chosen in such a way that $\lim_{c\rightarrow 0}\mu_{i}^{\mathrm{EX}}=0$, where $c$ is the salt concentration [65]. The salt concentration is defined as $c=c_{+}/\nu_{+}=c_{-}/\nu_{-}$ with $\nu_{+}$ and $\nu_{-}$ being the stoichiometric coefficients of the cation and the anion in a simple electrolyte with the stoichiometry

[TABLE]

where C and A refer to cations and anions, while $z_{+}$ and $z_{-}$ are the valences of the ions. The mean activity coefficient is defined as

[TABLE]

where $\nu=\nu_{+}+\nu_{-}$. Accordingly, the mean excess chemical potential is computed as

[TABLE]

The mean quantities, $\gamma_{\pm}$ and $\mu_{\pm}^{\mathrm{EX}}$, can be measured accurately [2, 3, 4].

2.1 Calculation of the II term

The II term is calculated on the basis of the Primitive Model (PM) of electrolytes, where the ions are modeled as charged hard spheres, while the solvent is modeled as a dielectric continuum with a dielectric constant $\epsilon(c)$. The pair-potential describing the II interaction can be given as

[TABLE]

where $R_{i}$ is the radius of ionic species $i$, $\epsilon_{0}$ is the permittivity of vacuum, and $r$ is the distance between the ions. The radius is the Pauling radius in this work, although other choices are also possible [23]. The Pauling radii and other data are collected in Table 1.

The experimental dielectric constant, $\epsilon(c)$, is shown in Fig. 1 as a function of concentration. The equations fitted to the experimental data are found in the caption of the figure. Note that these data are also to be considered with restrictions, because they are obtained from extrapolations to zero frequency using impedance measurement data. For example, the results reported by Barthel et al. [29, 30] from the 1970s for NaCl are quite different from the data of Buchner et al. [37] from 1999, who used lower frequencies (200 MHz vs. 1 GHz). In our previous work, however, we showed that small changes in the $\epsilon(c)$ function has little effect on the activities, because these changes have effects of opposite signs in the II and IW terms that tend to cancel each other. In any case, this effect is smaller than the accuracy of the theory.

We used the the Adaptive-GCMC simulation method of Malasics et al. [39, 40] to determine the II excess chemical potentials of individual ions that correspond to prescribed concentrations. It is an iterative algorithm based on repeated GCMC simulations. The core of the algorithm uses the assumption that the $kT\ln c_{i}$ term changes more than the $\mu_{i}^{\mathrm{EX}}$ term over a little change in the state point. In the case of large ($c>1$ M) concentrations of LaCl3, this assumption does not hold, because of the strong correlations between ions (due to the small dielectric constant and the presence of the trivalent cation that correlates strongly with the anions). For larger concentrations, therefore, we used the simple method of performing several GCMC simulations on a grid and determining the chemical potentials from interpolation.

2.2 Calculation of the IW term

The IW term is computed from

[TABLE]

where $\epsilon_{\mathrm{w}}=78.37$ is the dielectric constant of water. This equation contains only experimental data ($\Delta G_{i}^{\mathrm{s}}$ and $\epsilon(c)$) and is obtained from Born’s treatment of solvation [38]. It is equivalent with equation

[TABLE]

Equation 8 is obtained by writing up the Born-equation for infinite dilution

[TABLE]

and eliminating $R^{\mathrm{B}}_{i}$ (the Born-radius, see Table 1) from Eqs. 9 and 10. The main point of Eq. 8 is to describe the $\epsilon(c)$-dependence of the IW term in a way that it reproduces the experimental hydration free energies in the $c\rightarrow 0$ limit.

The fact that the ions have different radii in the II and IW calculations should not bother the reader. The Pauling radii are actual physical radii of the charged hard spheres present in the particle simulations. The Born radius, on the other hand, is not a physical radius, but rather an effective parameter with which we can make the results of a simple theory match the experimental hydration free energy.

3 Results

It is reasonable to present the results for the 3:1 electrolyte together with data for 2:1 and 1:1 systems. We chose to present results for NaCl and CaCl2 taken from our previous paper [24]. Results for other 1:1 and 2:1 electrolytes can be found in that work. The data for NaCl and CaCl2 represent these systems well as far as the order-of-magnitude comparison between 1:1, 2:1, and 3:1 systems is concerned.

Figure 2 shows the results for the mean activity coefficients. Note that we plot the logarithm that is just the excess chemical potential in $kT$ unit: $\ln\gamma_{\pm}=\mu_{\pm}^{\mathrm{EX}}/kT$. The II+IW theory reproduces the non-monotonic behavior qualitatively for all cases without fitting.

The minimum of the $\ln\gamma_{\pm}(c)$ curve is deeper for electrolytes containing multivalent ions. It is deeper for CaCl2 and even deeper for LaCl3. This behavior can be understood by looking at the activity coefficients of the individual ions, especially the cations. Panels A and B of Fig. 3 show the results for cations and anions, respectively.

The excess chemical potentials of Ca2+ and La3+ vary in a wider range than those of Cl*-* due to their larger valence. The cation term is responsible for the more negative values of the mean excess chemical potentials. The minimum for La3+ is about $-4kT$. Although the excess chemical potential of Cl*-* appears in the mean with a larger weight (see Eq. 6), this term is restricted to a more narrow range, typically between 0 and $-1kT$.

In what case can we consider the quantitative agreement good? Looking at the figures (Figs. 2 and 3), the picture is quite diverse (it gets even more diverse if we look at all the electrolytes published in Ref. [24]). In some cases, the agreement is excellent (Ca2+ and Cl*-* in CaCl2, Na*+* in NaCl). In some cases, there is a systematic deviation (Cl*-* in NaCl and in LaCl3). In some cases, the relative deviation is considerable, but we feel the qualitative agreement appealing (La3+).

Interestingly, the cases, when the absolute value of $\ln\gamma_{i}$ is small (below $1$) seem to be the difficult cases. These are the cases of the monovalent ions (Na*+* and Cl*-). The two cases that are the most interesting for us are those of the multivalent cations (Ca2+* and La3+). In these cases the energies involved are so large that we value the qualitative agreement and we are not bothered by the quantitative disagreement too much.

The quantitative explanation of all these “feelings” and “considerings” is that the non-monotonic behavior is the result of the balance of two competing effects. In our model, these competing effects are the II and IW interactions. If we increase concentration, $\epsilon(c)$ decreases (Fig. 1) and the Coulomb interactions get less screened (Eq. 7). At larger concentrations, the ions are closer to each other on average. Therefore, the II term becomes more and more negative as $c$ increases. The IW term, on the other hand, increases with increasing $c$, because the dielectric environment gets less and less favorable for the ion as $\epsilon(c)$ decreases. At small to moderate concentrations, the II term decreases faster, while at larger concentrations above the minimum point, the IW term catches up and increases faster; hence the non-monotonic behavior.

These two balancing terms are plotted in panels A and B of Fig. 4 for La3+ and Cl*-, respectively, in LaCl3*. This figure explains why the qualitative agreement for La3+ is so valuable (the corresponding figure for Ca2+ can be found in our previous papers [24, 25]). The II term goes down to $-20kT$ and the IW term goes up to $20kT$ at concentrations close to saturation. Their sum (EX) is so close to the experimental value (within $0.5kT$) that we consider the agreement remarkable especially in the light of the fact that our model does not contain any adjustable parameters. Looking at the results for many 1:1 and 2:1 electrolytes [24], deviations are usually in the range of a few tenth of $kT$s. That seems to be a natural limit to the accuracy of our very simplified model.

The situation with Cl*-* is the opposite. Here, the II and IW terms are within $3kT$ and the sum is smaller than $1kT$, so there is a more delicate balance in this case. Any error in the II and IW terms separately has a larger impact here, but its effect on the mean is less, because the cation dominates the mean.

For LaCl3, results are shown only for concentrations up to 1.2 M. At larger concentrations, close to 1.2 M, the interactions between ions become stronger due partly to the reduced dielectric constant, partly to the shorter average distance between the ions. Because of the strong interactions, successful insertions and deletions of ions (especially, La3+) are very rare, and sampling in GCMC simulations becomes inadequate above 1.2 M. Recent studies using molecular dynamics simulations, spectral graph theory, and various experimental techniques imply the presence of ion aggregates from ion pairs to networks in NaCl and KSCN [71, 72, 73, 74, 75, 76]. There is no reason to be surprised that strongly correlated ions in ordered structures appear as we approach the solubility limit [77, 78]. The ions of the electrolyte, however, can be considered solvated even if these structures are present. The radial distribution functions of Fig. 5A for $c=1.2$ M show a fluid-like behavior at long distance, while large first peaks indicate the increased probability of ions forming aggregates. Their nature is revealed by the snapshot shown in Fig. 5B. These aggregates, however, are not solid structures, but momentary groups of ions that are forming and breaking up during thermal motion. The formation of these is closely related to the disruption of the hydrogen-bond network of water, so the simulation study of the phenomena is probably beyond the scope of the implicit solvent framework used here. The solubility of electrolytes, a closely related problem, has been thoroughly investigated by several research groups [79, 80, 81, 82, 83, 84] using explicit-water models.

Summarized, we feel that the relatively good agreement of the individual excess chemical potential with experiments for multivalent cations (results for Mg2+ and Ba2+ are found in our previous work [24]) is a strong evidence that supports our theory. In these cases, the absolute values of the competing II and IW terms are so large that the error in the sum is dwarfed by them.

Acknowledgements

We gratefully acknowledge the financial support of the Hungarian National Research Fund (OTKA NN113527) in the framework of ERA Chemistry.