Open Access Article

This Open Access Article is licensed under a Creative Commons Attribution-Non Commercial 3.0 Unported Licence

Roland
Kjellander

Department of Chemistry and Molecular Biology, University of Gothenburg, SE-412 96 Gothenburg, Sweden. E-mail: roland.kjellander@gu.se; Tel: +46 31 7869027

Received
8th April 2019
, Accepted 4th June 2019

First published on 10th June 2019

A general, exact theory for the decay of interactions between any particles immersed in electrolytes, including surface forces between macroscopic bodies, is derived in a self-contained, physically transparent manner. It is valid for electrolytes at any density, including ionic gases, molten salts, ionic liquids, and electrolyte solutions with molecular solvent at any concentration. The ions, the solvent and any other particles in the system can have any sizes, any shapes and arbitrary internal charge distributions. The spatial propagation of the interactions in electrolytes has several decay modes with different decay lengths that are given by the solutions, κ_{ν}, ν = 1, 2,…, to a general equation for the screening parameter κ; an equation that describes the dielectric response. There can exist simultaneous decay modes with plain exponential decay and modes with damped oscillatory exponential decay, as observed experimentally and theoretically. In the limit of zero ionic density, the decay length 1/κ_{ν} of the mode with the longest range approaches the Debye length 1/κ_{D}. The coupling between fluctuations in number density and charge density, described by the density–charge correlation function H_{NQ}(r), makes all decay modes of pair correlations and interaction free energies identical to those of the screened electrostatic potential, and hence they have the same values for the screening parameters. The density–density and charge–charge correlation functions, H_{NN}(r) and H_{QQ}(r), also have these decay modes. For the exceptional case of charge-inversion invariant systems, H_{NQ}(r) is identically zero for symmetry reasons and H_{NN}(r) and H_{QQ}(r) have, instead, decay modes with different decay lengths.

Oscillatory pair distribution functions and, consequently, oscillatory free energy of intermolecular interaction (potential of mean force) are well-established features of dense electrolytes like molten salts. The recent discovery^{1,2} of monotonic, exponentially decaying long-range forces with decay lengths of 4–11 nm in ionic liquids was therefore quite surprising, not least because traditional theories of electrolytes give screening lengths that are orders of magnitude shorter under the conditions in question. Such forces with long decay lengths have been observed for various systems with high ionic densities, like ionic liquids and concentrated electrolyte solutions.^{3–10} Simultaneously, there are often oscillatory contributions with shorter decay lengths simultaneously present in the force curves for such systems.

These observations are conceptually important because exponentially decaying forces between surfaces in electrolyte solutions have often been taken as an experimental verification of the correctness of the Poisson–Boltzmann (PB) approximation for surface interactions, which predicts such forces for large distances. The PB expressions that are commonly used contain adjustable parameters like surface charge density, surface potential etc. that can be fitted to the experimental data. However, any reasonable theory gives exponentially decaying forces in electrolytes at least for low ionic densities and exact analysis of statistical mechanics of fluids says that the forces must be exponential under such conditions – apart from contributions due to dispersion interactions that have power law decay and therefore ultimately must dominate the forces for sufficiently large distances. Therefore the mere existence of an exponential decay says nothing about the validity of the PB approximation. In this connection we should note that an exponential decay of the interactions in planar geometry, i.e., proportional to e^{−κℓ}, where ℓ is the distance and κ is the screening parameter, translates for spherical geometry into a Yukawa function decay, e^{−κr}/r, where r is the radial distance. Furthermore, it is a statistical mechanical fact that the exponential decay of surface forces for large separations has the same parameter κ as the pair correlations in the bulk fluid in equilibrium with the fluid in the slit between the surfaces. The same applies for exponentially damped oscillatory forces.

Another test of whether the PB approximation is applicable or not is the magnitude of the decay length, which in this approximation is given by the Debye length 1/κ_{D} where the Debye screening parameter κ_{D} is defined for electrolyte solutions from

(1) |

The PB approximation is quite accurate for low ionic densities. The actual decay length for the electrolyte approaches 1/κ_{D} when the ion density goes to zero and the importance of the ion–ion correlations decreases. A pertinent question is how low the ion density must be for the correlations to be unimportant. In fact, long-range electrostatic ion–ion correlations can give substantial deviation from the Debye length also for low ion densities provided the electrostatic coupling is sufficiently strong, for instance for multivalent ions in dilute aqueous solution as we will see later.

Furthermore, in contrast to the PB prediction of monotonic exponential decay with a single decay length, there appear in general more than one decay length and/or exponentially damped oscillatory decay. This can be illustrated by a simple approximation^{11,12} that is valid for simple ionic fluids with rather low electrostatic coupling, for example ionic fluids at very high temperatures or electrolyte solutions in the primitive model at high ε_{r} for the dielectric continuum solvent; typical examples are monovalent ions in aqueous solution at room temperature and such ions in vacuum at T = 23400 K (these systems have the same value of ε_{r}T). For a symmetric electrolyte with ionic diameter d this approximation gives the following expression for the screening parameter κ

(2) |

(3) |

Fig. 1 A plot of the screening parameters divided by the Debye parameter κ_{D} for a 1:1 aqueous electrolyte solution at room temperature as functions of κ_{D}d, which is proportional to the square root of the electrolyte concentration. The decay length is smaller than the Debye length when κ/κ_{D} > 1. The full curves show the screening parameter κ that gives the leading decay length 1/κ for low concentrations. The other curves show the screening parameters for other decay modes of the system (see text). The thick curves are the results of the simple approximation in eqn (2), while the thin curves are from HNC calculations^{18–20} that are very accurate for this system. The open symbols show Monte Carlo data^{18} for the same system. The filled symbols show the cross-over point between monotonic exponential and exponentially damped oscillatory decays, i.e., real and complex screening parameters, respectively, as calculated in the HNC approximation and from eqn (2). |

The simple expression in eqn (2) is surprisingly accurate, considering its humble origin; in the figure its predictions are compared with the results from Monte Carlo (MC) simulations^{18} and Hypernetted Chain (HNC) calculations^{18–20} and it is seen that the agreement is nearly quantitative. For the oscillatory decay, the decay length 1/κ_{ℜ} increases quite rapidly with electrolyte concentration after the cross-over point. The wavelength 2π/κ_{ℑ} immediately after the cross-over is very long (κ_{ℑ} = 0 at the cross-over point), much longer than the decay length, so in practice the oscillations would not be observable. However, the wavelength decreases very rapidly with increased concentration and is quite soon comparable to the decay length. Therefore, for the forces between particles or between surfaces at increasing separations, some oscillations should be seen before their magnitude is too small.

The important points here are that there exist more than one decay mode, one with decay length 1/κ and one with 1/κ′ and, furthermore, that the decay modes predicted by the equation for κ can be oscillatory. The objective of the simple expression (2) is to give a simple illustration of these facts. In the general case, other decay modes are also simultaneously present in electrolytes. These matters constitute the main theme of the current work. We will see how the decay modes can be treated in formally exact theory that is valid for much more complex systems. Such modes have a major role for the interactions in electrolytes for all conditions from low to high ionic densities and all temperatures (provided that the system remains fluid).

The approximation behind eqn (2) is not sufficient for higher electrostatic coupling than that of the example in the figure, i.e., lower temperatures, lower ε_{r} and/higher ionic valencies. For divalent ions in aqueous solution at room temperature, results from HNC calculations and MC simulations^{18–20} show that the decay length at low ion densities is appreciably longer than the Debye length. In a plot as that in Fig. 1, the curve for κ/κ_{D} then lies below the value 1 for small κ_{D}d; for example, for divalent ions with diameter d = 4.6 Å, κ/κ_{D} reaches a minimal value ≈0.84 before it increases to values above 1, whereafter the qualitative behavior is similar to that for the monovalent case. In fact, the results from the HNC approximation show that κ/κ_{D} lies below the value 1 for very small κ_{D}d also for the monovalent ions,^{12,20} but this is not visible for monovalent electrolytes on the scale in Fig. 1.

It is quite significant that deviations in decay length from the Debye length occur not only for dense electrolytes. This is a fundamental property of electrolytes because in the limit of low ionic density, the ratio κ/κ_{D} for a z:z electrolyte (z is the valency) satisfies the exact limiting law^{20,21}

(4) |

This law was recently verified experimentally by Smith et al.^{23} for dilute 2:2 and 3:3 electrolytes in aqueous solution. Large negative deviations of κ/κ_{D} from the value 1 were observed for electrolyte concentrations in the interval 0.1–10 mM and the agreement with eqn (4) was nearly quantitative for a large part of this interval. For aqueous 1:1 electrolytes, however, the deviation from the value 1 was very small.

The effect described by eqn (4) is not included in the approximation given in eqn (2) or in any other linear theory like MSA and LMPB, so κ/κ_{D} > 1 for small κ_{D}d in these approximations. However, if one includes nonlinear terms,^{12} one obtains agreement with eqn (4) in the limit of small Λ. Since the HNC approximation is nonlinear, it is in agreement with eqn (4). The difference between κ/κ_{D} from eqn (2) and from the HNC approximation for small κ_{D}d is very small for the monovalent electrolyte in Fig. 1. Due to the factor of z^{4} in eqn (4) the deviation from the Debye length increases rapidly with ionic valency for systems at low ion densities. For systems with much lower ε_{r} and/or lower temperatures, the deviation will be substantial also for monovalent ions.

Other very relevant illustrations of the subject that is studied formally in the current work can be taken from the computer simulation work by Keblinski et al.,^{24} where they investigate NaCl for a large variety of conditions from thin gases to molten salt for various temperatures. They used an empirical, realistic model for the pair interaction potential for NaCl that has been shown to reproduce quite well the thermodynamic, structural and other properties of molten NaCl.

Keblinski and coworkers calculated the charge–charge and density–density correlation functions, H_{QQ}(r) and H_{NN}(r) respectively,¶ and determined the decay lengths, wavelengths and other characteristics of the dominant contributions to these functions for various conditions. They found that the leading term of the asymptotic decay for large r of the correlation functions can dominate also for shorter distances – in some cases down to a distance of one or two particle diameters. Their Fig. 1 in ref. 12 shows that in molten NaCl at 1000 K, the leading asymptotic term in H_{QQ}(r), which decays like e^{−κℜr}cos(κ_{ℑ}r + γ)/r, is practically indistinguishable from H_{QQ}(r) down to r ≈ 2d, where d is the average ion diameter. Furthermore, it gives a very good description of H_{QQ}(r) down to r ≈ d. This means that an asymptotic analysis of the correlation function gives important information not only for very large r, but also for a quite wide range of shorter distances. It can in many cases be sufficient to include a couple of leading decay modes in the asymptotic analysis. Such a dominance of the leading asymptotic terms has also been found in earlier studies of other electrolyte systems.^{18,20,25}

The simulations by Keblinski and coworkers also included a study of the Kirkwood cross-over and they showed explicitly the presence of two exponentially decaying modes in in H_{QQ}(r) with decay lengths 1/κ and 1/κ′ (in our notation) for densities lower than the cross-over point and the leading oscillatory mode after the cross-over (see Fig. 5 in ref. 24). Furthermore, they showed the existence of two simultaneous decay modes of different types in the pair distribution function g_{ij}(r), an oscillatory mode decaying as e^{−κℜr}cos(κ_{ℑ}r + γ)/r and a monotonic one decaying as e^{−κ″r}/r. The decay lengths 1/κ_{ℜ} and 1/κ″ vary with ion density so that for some densities 1/κ_{ℜ} is smaller than 1/κ″ and for other densities the reverse is true (see Fig. 7 in ref. 24, where 1/κ_{ℑ} is denoted by λ_{Q} and 1/κ″ is denoted by λ_{G}). The density where the decay lengths are equal is a so-called Fisher–Widom cross-over point,^{13} where changes in the decay lengths of g_{ij}(r) make an oscillatory term to become leading for large distances instead of a monotonic term or vice versa. For NaCl at T = 3000 K, which is close to the critical temperature, this cross-over occurs at the density n^{b}_{tot} = n^{b}_{+} + n^{b}_{−} ≈ 0.1d^{−3}. The monotonic term has the largest decay length for n^{b}_{tot} < 0.1d^{−3}, which is due to the proximity of criticality where the density–density correlations have a long range, and the oscillatory term has the largest decay length above this density.

In their analysis of the simulation data, Keblinski and coworkers extracted the decay length for the oscillatory decay mode from H_{QQ}(r) and that for the monotonic mode from H_{NN}(r). Therefore they denoted the former as the “screening length” λ_{Q} and the latter as the “density–density correlation length” λ_{G}. A very important point to be made here is that there are terms in both correlation functions with these decay lengths, i.e., there is an oscillatory term in both H_{QQ}(r) and H_{NN}(r) with decay length λ_{Q} and likewise a monotonic term in both functions with decay length λ_{G}. The magnitude of the prefactor for the term with decay length λ_{Q} in H_{NN}(r) is, however, small for the NaCl system and likewise the term with decay length λ_{N} is small in H_{QQ}(r), so these contributions are rather insignificant numerically. As we will see, this is due to the fact that g_{++}(r) ≈ g_{−−}(r) in this system. For other systems the magnitudes of these kinds of terms in H_{QQ}(r) and H_{NN}(r) do not differ that much in general. In some cases they can have similar magnitudes.

In general, the leading terms of H_{QQ}(r) and H_{NN}(r) for large r have the same decay length. This has been shown for binary electrolytes with spherical ions of different sizes or different valencies in ref. 13, 25 and 26. We will show in the current work that this applies to all decay modes in electrolytes of almost any kind, so the decay lengths are therefore both screening lengths and density–density correlation lengths – in this work they are simply called “screening lengths” and the κ parameters are called “screening parameters.”

As we will see, there exist a few exceptions, for example model systems where the anions and cations are identical apart from the sign of their charges. The restricted primitive model (RPM) is such a model because all ions are charged hard spheres of equal diameter and have the same absolute valency. Then g_{++}(r) = g_{−−}(r) and, as a consequence of this symmetry, the density–charge correlation function H_{NQ}(r) is identically equal to zero and H_{QQ}(r) and H_{NN}(r) have different decay lengths. It is quite astounding that the most common model for electrolytes is an exception as regards its screening behavior! However, as we saw from the case of NaCl where g_{++}(r) ≈ g_{−−}(r), terms in H_{QQ}(r) and H_{NN}(r) with the same decay length have very different magnitudes and one dominates over the other, which means that results from models with g_{++}(r) = g_{−−}(r) may be quite reasonable as an approximation. These matters will be investigated for the general case in the current work and we will see that the class of systems that constitute this kind of exception consists of systems that are invariant when one reverses the sign of all charges in the system (including those in polar molecules). They will be called “charge-inversion invariant systems.”

The normal, realistic cases are, of course, systems that do not have such an invariance. Anions and cations virtually always differ by much more than the sign of their charges and, in addition, the positive and negative parts of the solvent molecules (if present) normally differ a lot apart from the sign of the charge. Thus H_{QQ}(r) and H_{NN}(r) normally have contributions with exactly the same set of decay lengths as the electrostatic potential and so do g_{ij} and the free energy of interaction. It may appear puzzling that H_{QQ}(r) and H_{NN}(r) have the same decay length also in cases where long-range density fluctuations arise upon approach to a critical point. Since both functions have equal range, H_{QQ}(r) is also an equally long-range function but we will see that the ratio H_{QQ}(r)/H_{NN}(r) for large r decreases quickly with increasing decay length, so the influence of H_{QQ}(r) vanishes when the decay length diverges. We must, however, emphasize that the present analysis does not cover critical phenomena, so critical points are not included.

As mentioned earlier, simultaneous decay modes with long-range monotonic and shorter-range oscillatory exponential decays have been observed in ionic liquids and concentrated electrolyte solutions. The presence of oscillatory, exponentially damped contributions to the correlation functions reflect in many cases the granularity of the system at the molecular level. The wavelength can have about the same magnitude as the average molecular size or, for ionic systems, a combination of the anion and cation diameters. This is, however, not always the case, for example when the particle sizes are very different or when the electrostatic correlations dominate.

It is important to note that oscillatory surface forces that extend several oscillations, like those observed in electrolyte solutions and ionic liquids, have decay lengths and wavelengths that are the same as those in the bulk fluid in equilibrium with the fluid between the surfaces. It is common to think of such oscillations as being specifically caused by a layered structure close to the surfaces, but the corresponding structuring is present also for the pair distributions in the bulk. A specific structure may be induced by the surfaces that causes a few oscillations in the surface force for very short surface separations, but in order for the oscillations to continue when the separation is increased, they must be supported by a decay mode of the bulk phase. As we have seen, the leading modes of the decay of the correlation functions often dominate not only for very large distances, but also for surprisingly short ones. This applies also for surface forces. Therefore, the bulk properties of electrolytes that are analyzed in the current work have great significance also for surface forces.

Oscillatory contributions are expected for most dense liquid systems, including dilute electrolyte solutions with molecular solvent. The oscillations are then dominated by the structure of the solvent. This has been observed experimentally by Smith et al.^{9} in surface force experiments for mixtures of an ionic liquid and a polar solvent. Simultaneously, there was a monotonic long-range decay mode. Their findings regarding the oscillatory contributions (but not the long-range monotonic one) have been illustrated theoretically in calculations by Coupette et al.^{27,28} for a semi-primitive model of electrolyte solutions, where the monovalent ions and the solvent are hard spheres. They studied systems with equally-sized ions,^{27} where the model is a charge-inversion invariant system, and systems with unequally-sized ions.^{28} Several decay modes were found and analyzed. For the case of equally-sized ions, the charge–charge and the density–density correlations have different decay lengths and different wavelengths, while for unequally-sized ions all correlation functions have common lengths. The wavelengths and the decay lengths were determined by the bulk phase, as always when the oscillations are maintained for increased surface separation. These observations are in agreement with the findings in the current work and in ref. 29 and 31. The presence of multiple decay modes and the behavior of the decay lengths for this system are also in line with previous studies of the primitive model (without discrete solvent) in ref. 13, 18–20 and 25.

The oscillatory contributions due to solvent molecules must be present also for dilute aqueous solutions of simple salts and, as argued in ref. 29, such oscillatory forces determined by the bulk phase have been observed in both surface force experiments and in theoretical investigations of bulk systems with discrete water molecules; the oscillatory surface forces are therefore not caused by “hydration forces” specifically associated with the surfaces. Thus, for dilute electrolyte solutions in general, there are simultaneous decay modes with long-range monotonic decay (i.e., the traditional decay mode with κ ≈ κ_{D}) and a shorter-range oscillatory exponential decay mode with a complex-valued screening parameter – the latter mode due primarily to the structure of the solvent. Most importantly, both of these decay modes are present in the correlation functions, the mean electrostatic potential and the interaction free energy.

The fact that the coupling between fluctuations in number density and fluctuations in charge density leads to equal decay lengths for the screened electrostatic potential, and all correlation functions, including H_{NN}(r), are analyzed thoroughly in the present work. We will prove that all decay modes in electrolytes, both the monotonic and the oscillatory exponential ones, originate in the general case from the same fundamental equation, which describes the relevant dielectric response of the system and involves the static dielectric function (k), where k is the wave number. This equation, which constitutes the general exact equation for the screening parameter κ, can be written in several equivalent manners. Perhaps the most appealing manner is the following expression^{30,31} that is similar to eqn (1) for the Debye parameter

(5) |

The exact theory presented in the present paper is very general and is valid for electrolytes from low to high ion densities: from thin gases to dense liquids and from dilute to concentrated electrolyte solutions with molecular solvent. The ions, the solvent and any other particles can have any sizes, arbitrary shapes, any internal charge distributions and can interact with any reasonable nonelectrostatic pair interactions. For simplicity, the particles are, however, assumed to be rigid and unpolarizable. The theory is a generalization of the Dressed Ion Theory (DIT)^{21,25,32} and the Dressed Molecule Theory (DMT)^{33–35} and it extends previous related work on the interactions in electrolytes.^{29–31} The basis of the theory is derived in the present paper in a new self-contained manner, which allows it to be done in a physically intuitive, yet correct way without use of advanced statistical mechanics. This derivation should allow a much broader readership than otherwise would be possible. It also allows a gradual introduction of the principally new results of the present work in an equally transparent manner. Although this theory treats quite complex phenomena and therefore is intrinsically complicated, the introduction of the concept of dressed particles, which is a key element, allows the general exact features of electrolytes to be cast in a form that is much simpler than otherwise is possible.

The contents of the paper are as follows. Section 2 contains a treatment of the linear polarization response of electrolytes due to the effect of weak external electrostatic fields. This is done in a special manner which focuses on the free energy of interaction for each constituent particle in the fluid and is formulated such that it is suitable as a basis for the general treatment of interparticle interactions in electrolytes. The nonlocal nature of electrostatic interactions in electrolytes is thereby demonstrated in a straightforward manner and is contrasted to the treatment of electrolytes in the PB approximation. The relationship to the static dielectric function (k) is given. In Section 3, the concept of a dressed particle is defined from an elementary analysis of the polarization response of electrolytes. The charge density of a dressed particle of species i is denoted by ρ_{i}* and it is an entity that is involved in all major aspects of screened interactions in electrolytes. The screened Coulomb potential ϕ_{Coul}*(r), which describes the spatial propagation of electrostatic interactions in electrolytes, is defined for the general case and can be expressed in terms of (k) in a simple manner. It is shown that one can use ϕ_{Coul}*(r) and ρ_{i}* in Coulomb's law to calculate the mean electrostatic potential ψ_{i} due to the particle. The general exact eqn (5) for the screening parameters κ_{ν} of the decay modes of ϕ_{Coul}*(r) is derived. The potential ψ_{i} has the same set of decay modes as ϕ_{Coul}*(r) and the different magnitudes of the modes of ψ_{i} are determined by projections of ρ_{i}* on each mode. The projections can be interpreted in terms of “effective” charges, as we will see. Section 4 deals with the pair distribution function g_{ij} and the free energy of interaction w_{ij} between particles in electrolytes, i.e., the pair potential of mean force. The distribution functions g_{ij}* of the dresses of the particles are also obtained. It is shown that all decay modes of w_{ij} and g_{ij} are, in general, the same as those of ϕ_{Coul}*(r) and are hence determined by the dielectric response of the electrolyte. This is caused by the coupling between charge density and number density fluctuations and implies that H_{QQ}(r), H_{NN}(r) and H_{NQ}(r) also have the same set of decay modes. The exception is charge-inversion invariant systems where the modes of H_{NN}(r) differ from those of H_{QQ}(r) and of ϕ_{Coul}*(r) since charge and density fluctuations in this case are independent of each other for symmetry reasons. Such invariant systems are investigated in some detail. Finally, the decays of H_{QQ}(r), H_{NN}(r) and H_{NQ}(r) for large r are investigated for the general case. Section 5 contains a summary of the conclusions of this work and also some perspectives on the use of the results in experimental and theoretical investigations.

The total electrostatic potential in the system is equal to Ψ(r) given by

δΨ(r) = δΨ^{ext}(r) + δΨ^{pol}(r). | (6) |

(7) |

The density distribution n_{i}(r′) can be expressed in terms of the potential of mean force W_{i}(r′) as the Boltzmann relationship

n_{i}(r′) = n^{b}_{i}e^{−βWi(r}′^{)}, | (8) |

δn_{i}(r′) = −βn_{i}(r′)δW_{i}(r′), | (9) |

For the special case of a bulk electrolyte, which is the main subject in this work, the density of species i is equal to n^{b}_{i} and we have ρ(r) = 0. Then, the electrostatic potential Ψ is constant and we set it equal to zero by convention. Likewise we initially have Ψ^{ext} = 0. Let us expose the system to an external electrostatic potential δΨ^{ext}(r) that is small everywhere in the electrolyte. This potential polarizes the bulk electrolyte and gives rise to a polarization charge density δρ^{pol} and hence to nonzero δΨ^{pol} and δΨ given by eqn (6) and (7). Henceforth, when we consider the polarization of a bulk electrolyte, all functions with a δ, like δρ^{pol}, δΨ^{pol}, δΨ and δΨ^{ext}, denote entities that are weak everywhere in the electrolyte; in principle they are infinitesimally small.

For a bulk fluid eqn (9) becomes

δn_{i}(r′) = −βn^{b}_{i}δW_{i}(r′) | (10) |

As we will show below

(11) |

(12) |

To show eqn (11) we will make use of an exact relationship in statistical mechanics for fluids called Yvon's first equation. It says that for a bulk fluid mixture that we expose to a weak external potential δv_{j}(r_{2}) acting on the various species j, we have

(13) |

(14) |

We can write the central charge q_{i} of the ion as a charge density q_{i}δ^{(3)}(r), where δ^{(3)}(r) is the three-dimensional Dirac delta function. By introducing ρ^{tot}_{i}(r) = q_{i}δ^{(3)}(r) + ρ_{i}(r), which is the total charge density of the ion itself and the surrounding ion cloud, we can write eqn (14) as

(15) |

(16) |

From eqn (15) it follows that the polarization charge density for a bulk electrolyte exposed to the weak electrostatic potential Ψ^{ext} is given by

(17) |

(18) |

We can express this as χ

Let us now turn to electrolytes consisting of nonspherical ions and other particles, for example solvent molecules. We will for simplicity solely treat the case of rigid, unpolarizable particles, but for each species i, the particle size, shape and internal charge density σ_{i} are arbitrary. For a particle with center of mass at r_{3} and orientation ω_{3}, the internal charge density at point r_{1} is given by σ_{i}(r_{1}|r_{3},ω_{3}). We use a normalized orientation variable ω so that where the integral is taken over all orientations.|| The charge of the particle is , which is independent of r_{3} and ω_{3}. Note that σ_{j}(r_{1}|r_{3},ω_{3}) for a given ω_{3} is a function of only r_{31} = r_{1} − r_{3}, where the vector r_{31} starts at the center of the particle. To simplify the notation we will henceforth write (r_{ν},ω_{ν}) ≡ R_{ν} whereby we have σ_{i}(r_{1}|r_{3},ω_{3}) ≡ σ_{i}(r_{1}|R_{3}), which is the charge density at r_{1} for a particle with coordinates R_{3}.

The pair interaction potential is u_{ij} = u^{el}_{ij} + u^{ne}_{ij}, which is the sum of the electrostatic (el) and nonelectrostatic (ne) interaction. The former is given by Coulomb's law as

(19) |

The number density of particles with center of mass at r_{3} and orientation ω_{3} is equal to n_{i}(r_{3},ω_{3}) ≡ n_{i}(R_{3}) and we have n_{i}(R_{3}) = n^{b}_{i}e^{−βWi(R3)}. The charge density of the system is

(20) |

As shown in Appendix A, the equation that corresponds to eqn (15) for a bulk fluid of nonspherical particles exposed to a weak external electrostatic potential δΨ^{ext} is

(21) |

ρ^{tot}_{i}(r_{2}|R_{1}) = σ_{i}(r_{2}|R_{1}) + ρ_{i}(r_{2}|R_{1}) | (22) |

(23) |

The physical interpretation of this formula is identical to that for spherical particles.

For nonspherical particles the expression for the polarization response function χ^{ρ} is somewhat more complicated than eqn (18). By inserting eqn (21) into eqn (20) we obtain

(24) |

We can write eqn (24) as

(25) |

The fact that the functions δΨ^{ext}, δρ^{pol}, and δΨ are linearly related to each other implies that there exists a function χ*(r) that expresses the linear relationship between δρ^{pol} and δΨ as

(26) |

In fact, for a given χ^{ρ}(r), the function χ*(r) is the solution to the equation

(27) |

(28) |

Eqn (27) is easily derived as follows. By inserting (6) into eqn (26) and then using eqn (7) we obtain

Since the left-hand side (lhs) is equal to and δΨ^{ext}(r_{4}) is arbitrary, eqn (27) follows.

It is illustrative to see what the PB approximation says about the polarization response. In this approximation we have for spherical ions δW_{i}(r_{1}) = q_{i}δΨ(r_{1}) and hence from eqn (10) we obtain δn_{i}(r_{1}) = −βn^{b}_{i}q_{i}δΨ(r_{1}). Therefore

(29) |

and by comparing with the expression (26) for χ*, we can identify

(30) |

In reality, we have nonlocal electrostatics in the sense that the polarization at one point r_{1} is influenced by the electrostatic potential at other points in the surroundings. This can be seen in eqn (26). The polarization response function χ*(r_{12}) is nonzero for r_{12} ≠ 0 so δρ^{pol}(r_{1}) is influenced by δΨ(r_{2}) for all points r_{2} in the neighborhood of r_{1}. This feature is caused by the correlations between the particles in the electrolyte. The nonlocal nature of the response can be understood in the following manner. Any ion located at r_{2} interacts with the electrostatic potential and since this ion correlates with other ions it will affect the probability for ions to be at r_{1}. The density of ions at r_{1} accordingly depends on the potential elsewhere. In other words, the potential of mean force δw_{i}(r_{1}) and hence the charge density δρ^{pol}(r_{1}) depend on the values of δΨ(r_{2}) for all points r_{2} in the vicinity of r_{1}. This fact is expressed by the nonlocality in the general exact relationships we have obtained.

The response functions χ^{ρ}(r) and χ*(r) are closely related to the dielectric response of the electrolyte in terms of the static dielectric function (k), where k is a wave number. To see this we introduce the electrostatic fields corresponding to Ψ(r) and Ψ^{ext}, which are E(r) = −∇Ψ(r) and E^{ext}(r) = −∇Ψ^{ext}(r), respectively. E is sometimes called the Maxwell field and E^{ext} can be expressed in terms of the displacement field D given by D(r) = ε_{0}E^{ext}(r), but we will use E^{ext} in the present work. As we have seen E^{ext} is the field from the external source in the absence of the fluid. The reasoning is performed most easily in Fourier space, so we will consider Ẽ(k) and Ẽ^{ext}(k).

The static (longitudinal) dielectric function (k) relates these electrostatic fields when they are weak and therefore linearly related to each other. In general we have for a homogeneous and isotropic fluid

(31) |

In line with the notation above, when we expose the electrolyte to a weak field we put a δ on the potentials and ^{ext}. Therefore we write δ and δ^{ext} instead of and ^{ext}, whereby eqn (31) implies

(32) |

(33) |

By comparing the last two equations with eqn (32) we can identify

(34) |

The nonlocal nature of electrostatics in the microscopic domain, which is described by the nonlocal response functions χ^{ρ}(r) and χ*(r), is hence included in the k dependence of the dielectric function (k). In the PB approximation, where electrostatics is local, we have from eqn (30), where κ_{D} is the Debye screening parameter for ions in vacuum given by

(35) |

(36) |

For completeness, we note that in electrostatics it is common to use the electric susceptibility χ_{e}, which gives the polarization density P in terms of the total field as = ε_{0}_{e}Ẽ for weak fields.†† In our notation we have P = −ε_{0}E^{pol}, where E^{pol}(r) = −∇Ψ^{pol}(r). For weak fields we thus have −ε_{0}δẼ^{pol} = ε_{0}_{e}δẼ or in terms of potentials δ^{pol} = −_{e}δ. By comparing with eqn (33) we see that

To minimize the number of symbols we use, we will refrain ourselves from using _{e} in what follows. Instead we use *(k) [in ordinary space χ*(r)], which plays a central role in the present work. Obviously, * has a prominent role also in ordinary electrostatics.

By inserting δΨ^{ext}(r) = δΨ(r) − δΨ^{pol}(r) into eqn (23) we obtain

(37) |

where

(38) |

(39) |

(40) |

(41) |

(42) |

We will now interpret the physical meaning of ρ_{i}*. Let us immerse a particle of species i with orientation ω_{1} into a bulk electrolyte at the point r_{1}. Initially the charge density is zero and after the immersion the density at r_{2} is ρ^{tot}_{i}(r_{2}|R_{1}). We have ρ^{tot}_{i} = σ_{i} + ρ_{i} [eqn (22)] and ρ_{i} can be regarded as the total polarization charge density that the particle induces in the surrounding electrolyte due to all kinds of interactions between this particle and the other particles (not only the polarization due to electrostatic interactions). Since these interactions are strong close to the particle, one cannot treat the polarization by linear response. If, however, the total electrostatic potential ψ_{i} due to the particle were weak, we would according to eqn (26) have the polarization charge density due to this potential

(43) |

ρ_{i}*(r_{2}|R_{1}) = ρ^{tot}_{i}(r_{2}|R_{1}) − ρ^{lin}_{i}(r_{2}|R_{1}), | (44) |

ρ_{i}*(r_{2}|R_{1}) = σ_{i}(r_{2}|R_{1}) + ρ^{dress}_{i}(r_{2}|R_{1}), | (45) |

The density ρ_{i}* plays a key role in what follows. For instance, as we have seen from eqn (23) and (41), the dressed charge density ρ_{i}* has the same role vis-à-vis the total potential δΨ as the total charge density ρ^{tot}_{i} has vis-à-vis the external potential δΨ^{ext}. Eqn (41) says that in the linear domain, the free energy of interaction δW_{i} of the i-ion with the total electrostatic field is equal to the electrostatic interaction energy between δΨ and the dressed particle charge–density ρ_{i}* associated with the ion. As we will see ρ_{i}* has very important roles in the interparticle interactions in electrolytes.

Note that ρ_{i}* for each i can be expressed in terms of ρ^{tot}_{j} for all j viaeqn (40) where ψ_{i} is given in terms of ρ^{tot}_{i} by eqn (38) and where χ* can be expressed in terms of χ^{ρ}viaeqn (28) and finally in terms of ρ^{tot}_{j}viaeqn (25). Thus, given ρ^{tot}_{j} for all j, one can calculate ρ_{i}* at least in principle.

Let us now return to the case when a bulk electrolyte is polarized by a weak electrostatic potential. Then, δn_{i} is given by eqn (42) and we can readily derive

(46) |

We can write eqn (46) in Fourier space by introducing the notation

f(r_{12},ω_{1}) ≡ f(r_{2}|r_{1},ω_{1}) = f(r_{2}|R_{1}) | (47) |

(48) |

(49) |

(50) |

Let us yet again connect with the PB approximation. For this purpose we apply the general, exact equations above to the special case of spherical ions with a charge q_{i} at the center, whereby eqn (43) becomes

(51) |

ρ_{i}*(r_{12}) = σ_{i}(r_{12}) + ρ^{dress}_{i}(r_{12}) = q_{i}δ^{(3)}(r_{12}) + ρ^{dress}_{i}(r_{12}), | (52) |

(53) |

(54) |

Let us consider a single nonspherical j-particle immersed in an electrolyte consisting of spherical ions of equal sizes. In this case the PB approximation says that

(55) |

From eqn (43) with χ*(r) given by eqn (30) we have

so from eqn (44) it follows that

(56) |

As noted earlier, these results in the PB approximation apply only to a single particle in an electrolyte. The surrounding ions are treated as bare point ions that do not correlate with each other, so they do not have any ion clouds of their own and no dresses. For these ions the potential of mean force is q_{i}ψ_{j} as assumed in the exponent of eqn (55). When the distribution of ions around an ion is calculated in the PB approximation, i.e., the pair distributions, one accordingly treats the latter ion differently than the rest. This unequal treatment leads to an infamous feature that g_{ij} ≠ g_{ji} in the (nonlinear) PB approximation.

As we saw in the simple example in Fig. 1 based on the approximation in eqn (2), if we instead treat all ions on an equal basis so all have dresses and, as we will see, therefore have effective charges different from the bare charges, the consequences are quite dramatic with the appearance of multiple screening parameters and oscillatory decay instead of a simple screening parameter κ_{D}.

3.2.1 The general Green's function for screened Coulomb interactions.
So far we have seen that the dressed particle charge–density ρ_{i}* plays an important role in the free energy of interaction of a particle with the total electrostatic potential in the linear domain [eqn (41)]. Furthermore ρ_{i}* determines the polarization response function χ* [eqn (46) and (54)] and the dielectric function (k) [eqn (50)]. We will now see that the dressed particle charge–density has yet another important role, namely as the source of the screened electrostatic potential from a particle.

with the boundary condition ψ_{i}(r|R_{1}) → 0 when r → ∞. By subtracting ρ^{lin}_{i}(r|R_{1}) from both sides of Poisson's equation we obtain

so we have, using eqn (43),

The solution ψ_{i} of this equation, which is linear in ψ_{i}, can be written in terms of the Green's function of the equation. The Green's function, which we denote by ϕ_{Coul}*(r), is by definition the solution of

and the solution of eqn (58) is thereby given by

as can be verified by inserting this expression into eqn (58) and using eqn (59). The Green's function ϕ_{Coul}*(r), which accordingly is defined by eqn (59), is called the (unit) screened Coulomb potential for the general case. It has a key role in the spatial propagation of electrostatic interactions in the electrolyte.

which has the solution

a monotonic Yukawa function. Therefore

in this approximation.

where we have used eqn (34). If we expose the electrolyte to an external electrostatic potential δΨ^{ext}(r) = δqϕ_{Coul}(r) ≡ Ψ^{ext}_{{δq}}(r), which is the potential at distance r from a small point charge δq, we see from eqn (32) with δ^{ext}(k) = δq_{Coul}(k) that the resulting total potential is δΨ(r) = δqϕ_{Coul}*(r) ≡ Ψ_{{δq}}(r) at least if r is sufficiently large and δq is sufficiently small, so the linear response is adequate. In this sense we may say that

for r > 0.

By using _{Coul}*(k) = _{Coul}(k)/(k) = 1/|ε_{0}k^{2}(k)] we see from the last equality that

when the factor 1/|ε_{0}k^{2}] has been removed. This simple relationship can also be derived from eqn (40).

which is the Fourier transform of eqn (62). The feature that the decay behavior of ϕ_{Coul}*(r) in the PB approximation is proportional to e^{−κDr}/r is associated with the fact that the denominator in eqn (67) is zero when k = ±iκ_{D}, that is, _{Coul}*(k) has simple poles at k = ±iκ_{D}. This holds in general,‡‡ so a pair of simple zeros k = ±iκ of the denominator of eqn (64) corresponds to a term in ϕ_{Coul}*(r) that decays as e^{−κr}/r. At such zeros, i.e., poles of _{Coul}*(k), we have

where the latter holds for κ ≠ 0. Since _{Coul}*(k) is the Fourier transform of a real function of r, _{Coul}*(k) is an even function of k and therefore it is sufficient to consider a pole at iκ, whereby the pole at −iκ follows.

where we have used the Fourier transform of eqn (52) and the definition (35) of κ_{D}. It follows from eqn (69) that the dresses of the ions, which describe the effects of ion–ion correlations, make κ ≠ κ_{D}.

where A* is a constant that will be determined later [eqn (70) holds at least when the ion density is not too high]. In contrast to the PB result (62), this decay formula is only valid asymptotically for large r because, as we will see, there are other terms in ϕ_{Coul}*(r) that decay faster to zero than the contribution on the rhs. They give non-negligible contributions for small r.

The mean electrostatic potential ψ_{i}(r|R_{1}) from a particle with coordinates R_{1} can be obtained from the total charge density ρ^{tot}_{i} associated with the particle via Coulomb's law as expressed in eqn (38). This potential satisfies Poisson's equation

−ε_{0}∇^{2}ψ_{i}(r|R_{1}) = ρ^{tot}_{i}(r|R_{1}) | (57) |

−ε_{0}∇^{2}ψ_{i}(r|R_{1}) − ρ^{lin}_{i}(r|R_{1}) = ρ^{tot}_{i}(r|R_{1}) − ρ^{lin}_{i}(r|R_{1}) ≡ ρ_{i}*(r|R_{1}) |

(58) |

(59) |

(60) |

Incidentally we note that in the PB approximation, where χ*(r) given by eqn (30), eqn (59) becomes

−ε_{0}[∇^{2}ϕ_{Coul}*(r) − κ_{D}^{2}ϕ_{Coul}*(r)] = δ^{(3)}(r) (PB), | (61) |

(62) |

(63) |

In the general case, by taking the Fourier transform of eqn (59) we obtain

(64) |

for r > 0.

The right hand side (rhs) of eqn (60) has the same form as Coulomb's law in eqn (38), which is no surprise since the usual (unscreened) Coulomb potential ϕ_{Coul}(r) is a Green's function of Poisson's equation (57). Note that ψ_{i}(r|R_{1}) given by eqn (38) and (60) is exactly the same for all r; what differs in the two equations is the following: when ψ_{i} is expressed in terms of ϕ_{Coul}(r) the source is ρ^{tot}_{i} while when it is expressed in terms of ϕ_{Coul}*(r) the source is ρ_{i}*.

We can write eqn (38) and (60) in Fourier space by using the notation in eqn (47), whereby we obtain

_{i}(k,ω_{1}) = ^{tot}_{i}(k,ω_{1})_{Coul}(k) = _{i}*(k,ω_{1})_{Coul}*(k). | (65) |

(66) |

In the PB approximation where (k) is given by eqn (36) we have

(67) |

[ε_{0}k^{2} − *(k)]_{k=±iκ} = 0, (±iκ) = 0, | (68) |

In order to see what these facts lead to, let us first consider an electrolyte consisting of spherical simple ions in vacuum. From the Fourier transform of eqn (54) we see that _{Coul}*(k) has a pole k = iκ when

(69) |

When the density of the ions goes to zero, the charge density ρ_{i}(r) of the ion cloud goes to zero for each r and likewise the linear part ρ^{lin}_{i}(r), so ρ^{dress}_{i} ≡ ρ_{i} − ρ^{lin}_{i} goes to zero. Thus κ/κ_{D} → 1 in this limit and the PB result is approached. Likewise, in the same limit we have δW_{i}(r_{1}) ≈ q_{i}δΨ(r_{1}) from eqn (53), so the PB approximation is reasonable at least where the electrostatic potential is sufficiently small.

At finite densities the pole at k = iκ gives, as we will see, the leading term in the decay of ϕ_{Coul}*

(70) |

3.2.2 The general equation for the screening parameter κ vs. the dielectric function.
Let us now focus on the dielectric function (k) and use the condition (iκ) = 0 in eqn (68) for the pole of _{Coul}*(k). Since (k) = 1 − *(k)_{Coul}(k) we can write

where we have split (k) into two parts

which is singular when k → 0 (provided that *(0) ≠ 0), and

which is regular (non-singular) at k = 0 since its value there is

These integrals converge provided χ*(r) decays to zero sufficiently rapidly with increasing r.

where

is the dressed particle charge of a particle of species i, that is, the total charge of the dressed particle charge–density ρ_{i}*. This charge consists of the charge q_{i} of the ion and the total charge of its dress. Note that

because the total charge of ρ^{tot}_{i}(r) is zero due to the local electroneutrality condition. Since ρ^{lin}_{i} is the linear part of the polarization response due to ψ_{i} from the ion, q_{i}* is normally nonzero and has the same sign as q_{i}. It is, however, possible for q_{i}* for a species to change sign and thereby attain the opposite sign to q_{i}, so in rare cases q_{i}* for a species can fortuitously be zero for certain parameter values of the system. One can show§§ that for electrolytes, so it is not possible to have q_{i}q_{i}* < 0 for all ionic species simultaneously. In Appendix A of ref. 31 it is described how ρ_{i}*(r) and q_{i}* can be calculated from ρ^{tot}_{j}(r) for all j or from the pair distribution functions g_{ij}(r).

where we have used

which is independent of the orientation ω_{3}, and the corresponding relationship for σ_{i}. Again, q_{i}* is the total charge of the dressed particle charge–density ρ_{i}*. Thus we always have , where q_{i}* is the dressed particle charge.

Due to the presence of _{sing}(k), the dielectric function (k) for electrolytes diverges to infinity when k → 0. This divergence is commonly called perfect screening.

where

is the dielectric factor. Eqn (78) is the general equation for the screening parameter κ, eqn (5). It is obviously equivalent to (iκ) = 0. Note that eqn (78), which is an exact equation for κ with general applicability, has an appearance very similar to the definition of the Debye parameter κ_{D}; only the factor q_{i}* and the presence of _{r}*(κ) in the denominator differ. From eqn (73) it follows that

As we will see, the fact that _{r}*(κ) is a function of κ is crucial for the understanding of the properties of electrolytes. When the particle density goes to zero, we have q_{i}* → q_{i} and _{r}*(κ) → 1 with screening parameter κ near zero, so eqn (78) approaches the expression for κ_{D} in eqn (35) which implies that κ/κ_{D} → 1 in this limit. For electrolyte models with a dielectric continuum solvent, the initial 1 on the rhs of eqn (80) should be replaced by ε_{r} for the solvent.

which follows from eqn (35) and (78) and the fact that n^{b}_{+}q_{+} = −n^{b}_{−}q_{−}. Note that κ_{D} used here is calculated for particles in vacuum.

(71) |

(72) |

(73) |

These integrals converge provided χ*(r) decays to zero sufficiently rapidly with increasing r.

For the case of spherical simple ions we have from the Fourier transform of eqn (54)

(74) |

(75) |

because the total charge of ρ

For nonspherical particles eqn (49) yields for k = 0

(76) |

which is independent of the orientation ω

By inserting this result into eqn (72) we obtain

(77) |

Using eqn (77) we see that the condition for a pole, (iκ) = _{reg}(iκ) + _{sing}(iκ) = 0, is

(78) |

_{r}*(κ) = _{reg}(iκ) | (79) |

(80) |

For a binary electrolyte, the deviation in κ from κ_{D} can be obtained from

(81) |

3.2.3 Decay modes, multiple screening parameters and effective dielectric permittivities.
Let us now determine the coefficient A* in eqn (70). From eqn (64) and the results above we see that

Since the denominator is zero for k = ±iκ we can write

where A* is given by

Here_{reg}′(k) = d_{reg}(k)/dk and we have used l'Hospital's rule to obtain the last equality. If we define

it follows that A* = 1/[ε_{0}^{eff}_{r}(κ)] and hence we have

By comparing with the PB result in eqn (62) we see that a correct treatment of the correlations between all particles in the system has given rise to a change in magnitude of ϕ_{Coul}*(r) for large r by a factor 1/^{eff}_{r}(κ). In addition there is a change in the value of the screening parameter from κ_{D} to κ. Note that ^{eff}_{r}(κ) differs from _{r}*(κ) by the last term in eqn (82).

where the first expression yields eqn (82) after the differentiation and where we have used (iκ) = 0 to obtain the last equality. From eqn (73) and the first equality in eqn (84) we readily obtain

which can be compared with eqn (80). When the density of the particles in the electrolyte goes to zero, we have^{eff}_{r}(κ) → 1. Since κ/κ_{D} → 1 in this limit, eqn (83) approaches the PB result for ϕ_{Coul}*(r) in eqn (62), at least for sufficiently large r. (For electrolyte models with a dielectric continuum solvent, the initial 1 on the rhs should be replaced by ε_{r} for the solvent.)

where we have used the fact that q_{i}* ≈ q_{i} for a dilute electrolyte. In the last equality we have also identified κ_{D}^{2} for the electrolyte solution when the solvent is a dielectric continuum with dielectric constant ε_{r}, eqn (1).
_{D}, that is, the same as in the PB approximation when the pure solvent has dielectric constant ε_{r}.

where the “other terms” are analogous Yukawa function terms with shorter decay lengths [i.e., with other κ values that are solutions to eqn (78)] and functions with different functional dependences of r (more about this later). Note that each Yukawa term has its own value of ^{eff}_{r}.

where we have written ^{eff}_{r}(κ_{ℜ} + iκ_{ℑ}) = |^{eff}_{r}|e^{−iϑ} with a real ϑ_{}. The two poles hence give rise to an exponentially decaying, oscillatory term with decay length 1/κ_{ℜ}, wavelength 2π/κ_{ℑ} and phase shift −ϑ_{}. We will denote such a function as an “oscillatory Yukawa function.” The function e^{−κr}/r with real κ will be designated as a “monotonic Yukawa function.”

with wavelengths 1/κ > 1/κ′. For low ionic densities we have seen that κ ≈ κ_{D} and when the ionic density is increased, the two solutions κ and κ′ of eqn (78) approach each other as in the example of Fig. 1 and then merge when the density reaches the cross-over point. For even higher densities, κ and κ′ become two complex conjugate solutions, so the leading term for large r is oscillatory as shown in eqn (87).

Since the denominator is zero for k = ±iκ we can write

where A* is given by

Here

(82) |

(83) |

Alternative expressions for ^{eff}_{r} are

(84) |

which can be compared with eqn (80). When the density of the particles in the electrolyte goes to zero, we have

These results can be applied to the case of an electrolyte solution with a molecular solvent consisting of polar, uncharged molecules. Since such molecules have q_{i} = 0 they do not contribute to _{sing}(k), so the sum in eqn (77) runs in practice over the ionic species only. Furthermore, since these molecules have no net charge they do not contribute to q_{i}* for any species i. In the expression for κ, eqn (78), the solvent molecules solely contribute to the dielectric factor _{r}*(κ) [apart from their indirect influence on the distribution functions and other entities, of course].

For a pure polar solvent without ions, _{sing}(k) vanishes so (k) = _{reg}(k). As mentioned earlier the dielectric constant of a pure polar medium is defined microscopically as , so we have ε_{r} = (0) = _{reg}(0). In a very dilute electrolyte solution κ ≈ 0 and hence _{r}*(κ) ≈ _{r}*(0) ≈ ε_{r}. Eqn (78) then becomes

(85) |

In dilute solutions we have ^{eff}_{r}(κ) ≈ ^{eff}_{r}(0) ≈ ε_{r} since the last term in eqn (82) vanishes when κ → 0 because of the prefactor iκ. This implies that we approximately have when r → ∞

For electrolyte solutions in general, ^{eff}_{r}(κ) and _{r}*(κ) contain, as we have seen, contributions from the ions. Since 1/^{eff}_{r}(κ) determines the magnitude of the screened Coulomb potential for large r, ^{eff}_{r}(κ) can be designated as the effective relative dielectric permittivity of the entire electrolyte solution. Likewise, the dielectric factor _{r}*(κ) acts as a kind of relative dielectric permittivity of the solution since it takes the role that the dielectric constant of the pure solvent has in the expression for κ_{D}. Remember that ^{eff}_{r}(κ) and _{r}*(κ) are different from each other in general.

The only difference between the general equation for the screening parameter κ, eqn (78), and the definition (35) of the Debye parameter κ_{D} is, apart from the factor _{r}*(κ), that the former contains the factor q_{i}q_{i}* instead of q_{i}^{2}. While κ_{D} is given solely in terms of the system parameters n_{i}, q_{i}, and T, the actual screening parameter κ depends on the state-dependent entities q_{i}* and _{r}*(κ). The latter are, as we have seen, defined in terms of dressed charge densities or, equivalently, in terms of ρ^{tot}_{i} for all i via its relationship to ρ_{i}* as explained earlier.

In the expression (78) for κ, the dressed ion charge q_{i}* is a constant for each system with given system parameters, but the dielectric factor _{r}*(κ) is a function of κ. The latter fact makes a huge difference compared to the predictions of the PB approximation. While this approximation predicts a unique screening parameter κ_{D} from eqn (35), the exact equation, eqn (78), is an equation for κ. Apart from the solution κ, which gives the longest decay length, eqn (78) has in general several other solutions κ′, κ″ etc. (the three solutions can alternatively be denoted κ_{ν}, for ν = 1, 2 and 3). Each solution gives rise to a term in ϕ_{Coul}*(r) like the rhs of eqn (83), so we have¶¶

(86) |

Furthermore, solutions to eqn (78) can be complex-valued, in the case of which there are always two solutions that are complex conjugates to each other, say, κ = κ_{ℜ} + iκ_{ℑ} and κ′ = κ_{ℜ} + iκ_{ℑ}, where κ_{ℜ} and κ_{ℑ} are real. Since the sum of a complex number Z and its complex conjugate is given by Z + = 2ℜ(Z), where ℜ(Z) stands for the real part of Z, we then have

(87) |

Thus there exist several decay modes for the screened electric potential, some that give monotonic and others that give oscillatory decay. In the limit of low ionic densities, one of the modes approaches the single mode that is included in the PB approximations. Henceforth, we will use the term “screening parameters” to denote a set like κ, κ′, κ″ etc. that are solutions to eqn (78), so this term includes the inverse decay lengths κ or κ_{ℜ} and the inverse wavelengths κ_{ℑ}/2π of all monotonic and oscillatory Yukawa function contributions to ϕ_{Coul}*. More generally, the concept of “decay parameters” denotes the inverse decay lengths and inverse wavelengths of all kinds of contributions to the various functions.

A well-known example of the occurrence of an oscillatory term as in eqn (87) is the Kirkwood cross-over^{13} mentioned in the Introduction, where two real solutions turn into two complex-valued solutions when a system parameter like the ion density is changed. We take the example of two real solutions κ and κ′, where κ is the smallest and κ′ is the second smallest solution, i.e., for large r they give the two leading terms in ϕ_{Coul}*(r) as given by eqn (86) and we have

(88) |

Recall that eqn (78) is equivalent to (iκ) = 0, so before the cross-over κ and κ′ are two consecutive zeros of (iξ) as a function of the real variable ξ. As we have seen, ^{eff}_{r}(κ) > 0 for low ion densities, so from the rhs of eqn (84) we see that ′(iξ) > 0 for ξ = κ. The next zero of (iξ) must have an opposite derivative, so ′(iξ) < 0 for ξ = κ′, which implies that ^{eff}_{r}(κ′) < 0. Thus the two terms in eqn (88) have opposite signs. At the the cross-over point, where the two zeros of (iξ) merge, we must have ^{eff}_{r}(κ) = ^{eff}_{r}(κ′) = 0 but the sum of the two terms remains finite there.

There may also appear another kind of cross-over between a monotonic and an oscillatory Yukawa function decay, namely the Fisher–Widom cross-over^{13} mentioned in the Introduction. Say that the term with screening parameter κ″ in eqn (86) initially has a shorter decay length 1/κ″ than the oscillatory term we have just discussed. When the ionic density is increased it is possible that the former term becomes the leading term because 1/κ″ becomes larger than the decay length 1/κ_{ℜ} of the latter. This means that the decay behavior of ϕ_{Coul}*(r) for large r changes from oscillatory to monotonic at the density value where 1/κ″ and 1/κ_{ℜ} are equal. This kind of cross-over may, of course, also occur in the reverse direction, i.e., the decay changes from monotonic to oscillatory.

3.2.4 Multipolar effective charges.
The mean electrostatic potential ψ_{i} due to an i-particle is given by eqn (60) and for each Yukawa function term in ϕ_{Coul}*(r) there is a term in ψ_{i} with the same screening parameter, say κ_{ν},

Each of these contributions are like the single term in the PB result, eqn (63). The screening parameter κ_{ν} of some contributions may, as we have seen, be complex-valued and give rise to an oscillatory contribution to ψ_{i}.
_{12} and r′ = r_{13} (note that r − r′ = r_{2} − r_{3}) and using the notation introduced in eqn (47), we can conclude from eqn (89) in the limit r_{12} → ∞

provided that ρ_{i}* decays sufficiently rapidly with distance. Thus each contribution decays like a Yukawa function with distance r but has a magnitude that is different depending on the direction of the vector r_{12}, whereby the integral contains the direction dependence as expressed via_{12} in the exponent.

and we obtain from the terms in eqn (86) in the limit r_{12} → ∞

where q^{eff}_{i} is an “effective charge” of the i-particle. Note that each decay mode has its own value of the effective charge since q^{eff}_{i}(κ_{ν}) depends on κ_{ν}. Incidentally we also note that q^{eff}_{i} is different from the dressed particle charge q_{i}*, which is a constant that is independent of κ_{ν}.

where

is a direction dependent entity that has the unit of charge. We may call Q^{eff}_{i} a “multipolar effective charge” of the particle, where “multipolar” indicates that it is direction dependent. Each decay mode has its own value since Q^{eff}_{i} depends on κ_{ν}. The orientation independent part of the potential, the monopolar part, is given by the average

where

is the orientation average of the effective multipolar charge. ^{eff}_{i} can be described as a kind of monopolar effective charge of the particles. The rest of Q^{eff}_{i}(_{12},ω_{1},κ_{ν}) has an orientational angle dependence with a combination of dipolar, quadrupolar and higher multipolar characteristics.^{34} Note that electroneutral particles, like solvent molecules, acquire nonzero ^{eff}_{i} due to the presence of ions in their neighborhood, so they have nonzero effective charges in general.
_{ν} of _{Coul}*(k), that is, zero of (k), we obtain a contribution to ψ_{i}(r_{12},ω_{1}) that decays as e^{−κνr}/r with the same prefactor 1/[4πε_{0}^{eff}_{r}(κ_{ν})] as in eqn (86) times the factor

where we in ordinary space apply this result with = _{12} since the origin is selected at the particle center. We see that Q^{eff}_{i} for the mode with screening parameter κ_{ν} is the projection of ρ_{i}* on this mode.

where we have defined

Eqn (97) is an equation for κ that is equivalent to eqn (78).

which is similar to the general eqn (78). The differences are that q^{eff}_{i}(κ) depends on κ while q_{i}* is constant and that eqn (99) lacks _{r}*(κ) in the denominator. Incidentally, we may note that for spherical ions we have , but in general q^{eff}_{i}(κ) ≠ q_{i}*/_{r}*(κ). For a binary electrolyte of spherical ions, the deviation in κ from κ_{D} can be obtained from

which should be compared to eqn (81).

For an RPM electrolyte in the linearized PB approximation, the mean electrostatic potential due to an ion of species i, the “central” ion located at the origin, is

where κ = κ_{D} and prefactor for ψ_{i}(r) is dictated by local electroneutrality, . We can identify the effective charge as q^{eff}_{i} = q_{i}e^{κd}/(1 + κd) whereby we have ψ_{i}(r) = q^{eff}_{i}e^{−κr}/[4πε_{r}ε_{0}r]. In the PB approximation, only the central ion has q^{eff}_{i} ≠ q_{i} while all ions in its ion cloud are assumed to be bare, so for them q^{eff}_{i} = q_{i} and eqn (101) yields κ/κ_{D} = 1. However, since all ions should be treated on the same basis, they should all have q^{eff}_{i} ≠ q_{i}. If we in eqn (101) set q^{eff} equal to the value from the LPB case but with κ ≠ κ_{D}, we obtain the approximation in eqn (2).

(89) |

Let us investigate some consequences of eqn (89) and we start with real κ_{ν}. We will use the fact that

(90) |

For a spherically symmetric particle, the integral becomes

(91) |

(92) |

In the general case of nonspherical particles eqn (90) implies that

(93) |

(94) |

where

^{eff}_{i}(κ) = 〈Q^{eff}_{i}(,ω,κ)〉_{ω} | (95) |

For a pair of complex-valued screening parameters, there is an oscillatory Yukawa term in ψ_{i} with a direction dependent coefficient that can be determined from eqn (93) in a similar manner as in eqn (87).

Since the theory is valid for particles of any size and shape, the result in eqn (93) is also applicable to macroscopic particles. Therefore the decay of the potential from, for example, a planar surface is also given by the same decay modes. In such cases the effective particle charges can be translated into effective surface charge densities by dividing by the surface area.^{30}

The results in eqn (93) and (94) can alternatively be obtained in Fourier space where ψ_{i} is given by [eqn (65)]

(96) |

There exists an intimate relationship between the screening parameters κ_{ν} and the effective charges, i.e., q^{eff}_{i} for spherical and Q^{eff}_{i} for nonspherical particles. The parameter κ is a solution to (iκ) = 0 and, equivalently, to the general equation for κ, eqn (78). By inserting k = iκ into the expression (50) for (k) we obtain

(97) |

(98) |

For a spherical ion with charge q_{i} at the center we have _{i}(k) = q_{i} so we have Q_{i}(−,ω,κ) = q_{i} and Q^{eff}_{i}(,ω,κ) = q^{eff}_{i}(κ). When we deal with a system consisting solely of such ions, eqn (97) becomes [cf.eqn (69)]

(99) |

(100) |

For the special case of a binary symmetric electrolyte n^{b}_{+} = n^{b}_{−} ≡ n^{b} and q_{+} = −q_{−} ≡ q. If the spherical anions and cations differ only by the sign of their charges, as in the restricted primitive model, we have q^{eff}_{+}(κ) = −q^{eff}_{−}(κ) ≡ q^{eff} and q_{+}* = −q_{−}* ≡ q*. In this case q^{eff}(κ) = q*/_{r}*(κ) and from eqn (100) it follows that

(101) |

(102) |

We finally note that each Yukawa function term for nonspherical particles is always accompanied by terms that decay as e^{−κr}/r^{l} with l = 2, 3,… and have different orientational dependencies.^{34} The term with l = 2 has an orientational angle dependence with a combination of dipolar, quadrupolar and higher multipolar characteristics, but lacks a monopolar part, while the term with l = 3 has quadrupolar and higher multipolar orientational characteristics, but no monopolar and dipolar parts. This applies for all screening parameter values, κ, κ′, κ″ etc. These higher order terms are included in “other terms” above.

When κ → 0, the term with l = 2 goes over to a purely dipolar term that decays with distance as 1/r^{2} and the term with l = 3 goes to a purely quadrupolar term that decays as 1/r^{3}. Thereby the usual multipole expansion is obtained, which applies to the electrostatic potential from a fixed charge distribution immersed in a pure polar liquid.

provided that the screened electrostatic interactions between the two particles dominate in w

(103) |

(104) |

A very important task is to relate the decay behavior of ϕ_{Coul}*(r) to that of w_{ij} and the pair distribution functions g_{ij} = 1 + h_{ij} = e^{−βwij}. We have the expansions

(105) |

To find the connection between the decays of ϕ_{Coul}* and h_{ij}, let us introduce another pair correlation function h_{ij}*, which for spherical simple ions is defined by

(106) |

We will now show that the dressed ion charge density is given by

(107) |

where we have used eqn (54) to obtain the last equality. The rhs of this equation is equal to ρ

In the general case we define in an analogous manner

(108) |

(109) |

(110) |

(111) |

By using the definition (103) of w^{el}_{ij} we see that eqn (109) can be written as

h_{ij}(R_{1},R_{2}) = h_{ij}*(R_{1},R_{2}) − βw^{el}_{ij}(R_{1},R_{2}). | (112) |

The pair correlation function h_{ij} has accordingly been split into two parts: the function h_{ij}* of the dresses and the part −βw^{el}_{ij} that contains the screened Coulomb interaction. As we will see, the electrostatic term w^{el}_{ij} determines the decay behavior of pair correlation function h_{ij} in terms of Yukawa functions. More precisely, in the overwhelming number of cases

(i) the term −βw^{el}_{ij} in eqn (112) gives rise to all contributions to h_{ij} (and thereby to −βw_{ij}) that decays exponentially like a monotonic Yukawa function e^{−ar}/r or an oscillatory one e^{−aℜr}cos(a_{ℑ}r + ϑ)/r [with a_{ℜ} = ℜ(a) and a_{ℑ} = ℑ(a), the imaginary part] and

(ii) the decay parameter a of each such contribution satisfies (ia) = 0, which means that a is a solution to the general eqn (78) for κ, so a is a screening parameter. This implies that h_{ij}has the same set of screening parameters as ϕ_{Coul}*(r) and that the various a are equal to κ, κ′ etc. in eqn (86).

Thus, for each Yukawa term in eqn (86) there is a corresponding contribution in h_{ij} with the same screening parameters (but with different prefactor and phase shift, if any). The first term in eqn (112), the function h_{ij}*, also has contributions that decay like Yukawa functions (with other values of the decay parameters), but as shown below they do not give any contribution to h_{ij} (apart from some rare exceptions to be described later). Thus, the screened Coulomb potential and the pair correlation function normally have the same decay modes, each mode having its own values of the screening parameter, dielectric factor, effective relative dielectric permittivity and effective charges.

As we have seen, a contribution that decays like a monotonic or oscillatory Yukawa function corresponds to a simple pole in complex k-space. Let us therefore investigate the functions in Fourier space. Since the pair correlation function h_{ij}(R_{1},R_{2}) in the bulk phase depends on the separation vector r_{12} = r_{2} − r_{1} we can write h_{ij}(R_{1},R_{2}) = h_{ij}(r_{12},ω_{1},ω_{2}), so in Fourier space we obtain from eqn (103) and (112)

(113) |

(114) |

An important result of the current work is that the coupling between fluctuations in charge density and in number density, which normally takes place, makes all poles of _{ij} to be given by the zeros of (k), whereby they coincide with the poles of _{Coul}*(k). This follows from the fact, shown in Appendix C, that _{ij} and _{ij}* cannot have poles for the same k values (apart from a few exceptional cases that will be treated in Section 4.2). This result is a consequence of the fact that each pole of _{ij}* is cancelled by a corresponding pole in the last term in eqn (113), as shown explicitly in Appendix C. For binary simple electrolytes this cancellation in eqn (114) has previously been found.^{25} Since _{i}* and _{j}* are given by linear combinations of _{ij}*, their poles are also poles of _{ij}* and cannot coincide with any pole of _{ij}. Therefore, all poles of the rhs of eqn (113) [for spherical ions eqn (114)] originate from the zeros of (k) in the denominator and the assertions in (i) and (ii) above follow. The dielectric function (k) is also a linear combination of _{ij}* since _{i}* in eqn (50) is such a linear combination. The poles and zeros of (k) occur, of course, for different k values.

For spherical particles, these results and eqn (114) imply that [cf.eqn (92)]

(115) |

For nonspherical particles, an expression for h_{ij}(r_{12},ω_{1},ω_{2}) analogous to that in eqn (115) applies when r_{12} → ∞

(116) |

When κ is complex, the κ and κ′ terms combine and form an oscillatory term [cf.eqn (87)]. By writing Q^{eff}_{l}(,ω,κ) = |Q^{eff}_{l}(,ω,κ)|e^{−iγl(,ω,κ)} with a real-valued γ_{l} for l = i, j, we then have for r_{12} → ∞

(117) |

As noted in Section 3.2.4, the theory is valid for particles of any size and shape. Therefore results in eqn (116) and (117) are applicable for the correlations between, for example, a macroparticle and an ion. Thereby the potential of mean force acting on the ion as a function of distance from the macroparticle decays as w_{ij} ∼ −β^{−1}h_{ij}. Likewise, the results can be applied for the interaction between two macroparticles, for example, the surface forces between two macroscopic surfaces,^{29,30} which hence have the decay modes that are determined by the bulk electrolyte that the fluid phase between the surfaces is in equilibrium with.

The density–density, charge–charge, density–charge correlation functions, H_{NN}(r), H_{QQ}(r) and H_{NQ}(r) defined in Appendix B, have the same poles in Fourier space as _{ij}, which are in general determined by the zeros of (k). This can be realized as follows. _{NN}(k), given in Appendix B, eqn (138), is a linear combination of _{ij},

(118) |

(119) |

(120) |

As mentioned earlier there always exist terms in eqn (115) that have a different r dependence than Yukawa functions. They appear because when h_{ij}(r) and therefore w_{ij}(r) contain a term that decays as e^{−κr}/r with real κ, it follows from eqn (105) that h_{ij}(r) contains terms that decay as [e^{−κr}/r]^{2}, [e^{−κr}/r]^{3}etc. and that this also applies to w_{ij}(r). These higher order terms that appear because the system is intrinsically nonlinear have decay lengths (2κ)^{−1}, (3κ)^{−1}etc. so they decay faster than the “original” Yukawa function term, with the same κ, that has generated them.|||| Therefore they give important contributions mainly for small r. For large r some of them can, however, dominate over Yukawa function terms with larger κ values in eqn (115), for example the term with e^{−κ′r}/r provided κ′ > 2κ. There also exist terms that decay as ζ(r)[e^{−κr}/r]^{l}, where ζ(r) is a slowly varying function and l ≥ 2.^{21} Furthermore, there are cross-terms from two or more Yukawa functions with different decay lengths. All these higher order terms are included in “other terms” in eqn (115) and it is not worthwhile to consider more than the leading ones individually and explicitly. For complex-valued κ, similar conclusions are valid for exponentially decaying oscillatory terms and the same applies in eqn (116) for nonspherical particles.

All of these contributions are generated by the set of fundamental decay modes with screening parameters κ_{ν} that are solutions to the general eqn (78) for κ and, equivalently, solutions to (iκ_{ν}) = 0 and to eqn (97). Thus, the decay behaviors of both the screened electrostatic potential and the correlation functions originate from fundamental decay modes.

In cases where the non-electrostatic pair interactions u^{ne}_{ij}(r) do not have a short range (contrary to our assumptions earlier), but instead have a power law decay like the r^{−6} dispersion interactions, the functions h_{ij}(r), w_{ij}(r) and ϕ_{Coul}*(r) ultimately decay like a power law when r → ∞.^{35} These functions still have Yukawa function terms*** that are given by the present formalism and the power law terms are included in “other terms” in the various equations. The Yukawa terms can give dominant contributions for short to intermediate r values, but they can never dominate for very large r because they decay faster than any power law. In many cases the Yukawa function terms have a dominant influence in h_{ij}(r) for most r. This is, for example, seen in the simulation results by Keblinski et al.^{24} mentioned in the Introduction for the realistic model for NaCl that includes dispersion interactions. Their calculations show that the monotonic and oscillatory Yukawa function terms give the dominant contributions.

The other exceptional case is a category of systems where _{ij}*, in contrast to the main case, always gives Yukawa function contributions to h_{ij}, namely model systems that are invariant when we invert the sign of all charges of the particles. Such systems, which we will call charge-inversion invariant systems, remain exactly the same when one does such a charge inversion, whereby the internal charge distribution σ_{i}(r|R_{1}) for each particle changes to −σ_{i}(r|R_{1}) (positive regions become negative and vice versa without change in the absolute value of σ_{i} for each r). This means, for example, that the anions become cations and vice versa during charge inversion. For an invariant system, for each cation species there must exist an anion species that is identical in all respects apart from the sign of σ_{i}(r|R_{1}), like the same size, same shape and same nonelectrostatic interactions with other particles, and have the same number density. Examples of such systems include the restricted primitive model, where the anions and cations of the same absolute valency are charged hard spheres of equal size. As soon as there is any difference, however small, between anions and cations apart from the sign of their charges, the main result applies, so all Yukawa function terms in h_{ij} originate from the last term in eqn (113) and all poles of _{ij} arises from the zeros of (k). For charge-inversion invariant electrolyte systems, all electroneutral species present must also satisfy invariance conditions. Examples include electrolyte models with explicit solvent where the solvent molecules turn into themselves during a charge inversion, for instance spherical particles with a dipole at the center.

The reason why charge-inversion invariant systems are exceptions is that the extreme symmetry forces the density–charge correlation function H_{NQ}(r) to be identically zero, so fluctuations in charge density and in number density are uncoupled from each other. It is simple to realize that H_{NQ}(r) must be identically equal to zero in such systems. Suppose that H_{NQ}(r) is, say, positive for a certain r value and we invert all charges in the system, H_{NQ}(r) would become negative but since the system is charge-inversion invariant H_{NQ}(r) must remain the same, which implies that H_{NQ}(r) must be zero. Alternatively, this can be realized when we consider correlations between fluctuations in density and fluctuations in charge at two points separated by distance r, because for a given fluctuation in density, the probability for positive and negative fluctuations in charge must be equal, so the fluctuations will average to zero. Mathematically, the fact that H_{NQ}(r) ≡ 0 can be realized from the rhs of the definition of H_{NQ}(r) in Appendix B [eqn (139)]. For each positive contribution on the rhs there must exist an equally large negative contribution due to the charge-inversion invariance. This can likewise be realized from eqn (140). Exactly the same argument applies to eqn (120) because during a charge inversion _{i}*(k,ω_{1}) for each particle changes to −_{i}*(k,ω_{1}), so we must have

(121) |

Let us now consider the density–density correlation function for charge-inversion invariant systems. By inserting _{ij} from eqn (113) into eqn (118) we see that the contribution from the last terms in eqn (113) cancels identically in _{NN}(k) because of eqn (121). Therefore we have

(122) |

Consider two species of ions that swap their identities as anions and cations during a charge inversion. Let us call them A^{+} and A^{−}, so we have σ_{A+}(r|R_{1}) = −σ_{A−}(r|R_{1}). The common notation A indicates that they are identical apart from the sign of their internal charge distributions. We likewise consider two other species of ions B^{+} and B^{−} with σ_{B+}(r|R_{1}) = −σ_{B−}(r|R_{1}). They are also identical to each other apart from the sign of σ_{i}. The charge-inversion invariance implies that we have the symmetry

We now define h_{S,AB}(r,ω_{1},ω_{2}), h_{D,AB}(r,ω_{1},ω_{2}) and the corresponding h* functions from

(123) |

Due to the charge-inversion invariance we have _{A+}* = −_{A−}* ≡ _{A}*, which defines _{A}*, and likewise for _{B}*. Therefore we have from eqn (113)

_{S,AB}(k,ω_{1},ω_{2}) = _{S,AB}*(k,ω_{1},ω_{2}) | (124) |

(125) |

In fact, H_{NN}(r) is a linear combination of the S-functions and H_{QQ}(r) is a linear combination of the D-functions (including those for any electroneutral particles, if present). Furthermore, the total particle density around each of the particles in charge-inversion invariant systems is determined by S-functions and the total charge density ρ^{tot}_{i} is determined by D-functions, so these two entities have different decay parameters.

All these facts are well-known for the restricted primitive model where there are only two species A^{+} and A^{−} in a dielectric continuum solvent, but the new findings here are that this applies in general for charge-inversion invariant systems with any number of components including solvent molecules and other electroneutral particles. All particles can have any shape, size, internal charge density and nonelectrostatic interactions subject to the conditions of charge-inversion invariance.

In the RPM, where we have , it can happen that h_{S,AA}(r) has a longer decay length than h_{D,AA}(r) when the ion density is high.^{19,20,39} Then, the tail of h_{ij}(r) for large r has the same sign irrespective of the signs of the i and j-ions. This has been denoted as “core dominance,”^{39} because h_{S,AA}(r) is decoupled from electrostatics and is mainly determined by the hard core packing of the ions. For low ion densities, where h_{D,AA}(r) has the longest decay length, the tails of h_{++}(r) and h_{+−}(r) have different signs and there is “electrostatic dominance.” Such strict distinction between core dominance and electrostatic dominance exists only in charge-inversion invariant systems.

It is interesting to consider systems that are close to being charge-inversion invariant. This is the case, for example, in the primitive model of binary symmetric electrolyte solutions when anions and cations have equal absolute valency but differ slightly in size. Then, there is no longer a decoupling between electrostatic and nonelectrostatic correlations, so there is electrostatic coupling in the decay modes where hard core packing of the ions dominates and vice versa. Systems close to charge-inversion invariance also occur in more realistic models of symmetric electrolytes where the anions and cations have nearly the same nonelectrostatic pair interactions with their own species u^{ne}_{++}(r) ≈ u^{ne}_{−−}(r), as in the model for NaCl mentioned in the Introduction. Another example is a model with explicit solvent where anions and cations are identical apart from their signs but the solvent molecules are modeled as spheres with a radially aligned dipole that lies somewhat off-center.

For these systems all poles of _{ij} are given by the zeros of (k), while for charge-inversion invariant systems the decay parameter values of h_{ij} belong, as we have seen, to two distinct sets: those that are due to poles of _{ij}* and those that are due to zeros of (k). Since the two kinds of systems differ very little from each other, their decay parameters must be closely related and vary continuously when one changes a system from being nearly charge-inversion invariant to being invariant or vice versa. This must apply not only for the decay parameters that are given by zeros of (k) all the time, but also for those that initially are due to such zeros and then, for the invariant system, belong to the set that are due to poles of _{ij}* but not zeros of (k). In Appendix C it is shown why and how this happens. One can visualize the findings by following the “trajectories” of the poles of _{ij} and _{ij}* as functions of the system parameters like ion sizes etc. For a system that is nearly charge-inversion invariant, it is found in the appendix that there exist one set of poles of _{ij} [zeros of (k)] where each pole is close to a pole of _{ij}* and where the trajectories of the poles of _{ij} and _{ij}* cross each other at the point where the system become invariant. At the crossing point the pole of _{ij} ceases to be a zero of (k) and becomes instead a pole of both _{ij}* and _{ij}. The second set of poles of _{ij} do not have such crossings and these poles remain zeros of (k) throughout. These two sets correspond to the sets for the invariant system introduced above.

More precisely, it is shown in the appendix that there exist zeros of (k) for the nearly invariant system that lie close to the poles of _{ij}* for the same system. The latter are not poles of _{ij} while the zeros of (k) are such poles. Each of these zeros moves continuously with the system parameters so that when the system is turned into a charge-inversion invariant one, it merges with the corresponding pole of _{ij}* and ceases to be a zero of (k) but remains as a pole of _{ij}.

Let us consider cases where the leading term in h_{ij} is given by the monotonic Yukawa term shown explicitly in eqn (116). This term originates from the last term in eqn (113) evaluated at the leading zero of (k), i.e., k = iκ and we assume that κ is real. We will investigate the leading term when r → ∞ for H_{NN}(r), H_{QQ}(r) and H_{NQ}(r), respectively. The decay length is 1/κ for all three functions and we will find that the magnitudes of the leading term of these functions are interdependent via common prefactors, which can be obtained from the effective multipole charges Q^{eff}_{i}. We will also show that the intimate relationship between the screening parameter κ and Q^{eff}_{i} found in Section 3.2.4 has direct consequences for the relative importance of the different correlation functions.

These results can be obtained from eqn (118)–(120) together with eqn (116) and we obtain for r → ∞

where we can select = [cf.eqn (94) and (96)] and where Q

(126) |

(127) |

(128) |

(129) |

H
_{NN}(r), H_{NQ}(r) and H_{QQ}(r) have the same decay length, but depending on the values of the system parameters the leading term of H_{QQ}(r) can be larger than that of H_{NN}(r) or vice versa. Since κ^{2} = βn^{b}_{tot}q_{e}^{eff}_{Q}(κ)/ε_{0}, the decay length 1/κ is directly linked to the magnitude of H_{QQ}(r), while the magnitude of H_{NN}(r) does not have such a direct link. For a situation with long-range density–density correlations, the screening parameter κ is small and the charge–charge correlations must be small and become even smaller when the decay length increases. This is, for example, relevant when a critical point is approached, whereby the charge–charge correlations become less and less significant although they have the same decay length as the density–density correlations. In the present analysis we do, however, avoid the immediate neighborhood of critical points, where other considerations have to be made. We may, however, note that in the limit of 1/κ → ∞ the charge–charge correlations vanish as they must do.

For charge-inversion invariant systems ^{eff}_{N} ≡ 0 and H_{NQ}(r) ≡ 0, while the decay of H_{NN}(r) can be obtained from

Returning to the general case, the conclusions above expressed in eqn (127)–(129) are valid for the leading term in H_{NN}(r), H_{NQ}(r) and H_{QQ}(r) when κ is real. Analogous conclusions are valid for the Yukawa terms with other κ_{ν} values in these correlation functions, so similar relationships to these equation are valid for all decay modes (including the oscillatory case where there is also a cosine factor). The relative magnitudes of |^{eff}_{N}| and |^{eff}_{Q}| varies for the different modes, so the contributions with some decay lengths may have |^{eff}_{N}| > |^{eff}_{Q}| while the reverse can be true for those with other decay lengths. In this regard it is particularly illustrative to consider systems that are close to being charge-inversion invariant, but let us start with the corresponding invariant systems.

For charge-inversion invariant systems we have seen that there are two distinct sets of poles for _{ij}: the poles of _{NN}(k) [poles of _{ij}*] and those of _{QQ}(k) [zeros of (k)]. This implies that the Yukawa function terms in H_{NN}(r) and H_{QQ}(r) have distinct decay parameter values; a Yukawa term that is present in H_{NN}(r) is absent in H_{QQ}(r) and vice versa. A system that is close to being charge-inversion invariant has, however, all such terms present in both H_{NN}(r) and H_{QQ}(r). These two functions then have the same decay parameters and, as we saw at the end of Section 4.2, when one breaks the charge-inversion invariance all decay parameters vary continuously with the system parameters. We also saw that the two sets of poles for _{ij} remain in the sense that each pole in one set is close to a pole of _{ij}* (and merge with it for the invariant system) and the other set consists of poles that do not behave in this manner and are poles of (k) throughout. For continuity reasons, the magnitudes of the Yukawa terms in H_{QQ}(r) with decay parameters in the first set must be close to zero (they are exactly zero for the invariant system), while the magnitudes of the terms in H_{NN}(r) belonging to the second set must be close to zero for the same reason. Thus |^{eff}_{N}| ≫ |^{eff}_{Q}| for the first set and |^{eff}_{N}| ≪ |^{eff}_{Q}| for the second. This is the situation for the NaCl system mentioned in the Introduction.

Finally, we will turn to dilute electrolyte solutions with molecular solvent. The leading term in h_{ij} for large distances then has a decay parameter κ that approaches the Debye parameter κ_{D} at high dilution, as we saw in eqn (85). In a dilute solution where κ is small we have ^{eff}_{i}(κ) ≈ q_{i}, which implies that n^{b}_{tot}^{eff}_{N}(κ) is very small because due to electroneutrality (q_{solvent} = 0 since the solvent molecules are uncharged). Thus, the leading term in H_{NN}(r), which is proportional to [n^{b}_{tot}^{eff}_{N}(κ)]^{2}, is very small. H_{QQ}(r) is much larger because its coefficient is proportional to the square of

As in other systems, there exist other decay modes with terms of different decay lengths in the correlation functions, for example contributions where the density–density fluctuations of the solvent are prominent. In contrast to the very small terms in H_{NN}(r) with decay length 1/κ that we investigated above, these contributions to H_{NN}(r) are accordingly substantial. Such solvent-dominated contributions are given by Yukawa functions with screening parameters κ_{ν} that satisfy eqn (78) and each mode yields terms in H_{NN}(r), H_{NQ}(r) and H_{QQ}(r) with the same decay length. As we saw from the examples in the Introduction, for dilute electrolyte solutions at least one pair of these screening parameters is complex-valued and gives rise to an oscillatory component of the correlation functions that reflects the structure of the molecular solvent.

In the limit of zero electrolyte concentration, the screened Coulomb potential decays as ϕ_{Coul}*(r) ∼ 1/(4πε_{0}ε_{r}r) which originates from the leading term e^{−κr}/(4πε_{0}^{eff}_{r}(κ)r) when κ → 0. However, the solvent-dominated decay modes make ϕ_{Coul}*(r) oscillatory for small to intermediate r values as discussed in more detail for aqueous systems in ref. 29. These oscillations dominate in ϕ_{Coul}*(r) for pure water up to quite large r values where 1/r tail eventually takes over.

At higher concentrations of electrolyte, one cannot clearly distinguish between contributions to the correlation functions that are solvent-dominated and those that are electrolyte-dominated because there is an intricate coupling between the different constituents of the solution. In all cases there are several decay modes with different screening parameters that are simultaneously present.

Perhaps the most important finding in this work is the fact that all decay modes in electrolytes, with the exception of charge-inversion invariant electrolyte systems, are governed by the dielectric response and that the decay lengths hence are determined by the parameters κ_{ν} that are solutions to eqn (5), which are the screening parameters for the electrostatic potential. Accordingly, all decay modes of the pair potential of mean force w_{ij} for the particles are the same as for the screened Coulomb potential ϕ_{Coul}*(r); the latter modes are shown explicitly in eqn (86). This applies for both the monotonic and the oscillatory exponential modes, including those that are dominated by “packing” of molecules in dense liquids. It also applies to the long-range monotonic modes at high ionic densities and the long-range density–density correlations on approach of criticality (but not at the critical point itself). All these facts are a consequence of the coupling between fluctuations in charge density and fluctuations in number density.

Another very significant result is that each decay mode has its own value of the effective relative dielectric permittivity ^{eff}_{r}(κ_{ν}), which determines the magnitude of the mode's contribution to ϕ_{Coul}*(r) [eqn (86)] and, for oscillatory modes, determines also the phase. ^{eff}_{r}(κ_{ν}) is the physical entity for electrostatic interactions in electrolytes that corresponds to the dielectric constant ε_{r} for pure polar fluids and other nonelectrolytes. In contrast to ε_{r} it is not given by the infinite wavelength dielectric response (k → ∞) and does not solely contain the dipolar orientational dielectric polarization. As discussed in more detail in ref. 31, the values of the effective dielectric permittivities result from all kinds of polarizations: orientational polarizations (dipolar, quadrupolar, octupolar etc.) and from changes in ion distributions, including transient aggregate formations like ion pairing. The value of ^{eff}_{r}(κ_{ν}) reflects the polarization response of the electrolyte to an exponentially decaying disturbance with decay parameter κ_{ν}.

Furthermore, the dielectric factor _{r}*(κ_{ν}) that appears in the expression for the screening parameter κ_{ν}, eqn (5), has different values for the various modes. As we have seen, _{r}*(κ_{ν}) replaces the dielectric constant ε_{r} that appears in the expression for the Debye parameter κ_{D} [eqn (1)]. It is not correct to use a dielectric permittivity obtained from the dielectric polarization response at infinite wavelength instead of ^{eff}_{r}(κ_{ν}) and _{r}*(κ_{ν}), except possibly as an approximation for modes with long decay lengths. Thus, when κ_{ν} is not close to zero, one cannot use experimental “dielectric constants” for electrolytes evaluated macroscopically at k = 0 from measurements at low frequencies. Furthermore, theoretical estimates of such dielectric constants based on dipolar polarizations including the effect of ions on dipolar orientations are clearly inadequate, since they leave out the polarization from changes in ion distributions and multipole orientations.

It is shown in Section 4.1 that the decay modes of h_{ij} and w_{ij} are determined by the screened electrostatic pair interaction w^{el}_{ij}, which is a part of w_{ij}. The entity w^{el}_{ij} is equal to the electrostatic interaction energy of the dressed particles of species i and j (with charge densities ρ_{i}* and ρ_{j}*) as mediated by the screened Coulomb potential ϕ_{Coul}*(r) [eqn (103) and (104)], that is, w^{el}_{ij} = ρ_{i}* ⊛ ϕ_{Coul}* ⊛ ρ_{j}*, where ⊛ denotes a convolution integral. Alternatively this can be expressed as follows: w^{el}_{ij} is equal to the interaction energy between the mean electrostatic potential ψ_{i} due to the particle of species i and the charge distribution ρ_{j}* of a dressed particle of species j (or vice versa). Thereby we have used the fact that ψ_{i} = ρ_{i}* ⊛ ϕ_{Coul}*, that is, ψ_{i} is obtained via Coulomb's law for the screened potential, eqn (60), where ρ_{i}* is the source for ψ_{i}.

The screened Coulomb potential ϕ_{Coul}*(r) is a Green's function that describes the spatial propagation of electrostatic interactions in electrolytes. It is closely related to the dielectric function (k); in Fourier space we have _{Coul}*(k) = _{Coul}(k)/(k), where _{Coul}(k) = 1/(ε_{0}k^{2}) is the Fourier transform of the ordinary unscreened Coulomb potential ϕ_{Coul}(r). The decay modes of ϕ_{Coul}*(r) are given by the poles of _{Coul}*(k) in complex k-space, i.e., the zeros of (k), whereby the screening parameter κ satisfies (iκ) = 0, where the imaginary unit i appears because the dielectric response is evaluated for an exponentially decaying field. The condition (iκ) = 0 is equivalent to the general exact equation for the screening parameter, eqn (5).

In analyses of experimental measurements of interparticle interactions in electrolytes, including surface force measurements, it is important to realize the existence of the various decay modes with different values of the screening parameter κ_{ν}. It is therefore suitable to analyze the data using several decay modes with various decay lengths and, when appropriate, wavelengths. Thereby, the various modes may dominate in different distance intervals.

In both experimental and theoretical work, one can make curve fits to force curves and similar data plots on a log scale in order to find the decay lengths and the wavelengths. In cases where the decay modes have quite different decay lengths, it should be feasible to identify the modes, but if, say, two modes have nearly the same decay lengths, it is in general difficult to distinguish them. It can also be difficult to know if one has determined the functions for sufficiently large distances so the results cover the ultimate decay where the mode with the largest decay length dominates. Furthermore, if the magnitude of a mode is small it may not be possible to detect it.

In theoretical work, including simulations, one should, if practically feasible, solve eqn (5) for the decay parameter or the equivalent equation (iκ) = 0 in order to determine the various κ_{ν}. One can then also determine modes that are difficult to detect by fitting. For spherical ions, eqn (99) is an option. The solutions should then be compared with the decay parameters obtained by fitting. Thereby, one has two independent ways to determine the decay modes. Ref. 18 gives examples of how this can be done in practice in simulations of electrolytes with spherical ions.

For binary electrolytes, in particular, it can be useful to plot κ_{ν}/κ_{D} as a function of, for example, temperature or concentration since this ratio has a simple relationship to q_{i}* and _{r}*(κ_{ν}), see eqn (81). For spherical ions there is also a simple relationship to the effective ionic charges, see eqn (100).

A central theme of the present work is the demonstration of the fact that when one deals with the screened electrostatic interactions in electrolytes, it is appropriate and often advantageous to use the concept of a dressed particle instead of the entity composed of the particle together with the entire ion/solvent cloud that surrounds it. This cloud consists of the charge distribution ρ_{i} due to excess ions and solvent molecules close to the particle, including the effects of their orientational ordering due to interactions. A dressed particle is the particle together with its dress, which has a charge distribution ρ^{dress}_{i} that is different from ρ_{i}. The total charge density of the particle and its surrounding cloud is ρ^{tot}_{i} = σ_{i} + ρ_{i}, where σ_{i} is the internal charge density of the particle. Likewise, we have ρ_{i}* = σ_{i} + ρ^{dress}_{i}. While the ion/solvent cloud ρ_{i} is described by the pair distribution functions g_{ij} = h_{ij} + 1, the dress ρ^{dress}_{i} is described by the distribution functions g_{ij}* = h_{ij}* + 1 that are related to the former by the simple relationship g_{ij}* = g_{ij} + βw^{el}_{ij} [eqn (112)]. The dressed particle charge q_{i}* that occurs in eqn (5) for the screening parameter is the total charge of ρ_{i}*.

The reason why it is appropriate to use dressed particles is apparent already when considering the linear polarization response of the bulk electrolyte due to a weak external electrostatic potential δΨ^{ext}, as discussed in Section 2. This response can be described microscopically from the point of view of the free energy of interaction for each constituent particle in the fluid (ion, solvent molecule or any other particle). The interactional free energy δW_{i} for a particle of species i is simply equal to the electrostatic interaction energy δW_{i} = ρ^{tot}_{i} ⊛ δΨ^{ext} between δΨ^{ext} and ρ^{tot}_{i} [eqn (16) and (23)]. However, δΨ^{ext} is the potential from the external charges in the absence of the electrolyte, while the total potential δΨ, which also includes the potential δΨ^{pol} from the induced polarization charge density, is the actual potential in the presence of the electrolyte. The latter gives the real situation for the particles in the system and, as shown in Section 3, the free energy δW_{i} is equal to the electrostatic interaction δW_{i} = ρ_{i}* ⊛ δΨ between δΨ and ρ_{i}* [eqn (41)]. Thus, the use of dressed particles reflects the actual conditions for the particles in the weakly polarized electrolyte.

The strong, nonlinear polarization in the immediate neighborhood of the i-particle, which is due to the interactions between the particle and its surroundings, is included in the dress. The linear part of the electrostatic polarization due to the particle, ρ^{lin}_{i}, is included in ρ^{tot}_{i} but not in ρ_{i}*, and we have ρ^{tot}_{i} = ρ_{i}* + ρ^{lin}_{i} [eqn (44)]. The contribution from this linear part to the potential of mean force δW_{i} = ρ_{i}* ⊛ δΨ is instead included in the factor δΨ via the collective response to the external potential from all particles in the electrolyte.

A similar decomposition in linear and nonlinear parts of the electrostatic polarization due to a particle is done in the calculation of ψ_{i}via Coulomb's law. The mean potential ψ_{i} can be obtained from ρ^{tot}_{i} by using the usual unscreened Coulomb potential, ϕ_{Coul}(r) = 1/(4πε_{0}r) in this law [eqn (38)]. As mentioned above, exactly the same potential ψ_{i} can alternatively be obtained from ρ_{i}* by using ϕ_{Coul}*(r) in Coulomb's law [eqn (60)]. The difference is that the linear part of the electrostatic polarization in the latter case is included via ϕ_{Coul}*(r), which contains an electrostatic linear response via (k), while in the former case the linear part is included in ρ^{tot}_{i}.

The decay modes of ψ_{i} are always the same as those of ϕ_{Coul}*(r). The magnitude of each of these decay modes of ψ_{i} is proportional to a kind of effective charge Q^{eff}_{i} of the dressed i-particle, which has different values for the different modes, and inversely proportional to ^{eff}_{r} [eqn (93)]. The values of Q^{eff}_{i} are determined by the projections of ρ_{i}* on the various modes, which define Q^{eff}_{i} for each mode [eqn (94)]. For the oscillatory modes there is also a phase shift determined by ρ_{i}*. In this case, Q^{eff}_{i} is complex-valued and its absolute value gives the magnitude of the mode and its phase gives the phase shift. It is clear that the concept of dressed particle has a fundamental role since ρ_{i}* directly determines the magnitudes and the phase shifts of the modes of ψ_{i}. For a nonspherical particle, the effective charge and the phase shift are direction dependent quantities which reflect the anisotropy of the mean potential ψ_{i}. The former is therefore denoted as the “multipolar effective charge” [Section 3.2.4].

The fundamental role of the dressed particles is further accentuated from their appearance in w^{el}_{ij} as mentioned earlier and the fact that w^{el}_{ij} in general determines all decay modes of h_{ij} and w_{ij}. The magnitude of each decay mode of h_{ij} and w_{ij} is proportional to the product Q^{eff}_{i}Q^{eff}_{j} as evaluated for each mode [eqn (116)]. For oscillatory modes the phase shift contains the sum of the phases of Q^{eff}_{i} and Q^{eff}_{j} [eqn (117)].

The correlation functions H_{NN}(r), H_{NQ}(r) and H_{QQ}(r) have in general the same set of decay modes as w^{el}_{ij} and all three functions therefore have the same decay lengths. Their magnitudes are proportional to a set of prefactors ^{eff}_{N} and ^{eff}_{Q} that constitute averages of the effective charges of the particles in the electrolyte [eqn (126)]. The modes of H_{NN}(r) are proportional to (^{eff}_{N})^{2}, those of H_{NQ}(r) proportional to ^{eff}_{N}^{eff}_{Q} and those of H_{QQ}(r) proportional to (^{eff}_{Q})^{2} [eqn (127)–(129)].

For charge-inversion invariant systems, where H_{NQ}(r) is identically equal to zero, the decay modes of H_{QQ}(r) are determined by w^{el}_{ij}, but those for H_{NN}(r) are different. The latter are instead determined by the correlation functions h_{ij}* of the dresses [eqn (122)], so the concept of dressed particles has in this case yet another fundamental role. The decay modes of H_{NN}(r) are here given by the poles of _{ij}*, which are different from the poles of _{Coul}*(k). Only in charge-inversion invariant systems there is a strict distinction between “electrostatic dominance” and so-called “core dominance” (in the latter, packing effects dominate), for example in the restricted primitive model.^{39} The former kind of dominance occurs when H_{QQ}(r) has the longest decay length and the latter when H_{NN}(r) has the longest one.

There exist electrolytes that are nearly charge-inversion invariant, for example the NaCl system mentioned in the Introduction and any binary system in the primitive model that have anions and cations of the same valency but with slightly different sizes. For such systems, the conclusions for the main case apply and h_{ij}, w_{ij}, ϕ_{Coul}*(r), H_{NN}(r), H_{QQ}(r) and H_{NQ}(r) have the same set of decay modes. In this case the poles of _{Coul}*(k) [i.e., zeros of (k)] can be divided into two categories depending on what takes place when the system is converted into one that is charge-inversion invariant; in the given primitive model example this happens when the ion sizes are made equal. Each pole in the first category lies close to a pole of _{ij}* and the two poles move towards each other when the system is converted into an invariant one. When the system becomes invariant the pole of _{Coul}*(k) vanishes. The decay mode remains but is instead determined by the pole of _{ij}* for the invariant case. Each pole of _{Coul}*(k) in the second set continues to be such a pole throughout. For the decay modes given by first set of poles |^{eff}_{N}| ≫ |^{eff}_{Q}| and ^{eff}_{Q} becomes zero when the system becomes charge-inversion invariant, while for the second set |^{eff}_{N}| ≪ |^{eff}_{Q}| and ^{eff}_{N} becomes zero for the invariant case. Thereby, the various modes are smoothly changed when the electrolyte is turned into a charge-inversion invariant one.

and for a bulk fluid perturbed by the weak external potential δv

When the external potential is solely due to the electrostatic potential δΨ

and hence we obtain

(130) |

(131) |

where ρ

(132) |

(133) |

(134) |

(135) |

(136) |

The density–density correlation function H_{NN} is defined as

(137) |

(138) |

Finally, the density–charge correlation function H_{NQ} ≡ H_{QN} is defined as

(139) |

(140) |

The existence of a simple pole for _{ij}* at k = ±iα means that _{ij}* diverges like

(141) |

when r → ∞.

The factorization in eqn (141) is crucial for the cancellation of poles as we now will see. The Fourier transform of eqn (111) is given by the first line in the following equation and by inserting eqn (141) we obtain the second line

(142) |

(143) |

(144) |

We now insert eqn (142) and (144) into the last term in eqn (113) and since _{i}*, _{j}* and (k) are dominated by the diverging contributions for k values close to iα, it follows that we have when k → iα

Charge-inversion invariant systems, which are treated in Section 4.2, are, however, different. For such systems E(k,α) is identically equal to zero because for each positive contribution on the rhs of eqn (143) there is an equally large negative contribution. Therefore, the pole of _{ij}* at k = iα is not cancelled in _{ij}. For other systems it may happen at exceptional points in the system's parameter space that E(k,α)|_{k=iα} fortuitously happens to be equal to zero, which means that the cancellation may not take place at such points. Otherwise the cancellation always occurs and the poles of _{ij} are given by the zeros of (k).

For systems that are close to being charge-inversion invariant, E(k,α) is small but nonzero. Such systems should have properties that are very similar to the corresponding charge-inversion invariant system. For the latter, it is shown in Section 4.2 that the decay parameter values of h_{ij} belong to two groups: those that are poles of _{ij}* and those that zeros of (k). When the system deviates only slightly from being charge-inversion invariant and all poles of _{ij} are given by the zeros of (k), the former group must have turned into zeros of (k) because the poles of _{ij}* are not poles of _{ij}. This is indeed the case as we now are going to see.

When E(k,α) is small but nonzero and _{ij}* has a pole at k = iα, the function _{ij} does not have a pole there, but instead there is, in fact, a pole of _{ij} close to this k value, that is, there exist a zero of (k) close to k = iα, say, at k = iκ_{ν}. This can be realized as follows. Since (iκ_{ν}) = 0 we obtain from eqn (144)

(145) |

(146) |

- M. A. Gebbie, M. Valtiner, X. Banquy, E. T. Fox, W. A. Henderson and J. N. Israelachvili, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 9674 CrossRef CAS PubMed.
- M. A. Gebbie, H. A. Dobbs, M. Valtiner and J. N. Israelachvili, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 7432 CrossRef CAS PubMed.
- R. M. Espinosa-Marzal, A. Arcifa, A. Rossi and N. D. Spencer, J. Phys. Chem. Lett., 2014, 5, 179 CrossRef CAS PubMed.
- T. Baimpos, B. R. Shrestha, S. Raman and M. Valtiner, Langmuir, 2014, 30, 4322 CrossRef CAS PubMed.
- A. M. Smith, A. A. Lee and S. Perkin, J. Phys. Chem. Lett., 2016, 7, 2157 CrossRef CAS PubMed.
- 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 RSC.
- A. A. Lee, C. S. Perez-Martinez, A. M. Smith and S. Perkin, Faraday Discuss., 2017, 199, 239 RSC.
- N. Hjalmarsson, R. Atkin and M. W. Rutland, Chem. Commun., 2017, 53, 647 RSC.
- A. M. Smith, A. A. Lee and S. Perkin, Phys. Rev. Lett., 2017, 118, 096002 CrossRef PubMed.
- A. A. Lee, C. S. Perez-Martinez, A. M. Smith and S. Perkin, Phys. Rev. Lett., 2017, 119, 026002 CrossRef PubMed.
- D. G. Hall, Z. Phys. Chem., 1991, 174, 89 CrossRef CAS.
- R. Kjellander, J. Phys. Chem., 1995, 99, 10392 CrossRef CAS.
- R. J. F. Leote de Carvalho and R. Evans, Mol. Phys., 1994, 83, 619 CrossRef.
- J. G. Kirkwood, Chem. Rev., 1936, 19, 275 CrossRef CAS.
- J. G. Kirkwood and J. C. Poirier, J. Phys. Chem., 1954, 58, 591 CrossRef CAS.
- C. W. Outhwaite, Chem. Phys. Lett., 1970, 5, 77 CrossRef CAS.
- C. W. Outhwaite, in Statistical Mechanics: A Specialist Periodical Report, ed. K. Singer, The Chemical Society, London, 1975, vol. 2, p. 188 Search PubMed.
- J. Ulander and R. Kjellander, J. Chem. Phys., 2001, 114, 4893 CrossRef CAS.
- J. Ennis, PhD thesis, Australian National University, 1993.
- J. Ennis, R. Kjellander and D. J. Mitchell, J. Chem. Phys., 1995, 102, 975 CrossRef CAS.
- R. Kjellander and D. J. Mitchell, J. Chem. Phys., 1994, 101, 603 CrossRef CAS.
- D. J. Mitchell and B. W. Ninham, Chem. Phys. Lett., 1978, 53, 397 CrossRef CAS.
- A. M. Smith, P. Maroni, G. Trefalt and M. Borkovec, J. Phys. Chem. B, 2019, 123, 1733 CrossRef CAS PubMed.
- P. Keblinski, J. Eggebrecht, D. Wolf and S. R. Phillpot, J. Chem. Phys., 2000, 113, 282 CrossRef CAS.
- J. Ulander and R. Kjellander, J. Chem. Phys., 1998, 109, 9508 CrossRef CAS.
- G. Stell, J. Stat. Phys., 1995, 78, 197 CrossRef.
- F. Coupette, A. A. Lee and A. Härtel, Phys. Rev. Lett., 2018, 121, 075501 CrossRef CAS PubMed.
- ESI for ref. 27 at http://link.aps.org/supplemental/10.1103/PhysRevLett.121.075501.
- R. Kjellander, J. Chem. Phys., 2018, 148, 193701 CrossRef PubMed.
- R. Kjellander, Phys. Chem. Chem. Phys., 2016, 18, 18985 RSC.
- R. Kjellander, J. Chem. Phys., 2016, 145, 124503 CrossRef PubMed.
- R. Kjellander and D. J. Mitchell, Chem. Phys. Lett., 1992, 200, 76 CrossRef CAS.
- R. Ramirez and R. Kjellander, J. Chem. Phys., 2003, 119, 11380 CrossRef CAS.
- R. Ramirez and R. Kjellander, J. Chem. Phys., 2006, 125, 144110 CrossRef PubMed.
- R. Kjellander and R. Ramirez, J. Phys.: Condens. Matter, 2008, 20, 494209 CrossRef.
- J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids, Academic, Amsterdam, 2013 Search PubMed.
- M. Ohba and K. Arakawa, J. Phys. Soc. Jpn., 1981, 50, 743 CrossRef CAS.
- A. Chandra and B. Bagchi, J. Chem. Phys., 1989, 91, 3056 CrossRef.
- P. Attard, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1993, 48, 3604 CrossRef CAS.

## Footnotes |

† Dedicated to the memory of Per Linse, Lund University, Sweden. |

‡ As shown in Section 3.2.4, for a binary symmetric electrolyte with ions of equal diameters we have the exact relationship κ^{2}/κ_{D}^{2} = q^{eff}/q [eqn (101)], where q^{eff} is an effective ionic charge and q = q_{+} = −q_{−}. Using the LPB approximation as a guide for an approximative value of q^{eff} for all ions (i.e., not only for the ion at the origin), we set q^{eff} = qe^{κd}/(1 + κd) and obtain eqn (2); for details see Section 3.2.4. |

§ The general limiting law at high dilution for the decay parameter κ of the leading decay mode is from ref. 21 where z |

¶ In the work by Keblinski et al. they use the notation Q(r) and G(r), where Q(r) is proportional to H_{QQ}(r) and the function G(r) − 4n^{b}_{tot} is proportional to H_{NN}(r) [n^{b}_{tot} is the total ion density]. The definitions of H_{QQ}(r) and H_{NN}(r) can be found in Appendix B. |

|| We can, for example, take where (φ,θ,η) are the Euler angles of the particle or, for a linear molecule, since η is redundant in the latter case. |

** The formalism has a general validity, but in the presence of nonelectrostatic interactions that decay like a power law (like dispersion interactions) there are terms that are not explicitly included in the asymptotic expressions for large distances obtained in this work; see the comments at the end of Section 4.1. |

†† In macroscopic electrostatics we have P = ε_{0}χ_{e}E in ordinary space, while in the present microscopic case for bulk fluids we have the corresponding relationship in Fourier space. |

‡‡ This can be shown by contour integration and residue calculus in complex k-space. |

§§ This is a consequence of the Stillinger–Lovett second moment condition. |

¶¶ This can be shown by contour integration and residue calculus for _{Coul}*(k) in complex k-space. |

|||| The higher order terms give singularities in complex Fourier space that are different from simple poles, for instance e^{−2κr}/r^{2} gives a logarithmic branch point at k = 2iκ. |

*** At least for low ionic densities, the leading term Yukawa function term in the presence of dispersion interactions is oscillatory with a wavelength that is much larger than the decay length,^{35} so in practice it appears like monotonic Yukawa functions. |

††† A full treatment of the theory that uses and generalizes the approach in Appendix D of ref. 29 will be published in a separate paper. |

‡‡‡ The details will be published elsewhere. |

§§§ The factorization is a consequence of a general result for matrices that turn singular at some parameter value because the determinant becomes zero, which in the present case occurs at k = ±iα. In the notation of Appendix D in ref. 29 the matrix H*(k) = {_{ij}*(k,ω_{1},ω_{2})} is the solution of the Ornstein–Zernike (OZ) equation (N + NH*N) ★ (N^{−1} − C*) = 1, where C*(k) = {_{ij}*(k,ω_{1},ω_{2})} is the matrix of the short-range part of the direct correlation function c_{ij}* ≡ c_{ij} + βu^{el}_{ij}, N is the diagonal matrix of the densities n^{b}_{i} and 1 is the unit matrix (note that there is a misprint in ref. 29 for the OZ equation). The inverse matrix (N^{−1} − C*)^{−1} = Adj(N^{−1} − C*)/D*, where Adj denotes the adjugate matrix and D* = Det(N^{−1} − C*) is the determinant, becomes singular because D* turns zero at k = ±iα. The adjugate matrix can in this situation be factorized at k = ±iα in the manner indicated. These and other matters regarding the direct correlation function approach to the general exact theory will be dealt with in more detail in a forthcoming publication by the author. |

This journal is © The Royal Society of Chemistry 2019 |