Symmetry breaking of interacting charge-regulated planar macroions at low salt concentration

Abstract

The charge of macroions such as colloids and membranes can often be regulated by reversible ion adsorption or dissociation reactions. When two chemically identical macroions interact across an electrolyte, charge regulation will typically lead to the same charge density on both macroion surfaces. However, it was recently demonstrated that a non-ideal, Frumkin-like adsorption behavior of the ions can lead to a spontaneous symmetry breaking, where the coupling of charge regulation leads to different surface charge densities on both macroions in thermal equilibrium. While previous modeling was carried out numerically in the presence of salt for selected systems, we argue that the no-salt limit captures experimentally relevant situations at low salt concentration. We therefore focus on macroions with dissociable surface groups in the absence of salt, where the Poisson–Boltzmann equation can be solved analytically and the spontaneous emergence of uniform but different degrees of dissociation on both macroions signifies a symmetry breaking. Using a perturbation approach we derive analytic results for the local stability of the symmetric state. This not only provides a complete thermodynamic characterization of the symmetry breaking as function of all parameters, it also uncovers previously unrecognized features. First, depending on the degree of non-ideality, the symmetry breaking transition can be continuous or discontinuous. Second, metastable states do exist far away from critical points but not in their vicinity. And third, electrostatic interactions generally act toward weakening or suppressing the occurrence of symmetry breaking.

---

1 Introduction

Experimental studies suggest the complex surface chemistry of colloidal particles¹ can alter mutual interactions² and even trigger self-assembly via long-range forces, especially under low salt conditions.^3–5 Electrostatic modeling through Poisson–Boltzmann theory generally predicts repulsion between identical macroions,⁶ but when adjusting solvation properties^1,7 or charge regulation^2,8,9 at the surfaces of the involved macroions introduce additional order parameters, unexpected new properties can emerge even on the level of mean-field theory,^10,11 including non-monotonic screening lengths, buffering of the background electrolyte, and the possibility of long-range attraction between nominally like-charged objects. Not surprisingly, there has been considerable interest in modeling interacting macroions that are subject to a charge regulation mechanism using continuum electrostatics and computer simulations.^12–18

Charge regulation, the adjustment of a surface's effective charge through reversible association and dissociation of ions,¹⁹ is based on an underlying adsorption isotherm. Of particular interest is charge regulation through Frumkin-type adsorption,²⁰ which goes beyond the Langmuir adsorption model of independent adsorption sites by accounting for lateral interactions between already adsorbed ions. The Frumkin isotherm, also referred to as Frumkin–Fowler–Guggenheim isotherm, is frequently used to model the adsorption of contaminants such as heavy metals, dyes, and organics on porous carbons, clays, and oxides^21–24 but also for the adsorption of ionic surfactants and other amphiphiles on fluid interfaces,^25–28 the adsorption of hydrogen on platinum nanoparticles,²⁹ or the intercalation of lithium into inorganic host materials.³⁰ The nature and magnitude of lateral interactions among adsorbents has been estimated based on Frumkin-type adsorption isotherms for a multitude of systems such as the adsorption of surfactin on the liquid–air interface,³¹ ion-exchange and charge-regulation phenomena at surfaces of metal oxide particles immersed in aqueous electrolyte solutions,³² and the competitive adsorption of ions on mineral surfaces.³³

In a series of studies, Podgornik and coworkers^34–37 have proposed the occurrence of a spontaneous symmetry breaking for interacting, charge-regulated surfaces of chemically identical macroions. The symmetry breaking is driven by an underlying Frumkin–Fowler–Guggenheim isotherm that favors high and low effective charge densities over intermediate ones. Models for both planar³⁴ and spherical³⁷ macroions were analyzed in the presence of monovalent salt, and numerical solutions of the Poisson–Boltzmann equation were used to predict when symmetry breaking occurs. Simple analytic solutions of interacting macroions are, in general, not available in the presence of salt. This makes the complete characterization of the symmetry breaking as function of all parameters challenging. However, we argue that experimental studies of interacting, charge-regulated macroions such as colloidal particles or charged lipid membranes are often carried out at sufficiently low salt concentrations, at which the presence of salt can be ignored altogether and interacting macroions can be modeled under no-salt conditions. For monovalent salt, this is the case when the Debye screening length l_D ≫ a/l_B is much larger than the ratio of average cross-sectional area per charge a divided by the Bjerrum length l_B. The advantage of the no-salt regime is the availability of simple analytic solutions of the non-linear Poisson–Boltzmann equation for interacting planar surfaces. The opportunity to completely and accurately characterize the stability of a highly nonlinear system that previously was amenable only to numerical analysis has motivated the present study.

In this work, we model the interaction between planar charge-regulated macroions in the absence of salt, where the Poisson–Boltzmann equation offers an analytic solution. As in preceding models by Podgornik and coworkers,^10,34,37 we consider a charge regulation mechanism based on the Frumkin–Fowler–Guggenheim isotherm, which is equivalent to the mean-field Bragg–Williams free energy³⁸ in the way to account for non-electrostatic interactions between adsorbed ions. In contrast to the preceding studies, we assume surface groups become charged upon dissociation without the presence of additional adsorption sites for the counterions. Solutions of the Poisson–Boltzmann equation entail boundary conditions in the form of non-linear algebraic relationships. Using a perturbation approach we carry out a stability analysis, from which we obtain analytic predictions for the appearance of symmetry breaking. We demonstrate that broken symmetry exists only for intermediate distances between the macroions, but not for sufficiently small or large distances. In addition, we predict the existence of metastable symmetric states and present a stability diagram that fully describes the thermodynamic behavior of our model as function of all parameters. Finally, we discuss the applicability of the no-salt assumption to experimental situations, assess the role of electrostatic interactions for the appearance of spontaneous symmetry breaking, characterize the stability limit of an isolated macroion, and calculate the disjoining pressure between the two macroions.

2 Theory

We represent two interacting macroions by two very large planar and parallel surfaces of lateral area A, indexed “1” and “2” and separated by a distance d. Each macroion carries a fixed bare surface charge density σ = −e/a of negative sign, where a denotes the average cross-sectional area per elementary charge e. The net surface charge density of the macroion is equal to σ only when all counterions are released and none is located at the macroion surface. Sandwiched between the two planar macroions is an electrolyte that contains counterions only. The total number of counterions is N = 2A/a. Allowing for charge regulation, we assume counterions can be associated with or dissociated from the macroion surfaces. If N₁ = ϕ₁A/a and N₂ = ϕ₂A/a counterions are associated with macroions “1” and “2”, respectively, the remaining N_f = N − N₁ − N₂ counterions are mobile, thus forming a diffuse counterion cloud between the two surfaces. Counterion association reduces the net surface charge densities from σ to the effective values σ₁ = −(1 − ϕ₁)e/a and σ₂ = −(1 − ϕ₂)e/a, where ϕ₁ and ϕ₂ describe the degree of counterion association on each macroion, with 0 ≤ ϕ₁ ≤ 1 and 0 ≤ ϕ₂ ≤ 1. If ϕ₁ = ϕ₂ = 0, all counterions reside in the electrolyte and σ₁ = σ₂ = σ has the largest possible magnitude. If ϕ₁ = ϕ₂ = 1, all counterions are associated with the two macroions and σ₁ = σ₂ = 0 vanishes. Fig. 1 shows a schematic representation of the two planar macroions, with some counterions being macroion-associated and the remaining ones dissociated, forming a diffuse cloud of mobile counterions.

Fig. 1 Two planar macroions (shown in blue) carry a bare surface charge density σ, with counterions either associated with the surfaces (shown in dark red) or dissociated and thus being mobile in solution between the two macroions (shown in red). No salt is added. The macroion surfaces are located at positions x = 0 and x = d, with the x-coordinate pointing normal to the surfaces. The fractions ϕ₁ and ϕ₂ characterize the degrees of counterion association on macroions “1” and “2”, respectively. The local concentration of the intervening mobile counterions is denoted by n(x), and the dimensionless electrostatic potential Ψ(x) is shown schematically.

When the distance d between the planar macroions changes, ϕ₁ and ϕ₂, as well as the local concentration of the intervening mobile counterions (denoted below by n), can adjust. Characterization of this adjustment is the goal of the present work. We point out an important assumption: while the degrees of counterion association ϕ₁ and ϕ₂ can be different on each macroion, they are assumed to be uniform along each individual macroion surface. At least for isolated macroions we will make sure in this work that the degrees of counterion association remain uniform. For interacting macroions, however, the uniformity stands as an assumption.

We assume counterion association with a macroion surface is driven by a fixed association energy λ per counterion. In this work, all energies are expressed in units of the thermal energy k_BT, Boltzmann constant times absolute temperature. Hence, λ is a dimensionless quantity. We also consider a non-ideal contribution to the counterion association/dissociation isotherm and describe this by using the Bragg–Williams free energy³⁸ per charged site on each macroion surface,

f_BW(ϕ) = ϕ ln ϕ + (1 − ϕ)ln(1 − ϕ) + χϕ(1 − ϕ), (1)

where χ is a non-ideality parameter. Positive χ amounts to an effective non-electrostatic attraction between surface-associated counterions. The present work studies the consequences of χ > 0 without investigating the molecular origin that could give rise to attractive interactions.¹ We note that our counterion association/dissociation model based on the fixed association energy λ per counterion together with the non-ideality parameter χ in eqn (1) is equivalent to the Frumkin–Fowler–Guggenheim isotherm²⁸ as it was used in preceding work.^34,37 For χ = 0, counterion association/dissociation follows the simple Langmuir adsorption model of an electrified surface.

We express the total free energy of the two macroions with their counterions either being macroion-associated or dissociated and thus residing in the electrolyte as

The first line describes the non-ideal mixing free energies of the macroion-associated counterions together with their association energies. This involves the two parameters χ and λ. The second line expresses the free energy of the electrolyte, which occupies the space along the x-axis of a Cartesian coordinate system from the surface of macroion “1” (located at x = 0) to macroion “2” (located at x = d). All properties of our mean-field model depend only on the x-coordinate, including the dimensionless electrostatic potential Ψ = Ψ(x) and the local counterion concentration n = n(x). The term Ψ′²/8πl_B in the integrand characterizes the volume density of the electrostatic energy and l_B the Bjerrum length, thereby assuming a uniform dielectric constant in the solvent space between the macroions. Note that a prime denotes the derivative with respect to the argument; i.e. Ψ′ = dΨ(x)/dx. The term n ln(nν) − n accounts for the local free energy density of an ideal gas, where the “thermal” volume ν depends only on the temperature T. Keeping T fixed renders ν a constant. If we define the average counterion density in the electrolyte by , the conservation constraint of the total number of counterions, ϕ₁ + ϕ₂ + adn̄ = 2, can be incorporated into the free energy by defining

where Λ is a Lagrangian multiplier. Accounting for Poisson's equation Ψ″ = −4πl_Bn allows us to calculate the first variation of F̃ with respect to ϕ₁, ϕ₂, and n,

Stationarity demands δF̃ = 0. The variations δϕ₁, δϕ₂, and δn are arbitrary, thus implying the counterion association/dissociation equilibrium conditions

at the two macroion surfaces (recall that the prime in denotes the derivative with respect to the argument) as well as the Boltzmann distribution n = e^−Ψ−Λ/ν.

Note that the association/dissociation reaction S⁻ + A⁺ ⇌ SA, which we may consider taking place on macroion surface “1” (located at x = 0), implies surface fractions of neutral and anionic sites [SA] = ϕ₁ and [S⁻] = 1 − ϕ₁, respectively, as well as a local counterion concentration at the surface [A⁺] = n(0). Inserting Λ = −Ψ(0) − ln[νn(0)] from the Boltzmann distribution at x = 0 into eqn (5) and assuming χ = 0, yields the equilibrium constant

of the association/dissociation reaction. Eqn (6) is valid in the presence or absence of electrostatic interactions, and their absence entails n(0) = (1 − ϕ)/(ad) is uniform. If non-ideal behavior is present (χ ≠ 0), the surface fractions [SA] = ϕ₁ and [S⁻] = 1 − ϕ₁ are to be multiplied by corresponding activity coefficients.

Inserting the Boltzmann distribution n = e^−Ψ−Λ/ν into Poisson's equation yields the Poisson–Boltzmann equation Ψ″ = −(κ²/2) × e^−Ψ, where we have defined κ through κ² = 8πl_Be^−Λ/ν. The solution of the Poisson–Boltzmann equation for planar geometry in the absence of salt can be expressed as

The length δ will vanish for a symmetric state where ϕ₁ = ϕ₂, but not when the symmetry is broken and ϕ₁ ≠ ϕ₂. From eqn (7) we can calculate Ψ(0) and Ψ(d), which appear in eqn (5) as well as Ψ′(0) and Ψ′(d), which enter the electrostatic boundary conditions implied by Poisson's equation, Ψ′(0)/(4πl_B) = (1 − ϕ₁)/a and Ψ′(d)/(4πl_B) = −(1 − ϕ₂)/a. To proceed, we introduce the convenient notation B = a/(4πl_B) and switch from λ to the new quantity d_r = (ν/a) × e^λ. Both B and d_r carry a positive sign and have the dimension of a length. Note that B = (1 − ϕ)l_GC/2 is related to the Gouy–Chapman length³⁹ l_GC of a macroion surface with a degree of counterion association ϕ. The two electrostatic boundary conditions together with the counterion association/dissociation equilibria in eqn (5) then give rise to four algebraic equations

which define the four unknown quantities ϕ₁, ϕ₂, κ, and δ as function of the parameters B, d, d_r, and χ. Our goal is to obtain a complete understanding of the solutions eqn (8) offer. To this end, it is convenient to measure lengths in units of B. That is, we define the scaled quantities d̄ = d/(4B) = πl_Bd/a and d̄_r = d_r/(4B) = πl_Bνe^λ/a². In addition, we define the auxiliary dimensionless quantitiesallowing us to express eqn (8) in terms of d̄, d̄_r, C₁, C₂ throughwith

Inserting C₁ and C₂ into eqn (10) enables us to compute ϕ₁ and ϕ₂ based on solving only two algebraic equations, f₁(ϕ₁,ϕ₂) = 0 and f₂(ϕ₁,ϕ₂) = 0. Each equation defines a relationship, ϕ₁ = ϕ₁(ϕ₂), that can be viewed as a curve in a ϕ₁,ϕ₂-diagram. Because of the symmetry f₁(ϕ₁,ϕ₂) = f₂(ϕ₂,ϕ₁), the two curves are inverse to each other; that is, their graphs are mirror images across the line ϕ₁ = ϕ₂. Once ϕ₁ and ϕ₂ have been calculated, C₁ and C₂ follow from eqn (11), allowing us to determine from eqn (9)

This fully defines the solution Ψ(x) as function of d in eqn (7) and thus allows us to compute all other system properties, including the free energy F. An example is displayed in Fig. 2 for d̄_r = 0.8, χ = 3.0 and the two different values d̄ = 0.75 (left diagram) and d̄ = 0.95 (right diagram). We show solutions of eqn (10) (the red line for the upper and blue line for the lower equation), with white circles marking their intersections and thus the stationary points of the free energy. Also shown are contour lines that represent the free energy F × a/A, with darker colors indicating lower free energy. Clearly, the symmetric solution ϕ₁ = ϕ₂ represents a minimum of the free energy for d̄ = 0.75 whereas it has turned into a saddle point for d̄ = 0.95. The minima of the free energy are now located at asymmetric states ϕ₁ ≠ ϕ₂, indicating a symmetry breaking transition has taken place. Let us analyze that transition systematically.

Fig. 2 Degrees of counterion association, ϕ₁ and ϕ₂, depend on d̄_r, χ, d̄. The two diagrams show results for d̄_r = 0.8 and χ = 3.0, with d̄ = 0.75 (left diagram) and d̄ = 0.95 (right diagram). Each diagram displays solutions of eqn (10), the red line for the upper and blue line for the lower equation, with intersections locating the stationary points (marked by white circles). Also shown is the corresponding free energy F × a/A in eqn (2) by a contour plot, with darker colors for lower free energy values. A thin gray line marks the main diagonal where ϕ₁ = ϕ₂. Red and blue lines in each diagram are mirror images of each other with respect to that line. The stationary point on the line ϕ₁ = ϕ₂ corresponds to a minimum in the left diagram and to a saddle point in the right diagram.

3 Perturbation approach

For ϕ₁ = ϕ₂ both macroions carry the same effective charge. The solution Ψ(x) in eqn (7) is then mirror-symmetric with respect to the plane x = d/2, implying δ = 0. Eqn (10) always possess a symmetric solution ϕ_s = ϕ₁ = ϕ₂ but the corresponding stationary point may locate a minimum or a saddle point of the free energy. The transition from the former to the latter marks an instability that we can identify by a perturbation approach. We insert ϕ₁ = ϕ_s + Δϕ and ϕ₂ = ϕ_s − Δϕ into eqn (10)–(11). Sufficiently small Δϕ then gives rise to the two relations
which fully characterize the transition between symmetric states (where ϕ₁ = ϕ₂) corresponding to a local minimum and a saddle point. Eqn (13) represents a central result of this study, as the two relationships fully characterize the local stability of the symmetric state.

The upper diagram in Fig. 3 shows d̄_r and d̄ as function of ϕ_s for fixed χ = 3.0. For any given value d̄_r < 1 (which is marked by the blue circle), two distinct degrees of counterion association ϕ^A_s and ϕ^B_s with corresponding scaled distances d̄^A = d̄(ϕ^A_s) and d̄^B = d̄(ϕ^B_s) between the two macroions are predicted. Shown in the diagram is the specific choice d̄_r = 0.8, entailing ϕ^A_s = 0.3442 and ϕ^B_s = 0.6771 as well as d̄^A = 1.2569 and d̄^B = 0.8475. Recall that the two contour plots in Fig. 2 were calculated for χ = 3.0 and d̄_r = 0.8, with a minimum of the symmetric state for d̄ = 0.75 and a saddle point for d̄ = 0.95. Fig. 3 predicts the transition from local minimum to saddle point to take place at d̄ = d̄^B = 0.8475. Moreover, another transition back to a local minimum occurs at d̄ = d̄^A = 1.2569.

Fig. 3 The functions d̄_r(ϕ_s) (solid lines) and d̄(ϕ_s) (dashed lines) according to eqn (13). The upper diagram focuses on χ = 3.0, showing specific values of ϕ^A_s and ϕ^B_s at positions where d̄_r = 0.8. Corresponding scaled distances d̄^A = d̄(ϕ^A_s) and d̄^B = d̄(ϕ^B_s) are marked. The small circle locates the critical point ϕ_s = 1/2 and d̄_r = 1. The lower diagram includes five different values of χ as indicated in the legend. Critical points are marked by color-matching small circles. The gray solid line shows the progression of the critical points according to eqn (14).

The critical point, ϕ_s = 1/2 and d̄_r = 1, in the upper diagram of Fig. 3 is marked by a blue circle, which applies to χ = 3. No asymmetric equilibrium states exist for values d̄_r > 1. Because of d̄_r ∼ e^λ, it follows that sufficiently large counterion association energies λ prevent the formation of asymmetric equilibrium states ϕ₁ ≠ ϕ₂. Different values of χ lead to different critical points, which is shown in the lower diagram of Fig. 3. The gray line connecting the critical points for different χ follows the relationship

with the plus and minus signs applying to values of χ larger and smaller than , respectively. The corresponding marks the largest possible value for the onset of the asymmetry transition, and this is accompanied by d̄_r = 0.5580 and . No critical point can have a degree of counterion association larger than .

We also note the two critical values χ = 3 and χ = 2 for ϕ_s = 1/2. The former, χ = 3, implies d̄_r = 1 and d̄ = π/2. The latter, χ = 2, is the lowest value above which asymmetric states can exist. This occurs in the limit where both d̄_r → 0 and d̄ → 0. Here, electrostatic interactions are irrelevant, and since half the counterions are adsorbed, the critical point ϕ_s = 1/2 and χ = 2 is predicted entirely by the Bragg–Williams free energy in eqn (1).

4 Characterization of symmetry breaking

The perturbation approach above informs us that asymmetric states must exist if the symmetric state corresponds to a saddle point of the free energy. Yet, asymmetric states may even exist if the symmetric state falls into a local minimum of the free energy. In this case, a contour plot such as in Fig. 2 would have five stationary points, three local minima of the free energy separated by two saddle points. This can indeed be observed and has thus been studied by identifying all local free energy minima for choices of the three parameters, d̄_r, d̄, and χ. In Fig. 4 we show diagrams of ϕ₁ and ϕ₂ as function of d̄ for fixed d̄_r = 0.8 and four different values of χ.

Fig. 4 Degrees of counterion association, ϕ₁ and ϕ₂, as function of the scaled macroion-to-macroion distance d̄ for fixed d̄_r = 0.8 and four different values of χ as stated on top of diagrams A–D. States residing in local minima and at saddle points of the free energy are displayed by solid lines and dashed lines, respectively. Onsets of discontinuous and continuous transitions are marked by green and red bullets, respectively. Corresponding thin vertical lines assist in locating the positions d̄^A_b, d̄^A, d̄^B, and d̄^B_b. The four blue circles in diagram A mark the stationary points at d̄ = 0.75 and d̄ = 0.95 that correspond to the free energy contours in Fig. 2.

Fig. 4A, which is calculated for χ = 3.0, demonstrates the existence of a region with a metastable symmetric minimum and a stable asymmetric minimum between the values d̄^A_b = 0.7992 and d̄^A = 0.8475 and then again between d̄^B = 1.2569 and d̄^B_b = 1.2806 as the scaled macroion-to-macroion distance d̄ is increased. Two local minima of the free energy are always separated by a saddle point. The locations of all saddle points in Fig. 4 are shown as dashed lines. The four blue circles in Fig. 4A mark the stationary points for d̄ = 0.75 (one point) and d̄ = 0.95 (three points), which are consistent with the two diagrams in Fig. 2. We emphasize that the transitions from stable symmetric to stable asymmetric states in Fig. 4A at d̄ = d̄^A_b and d̄ = d̄^B_b are both discontinuous. They can be located by finding solutions of f₁(ϕ₁,ϕ₂) = 0 and f₂(ϕ₁,ϕ₂) = 0 such that the additional condition

ensures that two inverse curves intersect each other tangentially. The three equations specify d̄, ϕ₁, ϕ₂ for any given d̄_r and χ. For χ = 3.0, two such solutions exist, denoted by d̄^A_b, ϕ₁(d̄^A_b), ϕ₂(d̄^A_b), and d̄^B_b, ϕ₁(d̄^B_b), ϕ₂(d̄^B_b) and marked in Fig. 4A by green bullets.

The discontinuous nature of the symmetry breaking transition persists upon lowering χ until at χ = 2.9149 the locations d̄^B_b and d̄^B merge. Below χ = 2.9149, the point d̄^B_b does not exist, leaving only d̄^B to mark the onset of the transition, which is now continuous. For χ = 2.90, this is displayed in Fig. 4B. The red bullet at d̄^B marks the continuous transition. Note that the two values d̄^A_b = 0.8851 and d̄^A = 0.8923 remain distinct, implying one of the two transitions is still discontinuous. Below χ = 2.8562 both transitions become continuous and are predicted by eqn (13) (see also Fig. 3). The values of χ(d̄_r), where the transition from a discontinuous jump to a smooth symmetry breaking takes place, cannot be determined from eqn (13). We have computed them numerically. Fig. 4C was calculated for χ = 2.85, where both symmetry-breaking transitions occur in a continuous manner.

Decreasing χ even further renders the range of transition, from d̄^A to d̄^B, smaller until at χ = 2.7958 it ceases to exist. The critical point χ = 2.7958 is predicted by eqn (13). In the lower diagram of Fig. 3 it would be the critical point for the curve at χ = 2.7958, leading to d̄_r = 0.8 with ϕ_s = 0.5260 and d̄ = 1.0756. That curve falls between the displayed curves for χ = 2.7 and χ = 3.0. Fig. 4D shows the case χ = 2.797, barely above the critical point and thus leading to a narrowly confined region in which symmetry breaking takes place. No metastable regions are present in Fig. 4C and D.

5 Stability diagram

Recall that our model is completely defined by the three dimensionless parameters d̄_r = πl_Bνe^λ/a², d̄ = πl_Bd/a, and χ, expressing the counterion association strength λ, the macroion-to-macroion distance d, and the non-ideality parameter χ. Hence, our model's predictions can be cast into a three-dimensional stability diagram, which includes regions of local instability and metastability. This diagram is displayed in Fig. 5. Regions enclosed by solid lines have locally unstable symmetric states, implying symmetry breaking must occur. They are fully defined by eqn (13). Regions between color-matching solid and dashed lines are metastable: their symmetric states reside in a local but not in a global free energy minimum. The largest value of d̄_r = d̄_r(χ) that leads to an instability gives rise to a critical point, marked by circles in Fig. 5. Close to the critical point, no metastable regions exist. They appear only below values for d̄_r and d̄ marked by small bullets in Fig. 5. The locations of these small bullets explain when diagrams such as Fig. 4B exist, where one transition is continuous and the other discontinuous. The thin vertical blue line in Fig. 5 indicates the pathway, from d̄ = 0.5 to d̄ = 1.5 at fixed d̄_r = 0.8 and χ = 3.0, along which Fig. 4A was calculated. Here, stable and unstable symmetric states along this path are separated by two metastable regions, bounded by d̄^A_b, d̄^A and d̄^B, d̄^B_b. Using the stability diagram in Fig. 5, it is straightforward to predict if a given path intersects with zero, one, or two metastable regions.

Fig. 5 Stability diagram, indicating regions of spontaneous symmetry breaking as function of d̄_r and d̄, for different values of χ as specified in the legend. Regions enclosed by solid lines have locally unstable symmetric states. Regions between solid and dashed lines have metastable symmetric states. Open circles mark critical points, beyond which asymmetric states (with larger d̄_r) do not exist. No metastable states are present near critical points. Small bullets (for clarity shown only for χ = 2.7, 3.0, and 3.3) mark the beginning of the metastable regions. The thin blue vertical line shows the path along which the degrees of counterion association, ϕ₁ and ϕ₂, are shown in Fig. 4A.

Fig. 5 indicates that the region in which symmetry breaking occurs is always confined to a finite range, 0 < d < ∞, and does not extend to either contact (d = 0) or infinite separation (d → ∞). In the limit d → 0, the degree of counterion association approaches ϕ → 1 for any finite association energy −∞ < λ < ∞, driving the system close to full association (or, equivalently, zero dissociation). This proximity to ϕ = 1 suppresses the formation of asymmetric states due to the high entropy penalty. In the opposite limit d → ∞, symmetry breaking is absent because isolated, non-interacting macroions are assumed to be stable, an assumption that will be examined further in Section 6.4.

6 Discussion

Our theoretical model provides a description of spontaneous symmetry breaking of two interacting, chemically identical, planar macroions that are subject to charge regulation based on a non-ideal counterion association/dissociation isotherm according to the mean-field Bragg–Williams free energy in eqn (1). Each charge on the macroion is the result of a dissociation reaction that produces a mobile counterion. No salt is present, implying that only counterions reside between the two macroion surfaces. In the following, we discuss the applicability of our no-salt assumption, specify experimental conditions to which our model applies, investigate the role of electrostatic interactions for the spontaneous symmetry breaking, identify a maximal value of χ below which an isolated macroion remains stable, calculate the disjoining pressure between the interacting macroions, and compare our model with preceding work by Majee et al.³⁴

6.1 Applicability of the no-salt assumption

A crucial assumption of our model is the absence of salt. Salt, which is always present to some extent, gives rise to a characteristic length, the Debye length l_D, beyond which interactions are screened. The presence of salt can be neglected if interactions between macroions occur on length scales much smaller than l_D. To illustrate this point, we consider a single charged planar surface carrying a fixed negative surface charge density σ = −e/a in contact with a symmetric 1 : 1 electrolyte. The Debye length reflects the presence of salt with bulk concentration n₀. We solve the corresponding Poisson–Boltzmann equation l_D²Ψ″(x) = sinh Ψ(x) subject to the boundary conditions Ψ(0) = 0 and Ψ′(0) = 4πl_B/a = 1/B at the macroion surface x = 0. Inside the bulk of the electrolyte, for x → ∞, the potential will be constant and equal to the potential difference produced by the electric double layer. It is convenient to use the length B as unit length. The gray curves in Fig. 6 show Ψ as function of x/B for different values of the scaled Debye length l_D/B. As l_D/B grows, the potential Ψ(x) = 2 ln[1 + x/(2B)] is increasingly better approximated inside the region 0 ≤ x ≪ l_D. Yet, Ψ(x) = 2 ln[1 + x/(2B)] is exactly the potential produced by the Poisson–Boltzmann model in the absence of salt, solving Ψ″(x) = −(κ²/2) × e^−Ψ(x) with Ψ(0) = 0 and Ψ′(0) = 1/B. Hence, at distances much smaller than l_D away from a single charged surface, the use of the Poisson–Boltzmann model in the absence of added salt is justified.

Fig. 6 Log–log plot of the potential Ψ(x/B) as function of the scaled distance x/B (with B = a/(4πl_B)) away from a single planar macroion surface in contact with a symmetric 1 : 1 electrolyte of scaled Debye length l_D/B. The four gray curves labeled by the arrows in the diagram correspond to l_D = 10B, l_D = 30B, l_D = 100B, and l_D = 300B. For distances x ≪ l_D, the potential for the Poisson–Boltzmann model in the absence of salt, Ψ(x) = 2 ln[1 + x/(2B)] (shown as a black solid line) reproduces the potential for the Poisson–Boltzmann model in the presence of salt.

These considerations can be extended to the interaction between two planar surfaces with charge densities σ = −e/a. For separations d ≪ l_D, the potential Ψ(x) = Ψ(d/2) − (x − d/2)²/(Bd) with Ψ(d/2) = −ln[4l_D²/(Bd)] satisfies the Poisson–Boltzmann equation l_D²Ψ″(x) = sinh Ψ(x) together with the boundary conditions Ψ′(0) = −Ψ′(d) = 1/B at the two macroion surfaces located at x = 0 and x = d. Apart from the irrelevant constant Ψ(d/2), the same potential is obtained from eqn (7) for the salt-free model in the limit d ≪ l_D. This limit corresponds to the so-called ideal-gas regime,⁴⁰ which yields κ² = 4/(Bd).

6.2 Discussion of experimental relevance

The stability diagram in Fig. 5 indicates that if spontaneous symmetry breaking occurs, the dimensionless separation d̄ ≈ 1 will be of order unity. As discussed below in Section 6.4, for an isolated macroion becomes unstable with respect to internal phase separation, resulting in spatially heterogeneous counterion association along the macroion surface. Only for values of χ slightly below this threshold can spontaneous symmetry breaking occur at large, but not infinitely large, separations 1 ≪ d̄ < ∞. For example, at χ = 3.7, the largest separation at which symmetry breaking can occur is d̄ = 8. The characteristic value d̄ = πl_Bd/a ≈ 1 implies, roughly, a physical separation d ∼ a/l_B for the onset of symmetry breaking. Consequently, we expect our present salt-free theoretical predictions to become experimentally testable for salt concentrations with corresponding Debye length l_D ≫ a/l_B. Typical charge densities of macroions have magnitudes |e/a| with a varying by two orders of magnitude, from a ≈ 100 nm² for silica particles near neutral pH, weakly functionalized latex particles, weakly charged lipid membranes, and microgels under weak ionization, to a ≈ 1 nm² for DNA, cytoskeletal filaments, viral capsids, and highly charged lipid membranes.³⁹ Assuming a Bjerrum length of order l_B = 1 nm we arrive at salt concentrations significantly smaller than 10 µM for weakly charged macroions and significantly smaller than 100 mM for highly charged macroions. We conclude that for a wide range of macroions experimentally accessible salt concentrations exist that can test our present model. For example, recent experiments by Wang et al.⁴ employed charged silica particles and aminated silica particles at salt concentrations starting from 5 µM. Under these conditions our model would be applicable if a ≲ 100 nm². Experimentally supported charge densities of silica particles, including functionalized silica particles, fall well into this range.^41–43

Let us connect the dimensionless parameters d̄_r and d̄ to values for λ and d using a characteristic experimental scenario. Recall first our definitions d_r = ν/a × e^λ, d̄_r = d_r/(4B), d̄ = d/(4B), and B = a/(4πl_B), which yield d̄_r = πl_Bνe^λ/a² and d = d̄a/(πl_B). Let us choose a = 10 nm² for the area per dissociable group on the macroion surface, l_B = 1 nm for the Bjerrum length, χ = 3 for the non-ideality parameter, and λ = 0.24 for the adsorption strength. The volume ν sets a reference for measuring λ. Specifically λ = 0 and χ = 0 yield a 50%-degree of association/dissociation at a distance d = ν/(2a) in the absence of electrostatic interactions. Choosing ν = 20 nm³ sets that distance to be d = 1 nm. With that we obtain d̄_r = 0.8. From Fig. 4A we recall the region 0.8 < d̄ < 1.28 for spontaneous symmetry breaking. This then translates into a range 2.5 nm < d < 4.1 nm of macroion-to-macroion distances at which our model predicts symmetry breaking. Applicability of this example would require a Debye screening length of l_D ≫ a/l_B = 10 nm and thus a salt concentration of much less than 1 mM.

We finally note that while our model refers to the planar macroion geometry, it can also be applied to spherical colloids using the Derjaguin approximation,^39,44 given the radius of the colloid R ≫ l_D is much larger than the Debye screening length l_D.

6.3 Role of electrostatic interactions

In the following, we aim to understand the role electrostatic interactions play for the spontaneous symmetry breaking. Would “switching off” the Coulomb interaction between all involved charges suppress or enhance the occurrence of symmetry breaking? This can be investigated by ignoring the electrostatic energy term Ψ′²/(8πl_B) in eqn (2). The counterions, which now do no longer carry an electric charge, distribute uniformly between the two uncharged plates. Changing the plate-to-plate distance d still affects the degree of adsorption ϕ₁ and ϕ₂. In fact, the two equations and specify ϕ₁ and ϕ₂ for any stationary state of the free energy. The absence of electrostatics implies l_B → 0. The length B = a/(4πl_B) therefore diverges and thus can no longer be used to scale other quantities. Instead, the ratio d/d_r emerges as a convenient dimensionless quantity to characterize the onset of spontaneous symmetry breaking. A perturbation calculation analogous to that in Section 4 shows that in between the two values d/d_r = 2(χ − 1 + s)/(χ − s)e^−s with , the symmetric state is unstable. The two values of d/d_r are plotted in the inset of Fig. 7. To directly compare the instability regions in the absence and presence of electrostatic interactions, we re-plot the regions inside which symmetric states are locally unstable from Fig. 5 as function of d/d_r instead of d̄. This is shown in the main diagram of Fig. 7. Recall the horizontal axis shows d̄_r = d_r/(4B) with B = a/(4πl_B), implying d̄_r = πl_Bd_r/a. The strength of the Coulomb potential is proportional to the Bjerrum length l_B, implying d̄_r → 0 in the absence of electrostatic interactions. This corresponds to the instability region on the left axis of Fig. 7. For example, the red curve, for χ = 2.1, implies d/d_r = 1.200 and d/d_r = 0.7933 in the limit d̄_r = 0, where electrostatic interactions are absent. These two points appear in the inset of Fig. 7 at the same positions and are marked by the two red bullets. Instability points for other values of χ are marked by color-matching bullets. Fig. 7 reveals the behavior implied by “switching on” electrostatic interactions corresponds to moving horizontally through the stability diagram, starting at d̄_r = 0. This assumes all other parameters (χ and λ) remain fixed. If a system is unstable at d̄_r = 0, it will always move into a stable region. If it is not stable at d̄_r = 0, it may pass through a region of instability, but eventually end up in a stable region. In this sense, electrostatic interactions to not tend to induce spontaneous symmetry breaking. Instead, they stabilize the symmetric state, ϕ₁ = ϕ₂.

Fig. 7 Regions inside which symmetric states are locally unstable, re-plotted from Fig. 4 as function of d/d_r and d̄_r = πl_Bd_r/a for different values of χ as specified in the legend. The absence of electrostatic interactions (l_B → 0) corresponds to the line d̄_r = 0. For an isolated planar macroion becomes unstable, with the instability occurring first at , indicated by the dashed vertical black line. The inset shows d/d_r as function of χ in the absence of electrostatic interactions. The colored bullets match the curves shown in the main diagram at position d̄_r = 0.

We exemplify the predictions of Fig. 7 by a comparison of the degrees of counterion association, ϕ₁ and ϕ₂, in the absence and presence of electrostatic interactions for two specific cases shown in Fig. 8, with χ = 2.4 (left diagram) and χ = 3.0 (right diagram). The black lines show ϕ₁ and ϕ₂ in the presence of electrostatics with d̄_r = 0.8 as in Fig. 4. The gray lines refer to the absence of electrostatic interactions, where B → ∞. Note that Fig. 7 and 8 make consistent predictions: for χ = 2.4, the presence of electrostatic interactions inhibits symmetry breaking, whereas for χ = 3.0 the region of symmetry breaking is shifted to somewhat larger d/d_r.

Fig. 8 Degrees of counterion association, ϕ₁ and ϕ₂, as function of the scaled macroion-to-macroion distance d/d_r for fixed d̄_r = 0.8 as well as χ = 2.4 (left diagram) and χ = 3.0 (right diagram). Black and gray lines refer to the presence and absence of electrostatics, respectively. Solid lines correspond to states that are locally stable, dashed lines characterize stationary points of the free energy that are unstable. Green and red bullets mark discontinuous and continuous transitions between symmetric (ϕ₁ = ϕ₂) and asymmetric (ϕ₁ ≠ ϕ₂) degrees of counterion association. The black lines in the right diagram, which are computed for d̄_r = 0.8 and χ = 3.0 in the presence of electrostatic interactions, represent the same data as in Fig. 4A.

6.4 Stability of an isolated macroion

Throughout this work we have considered values of χ in the range from 2 < χ ≤ 3.3. If χ is significantly larger, even an isolated macroion will undergo an instability and phase separate into regions with different degrees of counterion association. Let us find the largest χ for which a single macroion is still stable. Repeating the theoretical approach of the present work for an isolated macroion with a degree of counterion association ϕ leads to the equationwhere we recall d̄_r = πl_Bνe^λ/a². Eqn (16) predicts a critical point , , and . At the critical point, two solutions for ϕ will emerge from eqn (16) upon any small increase of χ. We conclude that only for χ < 3.732 does an isolated macroion remain stable against an internal phase transition that would result in spatially inhomogeneous values of ϕ. To place this prediction into the context of the present work, we have added a curve (the gray line) for χ = 3.732 to the stability diagram in Fig. 7. This is the smallest χ for which the instability region grows to infinitely large values of d/d_r, and the divergence occurs at d̄_r = 1.952 (marked by the dashed vertical black line in Fig. 7), exactly as predicted by eqn (16). At this point, an isolated macroion surface becomes unstable with respect to the separation into domains with different degrees of counterion association.

Yet another line of context can be established by considering a single planar macroion surface of fixed degree of counterion association ϕ in the presence of salt with corresponding Debye length l_D, as in Fig. 6. Such a surface with fixed charge density σ = −(1 − ϕ)e/a gives rise to a spinodal line^45,46

The critical point predicted by eqn (17) increases from χ = 2 and ϕ = 1/2 in the limit of large salt content (l_D → 0) to and in the limit of no salt content (l_D → ∞), exactly as the critical point that emerges from eqn (16) in the presence of charge regulation. The critical point χ = 3.732 for a non-interacting planar surface in the absence of salt was first noted by Gelbart and Bruinsma.⁴⁷ In summary, the range of χ values considered in this work is well justified, since only for χ < 3.732 is an isolated, non-interacting planar macroion surface guaranteed to remain stable with respect to phase separation.

We point out that interacting macroions may adopt non-homogeneous degrees of counterion association even for 2 < χ < 3.732. The characterization of these non-homogeneous degrees is outside the scope of the present work. In fact, because the electrostatic contribution to the line tension of two coexisting domains of different surface charge density on a macroion is negative,⁴⁸ the consequence may be the formation of interlocked modulated phases with periodically changing degrees of dissociation along each macroion. Testing this hypothesis may be the subject of a future study.

6.5 Disjoining pressure

The disjoining pressure P = n(x) − Ψ′(x)²/(8πl_B) that acts between the two planar surfaces consists of two contributions, the osmotic pressure of the mobile counterions and the Maxwell stress. Indeed, Beltrami's identity⁴⁹ ensures P emerges as an invariant, independent of x, upon integrating the Poisson–Boltzmann equation Ψ″ = −(κ²/2) × e^−Ψ once. Using n(x) = e^−Ψ(x) × κ²/(8πl_B) together with Ψ(x) in eqn (7) leads to P = κ²/(8πl_B). For the scaled (dimensionless) pressure P̄ = PaB we thus obtainwhere C₁ and C₂ follow from eqn (11) after identifying ϕ₁ and ϕ₂ from solving f₁(ϕ₁,ϕ₂) = 0 and f₂(ϕ₁,ϕ₂) = 0 according to eqn (10). Note that the pressure is always positive, indicating the two macroion surfaces always repel each other.

Let us consider the example d̄_r = 0.8 and χ = 2.85, discussed in Fig. 4C. Recall there are continuous transitions from the symmetric to the asymmetric state at d̄ = d̄^A and then back from the asymmetric to the symmetric state at d̄ = d̄^B as d̄ increases. For d̄ = d̄^A and d̄ = d̄^B we obtain P̄ = 0.1850 and P̄ = 0.1889, respectively. That is, the pressure increases when increasing the distance between the macroions across the region where ϕ₁ ≠ ϕ₂. The full behavior of the scaled pressure as function of the scaled distance P̄ = P̄(d̄) is shown by the black line (corresponding to χ = 2.85) in Fig. 9 together with red and blue lines for χ = 2.9 and χ = 3.0. Solid and dashed lines correspond to locally stable and unstable states, respectively. The inset shows a magnified view of the region where spontaneous symmetry breaking takes place. Recall from Fig. 4A and B that for χ = 2.9 only one of the two transitions is continuous, the other discontinuous, and for χ = 3.0 both transitions are discontinuous. When the transition is discontinuous, the pressure jumps. Corresponding d̄-values are marked in the inset of Fig. 9 by open circles. In the d̄-region where the macroions have asymmetric degrees of dissociation, the pressure is a non-monotonic function of the distance d. This non-monotonic and in some cases even discontinuous behavior of the pressure is a consequence of the two macroions each being assumed to maintain a spatially uniform degree of counterion association when they interact.

Fig. 9 Scaled disjoining pressure P̄ = PaB as function of the scaled distance d̄ = d/(4B) = d × πl_B/a for d̄_r = 0.8 and χ = 3.0 (blue), χ = 2.9 (red), and χ = 2.85 (black). The three curves correspond to the degrees of counterion association shown in Fig. 4A–C. Solid and dashed lines correspond to locally stable and unstable states, respectively. The inset shows a magnified view of the region where spontaneous symmetry breaking takes place. Dashed lines between bullets mark unstable symmetric states. Open circles locate transitions to the asymmetric state that are discontinuous.

To understand the origin for the non-monotonic pressure, it is important to distinguish between spontaneous symmetry breaking and phase separation. In the former case, each of the two macroions may adopt a different and yet uniform degree of counterion association over the entire surface. In the latter case, regions with distinct degrees of counterion association grow or shrink, similar to coexisting phases in a phase separation. In the absence of electrostatics, it is straightforward to compare both cases. Recall our discussion of “switching off” electrostatic interactions in Section 6.3, where we have presented in Fig. 8 the degrees of counterion association ϕ₁ and ϕ₂ for the two values χ = 2.4 and χ = 3.0. The gray curves in Fig. 8 refer to the absence of electrostatic interactions. Fig. 10 shows the corresponding (conveniently scaled) disjoining pressure P̄ = Pad_r = (2 − ϕ₁ − ϕ₂)d_r/d in the absence of electrostatic interactions for the same parameters as in Fig. 8. Here again, dashed lines mark locally unstable states, corresponding to saddle points of the free energy. Gray solid lines refer to stable asymmetric states. The qualitative behavior of P̄ is similar in the absence (Fig. 10) and presence (Fig. 8) of electrostatic interactions, changing non-monotonically in the d-region where spontaneous symmetry breaking occurs. For discontinuous symmetry breaking (at positions marked by open circles in Fig. 10), the pressure jumps discontinuously. If we replace the occurrence of symmetry breaking by the ability to phase separate, two phases with different degrees of “counterion” (the counterions are now uncharged) association will coexist on the two surfaces, and the pressure will stay constant when d is changed. The corresponding pressure is determined by the well-known Maxwell construction and indicated in both diagrams of Fig. 10 by a brown horizontal line with two brown triangles locating the d-region where phase separation takes place. Thus, the pressure does not exhibit non-monotonic behavior during a phase transition if electrostatic interactions are absent.

Fig. 10 Scaled disjoining pressure P̄ = Pad_r = (2 − ϕ₁ − ϕ₂)d_r/d as function of the scaled distance d/d_r for χ = 2.4 (left diagram), and χ = 3.0 (right diagram). No electrostatic interactions are present (l_B → 0). Recall that and specify ϕ₁ and ϕ₂ for any stationary state of the free energy. Dashed lines correspond to locally unstable states. Black solid lines mark locally stable symmetric states, gray solid lines stable asymmetric states. Black bullets separate locally unstable and stable symmetric states. Gray open circles mark discontinuous transitions of the pressure, occurring in conjunction with discontinuous symmetry breaking; see Fig. 8. Brown horizontal lines terminating in two brown triangles locate a constant pressure according to the Maxwell construction, as is the case when two phases with different degrees of counterion association ϕ₁ and ϕ₂ coexist.

Our consideration of electrostatic interactions being absent makes the distinction between spontaneous symmetry breaking and phase separation straightforward. When electrostatic interactions are present, the situation is less clear. First, any macroscopic phases that form interact electrostatically with each other, preventing the formation of a d-region with constant pressure and thus the application of the Maxwell construction. Second, macroscopic phases may not form in the first place because the electrostatic contribution to the line tension between coexisting phases is negative.⁴⁸ This may lead to modulated phase structures that are characterized by a certain length scale.⁵⁰ Modulated phases on different macroions may then correlate. Third, a term penalizing gradients in the degree of counterion association may be present in the free energy. If that penalty is sufficiently high, it will lead to different and yet homogeneous degrees of counterion association on each macroion. In this case, the present model will apply. To sum up, to systematically account for heterogeneous degrees of counterion association along each macroion would be an interesting subject but is beyond the scope of the present study.

6.6 Comparison with model proposed by Majee et al.³⁴

Majee et al.³⁴ previously studied two interacting planar macroions with a charge regulation model governed by the Frumkin–Fowler–Guggenheim isotherm. They employ the notation χ_M = 2χ and λ_M = λ − χ in a free energy expression that is formally identical to eqn (2), yet with two differences. First, salt is present at a fixed chemical potential and corresponding Debye length l_D. Second, the number density of surface charges is (1/2 − ϕ)e/a, in contrast to (1 − ϕ)e/a used in the present work. This choice introduces additional adsorption sites, allowing the surface charge density to change sign and enabling attractive interactions between macroions with asymmetric degrees of dissociation ϕ₁ ≠ ϕ₂. Owing to the absence of simple analytic solutions of the Poisson–Boltzmann equation for interacting surfaces in the presence of salt and to the five-dimensional parameter space of the model, their study is restricted to selected examples that demonstrate symmetry breaking, without a systematic characterization of metastable regions, critical behavior, or the role of electrostatic interactions for the spontaneous symmetry breaking.

The consideration of the stability of an isolated planar macroion clearly illustrates the differences of the model proposed by Majee et al.³⁴ to our present work. Using our present notation, their free energy for a single planar isolated macroion reads

with p = (1/2 − ϕ) × l_D/(2B) and . Vanishing second derivative of aF/A with respect to ϕ yields a spinodal linethat is different from eqn (17) and implies a critical point χ = 2 at ϕ = 0.5 in the large salt limit (l_D → 0) and χ = 6.458 at ϕ = 0.203 in the no-salt limit (l_D → ∞). Majee et al.³⁴ have explored non-ideality parameters up to χ = 20 (that is, χ_M = 40 for their definition of the non-ideality parameter), suggesting the isolated macroions would not necessarily be stable with respect to separating into different phases with distinct values of ϕ. For our present model we ensure isolated macroions are stable by choosing χ < 3.732, as discussed in Section 6.4.

7 Conclusions

In this work we have analyzed a theoretical model for two interacting charge-regulated macroions, represented by two planar surfaces of the same bare charge density. The corresponding counterions can either be dissociated, thus residing in solution between the two macroions, or associated with the macroion surfaces through a Frumkin–Fowler–Guggenheim isotherm. This isotherm is structurally identical to the mean-field Bragg–Williams free energy of a lattice gas with interactions. That is, macroion charges that emerge from counterion dissociation are assumed to be subject of non-electrostatic interactions that can be attractive. No added salt is present.

We have shown that a spontaneous symmetry breaking can occur, leaving the two macroion surfaces with unequal degrees of dissociation and thus unequal effective charge densities. Analytic results show that, if present, symmetry breaking occurs only at intermediate distances between the two macroions, but not for sufficiently small or large distances. A three-dimensional stability diagram describes all parameter combinations where symmetry breaking is observed, either in the presence of metastable symmetric states or in their absence. Regions with metastable states do not exist near critical points.

We have discussed that our assumption of modeling the interaction between planar macroions in the complete absence of salt will become applicable if the Debye screening length l_D ≫ a/l_B is much larger than the typical cross-sectional area per charge a on the macroion divided by the Bjerrum length l_B. For example, a = 10 nm² and l_B = 1 nm would require a salt concentration significantly smaller than 1 mM. Many experimental systems carried out at low salt content will satisfy the condition l_D ≫ a/l_B. Another important assumption of our work is the uniformity of the degree of association/dissociation on each macroion. To ensure this is always the case at least for an isolated macroion, we have imposed the upper limit χ < 3.732 on our non-ideality parameter χ. The occurrence and structure of non-uniformly distributed degrees of counterion association for interacting macroions with χ < 3.732 remains an open question. Finally, we have studied the role of electrostatic interactions for the symmetry breaking. Interestingly, the hypothetical process of switching on the electrostatic interaction (that is, increasing the Bjerrum length l_B) reveals a tendency to stabilize the symmetric counterion association state ϕ₁ = ϕ₂. It is thus appropriate to state that, if present, symmetry breaking occurs not because of electrostatic interactions but in spite of them.

The present work offers a comprehensive understanding of spontaneous symmetry breaking for planar macroions in the absence of salt, based on the Frumkin–Fowler–Guggenheim isotherm. It complements preceding case studies in the presence of salt and for non-planar geometries,^34–37 where analytic results are not available.

Conflicts of interest

There are no conflicts to declare.

Data availability

All data that support the findings of this study are included within the article.

Acknowledgements

G. V. B. thanks Fondecyt (grant no. 11240064) for support.

Notes and references

1. G. Gonella, E. H. G. Backus, Y. Nagata, D. J. Bonthuis, P. Loche, A. Schlaich, R. R. Netz, A. Kühnle, I. T. McCrum and M. T. M. Koper, et al., Nat. Rev. Chem., 2021, 5, 466–485.
2. I. Popa, P. Sinha, M. Finessi, P. Maroni, G. Papastavrou and M. Borkovec, Phys. Rev. Lett., 2010, 104, 228301.
3. A. Kubincová, P. H. Hünenberger and M. Krishnan, J. Chem. Phys., 2020, 152, 104713.
4. S. Wang, R. Walker-Gibbons, B. Watkins, M. Flynn and M. Krishnan, Nat. Nanotechnol., 2024, 19, 485–493.
5. S. Wang, R. Walker-Gibbons, B. Watkins, B. Lin and M. Krishnan, Nat. Commun., 2025, 16, 1–16.
6. J. C. Neu, Phys. Rev. Lett., 1999, 82, 1072.
7. A. Behjatian, R. Walker-Gibbons, A. A. Schekochihin and M. Krishnan, Langmuir, 2022, 38, 786–800.
8. J. E. Hallett, D. A. J. Gillespie, R. M. Richardson and P. Bartlett, Soft Matter, 2018, 14, 331–343.
9. Y. Avni, D. Andelman and R. Podgornik, Curr. Opin. Electrochem., 2019, 13, 70–77.
10. T. Markovich, D. Andelman and R. Podgornik, Europhys. Lett., 2018, 120, 26001.
11. A. Behjatian, R. Blossey and M. Krishnan, J. Chem. Phys., 2025, 162, 064901.
12. S. H. Behrens and M. Borkovec, Phys. Rev. E:Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 60, 7040.
13. R. Pericet-Camara, G. Papastavrou, S. H. Behrens and M. Borkovec, J. Phys. Chem. B, 2004, 108, 19467–19475.
14. A. Bakhshandeh, D. Frydel, A. Diehl and Y. Levin, Phys. Rev. Lett., 2019, 123, 208004.
15. D. Frydel, J. Chem. Phys., 2019, 150, 194901.
16. A. Bakhshandeh, D. Frydel and Y. Levin, Phys. Chem. Chem. Phys., 2020, 22, 24712–24728.
17. T. Curk and E. Luijten, Phys. Rev. Lett., 2021, 126, 138003.
18. D. Izzo, A. Bakhshandeh and Y. Levin, J. Chem. Phys., 2025, 163, 064703.
19. G. Trefalt, S. H. Behrens and M. Borkovec, Langmuir, 2016, 32, 380–400.
20. K. H. Chu and B. C. Tan, Colloid Interface Sci. Commun., 2021, 45, 100519.
21. H. Tamura and R. Furuichi, J. Colloid Interface Sci., 1997, 195, 241–249.
22. C. A. Başar, J. Hazard. Mater., 2006, 135, 232–241.
23. K. Y. Foo and B. H. Hameed, Chem. Eng. J., 2010, 156, 2–10.
24. W. Rudzinski and D. H. Everett, Adsorption of gases on heterogeneous surfaces, Academic Press, 2012.
25. V. Kalinin and C. Radke, Colloids Surf., A, 1996, 114, 337–350.
26. V. Kolev, K. Danov, P. Kralchevsky, G. Broze and A. Mehreteab, Langmuir, 2002, 18, 9106–9109.
27. M. Volkova-Gugeshashvili, A. Volkov and V. Markin, Russ. J. Electrochem., 2006, 42, 1073–1078.
28. N. S. Mousavi, A. Javadi, V. I. Kovalchuk, E. V. Aksenenko, G. Loglio, D. Vollhardt, E. Schneck and R. Miller, Z. Phys. Chem., 2025, 239, 1705–1724.
29. A. Fomkin, A. Pribylov, I. Men'shchikov, A. Shkolin, O. Aksyutin, A. Ishkov, K. Romanov and E. Khozina, Reactions, 2021, 2, 209–226.
30. M. Levi and D. Aurbach, Electrochim. Acta, 1999, 45, 167–185.
31. S. A. Onaizi, M. Nasser and N. M. Al-Lagtah, J. Surfactants Deterg., 2016, 19, 645–652.
32. H. Tamura, N. Katayama and R. Furuichi, Environ. Sci. Technol., 1996, 30, 1198–1204.
33. S. J. Brugman, B. L. Werkhoven, E. R. Townsend, P. Accordini, R. van Roij and E. Vlieg, J. Colloid Interface Sci., 2020, 559, 291–303.
34. A. Majee, M. Bier and R. Podgornik, Soft Matter, 2018, 14, 985–991.
35. A. Majee, M. Bier, R. Blossey and R. Podgornik, Phys. Rev. Res., 2020, 2, 043417.
36. P. Khunpetch, A. Majee and R. Podgornik, Soft Matter, 2022, 18, 2597–2610.
37. H. Ruixuan, A. Majee, J. Dobnikar and R. Podgornik, Eur. Phys. J. E:Soft Matter Biol. Phys., 2023, 46, 115.
38. H. T. Davis, Statistical Mechanics of Phases, Interfaces, and Thin Films, Wiley-VCH, New York, 1996.
39. J. N. Israelachvili, Intermolecular and Surface Forces, Academic Press, London, 3rd edn, 2011.
40. D. Andelman, Handbook of biological physics, Elsevier, 1995, vol. 1, pp. 603–642.
41. C. J. Van Oss, M. K. Chaudhury and R. J. Good, Chem. Rev., 1988, 88, 927–941.
42. S. H. Behrens and D. G. Grier, J. Chem. Phys., 2001, 115, 6716–6721.
43. H. E. Bergna and W. O. Roberts, Colloidal silica: fundamentals and applications, CRC Press, 2005.
44. S. Oversteegen and H. Lekkerkerker, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2003, 68, 021404.
45. S. May, D. Harries and A. Ben-Shaul, Phys. Rev. Lett., 2002, 89, 268102.
46. C. L. Baciu and S. May, J. Phys.:Condens. Matter, 2004, 16, S2455.
47. W. M. Gelbart and R. Bruinsma, Phys. Rev. E:Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 55, 831.
48. G. V. Bossa, M. A. Brown, K. Bohinc and S. May, Int. J. Adv. Eng. Sci. Appl. Math., 2016, 8, 101–110.
49. G. B. Arfken, H. J. Weber and F. E. Harris, Mathematical methods for physicists: a comprehensive guide, Academic Press, 2011.
50. D. Andelman and R. E. Rosensweig, J. Phys. Chem. B, 2009, 113, 3785–3798.

---