# Charge symmetry broken complex coacervation

Arghya Majee, Markus Bier, Ralf Blossey, Rudolf Podgornik

► To cite this version:

Arghya Majee, Markus Bier, Ralf Blossey, Rudolf Podgornik. Charge symmetry broken complex coacervation. Physical Review Research, 2020, Physical Review Research, 2 (4), pp.043417. 10.1103/phys-revresearch.2.043417. hal-03179029

HAL Id: hal-03179029

https://hal.univ-lille.fr/hal-03179029

Submitted on 24 Mar 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Charge symmetry broken complex coacervation

Arghya Majee ●,¹,² Markus Bier ●,¹,²,³ Ralf Blossey ●,⁴ and Rudolf Podgornik ●,⁵,⁶

¹Max Planck Institute for Intelligent Systems, Stuttgart, Germany

²IV. Institute for Theoretical Physics, University of Stuttgart, Germany

³Fakultät Angewandte Natur- und Geisteswissenschaften, Hochschule für Angewandte Wissenschaften Würzburg-Schweinfurt, Germany

⁴Université de Lille, CNRS, UMR8576 Unité de Glycobiologie Structurale et Fonctionnelle (UGSF), 59000 Lille, France

⁵School of Physical Sciences and Kavli Institute for Theoretical Sciences, University of Chinese Academy of Sciences, Beijing, China

⁶CAS Key Laboratory of Soft Matter Physics, Institute of Physics, Chinese Academy of Sciences, Beijing, China

<figure><img src="image_2.png" /></figure>

(Received 1 October 2020; accepted 30 November 2020; published 24 December 2020)

Liquid-liquid phase separation has emerged as one of the important paradigms in the chemical physics as well as biophysics of charged macromolecular systems. We elucidate an equilibrium phase separation mechanism based on charge regulation, i.e., protonation-deprotonation equilibria controlled by pH, in an idealized macroion system which can serve as a proxy for simple coacervation. First, a low-density density functional calculation reveals the dominance of two-particle configurations coupled by ion adsorption on neighboring macroions. Then a binary cell model, solved on the Debye-Hückel as well as the full nonlinear Poisson-Boltzmann level, unveils the charge symmetry breaking as inducing the phase separation between low- and high-density phases as a function of pH. These results can be identified as a *charge symmetry broken complex coacervation* between chemically identical macroions.

DOI: 10.1103/PhysRevResearch.2.043417

I. INTRODUCTION

The importance of *complex coacervation* in polymers, colloids, and particularly proteins that exhibit an associative liquid-liquid phase separation (LLPS), driven by electrostatic interactions between oppositely charged macroions, has been recognized for about a century [1,2], though its fundamental role in compartmentalization and intracellular phase transitions in biological systems has been identified only recently [3]. The electrostatically driven attractions, as already hypothesized in the early Overbeek-Voorn theory [4], and later developed within more sophisticated theoretical frameworks [5], have been shown to result in LLPS, thus being recognized as the defining feature of complex coacervation [6]. On the other hand, for like-charged macroions with monovalent counterions, it is the variation of the solvent conditions, such as temperature, pH, and ionic strength [7], that modifies the electrostatic repulsion which would otherwise prevent coacervation, except when countered by the like-charge attraction mediated by multivalent counterions [8]. Studies of adhesive proteins [9] as well as several proteins involved in some protein aggregation diseases (e.g., Alzheimer's disease and amyotrophic lateral sclerosis) made it clear, however, that *simple coacervation* involving only sim-

ilarly charged macroions can also lead to LLPS, presumably because of short-range specific interactions of nonelectrostatic nature [10].

The proper understanding of the mechanisms of oppositely charged (complex) and similarly charged (simple) coacervations, interesting in the context of functional biomimetic and adhesive materials of the chemical, pharmaceutical, textile, and food industries [11], and particularly relevant in the biophysical milieu, where different facets of protein chemistry [12] can lead to coexisting liquidlike states, has been claimed to be one of the most important problems in the physical chemistry of the cytoplasm [13].

While experimentally well documented, the dependence of the associative LLPS on the bathing environment conditions, such as the solution pH [14], has lacked a comprehensive theoretical elucidation based on relevant microscopic models. That these effects are particularly important in protein solutions [12] is clear from the fact that the protein charge is not fixed, but is a result of the proton-mediated dissociation of amino-acid (AA) groups at the solvent accessible surface [15], whose chemical equilibrium then depends on the bathing environment parameters such as the solution pH [16]. The physical basis of the protein charging is consequently understood as the *charge-regulation (CR) mechanism*, i.e., an association and dissociation process that couples the local electrostatic potential with the local charge, leading to a self-consistent partitioning of the protein charge states with pronounced effects also on the properties of other macroions such as weak polyelectrolyte solution and gel conformational as well as charge properties (see Ref. [17] for details).

Theoretical analyses of the CR effects in the formation of macroion condensates, that depend explicitly on the solution

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Open access publication funded by the Max Planck Society.</footnote>
pH, have been scant. A simple cell model approach was used to analyze the CR macroions in solution [18], together with their effective charge [19], and the corresponding phase behavior [20]. A thermodynamic minimal model analysis was proposed to study LLPS in a fixed pH ensemble based on a set of reactions describing the protonation-deprotonation reactions of the solution macroions, conducive to multiple charge states [21]. The equilibrium charge state and critical behavior of CR macroions was studied based on a collective description of a solution composed of CR macroions and simple salt ions in the bulk [22]. Within the mean-field approximation it was found that above a critical concentration of salt, a nontrivial distribution of coexisting charge states leads to a liquid-liquid phase separation, similar to the behavior of micellar solutions close to the critical micelle concentration [23]. Extreme CR, implying a constant surface potential, was also invoked in binary suspensions of charged colloids [24], possibly leading to charge-alternating linear strings.

In what follows, we will present a detailed analysis not only of the liquid-liquid transition in CR macroion systems, but also the corresponding spatial charge distribution that is at its origin. The central idea, as depicted in Fig. 1, stems from the striking observation [25] that a pair of chemically identical interacting charge-regulated planar macroions are not necessarily equally charged and that the electric field at the midplane of the setup does not necessarily vanish. In order to provide a firm basis to the intuitive expectations on the charge symmetry breaking transition for a spherical macroion system, we present arguments based on a density functional theory (DFT) as well as on a binary cell model (BCM). Moreover, we show that the LLPS is based on a symmetry breaking transition of the macroion charge distribution, characterized by a spatially alternating sign of the macroion charge. In that respect this CR system driving a complex coacervation behaves not unlike the alternating multilayer structure of the electrical double layer in ionic liquids [26], except that here the charge alternation is driven by CR and not by the presence of different ion species. We identify this spatial charge layering, stemming from a symmetry broken charge distribution and leading to phase behavior that exhibits features of complex coacervation phenomenology, as charge symmetry broken complex coacervation between chemically identical macroions.

## II. CHARGE-REGULATION MODEL

As shown in Fig. 1(a), consider spherical macroions (e.g., proteins, polyelectrolytes, colloids, etc.) of radius $R_0$, whose surface charge is regulated following a mechanism identical to the charge-regulation model introduced in [25,27,28], that are suspended in a univalent salt solution. In short, each macroion surface contains a fixed number of negative charges and twice as many neutral sites where adsorption or desorption of protons can take place.

The surfaces are charge regulated through this adsorption and desorption, and the fraction $\eta$ of filled sites on a surface is a degree of freedom within our model. By construction, $\eta \in [0, 1]$. If the area per site is $a^2$, then the charge density is given by $\sigma = \frac{e}{a^2}(\eta - \frac{1}{2})$ with $e > 0$ being the elementary charge, so that $-\frac{1}{2}\frac{e}{a^2} \le \sigma \le \frac{1}{2}\frac{e}{a^2}$. The surface number density $1/a^2$

<figure><img src="image_3.png" /><figcaption>FIG. 1. Macromolecular solution (a) and the magnified view (b) of a small portion of it. Within a binary cell model, each macroion (indicated by red circles) of radius $R_0$ is surrounded by a cell of radius $R$ and the interaction between a pair such cell-surrounded spheres is considered. The electric field $E_R$ at the cell boundary is assumed to be uniform and as it is the case for interacting planar surfaces, the value of $E_R$ depends on the charge states of the neighboring macroions forming the pair. The two cells of the binary cell model allow for asymmetric charge configurations ($E_R \neq 0$), which are excluded from the standard symmetric charge (single-cell) cell mode ($E_R = 0$). In the latter case the binary cell model then reduces to the standard cell model.</figcaption></figure>

of adsorption sites is related to the number $K = 4\pi R_0^2/a^2$ of adsorption sites on a single colloidal macroion.

As in Refs. [25,27,28] we base our macroion CR model on the Frumkin-Fowler-Guggenheim isotherm [29] of the macroion surface defined with the phenomenological free energy of a single adsorption site in the units of thermal energy $\beta = 1/k_B T$ as

$$ \beta\hat{\Omega}^s(\eta) = -\alpha\eta - \frac{\chi}{2}\eta^2 + \eta \ln(\eta) + (1 - \eta)\ln(1 - \eta). \quad (1) $$

The parameters $\alpha$ and $\chi$ are phenomenological and describe the nonelectrostatic part of the proton-macroion and the proton-proton interactions at the macroion surface. In the case of (de)protonation reaction, the dependence of $\alpha$ on the bulk pH is model specific [22], but one can explicitly identify $\alpha = (pK-pH) \ln 10$, where pH = $-\log_{10}[H^+]$, with $[H^+]$ being the proton concentration in the bulk and pK is the dissociation constant of the deprotonation reaction, in the case of the Langmuir adsorption model [30]. Furthermore, $\chi$, as in the related lattice regular solutions theories (e.g., the
Flory-Huggins theory [31]), describes the short-range interactions between nearest-neighbor adsorption sites on the macroion surface [23]. An increase in the parameter value $\alpha$ encodes a favorable adsorption free energy between protons and the macroion surface, while $\chi \ge 0$ represents the tendency of protons on the macroion surface adsorption sites to phase separate into domains.

In what follows, we will use both $\alpha$ as well as $\chi$ as purely phenomenological interaction parameters, quantifying the adsorption energy in the surface (de)protonation reactions and the nearest-neighbor surface energy of filled surface adsorption sites.

The model (1) was applied to lamellar-lamellar phase transition in a charged surfactant system [28] and a good correspondence with experiments was obtained for the didodecyldimethylammonium chloride (DDACI) data with $\alpha = -3.4$, $\chi = 14.75$, and for the didodecyldimethylammonium bromide (DDABr) data with $\alpha = -7.4$ and $\chi = 14.75$ (see Ref. [28] for details). The same model was successfully applied also to other systems (see, e.g., [32,33]).

### III. DENSITY FUNCTIONAL THEORY IN THE LOW-DENSITY LIMIT

#### A. Formalism

As the configuration of a single colloidal macroion is described by the position $\mathbf{r} \in \mathcal{V}$ of the center of mass and the average degree of protonation $\eta \in [0, 1]$ on its surface, the whole suspension can be described by the number density $n(\mathbf{r}, \eta)$. The equilibrium number density minimizes a grand canonical density functional $\Omega[n]$ (see Ref. [34]), which is approximated in the low-density limit by

$$ \beta\Omega[n] = \int_{\mathcal{V}} d^3\mathbf{r} \int_0^1 d\eta n(\mathbf{r}, \eta) \times \left[ \ln\left(\frac{n(\mathbf{r}, \eta)}{\zeta}\right) - 1 + K\beta\hat{\Omega}^s(\eta) \right] + \beta F_{\text{hc}}^{\text{ex}}[n] + \beta F_{\text{el}}^{\text{ex}}[n]. \quad (2) $$

Here, $\zeta$ is the fugacity, $F_{\text{hc}}^{\text{ex}}$ represents the excess free energy due to the hard-core interaction between two colloidal macroions, and $F_{\text{el}}^{\text{ex}}$ describes the excess free-energy contribution of the electrostatic interaction. In the following, the hard-core excess free energy $F_{\text{hc}}^{\text{ex}}$ is based on the Percus-Yevick (PY) closure and the corresponding equation of state via the compressibility route is used [35,36].

The colloidal macroions are assumed to be suspended in an electrolyte solution with relative permittivity $\epsilon_r$ and Debye length $1/\kappa$. For not too highly charged macroions in a sufficiently dilute suspension one can use the Debye-Hückel (DH) approximation [36,37] for the electrostatic two-particle interaction potential

$$ \beta U_{\text{el}}(r, \eta, \eta') = \sigma^*(\eta)\sigma^*(\eta')K^2 \frac{\ell_B \exp[-\kappa(r-R_0)]}{(1+\kappa R_0)r}, \quad (3) $$

where here and below $r = |\mathbf{r} - \mathbf{r}'|$, while the dimensionless surface-charge density

$$ \sigma^*(\eta) = \frac{\sigma(\eta)a^2}{e} = \eta - \frac{1}{2} \quad (4) $$

of a colloidal macroion with average degree of protonation $\eta$ and the Bjerrum length $\ell_B = \beta e^2/(4\pi\epsilon_0\epsilon_r)$ of the solvent with the vacuum permittivity $\epsilon_0$ are introduced. Considering the electrostatic interaction $U_{\text{el}}$ as a perturbation of the hard-core interaction $U_{\text{hc}}$ with

$$ \beta U_{\text{hc}}(r) = \begin{cases} 0 & \text{for } r \ge 2R_0, \\ \infty & \text{for } r < 2R_0, \end{cases} $$

one obtains in the low-density limit [34]

$$ \beta F_{\text{el}}^{\text{ex}}[n] = \frac{1}{2} \int_{\mathcal{V}} d^3\mathbf{r} \int_{\mathcal{V}} d^3\mathbf{r}' \int_0^1 d\eta \int_0^1 d\eta' n(\mathbf{r}, \eta) n(\mathbf{r}', \eta') \\ \times \exp[-\beta U_{\text{hc}}(\mathbf{r} - \mathbf{r}')] \\ \times \{1 - \exp[-\beta U_{\text{el}}(\mathbf{r} - \mathbf{r}', \eta, \eta')]\}. \quad (5) $$

The considered model is then specified by the following five parameters: $\alpha$, $\chi$, $\kappa R_0$, $\kappa \ell_B$, $K$, among which $\alpha$ and $\chi$ describe the charge regulation (according to Sec. II). The values for the parameters $\alpha$ and $\chi$ are chosen keeping in mind that for $\chi = -2\alpha$ the surfaces remain charge neutral for $\chi < \chi_c$ below a certain critical value $\chi_c > 0$, whereas they can be oppositely charged for $\chi > \chi_c$ [25]. Assuming spherical colloidal macroions of radius $R_0 = 10$ nm in an aqueous electrolyte solution with ionic strength 1 mM, i.e., with Bjerrum length $\ell_B \approx 0.7$ nm and Debye length $1/\kappa \approx 10$ nm, leads to the values $\kappa R_0 \approx 1$ and $\kappa \ell_B \approx 0.07$. Finally, within the present DFT approach, $K \in \{20, 40, 45, 50\}$ adsorption sites per colloidal macroion are considered. This corresponds to surface areas per adsorption site $a^2 = 4\pi R_0^2/K \in \{63, 31, 28, 25\}$ nm$^2$, i.e., average distances between neighboring adsorption sites of $a \in \{7.9, 5.6, 5.3, 5\}$ nm, respectively.

#### B. Charge-regulation-induced phase separation

In the bulk of the colloidal suspension, no position dependence occurs for the equilibrium density profile as well as for the total bulk packing fraction, i.e., $n(\mathbf{r}, \eta) = n_b(\eta)$ and $\bar{\Phi}(\mathbf{r}) = \bar{\Phi}_b$. Note that uniformity of the bulk profiles of a fluid is not an assumption or an approximation but the necessary consequence of translational invariance. Upon solving the Euler-Lagrange equation (A4) given in Appendix A one obtains the bulk packing fraction profile

$$ \Phi_b(\eta) = \frac{4\pi}{3} R_0^3 n_b(\eta), \quad (6) $$

which provides the distribution of the average degree of protonation $\eta$ or, equivalently, of the surface-charge densities $\sigma^*(\eta)$ [see Eq. (4)]. Figure 2 displays this distribution for a suspension with bulk packing fraction $\bar{\Phi}_b = 0.1$ of colloidal macroions with $K = 20$ adsorption sites per macroion. For small values of $\chi > 0$, e.g., $\chi = 2$ (see the blue and the green curves in Fig. 2), the surface-charge distribution is unimodal, i.e., the colloidal macroions are essentially equally charged. If the charge-regulation parameters $\alpha$ and $\chi$ [see Eq. (1)] fulfill the relation $\chi = -2\alpha$ the peak is at $\sigma^* = 0$ (see the blue curve in Fig. 2), whereas for $\alpha \ge -\chi/2$ the majority of colloidal macroions carry a surface charge $\sigma^* \ge 0$ (see the green curve in Fig. 2). Upon increasing the value of the charge-regulation parameter $\chi$ the surface-charge distribution becomes bimodal with the peaks being located at increasingly
<figure><img src="image_4.png" /><figcaption>FIG. 2. Distribution of surface charges σ* and degrees of protonation η as well as the associated bulk packing fraction profile Φb(η) [as defined in Eq. (6)] in a suspension with packing fraction Φb = 0.1 of colloidal macroions with K = 20 adsorption sites per macroion. For small values of χ > 0 a unimodal distribution is observed, whereas for sufficiently large values of χ the surface-charge distribution becomes bimodal with increasingly large magnitudes of the peak surface charges.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 3. Binodals of the charge-regulation-induced phase separation of colloidal macroions with K = 40 (blue, light shade), 45 (green, intermediate shade), and 50 (purple, dark shade) adsorption sites per macroion. The charge-regulation parameter χ = 20 is chosen arbitrarily and the pH-sensitive parameter α ∈ [-13, -7] is tuned around the value -χ/2 = -10, where oppositely charged colloidal macroions are expected to occur. The interior of the loops corresponds to the two-phase regions, where phase separation into a low- and a high-density phase occurs at the given value of α. The two-phase region widens upon increasing the number K of adsorption sites per macroion as a result of an increasing magnitude of the electrostatic interaction.</figcaption></figure>

large magnitudes of the surface charges (see the yellow and the purple curves with the lightest and the darkest shades, respectively, in Fig. 2). For α = -χ/2 both peaks represent the same number of macroions, but of opposite charge. The presence of equal amounts of oppositely charged colloidal macroions is expected to lead to compact structures, i.e., to a high-density phase.

Recently, Avni et al. [23] described a two-phase (or even multiple-phase) coexistence region(s), where macroions with low-adsorption site occupation coexist with macroions with high site occupation, akin to the case presented in Fig. 2. However, the model in Ref. [23] differs from the present one in that there are two types of adsorption sites, one charging positively and one charging negatively, on the macroions, whereas here the negative surface charges are fixed and only the positively charging sites are charge regulated.

In order to illustrate the occurrence of a phase separation into a high- and a low-density phase for α ≈ -χ/2, the case χ = 20 for various numbers K of adsorption sites per colloidal macroion is considered. Figure 3 displays the binodals of the charge-regulation-induced phase separation transition for K = 40 (blue curve, light shade), 45 (green curve, intermediate shade), and 50 (purple curve, dark shade). The interior of the loops corresponds to the two-phase regions, where phase separation into a low- and a high-density phase at the given value α occurs. The two-phase region widens upon increasing the number K of adsorption sites per macroion as a result of an increasing magnitude of the electrostatic interaction.

It is well known that there are no fluid phases with packing fractions above Φb ≈ 0.5 as then crystallization sets in. This phenomenon is not covered within the present framework so that values of Φb ≳ 0.5 here are indicative of colloidal aggregation.

The binodals of the charge-regulation-induced phase separation presented in Fig. 3 are quite similar to those calculated by Adame-Arana et al. [21], even if the calculational details

differ, the main difference being that we include the electrostatic interactions explicitly via the two-body DH interaction, Eq. (3), while in Ref. [21] the charge-charge interaction is characterized by a phenomenological Flory-type parameter.

### C. Fluid structure

The bulk structure of the considered suspensions of charge-regulated colloidal macroions described by the density functional in Eq. (2) can be inferred from the partial structure factor S(q, η, η') (see Appendix A for details). It can be conveniently analyzed in terms of the number-number structure factor

$$S_{NN}(q) = \int_{0}^{1} d\eta \int_{0}^{1} d\eta' S(q, \eta, \eta'), \quad (7)$$

which describes the relative distribution of colloidal macroions irrespective of their charge, and the charge-charge structure factor

$$S_{ZZ}(q) = \int_{0}^{1} d\eta \int_{0}^{1} d\eta' \sigma^{*}(\eta) \sigma^{*}(\eta') S(q, \eta, \eta'), \quad (8)$$

which describes the relative distribution of charge within the fluid.

Figure 4 displays the structure factors SNN(q) and SZZ(q) for suspensions with packing fractions Φb ∈ {0.1, 0.2, 0.3} of colloidal macroions with K = 20 adsorption sites per macroion and charge-regulation parameters α = -5, χ = 10. The functional form of SNN(q) in Fig. 4(a) indicates a fluid structure of the colloidal suspension with an increasingly pronounced neighbor shell structure upon increasing the packing fraction Φb. In parallel, the form of SZZ(q) in Fig. 4(b) indicates a spatially alternating arrangement of oppositely charged
<figure><img src="image_6.png" /><figcaption>FIG. 4. Structure factors $S_{NN}(q)$ (a) and $S_{ZZ}(q)$ (b) of suspensions with packing fractions $\bar{\Phi}_b = 0.1$ (blue, light shade), 0.2 (green, intermediate shade), and 0.3 (purple, dark shade) of charge-regulated colloidal macroions of radius $R_0$ with $K = 20$ adsorption sites per macroion and charge-regulation parameters $\alpha = -5$, $\chi = 10$. As defined in Eqs. (7) and (8), the structure factor $S_{NN}(q)$ in (a) describes the relative distribution of colloidal macroions irrespective of their charge, whereas $S_{ZZ}(q)$ in (b) describes the relative distribution of charge within the fluid (see Table I for details). The location of the main peak of $S_{ZZ}(q)$ relative to that of $S_{NN}(q)$ indicates an alternating charge structure.</figcaption></figure>

colloidal macroions with a periodicity of approximately twice the nearest-neighbor distance.

A hypothetical macrophase separation of charges would lead to a peak of $S_{ZZ}(q)$ at $q = 0$, which does obviously

<table><caption>TABLE I. Characteristics (see main text) of suspensions of packing fractions $\bar{\Phi}_b = \{0.1, 0.2, 0.3\}$ of charge-regulated spherical colloidal macroions of radius $R_0$ with $K = 20$ adsorption sites per macroion and charge-regulation parameters $\alpha = -5$, $\chi = 10$ inferred from the structure factors $S_{NN}(q)$ [see Fig. 4(a)] and $S_{ZZ}(q)$ [see Fig. 4(b)]. Upon increasing the packing fraction $\bar{\Phi}_b$ the mean nearest-neighbor distance $\lambda_{NN}$ decreases to almost the close-contact distance $2R_0$ of two hard spheres of radius $R_0$, whereas the periodicity $\lambda_{ZZ}$ of the charge distribution is slightly larger than twice that distance, $\lambda_{ZZ} \gtrsim 2\lambda_{NN}$. In parallel, the coordination number $N_1$ of colloidal macroions in the nearest-neighbor shell increases.</caption><thead><tr><th>Φ<sub>b</sub></th><th>λ<sub>NN</sub>/R<sub>0</sub></th><th>ξ<sub>NN</sub>/R<sub>0</sub></th><th>λ<sub>ZZ</sub>/R<sub>0</sub></th><th>ξ<sub>ZZ</sub>/R<sub>0</sub></th><th>N<sub>1</sub></th></tr></thead><tbody><tr><td>0.1</td><td>2.40</td><td>0.66</td><td>5.80</td><td>0.87</td><td>3.4</td></tr><tr><td>0.2</td><td>2.21</td><td>0.98</td><td>4.82</td><td>1.11</td><td>6.1</td></tr><tr><td>0.3</td><td>2.06</td><td>1.44</td><td>4.39</td><td>1.34</td><td>8.3</td></tr></tbody></table>

not occur here. From the position of the major peaks in Figs. 4(a) and 4(b) one obtains the mean nearest-neighbor distance $\lambda_{NN}$ between two colloidal macroions and the periodicity $\lambda_{ZZ}$ of the charge distribution, respectively; the values are displayed in Table I. As is expected from the nature of the hard-core repulsion between two colloidal macroions, the nearest-neighbor distance $\lambda_{NN}$ is never smaller than the close-contact distance $2R_0$. Upon increasing the packing fraction $\bar{\Phi}_b$ both $\lambda_{NN}$ and $\lambda_{ZZ}$ decrease, but the relation $\lambda_{ZZ} \gtrsim 2\lambda_{NN}$ holds for all cases. Hence, the charge distribution oscillates on a length scale $\lambda_{ZZ}$ which is approximately twice the nearest-neighbor distance $\lambda_{NN}$. This is again showing the alternating charge structure within the suspension of charge-regulated spherical macroions. The widths of the major peaks in Figs. 4(a) and 4(b) yield, respectively, the decay lengths (correlation lengths) $\xi_{NN}$ and $\xi_{ZZ}$ of the structural features, which are of the order of the radius $R_0$ of the colloidal macroions, as can be observed in Table I. Finally, the coordination number $N_1$, i.e., the number of colloidal macroions in the nearest-neighbor shell, can be calculated from $S_{NN}(q)$; the corresponding values are displayed in Table I.

## IV. BINARY CELL MODEL

The DFT calculations presented so far work the best for dilute suspension of charged colloidal macroions. While it can also be expected to work well for an aggregating system of weakly charged macroions, a dense suspension of strongly charged macroions surely needs to be treated differently. Accordingly, in this section we present a variation on the standard cell model that we refer to as the binary cell model (BCM) that can be invoked to describe a dense suspension of identical charge-regulated macroions irrespective of their surface-charge densities and is consequently valid for symmetric as well as asymmetric charge configurations (see Fig. 1). Contrary to the standard cell model as well as the cell model in binary colloidal mixtures [38], the building blocks of our variety of the cell model are two adjacent cells of radius $R$ each of which encloses a charged particle of radius $R_0$ [see Fig. 1(b)]. The macroions as well as the wrapping cells are considered to be fixed in space. Such a model inherently assumes a certain infinite-range crystalline order and the validity of such assumption clearly improves the closer the system is to a high-density crystalline state.

The model medium built from this elementary cell construct thus has a particle-cell volume ratio $\Phi_v = (R_0/R)^3$ which is related to the bulk packing fraction $\bar{\Phi}_b$ defined earlier via the relation $\Phi_v = \bar{\Phi}_b/\Phi_{cp}$ where $\Phi_{cp} \approx 0.74$ is the packing fraction corresponding to the face-centered-cubic or hexagonal-close-packed arrangement of the cells. Clearly, $\Phi_v$ is inversely proportional to the cube of the intermacroion separation. The charges on the macroion surfaces are again regulated according to the mechanism introduced in Sec. II [see Eq. (1)].

The following considerations are based on a Poisson-Boltzmann (PB) theory of the binary cell system. The grand potential corresponding to a single cell, in units of the thermal energy $\beta = 1/k_B T$, can be written as

$$ \beta\hat{\Omega}(\eta, E_R) = \beta\hat{\Omega}^{\text{el}}(\eta, E_R) + K\beta\hat{\Omega}^{\text{s}}(\eta), \quad (9) $$
where $\beta\hat{\Omega}^s(\eta)$ is the energy contribution stemming from the chemical processes driving the (de)protonation reaction [as defined in Eq. (1)] at the macroion surface and

$$ \frac{\beta\hat{\Omega}^{\text{el}}(\eta, E_R)}{4\pi R_0^2} = \frac{-\epsilon}{\beta e^2} \int_{R_0}^{R} dr \frac{r^2}{R_0^2} \left[ \kappa^2 \cosh[\phi(r)] + \frac{1}{2}[\phi'(r)]^2 \right] + \frac{\phi(R_0)}{a^2} \left( \eta - \frac{1}{2} \right) - \frac{\epsilon}{e} \frac{R^2}{R_0^2} E_R \phi(R) \quad (10) $$

with the permittivity $\epsilon = \epsilon_r\epsilon_0$ of the embedding medium is the electrostatic part of the grand potential expressed per unit surface area of the macroion [39]. Herein, $\phi(r)$ is the dimensionless electrostatic potential expressed in the units of $\beta e$ which fulfills the PB equation in spherical symmetry,

$$ \frac{1}{r^2}(r^2\phi'(r))' = \kappa^2 \sinh(\phi(r)), $$

subject to the boundary condition of a charge density $\sigma$ at the macroion surface, i.e., at $r = R_0$,

$$ \phi'(R_0) = -\frac{\beta e^2}{\epsilon a^2} \left( \eta - \frac{1}{2} \right) = -\frac{\beta e \sigma}{\epsilon}, \quad (11) $$

and of a given radial component $E_R$ of the electric field at the cell boundary, i.e., at $r = R$,

$$ \phi'(R) = -\beta e E_R. \quad (12) $$

For two coupled cells the grand potential of the binary cell system can then be written as

$$ \beta\Omega(\eta_1, \eta_2, E_R) = \beta\hat{\Omega}(\eta_1, E_R) + \beta\hat{\Omega}(\eta_2, -E_R), \quad (13) $$

where the subscripts “1” and “2” are used to indicate the two adjacent cells in the binary cell model. The energy contribution $\beta\hat{\Omega}$ for the second cell is identical to the first one and is given by Eq. (9), albeit with the electric field $E_R$ replaced by $-E_R$ as the centers of the two adjacent cells imply a different unit normal at the boundary (see Fig. 1).

### A. Debye-Hückel case

First, we consider the Debye-Hückel case, i.e., the linearized PB equation, which renders the problem analytically tractable in part and also straightforwardly allows to scan the phase space spanning over the whole ranges of $(\eta_1, \eta_2)$ or equivalently $(\sigma_1^*, \sigma_2^*)$.

The exact solution of the linear electrostatic problem implies a nonlinear function

$$ \beta\Omega(\eta_1, \eta_2) = \min_{E_R} \beta\Omega(\eta_1, \eta_2, E_R) $$

of the degrees of protonation $\eta_1, \eta_2$, of the binary cell model, which subsequently needs to be minimized numerically. In order to achieve this, one can proceed as described in the Supplemental Material of [25], where the linearized PB equation is discussed in a planar geometry. In the present case we can make use of the known solution of the linearized PB equation in spherical geometry [40]; details are described in Appendix B. After inserting the equilibrium $E_R$ that minimizes the $\beta\Omega(\eta_1, \eta_2, E_R)$ obtained by expanding Eq. (10) up to second order in $\phi$ and replacing $\phi(r)$ by the known solutions, one finally arrives at the following expression for

the grand potential defined in Eq. (13) as a function of $\eta_1$ and $\eta_2$ only:

$$ \beta\Omega(\eta_1, \eta_2) = \frac{K}{2} \frac{\ell_B}{a^2 |\det \mathbf{M}|} \left[ \gamma \left( \eta_1 - \frac{1}{2} \right)^2 + \gamma \left( \eta_2 - \frac{1}{2} \right)^2 - \frac{\tau^2}{2|\nu|} (\eta_1 - \eta_2)^2 \right] + K \beta [\hat{\Omega}^s(\eta_1) + \hat{\Omega}^s(\eta_2)]. \quad (14) $$

As before, $K = 4\pi R_0^2/a^2$, and $\ell_B$ is the Bjerrum length. The factors $|\det \mathbf{M}|$, $\gamma$, $\tau$, and $\nu$ involve the three length scales of the problem: the radius of the macroion, the Debye length, and the cell size; the analytic expressions are given in Appendix B.

In line with the DFT calculation and the inherent approximations of DH theory, we consider the limit of small $K$ only. In this regime the binary cell model allows for asymmetric charge configurations already at smaller values of $(\alpha, \chi)$ parameters than within the DFT approach in Sec. III. The reason for this is a weaker electrostatic coupling between two colloidal macroions in Sec. III, which is based on a superposition approximation of the electrostatic interaction, as compared to the stronger coupling via the electric field $E_R$ at the common boundary between two adjacent cells within the binary cell model [41,42].

The phase coexistence around the symmetry axis $-2\alpha = \chi$ is brought about by an exchange of stability of the minima of the grand potential, as illustrated in Fig. 5. For a given particle-cell volume ratio $\Phi_v$, the charge state with minimum energy shifts from a symmetric to an asymmetric one as one moves away from the symmetry axis $-2\alpha = \chi$ [compare Figs. 5(a) and 5(b)]. Moreover, a comparison Figs. 5(a) and 5(c) suggests that for given $(\alpha, \chi)$, as the particle-cell volume ratio $\Phi_v$ diminishes, the difference in the grand potential between the symmetric and asymmetric configurations gradually diminishes too, leading ultimately to a symmetric equilibrium state.

It is worth noting that the linear theory ceases to be valid not only for high-$K$ values, but also for higher values of the packing fraction $\bar{\Phi}_b$ (or, equivalently, the particle-cell volume ratio $\Phi_v$) where the following nonlinear PB theory needs to be applied.

### B. Full Poisson-Boltzmann case

Within the nonlinear PB theory, the equilibrium $\eta$ values [or equivalently $\sigma^*$ values; see Eq. (4)] at the two surfaces along with the electric field at the cell boundary $E_R$ are obtained via a numerical minimization of the grand potential following the scheme described in the beginning of Sec. IV. We consider a system of macroions dispersed in an aqueous electrolyte solution with ionic concentration $I = 1$ mM, relative permittivity $\epsilon_r \approx 80$ at temperature $T = 300$ K.

The resulting variations of the degrees of protonations $\eta_1$ and $\eta_2$ leading to charge densities $\sigma_1^*$ and $\sigma_2^*$ at the two macroion surfaces are shown in Fig. 6 as functions of the interaction parameter $\alpha$ and the particle-cell volume ratio $\Phi_v$ for two different values of the parameter $K \in \{50, 1257\}$ corresponding to $a \in \{5, 1\}$ nm and $\chi = 30$. Note that unlike in the case of the linear DH theory, we are not constrained by any upper limit for $K$ here. Nevertheless, as one can infer from
<figure><img src="image_7.png" /><figcaption>FIG. 5. Top panel shows the variation of the grand potential $\beta\Omega(\sigma_1^*, \sigma_2^*)/K$ given by Eq. (14) as functions of the surface-charge density variables $\sigma_1^*$ and $\sigma_2^*$ for different values of the parameters $\alpha$ and $\Phi_v$. Bottom panel shows cuts through the potential surface at $\sigma_1^* = \sigma_2^*$ (in yellow, light shade) and $\sigma_1^* = -\sigma_2^*$ (in blue, dark shade). (a) Shows the results for $\alpha = -3$, $\chi = 6$, and $\Phi_v = 0.6$ which fall inside the phase separation region. Consequently, the asymmetric configuration is the minimum-energy configuration. (b) At the same value of $\Phi_v$, but with $\alpha$ changing from $-3$ to $-2.8$. As it is evident from the curves, the stability has changed toward a configuration with uncoupled cells, thus leaving the phase coexistence region. (c) Again, at $\alpha = -3$ and $\chi = 6$ but going down in the particle-cell volume ratio to $\Phi_v = 0.05$. As one can see, the grand potential of the symmetric configuration approaches that of the asymmetric configuration. For all the plots, $R_0 = 10$ nm, $\kappa = 0.1$ nm$^{-1}$, $\ell_B = 0.7$ nm, and $K = 20$ are used.</figcaption></figure>

Fig. 6, the results are qualitatively similar to those obtained within the linearized PB theory in the previous section and are

<figure><img src="image_8.png" /><figcaption>FIG. 6. Map showing the charge states of the two macroions as functions of the interaction parameter $\alpha$ and the particle-cell volume fraction $\Phi_v$ for $\chi = 30$ and two different $K$ values. The solid line corresponds to $K = 1257$ whereas the dashed line refers to $K = 50$. In each case, one observes a conical-shaped region centered around $\alpha = -\chi/2 = -15$ with stretched opening at the top. Inside this region and on the boundary, the two macroion surfaces are oppositely charged ($|\eta_1 - \eta_2| = 1$ or, equivalently, $\sigma_1^* = -\sigma_2^*$) whereas outside this region they are identically charged ($\eta_1 = \eta_2$ or, equivalently, $\sigma_1 = \sigma_2$). With increasing $K$ value, this region featuring charge asymmetry widens.</figcaption></figure>

consistent with the findings reported in Fig. 5. For any given $K$ value, one obtains a triangularly shaped region centered around $\alpha = -\chi/2$ with stretched opening at the top. Inside this region and on the boundary, the two macroion surfaces are oppositely charged ($\sigma_1^* \approx -\sigma_2^* \neq 0$) whereas outside this region they are identically charged ($\sigma_1^* = \sigma_2^*$). In accordance with the outcomes of the linear theory in Sec. III, within these identically charged regions, both the surfaces are negatively charged for $\alpha < \chi/2 = -15$, whereas for $\alpha > \chi/2 = -15$ they are positively charged.

With increasing $K$ value, the region with asymmetric charge configurations broadens as an increase in $K$ implies a higher surface-charge density, which in turn enhances the electrostatic attraction between the surfaces at the origin of the observed symmetry breaking. Although, in general, for a given $\alpha$ value the asymmetric configuration changes to the symmetric one with decreasing volume fraction (or, equivalently, increasing separation between the macroions), the asymmetric configurations observed on the line $\alpha = -\chi/2 = -15$ are very stable and persist down to $\Phi_v \approx 10^{-3}$ or lower. The specific choice of the parameter $\chi = 30$ is motivated by the observation that complete symmetry breaking, i.e., $\sigma_1^* \approx -\sigma_2^* \neq 0$, occurs above a critical value $\chi = \chi_c \approx 25$ for $K = 1257$. It is important to note that this $\chi_c$ value depends not only on $K$, but can be different within the linear and nonlinear PB theories. However, for any $\chi > \chi_c$, one observes the same qualitative features, i.e., transitions from symmetric to symmetry broken states as functions of $\alpha$ and $\Phi_v$ within both the theories.
V. CONCLUSIONS

In summary, we have described a charge-regulation-based mechanism that allows for pH-dependent phase separation in macroion solutions. A complex interplay of the different chemical interactions driving the charge regulation of individual macroion, as well as electrostatic interaction between them, leads to symmetry broken charge states of the macroions, thereby leading to aggregation. A density functional theory based approach, applicable for dilute suspension of macroions, fully accounting for the translational entropy of the macroions, indeed provides evidence of electrostatically driven macroion phase separation. A binary cell model Poisson-Boltzmann approach, applicable in the high-density limit, confirms the presence of electrostatic attraction, essential for the observed phase separation, stemming from a transition to asymmetrically charged states of nearest-neighbor macroion pairs.

An interesting aside to our calculation is the way it naturally ties together the well-studied complex coacervation between chemically oppositely charged macroions and the much less studied simple coacervation of chemically identical macroions in bathing electrolyte solutions, by introducing the actual chemical model of charging, as opposed to a priori chosen values of the surface charge (potential). We propose this charge symmetry broken complex coacervation between chemically identical macroions as a bridge between the two different types of coacervations or indeed liquid-liquid phase separations.

A systematic study based on the two contrasting approaches, the DFT and the BCM within the same charge-regulation model, therefore ensures the robustness of our results and allows us to conclude that the pH-dependent liquid-liquid phase separation in macroion solutions is a rule rather than exception, even in the case of chemically identical macroions. Further indications for the robustness

and generality of the described results are qualitative similarities with approaches based on alternative computational methodologies, such as the phenomenological Flory-type electrostatics [21] and the collective mean-field description [23]. It still remains to be seen which are the absolutely essential ingredients of the macroion surface-charge regulation promoting this liquid-liquid phase separation.

Finally, recent advances in the simulations of acid-base equilibria in systems coupled to a reservoir with a fixed pH, based either on a hybrid Monte Carlo method to resolve the charges of individual surface groups [43], on the grand-reaction method for coarse-grained simulations of acid-base equilibria with a fixed pH reservoir and salt concentration [44], or simulating the pH effects with classical coarse-grained molecular dynamics simulations [45], could in principle provide a proper background to different analytical approaches and hopefully elucidate the reality of the predicted phenomena. The comparison with experiments, in which other types of interactions and the notoriously difficult solvent effects come into play, poses another challenge that will have to be faced in the future.

ACKNOWLEDGMENTS

R.P. would like to acknowledge the support of the 1000-Talents Program of the Chinese Foreign Experts Bureau. He would also like to thank Y. Avni for her valuable comments.

APPENDIX A: DENSITY FUNCTIONAL THEORY

1. Bulk packing fraction

Upon expanding the term $1 - \exp(-\beta U_{el})$ in Eq. (5) of the main text in powers of $\beta U_{el}$ and using Eq. (3) of the main text one can perform both integrations over $V$ in Eq. (5) to obtain

$$ \beta F_{el}^{ex}[n_b] = \frac{3V}{4\pi R_0^3} \sum_{k=1}^{\infty} B_k \left( \int_0^1 d\eta \sigma^*(\eta)^k \Phi_b(\eta) \right)^2 \quad (A1) $$

with the system volume $V = |V|$, the constants

$$ B_k := \frac{3}{2} \frac{(-1)^{k-1} k^{k-3}}{k!(\kappa R_0)^3} \Gamma(3-k, 2k\kappa R_0) \left( K^2 \frac{\kappa \ell_B}{1 + \kappa R_0} \exp(\kappa R_0) \right)^k, $$

where $\Gamma(\nu, z)$ denotes the incomplete $\Gamma$ function [46], and the bulk packing fraction profile

$$ \Phi_b(\eta) = \frac{4\pi}{3} R_0^3 n_b(\eta) \quad (A2) $$

with $\bar{\Phi}_b = \int_0^1 d\eta \Phi_b(\eta)$. Using Eqs. (A1) and (A2) in Eq. (2) of the main text one obtains a scaled density functional in terms of the bulk packing fraction profiles $\Phi_b$:

$$ \Omega_b^*[\Phi_b] := \beta \Omega[n_b] \frac{4\pi R_0^3}{3V} = \int_0^1 d\eta \Phi_b(\eta) (\ln[\Phi_b(\eta)] + K\beta\hat{\Omega}^s(\eta)) + \bar{\Phi}_b[-1 - \mu^* + h_{PY}(\bar{\Phi}_b)] \\ + \sum_{k=1}^{\infty} B_k \left( \int_0^1 d\eta \sigma^*(\eta)^k \Phi_b(\eta) \right)^2 \quad (A3) $$
with the scaled chemical potential

$$ \mu^* := \ln \left( \frac{4\pi}{3} R_0^3 \zeta \right) $$

and

$$ h_{\text{PY}}(\bar{\Phi}_b) = -\ln(1 - \bar{\Phi}_b) + \frac{6\bar{\Phi}_b - 2\bar{\Phi}_b^2}{2(1 - \bar{\Phi}_b)^2}. $$

The equilibrium bulk packing fraction profile $\Phi_b^{eq}$ minimizes the scaled density functional $\Omega_b^*$ in Eq. (A3). Hence, it solves the Euler-Lagrange equation [34]

$$ \begin{aligned} 0 &= \frac{\delta \Omega_b^*}{\delta \Phi_b(\eta)} \Big|_{\Phi_b^{eq}} \\ &= \ln[\Phi_b^{eq}(\eta)] + K \beta \hat{\Omega}^s(\eta) - \mu^* + h_{\text{PY}}(\bar{\Phi}_b^{eq}) + \bar{\Phi}_b^{eq} h'_{\text{PY}}(\bar{\Phi}_b^{eq}) \\ &\quad + 2 \sum_{k=1}^{\infty} B_k \int_0^1 d\eta' \sigma^*(\eta')^k \Phi_b^{eq}(\eta') \sigma^*(\eta)^k. \end{aligned} \quad (A4) $$

This expression can be rewritten in the form

$$ \Phi_b^{eq}(\eta) = \exp \left( \mu^* - K \beta \hat{\Omega}^s(\eta) - h_{\text{PY}}(\bar{\Phi}_b^{eq}) - \bar{\Phi}_b^{eq} h'_{\text{PY}}(\bar{\Phi}_b^{eq}) \right. \\ \left. - 2 \sum_{k=1}^{\infty} B_k \bar{\Phi}_b^{eq} \langle \sigma^{*k} \rangle \sigma^*(\eta)^k \right) $$

with the kth moment of the surface-charge distribution

$$ \langle \sigma^{*k} \rangle := \frac{1}{\bar{\Phi}_b^{eq}} \int_0^1 d\eta \sigma^*(\eta)^k \Phi_b^{eq}(\eta). $$

## 2. Partial structure factor

The partial structure factor of the system under consideration can be written in the form

$$ S(q, \eta, \eta') = \frac{1}{n_b} \tilde{G}(q, \eta, \eta') := \frac{\sqrt{\bar{\Phi}_b(\eta)\bar{\Phi}_b(\eta')}}{\bar{\Phi}_b} \mathcal{G}(q, \eta, \eta'), $$

where $\tilde{G}(q, \eta, \eta')$ is the three-dimensional Fourier transform of the bulk correlation function $\mathcal{G}(r, \eta, \eta')$ of the number densities of colloidal spheres with degrees of protonation $\eta$ and $\eta'$ at a distance $r$. The auxiliary function $\mathcal{G}(q, \eta, \eta')$ fulfills the Ornstein-Zernike equation

$$ \mathcal{G}(q, \eta, \eta') = \int_0^1 d\eta'' \mathcal{C}(q, \eta, \eta'') \mathcal{G}(q, \eta'', \eta') + \delta(\eta - \eta') $$

with $\mathcal{C}(q, \eta, \eta') = \sqrt{n_b(\eta)n_b(\eta')}\tilde{c}(q, \eta, \eta')$, where $\tilde{c}(q, \eta, \eta')$, $q = |\mathbf{q}|$, is the Fourier transform of the bulk direct correlation function

$$ c(\mathbf{r}, \eta, \eta') = c_{\text{PY}}(\mathbf{r}, \bar{\Phi}_b) - \frac{\delta^2 \beta F_{\text{el}}^{\text{ex}}}{\delta n(\mathbf{r}, \eta) \delta n(\mathbf{0}, \eta')} \Big|_{n_b} $$

composed of the Percus-Yevick hard-core contribution $c_{\text{PY}}$ (see Ref. [35]) and the electrostatic contribution obtained from Eq. (5) of the main text (see Ref. [34]).

## APPENDIX B: DEBYE-HÜCKEL THEORY FOR BINARY CELL MODEL

In this Appendix we show how Eq. (14) of the main text can be derived. First, the linearized PB-term integral

$$ \int_{R_0}^{R} dr \, r^2 [(\phi'(r))^2 + \kappa^2 (\phi(r))^2], $$

can be simplified by partial integration since one has

$$ \begin{aligned} & \int_{R_0}^{R} dr (r^2 \phi'(r)) \phi'(r) \\ &= r^2 \phi'(r) \phi(r) \Big|_{R_0}^{R} - \int_{R_0}^{R} dr \, \partial_r (r^2 \partial_r \phi(r)) \phi(r), \end{aligned} $$

whereby the latter term can be transformed via the linearized PB equation and hence cancels out against the term in the original integral. Putting all constants and the boundary conditions in, one is left with the expression

$$ \frac{\beta \hat{\Omega}^{\text{el}}}{4\pi R_0^2} = -\frac{\varepsilon}{2\beta e^2} \left(\frac{R}{R_0}\right)^2 E_R \phi(R) + \frac{\sigma^*}{2a^2} \phi(R_0). $$

Coupling two cells, respecting electroneutrality,

$$ \Omega^{\text{el}} = \hat{\Omega}_{1}^{\text{el}}([\phi_1], E_R, \sigma_1^*) + \hat{\Omega}_{2}^{\text{el}}([\phi_2], -E_R, \sigma_2^*) $$

one has the expression

$$ \frac{\beta \Omega^{\text{el}}}{4\pi R_0^2} = -\frac{\varepsilon}{2\beta e^2} \left(\frac{R}{R_0}\right)^2 E_R[\phi_1(R) - \phi_2(R)] \\ + \frac{1}{2a^2}[\sigma_1^* \phi_1(R_0) + \sigma_2^* \phi_2(R_0)]. \quad (B1) $$

The boundary conditions at the cell particles and surfaces follow from the known exact solution of the electrostatic potential of a single cell with a solute of radius $R_0$, embedded in a spherical cell of radius $R$ which contains a salt solution, as given in [40]:

$$ \phi(r) = \frac{R_0 \phi(R_0)}{r} \cosh[\kappa(r - R_0)] \\ + \frac{R \phi(R) - R_0 \phi(R_0) \cosh[\kappa(R - R_0)]}{\sinh[\kappa(R - r_0)]} \\ \times \sinh[\kappa(r - R_0)]. $$

The boundary values $\phi(R_0)$ and $\phi(R)$ thus follow from the derivatives of $\phi(r)$ at $r = R_0, R$.

In the linear theory we discuss, denoting the vector $\phi' = (\phi'(R_0), \phi'(R))$ of derivatives, one needs to invert the matrix equation $\phi' = \mathbf{M} \cdot \phi$ which can be computed from the derivative of the solution. One finds

$$ \phi_i(R_0) = \frac{1}{\det \mathbf{M}} [\gamma \phi_i'(R_0) - \tau \phi_i'(R)] $$

and

$$ \phi_i(R) = \frac{1}{\det \mathbf{M}} [-\xi \phi_i'(R_0) + \nu \phi_i'(R)], $$
<figure><img src="image_9.png" /><figcaption>FIG. 7. Two-phase surface-charge density coexistence region around the symmetry axis with $\alpha = -3$, at $\chi = 6$ for $K = 20$.</figcaption></figure>

where the boundary values for the two cells are given by Eqs. (11) and (12) of the main text. The parameters in the equations are functions of the three characteristic lengths in the system (solute size, screening length,

cell radius):

$$ \nu = -\left[\frac{1}{R_0} + \frac{\kappa C(R)}{S(R)}\right], \quad \tau = \frac{R}{R_0} \frac{\kappa}{S(R)}, \\ \xi = -\frac{R_0}{R} \frac{\kappa}{S(R)}, \quad \gamma = -\frac{1}{R} + \frac{\kappa C(R)}{S(R)}, $$

where $C(R) = \cosh[\kappa(R - R_0)]$ and $S(R) = \sinh[\kappa(R - R_0)]$. Further,

$$ \det \mathbf{M} = \nu\gamma - \tau\xi. $$

One has $\nu < 0$, $\tau > 0$, $\xi < 0$, $\gamma > 0$, and $\det \mathbf{M} < 0$. After computation of the boundary conditions at the macroion radius $R_0$ and the binary cell radius $R$, the resulting expression (B1) needs to be minimized with respect to $E_R$, which is found to behave as $E_R \sim \sigma_1^* - \sigma_2^*$. Finally, collection of terms leads to Eq. (14) of the main text.

Within DH theory, the surface-charge density coexistence curve has a balloonlike shape (see Fig. 7). The upper part of the surface-charge density coexistence curve widens for increasing $\chi$. The location of the upper critical point of the phase coexistence curve moves from $(\chi, \Phi_v) = (3.6, 0.79)$ to $(\chi, \Phi_v) = (16, 0.97)$, which for $K = 20$ covers the interval in which surface-charge density coexistence exists in the DH limit.

[1] F. Tiebackx, Gleichzeitige Ausflockung zweier Kolloide, Z. Chem. Ind. Kolloide 8, 198 (1911).

[2] H. G. Bungenberg de Jong and H. R. Kruyt, Coacervation (Partial miscibility in colloid systems), Proc. Kon. Ned. Akad. Wet. 32, 849 (1929).

[3] C. Brangwynne, P. Tompa, and R. Pappu, Polymer physics of intracellular phase transitions, Nat. Phys. 11, 899 (2015).

[4] C. E. Sing and S. L. Perry, Recent progress in the science of complex coacervation, Soft Matter 16, 2885 (2020).

[5] P. Zhang, N. M. Alsaifi, J. Wu, and Z.-G. Wang, Polyelectrolyte complex coacervation: Effects of concentration asymmetry, J. Chem. Phys. 149, 163303 (2018).

[6] S. L. Perry, Phase separation: Bridging polymer physics and biology, Curr. Opin. Colloid Interface Sci. 39, 86 (2019).

[7] N. Chen, T. Nicolai, C. Chassenieux, and Y. Wang, pH and ionic strength responsive core-shell protein microgels fabricated via simple coacervation of soy globulins, Food Hydrocolloids 105, 105853 (2020).

[8] A. M. Smith, M. Borkovec, and G. Trefalt, Forces between solid surfaces in aqueous electrolyte solutions, Adv. Colloid Interface Sci. 275, 102078 (2020).

[9] B. P. Lee, P. B. Messersmith, J. N. Israelachvili, and J. H. Waite, Mussel-inspired adhesives and coatings, Ann. Rev. Mater. Res. 41, 99 (2011).

[10] S. Kim, J. Huang, Y. Lee, S. Dutta, H. Y. Yoo, Y. M. Jung, Y. Jho, H. Zeng, and D. S. Hwang, Complexation and coacervation of like-charged polyelectrolytes inspired by mussels, Proc. Natl. Acad. Sci. USA 113, E847 (2016).

[11] S. Srivastava and M. V. Tirrell, Polyelectrolyte complexation Advances in Chemical Physics, edited by S. A. Rice and A. R. Dinner, Advances in Chemical Physics, Vol. 161 (John Wiley & Sons, 2016), pp. 499–544.

[12] M. Lund and B. Jönsson, Charge regulation in biomolecular solution, Q. Rev. Biophys. 46, 265 (2013).

[13] A. A. Hyman, C. A. Weber, and F. Jülicher, Liquid-liquid phase separation in biology, Annu. Rev. Cell Dev. Biol. 30, 39 (2014).

[14] T. M. Franzmann, M. Jahnel, A. Pozniakovsky, J. Mahamid, A. S. Holehouse, E. Nüske, D. Richter, W. Baumeister, S. W. Grill, R. V. Pappu, A. A. Hyman, and S. Alberti, Phase separation of a yeast prion protein promotes cellular fitness, Science 359, eaao5654 (2018).

[15] M. Borkovec, B. Jönsson, and G. J. Koper, Ionization Processes and Proton Binding in Polyprotic Systems: Small Molecules, Proteins, Interfaces, and Polyelectrolytes, in Surface and Colloid Science, Vol. 16, edited by E. Matijevic (Springer, Boston, 2001).

[16] V. G. Taratuta, A. Holschbach, G. M. Thurston, D. Blankschtein, and G. B. Benedek, Liquid-liquid phase separation of aqueous lysozyme solutions: Effects of pH and salt identity, J. Phys. Chem. 94, 2140 (1990).

[17] Y. Avni, D. Andelman, and R. Podgornik, Charge regulation with fixed and mobile charged macromolecules, Curr. Opin. Electrochem. 13, 70 (2019).

[18] T. Gisler, S. F. Schulz, M. Borkovec, H. Sticher, P. Schurtenberger, B. D’Aguanno, and R. Klein, Understanding colloidal charge renormalization from surface chemistry: Experiment and theory, J. Chem. Phys. 101, 9924 (1994).

[19] N. Boon and R. van Roij, Charge reversal of moisturous porous silica colloids by take-up of protons, J. Colloid Interface Sci. 385, 66 (2012).

[20] R. J. Allen and B. W. Patrick, Complexation and phase behavior of oppositely charged polyelectrolyte/macroion systems, Langmuir 20, 1997 (2004).
[21] O. Adame-Arana, C. A. Weber, V. Zaburdaev, J. Prost, and F. Jülicher, Liquid phase separation controlled by pH, *Biophys. J.* **119**, 1590 (2020).

[22] Y. Avni, T. Markovich, R. Podgornik, and D. Andelman, Charge regulating macro-ions in salt solutions: Screening properties and electrostatic interactions, *Soft Matter* **14**, 6058 (2018).

[23] Y. Avni, R. Podgornik, and D. Andelman, Critical behavior of charge-regulated macro-ions, *J. Chem. Phys.* **153**, 024901 (2020).

[24] J. C. Everts, M. N. van der Linden, A. van Blaaderen, and R. van Roij, Alternating strings and clusters in suspensions of charged colloids, *Soft Matter* **12**, 6610 (2016).

[25] A. Majee, M. Bier, and R. Podgornik, Spontaneous symmetry breaking of charge-regulated surfaces, *Soft Matter* **14**, 985 (2018).

[26] M. V. Fedorov and A. A. Kornyshev, Ionic liquids at electrified interfaces, *Chem. Rev.* **114**, 2978 (2014).

[27] A. Majee, M. Bier, R. Blossey, and R. Podgornik, Charge regulation radically modifies electrostatics in membrane stacks, *Phys. Rev. E* **100**, 050601(R) (2019).

[28] D. Harries, R. Podgornik, V. A. Parsegian, E. May-Or, and D. Andelman, Ion induced lamellar-lamellar phase transition in charged surfactant systems, *J. Chem. Phys.* **124**, 224702 (2006).

[29] L. Koopal, W. Tan, and M. Avena, Equilibrium mono- and multicomponent adsorption models: From homogeneous ideal to heterogeneous non-ideal binding, *Adv. Colloid Interface Sci.* **280**, 102138 (2020).

[30] B. W. Ninham and V. A. Parsegian, Electrostatic potential between surfaces bearing ionizable groups in ionic equilibrium with physiologic saline solution, *J. Theor. Biol.* **31**, 405 (1971).

[31] I. Teraoka, *Polymer Solutions: An Introduction to Physical Properties* (Wiley, New York, 2002).

[32] L. Fink, J. Feitelson, R. Noff, T. Dvir, C. Tamburu, and U. Raviv, Osmotic stress induced desorption of calcium ions from dipolar lipid membranes, *Langmuir* **33**, 5636 (2017).

[33] L. Fink, A. Steiner, O. Szekely, P. Szekely, and U. Raviv, Structure and interactions between charged lipid membranes in the presence of multivalent ions, *Langmuir* **35**, 9694 (2019).

[34] R. Evans, The nature of the liquid-vapour interface and other topics in the statistical mechanics of non-uniform, classical fluids, *Adv. Phys.* **28**, 143 (1979).

[35] J.-P. Hansen and I. R. McDonald, *Theory of Simple Liquids* (Academic, San Diego, 1986).

[36] D. A. McQuarrie, *Statistical Mechanics* (University Science Books, Sausalito, 2000).

[37] P. Debye and E. Hückel, Zur Theorie der Elektrolyte, *Phys. Z.* **24**, 185 (1923).

[38] A. Torres, G. Téllez, and R. van Roij, The polydisperse cell model: Nonlinear screening and charge renormalization in colloidal mixtures, *J. Chem. Phys.* **128**, 154906 (2008).

[39] T. Markovich, D. Andelman, and R. Podgornik, Charged membranes: Poisson-Boltzmann theory, DLVO paradigm and beyond, in *Handbook of Lipid Membranes: Molecular, Functional, and Materials Aspects*, edited by C. R. Safinya and J. O. Rädler (CRC Press, Boca Raton, FL, 2021).

[40] P. M. Biesheuvel, S. Lindhoud, R. de Vries, and M. A. Cohen Stuart, Phase behavior of mixtures of oppositely charged nanoparticles: Heterogeneous poisson-boltzmann cell model applied to lysozyme and succinylated lysozyme, *Langmuir* **22**, 1291 (2006).

[41] A. Majee, M. Bier, and S. Dietrich, Electrostatic interaction between colloidal particles trapped at an electrolyte interface, *J. Chem. Phys.* **140**, 164906 (2014).

[42] A. Majee, T. Schmetzer, and M. Bier, Electrostatic interaction between dissimilar colloids at fluid interfaces, *Phys. Rev. E* **97**, 042611 (2018).

[43] T. Curk and E. Luijten, Charge-regulation effects in nanoparticle self-assembly, arXiv:2009.08544.

[44] J. Landsgesell, P. Hebbeker, O. Rud, R. Lunkad, P. Košovan, and C. Holm, Grand-reaction method for simulations of ionization equilibria coupled to ion partitioning, *Macromolecules* **53**, 3007 (2020).

[45] F. Grünewald, P. C. T. Souza, H. Abdizadeh, J. Barnoud, A. H. de Vries, and S. J. Marrink, Titratable Martini model for constant pH simulations, *J. Chem. Phys.* **153**, 024118 (2020).

[46] I. S. Gradshteyn and I. M. Ryzhik, *Table of Integrals, Series, and Products* (Academic, New York, 1980).