The intimate relationship between the dielectric response and the decay of intermolecular correlations and surface forces in electrolytes

Abstract

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.

---

1 Introduction

Electrolytes are prevalent in various systems of great importance in physics, chemistry, biosciences, surface and colloid sciences and many applied fields, including industrial applications. Such systems have been studied theoretically and experimentally for a very long time, but recently some unexpected experimental observations have been made, in particular the existence of long-range interactions in systems with high ion densities like ionic liquids and concentrated electrolyte solutions. These observations have renewed the interest in the basics of interactions in electrolytes.

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

where β = (k_BT)⁻¹, k_B is Boltzmann's constant, T is the absolute temperature, ε_r is the dielectric constant of the pure solvent modeled as a dielectric continuum (this is assumed in the primitive model of electrolyte solutions), ε₀ is the permittivity of vacuum, q_j is the charge of an ion of species j, and n^b_j is the number density in the bulk phase (superscript b stands for “bulk”). In the PB approximation one entirely neglects ion–ion correlations in the electrolyte adjacent to the object one considers. These ions are treated as point charges that do not correlate with each other. This applies also when the object one considers is an ion of the electrolyte itself (this is done when one deals with pair correlations in electrolytes), which means that this ion is treated in a different manner than the other ions of the same species. In a correct theory all ions should be treated on the same basis.

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 = 23 400 K (these systems have the same value of ε_rT). For a symmetric electrolyte with ionic diameter d this approximation gives the following expression for the screening parameter κ

which can be derived in a simple manner by making a very small extension of the linearized PB (LPB) approximation, whereby one requires that all ions in the systems should be treated on the same basis.‡ This expression can be written asso it is an equation for κ/κ_D as a function of κ_Dd. In the limit κ_Dd → 0 this equation has a solution with the property κ/κ_D → 1, but there is also another solution that we will denote as κ′/κ_D. These solutions are shown by the curves marked by κ and κ′ in Fig. 1 and we see that when κ_Dd is increased, the two solutions approach each other and at κ_Dd = 1.35 they merge (one can show that at this point), which for monovalent ions corresponds to a molar concentration of (4.10/d)² when d is measured in Å. For κ_Dd > 1.35, the solutions are complex-valued; the two solutions are then complex conjugates to each other: κ = κ_ℜ + iκ_ℑ and κ′ = κ_ℜ − iκ_ℑ, where i is the imaginary unit. The appearance of the complex solutions, marked by a filled symbol in the figure, means that there is a cross-over from a plain exponential decay to an oscillatory exponential one, that is, proportional to e^−κ_ℜℓ cos(κ_ℑℓ + γ) in planar geometry and e^−κ_ℜr cos(κ_ℑr + γ)/r in spherical geometry, where γ is a phase shift. This behavior is well-known and is called a Kirkwood cross-over,¹³ named after John G. Kirkwood who already in the 1930s showed^14,15 the existence of such a change to oscillatory behavior for electrolytes. Since then, several expressions for κ with the same behavior have been obtained, for example in the Mean Spherical Approximation (MSA), the Linearized Modified PB (LMPB) approximation by Outhwaite^16,17 and the generalized MSA.¹³ This behavior has also been verified many times in the past in numerical calculations, including simulations.

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 κ_Dd, 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¹⁸ 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¹⁸ 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 κ_Dd; 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 κ_Dd 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

where Λ = κ_Dℓ_B is the coupling parameter, ℓ_B = βq_e²/(4πε_rε₀) is the Bjerrum length, and q_e is the elementary (protonic) charge. The expression for ℓ_B used here is appropriate for the primitive model, otherwise ε_r = 1 both in ℓ_B and in eqn (1). The deviation of κ/κ_D from the value 1, as described by this law, is caused entirely by the long-range electrostatic correlations of the ions and is independent of other characteristics of the ions than their charges.§ Note that ln Λ is negative for small Λ, so formula (4) shows that the decay length is larger than the Debye length for low electrolyte concentrations in agreement with the numerical theoretical results mentioned above.

This law was recently verified experimentally by Smith et al.²³ 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 κ_Dd in these approximations. However, if one includes nonlinear terms,¹² 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 κ_Dd is very small for the monovalent electrolyte in Fig. 1. Due to the factor of z⁴ 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.,²⁴ 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,¹³ 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⁻³. The monotonic term has the largest decay length for n^b_tot < 0.1d⁻³, 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.⁹ 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,²⁷ where the model is a charge-inversion invariant system, and systems with unequally-sized ions.²⁸ 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 [small epsilon, Greek, tilde](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

where q_i* is a renormalized charge (called the dressed particle charge) for particles of species i and [scr E, script letter E]_r*(κ) is called the dielectric factor, which is obtained from [small epsilon, Greek, tilde](k) in a manner described later. Since the unknown variable κ appears on both sides in the expression, it is an equation for κ. The different solutions κ_ν, ν = 1, 2,…, to this equation comprise the set of screening parameters for the various decay modes. There exist in general both real and complex-valued solutions to the equation, which give rise to the monotonic and the oscillatory modes, respectively. An equivalent manner to write eqn (5) is, as we will see, [small epsilon, Greek, tilde](iκ) = 0, which explicitly involves the static dielectric function. The physical meaning of the appearance of the imaginary unit i in k = iκ is that the dielectric response is evaluated for an exponentially decaying field rather than a sinusoidal one (the latter corresponds to real k).

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 [small epsilon, Greek, tilde](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 [small epsilon, Greek, tilde](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.

2 Polarization response and nonlocal electrostatics in electrolytes

2.1 Polarization and external electrostatic potential

Consider an inhomogeneous electrolyte that is exposed to an external electrostatic potential Ψ^ext(r) from a source outside the system. The source may be a macroscopic body in contact with the electrolyte or a particle of any size and shape immersed in the electrolyte. The body or particle that is external to the system contains an arbitrary internal charge distribution that gives rise to Ψ^ext. Let us change the external potential somewhat by changing this charge distribution, so the system is instead exposed to the external potential Ψ^ext(r) + δΨ^ext(r) where δΨ^ext(r′) is assumed to be small everywhere in the electrolyte. Initially the charge distribution of the inhomogeneous fluid phase is ρ(r′), which can be obtained from the density distribution of constituent particles of the electrolyte. For the case of spherically symmetric ions with a point charge at the center, we have , where n_i(r′) is the number density distribution of ion species i (the corresponding formula for nonspherical ions and other particles will be given later). The change δΨ^ext in external potential makes the charge distribution become ρ(r′) + δρ(r′), where and δn_i is the variation in number density induced by the change in Ψ^ext. For reasons that soon will be apparent, this variation of ρ will be called a polarization charge density induced by δΨ^ext, so we have δρ^pol ≡ δρ.

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

where ϕ_Coul(r) = 1/(4πε₀r) is the (unit) Coulomb potential and the integral is taken over the whole space (this is the convention throughout the paper for integrals without limits). The variation in Ψ due to the influence of δΨ^ext is hence given by

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

whereis the potential from the polarization charge density δρ^pol induced by δΨ^ext.

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_ie^−βW_i(r′⁾, (8)

where n^b_i is the density of ions of species i in the bulk phase that is in equilibrium with the inhomogeneous phase. W_i is set to zero in the bulk phase. When n_i is changed by δn_i and W_i is changed by δW_i it follows from eqn (8) that δn_i(r′) = n^b_ie^−βW_i(r′⁾[−βδW_i(r′)], so we have

δn_i(r′) = −βn_i(r′)δW_i(r′), (9)

which will be of use later. In the Poisson–Boltzmann (PB) approximation W_i(r′) = q_iΨ(r′), while in reality there are also contributions to W_i that depend on ion–ion correlations, which are entirely neglected in the PB approximation.

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)

and the polarization charge density can be determined when we know δW_i. Thereby, a central question is: how can δW_i be determined from δΨ^ext or, in other words, what is the free energy of interaction δW_i(r′) of an i-ion at r′ when it interacts with the weak electrostatic potential δΨ^ext?

As we will show below

where ρ_i(r₁₂) is the charge density that surrounds an i-ion in the bulk phase, i.e., the charge density of the surrounding ion cloud which is given bywhere the sum is taken over all species, h_ij = g_ij − 1 is the total pair-correlation function and we have used the fact that , which follows from charge electroneutrality of the bulk phase. Eqn (11) says that the potential of mean force δW_i is given by the interactions of δΨ^ext with the i-ion itself, that is, with the charge q_i at the center of the ion [the first term on the right-hand side (rhs)], and with the ion cloud with density ρ_i [the integral], whereby the latter interaction affects the i-ion via the ion–ion correlations. Eqn (11) is an exact result, so it is precisely in this manner that ion–ion correlations enter into δW_i. For reasons that will be apparent in the following sections, we will, however, need to express the effects of such correlations in a more useful manner later.

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₂) acting on the various species j, we have

for spherical particles. For nonspherical particles a similar equation applies that we will consider later (see Appendix A). In our case of a weak external electrostatic potential δΨ^ext we have δv_j(r₂) = q_jδΨ^ext(r₂), which means thatwhere we have used the definition (12) of ρ_i. Since δn_i(r′) = −βn^b_iδW_i(r′) [eqn (10)] we can identify δW_i as the square bracket, which is equal to the rhs of eqn (11). Thereby the latter equation is demonstrated.

We can write the central charge q_i of the ion as a charge density q_iδ⁽³⁾(r), where δ⁽³⁾(r) is the three-dimensional Dirac delta function. By introducing ρ^tot_i(r) = q_iδ⁽³⁾(r) + ρ_i(r), which is the total charge density of the ion itself and the surrounding ion cloud, we can write eqn (14) as

and eqn (11) can be expressed asThis simple, exact relationship accordingly gives the free energy of interaction δW_i of the ion with the weak electrostatic potential δΨ^ext.

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

We can write this aswhereThe function χ^ρ(r₁₂) links the polarization charge density at one point, r₁, to the external electrostatic potential at another point, r₂. It is the polarization response function (a static, linear response function) that gives the contribution to δρ^pol at r₁ from the influence of the external electrostatic potential at r₂. The response function is a property of the unperturbed bulk fluid. By inserting the definitions of ρ^tot_i and ρ_i into eqn (18) it follows that
We can express this as χ^ρ(r₁₂) = −βq_e²H_QQ(r₁₂), where q_e is the elementary (protonic) charge and H_QQ is the charge–charge correlation function. The latter gives the correlation between fluctuations in charge densities at r₁ and r₂ irrespectively of which species the charges belong to and is for spherical ions equal to the square bracket divided by q_e². Its Fourier transform is equal to the charge–charge structure factor. These results for χ^ρ and H_QQ are well-known, but we have derived them here in a manner that is suitable for further development of the exact theory for electrolytes done later in this work.

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₃ and orientation ω₃, the internal charge density at point r₁ is given by σ_i(r₁|r₃,ω₃). 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₃ and ω₃. Note that σ_j(r₁|r₃,ω₃) for a given ω₃ is a function of only r₃₁ = r₁ − r₃, where the vector r₃₁ starts at the center of the particle. To simplify the notation we will henceforth write (r_ν,ω_ν) ≡ R_ν whereby we have σ_i(r₁|r₃,ω₃) ≡ σ_i(r₁|R₃), which is the charge density at r₁ for a particle with coordinates R₃.

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

and the latter is assumed to have a short range** (unless something else is stated explicitly) and to be strongly repulsive for small separations, but it is otherwise completely arbitrary.

The number density of particles with center of mass at r₃ and orientation ω₃ is equal to n_i(r₃,ω₃) ≡ n_i(R₃) and we have n_i(R₃) = n^b_ie^{−βW_i(R₃)}. The charge density of the system is

where dR₃ ≡ dr₃dω₃ and the integral is taken over the whole space and over all orientations. When a variation in external potential gives rise to a change δn_i in density, the resulting change in charge density iswhich is also applicable for the polarization charge density induced by the weak external potential δΨ^ext in a bulk phase with densities n^b_j for the various species j.

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

where

ρ^tot_i(r₂|R₁) = σ_i(r₂|R₁) + ρ_i(r₂|R₁) (22)

is the total charge density at r₂ around a particle with coordinates R₁ (i.e., it is located at r₁ and has orientation ω₁). This density consists of the internal charge density σ_i of the particle itself and the charge density ρ_i of the surrounding cloud of ions, solvent molecules and other particles present in the electrolyte. Electroneutral particles contribute to the distribution since they have orientational order due to the interactions with the i-particle. For simplicity, ρ_i will be denoted as the charge density of the ion/solvent cloud that surrounds the particle. In the same manner as before it follows that

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

which can be written as eqn (17) when we make the identificationwhere r₁₂ = |r₁₂| and r₁₂ = r₂ − r₁. As before, the response function can be expressed in terms of the charge–charge correlation function as χ^ρ(r₁₂) = −βq_e²H_QQ(r₁₂), where H_QQ for the case of nonspherical particles is defined in Appendix B [eqn (134)].

We can write eqn (24) as

where we have written the integral over ω₃ as an average since . Note that χ^ρ(r₁₂) depends only on the distance r₁₂ because the average over orientations is taken.

2.2 Polarization and total electrostatic potential

To continue we note that Ψ^ext is the electrostatic potential due to the external source in the absence of the electrolyte, while Ψ is the total potential in the presence of the electrolyte. This means that Ψ is the mean potential that actually exists in the electrolyte. Since this is the system that we are interested in, it is a relevant task to express δρ^pol in terms of δΨ rather than in terms of δΨ^ext as in eqn (17). Thereby we will obtain a relationship that can be used as an alternative to the latter equation.

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

which is similar to eqn (17). The function χ*(r) will also be denoted as a “polarization response function,” although a true response function normally gives the response of a system due to an influence that we can completely control by external means. The total potential δΨ contains δΨ^pol from the polarization charge density that we cannot control directly, so χ*(r) is a slightly different kind of function than χ^ρ(r).

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

which we will derive below. The solution is readily obtained in Fourier space aswhere the transform of the unit Coulomb potential is [small phi, Greek, tilde]_Coul(k) = 1/(ε₀k²). We define the Fourier transform f̃(k) of a function f(r) as For a function f(r) that only depends on r we have where k = |k|.

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

where we have changed the integration variable in the first integral on the rhs from r₂ to r₄. We now insert eqn (17) and obtain

Since the left-hand side (lhs) is equal to and δΨ^ext(r₄) 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₁) = q_iδΨ(r₁) and hence from eqn (10) we obtain δn_i(r₁) = −βn^b_iq_iδΨ(r₁). Therefore

where (PB) means that the equation is valid only in the PB approximation. This equation can be written as
and by comparing with the expression (26) for χ*, we can identifyIn eqn (29) we see that the polarization charge density δρ^pol(r₁) is proportional to the total mean electrostatic potential δΨ(r₁) at the same pointr₁. This fact is likewise expressed by the appearance of the Dirac delta function δ⁽³⁾(r₁₂) in χ*(r₁₂) for the PB case. Hence, electrostatics is local in the PB approximation: δρ^pol(r₁) is not influenced by the values of δΨ(r₂) at other points r₂ ≠ r₁.

In reality, we have nonlocal electrostatics in the sense that the polarization at one point r₁ is influenced by the electrostatic potential at other points in the surroundings. This can be seen in eqn (26). The polarization response function χ*(r₁₂) is nonzero for r₁₂ ≠ 0 so δρ^pol(r₁) is influenced by δΨ(r₂) for all points r₂ in the neighborhood of r₁. 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₂ interacts with the electrostatic potential and since this ion correlates with other ions it will affect the probability for ions to be at r₁. The density of ions at r₁ accordingly depends on the potential elsewhere. In other words, the potential of mean force δw_i(r₁) and hence the charge density δρ^pol(r₁) depend on the values of δΨ(r₂) for all points r₂ in the vicinity of r₁. 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 [small epsilon, Greek, tilde](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) = ε₀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 [small epsilon, Greek, tilde](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

where k̂ = k/k is a unit vector and where the field components along the wave vector k (the longitudinal components) are projected out. In electrostatics Ẽ and Ẽ^ext are parallel to k, which can be seen from the fact that in Fourier space we have Ẽ(k) = −ik[capital Psi, Greek, tilde](k) (the factor ik corresponds to ∇ in ordinary space) and likewise for Ẽ^ext. Hence we have k̂·Ẽ(k) = −ik̂·k[capital Psi, Greek, tilde](k) = −ik[capital Psi, Greek, tilde](k) and k̂·Ẽ^ext(k) = −ik[capital Psi, Greek, tilde]^ext(k). Therefore only the longitudinal components of the fields matter. In electrodynamics the transversal components, which are perpendicular to k̂, also matter, but in the present work we limit ourselves to the static equilibrium case. The frequency dependence of the dielectric function, that is commonly considered in the study of the polarization of fluids, is therefore not included in the treatment.

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

In Fourier space eqn (7), (17) and (26) areand since δ[capital Psi, Greek, tilde] = δ[capital Psi, Greek, tilde]^ext + δ[capital Psi, Greek, tilde]^pol we obtain
By comparing the last two equations with eqn (32) we can identifyNote that the last equality is equivalent to eqn (28). These expressions for [small epsilon, Greek, tilde](k) in terms of [small chi, Greek, tilde]^ρ(k) and [small chi, Greek, tilde]*(k) are equivalent to standard expressions for [small epsilon, Greek, tilde](k) given in ref. 36–38.

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 [small epsilon, Greek, tilde](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

and 1/κ_D is the Debye length. We see that [small chi, Greek, tilde]*(k) is independent of k in the PB case and we obtainIn this case, the entire k dependence of [small epsilon, Greek, tilde](k) originates from the Coulomb potential [small phi, Greek, tilde]_Coul(k).

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 P̃ = ε₀[small chi, Greek, tilde]_eẼ for weak fields.†† In our notation we have P = −ε₀E^pol, where E^pol(r) = −∇Ψ^pol(r). For weak fields we thus have −ε₀δẼ^pol = ε₀[small chi, Greek, tilde]_eδẼ or in terms of potentials δ[capital Psi, Greek, tilde]^pol = −[small chi, Greek, tilde]_eδ[capital Psi, Greek, tilde]. By comparing with eqn (33) we see that

so the electric susceptibility is related to [small epsilon, Greek, tilde] as [small epsilon, Greek, tilde](k) = 1 + [small chi, Greek, tilde]_e(k) as usual. In traditional treatments of polar media, P expresses dipolar polarization, while in the general case, including electrolytes, E^pol and hence P originate from all kinds of polarization, for example from changes in dipolar, quadrupolar, octupolar orientations and ion distributions.³¹ For a pure polar fluid (no ions), the dielectric constant (relative dielectric permittivity) is defined microscopically as and in this limit solely the dipolar polarization contributes. In the case of electrolytes this limit does not exist since [small epsilon, Greek, tilde](k) diverges at k = 0 and, as we will see, other k values matter a lot.

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

3 Dressed particles and screened electrostatic interactions

3.1 Definition of dressed particles

The potential of mean force δW_i due to a weak external electrostatic potential δΨ^ext was expressed in eqn (23) as the interaction energy between δΨ^ext and the total charge density ρ^tot_i associated with an i-particle, where ρ^tot_i is the sum of the internal charge densities σ_i of the particle and ρ_i of the surrounding ion/solvent cloud, i.e., the surrounding cloud of particles of any kind present in the system. We will now determine the corresponding relationship between δW_i and the total electrostatic potential δΨ. Thereby we are led to the definition of the concept of dressed particles and the corresponding charge density ρ_i*, which plays a fundamental role in the understanding of interactions in electrolytes.

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

and it remains to express the last term in terms of δΨ. By inserting δΨ^pol from eqn (7) into the last term, it can be written
whereis the electrostatic potential from the charge density ρ^tot_i, that is, the mean electrostatic potential due to the i-particle located at r₁ and with orientation ω₁. We can now express δρ^pol in terms of δΨ by using eqn (26) and by inserting the result into eqn (37) we obtainIf we definewe can write eqn (39) asBy comparing this expression with eqn (23) we see that when δW_i is expressed in terms of δΨ instead of δΨ^ext, the charge density ρ_i* takes the role corresponding to that of ρ^tot_i in eqn (23). Since δn_i(R₁) = − βn^b_iδW_i(R₁) [cf.eqn (10)] we obtain from eqn (41)which can be compared with eqn (21).

We will now interpret the physical meaning of ρ_i*. Let us immerse a particle of species i with orientation ω₁ into a bulk electrolyte at the point r₁. Initially the charge density is zero and after the immersion the density at r₂ is ρ^tot_i(r₂|R₁). 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

since ψ_i acts like the total electrostatic potential due to an external source, that is, the potential ψ_i from the immersed i-particle here acts like δΨ. Since ψ_i is not weak everywhere, this expression does not give the total polarization due to ψ_i, but instead the rhs gives the linear part of the electrostatic polarization response due to ψ_i and we therefore define in the general case this part asThe rest of the polarization of the bulk electrolyte caused by the interactions between the immersed i-particle and the other particles is contained in ρ^tot_i − ρ^lin_i. We now define the dressed particle charge–densityρ_i* of the particle as

ρ_i*(r₂|R₁) = ρ^tot_i(r₂|R₁) − ρ^lin_i(r₂|R₁), (44)

which is the same as eqn (40). Since ρ^tot_i = σ_i + ρ_i we can write

ρ_i*(r₂|R₁) = σ_i(r₂|R₁) + ρ^dress_i(r₂|R₁), (45)

where ρ^dress_i ≡ ρ_i − ρ^lin_i is the charge density of the “dress” of the i-particle. The dress expresses the effects of many-body correlations between the particles in the electrolyte in addition to those included in the linear response. Like the charge density ρ_i of the ion/solvent cloud, which is given by the pair distribution function g_ij, the dress is a distribution of particles of all kinds surrounding the i-particle. The pair distribution function g_ij* that gives the charge distribution ρ^dress_i will be investigated in Section 4.

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_jviaeqn (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

in the same manner as the derivation of eqn (25) from eqn (21). Note the similarity of the expression for χ* and the corresponding expression for χ^ρ. The only difference is that ρ_i* here takes the role of ρ^tot_i in eqn (25).

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

f(r₁₂,ω₁) ≡ f(r₂|r₁,ω₁) = f(r₂|R₁) (47)

for the functions σ_i and ρ_i*, where we have used the fact mentioned earlier that such functions only depend on the r₁₂ vector drawn from the center of the particle. We obtainand in Fourier space we therefore havewhere −k appears in [small sigma, Greek, tilde]_i since the integral over r₃ in eqn (48) becomes a convolution when r₃₁ is replaced by −r₁₃ = r₃₁. Since we take the average over ω₃, the direction of the unit vector k̂ = k/k is arbitrary. Eqn (49) implies that the dielectric function [small epsilon, Greek, tilde](k) = 1 − [small chi, Greek, tilde]*(k)[small phi, Greek, tilde]_Coul(k) can be expressed asThe reciprocal function 1/[small epsilon, Greek, tilde](k) = 1 + [small phi, Greek, tilde]_Coul(k)[small chi, Greek, tilde]^ρ(k) can likewise be expressed in terms of [small rho, Greek, tilde]^tot_i and thereby in terms of H_QQvia the Fourier transform of eqn (24). The latter is the usual path to [small epsilon, Greek, tilde](k), while we will see that eqn (50) is much more useful.

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

and we have ρ_i*(r₁₂) ≡ ρ^tot_i(r₁₂) − ρ^lin_i(r₁₂). In such cases eqn (45) becomes

ρ_i*(r₁₂) = σ_i(r₁₂) + ρ^dress_i(r₁₂) = q_iδ⁽³⁾(r₁₂) + ρ^dress_i(r₁₂), (52)

where ρ^dress_i(r₁₂) ≡ ρ_i(r₁₂) − ρ^lin_i(r₁₂). Eqn (41) can be written aswhere the first term on the rhs is the sole contribution to W_i included in the PB approximation and the last term describes the effects of ion–ion correlations. In the PB approximation we thus neglect that the ions have dresses when evaluating δW_i. For cases where all particles in the electrolyte are spherical eqn (46) reduces towhich should be compared with eqn (18). Note that χ*(r₁₂) for the the PB case, eqn (30), is obtained from eqn (54) when we set ρ_i*(r₁₂) = q_iδ⁽³⁾(r₁₂), that is, the charge density of a bare ion without its dress. This is in agreement with eqn (53) in the absence of its last term. In the PB approximation, all particles in the electrolyte are treated as being bare, except for the particle that causes the potential δΨ. That particle does have a dress in this approximation, as we will now see for the special case when the particle is an ion of species j. This is very illustrative because it gives some insights into the concept of a dressed particle in this simple approximation.

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

on the inside and outside, respectively, of the j-particle. Let us expand the exponential function in a Taylor series. For positions r outside the particle we get
From eqn (43) with χ*(r) given by eqn (30) we have
so from eqn (44) it follows thatsince the linear term in ρ^tot_j is canceled by −ρ^lin_i outside the particle. We see that ρ_j* is equal to the nonlinear terms of ρ^tot_j on the outside. Since ψ_j(r|R₁) and ρ^tot_j(r|R₁) in the PB approximation decay as e^−κ_Dr/r when r → ∞, we see that ρ_j*(r|R₁) decays as (e^−κ_Dr/r)². Thus ρ_j* has a considerably shorter range than ρ^tot_j; the former has half the decay length of the latter.

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 The screened Coulomb potential in the general case

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 [small epsilon, Greek, tilde](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.

The mean electrostatic potential ψ_i(r|R₁) from a particle with coordinates R₁ 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

−ε₀∇²ψ_i(r|R₁) = ρ^tot_i(r|R₁) (57)

with the boundary condition ψ_i(r|R₁) → 0 when r → ∞. By subtracting ρ^lin_i(r|R₁) from both sides of Poisson's equation we obtain

−ε₀∇²ψ_i(r|R₁) − ρ^lin_i(r|R₁) = ρ^tot_i(r|R₁) − ρ^lin_i(r|R₁) ≡ ρ_i*(r|R₁)

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 ofand the solution of eqn (58) is thereby given byas 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.

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

−ε₀[∇²ϕ_Coul*(r) − κ_D²ϕ_Coul*(r)] = δ⁽³⁾(r) (PB), (61)

which has the solutiona monotonic Yukawa function. Thereforein this approximation.

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

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 δ[capital Psi, Greek, tilde]^ext(k) = δq[small phi, Greek, tilde]_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.

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₁) 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

[small psi, Greek, tilde]_i(k,ω₁) = [small rho, Greek, tilde]^tot_i(k,ω₁)[small phi, Greek, tilde]_Coul(k) = [small rho, Greek, tilde]_i*(k,ω₁)[small phi, Greek, tilde]_Coul*(k). (65)

By using [small phi, Greek, tilde]_Coul*(k) = [small phi, Greek, tilde]_Coul(k)/[small epsilon, Greek, tilde](k) = 1/|ε₀k²[small epsilon, Greek, tilde](k)] we see from the last equality thatwhen the factor 1/|ε₀k²] has been removed. This simple relationship can also be derived from eqn (40).

In the PB approximation where [small epsilon, Greek, tilde](k) is given by eqn (36) we have

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, [small phi, Greek, tilde]_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 [small phi, Greek, tilde]_Coul*(k), we have

[ε₀k² − [small chi, Greek, tilde]*(k)]_k=±iκ = 0, [small epsilon, Greek, tilde](±iκ) = 0, (68)

where the latter holds for κ ≠ 0. Since [small phi, Greek, tilde]_Coul*(k) is the Fourier transform of a real function of r, [small phi, Greek, tilde]_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.

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 [small phi, Greek, tilde]_Coul*(k) has a pole k = iκ when

which implies thatwhere 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.

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₁) ≈ q_iδΨ(r₁) 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*

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.

3.2.2 The general equation for the screening parameter κ vs. the dielectric function. Let us now focus on the dielectric function [small epsilon, Greek, tilde](k) and use the condition [small epsilon, Greek, tilde](iκ) = 0 in eqn (68) for the pole of [small phi, Greek, tilde]_Coul*(k). Since [small epsilon, Greek, tilde](k) = 1 − [small chi, Greek, tilde]*(k)[small phi, Greek, tilde]_Coul(k) we can writewhere we have split [small epsilon, Greek, tilde](k) into two partswhich is singular when k → 0 (provided that [small chi, Greek, tilde]*(0) ≠ 0), andwhich 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.

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

whereis 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_iq_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).

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

where we have used
which is independent of the orientation ω₃, 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.

By inserting this result into eqn (72) we obtain

Due to the presence of [small epsilon, Greek, tilde]_sing(k), the dielectric function [small epsilon, Greek, tilde](k) for electrolytes diverges to infinity when k → 0. This divergence is commonly called perfect screening.

Using eqn (77) we see that the condition for a pole, [small epsilon, Greek, tilde](iκ) = [small epsilon, Greek, tilde]_reg(iκ) + [small epsilon, Greek, tilde]_sing(iκ) = 0, is

which can be writtenwhere

[scr E, script letter E]_r*(κ) = [small epsilon, Greek, tilde]_reg(iκ) (79)

is the dielectric factor. Eqn (78) is the general equation for the screening parameter κ, eqn (5). It is obviously equivalent to [small epsilon, Greek, tilde](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 [scr E, script letter E]_r*(κ) in the denominator differ. From eqn (73) it follows thatAs we will see, the fact that [scr E, script letter E]_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 [scr E, script letter E]_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.

For a binary electrolyte, the deviation in κ from κ_D can be obtained from

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.

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 [small epsilon, Greek, tilde]_reg′(k) = d[small epsilon, Greek, tilde]_reg(k)/dk and we have used l'Hospital's rule to obtain the last equality. If we defineit follows that A* = 1/[ε₀[scr E, script letter E]^eff_r(κ)] and hence we haveBy 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/[scr E, script letter E]^eff_r(κ). In addition there is a change in the value of the screening parameter from κ_D to κ. Note that [scr E, script letter E]^eff_r(κ) differs from [scr E, script letter E]_r*(κ) by the last term in eqn (82).

Alternative expressions for [scr E, script letter E]^eff_r are

where the first expression yields eqn (82) after the differentiation and where we have used [small epsilon, Greek, tilde](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 [scr E, script letter E]^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.)

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 [small epsilon, Greek, tilde]_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 [scr E, script letter E]_r*(κ) [apart from their indirect influence on the distribution functions and other entities, of course].

For a pure polar solvent without ions, [small epsilon, Greek, tilde]_sing(k) vanishes so [small epsilon, Greek, tilde](k) = [small epsilon, Greek, tilde]_reg(k). As mentioned earlier the dielectric constant of a pure polar medium is defined microscopically as , so we have ε_r = [small epsilon, Greek, tilde](0) = [small epsilon, Greek, tilde]_reg(0). In a very dilute electrolyte solution κ ≈ 0 and hence [scr E, script letter E]_r*(κ) ≈ [scr E, script letter E]_r*(0) ≈ ε_r. Eqn (78) then becomes

where we have used the fact that q_i* ≈ q_i for a dilute electrolyte. In the last equality we have also identified κ_D² for the electrolyte solution when the solvent is a dielectric continuum with dielectric constant ε_r, eqn (1).

In dilute solutions we have [scr E, script letter E]^eff_r(κ) ≈ [scr E, script letter E]^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 → ∞

with κ ≈ κ_D, that is, the same as in the PB approximation when the pure solvent has dielectric constant ε_r.

For electrolyte solutions in general, [scr E, script letter E]^eff_r(κ) and [scr E, script letter E]_r*(κ) contain, as we have seen, contributions from the ions. Since 1/[scr E, script letter E]^eff_r(κ) determines the magnitude of the screened Coulomb potential for large r, [scr E, script letter E]^eff_r(κ) can be designated as the effective relative dielectric permittivity of the entire electrolyte solution. Likewise, the dielectric factor [scr E, script letter E]_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 [scr E, script letter E]^eff_r(κ) and [scr E, script letter E]_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 [scr E, script letter E]_r*(κ), that the former contains the factor q_iq_i* instead of q_i². 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 [scr E, script letter E]_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 [scr E, script letter E]_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¶¶

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 [scr E, script letter E]^eff_r.

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 Z̲ is given by Z + Z̲ = 2ℜ(Z), where ℜ(Z) stands for the real part of Z, we then have

where we have written [scr E, script letter E]^eff_r(κ_ℜ + iκ_ℑ) = |[scr E, script letter E]^eff_r|e^{−iϑ_{[scr E, script letter E]}} with a real ϑ_{[scr E, script letter E]}. The two poles hence give rise to an exponentially decaying, oscillatory term with decay length 1/κ_ℜ, wavelength 2π/κ_ℑ and phase shift −ϑ_{[scr E, script letter E]}. 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.”

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¹³ 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

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).

Recall that eqn (78) is equivalent to [small epsilon, Greek, tilde](iκ) = 0, so before the cross-over κ and κ′ are two consecutive zeros of [small epsilon, Greek, tilde](iξ) as a function of the real variable ξ. As we have seen, [scr E, script letter E]^eff_r(κ) > 0 for low ion densities, so from the rhs of eqn (84) we see that [small epsilon, Greek, tilde]′(iξ) > 0 for ξ = κ. The next zero of [small epsilon, Greek, tilde](iξ) must have an opposite derivative, so [small epsilon, Greek, tilde]′(iξ) < 0 for ξ = κ′, which implies that [scr E, script letter E]^eff_r(κ′) < 0. Thus the two terms in eqn (88) have opposite signs. At the the cross-over point, where the two zeros of [small epsilon, Greek, tilde](iξ) merge, we must have [scr E, script letter E]^eff_r(κ) = [scr E, script letter E]^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¹³ 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.

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

where r̂ = r/r. By setting r = r₁₂ and r′ = r₁₃ (note that r − r′ = r₂ − r₃) and using the notation introduced in eqn (47), we can conclude from eqn (89) in the limit r₁₂ → ∞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₁₂, whereby the integral contains the direction dependence as expressed viar̂₁₂ in the exponent.

For a spherically symmetric particle, the integral becomes

and we obtain from the terms in eqn (86) in the limit r₁₂ → ∞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 κ_ν.

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

whereis 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

Q̄^eff_i(κ) = 〈Q^eff_i(r̂,ω,κ)〉_ω (95)

is the orientation average of the effective multipolar charge. Q̄^eff_i can be described as a kind of monopolar effective charge of the particles. The rest of Q^eff_i(r̂₁₂,ω₁,κ_ν) has an orientational angle dependence with a combination of dipolar, quadrupolar and higher multipolar characteristics.³⁴ Note that electroneutral particles, like solvent molecules, acquire nonzero Q̄^eff_i due to the presence of ions in their neighborhood, so they have nonzero effective charges in general.

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.³⁰

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

For each pole k = iκ_ν of [small phi, Greek, tilde]_Coul*(k), that is, zero of [small epsilon, Greek, tilde](k), we obtain a contribution to ψ_i(r₁₂,ω₁) that decays as e^−κ_νr/r with the same prefactor 1/[4πε₀[scr E, script letter E]^eff_r(κ_ν)] as in eqn (86) times the factorwhere we in ordinary space apply this result with k̂ = r̂₁₂ 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.

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 [small epsilon, Greek, tilde](iκ) = 0 and, equivalently, to the general equation for κ, eqn (78). By inserting k = iκ into the expression (50) for [small epsilon, Greek, tilde](k) we obtain

Inserting eqn (96) and using [small epsilon, Greek, tilde](iκ) = 0 we find thatwhere we have definedEqn (97) is an equation for κ that is equivalent to eqn (78).

For a spherical ion with charge q_i at the center we have [small sigma, Greek, tilde]_i(k) = q_i so we have Q_i(−k̂,ω,κ) = q_i and Q^eff_i(k̂,ω,κ) = q^eff_i(κ). When we deal with a system consisting solely of such ions, eqn (97) becomes [cf.eqn (69)]

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 [scr E, script letter E]_r*(κ) in the denominator. Incidentally, we may note that for spherical ions we have , but in general q^eff_i(κ) ≠ q_i*/[scr E, script letter E]_r*(κ). For a binary electrolyte of spherical ions, the deviation in κ from κ_D can be obtained fromwhich should be compared to eqn (81).

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*/[scr E, script letter E]_r*(κ) and from eqn (100) it follows that

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, iswhere κ = κ_D and prefactor for ψ_i(r) is dictated by local electroneutrality, . We can identify the effective charge as q^eff_i = q_ie^κd/(1 + κd) whereby we have ψ_i(r) = q^eff_ie^−κr/[4πε_rε₀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).

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.³⁴ 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² and the term with l = 3 goes to a purely quadrupolar term that decays as 1/r³. Thereby the usual multipole expansion is obtained, which applies to the electrostatic potential from a fixed charge distribution immersed in a pure polar liquid.

4 Screened interactions in electrolytes in the general case

4.1 The main case: all decay modes of the correlations are determined by the dielectric function

As we saw earlier, the electrostatic potential from a particle can be regarded as an external potential for the system even when the particle belongs to the same species as one of the constituent ones. The total mean potential ψ_j(r|R₂) from a j-particle with coordinates R₂ can then be treated in the same manner as Ψ(r) and when ψ_j is weak, the pair potential of mean force w_ij(R₁,R₂) for a particle of species i with coordinates R₁ satisfies according to eqn (41)
provided that the screened electrostatic interactions between the two particles dominate in w_ij. Here w_ij has taken the role of δW_i and ψ_i that of δΨ in eqn (41). With this in mind, we define in the general case the screened electrostatic part of the potential of mean force aswhere we have made use of eqn (60) for the potential ψ_j due to the j-particle. Note that the i-particle and the j-particle are treated in a symmetric manner and that we can also write . The integral on the rhs of eqn (103) has the same form as a Coulomb interaction energy between the charge densities ρ_i* and ρ_j*, but with the screened Coulomb potential ϕ_Coul*(r) instead of the usual unscreened one. As we will see w^el_ij has a central role for the interparticle interactions in electrolytes. For spherically symmetric particles we haveSince eqn (103) and (104) constitute definitions of w^el_ij, they can be used irrespectively of the magnitude of the potentials.

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^−βw_ij. We have the expansions

so for large separations, where h_ij and w_ij are small, these two functions decay in the same manner, that is, h_ij ∼ −βw_ij.

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

Recall that the total charge density ρ^tot_i associated with a spherical i-ion is given by
We will now show that the dressed ion charge density is given bywhich means that the functions h_ij* and g_ij* ≡ 1 + h_ij* have the same roles vis-à-vis ρ_i* as the usual pair functions h_ij and g_ij have vis-à-vis ρ^tot_i. In other words, the function g_ij* is the distribution function of the dress. Eqn (107) can easily be realized from the fact that eqn (106) implies that
where we have used eqn (54) to obtain the last equality. The rhs of this equation is equal to ρ_i*(r₁₂) according to its definition (40), so eqn (107) follows.

In the general case we define in an analogous manner

and hencewhere we have inserted eqn (60). The charge densities ρ^tot_i and ρ_i* are in this case given bywhere one can prove the expression for ρ_i* by using eqn (40) and (46) and analogous arguments as in the derivation of eqn (107). Also in this case, h_ij* gives ρ_i* in the same way as h_ij gives ρ^tot_i.

By using the definition (103) of w^el_ij we see that eqn (109) can be written as

h_ij(R₁,R₂) = h_ij*(R₁,R₂) − βw^el_ij(R₁,R₂). (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 [small epsilon, Greek, tilde](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_ijhas 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₁,R₂) in the bulk phase depends on the separation vector r₁₂ = r₂ − r₁ we can write h_ij(R₁,R₂) = h_ij(r₁₂,ω₁,ω₂), so in Fourier space we obtain from eqn (103) and (112)

where we have used eqn (64). For spherical particles this equation reduces toThese two results, which we have obtained here by physical reasoning, are key equations in DMT and DIT, respectively, that have been derived earlier^21,25,33 from the Ornstein–Zernike (OZ) equation. The entire theory of the present work, including the concept of dressed particles, can be formulated in terms of direct correlation functions, but in this work this has been avoided because the theory is then more accessible for a wider readership.

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 h̃_ij to be given by the zeros of [small epsilon, Greek, tilde](k), whereby they coincide with the poles of [small phi, Greek, tilde]_Coul*(k). This follows from the fact, shown in Appendix C, that h̃_ij and h̃_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 h̃_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.²⁵ Since [small rho, Greek, tilde]_i* and [small rho, Greek, tilde]_j* are given by linear combinations of h̃_ij*, their poles are also poles of h̃_ij* and cannot coincide with any pole of h̃_ij. Therefore, all poles of the rhs of eqn (113) [for spherical ions eqn (114)] originate from the zeros of [small epsilon, Greek, tilde](k) in the denominator and the assertions in (i) and (ii) above follow. The dielectric function [small epsilon, Greek, tilde](k) is also a linear combination of h̃_ij* since [small rho, Greek, tilde]_i* in eqn (50) is such a linear combination. The poles and zeros of [small epsilon, Greek, tilde](k) occur, of course, for different k values.

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

in the limit r₁₂ → ∞. The “other terms” in this equation are additional Yukawa function terms with different κ values (shorter decay lengths) that are solutions to eqn (78) and functions with different functional dependences of r₁₂ that normally decay faster than the leading term. The latter kind of functions always arises for reasons explained at the end of the current section.

For nonspherical particles, an expression for h_ij(r₁₂,ω₁,ω₂) analogous to that in eqn (115) applies when r₁₂ → ∞

where we have only shown the first term for κ_ν = κ, κ′ and κ″ [cf.eqn (93)]. Note that the vectors r̂₁₂ in Q^eff_i and −r̂₁₂ in Q^eff_j point towards the center of the other particle along the connecting line. As noted in Section 3.2.4, for these kinds of systems, each Yukawa function term is always accompanied by terms that decay as e^−κr/r^l with l = 2, 3,… They are included in “other terms” together with additional terms that will be discussed at the end of this section.

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

where γ_i = γ_i(r̂₁₂,ω₁,κ), γ_j = γ_j(−r̂₁₂,ω₂,κ) and where we likewise have suppressed the arguments of Q^eff_i and Q^eff_j. Both γ_l and |Q^eff_l| depend on orientation. This kind of oscillatory term is, of course, formed from any pair of complex conjugate screening parameters κ_ν that are solutions to eqn (78), say, κ_ν and .

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 ∼ −β⁻¹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 h̃_ij, which are in general determined by the zeros of [small epsilon, Greek, tilde](k). This can be realized as follows. H̃_NN(k), given in Appendix B, eqn (138), is a linear combination of h̃_ij,

so it cannot have any other poles than those of h̃_ij. By inserting eqn (66) into H̃_QQ(k) given in eqn (136) and H̃_NQ(k) given in eqn (140) we obtainand we see that the zeros of [small epsilon, Greek, tilde](k) are poles of these functions. Thus H_NN(r), H_QQ(r), H_NQ(r), h_ij, w_ij and ϕ_Coul*(r) have the same screening parameters, apart from in the exceptional cases mentioned earlier. They will be treated in the next section.

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]², [e^−κr/r]³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κ)⁻¹, (3κ)⁻¹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.²¹ 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 [small epsilon, Greek, tilde](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⁻⁶ dispersion interactions, the functions h_ij(r), w_ij(r) and ϕ_Coul*(r) ultimately decay like a power law when r → ∞.³⁵ 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.²⁴ 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.

4.2 Exceptions; charge-inversion invariant systems

Let us now turn to the exceptional cases, that is, cases where h̃_ij and h̃_ij* can have poles for the same k values and where Yukawa function contributions to h_ij therefore originate both from the last term in eqn (113) and from h̃_ij*. As shown in Appendix D of ref. 29 this can occur for most systems at exceptional points in the system's parameter space (for instance a critical point).††† This will not be dealt with any further because the systems behave in the way we have just described for all other parameter values.

The other exceptional case is a category of systems where h̃_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₁) for each particle changes to −σ_i(r|R₁) (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₁), 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 h̃_ij arises from the zeros of [small epsilon, Greek, tilde](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 [small rho, Greek, tilde]_i*(k,ω₁) for each particle changes to −[small rho, Greek, tilde]_i*(k,ω₁), so we must have

and hence H̃_NQ(k)≡0.

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

so the poles of H̃_NN(k) must be poles of h̃_ij*. Eqn (119) for H̃_QQ(k) remains, however, valid and the poles of H̃_QQ(k) are hence given by the zeros of [small epsilon, Greek, tilde](k) also in the present case. Thus there are two different sets of poles that are relevant, the poles of H̃_QQ(k) [zeros of [small epsilon, Greek, tilde](k)] and those of H̃_NN(k). The poles of the pair correlation function h̃_ij belong to both sets, but as we will see they enter in h̃_ij in different manners.

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₁) = −σ_A⁻(r|R₁). 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₁) = −σ_B⁻(r|R₁). They are also identical to each other apart from the sign of σ_i. The charge-inversion invariance implies that we have the symmetry

(we also have h_A⁺A⁻ = h_A⁻A⁺ and h_B⁺B⁻ = h_B⁻B⁺, as always). The first two lines are special cases of the last two (with A = B), so in the following we will only deal with the latter. The same symmetry relationships are valid for h_ij*. In the presence of more than four species of ions the corresponding relationships are valid for A, B, C, D, etc.

We now define h_S,AB(r,ω₁,ω₂), h_D,AB(r,ω₁,ω₂) and the corresponding h* functions from

(S stands for sum and D for difference) so we haveand likewise for h_ij*. Note that + sign on the rhs of eqn (123) applies to ions with equal sign of their charges (like A⁺B⁺) and − sign applies to ions with different signs (like A⁺B⁻). One can also define S and D functions for the electroneutral particles (like solvent molecules) in charge-invariant systems, but we will not enter into any details here.‡‡‡

Due to the charge-inversion invariance we have [small rho, Greek, tilde]_A⁺* = −[small rho, Greek, tilde]_A−* ≡ [small rho, Greek, tilde]_A*, which defines [small rho, Greek, tilde]_A*, and likewise for [small rho, Greek, tilde]_B*. Therefore we have from eqn (113)

h̃_S,AB(k,ω₁,ω₂) = h̃_S,AB*(k,ω₁,ω₂) (124)

Analogously to the arguments about the poles of h̃_ij* and h̃_ij, it follows that h̃_D,AB* cannot have poles common to h̃_D,AB so the poles of h̃_D,AB are due to the zeros of [small epsilon, Greek, tilde](k) [possibly apart from exceptional points in the system's parameter space]. The functions h_S,AB and h_D,AB have different decay parameters; the poles of h̃_S,AB belong to one of the sets mentioned earlier [the poles of H̃_NN(k)] and those of h̃_D,AB belong to the other set [the poles of H̃_QQ(k)]. Thus only h_D,AB has the same screening parameters as ϕ_Coul*(r) and h_D,AB/2 decays analogously to the rhs of eqn (115)–(117). The contribution from h_D,AB in eqn (123) appears with a plus sign for ions of the same sign of charge and with a minus sign for ions of opposite sign. The function h_S,AB (with different decay parameters than h_D,AB) appears with a plus sign in eqn (123) for all ions irrespectively of their signs of charge.

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,”³⁹ 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 h̃_ij are given by the zeros of [small epsilon, Greek, tilde](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 h̃_ij* and those that are due to zeros of [small epsilon, Greek, tilde](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 [small epsilon, Greek, tilde](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 h̃_ij* but not zeros of [small epsilon, Greek, tilde](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 h̃_ij and h̃_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 h̃_ij [zeros of [small epsilon, Greek, tilde](k)] where each pole is close to a pole of h̃_ij* and where the trajectories of the poles of h̃_ij and h̃_ij* cross each other at the point where the system become invariant. At the crossing point the pole of h̃_ij ceases to be a zero of [small epsilon, Greek, tilde](k) and becomes instead a pole of both h̃_ij* and h̃_ij. The second set of poles of h̃_ij do not have such crossings and these poles remain zeros of [small epsilon, Greek, tilde](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 [small epsilon, Greek, tilde](k) for the nearly invariant system that lie close to the poles of h̃_ij* for the same system. The latter are not poles of h̃_ij while the zeros of [small epsilon, Greek, tilde](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 h̃_ij* and ceases to be a zero of [small epsilon, Greek, tilde](k) but remains as a pole of h̃_ij.

4.3 The decay of H_QQ(r), H_NN(r) and H_NQ(r)

Let us now return to the general case of systems without charge-inversion symmetry, which is the normal, realistic case since anions and cations virtually always differ by much more than the sign of their charge and since the positive and negative parts of the solvent molecules normally differ a lot apart from the sign of the charge.

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 [small epsilon, Greek, tilde](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 the coefficient is independent of r̂ due to the average over the orientations,
where we can select k̂ = r̂ [cf.eqn (94) and (96)] and where Q_i is defined in eqn (98). The Yukawa function factor arises from 1/[small epsilon, Greek, tilde](k) = k²ε_r[small phi, Greek, tilde]_Coul*(k) evaluated in r space. By defining the average entitieswhere x^b_i = n^b_i/n^b_tot is the mole fraction of species i and where we have used eqn (97) to obtain the last equality, we can write these equations aswhen r → ∞. While Q̄^eff_N is simply the average over all particles of the orientation independent part of Q^eff_i [eqn (95)], Q̄^eff_Q contains a weighing with respect to Q_i that depends on the bare charge density σ_i of the particles [see eqn (98)].

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 κ² = βn^b_totq_eQ̄^eff_Q(κ)/ε₀, 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 Q̄^eff_N ≡ 0 and H_NQ(r) ≡ 0, while the decay of H_NN(r) can be obtained from

as follows from eqn (122). H_QQ(r), on the other hand, is given by the same expression as in the main case, eqn (128). In this case the behaviors of H_NN(r) and H_QQ(r) are independent of each other and the decay length of H_QQ(r) remains finite even in situations where the decay length of H_NN(r) becomes very large. This constitutes a considerable difference compared to the general case.

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 |Q̄^eff_N| and |Q̄^eff_Q| varies for the different modes, so the contributions with some decay lengths may have |Q̄^eff_N| > |Q̄^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 h̃_ij: the poles of H̃_NN(k) [poles of h̃_ij*] and those of H̃_QQ(k) [zeros of [small epsilon, Greek, tilde](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 h̃_ij remain in the sense that each pole in one set is close to a pole of h̃_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 [small epsilon, Greek, tilde](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 |Q̄^eff_N| ≫ |Q̄^eff_Q| for the first set and |Q̄^eff_N| ≪ |Q̄^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 Q̄^eff_i(κ) ≈ q_i, which implies that n^b_totQ̄^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_totQ̄^eff_N(κ)]², is very small. H_QQ(r) is much larger because its coefficient is proportional to the square of

where all terms in the sum are positive.

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πε₀ε_rr) which originates from the leading term e^−κr/(4πε₀[scr E, script letter E]^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.

5 Summary and concluding remarks

In this section we will give an overview of the key results and also give some perspectives on their usage in experimental and theoretical investigations.

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 [scr E, script letter E]^eff_r(κ_ν), which determines the magnitude of the mode's contribution to ϕ_Coul*(r) [eqn (86)] and, for oscillatory modes, determines also the phase. [scr E, script letter E]^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 [scr E, script letter E]^eff_r(κ_ν) reflects the polarization response of the electrolyte to an exponentially decaying disturbance with decay parameter κ_ν.

Furthermore, the dielectric factor [scr E, script letter E]_r*(κ_ν) that appears in the expression for the screening parameter κ_ν, eqn (5), has different values for the various modes. As we have seen, [scr E, script letter E]_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 [scr E, script letter E]^eff_r(κ_ν) and [scr E, script letter E]_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 [small epsilon, Greek, tilde](k); in Fourier space we have [small phi, Greek, tilde]_Coul*(k) = [small phi, Greek, tilde]_Coul(k)/[small epsilon, Greek, tilde](k), where [small phi, Greek, tilde]_Coul(k) = 1/(ε₀k²) is the Fourier transform of the ordinary unscreened Coulomb potential ϕ_Coul(r). The decay modes of ϕ_Coul*(r) are given by the poles of [small phi, Greek, tilde]_Coul*(k) in complex k-space, i.e., the zeros of [small epsilon, Greek, tilde](k), whereby the screening parameter κ satisfies [small epsilon, Greek, tilde](iκ) = 0, where the imaginary unit i appears because the dielectric response is evaluated for an exponentially decaying field. The condition [small epsilon, Greek, tilde](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 [small epsilon, Greek, tilde](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 [scr E, script letter E]_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 ψ_ivia Coulomb's law. The mean potential ψ_i can be obtained from ρ^tot_i by using the usual unscreened Coulomb potential, ϕ_Coul(r) = 1/(4πε₀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 [small epsilon, Greek, tilde](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 [scr E, script letter E]^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_iQ^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 Q̄^eff_N and Q̄^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 (Q̄^eff_N)², those of H_NQ(r) proportional to Q̄^eff_NQ̄^eff_Q and those of H_QQ(r) proportional to (Q̄^eff_Q)² [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 h̃_ij*, which are different from the poles of [small phi, Greek, tilde]_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.³⁹ 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 [small phi, Greek, tilde]_Coul*(k) [i.e., zeros of [small epsilon, Greek, tilde](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 h̃_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 [small phi, Greek, tilde]_Coul*(k) vanishes. The decay mode remains but is instead determined by the pole of h̃_ij* for the invariant case. Each pole of [small phi, Greek, tilde]_Coul*(k) in the second set continues to be such a pole throughout. For the decay modes given by first set of poles |Q̄^eff_N| ≫ |Q̄^eff_Q| and Q̄^eff_Q becomes zero when the system becomes charge-inversion invariant, while for the second set |Q̄^eff_N| ≪ |Q̄^eff_Q| and Q̄^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.

Conflicts of interest

There are no conflicts to declare.

Appendices

A Response due to variation in external electrostatic potential for the nonspherical case

For an inhomogeneous fluid of nonspherical particles Yvon's first equation is
and for a bulk fluid perturbed by the weak external potential δv_j(R) ≡ δv_j(r,ω) this becomes
When the external potential is solely due to the electrostatic potential δΨ^ext(r) we have
and hence we obtainWe now introduce the charge density ρ_i(r₃|R₁) at r₃ around a particle with coordinates R₁ (i.e., it is located at r₁ and has orientation ω₁)where we have obtained the last equality by using that follows from electroneutrality. Note that , where σ_j(r₂₃,ω₂) ≡ σ_j(r₃|r₂,ω₂) for given ω₂ is a function of r₂₃ only. Thus we can write eqn (130) as
where ρ^tot_i is the total charge density of the particle and its surroundings. The change in charge density, δρ(r₄) at a point r₄, is hence given bywhere H_QQ is the charge–charge correlation function given in eqn (134).

B The H_QQ, H_NN and H_NQ correlation functions

The charge–charge correlation function H_QQ is defined aswhich can be written asH_QQ(r₁,r₂) describes the correlations in fluctuations in charge densities at r₁ and at r₂. The first term on the rhs of eqn (133) originates from the self-correlations of the particles. Using the notation introduced in eqn (47) we can write eqn (134) asand in Fourier space this iswhich is independent of k̂ = k/k.

The density–density correlation function H_NN is defined as

where we have written h_ij(R₁,R₂) = h_ij(r₁₂,ω₁,ω₂). H_NN(r₁₂) describes the correlations in fluctuations in total number densities at r₁ and at r₂. Its Fourier transform iswhere .

Finally, the density–charge correlation function H_NQ ≡ H_QN is defined as

It describes the correlations in fluctuations in total number density at r₁ and the charge density at r₂ and vice versa. In Fourier space this is

C Relationships between the poles of h̃_ij* and h̃_ij

In this appendix we will investigate the relationships between the poles of h̃_ij and h̃_ij*. In particular we will show that each pole of h̃_ij*(k) in eqn (113) is normally cancelled by a corresponding pole in the last term of this equation. This implies that h̃_ij and h̃_ij* cannot have poles for the same k values (apart from some exceptional cases treated in detail in Section 4.2). We will also investigate some cases where poles of h̃_ij* can influence where the poles of h̃_ij appear.

The existence of a simple pole for h̃_ij* at k = ±iα means that h̃_ij* diverges like

where F̃_ij(kk̂,ω₁,ω₂) is a function that is finite at k = ±iα. One can show that F̃_ij factorizes§§§ at k = iα as , where τ_l for l = i, j is a function of orientation for a given α and k̂ and where K = K(α) is a constant, and we havewhere is a part of h̃_ij* that stays finite at k = iα. This means that this pole gives a term in h_ij* in ordinary space that decays as
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

whereand is the part of [small rho, Greek, tilde]_i* originating from . E(k,α) does not depend on the direction of k̂ due to the averaging over ω. The dielectric function [small epsilon, Greek, tilde](k) [given by eqn (50)] can be written as the first line in the following equation and by inserting eqn (142) we obtain the second linewhere [small epsilon, Greek, tilde]^rest is the part of [small epsilon, Greek, tilde](k) originating from . Since E(k,α)|_k=iα is normally nonzero, the functions h̃_ij*, [small rho, Greek, tilde]_i* and [small epsilon, Greek, tilde](k) diverge when k → iα due to k² + α² in the denominator. The remaining terms in the expressions stay finite for k = iα.

We now insert eqn (142) and (144) into the last term in eqn (113) and since [small rho, Greek, tilde]_i*, [small rho, Greek, tilde]_j* and [small epsilon, Greek, tilde](k) are dominated by the diverging contributions for k values close to iα, it follows that we have when k → iα

so on the rhs of eqn (113) this pole cancels the pole of h̃_ij* given by the first term on the rhs of eqn (141). This means that provided E(k̂,α) is nonzero, the only poles of the entire rhs of eqn (113) are those that originate from the zeros of [small epsilon, Greek, tilde](k) in the denominator.

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 h̃_ij* at k = iα is not cancelled in h̃_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 h̃_ij are given by the zeros of [small epsilon, Greek, tilde](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 h̃_ij* and those that zeros of [small epsilon, Greek, tilde](k). When the system deviates only slightly from being charge-inversion invariant and all poles of h̃_ij are given by the zeros of [small epsilon, Greek, tilde](k), the former group must have turned into zeros of [small epsilon, Greek, tilde](k) because the poles of h̃_ij* are not poles of h̃_ij. This is indeed the case as we now are going to see.

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

Provided that E(iκ_ν,α) ≠ 0 we can write this asWhen E²(iκ_ν,α) ≈ 0 this equation has a solution κ_ν ≈ ±α with the property κ_ν → ±α when E(iκ_ν,α) → 0. This is the solution that we are looking for. Its existence can also be inferred directly from eqn (145) because E²(iκ_ν,α)/(α² − κ_ν²) must remain finite in this limit since [small epsilon, Greek, tilde]^rest(iκ_ν) stays finite there. Likewise, there exist solutions of [small epsilon, Greek, tilde](k) = 0 that are connected in the same manner to other poles of h̃_ij*. The remaining solutions κ_v′ of eqn (146) are not interesting here because in the same limit they go to the solutions of 1 + [small epsilon, Greek, tilde]^rest(iκ_ν′) = 0, which is the same as [small epsilon, Greek, tilde](iκ_ν′) = 0 as seen from eqn (144) with E(k,α) = 0. Thus they belong to the second group mentioned above for charge-inversion invariant systems (these solutions are poles of h̃_D,AB in eqn (125) and not poles of h̃_ij*).

References

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

---

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 κ²/κ_D² = 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_j = q_j/q_e is the valency of species j and η_j is the stoichiometric coefficient. The term proportional to Λ, which was originally derived in ref. 22, vanishes for symmetric binary electrolytes. Both terms originate solely from the long-range, purely electrostatic correlations, while contributions that depend on other properties of the ions (like their sizes) affect the higher order terms, that here are included in [scr O, script letter O](Λ²) where Λ² is proportional to the total ion density.

¶ 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 = ε₀χ_eE 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 [small phi, Greek, tilde]_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² 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,³⁵ 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) = {h̃_ij*(k,ω₁,ω₂)} is the solution of the Ornstein–Zernike (OZ) equation (N + NH*N) ★ (N⁻¹ − C*) = 1, where C*(k) = {c̃_ij*(k,ω₁,ω₂)} 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⁻¹ − C*)⁻¹ = Adj(N⁻¹ − C*)/D*, where Adj denotes the adjugate matrix and D* = Det(N⁻¹ − 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.

---