Charge symmetry broken complex coacervation

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\"uckel 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.


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 similarly charged macroions can also lead to LLPS, presumably because of short range specific interactions of non-electrostatic 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 liquid-like 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/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 pH, have been scant. A simple cell model approach was used to analyse 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 non-trivial 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].
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 [24] that a pair of chemically identical interacting charge-regulated planar macroions are not necessarily equally charged and that the electric field at the mid-plane of the set-up 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 [25], 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 [24,26,27] -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-desorption, and the fraction η of filled sites on a surface is a degree of freedom within our model. By construction, η ∈ [0, 1]. If the area per site is a 2 , then the charge density is given by σ = e a 2 η − 1 2 with e > 0 being the elementary charge, so that − 1 2 e a 2 ≤ σ ≤ 1 2 e a 2 . The surface number density 1/a 2 of adsorption sites is related to the number K = 4πR 2 0 /a 2 of adsorption sites on a single colloidal macroion. As in Refs. [24,26,27] we base our macroion CR model on the Frumkin-Fowler-Guggenheim isotherm [28] of the macroion surface defined with the phenomenological free energy of a single adsorption site in the units of thermal energy β = 1/k B T as The parameters α and χ are phenomenological and describe the non-electrostatic part of the protonmacroion and the proton-proton interactions at the macroion surface. In the case of (de)protonation reaction, the dependence of α on the bulk pH is model specific [22], but one can explicitly identify α = (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 [29]. Furthermore, χ, as in the related lattice regular solutions theories (e.g., FIG. 1. Macromolecular solution (panel (a)) and the magnified view (part (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 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. the Flory-Huggins theory [30]) describes the short-range interactions between nearest neighbor adsorption sites on the macroion surface [23]. A parameter value α ≤ 0 encodes a favorable adsorption free energy between protons and the macroion surface, while χ ≥ 0 represents the tendency of protons on the macroion surface adsorption sites to phase separate into domains.
In what follows we will use both α as well as χ 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.

A. Formalism
As the configuration of a single colloidal macroion is described by the position r ∈ V of the center of mass and the average degree of protonation η ∈ [0, 1] on its surface, the whole suspension can be described by the number density n(r, η). The equilibrium number density minimizes a grand canonical density functional Ω[n] (see Ref. [33]), which is approximated in the low-density limit by Here, ζ is the fugacity, F ex hc represents the excess free energy due to the hard core interaction between two colloidal macroions and F ex el describes the excess free energy contribution of the electrostatic interaction. In the following, the hard core excess free energy F ex hc is based on the Percus-Yevick (PY) closure and the corresponding equation of state via the compressibility route is used [34,35].
The colloidal macroions are assumed to be suspended in an electrolyte solution with relative permittivity ε r and Debye length 1/κ. For not too highly charged macroions in a sufficiently dilute suspension one can use the Debye-Hückel (DH) approximation [35,36] for the electrostatic two-particle interaction potential where here and below r = |r − r |, while the dimensionless surface charge density of a colloidal macroion with average degree of protonation η and the Bjerrum length B = βe 2 /(4πε 0 ε r ) of the solvent with the vacuum permittivity ε 0 are introduced. Considering the electrostatic interaction U el as a perturbation of the hard core interaction U hc with one obtains in the low-density limit [33] βF ex The considered model is then specified by the following five parameters: α, χ, κR 0 , κ B , K, among which α and χ describe the charge regulation (according to Sec. II). The values for the parameters α and χ are chosen keeping in mind that for χ = −2α the surfaces remain charge neutral for χ < χ c below a certain critical value χ c > 0, whereas they can be oppositely charged for χ > χ c [24]. 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 B ≈ 0.7 nm and Debye length 1/κ ≈ 10 nm, leads to the values κR 0 ≈ 1 and κ B ≈ 0.07. Finally, within the present DFT approach, K ∈ {20, 40, 45, 50} adsorption sites per colloidal macroion are considered. This corresponds to surface areas per adsorption site a 2 = 4πR 2 0 /K ∈ {63, 31, 28, 25} nm 2 , i.e., average distances between neighbouring adsorption sites of a ∈ {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(r, η) = n b (η) and Φ(r) = Φ b . Upon solving the Euler-Lagrange Equation (S4) of the Supplementary Material [37] one obtains the bulk packing fraction profile which provides the distribution of the average degree of protonation η or, equivalently, of the surface charge densities σ * (η) (see Eq. (4)). charged. If the charge regulation parameters α and χ (see Eq. (1)) fulfill the relation χ = −2α the peak is at σ * = 0 (see the blue curve in Fig. 2), whereas for α ≷ −χ/2 the majority of colloidal macroions carry a surface charge σ * ≷ 0 (see the green curve in Fig. 2). Upon increasing the value of the charge-regulation parameter χ the surface charge distribution becomes bimodal with the peaks being located at increasingly 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-density and a low-density phase for α ≈ −χ/2, the case χ = 20 for various numbers K of adsorption sites per colloidal macroion are 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-density 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 characterised by a phenomenological Flory-like 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 the Supplementary Material [37]). It can be conveniently analyzed in terms of the  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. 3} of chargeregulated spherical colloidal macroions of radius R 0 with K = 20 adsorption sites per macroion and charge regulation parameters α = −5, χ = 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 Φ b the mean nearest-neighbor distance λ NN decreases to almost the close-contact distance 2R 0 of two hard spheres of radius R 0 , whereas the periodicity λ ZZ of the charge distribution is slightly larger than twice that distance, λ ZZ 2λ NN . In parallel the coordination number N 1 of colloidal macroions in the nearest-neighbor shell increases.
number-number structure factor which describes the relative distribution of colloidal macroions irrespective of their charge, and the charge-charge structure factor which describes the relative distribution of charge within the fluid. Figure 4 displays the structure factors S NN (q) and S ZZ (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 S NN (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 S ZZ (q) in Fig. 4(b) indicates a spatially alternating arrangement of oppositely charged 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 not occur here. From the position of the major peaks in Figs. 4(a) and 4(b) one obtains the mean nearest-neighbor distance λ NN between two colloidal macroions and the periodicity λ 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 λ NN is never smaller than the close-contact distance 2R 0 . Upon increasing the packing fraction Φ b both λ NN and λ ZZ decrease, but the relation λ ZZ 2λ NN holds for all cases. Hence the charge distribution oscillates on a length scale λ ZZ which is approximately twice the nearest-neighbor distance λ 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) ξ NN and ξ 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 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, the building block 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.
The model medium built from this elementary cell construct thus has a particle-cell volume ratio Φ v = (R 0 /R) 3 which is related to the bulk packing fraction Φ b defined earlier via the relation 74 is the packing fraction corresponding to the face-centered cubic or hexagonal close-packed arrangement of the cells. Clearly, Φ v is inversely proportional to the cube of the inter-macroion 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 β = 1/k B T , can be written as where β Ω s (η) is the energy contribution stemming from the chemical processes driving the (de)protonation reaction (as defined in Eq. (1)) at the macroion surface and with the permittivity ε = ε r ε 0 of the embedding medium is the electrostatic part of the grand potential expressed per unit surface area of the macroion [38]. Herein φ (r) is the dimensionless electrostatic potential expressed in the units of βe which fulfills the PB equation in spherical symmetry, 1 r 2 r 2 φ (r) = κ 2 sinh(φ(r)), subject to the boundary condition of a charge density σ at the macroion surface, i.e., at r = R 0 , and of a given radial component E R of the electric field at the cell boundary, i.e., at r = R, For two coupled cells the grand potential of the binary cell system can then be written as where the subscripts "1" and "2" are used to indicate the two adjacent cells in the binary cell model. The energy contribution β Ω 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 (η 1 , η 2 ) or equivalently (σ * 1 , σ * 2 ). The exact solution of the linear electrostatic problem implies a nonlinear function of the degrees of protonation η 1 , η 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 Supplementary Material of [24], 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 [39]; details are described in the Supplementary Material [37]. After inserting the equilibrium E R that minimizes the βΩ(η 1 , η 2 , E R ) obtained by expanding Eq. (10) up to second order in φ and replacing φ(r) by the known solutions, one finally arrives at the following expression for the grand potential defined in Eq. (13) as a function of η 1 and η 2 only: As before, K = 4πR 2 0 /a 2 , and B is the Bjerrum length. The factors |det M|, γ, τ and ν 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 [37].
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 (α, χ) 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. FIG. 5. Top panel shows the variation of the grand potential βΩ σ * 1 , σ * 2 /K given by Eq. (14) as functions of the surface charge density variables σ * 1 and σ * 2 for different values of the parameters α and Φ v . Bottom panel show cuts through the potential surface at σ * 1 = σ * 2 (in yellow, light shade) and σ * 1 = −σ * 2 (in blue, dark shade). Panel (a) shows the results for α = −3, χ = 6 and Φ v = 0.6 which falls inside the phase separation region. Consequently, the asymmetric configuration is the minimum-energy configuration. Panel (b): At the same value of Φ v , but with α changing from −3 to −2.8. As it is evident from the curves, the stability has changed towards a configuration with uncoupled cells, thus leaving the phase coexistence region. Panel (c): Again at α = −3 and χ = 6 but going down in the particle-cell volume ratio to Φ 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, κ = 0.1 nm −1 , B = 0.7 nm and K = 20 are used.
The phase coexistence around the symmetry axis −2α = χ 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 Φ v , the charge state with minimum energy shifts from a symmetric to an asymmetric one as one moves away from the symmetry axis −2α = χ (compare panels (a) and (b) of Fig. 5). Moreover, a comparison between panels (a) and (c) of Fig. 5 suggests that for given (α, χ), as the particle-cell volume ratio Φ 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 Φ b (or equivalently, the particle-cell volume ratio Φ v ) where the following nonlinear PB theory needs to be applied. In each case, one observes a conical shaped region centered around α = −χ/2 = −15 with stretched opening at the top. Inside this region and on the boundary, the two macroion surfaces are oppositely charged (|η 1 − η 2 | = 1 or equivalently, σ * 1 = −σ * 2 ) whereas outside this region they are identically charged (η 1 = η 2 or equivalently, σ 1 = σ 2 ). With increasing K-value, this region featuring charge-asymmetry widens.

B. Full Poisson-Boltzmann case
Within the nonlinear PB theory, the equilibrium η-values (or equivalently σ * -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 ε r ≈ 80 at temperature T = 300 K.
The resulting variations of the degrees of protonations η 1 and η 2 leading to charge densities σ * 1 and σ * 2 at the two macroion surfaces are shown in Fig. 6 as functions of the interaction parameter α and the particle-cell volume ratio Φ v for two different values of the parameter K ∈ {50, 1257} corresponding to a ∈ {5, 1} nm and χ = 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 Fig. 6, the results are qualitatively similar to those obtained within the linearized PB theory in the previous section and are consistent with the findings reported in Fig. 5. For any given K-value, one obtains a triangularly-shaped region centered around α = −χ/2 with stretched opening at the top. Inside this region and on the boundary, the two macroion surfaces are oppositely charged (σ * 1 ≈ −σ * 2 0) whereas outside this region they are identically charged (σ * 1 = σ * 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 α < χ 2 = −15 whereas for α > χ 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 α-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 α = −χ/2 = −15 are very stable and persist down to Φ v ≈ 10 −3 or lower. The specific choice of the parameter χ = 30 is motivated by the observation that complete symmetry breaking, i.e., σ * 1 ≈ −σ * 2 0, occurs above a critical value χ = χ c ≈ 25 for K = 1257. It is important to note that this χ c -value depends not only on K but can be different within the linear and nonlinear PB theories. However, for any χ > χ c , one observes the same qualitative features, i.e., transitions from symmetric to symmetry-broken states as functions of α and Φ v within both the theories.

V. CONCLUSIONS
In summary, we have described a charge regulation based mechanism that allows for pHdependent 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-like 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 liquidliquid 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 [40], on the grand-reaction method for coarse-grained simulations of acid-base equilibria with a fixed pH reservoir and salt concentration [41], or simulating the pH effects with classical coarse-grained molecular dynamics simulations [42], 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.
with the system volume V = |V|, the constants where Γ(ν, z) denotes the incomplete Γ-function [1], and the bulk packing fraction profile . Using Eqs. (S1) and (S2) in Eq. (2) of the main text one obtains a scaled density functional in terms of the bulk packing fraction profiles Φ b : with the scaled chemical potential The equilibrium bulk packing fraction profile Φ eq b minimizes the scaled density functional Ω * b in Eq. (S3). Hence it solves the Euler-Lagrange equation [2] This expression can be rewritten in the form with the k-th moment of the surface charge distribution B. Partial structure factor The partial structure factor of the system under consideration can be written in the form where G(q, η, η ) is the three-dimensional Fourier transform of the bulk correlation function G(r, η, η ) of the number densities of colloidal spheres with degrees of protonation η and η at a distance r. The auxiliary function G(q, η, η ) fulfills the Ornstein-Zernike equation with C(q, η, η ) = n b (η)n b (η ) c(q, η, η ), where c(q, η, η ), q = |q|, is the Fourier transform of the bulk direct correlation function c(r, η, η ) = c PY (r, Φ b ) − δ 2 βF ex el δn(r, η)δn(0, η ) n b composed of the Percus-Yevick hard-core contribution c PY (see Ref. [3]) and the electrostatic contribution obtained from Eq. (5) of the main text (see Ref. [2]).
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 [4]: The boundary values φ(R 0 ) and φ(R) thus follow from the derivatives of φ(r) at r = R 0 , R.
In the linear theory we discuss, denoting the vector φ = (φ (R 0 ), φ (R)) of derivatives, one needs to invert the matrix equation φ = M · φ which can be computed from the derivative of the solution. 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): where C(R) = cosh(κ(R − R 0 )) and S (R) = sinh(κ(R − R 0 )). Further, det M = νγ − τξ .
One has ν < 0, τ > 0, ξ < 0, γ > 0 and det M < 0. After computation of the boundary conditions at the macroion radius R 0 and the binary cell radius R, the resulting expression Eq. (S5) needs to be minimized with respect to E R , which is found to behave as E R ∼ σ * 1 − σ * 2 . Finally, collection of terms leads to Eq. (14) of the main paper.
Within DH-theory, the surface-charge density coexistence curve has a balloon-like shape. The upper part of the surface-charge density coexistence curve widens for increasing χ. The location of the upper critical point of the phase coexistence curve moves from (χ, Φ v ) = (3.6, 0.79) to (χ, Φ v ) = (16, 0.97), which for K = 20 covers the interval in which surface-charge density coexistence exists in the DH-limit.