A multiple decay-length extension of the Debye–Hückel theory: to achieve high accuracy also for concentrated solutions and explain under-screening in dilute symmetric electrolytes

Abstract

The Poisson–Boltzmann and Debye–Hückel approximations for the pair distributions and mean electrostatic potential in electrolytes predict that these entities have one single decay mode with a decay length equal to the Debye length 1/κ_D, that is, they have a characteristic contribution that decays with distance r like e^−κ_Dr/r. However, in reality, electrolytes have several decay modes e^−κr/r, e^−κ′r/r etc. with different decay lengths, 1/κ, 1/κ′ etc., that in general are different from the Debye length. As an illustration of the significance of multiple decay modes in electrolytes, the present work uses a very simple extension of the Debye–Hückel approximation with two decay lengths, which predicts oscillatory modes when appropriate. This approach gives very accurate results for radial distribution functions and thermodynamic properties of aqueous solutions of monovalent electrolytes for all concentrations investigated, including high ones. It is designed to satisfy necessary statistical mechanical conditions for the distributions. The effective dielectric permittivity of the electrolyte plays an important role in the theory and each mode has its own value of this entity. Electrolytes with high electrostatic coupling, like those with multivalent ions and/or with solvent of low dielectric constant, have decay lengths in dilute solutions that substantially deviate from the Debye length. It is shown that this is caused by nonlinear ion–ion correlation effects and the origin of under-screening, i.e., 1/κ > 1/κ_D, in dilute symmetric electrolytes is analyzed. The under-screening is accompanied by an increase in the effective dielectric permittivity that is also caused by these correlations. The theoretical results for the decay length are successfully compared with recent experimental data for simple electrolytes in various solvents. The paper includes background material on electrolyte theory and screening in order to be accessible for nonexperts in the field.

---

1 Introduction

1.1 Brief overview of electrolyte theories and screening

Electrolytes are ubiquitous and in many fields of science and industrial applications it is important to understand their properties and to be able to calculate, for example, various thermodynamical and structural quantities. The study of electrolyte systems has a long history and is still a very active field of science. This illustrates both the importance and the complexity of such systems.

There exist an abundance of statistical mechanical theories for electrolytes that are based on approximations of various degrees of sophistication and accuracy. A selection of the most prominent ones that are relevant for bulk systems are presented and discussed in ref. 1–6. A cornerstone is the Debye–Hückel (DH) approximation,⁷ which gave an answer to an important puzzle at the time of its development in the 1920s, namely experimentally observed deviations from ideality in very dilute electrolyte solutions. In solutions of nonelectrolytes such deviations disappear proportionally to the concentration, while for electrolytes they decay like the square root of concentration due to the long range 1/r decay of the Coulomb interactions as function of distance r. In the DH approximation, the mean electrostatic potential from a charge is screened and decays like a Yukawa function e^−κ_Dr/r, where κ_D is the is the Debye parameter and 1/κ_D the Debye length, which are obtained from

where β = 1/k_BT, k_B is Boltzmann's constant, T is the absolute temperature, ε_r is the dielectric constant of the solvent, ε₀ is the permittivity of vacuum, n_j is the number density and q_j the charge of ions of species j and the sum is taken over all species of ions in the system. The exponential screening leads to the Debye–Hückel limiting laws for the mean activity coefficient γ_± and the osmotic coefficient ϕ, which are for a binary electrolytein the limit of zero ionic density. These limiting laws and the existence of exponential screening for large r are, in fact, exact features of statistical mechanics of dilute electrolytes.

The DH approximation is the linearized version of the Poisson–Boltzmann (PB) approximation, which is a mean field theory where the correlations between ions in the ion cloud around each ion are neglected. It gives a quite good description of electrolytes for low electrolyte concentrations, but fails when the concentration is increased. We here primarily consider the primitive model of electrolyte solutions, where the ions are charged hard spheres and the solvent is modeled as a dielectric continuum with dielectric constant ε_r, and we will also treat classical plasmas of such ions for which ε_r = 1.

A large number of improvements of the DH approximation has been suggested over the years. One of the most well-known is the approach by Pitzer,^8–10 where various DH expressions for thermodynamical quantities are parametrized to fit experimental data for a wide concentration interval. In most of the improved DH approximations, the decay length of the mean electric potential ψ_i(r) from an ion of species i is equal to the Debye length. The actual decay length in electrolytes is, however, in general not equal to the Debye length; it is only approximately equal to the Debye length for sufficiently low concentrations. For higher concentrations the decay length depends significantly on other system parameters than those that appear in eqn (1), for example the ion sizes.

When the electrolyte concentration is increased and the decay length starts to deviate noticeably from the Debye length, the mean electrostatic potential ψ_i(r) decays for large r like e^−κr/r where κ ≠ κ_D. The radial charge density ρ_i(r) around an ion and, in general, the pair distribution function g_ij(r) also decay in this manner. The deviation of κ from κ_D is not described by the DH or PB approximations and neither is the fact that these functions become oscillatory at high ionic densities.^11,12 The latter feature was found by Kirkwood already in the 1930s^13,14 and later confirmed in other approaches by Martynov¹⁵ and Outhwaite.¹⁶ The presence of such oscillations is connected to the fact that there exist several decay modes in electrolytes, that is, contributions to the distribution functions and the electrostatic potential that decay like Yukawa functions e^−κr/r, e^−κ′r/r etc. where the decay parameters κ, κ′ etc. are different. For low ionic densities, the two leading decay parameters, which give the longest decay lengths, satisfy κ < κ′. When the ion density is further increased, κ and κ′ approach each other and at the so-called Kirkwood crossover point they become equal. Beyond that point κ and κ′ turn complex and we have κ = κ_ℜ + iκ_ℑ and , where i is the imaginary unit and κ_ℜ and κ_ℑ are real (underbar denotes complex conjugation). The pair of decay modes then gives rise to a contribution that decays like an “oscillatory Yukawa function” e^−κ_ℜrcos(κ_ℑr + α), where α gives the phase for the oscillations. This kind of oscillatory mode occurs in concentrated solutions and is the dominant one for molten salts.

The presence of several decay modes with different decay parameters has been numerically verified for simple electrolytes by computer simulations,^17–20 the Hypernetted Chain (HNC) approximation^18,21–23 and the Generalized Mean Spherical Approximation (GMSA).²⁴ The leading modes give rise to the oscillatory behavior at high ion densities in accordance with the scenery above. Many other approximate theories for electrolytes also predict such decay modes and the occurrence of the Kirkwood cross-over, for example the Mean Spherical Approximation (MSA),^25–29 the Linearized Modified Poisson–Boltzmann (LMPB) approximation by Outhwaite,^16,30,31 the Modified Debye–Hückel (MDH) approximation by Kjellander,³² the closely related “Local Thermodynamics” (LT) approximation by Hall,³³ the Generalized Debye–Hückel (GDH) approximation by Lee and Fisher,^34,35 the Modified MSA by Varela and coworkers,^36–39 the charge renormalization theory by Ding et al.⁴⁰ and the ionic cluster model approach by Avni and coworkers.⁴¹

These theories, from GMSA onwards, are linear approximations, meaning that ψ_i(r) and ρ_i(r) are proportional to the ionic charge q_i. They provide explicit equations for the decay parameter with multiple solutions κ, κ′ etc. including complex-valued ones when appropriate and so does also the original theory by Kirkwood and that by Martynov. All give decay parameters that are dependent on the ionic diameter a and qualitatively they give similar results for κ/κ_D and κ′/κ_D as functions of κ_Da. Most of these approximations have been used to obtain thermodynamical quantities, distribution functions and/or the mean electrostatic potential.^{15,16,27,28,31–33,37,39,42–47} Xiao and Song^48–50 have developed a linear theory that exploits all decay modes obtained in many of these approximations to calculate thermodynamical quantities for electrolytes. We will return to some of these linear, approximate theories later. All of them work in general well for low ionic densities but have shortcomings, especially for high densities.

When the electrostatic coupling between the ions is not weak, linear theories are inadequate. Many nonlinear approximations for electrolytes have been used to obtain pair distributions and thermodynamical quantities. They include the HNC approximation^51–55 and other integral equation theories.^54,56–61 There are also improved versions of the PB approximation that correct for the neglect of correlations in the ion cloud around each ion, like the Modified PB approximation by Outhwaite et al.^16,62–64 and the correlation-enhanced PB theory by Su and coworkers.⁶⁵ Furthermore, various field theories^4,66–70 have been developed for the study of bulk electrolytes. Pair distributions and thermodynamical quantities for bulk electrolytes have, of course, been calculated also by simulations.^{11,12,17–20,55,59,61,71–80}

The nonlinear theories do not in general provide any explicit equation for the decay parameter like the linear ones above do, but from these theories and from computer simulations, the values of κ, κ′ etc. have been extracted numerically for various electrolytes.^{17–23,81–84} The decay parameter κ of the leading decay mode can in many cases be obtained by curve fitting to the tails of the functions for large r, but it is in general not straightforward to obtain those of the other modes and one must use systematic means to extract them from the distribution functions.¹⁸

One approximate way to incorporate some effects of nonlinearity into linear theories is to introduce ion pairing as originally done by Bjerrum.⁸⁵ In this approach one assumes that anion–cation pairs exist in equilibrium with free ions in the electrolyte, so the concentration of the latter ions is lower than the total concentration. Since the pairs are electroneutral, they do not contribute to the Debye parameter in eqn (1) and the Debye length 1/κ_D becomes larger than in their absence – a quite trivial reason for a change in decay length. The pair formation is a nonlinear phenomenon because the electrostatic interaction between two paired ions is large. The effects of this nonlinearity is included in the value of the equilibrium constant that describes the equilibrium.

There exist various ways to construct theories of pairing. One kind of approach is to include ion pairs in a DH-type theory^86–88 for the electrolyte and another is to base the ion pair theory on the MSA;^88–92 see ref. 93 for a comparison of different variants of these approaches. Fisher and coworkers^34,88,93 stress the importance of including the interactions between the dipoles formed by the ion pairs and the free ions. In many cases, it is reasonable to consider a part of the effect from the presence of such dipoles as an augmentation of the dielectric permittivity of the electrolyte. We will treat some aspects of ion pairing and changes in dielectric permittivity later. Transient ion pairing is, in fact, an aspect of strong anion–cation correlations, which affect the decay length and the (effective) dielectric permittivity of the electrolyte.^94,95

Electrolyte systems have in general a complex behavior and most theories mentioned above are quite complicated to use, in particular those that yield accurate distribution functions. A fair question to ask is how thermodynamic properties and distribution functions possibly can be accurately obtained without undue complications and without the use of empirical fitting parameters. In order a get a perspective on this question, let us first consider some approaches that have been used for such a task.

A simple linear approximation that takes into consideration that the decay length depends on the ion diameter is the Corrected Debye–Hückel (CDH) approximation by Abbas, Nordholm and coworkers.^78,96 They make an Ansatz that the radial charge density ρ_i(r) around an ion (and therefore the electrostatic potential ψ_i(r)) is proportional to e^−κr/r, where the decay parameter κ is not the same as κ_D. The value of κ is then determined from an optimization of an approximate free energy of the system. By combining CDH with the Carnahan–Starling equation of state⁹⁷ for hard spheres, they obtain the activity coefficient for monovalent (1 : 1) electrolytes in aqueous solution in very good agreement with simulations for a wide range of concentrations, including concentrated solutions.⁷⁸ For 2 : 1 electrolytes they obtain equally good agreement except for small ion sizes, but the latter results are still good. As pointed out by Abbas et al.,⁷⁸ the thermodynamical quantities obtained from this Ansatz are in very good agreement with simulation results for high densities despite that the actual ρ_i(r) is oscillatory rather than proportional to e^−κr/r.

Likewise, for a parametrized DH approach with fitting parameters for thermodynamics like that of Pitzer mentioned earlier, the structural entities are qualitatively different from the accurate ones at high concentrations. The same is the case for approaches that use the PB approximation.

Attard⁸⁴ also makes the assumption that ρ_i(r) is proportional to e^−κr/r. By requiring that ρ_i(r) fulfills some necessary conditions that we will consider later, he obtains the Self-Consistent Screening Length (SCSL) approximation for the single decay parameter κ, which can be a real or complex number depending on the ion density. This is a linear theory that has been successfully used to calculate thermodynamical quantities for monovalent electrolytes solutions with low up to moderate concentrations. Another theory that uses only one decay mode in the calculation of thermodynamical properties is the MDH approximation.³² It is also successful for low up to moderate concentrations. These two approximations will be treated in quite some detail later.

Janecek and Netz⁸⁰ investigated whether properties of electrolytes can be described in terms of an effective screening length. They extracted a single effective screening parameter κ_eff by simulation of various electrolyte solutions (corresponding to 1 : 1 and 2 : 2 salt in water) and used it to calculate thermodynamic properties from a DH-like theory where κ_D is replaced by κ_eff. For high concentrations where the simulated radial charge density is oscillatory, a complex-valued screening parameter was extracted instead of κ_eff. This approach works very well for the thermodynamics of 1 : 1 electrolytes for the concentration range investigated and somewhat less well for the 2 : 2 case. The parameter κ_eff is in general not the same as the decay parameter κ of the asymptotic decay of ρ_i(r) discussed in the current work,† but κ_eff approaches the actual κ when the concentration is decreased towards zero. As pointed out by the authors, the effective screening parameter “does not provide an accurate description of the charge density distribution.” They also used the ion-pairing DH-type theory by Fisher and coworkers^88,93 to calculate activity coefficients, which worked very well for 2 : 2 electrolytes and less well for monovalent ones.

From these examples we see that even if an approximation reproduces thermodynamical data accurately, the structural entities like ρ_i(r) and the pair distribution function g_ij(r) obtained in the approximation do not need to be equal to the correct ones. Even when these functions are oscillatory, the thermodynamical entities can be accurately obtained from a single Yukawa decay mode with a real decay parameter. Likewise, when there are more than one decay mode, the thermodynamics can be accurately obtained from just one Yukawa mode. The value of the decay parameter used is in general different from that of the actual κ of the electrolyte. To judge whether an approximation gives a fully accurate description of a system, it is accordingly important to investigate both thermodynamical and structural quantities. A successful theory gives accurate results for both. An important aim of the current work is to develop a simple theory that accomplishes this. Thereby, the presence of more than one decay mode will be included in the treatment of simple electrolytes.

Incidentally, it should be noted that several simultaneous modes that decay like monotonic or oscillatory Yukawa functions have also a prominent role in dense ionic liquids, like room temperature ionic liquids. An analysis of these decay modes is important for the understanding of such systems – not only for structure and thermodynamics, but also for interactions in the systems including surface forces. The modes can, for example, be detected in surface force experiments. It has recently been shown^95,98 that all such decay modes, including those that are dominated by “packing” of molecules in the dense liquids, are in general also decay modes of the electrostatic interactions and are therefore governed by the dielectric response of the liquid.‡ The same is true also for the decay mode of long-range density–density correlation fluctuations with a decay length that diverges on approach to a critical point and the very long-range monotonically decaying mode that recently has been found⁹⁹ for dense electrolytes. Likewise, for dilute electrolyte solutions with discrete molecular solvent, oscillatory decay modes that are dominated by the structure of the solvent also constitute decay modes of the electrostatic interactions.^95,98 Such modes have often been designated as “solvation forces;” in the case of aqueous systems “hydration forces.”

1.2 Overview of the current work

We will limit ourselves in the current work to solutions of simple electrolytes in the primitive model and to classical plasmas. In order to illustrate the importance of the various decay modes, we will include two Yukawa modes and show that both modes are crucial in order to obtain very good description of both structural and thermodynamical quantities for a wide range of ion densities, including concentrated solutions and dense plasmas. As concluded in the previous section, a successful theory accurately predicts both kinds of quantities. The decay modes are physically distinguishable and the decay lengths 1/κ and 1/κ′ constitute characteristic properties of the propagation of screened electrostatic interactions in the electrolyte. For the oscillatory case, where κ and κ′ are complex, the decay length 1/κ_ℜ and the wave length 2π/κ_ℑ have the corresponding roles. The basic idea in this approach is simple: for a radial function f(r) like ψ_i(r), ρ_i(r) or the potential of mean force w_ij(r) between two ions, one makes the approximate Ansatz that f(r) is given byfor all r outside the hard core of the ion. The coefficients C and C′ and the decay parameters κ and κ′ are then determined from necessary statistical mechanical conditions for the radial function in question. Moreover, the coefficients are expressed in terms of physical entities like effective charges of ions and effective permittivities of the electrolyte – entities that are precisely defined in exact statistical mechanics as explained later. These entities are important for the understanding of screening phenomena and intermolecular interactions in electrolytes and the current paper also serves the purpose of illustrating their meaning and use. The simple structure of the approximations used in this work facilitates this objective.

An important further ingredient is the use of an exact relationship§ between the decay parameter κ and the effective charges, namely^100,101

where q^eff_j is the effective charge of a j-ion. This relationship, which is valid for a system with spherical ions, is similar to the definition of κ_D in eqn (1), but it has the product q_jq^eff_j instead of q_j² in the right hand side (rhs). In the next section we will see in a simple manner that eqn (4) is a direct consequence of the ion–ion correlations in the ion cloud around each ion. When such correlations are neglected like in the DH and PB approximations, we have q^eff_j = q_j for the ions in the cloud and eqn (1) for the Debye parameter is obtained.

In the present work, we will primarily investigate systems with sufficiently low electrostatic coupling so that electrostatic response can be linearized. The resulting mathematical formulas for both the structural and thermodynamic quantities are simple and straight-forward to apply in practice. They are obtained by using first-principle statistical mechanics and are free from any fitting parameters. The results obtained for conditions corresponding to monovalent ions in aqueous solutions at room temperature are in very good agreement with computer simulations and other accurate calculations for all concentrations. Such conditions also include classical plasmas at high temperatures. We will also see why it is insufficient with one single decay mode for the accurate description of distribution functions of electrolytes with medium to high densities.

Nonlinear correlational effects for high electrostatic coupling will also be investigated. In fact, the long-range nature of the Coulomb interactions gives such effects in highly diluted electrolyte solutions and thin plasmas. They can cause substantial deviations of κ from κ_D, which are given by the following exact limiting law¹⁰¹ – here written for the case of a binary electrolytes¶

where Λ = [small script l]_Bκ_D is a coupling constant, [small script l]_B = βq_e²/(4πε_rε₀) is the Bjerrum length, q_e is the elementary (protonic) charge and z_j = q_j/q_e is the ionic valency. The deviation of κ²/κ_D² from the value 1, as described by this law, arises solely from the tails of the long-range, purely electrostatic correlations between the ions when r → ∞ and does not depend on any other properties of the ions than their charges. These other properties, like the ion size, enter in the following term (not shown in eqn (5)) that is proportional to Λ², i.e., proportional to the ionic concentration.

For symmetric electrolytes (z₊ = |z₋|), the second term in the rhs of eqn (5) is identically equal to zero and we have κ < κ_D for sufficiently dilute solutions because ln Λ < 0 there, which means that the decay length is larger than the Debye length (“under-screening”). For dilute solutions of asymmetric electrolytes (z₊ ≠ |z₋|), the second term is dominant and for such electrolytes κ > κ_D in the limit of infinite dilution, so the decay length is smaller than the Debye length (“over-screening”).

These predictions have been verified experimentally. For a highly asymmetric electrolyte, z₊ = 12 and z₋ = − 1, the positive deviation of κ from κ_D has been experimentally measured in dilute solutions by Kékicheff and Ninham¹⁰³ in agreement with the leading term in eqn (5). For symmetric electrolytes, the negative deviation of κ from κ_D, has recently been measured in dilute solutions of simple electrolytes with high electrostatic coupling by Smith, Maroni, Trefalt and Borkovec¹⁰⁴ and Stojimirović, Galli and Trefalt¹⁰⁵ (henceforth referred to as “Trefalt and coworkers”) by surface force experiments. This latter case will be discussed in quite some detail in the current work.

While all linear approximations mentioned earlier give the correct limiting laws for the thermodynamical properties, eqn (2), in the limit of zero ionic density, they do not give the correct behavior of κ²/κ_D² in this limit as given by eqn (5). For example, the MSA, GMSA, LMPB, MDH, LT and GDH approximations predict that κ²/κ_D² decays proportionally to Λ² when the density goes to zero. On the other hand, the HNC approximation, which is a nonlinear approximation, gives a κ that agrees with eqn (5) for both symmetric²² and asymmetric⁸¹ electrolytes at low ionic densities. In the analysis presented in the current work, we will see the importance of including nonlinear ion–ion correlation effects in order to obtain the correct behavior of the decay parameter for dilute systems. Such effects influence the values the effective ionic charges and hence κ viaeqn (4).

The paper is organized as follows. In Section 2, some relevant parts of the theory of electrolytes are presented as a background for the rest of the paper and the main approximations made in the current work are introduced. A couple of approximations with a single decay-mode are treated in Section 3 and a simple equation for κ that is used in a large part of the paper is obtained. The Multiple-Decay Extended Debye–Hückel (MDE-DH) approximation, developed in this work, is described in Section 4, starting with the simplest version, which is sufficient provided that the ionic sizes are not too large so ionic core–core correlations are not too prominent. This version of the theory gives both radial distribution functions and thermodynamical quantities in very good agreement with simulations and Hypernetted Chain (HNC) calculations that are very accurate for the systems studied. The complete version of the MDE-DH approximation is used to successfully treat dense systems with large ions, whereby quite intricate details of the pair distribution functions are obtained in agreement with simulations. In Section 5, effects of nonlinear ion–ion correlations on the main decay mode are investigated. The screening properties of thin plasmas and dilute solutions of symmetric electrolytes are thereby studied, in particular the under-screening that occurs in such electrolytes with high electrostatic coupling. This under-screening is accompanied by an augmentation of the effective dielectric permittivity of the electrolyte caused by the strongly nonlinear ion–ion correlations. The various types of deviations of κ from κ_D in symmetric electrolytes and the reasons for them are studied in Section 6, whereby an approximate limiting law for κ²/κ_D² is derived for symmetric electrolytes that extends the exact law (5) by adding higher order contributions. An analysis is made of the experimental results by Trefalt and coworkers^104,105 for the deviations of κ from κ_D in dilute solutions of simple electrolytes with high electrostatic coupling. Different mechanisms for this phenomenon are discussed. Finally, in Section 7 the main results of the paper are summarized and concluding remarks are given. As an aid to the reader, this section contains references to some central equations of this work.

2 Background and basic approximations

We will consider electrolytes in bulk phase using the primitive model. This model is also applicable to classical plasmas where ε_r = 1. Each ion has a point charge q_j at its center and the total ion density in bulk is .

The average density of ions of species j at distance r from the center of an ion of species i is equal to n_jg_ij(r) = n_je^−βw_ij(r), where g_ij(r) is the pair distribution function and w_ij(r) is the pair potential of mean force. The mean charge density around the i-ion is given by the radial charge-distribution function

where h_ij(r) = g_ij(r) − 1 and the last equality follows from the electroneutrality condition for the bulk phase. When dealing with ρ_i(r) and other functions of r, we will denote this i-ion as the “central ion” and the surrounding charge density constitutes its ion cloud.

In the Poisson–Boltzmann approximation one assumes that

w_ij(r) = q_jψ_i(r) (PB) (7)

for r ≥ a, where ψ_i(r) is the mean electrostatic potential from the ionic charge q_i and the surrounding charge density ρ_i(r). The notation “(PB)” on the equation means that the equation is applicable in the PB approximation only. The potential ψ_i(r) decays like
where C_i is a constant. In the limit of infinite dilution κ_D → 0, C_i → q_i/(4πε_rε₀) and ψ_i(r) approaches the ordinary Coulombic potential q_i/(4πε_rε₀r) from the charge q_i. Therefore, for an electrolyte of any ionic density, it is reasonable to define an effective charge q^eff_i of the ion from C_i = q^eff_i/(4πε_rε₀) in the PB approximation, so we havewhere q^eff_i ≠ q_i. One can say that q^eff_i is the charge that the ion is perceived to have when seen from a distance, that is, as judged from the magnitude of the large r tail of the electrostatic potential. The magnitude of the tail is affected by the amount of charge in the ion cloud close to the central ion and the value of q^eff_i is determined by the distribution of charge there. Since ψ_i(r) varies in a nonlinear manner with q_i, the effective charge also varies nonlinearly with q_i. It follows from eqn (7) that in the PB approximation we havewhich does not obey the required symmetry w_ij = w_ji, a well-known deficiency of this approximation. The reason for this deficiency is that the ions in the ion cloud are treated in a different manner than the central ion. The former are treated as point ions that do not correlate with each other. They only correlate with the central ion, which is the only ion in the system that has an ion cloud of its own and therefore an effective charge q^eff_i ≠ q_i. All ions apart from the central ion enter in w_ij with their actual (bare) charge q_j. This applies also to those of the same species as the central ion.

In the linearized version of the PB approximation, the Debye–Hückel approximation, it is assumed that βw_ij(r) is sufficiently small so that e^−βw_ij(r) ≈ 1 − βw_ij(r) and the following well-known expressions are valid for r ≥ a

By making the identification of q^eff_i given by eqn (8) we obtainso we can write for r ≥ aNote that in this case, q^eff_i is proportional to q_i and so are ψ_i and ρ_i. In the DH approximation we have g_ij(r) = 1 − βw_ij(r) with w_ij(r) equal to the rhs of eqn (9) for all r ≥ a.

The total charge density associated with an i-ion, including the charge of the ion itself, is ρ^tot_i(r) = q_iδ⁽³⁾(r) + ρ_i(r), where δ⁽³⁾(r) is the three-dimensional Dirac delta function. The charge density must satisfy the condition of local electroneutrality, which is satisfied for ρ_i in eqn (11) and (14). The electrostatic potential is given by the solution to Poisson's equation −ε_rε₀∇²ψ_i(r) = ρ^tot_i(r) and can be written as

where ϕ_Coul = 1/(4πε_rε₀r) is the (unit) Coulomb potential and the integral is taken over the whole space (this is the convention throughout the paper for integrals without limits).

The (unit) screened Coulomb potential in electrolytes can be defined from the potential ψ_q(r) from a point charge q in the limit of zero q in the following manner||

In the PB and DH approximations we have
This function governs the spatial propagation of electrostatic potentials in the electrolyte. In the DH approximation we have ψ_i(r) = q^eff_iϕ_Coul*(r) for r ≥ a and in the PB case ψ_i(r) ∼ q^eff_iϕ_Coul*(r) for large r.

Also in the general (exact) case, ϕ_Coul*(r), w_ij(r) and the other functions decay like a Yukawa function if the density is not too high, but the decay parameter κ deviates from κ_D and the prefactors have other values. In exact statistical mechanics, all ions are treated on the same basis, which means that w_ij is symmetrical with respect to i and j and all ions have an effective charge q^eff_l ≠ q_l for l = i, j. Dressed Ion Theory (DIT),^{22,23,100,101} which is an exact reformulation of classical, equilibrium statistical mechanics of electrolytes in the fluid state, gives a formalism to handle these matters. Provided that the ion density is not too high we have from DIT

where is an effective dielectric permittivity of the electrolyte that differs from the dielectric constant ε_r of the solvent. Thus specifies the magnitude of the tail of ϕ_Coul*(r) for large r and gives a measure of the dielectric response of the electrolyte. In the PB and DH approximations one sets equal to ε_r, so the dielectric constant of the pure solvent is used in ϕ_Coul* of these approximations instead , which is of a property of the electrolyte.

The mean electrostatic potential due to an i-ion decays like

which is in analogy with the PB and DH approximations, but q^eff_i has a different value and we have in the denominator instead of ε_r. Note that ψ_i(r) ∼ q^eff_iϕ_Coul*(r) like in the PB approximation. Furthermore,which is symmetrical in i and j as it must. Note that the decay of h_ij(r) for large r is the same as for −βw_ij(r) because h_ij(r) = e^−βw_ij(r) − 1 ∼ −βw_ij(r) when βw_ij(r) is small. These asymptotic formulas constitute facts that are obtained from exact statistical mechanics for electrolytes. DIT provides methods to calculate q^eff_l and from the distribution functions of the system,** but we will not enter into any details here.

The appearance of the product q^eff_iq^eff_j in w_ij(r) of eqn (18) has some immediate, important consequences. Using eqn (17) we see that w_ij(r) ∼ q^eff_jψ_i(r) for large r, so the effective charge appears as the prefactor instead of the actual charge q_i that one has in the PB approximation w_ij(r) = q_jψ_i(r). Recall that the latter equation caused the incorrect symmetry of eqn (9). Since h_ij(r) ∼ −βw_ij(r) for large r we have from eqn (6)

Now, Poisson's equation and eqn (17) imply that for large r
By comparing the rhs of these two equations we see that the decay parameter κ is given bywhich is eqn (4). This is an exact equation that has been derived in a strict manner in DIT.^100,101 When the ion density goes to zero we have q^eff_j → q_j for all j, and κ → κ_D.

To conclude, the radial charge density decays like

The asymptotic expressions (16)–(18) and (20) do not depend on any particular approximation and are valid in general provided that the ion density is not too high. In order to treat higher densities, we must add other contributions to these expressions. This will be done next.

Eqn (18) only describes one decay mode of w_ij(r) of the electrolyte system. As mentioned in the Introduction, there exist several decay modes with different decay parameters, κ, κ′ etc. and we have^22,23,106 the likewise exact result

for r ≥ a, where κ < κ′ for low to moderate ion densities and “other terms” indicate terms from the additional decay modes as well as short-range terms with different decay behaviors.††Eqn (18) gives only the leading contribution with the longest decay length 1/κ. Each mode has its own values of the effective charges and the effective permittivity, like and . The decay parameter κ′ is given byanalogously to eqn (19). The screened Coulomb potential isand one may say that this function specifies the decay modes of the electrostatic interactions of an electrolyte.

As an alternative to eqn (19) and (22), one can express both of them as

whereby it is explicit that each mode has its own value of q^eff_j(κ). In fact, q^eff_j(κ) can be written as an integral that contains distribution functions for the electrolyte.^100,101 Since κ appears on both sides of eqn (24), this formula is an equation for κ that has several solutions, κ, κ′ etc., that constitute the decay modes.‡‡ This fact distinguishes eqn (24) from eqn (1) in a particularly important way. The latter equation gives κ_D uniquely and explicitly in terms of the system parameters q_j, n_j, T and ε_r, so there is one single decay mode in the PB and DH approximations.

The various decay modes in eqn (21) and (23) give the following principal contributions to the electrostatic potential and charge density for r ≥ a

andThe first term on the rhs of each equation gives the leading contributions (17) and (20), respectively, for large r when the ion density is sufficiently low.

For elevated ion densities, the leading decay of w_ij(r), ϕ_Coul*(r), ψ_i(r) and ρ_i(r) for large r become exponentially damped oscillatory at a Kirkwood cross-over point. This is included in eqn (21), (23), (25) and (26) and corresponds, as we have seen, to a pair of decay parameters that become complex-valued, κ = κ_ℜ+ iκ_ℑ and . This behavior cannot be obtained in the PB approximation since κ_D is a real number. The entities q^eff_i and also become complex-valued at the cross-over point and the corresponding primed quantities become the complex conjugates of the unprimed ones. Beyond the Kirkwood point, the two first terms in eqn (26) yield when r → ∞

where we have written , q^eff_i = |q^eff_i|e^−iη_i and κ = |κ|e^−iθ with real ϑ, η_i and θ and where we have defined α_i = η_i + 2θ − ϑ.

One can show^23,101 that before the Kirkwood crossover point we have and and that and go to zero at that point. Thus goes from ε_r to 0 when the ion density varies from zero to the Kirkwood point. The sum of the two first terms in the rhs of eqn (26) remains, however, finite at that point since the two infinities from and in the respective term have different signs and cancel each other.

In this paper we will explore approximations that express the electrostatic correlations in terms of the decay modes with different κ values. We will thereby make the following Ansatz for the electrostatic part w^el_ij of the potential of mean force

so the two leading modes in eqn (21) are included (we will also investigate cases where only one mode is included). This means that we assume that these contributions to w^el_ij are valid for all r down to contact, r ≥ a, and we neglect all other terms in w^el_ij, which decay faster than the leading mode. The main contributions of these other terms appear for small r, so we expect deviations from the exact w_ij(r) for small r values. Our Ansatz will work for cases where these deviations are sufficiently small. The complete w_ij(r) used in the current work also contains a contribution w^core_ij(r) from hard-core correlations among the ions. It will be specified in Section 4.2 and is important only for high densities. In the development of the theory we will gradually introduce various levels of refinements, whereby we will start with ρ_i(r) given by the first two decay mode terms of eqn (26) and later link this to w^el_ij given by eqn (28), which constitutes the basic approximation used in our approach.

Our approximation (28) is similar to the so-called exponential Debye–Hückel (DHX) approximation,^71,107 where the rhs of eqn (9) is used as an approximation for w_ij(r) in g_ij(r) = e^−βw_ij(r), that is,

whereby the DH value of q^eff_i from eqn (12) is used. The DHX approximation violates several necessary requirements of statistical mechanics for electrolytes, but it has nevertheless some merits.^{8,16,71,75,76}

Our Ansatz (28) is accordingly an extension of the DHX approximation to several decay modes. A major difference between it and the DHX approximation is that the effective entities q^eff_i, etc. in eqn (28) can be determined in a self-consistent manner from requirements that need to be fulfilled. In contrast, q^eff_i in the DHX approximation is explicitly determined by eqn (12) from the system parameters like density, temperature etc., which leads to the violations mentioned earlier.

In the main approximation of the current paper, the electrostatic response is linearized. This approximation will be denoted as the Multiple-Decay Extended Debye–Hückel, MDE-DH, approximation and is used in the first part of the current work, before Section 5. The linearization motivates the designation “an extended Debye–Hückel approximation,” since the DH approximation is a linear approximation for the electrostatic response. In the last part of the work, nonlinear contributions are studied and an extended DHX approximation is used.

Approximations that explicitly deal with several decay modes have been proposed earlier. Xiao and Song^48–50 have developed a linear theory for electrolytes that they designate as a molecular Debye–Hückel theory (they use the acronym MDH, which should not be confused with MDH used in the current work). This theory, which is based on DIT, focuses on the decay modes of ψ_i(r) as expressed in eqn (25) with several Yukawa function terms present, one for each decay mode with decay parameter κ_l. These modes are obtained from the static dielectric function [small epsilon, Greek, tilde](k), where k is the wave number, via an equation that is equivalent to eqn (24).§§ They write ψ_i(r) in terms of the Yukawa functions of all decay modes, for r ≥ a, where l in principle runs over all modes, whereby they neglect all other contributions to ψ_i. In practice the sum is limited to a few modes. The coefficients C_i,l are determined from equations that use the dielectric function [small epsilon, Greek, tilde](k) as input, whereby the latter can be taken from various theories like the MSA, LMPB, MDH or HNC approximations. The results are then used to calculate thermodynamical properties of the system. This approximation is closely related to the simplest version of the MDE-DH approximation of the current work, but it is not identical because the latter contains nonlinear contributions, which are important in the pair distributions and for thermodynamical consistency, as will be seen.

Another approach is that of Outhwaite and Bhuiyan,³¹ who have used the LMPB approximation to calculate ψ_i(r) in manner that also shares features with the simplest version of the MDE-DH approximation in the present work. They obtain κ and κ′ as the two smallest solutions of the LMPB equation for the decay parameters, whereby ψ_i(r) is approximated for r ≥ 2a like in eqn (3) and obtained for r < 2a by other means. The coefficients that correspond to C and C′ are determined via various alternative conditions on ψ_i, where one alternative is essentially the same as in the current work.

3 Modified Debye–Hückel approximations with a single decay-length

We saw in the Introduction that several approaches with a single decay lengths go quite a long way towards the objective to accurately obtain thermodynamical entities for electrolytes in a simple manner. In order to explicitly see why a single mode nevertheless is not sufficient in general, we will investigate two such approaches that are particularly relevant for the understanding of the approximations suggested in the current work.

As a preparation for the Multiple-Decay Extended DH approximation, let us consider the Modified Debye–Hückel (MDH) approximation by Kjellander,³² where one uses the following Ansatz with a single decay mode

which differs from the DH result in eqn (14) solely by the values of the decay parameter κ ≠ κ_D and the effective charge q^eff_i. This Ansatz means that the leading term in the large r decay formula, as given in eqn (20), is assumed to be valid for all r outside the ion and that , so the value of for low ion densities is assumed. Like the DH approximation, MDH is a linear approximation where the charge density of the ion cloud around an ion is proportional to the charge q_i. It can therefore be valid only for sufficiently low electrostatic coupling (high temperatures and/or high ε_r). By inserting eqn (30) into the condition of local electroneutrality we obtain
which gives(cf.eqn (12)). The mean electrostatic potential ψ_i(r) equalsas follows for Poisson's equation.

The deviation of q^eff_i from q_i expressed by eqn (31) is caused simply by the excluded volume that the ionic core gives rise to in the radial charge distribution. When all ions are treated on the same basis so all have q^eff_j ≠ q_j, this deviation has quite far-reaching effects as we now will see. By inserting eqn (31) into the exact eqn (19) and using eqn (1) we obtain the MDH result³²

This equation for κ has also been obtained in a different manner by Hall³³ in his LT approximation and recently by Ding et al. in their the charge renormalization theory.⁴⁰ It can be solved numerically to give κ as function of κ_D and when this solution inserted into eqn (31), one obtains q^eff as a function of κ_D. For low ion densities, where κ ≈ κ_D, eqn (30) agrees with the Debye–Hückel expression (11) for ρ_i, but it differs otherwise.

In order to illustrate how the solution to eqn (33) behaves, let us write eqn (33) as f(κa) = κ_Da, where f(x) = [x²(1 + x)/e^x]^1/2 and plot f as function of κa. This function is plotted as the full curve in Fig. 1 and the solution is the intersection point between the curve and a horizontal line at f = κ_Da; the figure shows an example with κ_Da = 0.477 for which the intersection (full circle) occurs at κa = 0.500. In the figure we see that there is a second intersection point, which is shown as an open circle. This point also corresponds to a solution of eqn (33), which hence has two solutions for κ which we will denote as κ and κ′, respectively. These two solutions exists for any κ_Da < 1.346, which is the maximum of the curve; the latter occurs for and is marked by a full square in the figure. The function f(x) goes to zero when x → ∞, so the right-hand intersection occurs for increasingly large κa when κ_Da goes to zero, which means that κ′ → ∞ in this limit. As we will see shortly, this solution is physically relevant and gives rise to a second decay mode.

Fig. 1 Plots of f(κa) for two different functions: the full curve shows f(x) = [x²(1 + x)/e^x]^1/2 and the dashed curve shows f(x) = [x²(1 + x)/exp₃(x)]^1/2, where exp₃(x) is an exponential sum function defined in eqn (37). The horizontal line shows f = κ_Da with κ_Da = 0.477. This line has two intersections with the full curve that are shown as a full circle and an open circle, respectively. The two intersections correspond to the two solutions κ and κ′ of eqn (33) when κ_Da = 0.477. The maximum of the full curve is marked with a full square.

We have thus shown that eqn (33) has two solutions κ and κ′ with κ < κ′, where κ goes to κ_D and κ′ goes to infinity when κ_Da goes to zero, or in other words, when the ion density goes to zero. The two solutions κ and κ′ appear despite the fact that the MDH theory was set up by using the ansatz in eqn (30) with only one κ parameter. This is an inconsistency that we will resolve later. When the ion density and hence κ_Da are increased, κ and κ′ approach each other and beyond κ_Da = 1.35, i.e., above the maximum in Fig. 1, the solutions κ and κ′ to eqn (33) become complex-valued, κ = κ_ℜ + iκ_ℑ and , which is in agreement with the general result mentioned above in Section 2.¶¶

The values of κ/κ_D, κ′/κ_D, κ_ℜ/κ_D and κ_ℑ/κ_D from the solutions to eqn (33) are plotted as functions of κ_Da in Fig. 2a and are compared with results from Monte Carlo (MC) simulations and Hypernetted Chain (HNC) approximation calculations of the decay of ψ_i(r) and ρ_i(r) for a 1 : 1 electrolyte in water at room temperature with a = 4.6 Å. These results also correspond to a classical 1 : 1 plasma in vacuum (ε_r = 1) when T = 23 400 K, which has the same value of k_BTε_r. The HNC approximation is very accurate for these systems. It is seen in the figure that the predictions from eqn (33) agree very well with the MC and HNC results, including the Kirkwood crossover point and beyond.|||| This is quite remarkable considering the humble origin of eqn (33). Obviously, the second root of eqn (33) and the cross-over to oscillatory decay have a great physical relevance. We will explore this in the next section. Fig. 2b shows a plot of κ_ℜ/κ_D and κ_ℑ/κ_D from eqn (33) for larger values of aκ_D. Other linear approximations like the MSA and LMPB theories, which are considerably more complicated, give predictions that are very similar to those shown in the figure.

Fig. 2 (a) Decay parameters divided by the Debye parameter κ_D for a 1 : 1 aqueous electrolyte solution at room temperature as functions of κ_Da, which is proportional to the square root of the ion density. The same results apply to a classical 1 : 1 plasma in vacuum at T = 23 400 K. Thick curves: obtained from the solutions of eqn (33); thin curves: HNC results and open symbols: MC simulation results. The Kirkwood crossover point is shown as a filled symbol for the HNC results and those from eqn (33). The MC data are taken from ref. 18 and the HNC data from ref. 21 and 22. (b) The real and imaginary parts of κ from eqn (33) divided by κ_D.

Apart from the local electroneutrality condition that the charge distribution has to fulfill, there is also the Stillinger–Lovett second moment condition¹⁰⁸ that must also be fulfilled. This condition expresses the fact that from the point of view of electrostatics, electrolytes behave like perfect conductors. It can be formulated as

which can alternatively be written in a more common manner as an equation involving the second moment of the pair correlation functions h_ij(r) = g_ij(r) − 1. This condition is not fulfilled for the MDH approximation [i.e., eqn (30) with κ given by eqn (33)], but this approximation does, however, fulfill it approximately for low to medium concentrations and is much better in this respect than the DH approximation.³³

We may strictly enforce the Stillinger–Lovett condition by refraining from the assumption that and select the value of so that eqn (34) is satisfied. Guided by eqn (20) we then take

If we insert this into eqn (34) and use the exact eqn (19) we obtainwhere we have defined the “exponential sum functions”(note that exp_∞(x) = exp(x) = e^x). Eqn (36) implies that for κa < 1, so it is a good approximation to set for low to medium ion densities as done in the MDH approximation.

By applying the local electroneutrality to eqn (35) we obtain

where we have inserted from eqn (36) to get the last equality. Using eqn (19), we obtainThis equation for κ is the Self-Consistent Screening Length (SCSL) approximation by Attard⁸⁴ mentioned earlier, which was set up as an approximation that satisfies both the local electroneutrality and the second moment conditions. In order to investigate κ obtained in this approximation, let us this time define f(x) = [x²(1 + x)/exp₃(x)]^1/2 whereby eqn (38) corresponds to f(κa) = κ_Da. In Fig. 1 the plot of this f as function of κa is shown as the dashed curve and we see that there is only one intersection with the horizontal line f(κa) = κ_Da. Therefore we see that this approximation does not predict a Kirkwood crossover point. Instead, the single real solution κ goes to infinity when κ_Da → 6, which follows from the fact that f(x) → 6 when x → ∞. For small κ_Da, however, the results of this approximation are accurate since the full and the dashed curves virtually coincide there. When κ_Da > 6, eqn (38) has only complex-valued solutions.

This is a clear example of the fact that one may not achieve an improvement of an approximate theory by enforcing a necessary condition. The problem is that it is too restrictive to assume that there only exist one decay mode e^−κr/r for the charge density, so the enforcement of the Stillinger–Lovett condition leads to a theory that gives qualitatively incorrect results for high ion densities. In order to obtain a Kirkwood crossover point we need two decay modes e^−κr/r and e^−κ′r/r. Thereby the oscillatory decay for high ion densities can be obtained in a correct manner. The MDH and SCSL approximations with a single decay mode can be used only for low to moderate ion densities. We will return to the SCSL approximation in Section 5; here we will continue with MDH approximation that is more accurate.

The MDH Ansatz (30) corresponds, as we will verify shortly, to the approximations w_ij(r) = q^eff_jψ_i(r) and g_ij(r) = e^−βw_ij(r) ≈ 1 − βw_ij(r), so by inserting ψ_i from eqn (32) we obtain

In the DH approximation, where ψ_i is given by eqn (13), we have insteadNote that h_ij(r) = g_ij(r) − 1 is proportional to q_iq_j in both cases since the effective charge is proportional to the actual (bare) change in these approximations.

The g_ij(r) in eqn (39) gives the charge density in eqn (30), as can be easily verified as follows. We have

where we have used the exact eqn (19) for κ and the fact that . This is the same as the second line in eqn (30). In the DH approximation we instead obtain the charge density given by eqn (14) for r ≥ awhere we have used the definition of κ_D in eqn (1). We can see that the appearance of a decay parameter κ ≠ κ_D in the MDH approximation is intimately linked to the treatment of all ions on the same basis, so all ions of species j have an effective charge q^eff_j ≠ q_j. As pointed our earlier, in the nonlinear Poisson–Boltzmann and the DH approximations, the central ion with effective charge q^eff_i ≠ q_i is treated differently from the ions in its surroundings. The latter are treated as point ions that do not correlate with each other, so for them we can say that q^eff_j = q_j and by inserting this into eqn (19) we obtain the DH expression (1) for κ_D.

The total ion density around an i-ion is equal to . In the DH and MDH approximations, this density is equal to the bulk value n_tot for r ≥ a because it follows from eqn (39) and (40) that .

4 Debye–Hückel extensions with multiple decay lengths

4.1 Decay parameters and effective permittivities

To resolve the inconsistencies of the single decay mode approaches, we use eqn (26) as basis for a better approximation. We take the first two terms as an Ansatz that we assume holds for all r ≥ awhere the second line is equal to ρ_i(r) for r ≥ a. We will, however, keep eqn (31) and hence eqn (33) as they are because the latter gives good results for κ and κ′, including the Kirkwood crossover, for monovalent ions in water at room temperature and for the classical plasma at high temperatures. These equations are also fulfilled by the primed quantities and the decay parameters used are those plotted in Fig. 2. Recall that the behavior of the decay modes as seen in this figure was obtained from an effect of excluded volume in the radial ion distribution, as expressed in terms of an effective charge in eqn (31) for each mode.

The Ansatz in eqn (43) constitutes an important ingredient of the simplest version of the MDE-DH approximation, which we will denote as the Simple MDE-DH approximation. We will later formulate the Complete MDE-DH approximation [see eqn (54), (55) and (97)], but eqn (43) is sufficient in many cases and works very well provided that the ion size is not too large and core–core correlations are not too prominent.

By applying the local electroneutrality to eqn (43) and using eqn (31) and its analogue for the primed quantities, we can readily deduce that

To determine and we need one more equation, which we can obtain by applying the Stillinger–Lovett condition to the ansatz in eqn (43) and use the exact eqn (19) and (22). This yieldsTogether, these two equations for and imply thatOne can show that these relationships imply that and , so is negative in agreement with the general results mentioned in Section 2. The value zero occurs at the Kirkwood point, that is,
as follows from eqn (46) and (47). This is also in agreement with the general exact results mentioned in Section 2. In the limit of infinite dilution where and , the second term in eqn (43) vanishes and eqn (30) becomes valid, but close to the Kirkwood cross-over point this second term is about equally important as the first term.

In the top frame of Fig. 3 the values of from eqn (46) are compared to results from MC simulations and the HNC approximation for the same system as in Fig. 2. The agreement is good; the curve from the present theory is somewhat shifted horizontally relative to that from HNC due the the slight difference in Kirkwood point as seen in Fig. 2. The prediction for is plotted in the bottom frame of Fig. 3.

Fig. 3 The ratios (top frame) and (bottom frame) as functions of κ_Da. The full curves show the result from eqn (46), the dashed curve show HNC results and the symbols show MC simulation results. The HNC and MC data are taken from ref. 18 for the same system as in Fig. 2 (see ref. 21 and 22 for further HNC data). Note the difference in ordinate scales of the two frames.

For the case where κ is complex-valued, κ = κ_ℜ + iκ_ℑ, we write and . One then obtains from eqn (44) that

and one can derive from eqn (46) thatwhere
In Fig. 4 we have plotted and ϑ as functions of κ_Da. At the Kirkwood crossover point, to the left in the figure, we have and ϑ = π/2.

Fig. 4 The ratio and the phase angle ϑ of from eqn (48) and (49) as functions of κ_Da beyond the Kirkwood crossover point where κ is complex-valued.

For complex-valued κ and κ′ the radial charge distribution is oscillatory and the second line of eqn (43) can be written as (cf.eqn (27))

where α_i = η_i + 2θ − ϑ with phases defined from , q^eff_i = |q^eff_i|e^−iη_i and κ = |κ|e^−iθ as before. Note that θ = −arctan(κ_ℑ/κ_ℜ).

4.2 Radial distribution functions

For binary symmetric electrolytes, the number densities of cations and anions are equal, n₊ = n₋ ≡ n, the ionic charges are q₊ = −q₋ ≡ q, the absolute valency z = q/q_e and n_tot = 2n. For ions of equal diameter (the restricted primitive model, RPM) we also have ρ₊(r) = −ρ₋(r) ≡ ρ(r) and q^eff₊ = −q^eff₋ ≡ q^eff. From eqn (1) and (19) followwhich is exact for the RPM.

In the RPM, the electrolyte bulk phase is completely specified by the two dimensionless parameters β_R ≡ βq²/(4πε_rε₀a) = z²[small script l]_B/a and τ = κ_Da; the latter is the abscissa of Fig. 2–4. Alternatively, any two independent combination of these two parameters can be used to specify an RPM system, like reduced total ionic density n_totR ≡ n_tota³ = τ²/(4πβ_R) and the product β_Rτ. The inverse of the latter is the so-called plasma parameter 4πn_tot/κ_D³.

We will test the predictions of the Simple MDE-DH Ansatz in eqn (43) by comparing with results from MC simulations and HNC calculations for an electrolyte with β_R = 1.55, which can be, for example, a 1 : 1 aqueous electrolyte solution at room temperature with ions of diameter a = 4.6 Å or a classical plasma of such ions in vacuum at T = 23 400 K. Depending on the ion density (concentration), the parameter τ has different values. The radial charge distribution ρ(r) is calculated using κ and κ′ obtained from the solutions of eqn (33), obtained from eqn (46) and from eqn (47). In Fig. 5 the function r|ρ(r)| as function of r for the cases of 0.5 and 0.7 M electrolytes is plotted on a semilogarithmic scale, whereby a Yukawa function decay becomes a straight line.

Fig. 5 The radial charge distribution ρ(r) for a 1 : 1 electrolyte plotted as rρ(r) as a function of r on a semilogarithmic scale. The ion diameter is a = 4.6 Å and the concentration is 0.5 M in frame (a) and 0.7 M in frame (b) [τ = κ_Da = 1.07 and 1.27, respectively, and β_R = 1.55]. The open squares in frame (a) show results from MC simulations¹⁸ and in both frames the full curves show HNC results,^18,21,22 the dotted curves show results from the Simple MDE-DH approximation, eqn (43), and the dashed curves show results from solely the first (leading) term in this equation. The unit for rρ(r) is q_ea⁻².

In Fig. 5a we see that the MC and HNC results virtually coincide with each other, so the HNC approximation is very accurate for the monovalent electrolyte. The prediction from eqn (43) (dotted curve) agrees very well with the MC and HNC results for all r except in a very short interval just outside the core region, 1 ≤ r ≲ 1.3a and the curves virtually coincide for r ≳ 2a. In Fig. 5b the agreement is equally good. The dashed curves are calculated from the first term in eqn (43), which is leading for large r. In Fig. 5b this curve coincides with the accurate results r ≳ 2.5a, but it deviates substantially from the dotted curve for smaller r values. This implies that the contribution from the second Yukawa term with decay parameter κ′ is important for the 0.7 M case. For the 0.5 M it is less important, but it still gives a noticeable contribution that makes the dotted curve to deviate from the dashed one. The difference in importance of the second term for the 0.7 and 0.5 M cases is due to the fact that the magnitude of in the denominator of this term increases rapidly when the concentration is decreased, see Fig. 3. We can conclude that the Ansatz (43) gives very good results for the 1 : 1 electrolyte systems; the deviation for small r has little consequence in most cases and is caused by terms in ρ(r) with shorter decay lengths that are not included in the Ansatz.

This conclusion is supported by Fig. 6, which shows results for 0.1, 0.7 and 1.0 M electrolytes (the data for 0.7 M are the same as in Fig. 5b). Here, ρ(r) is plotted on a linear scale as 4πr²ρ(r), which is proportional to the amount of charge at distance r. The curves for 0.1 M show that the first term in eqn (43) dominates for nearly all distances; the second term is very small because has very large magnitude when the concentration is low. The deviation from the HNC curve occurs, like for the other cases, mainly very close to the core region and has little consequence.

Fig. 6 The radial charge density ρ(r) for a 1 : 1 electrolyte plotted as 4πr²ρ(r) as function of r. The ion diameter is a = 4.6 Å and the concentration is 0.1 M (red curves), 0.7 M (black curves) and 1.0 M (blue curves) [τ = κ_Da = 0.48, 1.27 and 1.51, respectively, and β_R = 1.55]. For the cases 0.7 M (same as in Fig. 5b) and 0.1 M, the full curves show HNC results, the dotted curves show results from the Simple MDE-DH approximation, eqn (43), and the dashed curves show results from the first term in this equation. In the case of 1.0 M, the full curve shows the HNC data and the dotted curve show the results from eqn (50). The bottom frame shows a magnified view of the top one. The HNC data are from ref. 21 and 22.

For the 0.7 M case in Fig. 6a, we can very clearly see the importance of the second Yukawa term, since without it, the deviation from the accurate HNC curve is large for r < 2a as apparent from the dashed black curve. Thus it is very important to have the two decay modes present in the theory in order to describe the electrolyte in an appropriate manner. This conclusion is, of course, further emphasized when κ and κ′ and become complex-valued at the Kirkwood cross-over point and give rise to oscillations at higher concentrations. The conditions for the 1.0 M case in Fig. 6 are beyond this point, so ρ(r) has an exponentially damped oscillatory decay. We then have complex and the two terms in eqn (43) yield the expression for ρ in eqn (50). It is seen from the results from this equation, shown as the red dotted curve in the figure, that the agreement with the accurate HNC curve is about equally good as the results for the other concentrations. Note that the two red curves cross the abscissa axis at virtually the same point in Fig. 6b, which means that the phase α obtained from the Ansatz is very good.

Let us now turn to the pair distribution function g_ij(r) = e^−βw_ij(r). The Ansatz (43) can be obtained from eqn (28) for the electrostatic part, w^el_ij, of the potential of mean force. When βw_ij(r) is sufficiently small we can linearize the pair distribution function g_ij(r) ≈ 1 − βw_ij(r) and when the electrostatic contributions to w_ij dominate we can set w_ij(r) ≈ w^el_ij(r), so we have . By inserting eqn (28) and using eqn (19), we obtain eqn (43) [compare with the derivation of eqn (41)]. As we will see shortly, in the RPM the Ansatz (43) is valid also for somewhat more general conditions.

When the ion density is high, the approximation w_ij(r) ≈ w^el_ij(r) is insufficient because there are important steric contributions missing due to core–core correlations. The latter are in general coupled to the electrostatic correlations, so for simplicity, we will restrict ourselves to symmetric electrolytes in the RPM, where the density-charge correlation function is identically equal to zero and the coupling between electrostatic and core correlations is weak. We make the following Ansatz

and we approximate w^core_ij(r) with the potential of mean force obtained from a hard sphere fluid with pair distribution function g^hs(r) = e^−βwhs(r). The density of the hard sphere fluid is set equal to the total ionic density n_tot, that is, the hard sphere fluid has a reduced density n_tota³. A suitable choice of g^hs(r) in numerical calculations is to use the very accurate Verlet–Weis' semi-empirical expression¹⁰⁹ for g^hs(r). A Fortran code for this expression is publicly available.¹¹⁰ The g^hs(r) thus obtained gives contact values g^hs(a) in agreement with the very accurate Carnahan–Starling (CS) equation of state.⁹⁷ The contribution w^core_ij(r) expresses a further decay mode for the interactions.

Eqn (52) together with the two-mode Ansatz for w^el_ij(r) in eqn (28) constitute a starting point for the formal development of the approximations used in the current work. The details are given in Appendix A. Here we will only give the main points.

In the RPM, the Ansatz (52) implies that

g_+±(r) = e^{−β[±wel(r)+wcore(r)]} = e^∓βwel(r)g^core(r) (53)

with g^core(r) = g^hs(r) and we have for r ≥ awhere we have used eqn (51) to obtain the second line.

In the MDE-DH approximations we assume that βw^el(r) is sufficiently small so that we can approximate e^∓βwel(r) ≈ 1∓βw^el(r) + [βw^el(r)]²/2. We obtain

which can be written aswhere h^core(r) = g^core(r) − 1. The inclusion of the square term [βw^el(r)]²/2 makes the pair distribution to fulfill the necessary condition g_+±(r) ≥ 0, which can be violated in linear theories like MSA and the DH approximation. Eqn (54) and (55) constitute the Complete MDE-DH approximation for the RPM together with the eqn (31) for κ, which can be writtenand the corresponding relationship for the primed quantities.

The inclusion of a square term [βw^el_ij(r)]²/2 from electrostatics has been done for a long time in some kinds of theories, for example as an improvement of the DH approximation with w^el_ij(r) given by eqn (29).¹⁰ In the present case, the use of the square term of w^el_ij(r) with two decay modes and with coefficients that makes g_ij(r) satisfy the required electroneutrality and second moment conditions gives, as we will see in the next section, thermodynamic consistency to a considerable degree. In, for example, the GMSA approximation^24,47 one instead includes an empirical contribution with parameters that are selected to give thermodynamic consistency. The present approach does not contain any such parameters.

In systems where the core–core correlations are not very strong, it is a good approximation to neglect the cross-terms between h^core and w^el in eqn (56) since they are of higher order and are smaller than the other terms for most separations r in such cases. We then obtain for r ≥ a

Together with eqn (54) and (57), this equation constitutes the Simple MDE-DH approximation for the RPM. As shown in Appendix A, this is the approximation behind the Ansatz (43) in the RPM because eqn (58) yields ρ(r) used in the Ansatz (note that ρ(r) = qn[g₊₊(r) −g₊₋(r)] in the RPM). The square term and h^core do not contribute to ρ and therefore not to ψ. In this approximation all results in Section 4.1 are valid, in particular the expressions for and in eqn (44)–(49). The corresponding expressions in the Complete MDE-DH approximation are given in Appendices A and B.

Since we have retained solely the linear and quadratic terms in βw^el in both MDE-DH approximations, their range of validity is clearly limited to systems at sufficiently high temperatures and/or high ε_r so that β_R is sufficiently low. For monovalent ions of diameter 4.6 Å in aqueous solution at room temperature, which constitute the system for the MC simulations that we successfully have compared with in the previous examples, we have β_R = 1.55. The appropriate β_R values for the MDE-DH approximations apparently are of this order of magnitude or less.

As a test of the Complete MDE-DH approximation, the radial distribution function g_ij(r) has been calculated for a system with rather large ions, a = 6.6 Å, and high electrolyte concentration, where the core–core correlations are important and the charge density distribution ρ(r) is oscillatory. For this system β_R = 1.10. The results for g_ij and Δg = g_+ − − g_− − are presented in Fig. 7, where they are compared to MC simulations results of Zwanikken and coworkers.⁶¹ As seen in the figure, the MDE-DH approximation gives excellent results that nearly coincide with the MC ones. Note that the resolution of the ordinate scale in frame (d) of the figure is an order of magnitude larger than in frame (a). It is noteworthy that quite intricate details of the pair distribution functions are obtained in agreement with simulations. Our results are, in fact, slightly better than those of both the HNC approximation and Zwanikken's DHEMSA approximation,⁶¹ which is quite remarkable considering the simplicity of the present Ansatz.

Fig. 7 (a) The radial distribution function g_ij(r) for a 2 M 1 : 1 electrolyte with ions of diameter a = 6.6 Å (τ = κ_Da = 3.11 and β_R = 1.10) calculated in the Complete MDE-DH approximation (black curves) compared with data from MC simulations (symbols) taken from ref. 61. The red curve (short dashes) shows the pair distribution g^hs(r) for a hard sphere fluid of the same density and sphere diameter (packing fraction η = 0.37). The insert (b) shows a magnified view of the g_ij(r) curves in frame (a). (c) The difference Δg(r) = g₊₋(r) − g_− −(r), which is proportional to the radial charge-density distribution ρ₋(r) around a negative ion. Full curve is the MDE-DH and symbols are the MC results [the full circles are from Δg(r) data of ref. 61 and the full squares are calculated from g_ij(r) in frame (a) taken from the same reference]. The red curve (short dashes) shows Δg^el(r) = 2βw^el(r), which is the electrostatic part of Δg in the present approximation. The insert (d) shows a magnified view of the curves in frame (c). Note the large differences in ordinate scales of the various frames.

In order to investigate the effects of core–core correlations in this system, the pair distribution g^hs(r) of the hard sphere fluid is plotted in Fig. 7a and it is seen that the main oscillation of g_ij(r) with a wave length ≈a is caused by these correlations. It is also seen in frames (a) and (b) of the figure that there is another oscillatory component with longer wave length superimposed on the main oscillation. This component is equal to the term ∓βw^el(r) in eqn (56), which is oscillatory here.*** It makes the g₊₋ and g₋₋ curves to repeatedly intersect with one another and has a wave length 2π/κ_ℑ = 2.06a. In frames (c) and (d) this component gives the contribution Δg^el(r) = 2βw^el(r), which is the electrostatic part of Δg in the current approximation. For r ≳ 1.3a it is seen that Δg^el(r) gives nearly the entire Δg(r).

As seen in eqn (56), g_ij(r) has also contributions from [βw^el(r)]²/2 and from the two cross-terms involving h^core and w^el. All these contributions have shorter decay lengths than the other ones, so they mainly contribute for small r.

For the present system we have and ϑ = −0.828 radians as obtained from eqn (109) and (110) in Appendix B. This can be compared with the values obtained from the Simple MDE-DH approximation, eqn (48) and (49), which are 1.77 and −0.477, respectively. Eqn (58) of the latter approximation gives g_ij functions for this system that do not agree well with the MC data in Fig. 7. This shows the importance of using the Complete MDE-DH approximation here. All other numerical results obtained so far and those in Section 4.3 below are obtained for systems with smaller ions where it is sufficient to use the Simple MDE-DH approximation based on the Ansatz (43) and eqn (54), (57) and (58).

4.3 Thermodynamic quantities

Let us now turn to the thermodynamic properties and we start with the chemical potential μ_i of an ion of species i

μ_i = μ^ideal_i + μ^ex_i = k_BT ln(λ_i³n_i) + μ^ex_i,

where λ_i is the thermal de Broglie wave length for species i and μ^ex_i is the excess chemical potential. The latter has two parts μ^ex_i = μ^ex,core_i + μ^ex,el_i, where μ^ex,core_i is the contribution from the formation of a hole in the electrolyte that is large enough to fit an ion and μ^ex,el_i is the contribution from the electrostatic interactions. We define μ^ex,core_i as the change in free energy when an uncharged sphere of diameter a is inserted into the electrolyte and μ^ex,el_i as the free energy change when this sphere is charged from 0 to q_i. The average excess chemical potential μ^ex_± is given by where x_i = n_i/n_tot. The mean activity coefficient γ_± is defined from ln γ_± ≡ μ^ex_±/k_BT and we likewise define ln γ^el_± ≡ μ^ex,el_±/k_BT and ln γ^core_± ≡ μ^ex,core_±/k_BT. The latter is equal to μ^ex,core_i/k_BT for any i because all ions have the same size. Obviously, ln γ_± = ln γ^el_± + ln γ^core_±.

Let us start with the core term and consider the excluded volume hole in the electrolyte where an uncharged sphere of diameter a can fit. This volume, where the centers of the surrounding ions cannot enter, is equal to 4πa³/3. In the DH and MDH approximations

which can be understood from the fact that k_BT ln γ^core_± = k_BTn_tot × 4πa³/3 is the reversible pressure-volume work done when a hole of radius a and volume 4πa³/3 is formed in a low-density fluid and the density around the hole is equal to n_tot throughout. As we have seen, in these approximations the total ion density around an ion is equal to the bulk value n_tot for all distances from the ion and the same applies to an uncharged sphere or a hole.

For high ion densities it is, of course, not reasonable to evaluate ln γ^core_± as if the total density around an uncharged sphere is unaffected by the interactions with the latter. In the MDE-DH theories for the RPM, an uncharged sphere does not interact electrostatically with the ions because w^el is proportional to the charge of the central particle, but hard core interactions are included because of the presence of the g^core in the pair distribution function in eqn (53). Therefore, ln γ^core_± is equal to the excess chemical potential for the creation of a hole in a dense fluid of hard spheres of number density n^hs, which is set equal to the total ion density n_tot. Since we use the Verlet–Weis' semi-empirical expression for g^hs and set g^core = g^hs we have

where μ^ex_CS is given by the very accurate Carnahan–Starling expression for μ^ex,core obtained from the CS equation of state⁹⁷where η is the packing fraction η = πn^hsa³/6. The value of ln γ^core_± from eqn (60) approaches that of eqn (59) when the concentration is decreased to zero.

The electrostatic parts of the chemical potential and the activity coefficient can be calculated from the radial charge distribution functions of the electrolyte. In the DH approximation we have

and in the MDH approximation³²with κ from eqn (33). For a binary electrolyte where the ionic charges are q₊ and q₋, the prefactor of these expressions can alternatively be written as
where it is clearly seen that the prefactor is a constant that is independent of the ion density.

In the Simple MDE-DH approximation based on the Ansatz (43) and eqn (58), it is shown in Appendix C that

and we see that there is one contribution from each of the two decay modes with decay parameters κ and κ′. This expression is valid also for complex-valued κ and κ′, but for that case one can write it in a more convenient form given by eqn (114) in Appendix C.

In any approximate theory where the charge density ρ_i is proportional to q_i, like the DH, MDH and MDE-DH approximations, μ^ex,el_± is equal to the excess (interactional) average energy per ion, U^ex/N_tot (see eqn (115) ff. in Appendix C), where is the total number of ions and V is the volume of the system. We therefore have ln γ^el_± = βU^ex/N_tot in these approximations, so in the the Simple MDE-DH approximation we have

The equality between μ^ex,el_± and U^ex/N_tot is not true in general, but for conditions where linear electrostatic response is a good approximation we have ln γ^el_± ≈ βU^ex/N_tot.

In Fig. 8a, ln γ_± obtained in the Simple MDE-DH approximation is compared with results from MC simulations and we see that the agreement is very good for all concentrations; there are slight differences between the results for high concentrations that are hard to see in the plot. The red portions of the curves indicate the region where the radial charge distribution is oscillatory and the Kirkwood cross-over point is shown by red arrows. The DH prediction is also plotted and serves as a reference in the figure. The electrostatic contribution ln γ^el_± from these approximations is shown and these curves also give the predicted average energy plotted as βU^ex/N_tot. The latter quantity as obtained from MC simulations is plotted as blue triangles in the figure and the Simple MDE-DH approximation results are very good agreement with the MC data; the slight differences seen for high concentrations are actually nearly the same as for the ln γ_± results in the figure.

Fig. 8 (a) The logarithm of the mean activity coefficient γ_± and its electrostatic part γ^el_± plotted as functions of the square root of the salt concentration c_salt for 1 : 1 aqueous electrolyte solutions (a = 4.25 Å, T = 298 K, β_R = 1.68). The full and dotted curves show ln γ_± from the Simple MDE-DH approximation; the full curve is obtained directly and the dotted one is obtained via an integration of the osmotic coefficient. The symbols show the results from MC simulations taken from ref. 77. The alternatingly dashed-minidashed curve shows ln γ^el_± = βU^ex/N_tot from the Simple MDE-DH approximation (in this approximation μ^ex,el_± = k_BT ln γ^el_± is equal to U^ex/N_tot). The blue triangles show βU^ex/N_tot from the same MC simulations (ln γ^el_± is in general not equal to βU^ex/N_tot in simulations). The Kirkwood cross-over point (K) between monotonic and oscillatory decay is indicated by red arrows; for the red portions of the curves the radial charge distribution is oscillatory. For reference, we also show the Debye–Hückel (DH) predictions for ln γ_± and ln γ^el_± = βU^ex/N_tot as dashed curves. (b) The osmotic coefficient ϕ plotted in the same manner. The full curve shows Simple MDE-DH, symbols MC and dashes DH results.

Fig. 9a shows the low concentration region in more detail, where the results from the MDH approximation are also included. The latter cannot give any prediction at high ion densities (beyond the Kirkwood cross-over point) because there is only one term for the Ansatz in eqn (30) and, as seen in the figure, the MDH approximation works well only for concentrations up to about 0.25 M (c^1/2_salt ≈ 0.5 M^1/2). Again we see that it is very important to have the two decay modes present in the theory in order to describe the electrolyte in an appropriate manner. The dashed-dotted curve is calculated by using ln γ^core_± from eqn (59) rather than from eqn (60). By comparing the full and the dashed-dotted curves we can see that for concentrations below about 0.25 M (c^1/2_salt ≈ 0.5 M^1/2), it it does not matter which expression one uses, but for higher concentrations one must use the CS expression in order to obtain the correct values for ln γ_±.

Fig. 9 The quantities (a) ln γ_± and (b) ϕ for the same system as in Fig. 8, but displayed in an expanded view for low concentrations. The full, dotted and short-dashed curves and the open circles are the same as in Fig. 8. In frame (a), the crosses show MC simulation data from ref. 72 and the long-dashed curve shows the result of the MDH approximation. The dashed-dotted curves in both frames are based on the Ansatz (43), like the full and dotted curves, but they are calculated with the same hard core contribution for ln γ_± and ϕ as the DH and MDH approximations. The red portions of the curves indicate the region where the radial charge distribution is oscillatory.

The osmotic coefficient ϕ is defined from ϕ = P/P^ideal where P is the (osmotic) pressure of the fluid and P^ideal = k_BTn_tot is the ideal pressure. The excess pressure is P^ex = P − P^ideal and the excess osmotic coefficient is given by ϕ^ex = P^ex/P^ideal = ϕ − 1 (see eqn (116) ff. in Appendix C). In the DH approximation we have

and in the MDH approximation³²where the first term in each rhs is the contribution ϕ^cont = P^cont/P^ideal from the contact pressure P^cont at the hard core surfaces of the ions, the second term is ϕ^el = P^el/P^ideal and P^el is the pressure from the Coulombic interactions in the electrolyte.

As shown in Appendix C, in the Simple MDE-DH approximation we have

where ϕ^core = 2πa³n_totg^core(a)/3, andSince g^core(r) = g^hs(r) we have ϕ^core = ϕ^ex,hs = ϕ^hs − 1, that is, the excess osmotic coefficient for a pure hard sphere fluid of density n^hs = n_tot. The use of the Verlet–Weis' expression for g^hs implies that this part of the osmotic coefficient is given by the Carnahan–Starling equation of state, which says that
so we have

ϕ^core = ϕ^hs_CS|_n^hs=ntot − 1 (70)

in eqn (68). When the ion density goes to zero, this ϕ^core goes to ϕ^cont of the DH and MDH approximations in eqn (66) and (67). The latter only include the core contribution to ϕ^cont.

The prediction for ϕ obtained in the Simple MDE-DH approximation is shown as the full curve in Fig. 8b and the low concentration region is shown in Fig. 9b. It is seen that this prediction is in very good agreement with the MC data for all concentrations. Again, the curve nearly coincide with the MC data.

The activity coefficient can alternatively be obtained from the osmotic coefficient via the relationship

which can be derived from the Gibbs–Duhem equation (the selected integration variable is suitable since ϕ^ex(n_tot′) is proportional to for small densities). In the integrand, ϕ^ex(n_tot′) is evaluated for an electrolyte with total density n_tot′ that varies from 0 to the final value while keeping the mole fraction x_i constant, whereby one sets the ion density equal to n_i′ = x_in_tot′ when calculating ϕ^ex.

In approximate theories, there is no guarantee that ln γ_± calculated from ϕ^exviaeqn (71) is the same as ln γ_± obtained earlier. In order to have consistency, one would need to obtain the same value of each quantity irrespectively of how it is calculated, but this is guaranteed only in exact theory. A well-known fact is that the DH approximation shows very poor consistency for ln γ_± and the same is true for the MDH approximation.³² The use of eqn (71) is discussed in Appendix C, where we see that it is a common feature of all linear approximations (where ρ_i is proportional q_i) that the contribution ϕ^el in ϕ^ex integrated according to eqn (71) does not give ln γ^el_±. In the DH and MDH approximations, this is the cause of the inconsistency for ln γ_±.

In the Simple MDE-DH approximation, the contribution ϕ^el does not give ln γ^el_±via this integration for the same reason. However, there is also an electrostatic contribution from the square term in ϕ^cont given by eqn (68) that originates from the square term in eqn (58). This contribution makes the situation different. In Fig. 8a the dotted curve shows ln γ_± obtained via the integration (71) from the osmotic coefficient in Fig. 8b. In Fig. 9a the low concentration region is shown in more detail. The dotted curve agrees well with the full curve, so we see that this approximation shows a good thermodynamic consistency in this respect. Considering that the theory is based on the very simple Ansatz (43) and eqn (58), this is an accomplishment. In fact, it is not uncommon for more advanced theories to lack such a consistency.

For the system in Fig. 7 with large ions, the Complete MDE-DH approximation given by eqn (54) and (55) is, as we have seen, needed for the distribution functions because the cross-terms in eqn (56) give non-negligible contributions. However, when the thermodynamic quantities of this system are calculated in this approximation, ln γ_± deviates only by 0.4% and ϕ by 0.9% from the predictions by the Simple MDE-DH approximation, so the latter can still be used as a good approximation for these quantities. This is a further example of the fact that thermodynamical quantities sometimes are not overly sensitive to the details of the structural entities like the distribution functions. The corresponding deviations for the contributions ln γ^el_± and ϕ^el are 5.5% for both.††† If the latter quantities are needed with a higher accuracy than this, one has, of course, to use the Complete MDE-DH approximation. The formulas for such calculations are presented in Appendix C, eqn (125)–(127).

5 Influence of terms with higher powers of w_ij

As we have seen, the approximations in this paper are based on the basic Ansatz (52)
where is set equal to g^hs and w^el_ij is given by the Ansatz (28) with two Yukawa function terms. The electrostatic part of g_ij can be expressed as a Taylor expansionFor the RPM, we have w^el_+±(r) = ±w^el(r) with w^el given by eqn (54). In the MDE-DH approximations we only keep terms up to and including the quadratic one in the expansion (72). This leads to eqn (55), which contains a square term that does not contribute to ρ(r) so the latter is linear in w^el_ij. In fact, for all numerical results presented so far for the monovalent ions, the inclusion of all nonlinear terms changes the results only marginally, so the MDE-DH approximation used is fully adequate.

There remains, however, one issue to be considered. For dilute electrolytes where the screening is low, w^el_ij(r) has a long range and the terms with higher powers than two have, in fact, a qualitative effect on the screening behavior. In the limit of infinite dilution we have the exact limiting law in eqn (5) for the screening parameter κ, which takes the following form for symmetric binary electrolytes‡‡‡

where z = q/q_e. This asymptotic result originates, in fact, from the term in eqn (72) with w^el_ij to the third power.§§§ Later we will derive an approximate extension, eqn (81), of this limiting law.

As we have seen, since the logarithm in eqn (73) is negative for small Λ, the second term in the rhs is negative there, so κ < κ_D for dilute symmetric electrolytes and there is accordingly under-screening. However, the approximate formula (33) for κ in the MDE-DH approximation implies over-screening, κ > κ_D, and when κ_D → 0 it gives (κ/κ_D)² ∼ 1 + (κa)²/2 ∼ 1 + (κ_Da)²/2. Therefore, this κ is qualitatively incorrect for very dilute solutions in this approximation. Other linear approximations for the electrostatic response in electrolytes, like MSA, LMPB and GDH, share the same features; they all predict that κ ≥ κ_D and a limiting form proportional to (κ_Da)².

As mentioned in the Introduction, it has been experimentally observed that κ < κ_D in dilute solutions of symmetric electrolytes at high electrostatic coupling, so this phenomenon is relevant to investigate. Since eqn (73) has general validity as an exact limiting law, it is applicable also for low electrostatic coupling so it is theoretically interesting to start with this latter case.

In order to analyze these matters, let us investigate the predictions of the full nonlinear expression for ρ(r) = qn[g₊₊(r) − g₊₋(r)] with g_ij(r) given by eqn (53) and (54), that is,

for r ≥ a and zero otherwise. At high dilution, the contribution from the second Yukawa term with decay parameter κ′ can be neglected because is very large there and g^core_ij(r) ≈ 1, so ρ(r) is given to a very good approximation byApproximations that are similar to eqn (75) have been investigated earlier by Attard⁸⁴ and Kjellander.³² Here the task is to use this equation to evaluate the behavior of κ and in eqn (74) for low electrolyte concentrations.

Using the exact eqn (51), we can write the coefficient in the prefactor of the argument of sinh in eqn (75) as

The local electroneutrality condition and the Stillinger–Lovett second moment condition can be expressed in the RPM as
and
respectively. By inserting eqn (75) into these equations and using n_tota³ = τ²/(4πβ_R), we obtain the following set of equations for the two unknowns κ/κ_D and for given dimensionless system parameters τ and β_Randwhere we have made the variable substitution R = r/a. This set of equations can be solved numerically for κ/κ_D and . It constitutes an approximate nonlinear theory that is valid for dilute electrolytes. We will denote this approximation as the Extended DHX (EDHX) approximation. It is the limiting form for infinite dilution of the nonlinear (DHX) version of the MDE-DH approximation, where ρ(r) is given by eqn (74) and g_+±(r) = e^∓βwel(r)g^core(r) with w^el from eqn (54) [cf. the discussion in Section 2 in connection to eqn (28) and (29)].

Before proceeding we note that the linearized version of this approximation is equal to the self-consistent screening length, SCSL, approximation by Attard, which was discussed in Section 3. This follows because linearization of eqn (75) gives [cf.eqn (35)]

and the application of the local electroneutrality and Stillinger–Lovett conditions yields (κ/κ_D)² given by eqn (38), which is the SCSL approximation. Furthermore, this linearized expression for ρ(r) yields given by eqn (36). Since the EDHX and SCSL approximations contain only one decay mode, they can be useful only for dilute electrolytes.

The numerical solution of eqn (76) and (77) of the EDHX approximation gives the following results. Fig. 10 shows the ratios κ²/κ_D² and as functions of τ = κ_Da for a monovalent electrolyte with a = 4.6 Å in aqueous solution at room temperature, which corresponds to β_R = 1.55 (recall that κ_Da is proportional to the square root of concentration). The predictions of the EDHX approximation are compared in the figure with results from MC simulations, the HNC, MDE-DH and SCSL approximations. In the insert to Fig. 10a, which shows a magnified view of the dilute region, we see that the EDHX approximation indeed gives κ < κ_D for the dilute electrolyte in agreement with the exact asymptotic formula and the HNC approximation. However, the slight dip of κ²/κ_D² below 1 occurs for small ion densities where the deviation of κ from κ_D is negligible in practice. Despite that the curves for MDE-DH and SCSL approximations in the insert show the qualitatively incorrect behavior with κ > κ_D, it only on this greatly magnified scale that this can be seen. Likewise, in the insert to Fig. 10b we can see that for dilute solutions as shown by the HNC and EDHX approximations, while the MDE-DH and SCSL approximations give there, which is incorrect. Also in this case, the incorrect behavior occurs for low densities where the deviation of from ε_r is negligible for most practical purposes.

Fig. 10 The ratios (a) κ²/κ_D² and (b) plotted as functions of κ_Da for monovalent electrolytes in aqueous solution at room temperature (β_R = 1.55). The respective insert shows the same plot in an expanded view for low concentrations (low κ_Da). The dot-dashed curves are the result from the EDHX approximation, eqn (76) and (77), and the dotted curves are from the SCSL approximation; both are drawn in magenta color. The full curves are the Simple MDE-DH results, the dashed curves and crosses show the results from HNC calculations (the crosses occur only in b), and open circles are from MC simulations. The HNC and MC data are taken from ref. 18 and 22.

In the main frames of Fig. 10 we see, as before, that the MDE-DH approximation gives a very good over-all account of the properties of the monovalent electrolyte. Furthermore, we can see that the EDHX and SCSL approximations do not work for higher concentrations, because, as pointed out earlier, they only contain one decay mode. In the figure one can clearly see that it is very important to include two decay modes with different decay parameters κ and κ′, as in the MDE-DH approximation, in order to get agreement with the accurate calculations outside the dilute region.

The fact that κ < κ_D for low concentrations, as shown by the asymptotic formula eqn (73), is a general feature for symmetric electrolytes. Due to the factor z⁴ in this equation, the under-screening will become much more pronounced when the valency of the ions is increased. This has been verified for divalent (2 : 2) electrolytes in ref. 18 and 22 by HNC and MC calculations; these HNC and MC results are shown in Fig. 11 together with the corresponding result from the EDHX approximation. As seen in Fig. 11a, we have κ < κ_D for the divalent case in a wide concentration interval and the negative deviation of κ²/κ_D² from 1 reaches about 30% (shown by the MC data) before it rises again. This highly significant difference from the monovalent case is solely due to the difference in ionic charge by a factor of two; everything else is the same. In this case the HNC approximation (dashed curve) does not provide accurate values of κ, except for low concentrations. As pointed out earlier, the EDHX approximation is applicable only for sufficiently dilute electrolytes, but for higher concentrations it nevertheless has the correct qualitative behavior in a part of the concentration interval shown.

Fig. 11 Frames (a) and (b) show the same kind of plot as Fig. 10a and b, but for divalent (2 : 2) electrolytes in aqueous solution at room temperature. The symbols show the results from MC simulations, the dashed curves from HNC calculations, the dot-dashed curves from the EDHX approximation that is valid only for low concentrations, and the red dotted curves show the asymptotes given by eqn (81) and (82), respectively. The ions have diameter a = 4.6 Å and β_R = 6.22 except for the HNC results in frame (b), where a = 4.2 Å (short dashes) and 5.0 Å (long dashes); the corresponding HNC curve of a = 4.6 Å would lie between these curves. The HNC and MC data are taken from ref. 18 and 22.

The results for the effective dielectric permittivity for the divalent electrolytes is shown in Fig. 11b. Again, the HNC approximation provides accurate results only in the dilute region. In contrast to the monovalent case, we see that in the entire interval shown. The under-screening is accordingly accompanied by an augmentation of the effective dielectric permittivity of the electrolyte caused by the strongly nonlinear ion–ion correlations. For 1 : 1 electrolytes we have only in a small part of the very dilute region in the insert to Fig. 10b. Furthermore, we saw in Fig. 3 for the latter case, that quickly approaches ε_r from below when κ_Da is decreased from the Kirkwood cross-over point and then stays nearly constant when κ_Da → 0. In contrast, for the 2 : 2 case in Fig. 11, approaches ε_r from above in a slow manner when κ_Da → 0. Incidentally, we may note that for high values of κ_Da (outside the figure) there is a Kirkwood cross-over point at κ_Da ≈ 1.7 for the 2 : 2 case with a = 4.6 Å and decreases quickly to zero when that point is approached, see ref. 18 and 22.

The nonlinear electrostatic effects in ion–ion correlations are important for any system where the electrostatic coupling is appreciably stronger than the monovalent case we have investigated, that is, for lower temperatures, lower ε_r or higher valencies. This means that terms in higher powers of w^el_ij than the square one in eqn (72) cannot be neglected and the MDE-DH approximation is inadequate in these cases.

From the results from the EDHX approximation we see that apart from dilute systems it is clearly insufficient to have only one decay mode. We need to include at least two decay modes with different decay parameters κ and κ′ as shown in eqn (74) and the corresponding expression for g_+±(r) mentioned earlier. Then we have at least four parameters to determine κ, κ′, and . For the monovalent case we have the fortunate circumstance that κ and κ′ can be quite accurately be determined by solving eqn (33) and that and then can be determined from the application of the local electroneutrality condition and the Stillinger–Lovett second moment condition. This is not the case in general so one will need to apply further conditions that allow four parameters or more to be determined. This will not be pursued in the current paper.

6 The various types of deviations of κ from κ_D for symmetric electrolytes

The deviations of κ from κ_D can be analyzed in terms of the effective charges q^eff since we have the exact relationship eqn (19), , which can be written as [eqn (51)]for the RPM. For the monovalent electrolytes in aqueous solution we have seen that κ > κ_D for most concentrations, so the decay length 1/κ is shorter than the Debye length 1/κ_D and we have over-screening. The relationship to the effective charge given in eqn (78) implies that q^eff is larger than the actual (bare) ionic charge q. This is also apparent from our approximate formula (31) for q^eff.

The fact that q^eff > q for most concentrations in the monovalent case is simply a geometrical effect of the finite ion size. When the ion diameter a is increased from some small value, the surrounding ion cloud is pushed further out, as apparent already in the DH approximation, where eqn (11) for r ≥ a can be written as

with the factor (r − a) in the exponent. The range of the cloud is extended due to this factor when a is increased. This makes the electrostatic potential larger far from the ion, as apparent in eqn (10) for the DH case, so the ion appears to have a larger charge, that is, a larger q^eff. In the DH approximation this applies only for the central ion, but in the general case this applies to all ions because all of them are treated in the same manner, so if q^eff > q we have κ > κ_D.

For cases like divalent (2 : 2) electrolytes where κ < κ_D, there is, instead under-screening, so the decay length 1/κ is longer than the Debye length and q^eff < q. This lowering of q^eff is caused by strong anion–cation correlations; each ion strongly attracts ions of opposite charge causing large values of g₊₋(r) near contact, r = a. The large amount of opposite charge near the ion leads to a reduction of the potential for large r. The ion thus appear to have a smaller charge, that is, a smaller q^eff. The strong (nonlinear) anion–cation correlational attractions hence gives q^eff < q and κ < κ_D. Nothing more than strong anion–cation correlations are needed to give rise to the negative deviations of κ from κ_D in dilute solutions. This is apparent from the EDHX approximation, which shows this behavior, and when the electrostatic coupling is large, this gives, as we have seen, a large negative deviation.

Decay lengths 1/κ that are much longer than the the Debye length 1/κ_D have been experimentally obtained in surface force measurements of divalent (2 : 2) and trivalent (3 : 3) electrolytes in dilute aqueous solutions¹⁰⁴ and monovalent electrolytes in dilute isopropanol solutions¹⁰⁵ by Trefalt and coworkers. Isopropanol has dielectric constant ε_r = 17.9, so the electrostatic coupling for monovalent ions in isopropanol corresponds approximately to divalent ions in water. Their experimental results are shown in Fig. 12 together with theoretical data that will be discussed later.

Fig. 12 The ratio κ/κ_D as function of concentration on a logarithmic scale for monovalent (1 : 1), divalent (2 : 2) and trivalent (3 : 3) electrolytes in dilute aqueous solution and monovalent (1 : 1) electrolytes in dilute isopropanol solution. Experimental data from ref. 104 and 105 are shown as symbols and predictions of the asymptotic formula (81) for low concentrations are shown as curves. The values of the ionic diameter a, which is the single fitting parameter for the asymptote, are shown in the rectangles to the right in the plot, where also the ionic species are displayed.

They analyzed their results in terms of ion pairing, where anion–cation pairs are in equilibrium with dissociated ions, as described by a simple mass action law with equilibrium constant K, where the ionic activity coefficients are set to one. By only including the dissociated (free) ions with density n^free_i in the DH expression (1) for the decay parameter, they obtained an effective κ, defined from , which can be written

for symmetric electrolytes. The paired ions thereby do not contribute since the pair has zero net charge. This expression was successfully fitted by Trefalt and coworkers to the measured values of κ by selecting reasonable values for K. This approach leads to the asymptotic formula¹⁰⁴where a₁ is a negative constant that is proportional to K and a₂ = a₁β_R²/z⁴. This does not agree with the exact formula (73), but can nevertheless be a fair approximation within a finite concentration interval. Eqn (80) says that the deviation (κ^eff_D/κ_D)² from 1 is proportional to the ion concentration in the limit of infinite dilution. As we have seen, this is typical for approximations where nonlinear ion–ion correlations are not correctly included.

Ion pairing in simple electrolytes is an example of strong anion–cation correlations, where the contact value of g₊₋(r) is very large. The pairing is a transient state for the ions involved in the pairing. Fundamentally, all ions should be treated on the same basis irrespectively if they are paired or not, but as an approximation one can model the system such that paired ions are treated as a separate species as done in the pairing theories mentioned in the Introduction. The effective charge q^eff above for symmetric electrolytes, thereby correspond to an average over free ions and ions that are members of some pair,

where q^eff_paired = 0 in the RPM. Since the paired species does not contribute to κ² given by eqn (19), we obtain
In the approach by Trefalt and coworkers, they have accordingly tacitly assumed that the effective charge of the free ions is equal to the bare ionic charge, which gives eqn (79) without the factor q^eff_free/q. This can be a reasonable approximation for a system with high electrostatic coupling provided n^free ≪ n. However, this assumption, together with their additional assumption that the ionic activity coefficients are equal to one, set a limit on the applicability of their approach. The equilibrium condition used for paired and free ions thereby constitutes a simple model for the effects of ion–ion correlations. The reason why the asymptotic expression (80) does not agree with the exact formula (73) is the approximation q^eff_free = q.

Even if an ion pair is a transient state caused by strong anion–cation correlations rather than a kind of separate species, the two ions will act like a dipole in the interactions with the other ions including other transient pairs. There may also be other transient aggregates with more than two ions. Pairs and other transient aggregetes and their effects on the decay lengths and the effective dielectric permittivities of the electrolyte can be handled by the exact DIT formalism.^94,95 It is, for example, found that the strong anion–cation attractive correlations at short separations (transient pairs corresponding to high values of g₊₋(r) near contact) have qualitatively the same effect on the dielectric properties of the electrolyte as an actual pair. As we saw in Fig. 11, we have for 2 : 2 electrolytes of low to moderate concentrations, so there is an augmented dielectric permittivity of the electrolyte as expected from such a mechanism.

In order to investigate the dilute region further, let us extract the leading asymptotic terms for κ²/κ_D² and in the limit τ = κ_Da → 0 from eqn (76) and (77) of the EDHX approximation. We will thereby include terms up to order τ², that is, up to terms that are proportional to concentration. Since this approximation is designed to give accurate results for dilute electrolytes, the derived asymptotic formulas should be correct in the limit of zero concentration. As shown in Appendix D we obtain

andwhere the constant , is Euler's constant and s(x) = x/(5!·2) + x²/(7!·4) + … is defined in eqn (136) of Appendix D.

In the rhs of eqn (81), the second term, which is negative for small κ_Da, dominates when κ_Da → 0 due to the logarithm and decays to zero in this limit. Eqn (81) agrees with the exact limiting law (73) since β_Rτ = z²Λ and ln τ = ln(Λ) + ln(z²/β_R) ∼ ln(Λ) when the concentration and hence Λ and τ go to zero.¶¶¶

For divalent electrolytes in aqueous solution, the predictions of the asymptotic formulas (81) and (82) are plotted as red dotted curves in Fig. 11a and b, respectively. Despite that these formulas were derived from the EDHX approximation in the limit zero concentrations, each asymptote continues not far from the simulated results in a quite wide concentration interval and give therefore surprisingly reasonable estimates of the latter. This means that the asymptotes can be useful also outside the region of very high dilution, at least for the purpose of approximate assessments of ion–ion correlation effects in κ²/κ_D² and at high electrostatic coupling where the MDE-DH approximation is inadequate.

As a test of such an assessment, the asymptote from formula (81) has been compared with the experimental data by Trefalt an coworkers^104,105 mentioned earlier.|||||| Data are available for CsCl, NaCl, CuSO₄, MgSO₄ and LaFe(CN)₆ solutions in water and tetrabutylammonium bromide (TBAB) and LiCl solutions in isopropanol. Note that for a given T, ε_r and valency z of the symmetric electrolyte, the predicted asymptote as function of concentration contains only one parameter, namely the ionic diameter a. In the asymptotic formula, the diameter will correspond to an average value since the anions and cations in the experiments have somewhat different sizes.

The results of a comparison with the experimental data with a as the single fitting parameter is presented in Fig. 12, whereby all cases apart from monovalent ions in water (CsCl and NaCl) were fitted; the experimental uncertainties for the latter were too large to make a fit meaningful. For the other cases, the fits to the low-concentration part of the experimental data are quite good overall considering the experimental uncertainties. The values of the ionic diameters obtained are reasonable, but as we will see below, they are probably somewhat too large. Note that ionic sizes used in primitive model calculations fitted to experimental data are usually different from bare diameters and often include effects of hydration, in particular for multivalent ions. As regards the LiCl solutions in isopropanol there are to few experimental data points available in order to say something definite about the fit in this case. The first data point to the left in the plot for this system and for TBAB in isopropanol are outliers, which demonstrate the experimental difficulties in this kind of measurements.

We can conclude that strong anion–cation correlations constitute the explanation for the experimentally observed over-screening in dilute solutions at room temperature when z is high or ε_r is low, whereby the decay lengths are considerably longer than the Debye length. The same is true for thin ionic plasmas at moderately high temperatures. There is no need to assume actual ion pairing, although such pairing can be used as a simple way to model such correlations as mentioned in the Introduction and discussed above. The asymptotic formula (81) can accordingly be used as a simple means to approximately assess nonlinear ion–ion correlations effects in dilute solutions. The fact that the asymptote nearly quantitatively describes the variation of the decay length with electrostatic coupling strengths (varying z, ε_r or T) is thereby a great advantage.

Finally, it is of interest to investigate the asymptotes from eqn (81) and (82) in more detail in the dilute region. The ratios κ²/κ_D² and for divalent and monovalent symmetric electrolytes in water are plotted in Fig. 13 as functions of κ_Da for low values of the latter. Fig. 13a shows a magnified view of Fig. 11a in the dilute region of the 2 : 2 case and we see that the asymptote joins the other curves when κ_Da → 0 as it should. It thereby joins the HNC curve somewhat faster than the curve from the EDHX approximation in this limit.

Fig. 13 The ratios κ²/κ_D² and as functions of κ_Da for divalent (2 : 2) and monovalent (1 : 1) electrolytes in dilute aqueous solution compared with the asymptotes (red curves). (a) The same curves as in Fig. 11a for the divalent case shown the dilute region at large magnification. The red dotted curve is the asymptote for κ²/κ_D² from eqn (81). The ion diameter is a = 4.6 Å and β_R = 6.22. (b) An illustration of the ion size dependence of κ²/κ_D² for divalent electrolytes; from bottom to top a = 4.2, 4.6 and 5.0 Å. Black curves are from HNC calculations and red curves from the asymptotic formula (81). (c) The red dotted curve shows the asymptote for κ²/κ_D² from eqn (81) for monovalent electrolytes with β_R = 1.55. The full, dashed-dotted and dashed curves (MDE-DH, EDHX and HNC results respectively) are the same as the corresponding curves in Fig. 10a. (d) for divalent electrolytes (black curves and red dots, β_R = 6.22) and monovalent electrolytes (blue curves and red mini-dashes, β_R = 1.55) in aqueous solution. The red curves show from the asymptotic expression (82). The curves for the divalent case are the same as in Fig. 11b shown the dilute region at large magnification and the blue curves are the same as the corresponding ones in Fig. 10b. The HNC and MC data are taken from ref. 18 and 22.

Fig. 13b shows how the ion diameter affects both the HNC curves and the asymptotes when a varies between 4.2 to 5.0 Å for 2 : 2 electrolytes. The middle case, a = 4.6 Å, is the same as in frame (a) of the figure. We can see that the asymptote for a certain a value, say a = 5.0 Å, lies close to the HNC curve for a slightly smaller a value, say 4.6 Å, in a whole κ_Da interval. This implies that the fits made in Fig. 12 for the asymptotes most likely give somewhat too large values of a.

In Fig. 13c it is seen that the asymptote for the monovalent case accurately describes the slight dip of κ²/κ_D² below 1 for small κ_Da. When κ_Da is increased, the asymptote approaches and then crosses the full curve, that is, the decay parameter from the MDE-DH approximation, eqn (33), which goes like like κ²/κ_D² ∼ 1 + τ²/2 for small τ. We may note that eqn (81) contains a τ²/2 contribution, namely from the term 1/2 in the coefficient of the τ² term.

The ratio for both divalent and monovalent electrolytes is plotted in Fig. 13d. For the 2 : 2 case the asymptote lies quite closely to the the accurate results in a wide interval, as we saw previously in Fig. 11. The asymptote in the 1 : 1 case joins the accurate curves from above, , in the dilute region, where these latter curves lie slightly above 1, as seen more clearly in the inset to Fig. 10b.

7 Summary and concluding remarks

This work is based on the following observations. In electrolytes there are several decay modes for the electrostatic interactions, each having a decay length and a magnitude. In the DH and PB approximations there is, however, only one mode that gives a screened Coulomb potential equal to ϕ_Coul*(r) = e^−κ_Dr/(4πε₀ε_rr) with a decay length equal to the Debye length 1/κ_D and magnitude proportional to 1/ε_r. The mean electrostatic potential ψ_i(r) due to an ion of species i is in the DH approximation obtained for r > a by multiplying ϕ_Coul*(r) by the effective charge q^eff_i given by eqn (12), so we have ψ_i(r) = q^eff_ie^−κ_Dr/(4πε₀ε_rr). In the PB approximation ψ_i(r) decays for large r as the rhs of this expression.

A correct treatment of electrolytes provides a screened Coulomb potential ϕ_Coul*(r) with multiple decay modes, where the modes give contributions like (cf.eqn (23))

with magnitudes inversely proportional to and , respectively, which are effective dielectric permittivities of the electrolyte for the respective mode. They depend on the dielectric response of the electrolyte. We have here only included the two main modes with the largest decay lengths. The mean electrostatic potential ψ_i(r) has the following principal contributions from these modes (cf.eqn (25))where q^eff_i and are the effective charges of the i-ion, whereby each mode has its own value of this charge. The decay parameter κ can be expressed in terms of q^eff_j as [eqn (4) and (24)] and κ′ is given by the analogous expression with . The electrostatic part of the potential of mean force w_ij(r) for species i and j has contributions with yet another factor of effective charge, namely (cf.eqn (21))For elevated ionic densities these functions turn into exponentially damped oscillatory ones, whereby κ and κ′ are complex-valued, κ = κ_ℜ + iκ_ℑ and κ′ = κ_ℜ − iκ_ℑ. The contribution in eqn (85) from the second mode then is the complex conjugate of that from the first. The sum of the two modes then decays like e^−κ_ℜrcos(κ_ℑr + α), where α is a phase.

In the Multiple-Decay Extended Debye–Hückel, MDE-DH, approximation developed in the current work, one assumes as an approximation that solely the first two decay modes contribute to the electrostatics in electrolytes, so w^el_ij is given by the sum of the contributions in eqn (85). Furthermore one takes w_ij(r) = w^el_ij(r) + w^core_ij(r), where w^core_ij is the contribution from the ionic core–core correlations that is approximated by w^hs from a hard sphere fluid. The term w^core_ij expresses a further decay mode for the interactions. We have

and only the first few terms of the expansion in βw^el_ij are used as a further approximation, so the MDE-DH approximation is limited to sufficiently weak electrostatic coupling where βw^el_ij is small. The decay parameters κ and κ′ used are solutions to the equation κ²/κ_D² = e^κa/(1 + aκ) [see eqn (33)] originally derived in the LT and MDH approximations.^32,33 These solutions are presented in Fig. 2. The effective charges are obtained by inserting the respective κ value into q^eff_i = q_ie^κa/(1 + κa) [eqn (31)]. By applying the necessary local electroneutrality and the Stillinger–Lovett second moment conditions, explicit expressions for and are derived and the resulting theory does not contain any fitting parameters.

For symmetric electrolytes with ions of diameter a, which is the main case treated in the current work, ψ_i(r) is proportional to the ionic charge q_i, so the MDE-DH theory is a linearized approximation for the electrostatic response like, for example, the MSA, GMSA, LMPB, MDH and GDH theories. In contrast to these latter approximations, the MDE-DH theory contains the square term [βw^el_ij(r)]²/2, which contributes to the pair distributions but not to ψ_i (see Section 4.2). This term guarantees that g_ij(r) ≥ 0 for all r and leads to a good thermodynamic consistency.

The approximation comes in two versions, the Simple and the Complete MDE-DH approximations (the differences between the two are explained in Section 4.2, eqn (55) ff.). The simple version is applicable when w^core_ij is not too large, for example systems with small ions. Its results are in very good agreement with simulations and Hypernetted Chain calculations for all concentrations investigated, including high ones, provided the electrostatic coupling is not too high for instance monovalent electrolytes in water at room temperature. All numerical results in Section 4 apart from Fig. 7, where the ions are large, are calculated in this approximation. Furthermore, in the Simple MDE-DH approximation, there are simple analytical expressions for and [eqn (46) and (47)] and for the thermodynamic properties: the activity coefficient γ_± is given by ln γ_± = ln γ^el_± + ln γ^core_± with ln γ^el_± from eqn (64) and ln γ^core_± from eqn (60) and the osmotic coefficient ϕ = 1 + ϕ^el + ϕ^cont with ϕ^el from eqn (69) and ϕ^cont from eqn (68) [see also eqn (124)]. These expressions are valid also for complex-valued κ and κ′, but they are reformulated in Appendix C into equivalent formulas that are more easy to use for that case.

The Complete MDE-DH approximation is used to obtain the results in Fig. 7, which shows the pair distribution functions for a dense system with large monovalent ions. These results are in excellent agreement with simulations, including quite intricate details. It is, however, sufficient to use the formulas of Simple MDE-DH approximation to calculate the thermodynamic quantities for this case; the difference between the two versions of the approximation for these quantities is less than 1%. Thus, in practical applications, the simple version of the approximation is in many cases sufficient unless structural entities are required. This is agreement with the observation discussed in the Introduction that it is simpler to calculate thermodynamical entities than structural ones.

In the MDE-DH approximations and the other approximations with linearized electrostatic response mentioned above, we have κ ≥ κ_D for all concentration where κ is real, while the exact limiting law (5) says that there is under-screening, κ < κ_D, for sufficiently dilute solutions of symmetric electrolytes. This under-screening is caused by nonlinear electrostatic response of the electrolyte that is not included in the MDE-DH approximations. In order to investigate this matter, the nonlinear approximation given by the first equality in eqn (86) is used in the low density limit, where contributions from w^core_ij and the second decay mode in w^el_ij go to zero and therefore can be neglected. The approximation where these contributions to w_ij are not included is called the Extended DHX, EDHX, approximation since it is an extension of the exponential Debye–Hückel, DHX, approximation (cf. Section 2).

The EDHX approximation is used to calculate κ and for symmetric electrolytes and it is found that κ < κ_D for sufficiently dilute solutions in agreement with the exact limiting law. Furthermore, it is found that in the dilute region, while the linear approximations predict that . The EDHX results are thereby in agreement with the HNC approximation, which is considerably more complex. In the case of monovalent electrolytes in water at room temperature, it is found that this qualitatively incorrect behavior of the linear approximations has negligible consequences in practice because they occur for so low concentrations that the deviations of κ/κ_D and from 1 are miniscule.

For multivalent ions in water and monovalent ions in solvents with considerable lower ε_r, the situation is completely different and under-screening is very substantial for dilute solutions of symmetric electrolytes, as found by MC simulations, HNC calculations and experiments. Likewise, is considerably larger than ε_r for the dilute solutions due to the nonlinear ion–ion correlation effects. An approximate limiting law (81) for symmetric electrolytes that extends the exact law by including higher order terms is derived in the EDHX approximation and a corresponding law for , eqn (82), is also derived. These approximate laws contain the ion size dependence of κ and for dilute solutions. Despite that they are strictly valid only for very low concentrations, they can be used as a reasonable approximation in a wider concentration interval to describe the effect of nonlinear ion–ion correlation effects in dilute electrolytes. The approximate limiting law for κ/κ_D is used to analyze the experimental results by Trefalt and coworkers^104,105 for this ratio in aqueous multivalent electrolyte solutions and monovalent electrolytes dissolved in isopropanol. The law is successfully fitted to the experimental results for low concentrations whereby the ion diameter is the sole fitting parameter. This gives important insights into the reasons behind the deviations of κ from κ_D and mechanisms of the latter is discussed in some detail Section 6.

In order to give a further perspective on the decay modes and the associated entities in electrolytes, it is worthwhile to mention the following facts regarding more sophisticated models for electrolytes than the primitive model. Recall that for a classical plasma of spherical ions, the decay parameters κ, κ′ etc. satisfy the exact eqn (24) with ε_r = 1. For an electrolyte solution with molecular solvent, that is, in presence discrete polar molecules, we still have ε_r = 1 because the particles are in vacuum, but eqn (24) is not valid as it stands because it is, in fact, restricted to systems where all particles are spherically symmetric. The solvent molecules and any ion species that is not spherical give rise to additional terms in this equation that depend on the orientational degrees of freedom of these particles. The equation thereby is replaced by the following exact equation⁹⁸

where the first sum is over all spherical species and the second sum over the nonspherical ones. Here, Q_j is a particle-orientation dependent entity that originates from the internal charge distribution of a j-particle and Q^eff_j is also orientation dependent but originates from both the particle and its surrounding distribution of charge. The entity Q^eff_j plays a similar role as a prefactor for ψ_i for a nonspherical particle as q^eff_j does for a spherical one (for definitions and further details see ref. 98). The average 〈·〉 in eqn (87) is taken over all orientations of the j-particle. The last term contains, for example, effects of solvent-solvent correlations (both orientational and spatial ones). The solvent-ion correlations also contribute to this term and such correlations also influence the values of q^eff_j in the presence of solvent because oriented solvent molecules around an ion affects the charge distribution there.

Eqn (87) has a rather complicated form. It is, however, mathematically equivalent to a much simpler, exact equation^95,98

where q_j* is a renormalized ionic charge**** that is different from q^eff_j. , called the dielectric factor, is a kind of renormalized dielectric permittivity, which is not the same as . Both and can be defined in terms of static dielectric response function of the electrolyte.^95,98 Note that is a function of κ, while q_j* is independent of κ and has the same value for all decay modes. The decay parameters κ, κ′ etc. are solutions to this equation and, of course, also solutions to the equivalent eqn (87).†††† contains the major effects of solvent-solvent, solvent-ion and ion–ion correlations that are important for the decay mode with decay parameter κ. The value of q_j* is, of course, also affected by these correlations, but since it does not depend on κ, it is not specific for any mode. The general eqn (88) has a form that is very similar to the definition (1) of the Debye parameter.

Thus there are two different entities and that take the role that ε_r has in the PB and DH approximations. appears in the eqn (88) for κ and appears, as we have seen, in the magnitude of, for example, the screened Coulomb potential ϕ_Coul*, the mean electrostatic potential ψ_j and the charge density ρ_j. Each mode has its own values for these entities.

In some PB, DH and other primitive model approaches one replaces ε_r by an ion-concentration dependent entity ε^electrolyte_r, which is obtained from macroscopic measurements of the dynamic dielectric response of bulk electrolytes at finite but small frequency, whereby a zero frequency macroscopic response is extracted. This entity contains contributions from the orientational and translational degrees of freedom of the solvent molecules and from ion-solvent correlations that are not included in the primitive model. It is, however, not the same as or . When some mode has a small κ it may, however, be a reasonable approximation to set and equal to ε^electrolyte_r for that particular mode, but not for the other modes. When the ion concentration goes to zero, and for the mode with longest decay length, which goes to the Debye length in this limit.

For a system with spherical particles, eqn (88) is mathematically equivalent to eqn (24). In the primitive model both and solely contains contributions from the ionic degrees of freedom.

Conflicts of interest

There are no conflicts to declare.

Appendices

A Pair distribution functions for RPM, theoretical considerations

In this Appendix we will fill in the details regarding the pair distribution functions for RPM electrolytes in the MDE-DH approximation. As customary, we define the functionswhereby g_+±(r) = g_S(r)±g_D(r). Eqn (53) implies that g_S(r) = cosh[βw^el(r)]g^core(r) and g_D(r) = − sinh[βw^el(r)]g^core(r), so by inserting eqn (54) we haveandNote that

ρ(r) = qn2g_D(r) = qn_totg_D(r) (92)

and that the density of ions around an ion is given by n_totg_S(r). The Debye parameter is given by κ_D² = βn_totq²/ε_rε₀ and we have κ² = βn_totqq^eff/ε_rε₀.

When the electrostatic coupling is sufficiently weak, we can use the approximations cosh(x) ≈ 1 + x²/2 and sinh(x) ≈ x in these expressions and obtain as a good approximation

g_D(r) = −βw^el(r)g^core(r). (94)

Explicitly, this isandThis is the basis of the Complete MDE-DH approximation for the RPM. The radial charge distribution isin this approximation.

In the Simple MDE-DH approximation, where the cross-terms between h^core and w^el are neglected, we obtain for r ≥ a

g_D(r) = −βw^el(r) (99)

or explicitlyandfor r ≥ a. By inserting this g_D(r) into eqn (92) and using κ² = βn_totqq^eff/ε_rε₀, we obtain the radial charge density ρ(r) in Ansatz (43). The presence of the square term in g_S(r) in this approximation means that there is a contribution with the decay length 1/2κ in this function, which is in semi-quantitative agreement with the exact asymptotic analysis in ref. 101 (cf. footnote †† in Section 2).

When one uses the Complete MDE-DH approximation, the expressions for and are different from eqn (44)–(49). As shown in Appendix B, the local electroneutrality and the Stillinger–Lovett second moment conditions yield

whereand the corresponding expressions for M₁′ and M₃′ with κ′ inserted. When g^core(r) ≈ 1 one obtains M₁ ≈ 1 and M₃ ≈ 6e^−κa exp₃(κa) and likewise for M₁′ and M₃′, whereby the result in eqn (46) and (47) is obtained as a good approximation.

In the case of complex-valued decay parameters κ and , we have and , whereby eqn (102) and (eqn 103) yield and . Explicit expressions for and ϑ are given in Appendix B [eqn (109) and (110)]. From eqn (95) and (96) it follows that

andwhere η is defined from q^eff = |q^eff|e^−iη and the phase is α = 2η − ϑ = 4θ − ϑ since η = 2θ in the RPM, which follows from eqn (51).

B Equations for and

In this Appendix we will derive eqn (102) and (103) of Appendix A for and of RPM electrolytes in the Complete MDE-DH approximation when κ and κ′ are real-valued. We will also obtain the corresponding equations for the complex-valued case.

We have from eqn (97)

Let us introduce the unknowns
The local electroneutrality condition yields
and the Stillinger–Lovett second moment condition (34), which can be written as in the RPM, yields

This set of equations can be expressed as

where the matrix elements M₁ and M₃ are defined in eqn (104) and (105) in Appendix A and M₁′ and M₃′ are defined analogously with κ′ inserted. The solution to eqn (108) is Y = (M₃′ − 6M₁′)/D and Y′ = −(M₃ − 6M₁)/D, where D = M₁M₃′ − M₁′M₃ is the determinant. This solution can be written as in eqn (102) and (eqn 103) in Appendix A.

In the case of complex-valued decay parameters κ = κ_ℜ + iκ_ℑ = |κ|e^−iθ and we have , and . In this case we can write and , where ℑ(·) means the imaginary part. Furthermore we have

and
where G^cos_ν and −G^sin_ν are defined as the real and the complex parts of the integrals, respectively. We hence obtain
where we have used −i = e^i3π/2. By extracting |Y| and arg(Y) we obtainandwhen G^cos₁G^sin₃ − G^sin₁G^cos₃ ≥ 0 (otherwise there is a further term −π in ϑ).

C Theoretical considerations for U^ex, γ_± and ϕ

Here we will present the theoretical details regarding the thermodynamic entities. The electrostatic contribution μ^ex,el_i to the excess chemical potential is, as we have seen, the free energy change (= the reversible work) when the charge of a sphere inserted in the electrolyte is changed from 0 to q_i. It can be calculated via a coupling parameter integration
where u^el_ij(r) = q_iq_j/(4πε_rε₀r) is the electrostatic pair potential, ξ is the coupling parameter 1 ≤ ξ ≤ 1, g_ij(r;ξ) is the pair distribution function around a partially charged ion with charge ξq_i placed at the origin and ρ_i(r;ξ) is the charge density around this ion. All other ions are fully charged throughout. The corresponding contribution to the mean activity coefficient γ_± isFor approximations where the charge density ρ_i(r) is proportional to the ionic charge q_i, like the DH, MDH and MDE-DH approximations, we have ρ_i(r;ξ) = ξρ_i(r). Then, the integration over ξ contributes with a factor so we havefor such linear approximations. In the DH and MDH approximations this gives ln γ^el_± in eqn (62) and (63), respectively.

In the Simple MDE-DH approximation based on the Ansatz (43), we take ρ_i(r) from the second line in this equation and obtain after simplification

where we have made use of eqn (19) and its primed counterpart eqn (22). Using eqn (33), we can alternatively write this expression as eqn (64). When the charge density is oscillatory and κ and κ′ are complex-valued, these expressions for ln γ^el_± can be used as they stand, but they can alternatively be written in a more convenient manner aswhich can be explicitly expressed in terms of κ_ℜ and κ_ℑ as
Note that tan ϑ can be obtained in terms of κ_ℜ and κ_ℑ from eqn (49).

The average internal energy U of the system is

The second term on the rhs can be written aswhere we have used N_i/N_tot = n_i/n_tot. By comparing with eqn (112), which is valid for approximations where the charge density ρ_i(r) is proportional to the ionic charge q_i, we see that βU^ex/N_tot = ln γ^el_± in such cases. Consequently, in the Simple MDE-DH approximation βU^ex/N_tot is equal to the rhs of (eqn 113), so (eqn 65) follows when (eqn 33) is used. When κ and κ′ are complex, βU^ex/N_tot can be calculated from rhs of eqn (114).

Next, we treat the osmotic coefficient. The excess pressure of the electrolyte consists of two parts P^ex = P^cont + P^el, where P^cont is the contact pressure evaluated at the hard core surfaces of the ions

andThe excess osmotic coefficient, ϕ^ex = ϕ^cont + ϕ^el is accordingly given bywhere the first term is ϕ^cont = P^cont/P^ideal and the second is ϕ^el = P^el/P^ideal.

In the Simple MDE-DH approximation based on the Ansatz (43) we obtain

which can be written as eqn (69). For the case of complex-valued κ and κ′, this expression is still valid, but it can be written analogously to eqn (114) [the prefactor for ϕ^el has 24π instead 8π in the denominator].

As mentioned in Section 4.3, the DH and MDH approximations show very poor consistency for ln γ_± obtained from eqn (112) and that obtained from the osmotic coefficient viaeqn (71). The reason for this lies in the handling of the electrostatic contributions. Both approximations treat the core contribution consistently, because we have

which yields the contribution 4πa³n_tot/3 when inserted into eqn (71), that is, the same as ln γ^core_± given in eqn (59). Eqn (120) follows from the fact that for r ≥ a in the DH and MDH approximations, so we have . The inconsistency originates from the fact that ϕ^el inserted into eqn (71) does not give the same electrostatic contribution as the one obtained from eqn (111). The latter is, in fact, a general feature of linear approximations, where the charge density ρ_i(r) is proportional to q_i. This can be realized as follows.

For these kinds of approximations, ln γ^el_± is calculated from eqn (112) and a comparison with the second term in eqn (118) shows that the integral is the same in both equations and that the prefactors differ only by a factor of 1/3, which implies that ln γ^el_± = 3ϕ^el. This means that in eqn (71) one obtains one third of ln γ^el_± from the first term on the rhs and that the second term must yield two thirds of ln γ^el_± in order to have consistency. Thus the integral in this second term must be equal to one third of ln γ^el_±, that is, equal to ϕ^el. Consistency is achieved provided the function ϕ^ex(n_tot) satisfies , where we have made the substitutions n_tot′ = t² and n_tot = s². This implies that ϕ^el(s²) = Ks, that is, , where K is a constant. This is satisfied for the Debye–Hückel limiting law (2) for ϕ^ex = ϕ − 1 at infinite dilution, but not otherwise. Since the limiting law is an exact result of statistical mechanics, it must obey the consistency requirement, but no other linear approximation can be consistent in the sense that ϕ^el inserted into eqn (71) yields ln γ^el_± (except in the limit of infinite dilution).

Since the charge density ρ_i(r) is proportional to q_i in the MDE-DH approximations, they share the same feature, that is, ln γ^el_± is not obtained from ϕ^elvia the integration in eqn (71). However, ϕ^el is not the only contribution to ϕ^ex where electrostatics matter in these approximations; there is a further contribution due to the quadratic terms in w^el in eqn (55) and (58) in the complete and simple versions of this approximation, respectively. We will start by considering the Simple MDE-DH approximation and see that there is an electrostatic part in ϕ^contact due to the presence of the quadratic term [βw^el(r)]/2 in eqn (58) and hence in eqn (98). This is an important contribution that is necessary for the good agreement between the MDE-DH and MC data for ϕ in Fig. 8 and 9 and for the good consistency between the two different manners to calculate ln γ_±.

From eqn (118) we obtain for the RPM

where g_S and g_D are defined in eqn (89) and we have used (eqn 92). It follows from eqn (98) that in the Simple MDE-DH approximation we havewhere ϕ^core = 2πa³n_totg^core(a)/3. By inserting w^el(r) from eqn (54) into eqn (123) we obtainwith ϕ^core given by eqn (70). In the Simple MDE-DH approximation, we accordingly have ϕ = 1 + ϕ^cont + ϕ^el with ϕ^cont from eqn (124) and ϕ^el from eqn (119).

In the Complete MDE-DH approximation it follows from eqn (93) and (121) that

and we have from eqn (94) and (122)When w^el(r) from eqn (54) is inserted, ϕ^el can be written in terms of M₁ and M₁′, which are integrals involving g^core defined in eqn (104) in Appendix A. Furthermore, in this approximationwhich contains the same integral, and ln γ^core_± is given by eqn (60).

D Asymptotic formulas at high dilution

In this appendix we consider κ²/κ_D² and in very dilute solutions of symmetric electrolytes of ions with diameter a and derive approximate asymptotic formulas for these entities when τ = κ_Da → 0. These formulas will be extracted from eqn (76) and (77), which can be written aswhere X = κ/κ_D and . We will extract all contributions to order τ² or less, that is, up to and including contributions that are proportional to the total ionic density n_tot. We do this by expanding sinh x = x + x³/3! + x⁵/5! + … and expressing the integrals involving [e^−XτR/R]^ν in terms of known functions, namely
and
where is the incomplete Gamma function and where we have defined fromwhere E₁(x) is the exponential integral . The function for m > 5 has the limitbut for m ≤ 5 it diverges when x → 0. For m = 5 the divergence is logarithmic since when x → 0, where is Euler's constant.

Using these results in eqn (128), which is the local electroneutrality condition, we obtain

and likewise eqn (129), which is the Stillinger–Lovett condition, can be writtenSince X = κ/κ_D → 1 and when τ → 0, we have X = 1 + α_X(τ) and Y = 1 + α_Y(τ), where the as yet unknown functions α_X(τ) → 0 and α_Y(τ) → 0 when τ → 0. Hence

X^l ∼ 1 + lα_X(τ) and Y^l ∼ 1 + lα_Y(τ) when τ → 0. (134)

We start with the Stillinger–Lovett condition (133). We first note that when τ → 0. This implies that the term for ν = 5 in eqn (133) decays like −τ⁴ ln τ in this limit. This factor decays faster than τ² so the ν = 5 term can be skipped. Using eqn (131) we can see that the terms for ν > 5 decay like τ⁴, so they can also be skipped. Thus there remains only the terms for ν = 1 and 3.

The ν = 1 term in the left-hand side (lhs) of eqn (133) is

where we have used eqn (130), and the ν = 3 term is
We now insert these terms into eqn (133), whereby we skip contributions that decay faster to zero than τ² in order to extract the asymptote. In the ν = 3 term, the factor e^−3Xτ(1 + 3Xτ) decays when τ → 0 like

e^−3Xτ(1 + 3Xτ) ∼ 1 − (3Xτ)²/2  ∼  1 − (3τ)²/2 − α_X(τ)(3τ)² ∼ 1 − (3τ)²/2,

where we have used eqn (134). Furthermore, X¹⁰ ∼ 1 + 10α_X(τ) and Y³ ∼ 1 + 3α_Y(τ). This means that τ²X¹⁰Y³e^−3Xτ(1 + 3Xτ) decays like τ² and we can conclude that the ν = 3 term decays like β_R²τ²/324. By inserting this result and the ν = 1 term into eqn (133) and rearranging we obtain
Now, e^−Xτ exp₃(Xτ) ∼ 1 − (Xτ)⁴/24 ∼ 1 − τ⁴/24, so [e^−Xτ exp₃(Xτ)]⁻¹ ∼ 1 + τ⁴/24, so we can finally conclude thatto the order τ².

We now proceed with the local electroneutrality condition (132). The ν = 1 term in the lhs of this equation is

and the ν = 3 term is
Analogously to the previous arguments we find that τ²X¹²Y³ decays like τ², so the entire ν = 3 term decays like up
to the order τ². There are also contributions of the order τ² from the terms with ν ≥ 5. The factor τ²X^4νY^ν in these terms decays like τ² and by using eqn (131) we see that the terms decay when τ → 0 like
The leading contribution from all terms with ν ≥ 5 is accordingly equal to β_R²τ²s(β_R²), where we have definedwhere we have made the substitution ν = 2l + 3. By gathering these results and inserting them together with the ν = 1 term into eqn (132), we obtain after rearrangement
Now, from eqn (135) it follows that 1/Y decays like 1 + β_R²τ²/324 and we have [e^−Xτ(1 + Xτ)]⁻¹ ∼ 1 + (Xτ)²/2 ∼ 1 + τ²/2. Thus we can finally conclude thatto the order τ². Eqn (135) and (137) constitute the wanted result, which can be written as eqn (82) and (81), respectively, in the main text.

Acknowledgements

I want to thank Andrés González de Castilla for giving me the inspiration to do this work, for interesting discussions on practical applications of the theory of electrolytes and for helpful comments on the manuscript of this paper. Gregor Trefalt is cordially thanked for sending me preprints and raw data for the measurements of the decay lengths in dilute electrolytes and for useful comments on the manuscript. I also want to thank Sture Nordholm for his constructive remarks on the paper.

References

1. M. Su and Y. Wang, Commun. Theor. Phys., 2020, 72, 067601.
2. L. L. Lee, Molecular Thermodynamics of Electrolyte Solutions, World Scientific Publishing, Singapore, 2008.
3. L. M. Varela, M. García and V. Mosquera, Phys. Rep., 2003, 382, 1.
4. L. Lue, in Variational Methods in Molecular Modeling, ed. J. Wu, Springer, Singapore, 2017, p. 137.
5. B. Maribo-Mogensen, G. M. Kontogeorgis and K. Thomsen, Ind. Eng. Chem. Res., 2012, 51, 5353.
6. G. M. Kontogeorgis, B. Maribo-Mogensen and K. Thomsen, Fluid Phase Equilib., 2018, 462, 130.
7. P. Debye and E. Hückel, Phys. Z., 1923, 24, 185.
8. K. S. Pitzer, J. Phys. Chem., 1973, 77, 268.
9. K. S. Pitzer and G. Mayorga, J. Phys. Chem., 1973, 77, 2300.
10. K. S. Pitzer, Acc. Chem. Res., 1977, 10, 371.
11. J. C. Rasiah, D. N. Card and J. P. Valleau, J. Chem. Phys., 1972, 56, 248.
12. G. Hummer and D. M. Soumpasis, J. Chem. Phys., 1993, 98, 581.
13. J. G. Kirkwood, Chem. Rev., 1936, 19, 275.
14. J. G. Kirkwood and J. C. Poirier, J. Phys. Chem., 1954, 58, 591.
15. G. A. Martynov, Sov. Phys. - Usp., 1967, 10, 171.
16. C. W. Outhwaite in Statistical Mechanics: A Specialist Periodical Report, ed. K. Singer, The Chemical Society, London, 1975, vol. 2, p. 188.
17. P. Keblinski, J. Eggebrecht, D. Wolf and S. R. Phillpot, J. Chem. Phys., 2000, 113, 282.
18. J. Ulander and R. Kjellander, J. Chem. Phys., 2001, 114, 4893.
19. F. Coupette, A. A. Lee and A. Härtel, Phys. Rev. Lett., 2018, 121, 075501.
20. Supplementary material for ref. 19 at http://link.aps.org/supplemental/10.1103/PhysRevLett.121.075501.
21. J. Ennis, PhD thesis, Australian National University, 1993.
22. J. Ennis, R. Kjellander and D. J. Mitchell, J. Chem. Phys., 1995, 102, 975.
23. J. Ulander and R. Kjellander, J. Chem. Phys., 1998, 109, 9508.
24. R. J. F. Leote de Carvalho and R. Evans, Mol. Phys., 1994, 83, 619.
25. E. Waismann and J. L. Lebowitz, J. Chem. Phys., 1972, 56, 3086.
26. E. Waismann and J. L. Lebowitz, J. Chem. Phys., 1972, 56, 3093.
27. C. W. Outhwaite and V. C. L. Hutson, Mol. Phys., 1975, 29, 1521.
28. L. Blum. in Theoretical Chemistry, Advances and Perspectives, ed. H. Eyring and D. Henderson, Academic Press, New York, 1980, vol. 5, p. 1.
29. L. Blum, J. N. Herrera and J. F. Rojas-Rodríguez, Rev. Mex. Fis., 1995, 41, 618.
30. C. W. Outhwaite, Chem. Phys. Lett., 1970, 5, 77.
31. C. W. Outhwaite and L. B. Bhuiyan, Condens. Matter Phys., 2019, 22, 23801.
32. R. Kjellander, J. Phys. Chem., 1995, 99, 10392.
33. D. G. Hall, Z. Phys. Chem., 1991, 174, 89.
34. B. P. Lee and M. E. Fisher, Phys. Rev. Lett., 1996, 76, 2906.
35. B. P. Lee and M. E. Fisher, Europhys. Lett., 1997, 39, 611.
36. L. M. Varela, M. Perez-Rodriguez, M. Garcia, F. Sarmiento and V. Mosquera, J. Chem. Phys., 1998, 109, 1930.
37. L. M. Varela, M. Perez-Rodriguez, M. García and V. Mosquera, J. Chem. Phys., 2000, 113, 292.
38. L. M. Varela, J. M. Ruso, M. García and V. Mosquera, J. Chem. Phys., 2000, 113, 10174.
39. L. M. Varela, M. García and V. Mosquera, Physica A, 2003, 323, 75.
40. M. Ding, Y. Liang, B.-S. Lu and X. Xing, J. Stat. Phys., 2016, 165, 970.
41. Y. Avni, R. M. Adar and D. Andelman, Phys. Rev. E, 2020, 101, 010601.
42. R. Triolo, J. R. Grigera and L. Blum, J. Chem. Phys., 1976, 80, 1858.
43. L. Blum and J. S. Høye, J. Chem. Phys., 1977, 81, 1311.
44. T. Sun, J.-L. Lénard and A. S. Teja, J. Phys. Chem., 1994, 98, 6870.
45. B. Larsen, G. Stell and K. C. Wu, J. Chem. Phys., 1977, 67, 530.
46. M. C. Abramo, C. Caccamo and G. Pizzimenti, J. Chem. Phys., 1983, 78, 357.
47. G. Stell and S. F. Sun, J. Chem. Phys., 1975, 63, 5333.
48. T. Xiao and X. Song, J. Chem. Phys., 2011, 135, 104104.
49. T. Xiao, Electrochim. Acta, 2015, 178, 101.
50. T. Xiao and X. Song, J. Chem. Phys., 2017, 146, 124118.
51. J. Rasaiah, Chem. Phys. Lett., 1970, 7, 260.
52. J. C. Rasiah, J. Chem. Phys., 1972, 56, 3071.
53. B. Larsen, J. Chem. Phys., 1978, 68, 4511.
54. P. Sloth and T. S. Sørensen, J. Phys. Chem., 1990, 94, 2116.
55. E. Gutiérrez-Valladares, M. Lukšič, B. Millán-Malo, B. Hribar-Lee and V. Vlachy, Condens. Matter Phys., 2011, 14, 33003.
56. T. L. Croxton and D. A. McQuarrie, J. Phys. Chem., 1979, 83, 1840.
57. T. Ichiye and A. D. J. Haymet, J. Chem. Phys., 1990, 93, 8954.
58. D.-M. Duh and A. D. J. Haymet, J. Chem. Phys., 1992, 97, 7716.
59. C. S. Babu and T. Ichiye, J. Chem. Phys., 1994, 100, 9147.
60. T.-H. Chung and L. L. Lee, J. Chem. Phys., 2009, 130, 134513.
61. J. W. Zwanikken, P. K. Jha and M. Olvera de la Cruz, J. Chem. Phys., 2011, 135, 064106.
62. D. M. Burley, V. C. L. Hutson and C. W. Outhwaite, Mol. Phys., 1972, 23, 867.
63. C. W. Outhwaite, M. Molero and L. B. Bhuiyan, J. Chem. Soc., Faraday Trans., 1993, 89, 1315.
64. A. O. Quiñones, L. B. Bhuiyan and C. W. Outhwaite, Condens. Matter Phys., 2018, 21, 23802.
65. M. Su, Z. Xu and Y. Wang, J. Phys.: Condens. Matter, 2019, 31, 355101.
66. A. L. Kholodenko and A. L. Beyerlein, Phys. Rev. A: At., Mol., Opt. Phys., 1986, 34, 3309.
67. D. di Caprio, J. Stafiej and J. P. Badiali, J. Chem. Phys., 1998, 108, 8572.
68. R. R. Netz and H. Orland, Europhys. Lett., 1999, 45, 726.
69. R. R. Netz and H. Orland, Eur. Phys. J. E: Soft Matter Biol. Phys., 2000, 1, 67.
70. J.-M. Caillol, J. Stat. Phys., 2004, 115, 1461.
71. D. N. Card and J. P. Valleau, J. Chem. Phys., 1970, 52, 6232.
72. J. P. Valleau and L. K. Cohen, J. Chem. Phys., 1980, 72, 5935.
73. J. P. Valleau, L. K. Cohen and D. N. Card, J. Chem. Phys., 1980, 72, 5942.
74. W. van Megen and I. K. Snook, Mol. Phys., 1980, 39, 1043.
75. P. Sloth, T. S. Sørensen and J. B. Jensen, J. Chem. Soc., Faraday Trans. 2, 1987, 83, 881.
76. T. S. Sørensen, P. Sloth, H. B. Nielsen and J. B. Jensen, Acta Chem. Scand., Ser. A, 1988, 42, 237.
77. B. R. Svensson and C. E. Woodward, Mol. Phys., 1988, 64, 247.
78. Z. Abbas, E. Ahlberg and S. Nordholm, Fluid Phase Equilib., 2007, 260, 233.
79. Z. Abbas, E. Ahlberg and S. Nordholm, J. Phys. Chem. B, 2009, 113, 5905.
80. J. Janecek and R. R. Netz, J. Chem. Phys., 2009, 130, 074502.
81. R. Kjellander and J. Ulander, Mol. Phys., 1998, 95, 495.
82. A. McBride, M. Kohonen and P. Attard, J. Chem. Phys., 1998, 109, 2423.
83. B. Forsberg, J. Ulander and R. Kjellander, J. Chem. Phys., 2005, 122, 064502.
84. P. Attard, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1993, 48, 3604.
85. N. Bjerrum, Kgl. Dan. Vidensk. Selsk. Mat.-Fys. Medd., 1926, 7, 1.
86. H. Yokoyama and H. Yamatera, Bull. Chem. Soc. Jpn., 1975, 48, 1770.
87. M. E. Fisher and Y. Levin, Phys. Rev. Lett., 1993, 71, 3826.
88. Y. Levin and M. E. Fisher, Physica A, 1996, 225, 164.
89. W. Ebeling and M. Grigo, Ann. Phys., 1980, 37, 21.
90. T. Cartailler, P. Turq, L. Blum and N. Condamine, J. Phys. Chem., 1992, 96, 6766.
91. Y. Zhou, S. Yeh and G. Stell, J. Chem. Phys., 1995, 102, 5785.
92. S. Yeh, Y. Zhou and G. Stell, J. Phys. Chem., 1996, 100, 1415.
93. D. M. Zuckerman, M. E. Fisher and B. P. Lee, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1997, 56, 6569.
94. R. Kjellander, J. Chem. Phys., 2016, 145, 124503.
95. R. Kjellander, J. Chem. Phys., 2018, 148, 193701.
96. Z. Abbas, M. Gunnarsson, E. Ahlberg and S. Nordholm, J. Colloid Interface Sci., 2001, 243, 11.
97. N. F. Carnahan and K. E. Starling, J. Chem. Phys., 1969, 51, 635.
98. R. Kjellander, Soft Matter, 2019, 15, 5866.
99. M. A. Gebbie, A. M. Smith, H. A. Dobbs, A. A. Lee, G. G. Warr, X. Banquy, M. Valtiner, M. W. Rutland, J. N. Israelachvili, S. Perkin and R. Atkin, Chem. Commun., 2017, 53, 1214.
100. R. Kjellander and D. J. Mitchell, Chem. Phys. Lett., 1992, 200, 76.
101. R. Kjellander and D. J. Mitchell, J. Chem. Phys., 1994, 101, 603.
102. D. J. Mitchell and B. W. Ninham, Chem. Phys. Lett., 1978, 53, 397.
103. P. Kékicheff and B. W. Ninham, Europhys. Lett., 1990, 12, 471.
104. A. M. Smith, P. Maroni, G. Trefalt and M. Borkovec, J. Phys. Chem. B, 2019, 123, 1733.
105. B. Stojimirović, M. Galli and G. Trefalt, Phys. Rev. Res., 2020, 2, 023315.
106. R. Kjellander and R. Ramirez, J. Phys.: Condens. Matter, 2005, 17, S3409.
107. D. D. Carley, J. Chem. Phys., 1967, 46, 3783.
108. F. H. Stillinger and R. Lovett, J. Chem. Phys., 1968, 48, 3858.
109. L. Verlet and J.-J. Weis, Phys. Rev. A: At., Mol., Opt. Phys., 1972, 5, 939.
110. W. R. Smith, D. J. Henderson, P. J. Leonard, J. A. Barker and E. W. Grundke, Mol. Phys., 2008, 106, 3.

---

Footnotes

† In fact, κ_eff is evaluated in ref. 80 from the simulations via a route that implies that κ_eff satisfies Attard's approximate SCSL⁸⁴ equation for κ [eqn (38) in the current work]. This is the reason why κ_eff from the simulations agrees with κ from the SCSL approximation; see the authors' remark about this agreement in the text below eqn (47) in ref. 80. Any difference between the two values is due numerical noise in the simulation.

‡ An exception is charge-inversion invariant systems, for example, the restricted primitive model where anions and cations are identical apart from their sign of charge. In real systems, the anions and cations differ by more than the signs of their charges, so the systems are not charge-inversion invariant.

§ In this work “exact relationship,” “exact equation” etc. mean that they are derived without any approximations in statistical mechanics for a given model, that is, for a given Hamiltonian that defines the system.

¶ The term proportional to Λ in eqn (5) was originally derived in ref. 102.

|| The function ϕ_Coul*(r) is, in fact, a Green's function for spatial propagation of electrostatic field the electrolyte,⁹⁸ which is the proper way to define it.

** The notation used here differs from that in the original DIT papers ref. 22, 23, 100 and 101, where the product is denoted as E or and the present q^eff_j is denoted as q_j*, while q_j* in the present work denotes something different (see Section 7).

†† Examples of such terms include terms that decay like f(r)e^−2κr/r², where f(r) is a slowly varying function,^23,101 and cross-terms between the various modes (cf.ref. 98). Among these terms there are also core–core correlation contributions, which in the restricted primitive model can be distinguished from the electrostatic contributions. In the approximation used in the present work, the core–core contributions are approximated by a hard sphere correlation term w^core, see eqn (52).

‡‡ Eqn (24) can be formulated in several equivalent manners in DIT. One of them is an equation based on the static dielectric function [small epsilon, Greek, tilde](k), namely [small epsilon, Greek, tilde](iκ) = 0, where i is the imaginary unit. The decay modes are hence given by the zeros of [small epsilon, Greek, tilde](k) in complex Fourier space, which are poles of the Fourier transforms w̃_ij(k), [small psi, Greek, tilde]_i(k) and [small rho, Greek, tilde]_i(k).^23,100,101

§§ See footnote ‡‡ above.

¶¶ There also exist other complex-valued solutions to eqn (33), but here we are only concerned with those two that are connected continuously with the two real solutions κ and κ′.

|||| The HNC data for κ_Da ≳ 1.5 in ref. 22 were not reliable and have been ignored in the plot.

*** In the present case βw^el(r) is equal to the square bracket in eqn (107) in Appendix A.

††† In absolute values, ln γ^el_± give 8% of ln γ_± and ϕ^el 3% of ϕ in this system.

‡‡‡ Note that for any constant C > 0 we have Λ² ln(CΛ) = Λ² ln(Λ) + Λ² ln(C) ∼ Λ² ln(Λ) when Λ → 0. This means that any factor C can be included in the logarithm and the limiting law is still valid. For eqn (73) to be useful apart from at high dilution, one must therefore at least add the next order asymptotic term, which proportional to Λ², that is, proportional to the ion concentration. That term, which depends on the ionic characteristics like the ionic size, specifies an appropriate C. No exact simple expression for that term is, however, presently known, but we will later derive an approximate expression for it, eqn (81). We may note that z⁴Λ² ln(Λ) = (β_Rτ)²ln(β_Rτ/z²), which means that for two systems with different ionic valencies but with the same value of β_R, eqn (73) does not give the same asymptote as function of τ except in the limit of infinite dilution. This peculiarity occurs despite that the statistical mechanics in the RPM is determined solely by β_R and τ and it can easily be eliminated by selecting C = z². Eqn (73) is, however, derived for use in the infinite dilution limit only, where the value of C does not matter. A better selection of C is done in eqn (81), which is valid for a larger range of concentrations; see also footnote ¶¶¶.

§§§ The term −(βw^el_ij)³/3! corresponds to the low concentration limit of the term (h_ij^#)³/3! in eqn (103) of ref. 101 in the derivation of the exact limiting law for κ²/κ_D².

¶¶¶ Note that 1 + β_R²τ² ln(τ)/6 = 1 + z⁴Λ² ln(κ_Da)/6 contains the hard core diameter a and is not a universal limiting law for κ²/κ_D² valid for all symmetric electrolytes because it presumes that the ions are hard spheres. The universal law (73) instead contains the factor ln(κ_D[small script l]_B).

|||||| A favorable comparison between the experimental data and the exact asymptotic expression (73) was done in ref. 104 for the aqueous solutions. Such a comparison could not be done for the isopropanol solutions for reasons outlined in footnote ‡‡‡, but eqn (81) can be used for both kinds of solutions since it is expressed solely in terms of τ and β_R.

**** The notation used here differs from that in ref. 22, 23, 100 and 101, where q_j* denotes the effective charge q^eff_j of the present work.

†††† Both eqn (87) and (88) are mathematically equivalent to the equation [small epsilon, Greek, tilde](iκ) = 0, cf. footnote ‡‡ in Section 2.

---