Solvent-induced ion clusters generate long-ranged double-layer forces at high ionic strengths

Abstract

Recent experimental results by the surface force apparatus (SFA) have identified a dramatic deviation from previously established theories of simple electrolytes. This deviation, referred to as anomalous underscreening, suggests that the range of electrostatic interactions increase upon a further addition of salt, beyond some threshold concentration (usually about 1 M). In this theoretical work, we explore an extension of the restricted primitive model (RPM) wherein a short-ranged pair potential of mean force (sPMF) is added to the usual Coulombic interactions so as to mimic changes of the hydration as two ions approach one another. The strength of this potential is adjusted so that the modified RPM saturates at a realistic concentration level (within a range 4–7 M, typical to aqueous 1 : 1 salts). We utilise grand canonical simulations to establish surface forces predicted by the model and compare them directly with SFA data. We explore different sPMF models, which in all cases display significant clustering at concentrations above about 1 M. In these models, we find significant double-layer repulsion at separations that significantly exceed those expected from standard RPM predictions. We do not, however, observe an increase of the screening length with salt concentration, but rather that this screening length seemingly saturates at a (rather high) value. The simulated long-ranged interactions are shown to correlate with ion cluster formation, implicating the important role of accompanying cluster–cluster interactions. In particular, steric interactions between clusters (manifested in density–density correlations) are quite relevant in these systems.

---

1 Introduction

The theoretical study of interactions between large charged particles is important in almost all areas of colloid science.^1–3 For some time there had been a general consensus that the behaviour of charged particles in aqueous solutions of monovalent (1 : 1) electrolytes could be described by so-called primitive models (PM), at least qualitatively. In these models the aqueous solvent is implicitly treated, only asserting itself via a uniform relative dielectric constant, ε_r. The ionic species, on the other hand, are explicitly accounted for as charged hard spheres. In a mean-field treatment that ignores ion–ion correlations, this description can be even further simplified by neglecting the hard core interactions between the background salt ions. Meanwhile, colloidal particles can be treated as charged surfaces. Assuming additive van der Waals forces (of quantum origin) between the surfaces results in the famous DLVO theory,^4,5 which has proven very useful for the study of the stability of colloidal suspensions.

One generally finds that mean-field (MF) treatments of the PM predict that, at low electrolyte concentrations, the dominant non-quantum interaction between charged surfaces is electrostatic in origin and is exponential in form, with a range largely dictated by the Debye screening length, λ_D. This interaction arises from the direct forces between the charged surfaces that are screened by a counter-charge density sourced from the background salt. The Debye length measures the typical spatial length over which non-electroneutral fluctuations occur in the electrolyte solution. Structural variations in the overall density (rather than charge density) are also induced by the surfaces, but the forces that arise from this effect are usually negligible at large surface separations.^6,7 The Debye length decreases with electrolyte concentration, c, according to, . Corrections to MF theories require the inclusion of direct (short-ranged) ion correlations, which lead to minor modifications of the Debye length. For example, taking into account hard core interactions leads to an excluded volume corrected ‘effective’ value λ_eff, which is smaller than λ_D.^6,8 The qualitative correctness of MF theories using the PM, and the role played by the Debye length, is reasonably well-supported at low electrolyte concentrations by experiments, as well as simulations.¹

MF theories break down at high concentration where the Debye length becomes small. In this regime, ion correlation corrections must be included to provide even a qualitative description of the PM. Such corrections predict that the surface interactions of electrostatic origin will switch from exponential to oscillatory, as the Debye length approaches the diameter scale of the salt particles. This is the so-called Kirkwood transition.⁹ This transition has been confirmed by simulations of the PM and is due to the monotonically decreasing (essentially exponential) ionic double-layer becoming unstable to the formation of charged layers by the electrolyte.¹⁰ In this regime, the non-electrostatic structural forces may compete with the electrostatic interactions for relative importance as both are oscillatory and of comparable range.^6,7 However, recent experimental measurements using the surface force apparatus (SFA) have challenged this result.^11–14 That is, instead of transitioning to an oscillatory function, the asymptotic interaction between charged surfaces remains exponential, but with a screening length that suddenly begins increasing (rather than decreasing) with concentration, eventually becoming significantly larger than the nominal Debye length. Surprisingly, this occurs even in aqueous 1 : 1 salt solutions, where it was generally believed that the PM should be reasonably reliable. This phenomenon has been described as anomalous underscreening. One consequence of this is that suspensions of like-charged colloids may remain stable at high ionic strength due to significant repulsive interactions at colloidal separations where the attractive van der Waals interactions remain small. This is counter to the predictions of traditional DLVO theory. Despite considerable theoretical efforts,^15–21 there is at present no consensus as to the physical origin of interactions at such large inter-particle distances. One can find independent (non-SFA) experimental evidence for anomalous underscreening in literature.^22–25 It should also be noted that there is at least one recent experimental work that does not support this phenomenon.²⁶ Kumar et al. used an atomic force microscope to measure interactions between charged surfaces (mainly silica-silica) in aqueous salt solutions. They did not find any sign of anomalous underscreening, i.e. at high ionic strengths, the measured surface forces were quite short-ranged.

The formation of ionic clusters of low overall charge has been proposed as a possible mechanism for the putative increase in the Debye length observed in the underscreening concentration regime. Such a model implicitly assumes that the electrostatic interactions between the surfaces remains dominant and the formation of clusters will effectively reduce the concentration of screening charges.^12,27 For this to be a feasible explanation, clustering needs to be sufficient to reduce the free charge concentration by a factor of up to ∼10⁴. Recent simulations of the restricted primitive model (RPM), wherein ions have equal sizes, have investigated the degree of clustering that occurs in a 1 : 1 electrolyte.²¹ That work also took into account dielectric saturation by allowing the relative dielectric constant of the solvent, ε_r, to become concentration-dependent. In particular ε_r decreased with c to account for the diminished dielectric response of rotationally constrained water molecules in the ionic hydration shells.^28–34 It was seen that, while increased clustering was evident at high concentrations, and lower ε_r, it was not of a sufficient magnitude to explain anomalous underscreening.²¹ That is, the RPM (even accounting for dielectric saturation) would predict that forces between like-charged surfaces would be essentially completely screened by salt ions at the large distances where experiments still find a significant repulsion.

This failure of the RPM led to a recent work³⁵ in which we proposed a localised modification of the RPM by addition of an additional short-ranged potential of mean force (sPMF) between ions. The inclusion of the short-ranged potential induced significant ion clustering (as confirmed via simulations). In the present article, we study the consequences of the formed clusters on the interactions between two charged surfaces. These interactions are directly compared with experimental SFA data, and it is demonstrated that cluster formation indeed leads to remarkably long-ranged surface forces. We provide analyses of the separate contributions to the correlations in these systems in order to uncover the dominant physical mechanism leading to these long-ranged forces.

2 Model and methods

A weakness of the RPM is that it does not (with an aqueous solvent) saturate at reasonable concentrations. In this work, we explore the hypothesis that what the RPM is lacking is a short-ranged potential of mean force (sPMF), adjusted to induce solution instability at a concentration that is close to typical saturation values (4–7 M) for simple salts. While we do not have a detailed knowledge of this potential of mean force, it is plausibly a result of water restructuring as the ionic environment becomes crowded. In previous work,³⁵ we took an approach whereby we assigned the sPMF a range of about one water diameter (3 Å). Since our simulations are based on Coulomb interactions, a convenient option was to increase the strength of the Coulomb interaction in this separation range. This is achieved by varying the dielectric constant using linear extrapolation from the bulk down to a smaller contact value. Of course, this choice of additional short-ranged interaction is not unique and in this paper we will consider different functional forms of the sPMF, including the model we used earlier.³⁵ Due to its relatively narrow range, we will dub that model the “narrow” sPMF.

2.1 The “narrow” sPMF, ϕ_n

As described above, the explicit form for this sPMF uses the following linear ramp function for the spatial variation of the dielectric constant, ε_r, within its range,Here, d is the hard-sphere diameter of the ions, while Δ is a measure of the thickness of a hydration shell and chosen to be Δ = 3 Å, the diameter of a water molecule. The parameters, ε_c and ε_b are the contact and bulk values of the dielectric factor, respectively. At room temperature, ε_b = 78.3 and we have chosen, ε_c < ε_b. It could be argued that this inequality reflects the exclusion of water solvent as ions approach each other, together with the strengthening of the electric field leading to dielectric saturation. However, as we discuss later, these effects are really many-body in character and can only qualitatively be reproduced via an sPMF, which is a 2-body interaction. More pragmatically, this approach can be viewed as merely a convenient way to introduce a solvent-induced PMF at short range which promotes solution saturation and does not mean that the dielectric factor actually varies in this way. Indeed, we would argue that since the effects we attempt to incorporate are effective at the molecular scale, they are not describable by continuum electrostatics.

With this choice for the dielectric factor, the sPMF between ions i and j of valency |Z_i| = |Z_j| = 1 (ϕ^ij_n) can be written as:

where l_B(r) = βe₀²/(4πε₀ε_r(r)) and l_B(bulk) = βe₀²/(4πε₀ε_b). It should be noted that the nature of ϕ^ij_n is given by the valency of the interacting species, i.e. it is always attractive between unlike charges, and repulsive between charges of equal sign.

For a given choice of d, we can tune ε_c to a value where the solution phase separates at a concentration typical of the saturation limit for an aqueous monovalent electrolyte solution (about 4–7 M, at room temperature). In practice, we chose d = 3 Å and ε_c = 23 where the solution is stable at 3.45 M, but will phase separate at a marginally smaller value for ε_c. Our rationale for choosing these values in given in ref. 35. We also consider d = 4 Å, in which case ε_c = 20 will give a similar stability. We will use this as the main model of investigation in this paper, including looking at the effect of salt concentration on the interaction between charged surfaces.

We will also perform a more limited study on a so-called “wide” version of an sPMF, where we will also explore other different combinations of attractions and repulsions between ionic species. The details of this interaction are given below.

2.2 The “wide” sPMF, ϕ_w

Here, we adopt a sPMF with the same cut-off (Δ) as ϕ^ij_n, but with a broader range over which it maintains a substantial strength. Specifically, we define it as,One important difference, compared with the narrow version, is that the sign of A_ij can be arbitrarily assigned. This leads to various sPMF:s with a qualitatively different impact. The following nomenclature describes the choice of parameter values we used in this study. For instance, ϕ_w(a, r) indicates that A₊₋ is positive between unlike charges (attractive sPMF), but that A₊₋ = −A₋₋ = −A₊₊ = 2.5. The corresponding notation for an overall attractive sPMF is ϕ_w(a, a), for which A₊₊ = A₋₋ = A₊₋ = 0.5. Finally, ϕ_w(a, 0) signifies an attractive sPMF between unlike charges, but with no impact on the interaction between ions of the same sign, i.e. A₊₊ = A₋₋ = 0, and A₊₋ = 1.0. In all cases, the amplitudes used were regulated to be strong enough to “almost” induce phase separation in a bulk solution at 3.45 M. We have only investigated ϕ_w systems in which d = 3 Å.

The overall qualitative difference between ϕ_n(r) and ϕ_w(r) is illustrated in Fig. 1. The total potential of mean force between ions i and j, ϕ_ij(r) (separated a distance r) can thus be written as:

where α = w or n, depending on the choice of sPMF. Importantly, for r ≥ d + Δ, the pair potential is equivalent to the RPM pair potential, corroborating the local nature of our modification.

Fig. 1 A comparison between ϕ_n(r) and ϕ_w(r). The illustration is made for the +− sPMF, which is always attractive, and for d = 3 Å. Moreover, for the narrow sPMF, we have set ε_c = 23, whereas the wide sPMF is based on an amplitude A₊₋ = 2.5.

2.3 Surface force simulation details

Our simulation system is comprised two parallel flat and uniformly charged surfaces, separated at a distance H, and extending infinitely along the x and y directions. We imagine these surfaces to be immersed in a salt solution. The bulk region external to the space between the surfaces is characterised by the electro-chemical potential and the bulk osmotic pressure, p_b. A grand canonical simulation scheme ensured that the simulated fluid, contained between the surfaces, was in equilibrium with the bulk and that the surface charges were properly neutralised.³⁶ The surface charge density was chosen to be σ_s = 1/70 ≈ −0.014336 e Å⁻², which is a value typical to mica,³⁷i.e. the surfaces used in SFA experiments. All simulations were performed at 298 K. Periodic boundary conditions were applied along the (x, y) directions, and the “charged sheet” method¹⁰ was adopted to account for long-ranged interactions. Cluster moves were implemented, leading to crucial improvements of the statistical performance. The internal (osmotic) pressure, p, acting transverse to the walls was evaluated from the average normal force across the walls, separated by the expected distance between Gibbs dividing (solid–fluid) surfaces, h ≡ H−d. Subtracting p_b, we obtain the net pressure, p_net = p − p_b. Integrating this pressure, and utilising the Derjaguin approximationallows us to construct force per radius curves, F/R, facilitating direct comparisons with SFA data (R is the radius of the curved surface used in the experimental setup). To make comparisons with the SFA data, the simulated interaction curves were matched to the experimental results for a 2 M NaCl solution at a separation of 57 Å.

The non-electrostatic ion–wall interaction was modelled by an analytic and purely repulsive wall potential. We explored two different choices, decaying as δ⁻⁴ and δ⁻⁶, respectively, where δ is the transverse distance between the ion and the wall. The first model was adopted when investigating the effect of salt concentration, using ϕ_n as the sPMF. One might argue that this softer choice amounts to a coarse-grained description of molecularly rough surfaces. In another part of this study, which explores the various choices for the sPMF in concentrated solutions (1.6–2 M), we chose the steeper wall, decaying as δ⁻⁶. This represents a better model for a molecularly smooth surface, such as mica. Details are provided in the ESI,† including a comparison between results obtained with either wall, for a given system. The ESI,† also contains a more thorough description of our model and simulation methods.

3 Results and discussion

Our results comprise two parts. In the first part, we evaluate how surface forces respond to salt concentration changes, using the “soft” wall (δ⁻⁴) description (see ESI†), and ϕ_n as the sPMF. The second part compares surface interactions at a high salt concentration, using various choices for the sPMF. In that part, we have employed the “steep” wall (δ⁻⁶) description, which arguably is a better choice for mica surfaces. On the other hand, as we discuss in the ESI,† the difference in wall softness essentially amounts to a small shift (about 3 Å) of the surface separation.

3.1 Salt concentration effects, using ϕ_n

In previous work we performed simulations using the “narrow” ϕ_n potential as the sPMF in addition to the standard RPM interaction for bulk electrolyte solutions.³⁵ We will denote this as the sRPM(ϕ_n) model The findings important to the current investigation can be summarised as follows:

1. At low electrolyte concentration, charge–charge correlations displayed behaviours typical of classic MF theory, and the Debye length is the relevant length-scale.

2. At higher concentration (beyond 1 M) there was evidence of strong cluster formation, much larger than predicted by the RPM. Charge–charge correlations became short-ranged and oscillatory (the Kirkwood transition), while density–density correlations grew to become more dominant, displaying a monotonic decay, with a length-scale much longer than the Debye length.

3.1.1 Simulated surface forces. The presence of large clusters in the sRPM(ϕ_n) model prompted us to consider their consequence on the interaction between charged surfaces. Specifically, we report here simulations of charged surfaces immersed in the sRPM(ϕ_n) solution and compare our results with experimentally measured forces that show anomalous underscreening. Simulated net pressures acting on the charged surfaces for a range of different concentrations of the sRPM(ϕ_n) are shown in Fig. 2. We observe a significant repulsion for the simulated curves over a large range of surface–surface separations. The range of the repulsion initially decreases as the ionic strength increases, but then appears to level out and become essentially concentration-independent above some threshold value (1.6 M). For comparison, we have also included results for the 1.6 M RPM electrolyte in which the relative dielectric constant was uniform and set to ε_r = 23. We found no detectable net pressure for the RPM electrolyte at any separation. A short-ranged sPMF can induce clustering more efficiently than a global increase of the electrostatic coupling. This can be explained as follows. Consider linear cluster of alternating charges: +−+−+−. A global increase of the Debye length will generate a significant repulsion between next-nearest neighbours, but this is not so with a short-ranged sPMF that is attractive between ions of opposite charge but repulsive when the ions share the same valency. Similar results would also be obtained for an RPM model with ε_r = 78.3. That is, the sRPM(ϕ_n) predicts surface forces of a much longer range than the traditional RPM, even if a broad a range of dielectric constants is used in the latter. From our previous work, we found that the sRPM(ϕ_n) predicts a much higher degree of clustering than the RPM, again despite the value of ε_r. It follows that the much longer-ranged forces observed in the sRPM(ϕ_n) in Fig. 2 is related to this increased cluster formation. The F/R predictions that follow from the net pressures are given in Fig. 3. These are compared with SFA data on NaCl_(aq), at a concentration of 2 M. The simulation predictions at close to 2 M agree reasonably well with experiments for separations less than about 40 Å. To put this agreement into context, we note that the F/R curves predicted by the RPM would be essentially zero if included in Fig. 3. Given the discussion above, this prompts us to assert that the large range of the forces seen in experiments is due to significant ionic clustering in the electrolyte. This assertion implies that the relevant mechanism is to be found in the behaviour of the solution phase and is not, for example, a surface phenomenon. Moreover, given our previous findings of the dominance of density–density correlation in the high density regime (beyond 1 M), it is anticipated that the surface forces are due to structural changes in the overall density, rather than the charge density, as the surfaces approach. That is, clusters will adsorb onto the surfaces, primarily due to charge and multipole interactions. The adsorbed layers will begin to breakdown as the surfaces encroach each other, leading to a net repulsion with a length-scale of the typical cluster size. We expect that the average cluster size should in principle grow with the concentration of the electrolyte solution, perhaps explaining the apparent growth in screening length observed in anomalous underscreening. This notwithstanding, the sRPM(ϕ_n) predictions are not completely consistent with experiments. Fig. 3 indicates that at large concentrations the range of the interactions appears to reach a plateau where F/R becomes concentration-independent. This disagrees with the “anomalous underscreening” predictions by SFA, according to which there is an increase of the decay length as the ionic strength increases in the high concentration regime.

Fig. 2 Results from sRPM(ϕ_n) grand canonical simulations of net pressures, at various bulk salt concentrations, with d = 3 Å, ε_c = 23 and ε_b = 78.3. An exception is the dashed line, which displays results using the standard RPM, with a uniform ε_r = 23. The inset is a zoom-in at large separations.

Fig. 3 Simulated prediction of F/R, with d = 3 Å and d = 4 Å, at various salt concentrations, compared with SFA data with a 2 M NaCl_(aq) solution, obtained from prof. Susan Perkin (data previously published in ref. 13). Simulated curves have been shifted to roughly match the SFA data, at the largest investigated separation. (a) d = 3 Å (linear scale). (b) d = 3 Å (log scale). (c) d = 4 Å (linear scale).

For completeness, we include results obtained with a larger value for the hard-sphere diameter of the ions, d = 4 Å. Recall also that we chose ε_c = 20 at this diameter, as described in detail in ref. 35. Comparisons with raw SFA data are given in Fig. 3(c). We note that the difference between the interactions obtained at 0.40 M and 1.6 M is considerably smaller than we find by a similar comparison (0.50 M and 1.6 M), with d = 3 Å.

We also compare the simulated net pressures directly with a derivative of the experimentally observed free energies per unit area (taking into account an appropriate scaling factor), Fig. 4. In the latter case, we needed to smooth out the raw data by performing “running averages” prior to taking the derivative (via discrete difference). The extremely long-ranged tail found by the SFA is perhaps even more apparent in this presentation, and it appears to be an aspect of the experimental data that the sRPM(ϕ_n) is unable to reproduce.

Fig. 4 Simulated net pressures, and corresponding predictions from SFA data. The inset is a zoom-in at large separations.

3.1.2 Structural analyses. Here we illustrate how the long range forces are related to a slow density decay of the ionic solution, which results from cluster formation. We will only consider systems with d = 3 Å. We start by analysing ion concentration profiles, n(z), along the direction normal to the surfaces. In Fig. 5(a), we note how the ion concentrations rapidly approach the bulk value in RPM simulations. However, using the sRPM(ϕ_n), this approach is quite slow, and even at a surface separation of about 5 nm, both cation and anion mid plane concentrations are greater the bulk value (recall that these results are from grand canonical simulations). This is because a significant fraction of the anions are members of a cluster, many of which are expected to be net positively charged and hence are adsorbed to the negative surface. This leads to an anion concentration profile that approaches the mid plane from above. By comparison, the RPM would produce an anionic profile approaching its mid plane value from below.

Fig. 5 Ion concentration profiles along z. Solid lines depict cation density profiles, whereas the dashed lines show anion profiles. The dotted line indicates the bulk concentration. (a) Comparing concentration profiles at h = 48 Å, as obtained with the sRPM(ϕ_n) and the standard RPM. (b) Concentration profiles from sRPM(ϕ_n) simulations, at two different separations: h = 48 Å and h = 57 Å.

In Fig. 5(b) we see how the mid plane values of the ion profiles gradually decreases towards the bulk concentration, but we note that a significant difference persists even at quite large surface separations.

We also investigated the correlations parallel the surface, focusing on the region at the mid plane. We computed lateral (2D) species resolved pair correlation functions (presented in detail in the ESI†) for low (0.13 M) and high (3.4 M) concentrations, for all slit widths under investigation. No evidence of a continuous solid phase was found in any of these systems indicating the absence of a frozen state. In order to further investigate the dominant correlations in our system, we computed the 2D charge–charge, h_cc(ρ), and density–density, h_nn(ρ), total correlation functions from the species resolved pair correlation functions.⁷ Note that the total correlation functions are related to the pair distribution functions via h(ρ) = g(ρ) − 1. For brevity, the explicit dependence is omitted.

This allows us to essentially deconstruct the mean like-charge total correlation functions and observe the correlations dominating at the asymptotic limit via. The asymptotic behaviour of correlation functions, as r → ∞, is described by a Yukawa decay, or an oscillatory Yukawa decay⁷where μν ∈ {c, n, +, −} denotes the corresponding correlation type μν in a general sense. We define 1/ξ^(μν)₀ as the asymptotic decay correlation length (the effective screening length) and 2π/ξ^(μν)₁ as the asymptotic frequency correlation length, which describes the spatial range of electrolyte charged layer oscillations appearing after the Kirkwood transition. A^(μν) is used for the amplitude and Θ^(μν) for the phase shift. Interested readers are referred to ref. 6, 7 and 38 for detailed derivations. The results for the c ≈ 3.4 M system, obtained at the mid plane, is presented on Fig. 6. Further details including all results from various slit widths under investigation is presented in the ESI.†Fig. 6(a) and (b) presents the density–density and charge–charge total correlation functions on a linear scale for two different surface separations, h = 21 Å (top) and h = 57 Å (bottom). For the large surface separation we observe apparent monotonic decay of the density–density correlations at long range and a damped oscillatory decay for the charge–charge correlations for both surface separations. In order to better observe the long-ranged asymptotic behaviour, we plot ρ|h_μμ(ρ)|, μμ ∈ {nn, cc} using a log-linear scale in Fig. 6(c) and (d). Here we observe two major changes in the correlation functions as the surface separation increases. Firstly at the smaller surface separations the charge–charge correlations display a higher amplitude compared to larger separations. Secondly, the density–density correlations decay faster at the shorter separation, with a much longer range decay at larger separations. This is consistent with the cluster picture described above. Namely at shorter separations, typical clusters will be smaller with highly correlated charge–charge interactions, corresponding to more tightly packed ions. At larger separations, much larger clusters can be accommodated between the surfaces. Such clusters will have a more loose charge arrangement. It is probably the exclusion of larger clusters as the surface separation decreases that leads to the repulsion between the surfaces. Indeed, this is a major conclusion of this paper. This notwithstanding, our results have by no means proven that substantial ion clusters do indeed exist in reality for such systems. Our sRPM(ϕ_n) potential approximates an effect due to the interaction between salt ions and the solvent, which should in principle be modelled by many-body interactions in a solute-only approach. It should be noted that some all-atomistic molecular dynamics simulations do suggest that ionic clustering can be quite significant in the presence of explicit water models.^39,40 On the other hand, such findings are quite dependent upon the potential model employed.^41,42

Fig. 6 Correlation function analysis in the slit, obtained for the c ≈ 3.4 M sRPM(ϕ_n) system at the mid plane, for h = 21 Å and h = 57 Å respectively. Here, red is used for the charge–charge and blue for the density–density total correlation functions, while grey denotes the mean like-charge total correlation function. Note: h_mean = h_cc + h_nn. Graphs (a) and (b) present the 2D charge–charge, h_cc, and density–density, h_nn, total correlation functions on a linear scale. Graphs (c) and (d) provide the asymptotic analysis plots on a logarithmic scale.

In order to complement the investigation of lateral correlation functions obtained at the mid plane, we present a comparison with bulk correlation functions on Fig. 7, computed with the method described in our previous work.³⁵ We explicitly investigate three extreme cases: the smallest h = 21 Å surface separation, largest h = 57 Å, and a neutral wall h = 77 Å, by directly comparing with bulk correlation functions for a system of equal concentration. We plot all correlation functions on a linear scale on Fig. 7(a). We can observe a minimal difference between the full (h = 57 Å) and dashed (h = 77 Å neutral walls) lines, insert on Fig. 7(a). All of the correlations in the slit have amplitudes lower than the bulk. Fig. 7(b)–(d) demonstrate the similarity of the frequency correlation length 2π/ξ^(cc)₁ and decay correlation length 1/ξ^(cc)₀ between the bulk and slit systems. However, in both cases the amplitude of charge–charge correlations is smaller than the corresponding bulk value, demonstrated as a downwards shift on the y-axis. We also observe a slight phase shift between the bulk and slit system charge–charge correlations. The density–density correlations have a smaller decay correlation length for the short surface separation Fig. 7(b), and larger for the wider surface separation Fig. 7(c), when compared to the bulk values. In both cases, the amplitude is again smaller for the slit system. Fig. 7(d) demonstrates the correlation analysis for a wide surface separation with neutral walls. We can see that the difference between Fig. 7(c) and (d) is minimal, demonstrating that the charge of the surfaces has no effect on the lateral correlation decay at the mid plane at large surface separations. Importantly, the similarity of frequency and decay charge–charge correlation lengths between slit and bulk results at all surface separations (albeit with moderately large differences in amplitudes and a slight phase shift) demonstrate that the asymptotic behaviour lateral to the confining charged surfaces is dictated by bulk properties of the ionic fluid and not by the surfaces themselves. For completeness, the same analysis is also presented for the low concentration system in the ESI.† The ESI,† also contains illustrative configurational snapshots, at various separations, and bulk concentrations.

Fig. 7 Comparisons of correlation functions obtained for the c ≈ 3.4 M sRPM(ϕ_n) system at the mid plane for h = 21 Å, h = 57 Å charged wall systems, as well as the h = 77 Å neutral wall (n.w.) system, with bulk correlation functions. For the bulk, all correlation functions have a standard radial dependence i.e. h_μμ ≡ f(r). Subplots compare the bulk correlation functions with (a) all the correlation functions on a linear scale, (b) h = 21 Å slit correlation functions, (c) h = 57 Å slit correlation functions, and (d) h = 77 Å with neutral walls slit correlation functions.

3.2 Comparison of different sPMF models at high salt concentrations

Here, we will compare results obtained with the sRPM(ϕ_n), and the sRPM(ϕ_w) models (where the latter uses the wider form of the sPMF in addition to the RPM interactions). Recall that with this latter model we have different combinations of attractions and repulsions between ionic species. To reiterate, they are denoted as, ϕ_w(a, r) (attractive +−, repulsive ++, −−), ϕ_w(a, a) (attractive +−, attractive ++, −−) and ϕ_w(a, 0) (attractive +−, no sPMF between ++ and −−). In all cases below, we have set d = 3 Å, and used the “steep” wall description. Bulk salt concentrations were adjusted (via the chemical potential) to lie in the regime 1.6–2 M.

Net pressure curves in concentrated samples (1.6–2 M), for various choices of sPMF, are summarised in Fig. 8. A few important observations immediately emerge:

Fig. 8 Net pressure curves, as obtained with the sRPM(ϕ_n) (reference) as well as with various versions of ϕ_w. The bulk salt concentrations are within the regime 1.6–2 M. Graph (a) focuses on the large separation regime. The negative part of the ϕ_w(a, a) curve has been removed in the inset of graph (b), that displays the net pressures on a log scale.

• The long-ranged pressures that are obtained with models in which the additional sPMF is attractive between unlike charges, and repulsive between like charges, is large. That is, the net pressures obtained from sRPM(ϕ_w(a, r)) are qualitatively similar to those of sRPM(ϕ_n). The sRPM(ϕ_w(a, r)) model generates a somewhat stronger repulsion, but the difference is small.

• The qualitative nature of the sPMF is quite important. Specifically, the sRPM(ϕ_w(a, a)) model only produces a short-ranged repulsive regime, outside of which we even note a weak attractive interaction. The sRPM(ϕ_w(a, 0)) model is an intermediate case, where the repulsion clearly extends further than with the pure RPM, but it is not as long-ranged as with ϕ_w(a, r) or ϕ_n.

• We do not observe any force oscillations in the investigated separation regime (irrespective of the chosen sPMF). SFA measurements^11–13 have reported force “jumps”, indicative of such oscillations. On the other hand, such behaviours were not observed in a recent AFM²⁶ study. It should be noted that possible oscillations due to solvent packing will be unaccounted for by our implicit solvent treatment.

It is also of interest to compare the cluster tendency displayed by some of our investigated models. On Fig. 9 we plot cluster probability distributions, P_c(N_c), for a range of different models. Here, P_c measures the probability for a cluster of size N_c. A cluster is defined such that an ion must be within a distance δ_c or less of at least one other ion within a cluster, in order to be a member of that cluster. We have set δ_c = d + 2 Å. These distributions were obtained from grand canonical slit simulations, where the chemical potential is adjusted so that the corresponding bulk solution has a concentration of about 1.8 M. The distributions are based on at least 100 different configurations, separated by 10⁸ attempted moves to ensure statistical independence.

Fig. 9 Cluster probability distributions. (a) Comparing the cluster distribution (P_c(N_c)), at a rather wide surface separation (61 Å), of the sRPM(ϕ_n) with RPM:s with a uniform dielectric constant of 23 and 78.3. Also shown is the cluster distribution for the sRPM(ϕ_n) at a narrow separation, 21 Å. (b) The long-ranged tail of P_c(N_c), at a surface separation of 61 Å, for sRPM(ϕ_n), sRPM(ϕ_w(a, r)) and sRPM(ϕ_w(a, a)).

We note that the RPM approach does not generate large clusters, even if the (uniform) dielectric constant is set as low as 23. The sRPM:s, on the other hand, are able to produce quite large ion clusters. Interestingly enough, this is true also for the sRPM(ϕ_w(a, a)), even though we have seen that the corresponding surface forces are rather short-ranged. Nevertheless, the large cluster tail does not extend as far as with sRPM(ϕ_w(a, r)) or (in particular) sRPM(ϕ_n). Even so, these results imply that it is possibly not only the number of particles in each cluster that matters, but also the intrinsic structure of such clusters. This is supported by the long-ranged radial distribution tails that we observe with ϕ_n and ϕ_w(a, r), but not with ϕ_w(a, a) (see the ESI,† for details). In future work, we will scrutinise this further, making use of polymer classical density functional theory.

In graph (a) of Fig. 9 we also see that, as expected, very large clusters become improbable at short separations. It should be emphasised that this is not simply an effect of a diminishing number of ions within the simulation box in this regime. This is illustrated in the ESI,† where we find a very similar cluster probability distribution at short separations, but with an increased lateral (x, y) size of the simulation box.

Another interesting aspect is the average cluster charge, and how this varies with separation. We provide examples in Fig. 10. We note that the slopes of the curves, which provide the average net charge of an ion in a cluster, are similar between the models, but vary strongly with separation. At h = 21 Å, the slope is about +0.25e₀, dropping to about +0.075e₀ at h = 61 Å. The main reason is most likely a diminishing relative fraction of counterions, as the separation increases.

Fig. 10 The average net cluster charge, at two different separations, as obtained with the sRPM(ϕ_n), sRPM(ϕ_w(a, r)) and sRPM(ϕ_w(a, a)).

4 Conclusion

In this article we have demonstrated that a modified RPM with an added short-ranged PMF, adjusted to provide the model more realistic saturation properties, produces surface forces in much better agreement with experiments than can be obtained with the RPM model itself. This is related to the increased prevalence of clusters, brought about by the added short-ranged solvent-induced potential, that we, in the initial part of this work (concentration dependence), chose to model by a locally enhanced electrostatic coupling. Density–density correlations appear to dominate at high salt concentrations. Thus, we hypothesise that the experimentally observed anomalous screening lengths originate from cluster correlations rather than simple ion correlations. One apparent failure of the sRPM(ϕ_n) is its inability to predict a growth in the interaction screening length with electrolyte concentration. This is possibly due to the lack of many-body contributions to the ion–ion interaction, which potentially leads to an underestimation of the effects from concentration on particle interactions. A future aim is to implement a simple many-body scheme which will provide some estimate of this effect.

In the second part of this work, we established the qualitative nature of the sPMF is of crucial importance, even though its functional form is not. Our combined surface force and bulk solution cluster analyses suggest that the internal structure of the ion clusters may be as important as their average size. This will be investigated more closely in future work, employing classical polymer density functional theory.

Conflicts of interest

There are no conflicts to declare.

Data availability

All generated data, as well as all codes used to generate these data, are freely available upon reasonable request.

Acknowledgements

We thank professor Susan Perkin for sending us raw data from SFA measurements, as well as for fruitful discussions. Professors Sture Nordholm and Christian Holm are also acknowledged for enlightening discussions. J. F. acknowledges financial support by the Swedish Research Council, and computational resources by the Lund University computer cluster organisation, LUNARC.

References

1. J. N. Israelachvili, Intermolecular and Surface Forces, Academic Press, London, 2nd edn, 1991.
2. F. A. Evans and H. Wennerström, The colloidal domain: where Physics, Chemistry, Biology and Technology meet, VCH Publishers, New York, 1994.
3. C. Holm, P. Kekicheff and R. Podgornik, Electrostatic Effects in Soft Matter and Biophysics, Kluwer Academic Publishers, Dordrecht, 2001.
4. B. V. Derjaguin and L. Landau, Acta Physicochim. URSS, 1941, 14, 633–662.
5. E. J. W. Verwey and J. T. G. Overbeek, Theory of the Stability of Lyophobic Colloids, Elsevier Publishing Company Inc., Amsterdam, 1948.
6. P. Attard, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1993, 48, 3604–3621.
7. R. Leote de Carvalho and R. Evans, Mol. Phys., 1994, 83, 619–654.
8. J. Forsman, D. Ribar and C. E. Woodward, Phys. Chem. Chem. Phys., 2024, 26, 19921–19933.
9. J. G. Kirkwood, Chem. Rev., 1936, 19, 275–307.
10. G. M. Torrie and J. P. Valleau, J. Chem. Phys., 1980, 73, 5807–5816.
11. 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–9679.
12. M. A. Gebbie, H. A. Dobbs, M. Valtiner and J. N. Israelachvili, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 7432–7437.
13. A. M. Smith, A. A. Lee and S. Perkin, J. Phys. Chem. Lett., 2016, 7, 2157–2163.
14. G. R. Elliott, K. P. Gregory, H. Robertson, V. S. Craig, G. B. Webber, E. J. Wanless and A. J. Page, Chem. Phys. Lett., 2024, 843, 141190.
15. F. Coupette, A. A. Lee and A. Härtel, Phys. Rev. Lett., 2018, 121, 075501.
16. B. Rotenberg, O. Bernard and J.-P. Hansen, J. Phys.: Condens. Matter, 2018, 30, 054005.
17. R. Kjellander, J. Chem. Phys., 2018, 148, 193701.
18. S. W. Coles, C. Park, R. Nikam, M. Kanduc, J. Dzubiella and B. Rotenberg, J. Phys. Chem. B, 2020, 124, 1778–1786.
19. J. Zeman, S. Kondrat and C. Holm, Chem. Commun., 2020, 56, 15635–15638.
20. P. Cats, R. Evans, A. Härtel and R. van Roij, J. Chem. Phys., 2021, 154, 124504.
21. A. Härtel, M. Bültmann and F. Coupette, Phys. Rev. Lett., 2023, 130, 108202.
22. P. Gaddam and W. Ducker, Langmuir, 2019, 35, 5719–5727.
23. H. Yuan, W. Deng, X. Zhu, G. Liu and V. S. J. Craig, Langmuir, 2022, 38, 6164–6173.
24. R. J. E. Reinertsen, S. Kewalramani, F. Jiménez-Ángeles, S. J. Weigand, M. J. Bedzyk and M. Olvera de la Cruz, Proc. Natl. Acad. Sci. U. S. A., 2024, 121, e2316537121.
25. H. Robertson, G. R. Elliott, A. R. J. Nelson, A. P. Le Brun, G. B. Webber, S. W. Prescott, V. S. J. Craig, E. J. Wanless and J. D. Willott, Phys. Chem. Chem. Phys., 2023, 25, 24770–24782.
26. S. Kumar, P. Cats, M. B. Alotaibi, S. C. Ayirala, A. A. Yousef, R. van Roij, I. Siretanu and F. Mugele, J. Colloid Interface Sci., 2022, 622, 819–827.
27. K. Ma, J. Forsman and C. E. Woodward, J. Chem. Phys., 2015, 142, 174704.
28. J. B. Hasted, D. M. Ritson and C. H. Collie, J. Chem. Phys., 2004, 16, 1–21.
29. J. de Souza, A. A. Kornyshev and M. Z. Bazant, J. Chem. Phys., 2022, 156, 244705.
30. B. Conway and S. Marshall, Aust. J. Chem., 1983, 36, 2145–2161.
31. D. J. Bonthuis, S. Gekle and R. R. Netz, Langmuir, 2012, 28, 7679–7694.
32. I. Danielewicz-Ferchmin, E. Banachowicz and A. Ferchmin, J. Mol. Liq., 2013, 187, 157–164.
33. R. M. Adar, T. Markovich, A. Levy, H. Orland and D. Andelman, J. Chem. Phys., 2018, 149, 054504.
34. T. R. Underwood and I. C. Bourg, J. Phys. Chem. B, 2022, 126, 2688–2698.
35. D. Ribar, C. E. Woodward, S. Nordholm and J. Forsman, J. Phys. Chem. Lett., 2024, 15, 8326–8333.
36. S. Stenberg and J. Forsman, Langmuir, 2021, 37, 14360–14368.
37. A. R. Crothers, C. Li and C. Radke, Adv. Colloid Interface Sci., 2021, 288, 102335.
38. R. Evans, R. J. F. Leote de Carvalho, J. R. Henderson and D. C. Hoyle, J. Chem. Phys., 1994, 100, 591603.
39. L. Degreve and F. L. B. da Silva, J. Chem. Phys., 1999, 111, 5150–5156.
40. K. Komori and T. Terao, Chem. Phys. Lett., 2023, 825, 140627.
41. P. Auffinger, T. E. Cheatham and A. C. Vaiana, J. Chem. Theory Comput., 2007, 3, 1851–1859.
42. J. Tong, B. Peng, G. M. Kontogeorgis and X. Liang, J. Mol. Liq., 2023, 371, 121086.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5sm00027k

---