Strong coupling electrostatics for randomly charged surfaces: antifragility and effective interactions

Abstract

We study the effective interaction mediated by strongly coupled Coulomb fluids between dielectric surfaces carrying quenched, random monopolar charges with equal mean and variance, both when the Coulomb fluid consists only of mobile multivalent counterions and when it consists of an asymmetric ionic mixture containing multivalent and monovalent (salt) ions in equilibrium with an aqueous bulk reservoir. We analyze the consequences that follow from the interplay between surface charge disorder, dielectric and salt image effects, and the strong electrostatic coupling that results from multivalent counterions on the distribution of these ions and the effective interaction pressure they mediate between the surfaces. In a dielectrically homogeneous system, we show that the multivalent counterions are attracted towards the surfaces with a singular, disorder-induced potential that diverges logarithmically on approach to the surfaces, creating a singular but integrable counterion density profile that exhibits an algebraic divergence at the surfaces with an exponent that depends on the surface charge (disorder) variance. This effect drives the system towards a state of lower thermal ‘disorder’, one that can be described by a renormalized temperature, exhibiting thus a remarkable antifragility. In the presence of an interfacial dielectric discontinuity, the singular behavior of counterion density at the surfaces is removed but multivalent counterions are still accumulated much more strongly close to randomly charged surfaces as compared with uniformly charged ones. The interaction pressure acting on the surfaces displays in general a highly non-monotonic behavior as a function of the inter-surface separation with a prominent regime of attraction at small to intermediate separations. This attraction is caused directly by the combined effects from charge disorder and strong coupling electrostatics of multivalent counterions, which dominate the surface–surface repulsion due to the (equal) mean charges on the two surfaces and the osmotic pressure of monovalent ions residing between them. These effects can be quite significant even with a small degree of surface charge disorder relative to the mean surface charge. The strong coupling, disorder-induced attraction is typically much stronger than the van der Waals interaction between the surfaces, especially within a range of several nanometers for the inter-surface separation, where such effects are predicted to be most pronounced.

---

I. Introduction

The counter-intuitive electrostatic effects produced by mobile multivalent counterions in the vicinity of charged macromolecular surfaces (such as biopolymers like DNA, membranes, colloids, nano-particles and virus-like nano-capsids) have been studied extensively within the counterion-only models and in the case of non-disordered, homogeneous surface charge distributions (see, e.g., recent reviews in ref. 1–8 and references therein). These effects include, most notably, the so-called like-charge attraction induced by multivalent counterions between charged surfaces of the same sign. The like-charge attraction phenomena stand at odds with the weak coupling (WC) paradigm based traditionally on the mean-field Poisson–Boltzmann (PB) theory.^9–13 They require the framework of the strong coupling (SC) theories^{1–7,14–19} that incorporate strong ion–surface correlations, to the leading order, and ion–ion correlations, to subleading orders, ubiquitous in the situation where electrostatic interactions are strong; this is realized at low temperatures, solvents with low dielectric constant, and/or with highly charged macromolecular surfaces but, most prominently, with multivalent counterions. Usually, however, charged bio-soft systems with multivalent counterions include also a finite amount of monovalent salt ions that couple weakly to charged surfaces. A few experimental examples of such systems are furnished by the condensation of DNA by multivalent cations in the bulk^20–26 or in viruses and virus-like nano-capsids,^27–29 and formation of large aggregates (or bundles) of highly charged biopolymers such as DNA,^30,31 F-actin³² and microtubules.³³ In these situations, one deals with a difficult problem in which neither the WC nor the SC limiting laws can be applied without reservations. For such asymmetric Coulomb fluids, a generalized dressed multivalent-ion approach, that bridges the two limits in one single theoretical framework, has been introduced^7,34 and tested against extensive explicit- and implicit-ion simulations.^7,34–38

In counterion-only systems and in the case of uniformly charged surfaces with surface charge density −σe₀, the strength of the electrostatic interactions is described by the electrostatic coupling parameter Ξ = q²[small script l]_B/μ = 2πq³[small script l]_B²σ.¹⁴ It is the ratio of the rescaled Bjerrum length, q²l_B with [small script l]_B = e₀²/(4πεε₀k_BT), that measures the strength of electrostatic interactions between counterions of charge valency q, and the Gouy–Chapman length, μ = 1/(2πl_Bqσ), that measures the strength of electrostatic interaction between a counterion and the surface charge in units of thermal energy k_BT (conventionally, we take σ > 0 and q > 0). In the WC regime, Ξ ≪ 1, typically realized with monovalent counterions, thermal energy is dominant and, therefore, the counterions form a diffuse gas-like phase next to a charged surface that can be described by the mean-field PB theory. In the SC regime, Ξ ≫ 1, typically realized with multivalent counterions, electrostatic interactions are large enough to reduce the three-dimensional ionic cloud above the charged surface to a quasi-two-dimensional, strongly correlated liquid layer, or even a 2D crystal.^{1–7,14–16,19}

Subsequent studies also revealed that the WC–SC paradigm can be generalized to more complicated and realistic situations usually encountered in an experimental context, with an asymmetric Coulomb fluid comprising a mixture of mono- and multivalent ions in the vicinity of charged macromolecular interfaces. Neither the WC nor the SC paradigm can be consistently applied to this situation where the system is electrostatically part strongly and part weakly coupled. In this case the dressed multivalent-ion theory, based on piecewise application of the WC and SC formalism, has been shown to possess a robust regime of validity.^7,34–38 Other cases demanding an extension and/or modification of the WC–SC framework have been recently reviewed in ref. 7.

A separate issue, of the nature of the charge distribution on macromolecular surfaces, has also received progressively more focused attention. In many cases, these surfaces are not only inhomogeneous in terms of their charge configuration, but exhibit a fundamentally disordered distribution of charges^2,39–72 dictated either by the method of sample preparation and/or by the inherent structural properties of the surface materials.^39–48 These disordered systems, with either thermalized, annealed surface charges^{50,51,59,60,66–71} or structurally disordered, quenched surface charges^2,47–67 (or even partially annealed or partially quenched ones^58,72) have received a more rigorous attention, as it is becoming clear that the nature of charge disorder can interfere fundamentally with the behavior of the system approximated with fixed, homogeneous surface charges.

In the most realistic model of the interaction between charged macromolecules, one would thus need to consider not only an asymmetric confined Coulomb fluid between macromolecular surfaces, with a mixture of mono- and multivalent salt ions, but also surface charges that exhibit a disordered component. While off hand this seems to be like a formidable task, detailed calculations have shown that the final results are rather intuitive and revealing.^56–64 First of all, it became clear that, for quenched charge disorder, which is of interest in this paper, and in the absence of dielectric asymmetry between the solution subphase and the macromolecular surface, there is no disorder-induced effect in the WC regime.^56,58 However, in the SC limit and in the same general regime of parameters, one finds a long-ranged, disorder-induced attraction, which must be added to the repulsive force between equally charged surfaces.⁵⁶ It has also been shown that in the linear WC regime with Gaussian field-fluctuations around the mean-field solution, the quenched disorder would lead to repulsive or attractive interactions provided that the system has an inhomogeneous dielectric constant or monovalent salt distribution (in this case, repulsion arises when the medium is more polarizable than the surfaces and vice versa).^57,59–63

Just as the coupling parameter Ξ quantifies the electrostatic interactions between mobile counterions and fixed, homogeneous surface charges, the disorder coupling parameter, χ = 2πq²[small script l]_B²g, quantifies the strength of the disorder-induced effects in the case of disordered surfaces with g being the surface charge variance within the Gaussian Ansatz for its statistical distribution.^56–58 Obviously, the two coupling parameters, Ξ and χ, are very similar but the latter depends in a noticeably less drastic manner on the valency of the counterions. Just as the overall features of the behavior of a system with homogeneous charge distributions depend on the value of Ξ, the behavior of a disordered system depends in a very fundamental manner not only on the exact value of χ, but also on the presence of salt ions and dielectric inhomogeneities (“image charges”) in the system.^59–64

Here, we study the effective interaction mediated by strongly coupled Coulomb fluids between two plane-parallel dielectric surfaces carrying quenched, random monopolar charges with equal mean and variance, both when the Coulomb fluid consists only of mobile multivalent counterions and when it consists of an asymmetric ionic mixture. We analyze the consequences that follow from the interplay between surface charge disorder, dielectric and salt image effects, and the strong electrostatic coupling that result from multivalent counterions on the distribution of these ions and the effective interaction pressure they mediate between the surfaces. In a dielectrically homogeneous system, the multivalent counterions are found to be attracted towards the surfaces with a singular, disorder-induced potential that diverges logarithmically on approach to the surfaces, creating a singular but integrable counterion density profile that exhibits an algebraic divergence at the surfaces with an exponent given by the disorder coupling parameter. This effect drives the system towards a state of lower thermal ‘disorder’, one that can be described by a renormalized effective temperature, exhibiting thus a remarkable antifragility⁷³ reported in a previous work by the present authors.⁶⁴ We thus find that, in this general context of strongly coupled, disordered systems, anti-fragile behavior is ubiquitous and stems from the interplay of structural and thermal disorder, with the former engendering a decrease in the translational entropy of multivalent counterions as a consequence of the fact that the disordered charge distribution generates a finite degree of (non-thermal) configurational entropy. In the presence of an interfacial dielectric discontinuity, the singular behavior of the counterion density at the surfaces is removed but multivalent counterions are still accumulated much more strongly close to randomly charged surfaces as compared with uniformly charged ones. The interaction pressure between the surfaces displays in general a highly non-monotonic behavior as a function of the inter-surface separation with a prominent regime of attraction at small to intermediate separations. This attraction is caused directly by the combined effects from charge disorder and SC electrostatics of multivalent counterions, which dominate the surface–surface repulsion due to the (equal) mean charges on the two surfaces and the osmotic pressure of monovalent salt ions residing between them. These effects are quite significant even with a small degree of surface charge disorder (variance) relative to the mean surface charge.

The organization of the paper is as follows: we introduce our model in Section II and present our theoretical formalism in Section III. The distribution of multivalent counterions in the disordered counterion-only case is studied in Section IV A and the effects of salt ions and dielectric inhomogeneities are considered in Section IV B. The behavior of the effective inter-surface pressure is studied in Sections IV C–IV E, followed by the conclusion and discussion in Section V.

II. The model

We consider two plane-parallel dielectric slabs of infinite surface area S, finite thickness b and dielectric constant ε_p, placed perpendicular to the z axis with the inner bounding surfaces separated by a distance d (see Fig. 1). The outer bounding surfaces of the slabs are assumed to be neutral, while the inner ones carry a quenched, spatially uncorrelated, random charge distribution ρ(r), characterized by a Gaussian probability weight

Fig. 1 Schematic view of two infinite, plane-parallel dielectric slabs of thickness b and dielectric constant ε_p in a bathing ionic solution of dielectric constant ε_m, containing a monovalent and multivalent salt mixture (with multivalent counterions shown by large blue spheres and monovalent salt anions and cations shown by small orange and blue spheres). The multivalent counterions are confined in the slit between the randomly charged inner surfaces of the slabs.

Here, C is a normalization factor and ρ₀(r) and g(r) are the mean and the variance of the disordered charge distribution, respectively. For the specific model considered here, we assume that the charge distributions on the inner surfaces of the slabs are statistically identical and given by

ρ₀(r) = −σe₀[δ(z + d/2) + δ(z − d/2)], (2)

g(r) = ge₀²[δ(z + d/2) + δ(z − d/2)], (3)

where, with no loss of generality, we take g ≥ 0 and σ > 0 (in general, the surfaces may also have no net charge but only a finite disorder variance, which we shall consider in detail elsewhere⁷⁴). We develop the main formalism for an arbitrary shape of the charged boundaries, which we then apply to the specific example of the two planar slabs described above.

The slabs are assumed to be immersed in a solution of dielectric constant ε_m containing an asymmetric Coulomb fluid composed of a monovalent 1 : 1 salt of bulk concentration n₀ and a multivalent q : 1 salt of bulk concentration c₀ with q > 0 being the charge valency of multivalent counterions. The dielectric slabs are assumed to be impermeable to mobile ions. The multivalent counterions are assumed to be confined within the slit −d/2 ≤ z ≤ d/2. This particular constraint enables us to reproduce the standard counterion-only results⁵⁶ as a limiting case (practically, multivalent counterions can be prevented from entering the outer regions, |z| > d/2 + b, by enclosing these regions in semi-permeable membranes). This assumption has no impact on our results in a large part of the parameter space, e.g., for slab thicknesses larger than the Debye screening length, and especially for semi-infinite slabs that will be of main interest in this work.

III. The formalism

A. Dressed multivalent-ion theory

The problem of an asymmetric Coulomb fluid containing monovalent and multivalent ions next to charged macromolecular surfaces is a complicated one, mainly because different components of the Coulomb fluid couple differently to the surface charges: multivalent counterions tend to couple strongly, and therefore generate non-mean-field effects, whereas monovalent anions and cations couple weakly, and are thus expected to follow the standard mean-field paradigms as they would in the absence of multivalent ions.^9–13 Strong coupling effects due to multivalent counterions within the counterion-only models can be described well using the recent SC theories that have been studied thoroughly over the last several years (see, e.g., ref. 1–8 and references therein). In the case of an asymmetric Coulomb fluid, although it is not a priori clear if a single approximation can be introduced in order to describe the hybrid nature of the electrostatic couplings in the presence of charged boundaries, it turns out that a simple generalization of the counterion-only SC theory,^{5–7,14–18} dubbed dressed multivalent-ion theory, provides a very good approximation in a large part of the parameter space as discussed in detail in ref. 7 and 34–38. This approach is obtained by first tracing the partition function over the degrees of freedom associated with monovalent ions within a linearization approximation and then virial expanding it with respect to the fugacity of multivalent counterions. Both these steps can be justified systematically³⁴ only in the case of highly asymmetric mixtures with a sufficiently large counterion valency q and provided that the bulk concentration of multivalent counterions is relatively small (as is often the case in experimental systems containing asymmetric ionic mixtures^20–33), because otherwise the different ionic species in the solution cannot be treated on such different levels. The dressed multivalent-ion theory described and reviewed recently⁷ bridges the gap between the weak coupling DH theory and the counterion-only SC theory and reproduces them as two special limiting cases at large and small screening parameters, respectively.

B. General form of the interaction free energy

In the dressed multivalent-ion theory, the grand-canonical partition function of an asymmetric Coulomb fluid with a fixed realization of external charges ρ(r) can be written, in a functional-integral form (see, e.g.ref. 75–79 and references therein), as follows:³⁶

Here, β = 1/(k_BT), ϕ(r) is the fluctuating electrostatic potential and S[ϕ, ρ] is the effective “field-action”, i.e.,

where λ_c = c₀ is the fugacity (or bulk concentration) of multivalent counterions, Ω_c(r) is an indicator function that determines the region of space that is accessible to multivalent counterions (e.g., in the present model, we have Ω_c(r) = (θ(z + d/2) − θ(z − d/2)), and G⁻¹(r, r′) is the operator inverse of Green's function G(r, r′) that, in the dressed multivalent-ion theory, satisfies the DH equation

−ε₀∇·ε(r)∇G(r, r′) + ε₀ε(r)κ²(r)G(r, r′) = δ(r − r′). (6)

Here, κ(r) is the Debye (or salt) screening parameter which is non-zero only outside the region occupied by the dielectric slabs, i.e., κ² = 4πl_Bn_b, where n_b = 2n₀ + qc₀ is the bulk concentration due to all monovalent ions. Therefore, the role of monovalent ions is incorporated into the Green's function on the DH level.

As noted before, in the presence of multivalent counterions with sufficiently small bulk concentration, the partition function can be virial-expanded in terms of the counterion fugacity, λ_c, as

[scr Z, script letter Z][ρ] = [scr Z, script letter Z]₀[ρ] + λ_c[scr Z, script letter Z]₁[ρ] + [scr O, script letter O](λ_c²), (7)

with the first two terms corresponding to the dressed multivalent-ion theory on the leading order. The zeroth-order term corresponds to the electrostatic interaction of fixed charged objects in the absence of multivalent counterions,and the first-order term gives the single-particle contribution due to multivalent counterions

[scr Z, script letter Z]₁[ρ] = [scr Z, script letter Z]₀[ρ]∫drΩ(r)e^{−βu(r;[ρ])}.

Here, u(r;[ρ]) is the single-particle interaction energy that has the form

where G_im(r, r) is the generalized Born energy contribution, which is generated purely by the dielectric and/or salt polarization effects (or the “image charges”), i.e., G_im(r, r) = G(r, r) − G₀(r, r), with G₀(r, r) representing the formation (self-)energy of individual counterions in a homogeneous background; this can be calculated from the screened free-space Green's function defined through −ε₀ε_m(∇² − κ²)G₀(r, r′) = δ(r − r′).

In the case of counterion-only systems, the two leading terms in eqn (7) were shown to generate a finite contribution to the free energy, constituting the SC theory.^{5–7,14–18} The higher-order terms contain subleading contributions from multi-particle interactions that become important only in the crossover regime between the WC and SC regimes. In the more general context of the dressed multivalent-ion theory, the leading-order virial terms can generate both the SC and WC behaviors of the system in the limits of small and large screening parameters, respectively.³⁴ The predictions of this theory were analyzed for uniformly charged surfaces (and also for strictly neutral surfaces) and were compared with extensive numerical simulations elsewhere.^7,34–38

In the present study, the charge distribution of the dielectric surfaces, ρ(r), has a quenched, Gaussian disordered component with a finite variance around the mean charge density ρ₀(r). One therefore needs to average the thermodynamic quantities over different realizations of the charge disorder as well. Hence, for instance, the free energy follows from

β[scr F, script letter F] = −〈〈ln [scr Z, script letter Z][ρ]〉〉. (10)

The disorder average is given by 〈〈…〉〉 = ∫[scr D, script letter D]ρ(…)[scr P, script letter P][ρ] and the Gaussian weight by eqn (1). The averaged quantity 〈〈ln [scr Z, script letter Z][ρ]〉〉 can be calculated in general by employing the Edwards–Anderson's replica Ansatz.⁵⁶ Since we are interested only in the leading virial terms, we can also directly average the leading-order free energy terms over the quenched charge disorder weight, yielding

The averages in the above expression can be carried out straightforwardly and we find the (grand-canonical) dressed multivalent-ion free energy of the system as

where the first and the second terms are the contributions in the absence of multivalent counterions, representing the interaction free energies due to the mean surface charge density, ρ₀(r), and the disorder variance, g(r), respectively. The second term arises also in the analysis of the fluctuation-induced forces between disordered surfaces in vacuum or in a weakly coupled Coulomb fluid;^57,59–63 it gives a non-vanishing contribution only in inhomogeneous systems with a finite dielectric discontinuity at the bounding surfaces and/or a spatially inhomogeneous distribution of salt ions. The third term represents the contribution from multivalent counterions on the leading (single-particle) level, in which u(r) is the effective single-particle interaction energy⁶⁴

We note that the first term in eqn (13) originates from the interaction of multivalent counterions with the mean surface charge density, the second term from the self-interactions of individual counterions (with their own image charges) and the third term is due to the presence of surface charge disorder. This term is proportional to the disorder variance and shows an explicit temperature dependence and a quadratic dependence on the Green's function and the multivalent ion charge valency, q; it arises from the sample-to-sample fluctuations (or variance) of the single-particle interaction energy (9).⁶⁴

In the SC limit or within the multivalent dressed-ion theory,^14,16,34 the density profile of multivalent counterions can be derived solely in terms of the effective single-particle interaction energy as

c(r) = λ_cΩ(r)e^−βu(r). (14)

This concludes the recapitulation of the multivalent dressed-ion theory that forms the basic framework of our analysis of electrostatic coupling and quenched charge disorder in what follows.

C. Specific case of planar dielectric slabs

The free energy expression (12) is valid regardless of the shape of the boundaries. In the rest of this paper, however, we shall delimit ourselves to the specific example of two planar dielectric slabs (Section 2), in which case the Green's function G(r, r′) = G(ρ, ρ′;z, z′), with ρ = (x, y) and ρ′ = (x′, y′) being the transverse (in-plane) coordinates, is only a function of |ρ − ρ′|. Thus, its Fourier-Bessel transform Ĝ(Q; z, z′) can be defined by

Using standard methods, we find

where γ² = κ² + Q²,and

For semi-infinite slabs (b → ∞), one can recover the well-known expression

Note that the information about the image charge effects, which result from the inhomogeneous distribution of the dielectric constant or the bathing salt solution in the system, enters here through the parameter Δ_s. In the absence of salt ions (κ = 0), Δ_s reduces to the bare dielectric discontinuity parameter

In the cases treated below, relevant for aqueous solvents in the presence of bounding slabs of low dielectric constant, we delimit ourselves to Δ ≥ 0, which gives repulsive image interactions.

The effective single-particle interaction energy in the two-slab system can be written using eqn (13) as

The density profile of counterions in the slit region, −d/2 ≤ z ≤ d/2, follows as (see eqn (14))

c(z) = λ_ce^−βu(z), (22)

and the grand-canonical free energy (per k_BT and per unit area, S) can be written as

In the above expression, the different contributing terms appear in the same order as in eqn (12) and we have

where we have omitted irrelevant additive terms that are independent of the surface separation, d. The interaction free energy follows from eqn (23) by subtracting all such terms or, in other words, the reference free energy of the system for d → ∞.

D. Dimensionless representation

It is convenient to make use of a dimensionless set of quantities by rescaling the spatial coordinates with the Gouy–Chapman length μ = 1/(2πql_Bσ) as r̃ = r/μ. Other parameters are rescaled accordingly, e.g., the inter-surface separation, the salt screening parameter and an analogously defined length scale, χ_c² = 8πq²[small script l]_Bc₀, which is referred to as the rescaled bulk concentration of multivalent counterions, i.e.,

d̃ = d/μ, [small kappa, Greek, tilde] = κμ, [small chi, Greek, tilde]_c = χ_cμ. (25)

The rescaling of the Bjerrum length leads to the dimensionless electrostatic coupling parameter Ξ = q²[small script l]_B/μ = 2πq³[small script l]_B²σ, associated with the mean surface charge density,¹⁴ while the surface charge variance appears in the dimensionless disorder coupling (or disorder strength) parameter, χ = 2πq²[small script l]_B²g.⁵⁶ The rescaling of the density profile and the interaction pressure will be discussed later.

IV. Results

A. Distribution of counterions: the counterion-only model

Let us first consider the distribution of multivalent counterions in the special case of the counterion-only model, where the salt screening and image charge effects are set to zero by assuming a homogeneous system with κ = 0 and Δ = 0, or equivalently, ε_p = ε_m (this case was considered in a previous work,⁵⁶ which however focused only on the effective interaction between the surfaces and did not investigate the distribution of counterions). The two impermeable, randomly charged surfaces are placed at a separation distance, d, confining in the slit between them a fixed number of multivalent counterions, N, which is given by the mean charge on the two surfaces through the global electroneutrality condition, Nq = 2Sσ. Hence, the fugacity of counterions is given by^14,56

The Green's function of the counterion-only model reduces to the bare Coulomb interaction, G₀(r, r′) = 1/(4πε₀ε_m|r − r′|), and the effective single-particle interaction energy (up to an irrelevant additive constant and in rescaled units) follows straightforwardly from eqn (13) or (21) as

Note that while, because of the symmetric plane-parallel geometry of the model, there is no net attraction acting on individual counterions from the (oppositely signed) mean surface charges, the counterions experience an attractive potential due to the quenched surface charge disorder, which is given by the second term in eqn (27).

The rescaled density profile of counterions in the slit is obtained by using eqn (22) as

where we have defined

In the no-disordered case (χ = 0), the above expression reduces to the standard SC density profile c̃(z̃) = 2/d̃.^{3,5–7,14–16} As shown in Fig. 2, the density profile of multivalent counterions is strongly modified by the surface charge disorder and exhibits algebraic singularities close to the randomly charged boundaries at z = ±d/2 for any finite value of the disorder coupling parameter, in clear violation of the contact-value theorem established for uniformly charged surfaces.^80–83 This kind of behavior was discussed in our recent work on the distribution of multivalent counterions next to a single, randomly charged interface⁶⁴ and we shall not delve further into the details, making only a few remarks in what follows.

Fig. 2 Rescaled density profiles of multivalent counterions within the counterion-only model as a function of the rescaled normal position, z̃ = z/μ, between two randomly charged surfaces located at inter-surface separation d̃ = d/μ = 4. Different curves correspond to different values of the disorder coupling parameter as shown on the graph.

The singular behavior of the counterion density profile comes directly from the logarithmic disorder term in the single-particle energy in eqn (27), which drives the counterions towards the surfaces. The density of counterions away from the surfaces decreases as the disorder coupling parameter increases (Fig. 2) and, eventually, it tends to zero c̃(z̃) → 0 for any |z̃| < d̃/2 when χ → 2⁻ since C₀(χ) → ∞. Counterions are thus densely accumulated in the immediate vicinity of the surfaces. This behavior must be distinguished from the surface adsorption or counterion condensation phenomena⁷ as, in the present context, one can show systematically that the mean surface charge is not renormalized by the surface accumulation of counterions (see Appendix C in ref. 58).

The accumulation of counterions in the vicinity of the charged boundaries furthermore drives the system towards a state of lower thermal ‘disorder’, since, as one can show, the translational entropy of multivalent counterions decreases as a consequence of the disordered charge distribution, which generates a finite (non-thermal) configurational entropy; this latter type of entropy stems from the different realizations of the quenched disorder. This behavior was associated in ref. 64 with the anti-fragility⁷³ of the system, in which introducing an external (quenched) disorder source diminishes the intrinsic thermal disorder and drives the system towards a more ‘ordered’ state, which is also a thermodynamically more stable one as compared with the corresponding non-disordered case (see ref. 56 for the canonical free energy expression of the disordered counterion-only model).

The rescaled interaction pressure acting on each of the surfaces, P̃ = βP/(2π[small script l]_Bσ²), can be obtained as⁵⁶

which, by setting χ = 0, reduces to the standard SC pressure between two non-disordered like-charged surfaces¹⁴

The attractive (electrostatic) and repulsive (entropic) terms in the pressure give rise to an equilibrium bound-state separation of d̃_* = 2 between two identical, uniformly charged surfaces.^5,6,14–16 As is clear from eqn (30), the surface charge disorder gives an additive attractive contribution to the total interaction pressure, which renormalizes the entropic term and leads to a more closely packed bound state with the equilibrium inter-surface separation d̃_* = 2(1 − χ). This separation tends to zero for χ → 1⁻, beyond which the two surfaces collapse into a primary minimum. This is because the entropic contribution (second term in eqn (30)) changes sign at χ = 1 and, hence, the pressure becomes attractive at all separations (and diverges for d̃ → 0) in the regime χ > 1.

Another notable point is that the renormalization of the repulsive entropic term by a factor 1 − χ can be interpreted also as a renormalization of the effective temperature of the system to a lower value. However, this interpretation should be considered with caution because first, this particular form of the renormalization of the repulsive pressure is found only in the SC limit of the two-surface counterion-only model and, secondly, one can show, by inspecting the thermodynamic quantities of the system, that the renormalizing term, −2χ/d̃, in eqn (30) indeed comes from both electrostatic energy and entropic contributions.

B. Distribution of counterions: salt screening and image charge effects

In the next step, we assume that, in addition to the multivalent counterions that are introduced through an asymmetric q : 1 salt of bulk concentration c₀, the ionic mixture also contains a monovalent salt of bulk concentration n₀. As noted before, the Debye screening parameter is κ = (4πl_Bn_b)^1/2 with n_b = 2n₀ + qc₀. For the time being, we also assume that the system is dielectrically homogeneous, i.e., ε_p = ε_m and that the slabs are semi-infinite (b = ∞) and impermeable to all ions. The density profile of multivalent counterions can be calculated from eqn (21) and (22) and by making use of the appropriate expressions from eqn (16)–(19).

In the absence of surface charge disorder, the density profile of multivalent counterions shows a non-monotonic behavior with a peak at a small distance from each of the two bounding surfaces (see the black solid curve in Fig. 3a; note also that here the density profiles are rescaled with the bulk value c₀). This behavior is due to the interplay between two distinct factors, namely, the salt screening effect and the “salt image” effect. The former dominates at intermediate to large distances from the surfaces that are comparable to or larger than the Debye screening length [small kappa, Greek, tilde]⁻¹, while the latter dominates at small to intermediate distances from the surfaces, causing a partial depletion of multivalent counterions from the proximity of the surface boundaries. The salt image effect is generated because of the inhomogeneous distribution of salt ions that are not allowed to permeate into the wall regions (i.e., |z̃| > d̃/2), leading in turn to a discontinuous change in the polarization of the medium at the interfacial boundaries, with the slit region having a larger polarizability response than the slab region. This gives rise to repulsive “salt image” effects that show some qualitative similarities to the “dielectric image” effects.

Fig. 3 Rescaled density profiles of dressed multivalent counterions as a function of the rescaled normal position, z̃, between the randomly charged inner surfaces of two semi-infinite slabs placed at rescaled separation d̃ = 10 in (a) dielectrically homogeneous (Δ = 0) and (b) dielectrically inhomogeneous systems (Δ = 0.95) at fixed parameter values [small kappa, Greek, tilde] = 0.3 and Ξ = 50 and various disorder coupling parameters as shown on the graphs. (c) is the same as (b) but with fixed χ = 2 and varying slab thickness b̃ = 5, 10, 15 and ∞.

When the bounding surfaces carry a finite degree of quenched charge randomness (χ > 0), the multivalent counterions exhibit a strong attraction towards them and, again, generate a singular density profile with a diverging contact value (dashed curves in Fig. 3a). In fact, the counterion density profile exhibits the same algebraic singularity on approach to the surfaces as in the counterion-only case, eqn (28).⁶⁴ Thus, at small distances from the surfaces, the disorder-induced effects overcome both salt image and salt screening effects.

These features change qualitatively when the system is dielectrically inhomogeneous, exhibiting a finite dielectric discontinuity at the interfacial boundaries at z = ±d/2. In Fig. 3b, we show the counterion density profiles for Δ = 0.95 (corresponding to the dielectric discontinuity at the water/hydrocarbon boundary with ε_m = 80 and ε_p = 2) and a few different values of the disorder coupling parameter. Strong dielectric image repulsions lead to a zone of complete depletion near the surfaces with vanishing counterion density, followed by enhanced peaks at an intermediate distance from each surface. These features are qualitatively similar in the presence or absence of surface charge disorder, although the disorder generates larger densities, especially at the peak regions, by attracting a larger number of multivalent counterions from the bulk reservoir into the slit. The interfacial depletion zone is generated because (unlike the salt image effects) the dielectric image effects can be described in terms of point-like image charges (especially for large Δ ≃ 1);³⁶ these images repel the counterions with a singular image potential (second term in eqn. (13) and (21)) that behaves approximately as the first-image interaction potential βu_im ≃ ΞΔ/4(d̃/2 ± z̃) at small distances from the boundaries, thus overcoming the logarithmic, disorder-induced attraction experienced by the counterions.

We should also note that the finiteness of the slab thickness has typically only a small effect on the counterion distribution, especially when it is comparable to or larger than the screening length, κb ≳ 1, which is in fact often the case in realistic systems. For instance, in Fig. 3c, we show the results for the same parameters as in Fig. 3b but with χ = 2 and different slab thicknesses in the range b̃ ≥ 5 (which covers the range b > 1 nm in actual units, see the Discussion). The density of counterions in the slit is slightly increased but saturates quickly when the slab thickness is increased to infinity. The smaller counterion density found in the case of thinner slabs is due to the fact that the overall surface attraction experienced by counterions becomes smaller for smaller slab thicknesses and, at the same time, the salt ions in the outer region behind the slabs (see Fig. 1) also contribute more strongly to the screening effects.

An interesting effect seen in the above results is that the competition between disorder-induced attraction and image repulsion leads to a highly pronounced bimodal profile with two humps that correspond to two distinct counterion-populated regions in the slit. (Such bimodal profiles have also been found for the counterion density between heterogeneous but regularly patterned, planar charged surfaces⁸⁴ and also for the monomer density of polyelectrolyte chains between uniformly charged planar surfaces.⁸⁵) These humps are expected to appear at sufficiently large disorder variances or relatively large inter-surface separations as compared with the screening length. This kind of morphological change in the distribution of multivalent counterions can be quantified by defining the distance between the peaks, δz_p, as an analog of the ‘order parameter’ in the phase transition context. As seen in Fig. 4a (for [small kappa, Greek, tilde] = 0.3, d̃ = 8 and Δ = 0.95), this quantity shows a sharp, continuous change at a threshold value of χ_* ≃ 0.385 from a single-hump profile to a double-hump one. In Fig. 4b and c, we show the results obtained for the threshold value χ_* as a function of the dielectric discontinuity parameter, Δ (at fixed d̃ = 5), and as a function of the rescaled inter-surface distance, d̃ (at fixed Δ = 0.95), respectively. The region below (above) the curves in these figures corresponds to the parameter values for which the density profile is uni- (bi-)modal. As seen in Fig. 4b, at larger values of the dielectric discontinuity or at larger salt screening parameters, a larger disorder coupling parameter (disorder variance) is required in order to counteract the counterion-image repulsions from the boundaries and create a bimodal structure. The same is true for a system with a smaller inter-surface separation, see Fig. 4c.

Fig. 4 (a) Rescaled distance between the peaks (humps) in the density profile of multivalent counterions as a function of the disorder coupling parameter, χ, for [small kappa, Greek, tilde] = 0.3, d̃ = 8 and Δ = 0.95. (b) The threshold value χ_*, where the counterion density profile changes from a unimodal to a bimodal shape, as a function of the dielectric discontinuity parameter, Δ, at fixed d̃ = 5 and different values of [small kappa, Greek, tilde] as shown on the graph. (c) The threshold value χ_* as a function of the rescaled inter-surface distance, d̃, at fixed Δ = 0.95 and different values of [small kappa, Greek, tilde]. In (b) and (c), the region below (above) the curves correspond to the parameter values for which the density profile is uni- (bi-) modal. In all these cases, the slabs are semi-infinite.

The salt screening has a reverse effect at small or large inter-surface separations: whereas at small separations it increases the values of χ_*, at large separations it decreases them. Another point to be noted here is that, for the parameter values used in Fig. 4b and c, there is a plateau-like region with χ_* = 0, where the density profile is bimodal for any value of χ.

C. Interaction pressure

In the most general case, where the system is immersed in an asymmetric electrolyte bath, the interaction pressure (equivalent to the osmotic or disjoining pressure) acting on the slabs follows from the difference in the slit pressure and the bulk (electrolyte) pressure, i.e., P = P_S − P_b. The slit pressure is obtained by differentiating the free energy expression (23) with respect to the inter-surface separation as P_S = −∂[scr F, script letter F]/(S∂d), where all other parameters are kept fixed, and the bulk pressure is given by P_b = (n_b + c₀)k_BT with n_b = 2n₀ + qc₀ being the bulk concentration of the monovalent ions as defined before. It should be noted that the slit pressure obtained from the differentiation of the dressed multivalent-ion free energy with respect to d does not contain the contribution from the osmotic pressure of monovalent ions in the slit, P_mon, which can be calculated in terms of the mid-plane density of monovalent ions.³⁶

We decompose the rescaled interaction pressure, P̃ = βP/(2π[small script l]_Bσ²), into its different components as

P̃ = P̃_σ + P̃_dis + P̃_c + P̃_mon, (32)

where

The first term in eqn (32), P̃_σ, follows from the first term in eqn (23) and represents the repulsive pressure due to the mean charges on the inner surfaces of the slabs. The second term, P̃_dis, gives the contribution from the surface charge disorder (second term in eqn (23)) with f̃([small kappa, Greek, tilde], d̃, Δ) being defined using eqn (24) as

where Q̃ = Qμ and [small gamma, Greek, tilde] = γμ. This contribution can be attractive or repulsive and will be non-vanishing only in inhomogeneous systems with a finite dielectric discontinuity and/or an inhomogeneous distribution of salt ions (i.e., when Δ_s ≠ 0). These two contributions to the interaction pressure will be present regardless of the multivalent counterions; they have been analyzed in detail for semi-infinite slabs in ref. 57 and 59–63.

The contributions P̃_c and P̃_mon, on the other hand, represent the osmotic pressure components from multivalent counterions and monovalent salt ions, respectively. P̃_c is given in terms of the effective single-particle interaction energy, ũ(z̃), and follows from the third term in eqn (23) with ũ(z̃) obtained straightforwardly by rescaling the parameters in eqn (21).

Finally, P̃_mon is calculated from the contact-value theorem in terms of the total mid-plane density of monovalent ions, which can be estimated here through the relationship n(z) = λ₊exp[−βu(z)]|_q=1 + λ₋exp[−βu(z)]|_q=−1 as discussed in detail in ref. 36, where λ₊ = n₀ and λ₋ = n₀ + qc₀ are the bulk concentrations of monovalent cations and anions, respectively; the rescaled mid-plane density in eqn (36) is then defined as ñ(0) = n(z̃ = 0)/n_b.

D. Interaction of non-disordered surfaces

We first consider the interaction pressure between two non-disordered (uniformly charged) surfaces within the dressed multivalent-ion theory. The interaction pressure acting on the slabs in this case follows from eqn (32)–(36) by noting that in all the expressions involved we need to set χ = 0; hence, in particular we have P̃_dis = 0. We shall primarily focus on the case of semi-infinite slabs (b = ∞) and fix the electrostatic coupling parameter at Ξ = 50, which can be achieved with tetravalent counterions (q = 4) and the mean surface charge density σ = 0.24 nm⁻² in water (ε_m = 80) and at room temperature T = 293 K (see Table 1). Unless otherwise specified, we take Δ = 0.95, which is appropriate for water/hydrocarbon interfaces.

Table 1 A few typical examples for the actual values of the bulk concentrations n₀ and c₀, which can correspond to the typical values we have chosen for the rescaled parameters [small kappa, Greek, tilde] and [small chi, Greek, tilde]_c. Here, we have fixed the other parameter values as q = 4 (tetravalent counterions), l_B = 0.71 nm (corresponding to water with ε_m = 80 at room temperature, T = 293 K) and surface charge density σ = 0.24 nm⁻², which give μ = 0.23 nm and Ξ = 50. We also show the actual values of the disorder variance g, which can correspond to a few typical values of the disorder coupling parameter χ (see the text for definitions)

[small kappa, Greek, tilde] = 0.2	0.25	0.3	0.4			χ	g (nm^−2)
n 0 = 68 mM	108 mM	156 mM	279 mM	c 0 = 1.1 mM	[small chi, Greek, tilde] c = 0.1	0.5	0.01
65 mM	105 mM	153 mM	277 mM	2.5 mM	0.15	1.0	0.02
62 mM	101 mM	150 mM	273 mM	4.4 mM	0.20	2.0	0.04
57 mM	96 mM	145 mM	268 mM	6.9 mM	0.25	4.0	0.08

---

The results for the rescaled interaction pressure are shown in Fig. 5a as a function of the rescaled inter-surface separation for a few different values of the rescaled screening parameter, i.e., [small kappa, Greek, tilde] = 0.25, 0.27, 0.28, 0.29, 0.30, and [small chi, Greek, tilde]_c = 0.15 (in actual units, these parameter values can be obtained, for instance, by taking salt bulk concentrations n₀ ≃ 105, 123, 133, 143, 153 mM and c₀ ≃ 2.5 mM). At sufficiently large [small kappa, Greek, tilde] (e.g., high monovalent salt concentration), the pressure is repulsive and decays monotonically with the distance. In this regime, the salt screening effects are dominant and the SC effects due to the multivalent counterions are strongly suppressed; the repulsive pressure is in fact determined mainly by the surface–surface repulsion (P̃_σ) and the osmotic pressure of monovalent ions (P̃_mon) with the latter contribution being the larger of the two, as can be seen from the pressure components in Fig. 5b.

Fig. 5 (a) Rescaled interaction pressure as a function of the rescaled distance between the uniformly charged inner surfaces of two identical semi-infinite dielectric slabs with Ξ = 50, Δ = 0.95 and [small chi, Greek, tilde]_c = 0.15 and for a few different values of the rescaled screening parameter, [small kappa, Greek, tilde], as indicated on the graph. Panels (b) and (c) show the different components that contribute to the interaction pressure for [small kappa, Greek, tilde] = 0.3 and [small kappa, Greek, tilde] = 0.25, respectively (see the text for definitions).

As [small kappa, Greek, tilde] is decreased, the pressure develops a non-monotonic behavior with a pronounced local minimum at intermediate inter-surface separations within a range comparable to the Debye screening length (Fig. 5a). This local minimum eventually turns into a negative global minimum as [small kappa, Greek, tilde] is decreased further, indicating a strong attractive pressure induced between the slabs. For [small kappa, Greek, tilde] = 0.27 (red dashed curve), the minimum attractive pressure is approximately P̃ ≃ −2.7 (or, with the choice of the actual parameter values mentioned above, P ≃ −28 bar). For [small kappa, Greek, tilde] = 0.25 (black solid curve in Fig. 5a which is replotted as the red solid curve in Fig. 5c), the minimum attractive pressure is P̃ ≃ −75 (or, P ≃ −780 bar). Note that the net attraction here appears despite the fact that the surfaces are like charged. This is because of the SC effects, which are produced by the leading surface–counterion correlations^{5–7,14–16} and enter through the pressure component P̃_c. This contribution becomes quite large at small screening parameters (see Fig. 5c) and exhibits a non-monotonic behavior that we shall consider later in more detail. The non-monotonic behavior of the net pressure with the inter-surface separation (Fig. 5a) stems directly from the interplay between its different components, with P̃_c being the most essential one.

Similar behavior to those shown in Fig. 5a can be seen upon increasing the rescaled bulk concentration of multivalent counterions, [small chi, Greek, tilde]_c (Fig. 6), and upon decreasing the dielectric discontinuity parameter, Δ (not shown) with the minimum attractive pressure turning out to be quite sensitive to the exact values of these parameters. Fig. 6a shows the results for [small chi, Greek, tilde]_c = 0.09, 0.11, 0.13, 0.15 and [small kappa, Greek, tilde] = 0.27 (corresponding, for instance, to choosing the bulk concentrations as c₀ ≃ 0.9, 1.3, 1.9, 2.5 mM with n₀ ≃ 127, 126, 125, 123 mM, respectively). The pressure still remains strongly repulsive at very small separations because of the contributions P̃_σ and P̃_mon as discussed above (see Fig. 6b and also Fig. 5b and c). However, the multivalent counterions pressure component, P̃_c, again shows a non-monotonic behavior with both attractive and repulsive regions (Fig. 6b).

Fig. 6 (a) Rescaled interaction pressure as a function of the rescaled distance between the uniformly charged inner surfaces of two identical semi-infinite dielectric slabs with Ξ = 50, Δ = 0.95 and [small kappa, Greek, tilde] = 0.27 and for a few different values of the rescaled multivalent counterions concentration, [small chi, Greek, tilde]_c, as indicated on the graph. The inset shows a close-up view of the region around the shallow local minimum at large separations. Panels (b) and (c) show different components that contribute to the interaction pressure for [small chi, Greek, tilde]_c = 0.11 in the vicinity of the small-separation and large-separation minima, respectively.

The non-monotonic behavior of P̃_c can be understood by inspecting the average number of multivalent counterions in the slit between the slabs, , which, in rescaled units and per unit area, is given by

This quantity is shown in Fig. 7 for [small kappa, Greek, tilde] = 0.27 and [small chi, Greek, tilde]_c = 0.15 along with its corresponding pressure component P̃_c that in fact corresponds to the red dashed curve in Fig. 5a (note that P̃_c is directly related to ∂N̄/∂d̃ as one can see by comparing eqn (22), (35) and (38)). At large separations and upon decreasing the inter-surface distance, the number of multivalent counterions in the slit is increased due to a larger uptake of these ions from the bulk solution, which is caused by an increased counterion–surface correlation (attraction) that, in turn, enhances the attractive (negative) inter-surface pressure component P̃_c. As the separation is decreased further, the number of multivalent counterions in the slit reaches a maximum value. The multivalent counterions at smaller inter-surface separations are strongly repelled by their image charges and their number in the slit decreases; they are eventually completely ejected from the slit when the surfaces come close to contact. The pressure component due to multivalent counterions thus changes sign at the location where N̄ reaches a maximum and, eventually, tends to the bulk pressure P̃_c → −[small chi, Greek, tilde]_c²/4 (or, in actual units, P_c → c₀k_BT) when d̃ → 0 (this limiting value is not discernible at the range of scales shown in Fig. 7).

Fig. 7 Rescaled average number of multivalent counterions in the slit per rescaled surface area (blue curve) and the corresponding rescaled pressure component, P̃_c (red dashed curve), as a function of the rescaled distance between the uniformly charged inner surfaces of two identical semi-infinite dielectric slabs for [small kappa, Greek, tilde] = 0.27, [small chi, Greek, tilde]_c = 0.15, Ξ = 50 and Δ = 0.95.

We should also note that the intermediate attractive-pressure regime seen in Fig. 5a and 6a is followed by a weakly repulsive regime (hump) at larger separations but there is a very shallow local minimum at large separations that gives the pressure curves a weakly attractive long tail. This is shown in a close-up view in the inset of Fig. 6a for the curves that appear in the main set. Comparing the different pressure components around this large-separation minimum (Fig. 6c) shows that the contribution from the surface–surface repulsion (P̃_σ) nearly cancels the attractive contribution from multivalent counterions (P̃_c) and, thus, the total pressure is determined almost completely by the contribution from monovalent ions (P̃_mon), which itself shows a shallow minimum at large separations.

In general, the depth of this minimum, and the large-separation attraction regime, can be enhanced slightly also by P̃_c (not shown). A closer inspection of P̃_mon shows that the partial osmotic pressure that contributes to this component from monovalent salt anions is negative (because these ions have the same sign as the surface charges and are therefore depleted from the slit) and is slightly larger in magnitude than the (positive) partial osmotic pressure from monovalent salt cations, hence, causing the weakly attractive, large-separation behavior of P̃_mon. This behavior is robust and appears for a wide range of parameter values.

Another important point to be considered here is whether the van-der-Waals-like⁸⁶ non-monotonicity observed in the pressure–distance curves in Fig. 5a and 6a gives any indication of a phase coexistence, if the pressure curves are to be interpreted in the thermodynamic sense as, e.g., the pressure in a stack of plane-parallel like-charged membranes? Such a phase coexistence would imply a simultaneous existence of a dense phase and a swollen phase in the system. However, if one applies the Maxwell equal-area construction, the negative area under the pressure curves, which comes from the region around the attractive minimum at small separations, is typically much larger than the positive area coming from the region around the repulsive maximum (hump) at intermediate separations. This translates into a statement that the equal-area Maxwell construction cannot be fulfilled in general.

E. Interaction of disordered surfaces

We now consider the situation where the inner surfaces of the slabs carry a finite degree of charge disorder characterized by the dimensionless coupling parameter, χ. The rescaled interaction pressure again follows from eqn (32)–(36) and its behavior as a function of the rescaled inter-surface separation is shown in Fig. 8 for a few different sets of parameters in the case of semi-infinite slabs (b = ∞) with fixed Ξ = 50 and Δ = 0.95.

Fig. 8 Rescaled interaction pressure as a function of the rescaled distance between the randomly charged inner surfaces of two semi-infinite dielectric slabs with fixed Ξ = 50 and Δ = 0.95 and for (a) fixed [small kappa, Greek, tilde] = 0.3, [small chi, Greek, tilde]_c = 0.2 and different values of χ, (b) fixed [small kappa, Greek, tilde] = 0.3, χ = 2 and different values of [small chi, Greek, tilde]_c, and (c) fixed [small chi, Greek, tilde]_c = 0.2, χ = 2 and different values of [small kappa, Greek, tilde] as shown on the graphs. The inset in (a) shows a close-up view of the region around the shallow local minimum at large separations. Insets in (b) and (c) show the behavior of the rescaled interaction pressure over a wider range of scales.

In Fig. 8a, we show the results for a few different values of the disorder coupling parameter χ = 0, 1, 2 and 4 with [small kappa, Greek, tilde] = 0.3 and [small chi, Greek, tilde]_c = 0.2. With the choice of physical parameter values as q = 4 (tetravalent counterions), σ = 0.24 nm⁻² for water (ε_m = 80) at room temperature (T = 293 K), these rescaled parameter values correspond to surface charge disorder variances g = 0, 0.02, 0.04 and 0.08 nm⁻², respectively, and mono- and multivalent salt bulk concentrations of n₀ ≃ 150 mM and c₀ ≃ 4.4 mM (see Table 1).

For the non-disordered case with χ = 0 (black solid curve), the results show a repulsive, monotonically decaying interaction at small separations and a shallow attractive minimum at larger separations around d̃ ≃ 25 (inset) of the type that were analyzed in detail in Fig. 5 and 6. The presence of charge disorder on the inner surfaces of the slabs leads to significant qualitative differences in the effective interaction profile of the two surfaces (dashed curves). The disorder effects dominate at small to intermediate separations and turn the surface repulsion to a very strong attractive interaction. They diminish at larger separations as the curves with different values of χ (including χ = 0) converge. These features indicate an interplay between different contributions to the net interaction pressure that we shall examine later.

The attractive regime at small separations is followed by a repulsive regime with a pronounced hump at intermediate separations. For the set of parameter values used in Fig. 8a and for χ = 2 (blue dashed curve), the position of the hump and the maximum pressure are given by d̃ ≃ 11.2 and P̃ ≃ 0.05, which, for the set of the physical parameter values mentioned above, correspond to d ≃ 2.6 nm and P ≃ 0.5 bar. The shallow attractive minimum of the pressure curves at larger separations (d̃ ≃ 25 or d ≃ 5.7 nm) is only weakly influenced by the presence of disorder, as expected.

For sufficiently large disorder strengths (e.g., for χ = 4 shown by the green dashed curve), the attractive regime extends to the whole range of small to intermediate separations, even beyond the large-separation minimum, as shown in the inset, producing thus a long-ranged attractive interaction between the slabs. This follows as a consequence of the combined effects from the SC electrostatics of multivalent counterions and the surface charge disorder that couple to one another through the multivalent counterions contribution, P̃_c (see below), leading to features that are distinctly different from what we found in the case of uniformly charged surfaces (Section IV D).

A similar trend is found if the disorder strength is kept fixed and the rescaled bulk concentration of multivalent counterions, [small chi, Greek, tilde]_c, is increased (Fig. 8b). In this case both the range and the strength of the attractive pressure acting on the slabs are increased. We find a slightly different behavior when the Debye screening parameter (or equivalently, the bulk concentration of monovalent salt) is changed. For the parameters shown in Fig. 8c, we see that, as [small kappa, Greek, tilde] is decreased, the interaction pressure at intermediate separations first increases and becomes slightly more repulsive (compare green and blue dashed curves) and then turns to become attractive, in accordance with the intuitive expectation that the SC and disorder-induced effects become stronger at lower salt concentrations.

The disorder-induced attraction is enhanced and the repulsive hump and the large-separation minimum in the pressure curves disappear also when Δ → 0 (see Fig. 9), which clearly points to the key role of dielectric images, generating stronger repulsive interaction in systems with larger interfacial dielectric discontinuity. In the dielectrically homogeneous case (i.e., with Δ = 0 and [small kappa, Greek, tilde] > 0), we find long-ranged, monotonic attraction in the whole range of inter-surface separations with a diverging attractive pressure in the limit d̃ → 0 for all values of χ ≥ 0 (note that this behavior is in contrast with what we found in the canonical counterion-only case in Section IV A, where the pressure for d̃ → 0 becomes repulsive in the regime χ < 1). When the dielectric image effects are included (Δ > 0), the pressure turns repulsive again at small separations as, e.g., seen for χ = 1 in Fig. 8a and for [small kappa, Greek, tilde] = 0.32 in Fig. 8c; the other curves with larger χ or smaller [small kappa, Greek, tilde] in these figures, and also those shown in Fig. 8b, turn repulsive at extremely small separations, where the counterion–image repulsions eventually dominate and the multivalent counterions are fully depleted from the slit region (see also the discussion relating to Fig. 3b in Section IV B). However, the upturn of the pressure curves in these latter cases occurs at separations that are not physically meaningful (insofar as the experimental realizations of our model are concerned), e.g., around d̃ ≃ 0.02 (or d ≃ 0.005 nm) for χ = 2, the blue dashed curve, in Fig. 8a, and, thus, for the sake of presentation, it has not been shown in the graphs. This, on the other hand, means that, for sufficiently large χ and [small chi, Greek, tilde]_c and/or sufficiently small [small kappa, Greek, tilde], the stable equilibrium separation between the slabs (corresponding to the smallest point of zero pressure as, e.g., seen in Fig. 5a and 6a in the case of non-disordered surfaces and in Fig. 8a and c in the case of disordered surfaces) is pushed down to very small values, where the surfaces are nearly in contact. Indeed, in the case of disordered (but otherwise effectively like-charged) surfaces with no interfacial dielectric discontinuity and with counterions only, the disorder-induced attraction can become strong enough to cause a continuous transition to a collapsed state.^56,58

Fig. 9 Same as Fig. 8 but for fixed Ξ = 50, [small kappa, Greek, tilde] = 0.3, [small chi, Greek, tilde]_c = 0.2, χ = 2 and different values of Δ as shown on the graph.

The disorder-induced non-monotonicity found in the interaction profiles can be understood by examining the different components that contribute to the net pressure as defined in eqn (32)–(36). These are the mean surface–surface repulsion, P̃_σ, the interaction pressure due to the disorder variance, P̃_dis, the interaction pressure mediated by the multivalent counterions, P̃_c, and the monovalent ions contribution, P̃_mon. The effects of charge disorder, the dielectric and salt images are systematically included in all these components. In Fig. 10a, we show their behavior as a function of the rescaled inter-surface separation for Ξ = 50, Δ = 0.95, [small kappa, Greek, tilde] = 0.3, [small chi, Greek, tilde]_c = 0.2 and χ = 2. In order to enable a comparison between the non-disordered and disordered cases, we also show the corresponding curves of the non-disordered case with χ = 0 (marked in the legends with the superscript ‘0’). Note that P̃_σ is independent of the disorder strength and decays monotonically with the inter-surface distance. This contribution is comparable to that from the monovalent ions P̃_mon, which is also repulsive and shows a very weak dependence on the disorder strength with a slightly larger (more repulsive) value in the case of disordered surfaces as compared with non-disordered ones (P̃_mon > P̃⁰_mon). This is because, by increasing the disorder strength, more monovalent ions are attracted to the slit from the bulk, creating a larger osmotic pressure component. The effects of surface charge disorder on the multivalent counterions pressure component are quite substantial (compare P̃_c and P̃⁰_c) and most significant at small to intermediate separations, where the non-disordered contribution, P̃⁰_c (light-blue solid curve) shows a highly non-monotonic behavior with a repulsive (positive) hump as discussed in Section IV D, while the disorder contribution, P̃_c (red dashed curve), is attractive (negative) and increases in strength as the separation distance becomes small (until, at extremely small separations, where it turns repulsive and then decreases down to the bulk value for d̃ → 0; see below). Again, the disorder effects diminish and these two cases converge at separations larger than the Debye screening length.

Fig. 10 (a) Different components that contribute to the rescaled interaction pressure plotted as a function of the rescaled inter-surface distance between the slabs for the cases with disordered (χ = 2) and non-disordered (uniformly charged, χ = 0) surfaces. The curves that correspond to the non-disordered case are marked in the legends with the superscript ‘0’. The system parameters are Ξ = 50, Δ = 0.95, [small kappa, Greek, tilde] = 0.3 and [small chi, Greek, tilde]_c = 0.2. For comparison, we also show the vdW pressure between the slabs. (b) Different rescaled pressure components for the disordered case for the same parameter values as in (a) but with χ = 1, 2 and 4.

The other pressure component that needs to be considered here is P̃_dis that stems from the disorder variance and exists irrespective of the multivalent counterions in the system. This contribution is repulsive for Δ > 0 and becomes relevant only at very small separations. For the sake of comparison, we also show the inter-surface pressure generated by the van der Waals (vdW) interaction between two semi-infinite slabs immersed in an ionic mixture (dot-dashed curve). This interaction pressure can be calculated from the standard Lifshitz theory as (in actual units)⁸⁷

where the first term comes from the zero-frequency mode of the electromagnetic field-fluctuations and the second term comes from the higher-order Matsubara frequencies. A is the so-called Hamaker coefficient, which we take as A = 3 zJ (upper bound for the non-zero Matsubara modes in the case of hydrocarbon slabs interacting across an aqueous medium⁸⁸). For slabs of lower dielectric polarizability than the solvent, such as the cases considered here, the vdW interaction is attractive. In the figures, we show the rescaled quantity P̃_vdW = βP_vdW/(2π[small script l]_Bσ²) for the given parameter values. It is clear from Fig. 10a that the vdW component has a comparable range and strength as P̃_dis; these two contributions can thus compete at very small separations, and in the absence of multivalent counterions, generate a non-monotonic interaction profile between randomly charged dielectrics as shown in ref. 57, 59–63 and 74. These effects are however masked by the multivalent counterions contribution, P̃_c.

If χ is increased further (Fig. 10b), the contribution of monovalent ions P̃_mon changes only slightly (becoming more repulsive as noted before), but the effect of the increase of χ on both P̃_dis (which becomes more repulsive) and P̃_c (which becomes more attractive) is rather substantial. Note that P̃_dis depends linearly on the parameter χ according to eqn (34), while the dependence of P̃_c (and also P̃_mon) on χ is nontrivial and occurs through the effective single-particle energy, eqn (21).

We can thus conclude that the qualitative differences found between the interaction pressure curves in the non-disordered and disordered systems are closely connected with the behavior of the multivalent counterions contribution, P̃_c, that becomes significantly more attractive and dominant at small to intermediate separations when the disorder strength is increased, which explains the trends observed in Fig. 8a. We should emphasize that the increase in the strength of this attractive component is due to the fact that a larger amount of multivalent counterions are pulled into the slit region because of the stronger counterion–surface attractions in the presence of disorder, mediating also a stronger inter-surface attraction. This can be seen by inspecting the average number of multivalent counterions in the slit as shown in Fig. 11. This figure also shows that, for sufficiently large χ, the number of multivalent counterions in the slit increases by decreasing d̃, which implies that the attractive interaction between the surfaces increases in accord with Fig. 8a. However, as noted above, in systems with a finite dielectric discontinuity Δ > 0, the dielectric image repulsions dominate at small separations and the counterions are eventually ejected from the slit; therefore, in this case, N̄ → 0 for d̃ → 0 (not shown in the figure); consequently, we have P̃_c → −[small chi, Greek, tilde]_c²/4 (bulk value) and the total interaction pressure turns repulsive because of the other (repulsive) pressure components. (This is to be contrasted with the dielectrically homogeneous case with Δ = 0, where P̃_c and thus the total pressure diverge to large negative values for d̃ → 0, see Fig. 9.)

Fig. 11 Rescaled average number of multivalent counterions in the slit per rescaled surface area as a function of the rescaled distance between the randomly charged inner surfaces of two semi-infinite dielectric slabs for a few different values of the disorder coupling parameter, χ, for fixed [small kappa, Greek, tilde] = 0.3, [small chi, Greek, tilde]_c = 0.2, Ξ = 50 and Δ = 0.95.

A similar behavior is found for P̃_c when the rescaled bulk density of multivalent counterions is increased (Fig. 8b), and/or when the Debye screening length (Fig. 8c) or the dielectric discontinuity parameter (Fig. 9) are decreased. For instance, we show the behavior of the different pressure components for Δ = 0.25 (solid curves) and 0.95 (dashed curves) in Fig. 12. Only the two disorder-induced components, P̃_c and P̃_dis, show significant variations with the dielectric discontinuity parameter, Δ. The changes in the interaction pressure in Fig. 9 can be assigned primarily to the changes in the pressure component P̃_c as, e.g., at larger values of Δ, these ions are affected more strongly by image repulsions and are depleted more strongly from the slit region, giving a smaller attractive pressure, P̃_c.

Fig. 12 Same as Fig. 10a but for two different values of Δ = 0.25 and 0.95 (the total pressure for these cases is shown by black solid and black dashed curves, respectively).

V. Conclusion and discussion

In this paper, we have studied the effective interactions mediated by multivalent counterions between dielectric slabs that carry a quenched, random monopolar charge distribution on their juxtaposed, plane-parallel, inner surfaces. We have considered both a case in which the Coulomb fluid in the slit between the slabs consists only of mobile multivalent counterions, and also a case in which the Coulomb fluid is an asymmetric ionic mixture containing a q : 1 salt (with multivalent counterions of charge valency q > 1) and a monovalent 1 : 1 salt in equilibrium with a bulk reservoir. Our goal has been to elucidate the effects due to the coupling between the charge disorder and electrostatic correlations in the so-called strong coupling (SC) regime, realized experimentally when the system contains mobile multivalent counterions that give rise to strong surface–counterion correlations and also, to subleading orders, counterion–counterion correlations. The counter-intuitive phenomena in SC electrostatics, such as attraction between like-charged surfaces, have been well studied in the case of non-disordered (and, in most cases, uniformly charged) macromolecular surfaces (see, e.g., recent reviews in ref. 3–7 and references therein). The SC effects have, however, remained largely unexplored in the situation where the bounding (macromolecular) surfaces bear quenched, disordered charge distributions.

In the weak coupling (WC) regime, where, e.g., all ions in the Coulomb fluid are monovalent, the quenched charge disorder on bounding surfaces turns out to have either no or only very small effects and does not lead to any qualitatively new features in the behavior of the Coulomb fluid.^56,66,67 In the SC regime, the interaction of quenched random charge distributions across a Coulomb fluid has been considered only within the plane-parallel counterion-only model with no interfacial dielectric discontinuity.^56,58 It was shown that quenched surface charge disorder leads to a renormalization (reduction) of the entropic contributions to the interaction free energy (corresponding to a renormalized effective temperature), and that for sufficiently large disorder coupling parameter, i.e., χ ≥ 1, the surfaces undergo a continuous transition to a collapsed state.

In the present work, we have derived a generalized dressed multivalent-ion formalism that incorporates dielectric image effects as well as salt screening and salt image effects, being applicable to any arbitrary geometry of the bounding surfaces with disordered charge distributions. The charge disorder over the bounding surfaces is assumed to be Gaussian, without any spatial correlations, with possible generalizations to be discussed elsewhere.⁸⁹ The limiting case of this generalized model in the case of a single randomly charged dielectric slab has been discussed recently.⁶⁴

We have analyzed the prediction of the present theory in the specific example of two plane-parallel dielectric slabs with randomly charged inner surfaces immersed in an asymmetric Coulomb fluid. We have shown that, in a dielectrically homogeneous system (with or without additional salt), the multivalent counterions are attracted towards the surfaces with a singular, disorder-induced potential that diverges logarithmically in the vicinity of the surfaces; this creates a singular but integrable counterion density profile that exhibits an algebraic divergence at the surfaces with an exponent that depends on the disorder coupling parameter, χ. Remarkably, this behavior is in contrast with the contact-value theorem, which describes the behavior of counterions at uniformly charged surfaces and predicts a finite contact density.^80–82

Our results for the counterion-only case also shed further light on the previous findings for this system⁵⁶ that, as noted above, predicted a renormalized entropic contribution to the interaction free energy. We thus show here that such a behavior follows as a result of the singular accumulation of counterions in the immediate vicinity of the two surfaces, resulting in a renormalized temperature for the system. This notion of a renormalized temperature should be used with caution because first, this particular form of the renormalization of the system entropy is found only in the SC limit of the two-surface counterion-only model, and secondly, one can show that the renormalization in fact originates from both energetic and entropic sources. Nevertheless, the interplay between the translational entropy of strongly coupled counterions and the (non-thermal) configurational entropy, due to the averaging over different realizations of the quenched disorder, does result in the anti-fragile behavior⁷³ of multivalent counterions. Antifragility does not only result in a more ‘ordered’ state, characterized by a diminished intrinsic thermal ‘disorder’ induced by the externally imposed disorder, but also in a thermodynamically more stable state since quenched surface charge disorder engenders also a lower free energy as compared with a non-disordered system. The appearance of antifragility is possibly one of the most fundamental and perplexing features in the complicated Coulomb world, whose ramifications we are only beginning to unravel.

The singular behavior of the multivalent counterion density on approach to randomly charged surfaces persists also when the dielectrically homogeneous system contains an additional monovalent salt component. The mobile monovalent salt ions generate screening effects at intermediate to large distances (comparable to the Debye screening length), while at small to intermediate distances from the bounding surfaces, we find a narrow region where the multivalent counterions are partially depleted from the vicinity of the surfaces because of the salt image repulsions. This depletion effect arises because of the inhomogeneous distribution of the monovalent salt ions in the system (as the slabs are assumed to be impermeable to these ions) but its overall effect is quite weak and, when the surfaces are randomly charged, gives way to the singular attraction of counterions by disordered surface charges at small separations.

This salt image (or screening-induced) depletion is in contrast with the depletion due to dielectric images, which in dielectrically inhomogeneous systems always dominates in the vicinity of the dielectric surfaces and creates an interfacial zone of complete depletion, even when the bounding surfaces of the slabs are randomly charged, the reason being the repulsive ion-surface potential generated by the dielectric images that diverges at a dielectric surface (with Δ > 0) in a way that overcomes the disorder-induced attraction experienced by the multivalent counterions at small distances from the surfaces. The amount of multivalent counterions accumulated at some short distance away from the dielectric surface is, however, still overall enhanced because of the disorder-induced attraction.

We have also analyzed in detail the consequences that result from the interplay between charge disorder, dielectric and salt images, and the SC effects on the effective interaction pressure between the slabs. The interaction pressure shows in general a highly non-monotonic behavior as a function of the separation distance between the inner surfaces of the slabs. At small to intermediate separations (e.g., typically within a few nanometers), the SC effects can be quite significant. They contribute an attractive component to the net interaction pressure that opposes the repulsive osmotic pressure due to monovalent salt ions and the repulsive interaction between the mean surface charges on the two slabs. In the absence of disorder, this SC attraction between the two like-charged surfaces becomes dominant when the bulk concentration of multivalent counterions is increased and/or when the salt screening parameter or the dielectric discontinuity parameter is decreased. Nevertheless, when the inter-surface separation decreases, the pressure at very small separations becomes repulsive once again since counterions are completely depleted from the slit due to dielectric image repulsions. In the presence of disorder, counterions are attracted to the surface much more strongly than in the absence of disorder and a much larger number of counterions are sucked into the slab region from the bulk. As a result, the SC attraction mediated by multivalent counterions between the disordered charged surfaces becomes increasingly more enhanced, especially as the disorder coupling parameter and/or the bulk concentration of multivalent counterions are increased. The repulsive regime at small separations, which could stabilize the surfaces in the absence of disorder, is thus squeezed down to zero due to extremely large attractive pressures acting on the slabs.

We have presented our results in a rescaled (dimensionless) form, in terms of a few dimensionless parameters such as the rescaled screening parameter and the electrostatic and disorder coupling parameters. These parameters, for any given set of values, can be mapped to a wide range of values for the actual parameters, namely, the counterion and salt bulk concentrations, mean surface charge density, counterion valency, solvent dielectric constant, temperature, etc. A few typical examples of the physical parameter values that can correspond to some typical values of the rescaled parameters are shown in Table 1; other sets of physical parameter values are conceivable, e.g., by using divalent and trivalent counterions or other surface charge densities, etc. It should be noted that the range of values we have given in the table for the disorder coupling parameter, i.e., χ ≃ 0.5–4, corresponds to typically small degrees of charge disorder with disorder variances of around g ≃ 0.01–0.08 nm⁻² that can be achieved by relatively small surface density of impurity charges (≲0.1 nm⁻²) as compared with the mean number of surface charges (typically σ ≲ 1 nm⁻²). While larger values of the disorder variance are equally permissible, our results show that the effects due to the coupling between surface charge randomness and the SC electrostatics due to mobile multivalent counterions can be quite significant even at very small degrees of surface charge disorder.

The regime of applicability of the dressed multivalent-ion theory (which can be justified strictly only in the case of highly asymmetric Coulomb fluids³⁴) has been discussed in detail elsewhere^7,34–38 by making extensive comparison with implicit- and explicit-ion simulations in non-disordered systems, where it is found to be in the experimentally accessible parameter space. This includes the situations where electrostatic interactions are strong enough so that the effects of multivalent counterions next to an oppositely charged boundary are adequately included on the lowest-order single-particle level.^{5–7,14–18} The higher-order effects of multi-particle interactions between counterions are assumed to be sufficiently weak for moderate to high salt concentrations.^5,7,34–36 In particular, the present approach is expected to be applicable for relatively small (a few mM) bulk concentrations of the multivalent counterions, which is in fact the typical case in most experiments (see, e.g., ref. 20–31). We expect that the previously determined regimes of validity hold also for the present case with randomly charged surfaces, although this remains to be determined more systematically through extensive computer simulations that are still missing for the disordered systems.

Our results yield concrete predictions that can be tested against these future simulations. The fingerprints of disorder effects can show up in surface-force measurements as well (see, e.g., ref. 2, 47, 48, 53 and 59–63 and references therein for charge disorder effects in the context of Casimir force measurements and ref. 41–45 for force measurement between surfactant-coated surfaces in aqueous media). One should, however, note that the precise statistical distribution of charge disorder in real systems can be sample and material dependent and it can vary depending also on the method of preparation. These points need to be addressed first if the theoretical predictions are to be compared against experiments.

The present study is based on a primitive model of Coulomb fluids and makes a few simplifying assumptions that have been discussed elsewhere in detail (see, e.g., ref. 7, 34–36 and 64). For instance, we have neglected the solvent structure (see, e.g., ref. 1, 10 and 90–94 and references therein), the polarizability of mobile ions (see, e.g., ref. 95 and 96 and references therein), specific surface ion-adsorption effects,⁹⁷ and the size and internal structure of the counterions.⁹⁸ For multivalent counterions that can be modeled as spherical particles, their finite size can be incorporated in our approach in a straightforward manner as discussed in detail in ref. 64. It is important to note that most multivalent counterions possess an internal charge distribution that can introduce higher-order multipolar effects; these effects can be quite significant especially for multivalent counterions that have an extended structure (see, e.g., ref. 7 and 99 and references therein), including rod-like polyamines such as the trivalent spermidine and tetravalent spermine that have chain lengths of up to 1–1.5 nm.¹⁰⁰ We plan to address these issues in future publications. Other cases that can be studied in the present context include spatially correlated surface charge disorder^60,63,89 as well as surfaces with annealed (mobile) disordered charges,^{50,59,60,66–71} surfaces with partially quenched or partially annealed charges,^58,72 and also charge regulating surfaces.^71,101–105

We should also note that our approach is based on a continuum approximation for the surface charge distribution. It would be interesting to include other features of the detailed molecular nature of the surface charge such as the surface charge discreteness (see, e.g., ref. 66, 67, 84 and 106–108 and references therein). However, charge discreteness alone will not improve the reliability of results and concurrently, one would need to consistently include at least the discreteness of solvent, the interfacial dielectric constant profile, etc. Our approach is therefore based on the level that is presently tractable. Specifically, the charge discreteness effects alone introduce a new length scale into the system and are in some respects similar to the correlated disorder effects,^60,63 which we have analyzed elsewhere,⁸⁹ implying that while they would quantitatively change the details of our predictions they will not interfere with their qualitative features, which are after all the main focus of our endeavors.

Another point to be noted is that, in systems containing an added monovalent salt, our theoretical approach may become invalid when the mean electrostatic potential near the randomly charged surfaces becomes large, in which case the validity of the underlying DH approximation used for the monovalent ions can break down.³⁵ This can happen when the disorder strength is large and/or when the dielectric discontinuity parameter is small. Other cases that go beyond the present framework include the situation where nonlinear charge renormalization and/or Bjerrum pairing effects become relevant (see, e.g., ref. 109–113). These latter effects, however, turn out to be negligible in the regime of parameters that is of interest here.^7,34–36 These and other possible issues such as higher-order virial corrections,^{5–7,14–16,18} the intermediate coupling effects that will be relevant in a wide range of realistic systems (see, e.g., ref. 5–7, 16 and 114–117), the discrete nature of monovalent salt ions,³⁵ ion–ion excluded-volume repulsions,^118–122etc., that may be relevant especially at intermediate electrostatic couplings and/or large multivalent counterion concentrations, remain to be elucidated in future simulations.

Acknowledgements

A.N. acknowledges partial support from the Royal Society, the Royal Academy of Engineering, and the British Academy (UK). H.K.M. acknowledges support from the School of Physics, Institute for Research in Fundamental Sciences (IPM), where she stayed as a short-term visiting researcher during the completion of this work. R.P. acknowledges support from the ARRS through Grants no. P1-0055 and J1-4297. We thank E. Mahgerefteh for useful comments and discussions.

References

1. Electrostatic Effects in Soft Matter and Biophysics, ed. C. Holm, P. Kekicheff and R. Podgornik, Kluwer Academic, Dordrecht, 2001.
2. Electrostatics of Soft and Disordered Matter, ed. D. S. Dean, J. Dobnikar, A. Naji and R. Podgornik, Pan Stanford Publishing, Singapore, 2014.
3. A. Yu Grosberg, T. T. Nguyen and B. I. Shklovskii, Rev. Mod. Phys., 2002, 74, 329.
4. Y. Levin, Rep. Prog. Phys., 2002, 65, 1577.
5. H. Boroudjerdi, Y. W. Kim, A. Naji, R. R. Netz, X. Schlagberger and A. Serr, Phys. Rep., 2005, 416, 129.
6. A. Naji, S. Jungblut, A. G. Moreira and R. R. Netz, Phys. Stat. Mech. Appl., 2005, 352, 131.
7. A. Naji, M. Kanduč, J. Forsman and R. Podgornik, J. Chem. Phys., 2013, 139, 150901.
8. R. French, et al. , Rev. Mod. Phys., 2010, 82, 1887.
9. E. J. Verwey and J. T. G. Overbeek, Theory of the Stability of Lyophobic Colloids, Elsevier, Amsterdam, 1948.
10. J. Israelachvili, Intermolecular and Surface Forces, Academic Press, London, 1991.
11. S. A. Safran, Statistical Thermodynamics of Surfaces, Interfaces, and Membranes, Addison-Wesley Publishing Co., Massachusetts, 1994.
12. R. J. Hunter, Foundations of Colloidal Science, Oxford University Press, USA, 2nd edn, 2001.
13. Soft Condensed Matter Physics in Molecular and Cell Biology, ed. W. C. K. Poon and D. Andelman, Taylor & Francis, New York, London, 2006.
14. R. R. Netz, Eur. Phys. J. E, 2001, 5, 557.
15. A. G. Moreira and R. R. Netz, Phys. Rev. Lett., 2001, 87, 078301.
16. A. G. Moreira and R. R. Netz, Eur. Phys. J. E, 2002, 8, 33.
17. A. Naji and R. R. Netz, Eur. Phys. J. E, 2004, 13, 43.
18. M. Kanduč, M. Trulsson, A. Naji, Y. Burak, J. Forsman and R. Podgornik, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 061105.
19. L. Šamaj and E. Trizac, Phys. Rev. Lett., 2011, 106, 078301.
20. V. A. Bloomfield, Curr. Opin. Struct. Biol., 1996, 6, 334.
21. J. Pelta, D. Durand, J. Doucet and F. Livolant, Biophys. J., 1996, 71, 48.
22. J. Pelta, F. Livolant and J.-L. Sikorav, J. Biol. Chem., 1996, 271, 5656.
23. K. Yoshikawa, Adv. Drug Delivery Rev., 2001, 52, 235.
24. M. Takahashi, K. Yoshikawa, V. V. Vasilevskaya and A. R. Khokhlov, J. Phys. Chem. B, 1997, 101, 9396.
25. G. E. Plum and V. A. Bloomfield, Biopolymers, 1988, 27, 1045.
26. E. Raspaud, I. Chaperon, A. Leforestier and F. Livolant, Biophys. J., 1999, 77, 1547.
27. H. S. Savithri, S. K. Munshi, S. Suryanarayana, S. Divakar and M. R. N. Murthy, J. Gen. Virol., 1987, 68, 1533.
28. M. de Frutos, S. Brasiles, P. Tavares and E. Raspaud, Eur. Phys. J. E, 2005, 17, 429.
29. A. Leforestier, A. Siber, F. Livolant and R. Podgornik, Biophys. J., 2011, 100, 2209.
30. D. C. Rau and V. A. Parsegian, Biophys. J., 1992, 61, 246.
31. D. C. Rau and V. A. Parsegian, Biophys. J., 1992, 61, 260.
32. T. E. Angelini, H. Liang, W. Wriggers and G. C. L. Wong, Proc. Natl. Acad. Sci. U. S. A., 2003, 100, 8634.
33. D. J. Needleman, M. A. Ojeda-Lopez, U. Raviv, H. P. Miller, L. Wilson and C. R. Safinya, Proc. Natl. Acad. Sci. U. S. A., 2004, 101, 16099.
34. M. Kanduč, A. Naji, J. Forsman and R. Podgornik, J. Chem. Phys., 2010, 132, 124701.
35. M. Kanduč, A. Naji, J. Forsman and R. Podgornik, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 011502.
36. M. Kanduč, A. Naji, J. Forsman and R. Podgornik, J. Chem. Phys., 2012, 137, 174704.
37. L. Javidpour, A. Lošdorfer Božič, A. Naji and R. Podgornik, J. Chem. Phys., 2013, 139, 154709.
38. L. Javidpour, A. Lošdorfer Božič, A. Naji and R. Podgornik, Soft Matter, 2013, 9, 11357.
39. S. Manne, J. P. Cleveland, H. E. Gaub, G. D. Stucky and P. K. Hansma, Langmuir, 1994, 10, 4409.
40. S. Manne and H. E. Gaub, Science, 1995, 270, 1480.
41. E. E. Meyer, Q. Lin, T. Hassenkam, E. Oroudjev and J. N. Israelachvili, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 6839.
42. E. E. Meyer, K. J. Rosenberg and J. Israelachvili, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 15739.
43. S. Perkin, N. Kampf and J. Klein, Phys. Rev. Lett., 2006, 96, 038301.
44. S. Perkin, N. Kampf and J. Klein, J. Phys. Chem. B, 2005, 109, 3832.
45. G. Silbert, D. Ben-Yaakov, Y. Dror, S. Perkin, N. Kampf and J. Klein, Phys. Rev. Lett., 2012, 109, 168305.
46. H. T. Baytekin, A. Z. Patashinski, M. Branicki, B. Baytekin, S. Soh and B. A. Grzybowski, Science, 2011, 333, 308.
47. W. J. Kim, A. O. Sushkov, D. A. R. Dalvit and S. K. Lamoreaux, Phys. Rev. Lett., 2009, 103, 060401.
48. W. J. Kim, A. O. Sushkov, D. A. R. Dalvit and S. K. Lamoreaux, Phys. Rev. A, 2010, 81, 022505.
49. Y. Kantor, H. Li and M. Kardar, Phys. Rev. Lett., 1992, 69, 61.
50. I. Borukhov, D. Andelman and H. Orland, Eur. Phys. J. B, 1998, 5, 869.
51. S. J. Miklavic, D. Y. C. Chan, L. R. White and T. W. Healy, J. Phys. Chem., 1994, 98, 9022.
52. S. Panyukov and Y. Rabin, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 56, 7055.
53. C. C. Speake and C. Trenkel, Phys. Rev. Lett., 2003, 90, 160403.
54. D. B. Lukatsky, K. B. Zeldovich and E. I. Shakhnovich, Phys. Rev. Lett., 2006, 97, 178101.
55. D. B. Lukatsky and E. I. Shakhnovich, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 020901.
56. A. Naji and R. Podgornik, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 041402.
57. R. Podgornik and A. Naji, Europhys. Lett., 2006, 74, 712.
58. Y. Sh. Mamasakhlisov, A. Naji and R. Podgornik, J. Stat. Phys., 2008, 133, 659.
59. A. Naji, D. S. Dean, J. Sarabadani, R. Horgan and R. Podgornik, Phys. Rev. Lett., 2010, 104, 060601.
60. J. Sarabadani, A. Naji, D. S. Dean, R. R. Horgan and R. Podgornik, J. Chem. Phys., 2010, 133, 174702.
61. D. S. Dean, A. Naji and R. Podgornik, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 011102.
62. A. Naji, J. Sarabadani, D. S. Dean and R. Podgornik, Eur. Phys. J. E, 2012, 35, 24.
63. V. Rezvani, J. Sarabadani, A. Naji and R. Podgornik, J. Chem. Phys., 2012, 137, 114704.
64. A. Naji, M. Ghodrat, H. Komaie-Moghaddam and R. Podgornik, J. Chem. Phys., 2014, 141, 174704.
65. D. Ben-Yaakov, D. Andelman and H. Diamant, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 022402.
66. C. Fleck and R. R. Netz, Europhys. Lett., 2005, 70, 341.
67. C. Fleck and R. R. Netz, Eur. Phys. J. E: Soft Matter Biol. Phys., 2007, 22, 261.
68. A. Naydenov, P. A. Pincus and S. A. Safran, Langmuir, 2007, 23, 12016.
69. R. Brewster, P. A. Pincus and S. A. Safran, Phys. Rev. Lett., 2008, 101, 128101.
70. Y. S. Jho, R. Brewster, S. A. Safran and P. A. Pincus, Langmuir, 2011, 27, 4439.
71. Y. S. Velichko and M. Olvera de la Cruz, J. Chem. Phys., 2006, 124, 214705.
72. B. Hribar-Lee, M. Lukšič and V. Vlachy, Annu. Rep. Prog. Chem., Sect. C: Phys. Chem., 2011, 107, 14.
73. N. N. Taleb, Antifragile: Things that Gain from Disorder, Random House, 2012.
74. M. Ghodrat, A. Naji, H. Komaie-Moghaddam and R. Podgornik, “Counterion-mediated non-monotonic interactions between net-neutral randomly charged dielectrics”, manuscript in preparation.
75. S. F. Edwards and A. Lenard, J. Math. Phys., 1962, 3, 778.
76. R. Podgornik and B. Žekš, J. Chem. Soc. Faraday Trans. II, 1988, 84, 611.
77. R. Podgornik, J. Phys. A: Math. Gen., 1990, 23, 275.
78. R. R. Netz and H. Orland, Eur. Phys. J. E, 2000, 1, 203.
79. S. Tsonchev, R. D. Coalson and A. Duncan, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 041804.
80. D. Henderson and L. Blum, J. Chem. Phys., 1981, 75, 2025.
81. S. L. Carnie and D. Y. C. Chan, J. Chem. Phys., 1981, 74, 1293.
82. H. Wennerström, B. Jönsson and P. Linse, J. Chem. Phys., 1982, 76, 4665.
83. D. S. Dean and R. R. Horgan, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 061106.
84. W. Pezeshkian, N. Nikoofard, D. Norouzi, F. Mohammad-Rafiee and H. Fazli, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 061925.
85. R. Podgornik, J. Phys. Chem., 1991, 95, 5249.
86. This refers to the isotherms as opposed to van der Waals interactions to be considered later.
87. V. A. Parsegian, van der Waals Forces: A Handbook for Biologists, Chemists, Engineers, and Physicists, Cambridge University Press, 2005.
88. R. Podgornik, R. H. French and V. A. Parsegian, J. Chem. Phys., 2006, 124, 044709.
89. E. Mahgerefeth, A. Naji and R. Podgornik, “Strong coupling interaction between patchy charge patterns”, manuscript in preparation.
90. D. Ben-Yaakov, D. Andelman, D. Harries and R. Podgornik, J. Phys.: Condens. Matter, 2009, 21, 424106.
91. D. Ben-Yaakov, D. Andelman and R. Podgornik, J. Chem. Phys., 2011, 134, 074705.
92. D. Ben-Yaakov, D. Andelman, R. Podgornik and D. Harries, Curr. Opin. Colloid Interface Sci., 2011, 16, 542.
93. Y. Burak and D. Andelman, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2000, 62, 5296.
94. Y. Burak and D. Andelman, J. Chem. Phys., 2001, 114, 3271.
95. J. J. Bikerman, Philos. Mag., 1942, 33, 384.
96. V. Démery, D. S. Dean and R. Podgornik, J. Chem. Phys., 2012, 137, 174903.
97. J. Forsman, Langmuir, 2006, 22, 2975.
98. D. Frydel, “Mean-field electrostatics beyond the point-charge description”, e-print: arXiv:1411.7577.
99. M. Kanduč, A. Naji, Y. S. Jho, P. A. Pincus and R. Podgornik, J. Phys.: Condens. Matter, 2009, 21, 424103.
100. M. L. Di Paolo, R. Stevanato, A. Corazza, F. Vianello, L. Lunelli, M. Scarpa and A. Rigo, Biochem. J., 2003, 371, 549.
101. B. W. Ninham and V. A. Parsegian, J. Theor. Biol., 1971, 31, 405.
102. D. Chan, J. W. Perram, L. R. White and T. H. Healy, J. Chem. Soc., Faraday Trans. 1, 1975, 71, 1046.
103. H. H. von Grünberg, J. Colloid Interface Sci., 1999, 219, 339.
104. G. S. Longo, M. Olvera de la Cruz and I. Szleifer, Soft Matter, 2012, 8, 1344.
105. N. Adžić and R. Podgornik, Eur. Phys. J. E, 2014, 37, 49.
106. A. G. Moreira and R. R. Netz, Europhys. Lett., 2002, 57, 911.
107. M. O. Khan, S. Petris and D. Y. C. Chan, J. Chem. Phys., 2005, 122, 104705.
108. J. M. A. Grime, M. O. Khan and K. Bohinc, Langmuir, 2010, 26, 6343.
109. R. Kjellander, Colloid J., 2007, 69, 20.
110. S. Alexander, P. M. Chaikin, P. Grant, G. Morales, P. Pincus and D. Hone, J. Chem. Phys., 1984, 80, 5776.
111. N. Bjerrum, Kgl. danske Vid. Selsk. Math.-fys. Medd., 1926, 7, 1.
112. J. Zwanikken and R. van Roij, J. Phys.: Condens. Matter, 2009, 21, 424102.
113. M. E. Fisher and Y. Levin, Phys. Rev. Lett., 1993, 71, 3826.
114. Y. Burak, D. Andelman and H. Orland, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 70, 016102.
115. C. D. Santangelo, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 73, 041512.
116. M. M. Hatlo and L. Lue, Soft Matter, 2009, 5, 125.
117. J. M. Rodgers, C. Kaur, Y.-G. Chen and J. D. Weeks, Phys. Rev. Lett., 2006, 97, 097801.
118. R. Messina, E. González-Tovar, M. Lozada-Cassou and C. Holm, Europhys. Lett., 2002, 60, 383.
119. M. Quesada-Pérez, E. González-Tovar, A. Martin-Molína, M. Lozada-Cassou and R. Hidalgo-Álvarez, ChemPhysChem, 2003, 4, 234.
120. R. Kjellander, J. Phys.: Condens. Matter, 2009, 21, 424101.
121. G. I. Guerrero-Garcia, E. González-Tovar and M. Olvera de la Cruz, Soft Matter, 2010, 6, 2056.
122. D. Frydel and Y. Levin, J. Chem. Phys., 2012, 137, 164703.

---