Screening length for finite-size ions in concentrated electrolytes
The classical Debye-Hückel (DH) theory clearly accounts for the origin of screening in electrolyte solutions and works rather well for dilute electrolyte solutions. While the Debye screening length decreases with the ion concentration and is independent of ion size, recent surface-force measurements imply that for concentrated solutions, the screening length exhibits an opposite trend; it increases with ion concentration with a dependence on the ionic size. The screening length is usually defined by the response of the electrolyte solution to a test charge, but can equivalently be derived from the charge-charge correlation function. By going beyond DH theory, we predict the effects of ion size on the charge-charge correlation function. A simple modification of the Coulomb interaction kernel to account for the excluded volume of neighboring ions yields a non-monotonic dependence of the screening length (correlation length) on the ionic concentration, as well as damped charge oscillations for high concentrations.
Ionic solutions can be found in a wide range of biological and electrochemical systems VO (); Israelachvily (); David95 (); SafinyaBook (); SamBook () and are studied for both their fundamental properties and industrial applications. A key feature of ionic solutions, captured already within the seminal work of Debye and Hückel (DH) DH1 (); DH2 (), is the screening of electrostatic interactions between charged objects immersed in solution VO (); Israelachvily (); David95 (); SafinyaBook (); SamBook (); DH1 (); DH2 (), characterized by the Debye screening length, .
For a solution of monovalent cations and anions of bulk concentration in a solvent of dielectric constant , the Debye screening length is given by , where is the Bjerrum length (SI units), is the elementary charge and the thermal energy. The Debye length describes not only the screened potential due to external charges, but also the spatial decay of charge-charge correlations that arise from thermal fluctuations in the ionic concentration. Therefore, is also the DH charge-charge correlation length.
The original DH theory was derived for point-like ions and is independent of the ionic diameter. This is reasonable for dilute solutions, but not for more concentrated solutions with a non-negligible volume fraction of ions. Although it is possible to also obtain the DH result for finite-sized ions Attard 1993 (), the charge-charge correlation function according to this theory violates the Stillinger-Lovett second-moment condition SL1 (); SL2 (); Martin88 (), originating from the perfect long-range shielding of the electric field within a conducting medium.
Several theories that improve upon the DH result Attard 1993 (); Carvalho94 (); Caccamo96 (); Lee96 (); Lee97 (); Varela03 (); Janecek09 (); Kjellander2018 () have been suggested in order to surpass the above-mentioned difficulties with the charge-charge correlation function for finite-sized ions. For example, Lee and Fisher Lee97 () introduced a generalized DH (GDH) model based on a Debye charging process of charge oscillations. A common feature emerging from these theories is an oscillatory regime of the charge-charge correlation function at moderately high ionic concentrations. Such oscillations were predicted much earlier by Kirkwood Kirkwood34 (), and the crossover at the onset of oscillations is known as the Kirkwood line. We note that charge oscillations emerge from several Landau-Ginzburg-type theories that are used, for example, to describe the charge density of ionic liquids in the proximity of charged surfaces Bazant11 (); Fedorov14 (); Limmer15 (); Gavish16 (); Gavish18 ().
What rekindled the interest in this problem are recent surface-force apparatus (SFA) experiments on several concentrated electrolytes and ionic liquids Gebbie15 (); Smith16 (). The long-range forces in these experiments revealed an anomalously large screening length at high ionic concentrations. The measured forces follow a common scaling law Smith16 () for all the electrolytes and ionic liquids used in the experiments. The force decays as , where is the screening length. For low ionic concentrations, , where is the ionic diameter, it follows the DH result, . However, in the limit of , it scales differently as . Hence, increases linearly with the ionic concentration, .
The above scaling law, inferred from SFA experiments at high ionic concentrations, has motivated several recent theoretical works Rotenberg18 (); Kinjal18 (); Coupette18 (), but has not yet been fully understood. While these works rely on different assumptions and yield slightly different results, all of them introduce in a similar way some non-Coulomb short-range interactions between the ions (and possibly, the solvent). These interactions become significant at moderately high concentrations, where the average separation between ions is comparable with the ionic diameter.
In the present work we suggest an alternative approach. The ionic finite size is taken into account by a straight-forward modification of the Coulomb interaction. We introduce two possible modifications: either an interaction that is restricted to separations larger than the ionic diameter, or an interaction between finite-sized ions with a distribution of charge that we refer to as the internal charge-density Frydel13 (); Frydel16 (). A general expression for the charge-charge correlation function is derived analytically, by considering Gaussian charge-density fluctuations in the bulk electrolyte.
The outline of our paper is as follows. In Sec. II, we present our model and focus in Sec. II.1 on the modification of the Coulomb interaction for finite-size ions. A general result for the charge-charge correlation function is derived in Sec. III, followed by a description of the different regimes of the correlation length in Sec. III.1. The dilute and concentrated electrolyte limits are discussed, respectively, in Secs. III.2 and III.3. Finally, we discuss our results, relate our findings to experiments and provide some concluding remarks in Sec. IV.
Consider a monovalent electrolyte of bulk concentration . The solvent is modeled as a homogeneous medium with dielectric constant . The system is held at a constant temperature . We focus on properties of the bulk electrolyte, far from any charged objects and surfaces. Furthermore, it is assumed that the electrolyte is far from the liquid-liquid critical point, where the solution phase-separates into two electrolytes of different concentrations Lee96 (). Here, this assumption means that fluctuations in the ionic concentration can be neglected. This is in accordance with the experimental setup described in Ref. Smith16 (), where no inhomogeneity was observed.
The leading electrostatic contribution to the free energy originates from fluctuations in the charge density, , where are the number densities of the cations and anions. Namely, these fluctuations contribute an electrostatic energy term given by
where is the dimensionless electrostatic interaction kernel (units of ). Due to the finite size of the ions, is different than the standard Coulomb kernel, , where is the Bjerrum length. This modification of the Coulomb interaction lies at the heart of our work. Possible forms of are described in detail in Sec. II.1.
Similarly, it is possible to expand the free energy, where is the ions’ mixing entropy, to quadratic (Gaussian) order in such fluctuations, according to
The first term on the right-hand-side of Eq. (II) accounts for the free energy of a homogeneous solution with , while the second term corresponds to the entropic contribution of small fluctuations, and the third term is the electrostatic energy of Eq. (1).
ii.1 Modified electrostatic interaction
The electrostatic interaction between finite-size ions differs from that of point-like ions. This difference is especially important for small ionic separations, relevant for concentrated electrolytes. For simplicity, we assume that both ionic species have a diameter and are spherical-symmetric. Consequently, the interaction between two ions depends only on the ions’ relative position and not on their orientations. We consider two ways in which the ionic interaction can be modified: either a minimal distance cutoff for the interaction or an internal charge density of ions.
As the centers of two hard spheres of diameter must lie least a distance apart, it is possible to restrict the Coulomb interaction to ionic separations of with a cutoff (co), given by , where is the Heaviside function. It is useful to rewrite this function in Fourier space. We denote the Fourier transform of a function as , and find that
The minimal distance cutoff reduces to a cosine term, while the Coulomb interaction is restored in the limit.
An alternative approach is to consider ions with an internal charge-density, . The resulting interaction kernel in Fourier space is
The internal charge-density gives rise in Fourier space to the form factor, . As the ions are monovalent, the charge density satisfies , and the Coulomb interaction is restored in the limit .
We focus on the two simplest internal charge densities, , with spherical symmetry: a homogeneous spherical shell, , and a homogeneous sphere , given by
where is the Dirac delta function. The Fourier transform of these charge densities is given in terms of spherical Bessel functions, , according to
We note that is an even function, because of the assumed radial symmetry of the interaction, and it satisfies , reducing to the standard Coulomb interaction for point-like ions. In the three examples above, the different functions in Fourier space are given by
These three functions are plotted in Fig. 1.
Iii Charge-density correlation function
where is the Debye length. The equipartition over fluctuation modes yields the Fourier transform of the charge-charge correlation function, ,
The correlation function of Eq. (10) is our main result. It reduces to the DH result in the limit of vanishing ionic diameter ( and ). Furthermore, it coincides with the DH result up to quadratic order in with
This expansion to leading order satisfies two important conditions. The vanishing zeroth order () expresses electroneutrality, while the second-order coefficient satisfies the Stillinger-Lovett second-moment condition SL1 (); SL2 (); Martin88 ().
The real-space behavior of the correlation function is obtained from the inverse Fourier transform of Eq. (10). Its spatial dependence is given by exponential terms, determined by the poles of [Eq. (10)] that solve . These poles are generally complex numbers, and at large distances (), the correlation function can be approximated by the contribution of the pole whose imaginary part is closest to the real axis from above. We denote this pole as with , and find that
where is the amplitude and the phase. Both and can be obtained analytically, but we rather focus on , the value of the pole itself.
iii.1 Correlation length
As is evident from Eq. (12), the real part of the pole, , defines a wavenumber of charge-density oscillations, while the imaginary part, , defines the charge-density inverse screening length. A purely imaginary describes a gas-like phase with a finite correlation length, . A general complex corresponds to a liquid-like disordered phase, with finite-ranged charge density oscillations. Finally, a purely real infers a long-range order of alternating positive and negative charges having a wavelength . Numerical results for as function of for the different modified Coulomb interactions of Eq. (II.1) are presented in Figs. 2 and 3.
The inverse decay length and oscillation wavenumber that correspond to the cutoff interaction, are depicted in Fig. 2. For very dilute solutions , the DH result is restored with . For higher concentrations, the correlation is purely decaying, and the inverse decay length, , increases with and slightly differs from the DH result. This behavior persists up to the Kirkwood line, (left dashed line), where the pole has a non-zero real part, , and the finite-ranged charge-density fluctuations become oscillatory.
Beyond the Kirkwood line, the wavenumber increases as the inverse decay length decreases. For high values, the inverse decay length vanishes (right dashed line), and the pole is purely real, corresponding to a long-range order. This qualitative behavior, as is plotted in Fig. 2, restores the results of Lee and Fisher Lee97 (), obtained using the GDH theory. We note that such a long-range order is different than a solid salt crystal and is expected to be unstable at high concentrations. This unphysical high-concentration limit is resolved by considering different modified interactions.
The inverse decay length and oscillation wavenumber that correspond to the modified interaction of an internal charge-density distributed on a spherical shell, , is depicted in Fig. 3. Similarly to the cutoff interaction of Fig. 2, is purely imaginary for low ionic concentrations, and oscillations start to occur at the Kirkwood value, (dashed line), whose value is larger than in the cutoff case. However, unlike the case plotted in Fig. 2, here the disordered phase persists for arbitrarily high ionic concentrations, and no long-range order is formed. Rather, the inverse decay length gradually decays with . Finally, for the modified interaction due to a homogeneous spherical ionic charge-density () of each ion, the results are qualitatively similar to those plotted in Fig. 3 and are not presented explicitly.
It is evident that charge-charge correlations in dilute electrolytes exhibit the same dependence on , regardless of the exact form of the modified interaction. This is sensible, because the typical separation between ions in a dilute solution exceeds the ionic diameter, making finite-size effects weak. This global behavior in the dilute limit is analyzed below in Sec. III.2. In the opposite limit of concentrated electrolyte, on the other hand, Figs. 2 and 3 demonstrate qualitatively different results for the two types of modified interactions. This concentrated limit is further explored below in Sec. III.3.
iii.2 Dilute electrolyte limit
For low ionic concentrations, similarly to the DH result, the inverse decay length is expected to be small. Therefore, the poles of the correlation function in -space are small and can be found from an expansion of in powers of . We write , where and . Such second and fourth order terms are in accordance with the modified interactions of Eq. (II.1).
Substituting the small- expansion of in Eq. (10) yields the following equation for :
At the Kirkwood value, the inverse decay length obtains its maximal value, .
We note that the inverse decay length, obtained from the low-concentration of Eq. (13), agrees with the experimental results of Ref. Smith16 () for low and moderate concentrations. Figure 4 presents data from Ref. Smith16 () for the electrostatic decay length of an aqueous NaCl solution (red circles) and an ionic liquid in propylene carbonate (blue squares). The black curves correspond to the solution of Eq. (13), with the coefficients and treated as fit parameters. For the NaCl solution, the fitted values are , while for the ionic liquid, they are .
Negative values are expected from the Taylor expansion of the functions of Eq. (II.1), reflecting how the Coulombic interaction is diminished due to the finite ionic size. Smaller values correspond to larger values of the Kirkwood line, . The values, on the other hand, correspond to the peak of the inverse decay length, which satisfies . In addition, values were inferred from x-ray scattering experiments Santos11 () for the ionic liquid ( nm) and correspond to the average hydrated ion diameter for NaCl ( nm). We note that the inverse decay length of Eq. (13) vanishes at high concentrations.
iii.3 Concentrated-electrolyte limit
At high ionic concentrations, the oscillation wavenumber, , saturates, while the inverse decay length, , gradually decreases (Figs. 2 and 3). We assume that decays algebraically with the ionic concentration, according to , where and are positive numbers. In order to find the decay exponent, , this expansion is substituted in the denominator of the correlation function [Eq. (10)], and the pole of the correlation function is found by equating all orders of to zero.
The oscillation wavenumber is determined by the highest order term, . The next-order term vanishes for
We distinguish between two possible solutions. For , the inverse decay length is identically zero , as is the case for the cutoff interaction with . Alternatively, can be a degenerate root of , such that and decays gradually. This is the case for the internal charge densities, where . We focus on this latter case, and find the appropriate from the next order term,
resulting in . This result can be written in the form
The above relation is in accordance with our numerical calculations of for and , as is demonstrated in Fig. 5. We note that the above scaling does not coincide with the master curve produced in Ref. Smith16 () from experimentally measured decay lengths, suggesting that .
We propose a simple theory to account for the finite ion-size and its effect on charge-charge correlations. We show that a simple and physical modification of the electrostatic interaction is sufficient in order to capture the non-monotonic dependence of the electrostatic decay length on the ionic concentration, as well as charge-density oscillations at high ionic strength. A hardcore exclusion of the Coulomb interaction was shown to restore qualitatively the GDH result Lee97 () of undamped charge oscillations for high concentrations. Considering the internal charge densities of ions, on the other hand, yields a more realistic result of damped charge oscillations.
While other theoretical frameworks for concentrated electrolytes invoke additional short-range interactions, such as a hardcore repulsion or vdW-like attraction, the origin of the interaction in our present study is purely electrostatic. Some of the Coulomb interaction is excluded due to the finite ion size, resulting in a weaker interaction for small ionic separations . This short-range behavior can affect the charge-density correlation only if some ions are separated by distances smaller than . Therefore, the ions must be considered penetrable. We note that penetrable ions with an internal charge-density were considered in Refs. Frydel13 (); Frydel16 (). In addition, it was shown that the mean spherical approximation (MSA) can be interpreted in terms of the interaction between such spherical charged shells Roth16 ().
For high ionic concentrations, the modified interaction between ions with an internal charge density results in damped charge-density oscillations (Fig. 3). These oscillations decay according to . However, SFA experiments Smith16 () suggest that the electrostatic forces decay with an inverse decay length that decreases more rapidly with the ionic concentration. The relation was shown to be satisfied by several concentrated electrolytes and ionic liquids Smith16 () that fall on a master curve when is plotted as a function of , for all the chemical systems studied.
One of the possible reasons for the different scaling in our work and the experimental results lies in the different meaning of the three physical quantities in the above relation, , , and . First, the ionic diameter within our framework, , characterizes the ionic charge-density in solution, while in Ref. Smith16 (), the ionic diameter is the average bare ionic diameters measured for the two ionic species.
More importantly, the inverse decay length, , in our work is derived from the charge-density correlation length in the bulk electrolyte, while the measured quantity in SFA experiments is inferred from the decay of inter-surface electrostatic forces. Although the two are expected to coincide in the limit of large separations, surface effects may prevent identification of the measured screening length with the predicted bulk correlation length for charge fluctuations. For example, optical-tweezers experiments on charged colloids in nonpolar solvents Bartlett18 () suggest that surface charge-regulation may be the origin of the nonmonotonic behavior of the electrostatic decay length. We note that future scattering experiments may more directly measure the bulk correlation length with no complications due to the surface effects in SFA.
The inverse Debye length also deserves attention, due to the dielectric constant, . In our work, the dielectric constant is independent of the ionic concentration, and the pure solvent value is used. In Ref. Smith16 (), on the other hand, is given in terms of the dielectric constant of the electrolyte solution, . The static dielectric constant of electrolytes decreases with concentration for simple salts in water BarthelBook (); Hasted48 (); Adar18 (), due to excluded solvent volume and electrostatic correlations between the solvent and solute Adar18 (). For concentrated electrolytes, we have found that . The decrement of the dielectric constant, , thus suggests a further decrement of . However, this effect is not strong enough in order to solely explain the experimental scaling of .
We note that our theory cannot be extended to arbitrarily high ionic concentrations within the above concentrated electrolyte limit. At sufficiently high concentrations, the electrolyte approaches a critical point, where it phase-separates into two electrolytes of different concentrations Lee96 (). Close to this phase transition, large concentration fluctuations occur Lee96 (), and the present formulation for a homogeneous electrolyte must be refined. However, we note that such fluctuations were not observed in the discussed experiments.
For higher concentrations, ion pairing can also become significant Bjerrum (); Zwanikken09 (); Adar17 (); Huang18 (). Ion pairing was even suggested as a mechanism for the observed under-screening Huang18 (). However, for simple salts such as NaCl, pairing is not expected to be substantial for concentration of a few molars, as is indicated, for example, by dielectric data Adar18 ().
Finally, we mention possible future extensions of our theory. One can consider other modified interactions and, namely, other internal-charge form factors. For example, it is possible to describe charge densities of solvated ions rather than bare ions or to determine the form factor from a variation of the free energy. Performing such a variation for fixed total charge has resulted in qualitatively similar results. In addition, it is possible to model the solvent explicitly as charge dipoles in order to account for the dielectric decrement in the presence of ions. The compact analytical form of the correlation function [Eq. (10)] should be helpful for developing such extensions in future studies.
Acknowledgments. We thank Susan Perkin, Alpha Lee and Phil Pincus for fruitful discussions and suggestions. SAS is grateful for the support of the Israel Science Foundation and the US-Israel Binational Science Foundation. DA acknowledges partial support from the ISF-NSFC (Israel-China) joint program under Grant No. 885/15.
- (1) E. J. Werwey and J. Th. G. Overbeek, Theory of the Stability of Lyophobic Colloids (Elsevier, New York, 1948).
- (2) J. N. Israelachvili, Intermolecular and Surface Forces, 3rd ed. (Academic, New York, 2011).
- (3) D. Andelman, in Handbook of Physics of Biological Systems, edited by R. Lipowsky and E. Sackmann, Vol. I (Elsevier Science, Amsterdam, 1995), Chap. 12.
- (4) T. Markovich, D. Andelman, R. Podgornik, in Handbook of Lipid Membranes, edited by C. Safinya, J. Rädler (Taylor and Francis), to be published.
- (5) S. Safran, Statistical Thermodynamics of Surfaces, Interfaces, and Membranes (CRC Press, Boca Raton, 2003).
- (6) P. Debye and E. Hückel, Phyzik Z. 24, 185 (1923).
- (7) P. Debye and E. Hückel, Phyzik Z. 25, 97 (1923).
- (8) F. H. Stillinger and R. Lovett, J. Chem. Phys. 48, 3858 (1968).
- (9) F. H. Stillinger and R. Lovett, J. Chem. Phys. 49, 1991 (1968).
- (10) P. A. Martin, Rev. Mod. Phys. 60, 1075 (1988).
- (11) P. Attard, Phys. Rev. E 48, 3604 (1993).
- (12) R. J. F. L. de Carvalho and R. Evans, Mol. Phys. 83, 619 (1994).
- (13) C. Caccamo, Phys. Rep. 274, 1 (1996).
- (14) B. P. Lee and M. E. Fisher, Phys. Rev. Lett. 76, 2906 (1996)
- (15) B. P. Lee and M. E. Fisher, Europhys. Lett. 39, 611 (1997).
- (16) L. M. Varela, M. GarcÃa, and V. Mosquera, Phys. Rep. 382, 1 (2003).
- (17) J. Janeček and R. R. Netz, J. Chem. Phys. 130, 074502 (2009).
- (18) R. Kjellander, J. Chem. Phys. 148, 193701 (2018).
- (19) J. G. Kirkwood, J. Chem. Phys. 2, 767 (1934).
- (20) M. Z. Bazant, B. D. Storey, and A. A. Kornyshev, Phys. Rev. Lett. 106, 046102 (2011).
- (21) M. V. Fedorov and A. A. Kornyshev, Chem. Rev. 114 2978 (2014).
- (22) D. T. Limmer, Phys. Rev. Lett. 115, 256102 (2015).
- (23) N. Gavish and A. Yochelis, J. Phys. Chem. Lett. 7,1121 (2016).
- (24) N. Gavish, D. Elad, and A. Yochelis, J. Phys. Chem. Lett. 9, 36 (2018).
- (25) M. A. Gebbie, H. A. Dobbs, M. Valtiner, and J. N. Israelachvili, Proc. Natl. Acad. Sci. 112, 7432 (2015).
- (26) A. M. Smith, A. A. Lee, and S. Perkin, J. Phys. Chem. Lett. 7, 2157 (2016).
- (27) B. Rotenberg, O. Bernard, and J. P. Hansen, J. Phys.: Condens. Matt. 30, 054005 (2018).
- (28) N. B. Ludwig, K. Dasbiswas, D. V. Talapin, and S. Vaikuntanathan, J. Chem. Phys. 149, 164505 (2018).
- (29) F. Coupette, A. A. Lee, and A. Härtel, Phys. Rev. Lett. 121, 075501 (2018).
- (30) C. S. Santos, N. S. Murthy, G. A. Baker, and E. W. Castner, J. Chem. Phys. 134, 121101 (2011).
- (31) D. Frydel and Y. Levin, J. Chem. Phys. 138, 174901 (2013).
- (32) D. Frydel, J. Chem. Phys. 145, 184703 (2016).
- (33) R. Roth and D. Gillespie, J. Phys.: Condens. Matt. 28, 244006 (2016).
- (34) F. Waggett, M. D. Shafiq, and P. Bartlett, Colloids and Interfaces 2, 51 (2018).
- (35) J. Barthel, H. Krienke, and W. Kunz, Physical chemistry of electrolyte solutions: modern aspects Vol. 5 (Springer Science & Business Media, New York, 1998).
- (36) J. B. Hasted, D. M. Ritson, and C. H. Collie, J. Chem. Phys. 16, 1 (1948).
- (37) R. M. Adar, T. Markovich, A. Levy, H. Orland, and D. Andelman, J. Chem. Phys. 149, 054504 (2018).
- (38) N. Bjerrum, Kgl. Dan. Vidensk. Selsk. Mat. Fys. Medd. 7, 1 (1926).
- (39) J. Zwanikken and R. van Roij, J. Phys.: Condens. Matt. 21, 424102 (2009).
- (40) R. M. Adar, T. Markovich, and D. Andelman, J. Chem. Phys. 146, 194904 (2017).
- (41) J. Huang, J. Phys. Chem. C 122, 3428 (2018).