Role of electrostatic correlations in polyelectrolyte charge association

Polyelectrolytes have been extensively investigated in the past and continue to attract attention owing to both practical and fundamental interest.¹ On one hand, they are relevant to numerous applications, including solution property modifiers,^2,3 food processing technologies,⁴ pharmaceutical and therapeutic encapsulation vehicles,^5,6 and tissue engineering scaffolds.⁷ On the other hand, many fundamental questions in such systems remain poorly understood.⁸ Earlier scaling analyses^9,10 and field theoretical approaches¹¹ provide a basis for qualitatively rationalizing certain behaviors of polyelectrolyte solutions. Yet quantitative agreement between experimental results and theories regarding the structural, thermodynamic, and rheological behavior of these systems is scarce, unlike that for neutral polymeric systems, implying, to say the least, that our picture is incomplete.¹ 

Much recent activity in this area has been spurred by a resurgence of interest in the study of polyelectrolyte complex coacervation.¹² Coacervation refers to a demixing instability in polyelectrolyte solutions which results in a polymer-rich coacervate phase and a dilute supernatant phase. It is driven largely by associative electrostatic attractions among oppositely charged polyelectrolyte species. Complex coacervation is, in particular, ubiquitous within biology and has been shown to play a crucial role in cellular compartmentalization¹³ and the aggregation of intrinsically disordered peptides responsible for the onset of various neurodegenerative diseases.¹⁴ A series of experimental and theoretical studies¹⁵ have demonstrated that coacervation can be understood by considering the competition between the mixing entropy and various associative forces.¹² Given our comparatively poor understanding of simple polyelectrolyte solutions, it is in fact surprising that the presumably more complicated coacervation process can be modeled reasonably well. Therefore, many theoretical efforts have been devoted to this system, such as incorporating the effects of liquid-state correlations,^16,17 ion binding equilibria,^18,19 charge connectivity,²⁰ and adaptive chain structures,^21,22 in hopes of truly reconciling theory and experiments.

The focus of this work is ion binding equilibria. A basic question in polyelectrolyte solutions is understanding the long-range correlations of ionic species, including both the ionic polymers and their counterions.¹ On top of this, there exist a variety of short-range correlations that strongly affect polyelectrolyte solution thermodynamics. Ionized groups on oppositely charged chains can form specific ionic cross-links or ion-pairs, which contribute to the viscoelastic nature of coacervates.²³ The fraction of cross-links or bound pairs depends sensitively on the concentration of small ions through a doping-like equilibrium.^24,25 In addition, the activity of counterions in polyelectrolyte solutions deviates strongly from that of simple electrolyte solutions due to strong localization of counterions in the vicinity of the polymer chains.²⁶ The classical theory of Manning²⁷ posits that this localization helps to alleviate the strong, mutual repulsion among ionized groups on the backbones of polyelectrolytes. Counterion binding regulates the effective charge of polyelectrolytes and reduces the mobility of counterions, allowing this effect to be probed experimentally through a variety of conductometric and spectroscopic techniques.^28–33 

A multitude of theoretical models have been applied to describe the effects of ion binding or charge association. Poisson-Boltzmann approaches rely on a continuum description of the counterion charge density around a chain.^34–36 A transfer matrix formalism has been developed which describes the adsorption states along a test chain through a series of transition probabilities weighted by appropriate Boltzmann factors.^37,38 Beyond these, a general framework is to describe ion association through a set of coupled quasi-chemical reactions, featuring a balance of entropic effects favoring unbinding and energetic factors driving association.¹⁹ 

Given a form of solution free energy which includes, in general, the translational entropy of all species, entropic contributions from permuting bound and unbound states along the chains, and all excess interactions, the equilibrium extent of association is determined self-consistently by minimizing the total free energy with respect to the fraction of bound charges. The outcome is effectively a set of laws of mass action relating the bound fraction of ions to energetic considerations within the model, which have the following general form:

Here, K is an equilibrium constant and the temperature k_BT has been implicitly set to unity. The term μ is an exchange chemical potential derived from variations of the excess free energy. Contributions to this potential may arise from solvation effects,¹⁹ dispersion interactions,¹⁹ dipolar interactions between bound charges,^18,39,40 or electrostatic correlations between dissociated ionic species.^{18,19,41–44} The equilibrium constant depends as well on the prescribed binding free energy between oppositely charged species, ΔG, which may represent any local interactions not included in the excess free energy. Consequently, the expected association behavior depends strongly on the molecular interactions present in the system, which in turn depend on the concentrations, energetic parameters, and other solution properties.

The minimum force needed in the exchange chemical potential is the electrostatic correlation free energy. This contribution stems from the propensity for ionic species to preferentially arrange themselves around oppositely charged ions due to favorable electrostatic interactions. Mean field models cannot describe this effect in homogeneous, bulk solutions where the average charge density vanishes, and consequently fluctuation effects must be considered. A simple model to capture these effects, which has been applied previously,¹⁹ is the Debye-Hückel (DH) free energy.⁴⁵ DH theory applies to a constellation of disconnected point-like charges, which says nothing as to the physical connectivity of ionized groups, or the architecture, of the polyelectrolytes. Generalization of DH theory to polyelectrolytes has been worked out,^11,42,46 and the consequences on the thermodynamics of non-associating fluids investigated.²⁰ In those studies, charge connectivity was accounted for through specific reference to the intramolecular correlation function, i.e., structure factor, of the chain. More sophisticated treatments have allowed the structure factor (usually, the chain stiffness) to be adaptive and self-consistently determined.^21,22

Our work aims at examining the effects of collective charge correlations and chain connectivity on ion binding equilibria. This is studied within the context of counterion binding in a single polyelectrolyte solution, which we treat conceptually as a reversible adsorption process through the law of mass action approach described above. The effective association constant in this process receives a contribution from the correlation free energy, which describes the effects of long-range charge density fluctuations. Charge connectivity and the architecture of the polyelectrolyte are captured in the correlation free energy by modifying a previous theory developed for non-associating solutions.²⁰ Particular attention is paid to the fractal dimension and backbone stiffness of polyelectrolytes. In contrast to previous studies, we find that a consistent treatment of charge association behavior must incorporate the self-energy of dissociated ions in the correlation free energy with caution. In particular, we show that the standard DH expression for the correlation free energy cannot be directly applied to study reversible association behavior, and an alternative model is proposed.

This writing is organized as follows. In Sec. II, we introduce a general thermodynamic model for treating reversible counterion binding in polyelectrolyte solutions. In Sec. III, we demonstrate the unphysical behavior related to the application of the standard DH free energy to describe electrostatic correlations and show how it can be cured by smearing charges to account for the ion self-energy.⁴⁷ Section IV discusses how chain connectivity impacts charge association behavior. Various chain structures are considered, allowing us to reveal the effects of molecular length and compactness on ion binding. Finally, Sec. V summarizes our findings and implications for future modeling.

We consider an aqueous solution of polyanion chains, counterions, and coions. For simplicity, we assume all monomers are monovalent with identical molecular volume $v$⁠. The volume fractions of polyanion, cation, anion, and solvent are denoted as ϕ_A, ϕ₊, ϕ₋, and ϕ_W, respectively. In the absence of added salt and if polyanions are fully dissociated, overall charge neutrality demands the counterion fraction be equal to ϕ_A. If a total fraction of 2ϕ_S of salt is added, the cation and anion fractions become ϕ₊ = ϕ_A + ϕ_S and ϕ₋ = ϕ_S. At any given composition, however, the polyanion is only partially charged due to binding of small cations along the chain. Out of n_A polyanions of degree of polymerization N_A, N_A+ monomer units harbor a bound cation resulting in an associated fraction of α = N_A+/(n_AN_A) and an effective linear charge density along the chain of σ = 1 − α. Consequently, the fraction of free cations is $ϕ+f=ϕ+−αϕA=ϕS+σϕA$⁠, and the fraction of bound ions is $ϕ+b=αϕA$⁠.

The basis for our analysis of counterion association behavior is an expression for the solution free energy that depends on the degree of charge association α explicitly, which is to be minimized to yield the equilibrium fraction of bound ions. The free energy contains four additive contributions,

related to the translational entropy of free species, the entropy of permuting bound cations along the polyanion backbone, a binding energy for cation-polyanion pair formation, and electrostatic correlations between unpaired ionic species. Above, the free energy has been normalized by the reference (monomer) volume $v$ and the thermal energy k_BT.

The first term f_T comes from the mixing entropy and includes explicitly the contribution from solvent,

where the summation index γ runs over unbound species. We note that the term for cations containing a factor of σϕ_A imparts an α dependence to f_T, which proves important when minimizing the free energy to find equilibrium association extents.

The second term f_C is the entropic contribution due to permutations of bound ions along the polyanion backbone. For a system of n_A polyanions with degree of polymerization N_A, the total number of ways to arrange N_A+ bound sites on n_AN_A total monomer units is (n_AN_A)!/[(n_AN_A − N_A+)!(N_A+)!], resulting in a free energy contribution of the form,

This expression neglects any potential cooperativity effects due to chain connectivity and is effectively a mean-field treatment of ion adsorption along the chain.^18,19,42

The next term f_ad includes the binding energy ΔG_α for counterion-polyelectrolyte pair formation, scaled by the fraction of bound ions and polyanion composition. The resulting free energy contribution is

The binding energy captures any short-range, nonlinear interactions that contribute to the effective association constant for counterion localization along the chains.^18,19,44,48 In prior studies, this term has been correlated with a local dielectric mismatch and treated as a function of the relative permittivity.^18,39 Additionally, depending on the accuracy with which molecular interactions are captured by the excess free energy, ΔG_α may depend on factors such as molecular weight and specific chemical functionality. While aware of the importance of solvation effects and chemospecific interactions, we refrain from proposing a specific form for the dependence of the binding energy on these contributions. For the purpose of this work, it suffices to treat ΔG_α as a measure of the free energy gain per counterion bound along the backbone in a mean-field fashion. This parameter is prescribed a constant value of order k_BT and is varied in later calculations to illustrate its effect on the association equilibrium.

The last term f_corr contains the contribution of all fluctuation effects in the system. In our case where electrostatic interactions dominate, this term represents the free energy contribution due to long-range charge correlations and is treated with variants of the Debye-Hückel free energy or models derived through a field-theoretic Gaussian fluctuation analysis.²⁰ A general formalism for our Gaussian theory approach can be found in the  Appendix. The association constant for ion binding ultimately receives contribution related to this free energy term, in addition to the local binding energy mentioned above. For now, these various models are left unspecified and examined independently in detail in Secs. III and IV.

Minimizing the solution free energy with respect to α results in a law of mass action relating the fraction of bound ions to the binding energy (ΔG_α) and differentials of the correlation free energy (f_corr). For our system, the resulting expression is

Here K_α is the association constant for reversible ion binding. We note that even in the absence of a short-range interaction energy (ΔG_α = 0), K_α may be large depending on the relative sign and magnitude of the exchange potential μ_α, which in turn may yield a significant degree of ion association. The exchange chemical potential $μα=1ϕA∂fcorr∂α$ depends on the specific form of correlation free energy and represents the driving force for ion association due to (long-range) electrostatic fluctuations in the system. The form of μ_α is dictated by the structure of the polyelectrolytes. The focus of this work involves examining equilibrium predictions for the extent of binding α with various models for μ_α. These predictions are obtained by simply solving Eq. (6) as a function of composition or other energetic parameters, e.g., ΔG_α.

In the absence of any correlation effects (μ_α = 0), and assuming a constant binding energy, the right-hand side of Eq. (6) can be solved for α analytically and evaluated as a function of concentration. This simplified binding behavior is portrayed in Fig. 1, where we see that α is a monotonically increasing function of both polymer and salt concentrations with a rate dependent on the strength of binding energy. In this regime, the bound fraction ultimately approaches a constant depending on the value of ΔG_α. This trend can be interpreted as an entropic driving force for charge association solely due to the greater number of free charges available for binding at higher concentrations. On the other hand, in the limit of low ion concentration, the law of mass action without correlation effects predicts no bound charges, owing to the increased translational entropy of free ions in this regime. Deviations from this baseline behavior at both low and high concentrations due to correlation effects are of particular interest in Secs. III and IV.

The Debye-Hückel free energy has historically been used to describe charge correlations in polyelectrolyte systems, largely owing to its simplicity and ability to capture several qualitative features of experimental polyelectrolyte phase behavior.^15,49,50 This model captures the effect of charge density fluctuations due to favorable electrostatic interactions between oppositely charged species for a constellation of disconnected point-like ions, neglecting any notion of chain structure and connectivity.²⁰ The DH theory has been applied previously in charge association models to treat ion interactions for a background of dissociated point-like ions,^18,39 as well as for both chains and ions.¹⁹ In the analysis below, we adopt this latter approach and apply the DH model to treat both polyelectrolyte and salt correlations within our general charge association model. We show that the straightforward application of the DH free energy leads to unphysical binding behavior at high concentrations. This behavior is caused by an inconsistent treatment of the ion self-energy and can be cured by incorporating the self-energy explicitly in the correlation free energy using an approach developed in the literature for a different problem.⁴⁷ From this, we gain insight into a physical interpretation of the electrostatic correlation effect on ion binding. An analysis of chain structure and connectivity effects within the correlation free energy is deferred until Sec. IV.

We begin by examining the concentration dependence of the degree of ion association calculated using the classical Debye-Hückel theory. For our model system, DH theory yields the following form of the correlation free energy:

where the non-dimensional inverse Debye length is given by $κ¯2=4πl¯B∑γσγϕγ=4πl¯B(2ϕS+2σϕA)$ and the factor of two in each term stems from the free cation fraction $ϕ+f=σϕA+ϕS$⁠. The term $l¯B=lB/l$ is a reduced Bjerrum length, l = $v1/3$ is the cube root of the molecular volume, and l_B = e²/(4πϵk_BT) is the Bjerrum length. For the remainder of this work, we set the reference volume to $v$ = 3 × 10⁻²⁹ m³, corresponding to the volume of a water molecule. Then using the dielectric permittivity of water at room temperature, ∼80, gives a dimensionless Bjerrum length of $l¯B≈2.3$⁠. With this convention, the exchange potential is

Figure 2 shows the equilibrium values of α as a function of the concentrations of polymer and added salt, for a set of association free energies ΔG_α. In Fig. 2(a), the solution is salt-free, so ϕ₊ increases in proportion to ϕ_A to maintain charge balance. In the dilute regime, the association degree increases with ion concentration. In the more concentrated regime, surprisingly, the model predicts a vanishing α, i.e., completely dissociated ions. This effect is more pronounced when the binding energy is small, to the extent that no association is predicted when ΔG_α = 0.

Although this behavior at high concentration is not strictly prohibited from a thermodynamic perspective, a vanishing bound fraction at high concentrations is incompatible with a multitude of reports available in the literature. Simulation studies of polyion solutions have highlighted an increase (decrease) in the number of condensed counterions (effective charge of chains) with increasing polymer or salt concentration, lending notion to the idea that counterion binding persists through a reversible, concentration-dependent equilibrium.^51,52 Additionally, although we do not expect our simple model to capture experimental trends with quantitative accuracy, a variety of studies relying on conductiometric and spectroscopic measurements of polyion solutions also support the finding of a significant bound fraction in the concentrated regime.^28,30,32,33 We carry on in Sec. III B to show that the discrepancies between these prior studies and our model result from a thermodynamically inadmissible treatment of the electrostatic self-energy implicit within the DH free energy, and then propose a method to remedy this error.

We show here that the non-intuitive behavior observed above originates from an improper treatment of the ion self-energy for systems where association is considered. We first present an alternative derivation of the DH free energy that highlights the nature of this error, then propose a modified DH free energy to correct this point, which incorporates the self-energy explicitly by endowing point charges with a finite width.

It is well-known^20,21,42 that the classical DH free energy can be obtained from the perturbative Gaussian fluctuation or random phase approximation (RPA) approach for a system of electrostatically interacting point charges. The correlation free energy, in this case, adopts the form (the  Appendix),

This expression is presented in reduced units, where the wavevector q is normalized by the monomer length l. The integral over wavenumber exhibits an ultraviolet (UV) divergence due to insufficient decay of short-range interaction modes. This divergence is an artifact of the delta-function description of charged species commonly employed in field-theoretic models and would not occur in a more realistic model which accounts for the finite size of ions and molecules.^53,54 A common practice to regulate the theory is to subtract off the diverging contribution $κ¯2/q2$ from the integrand, which is associated with infinite, zero-point self-interactions through a Coulomb potential.^20,46,47 We refer to this term as an effective self-energy for the ions. In many cases, this subtraction is thermodynamically admissible due to the linear dependence of $κ¯2$ on composition and physically amounts to choosing a plasma as the reference state for the electrostatic correlation free energy.²⁰ In the point-like example considered above, this yields

which is identical to Eq. (7) for the Debye-Hückel free energy.

However, this form of regularization is no longer valid when reversible charge association is considered. In such cases, the self-energy is a function of the total charge density, which in turn depends on the bound fraction of ions α. After the free energy is minimized against α, σ and $ϕ+f$ appearing in $κ¯2$ become highly nonlinear functions of composition, which results in a non-constant shift in the reference chemical potential of the system. This thermodynamically inconsistent treatment manifests itself as the non-intuitive trend observed in Fig. 2. This type of phenomenon has been mentioned in a previous study on associating polyelectrolyte networks,⁴³ though no connection has been made to the widely used form of the Debye-Hückel free energy. We emphasize that this subtle observation regarding the DH free energy applies not only to treatments of charge association but rather to any system with an inhomogeneous interaction strength, such as ions in a heterogeneous dielectric medium.⁴⁷ 

We propose a remedy for this issue in accordance with more recent literature practices^{22,47,55–57} by smearing charges such that each ion acquires a finite width. This procedure damps the short wavelength modes sufficiently rapidly and regularizes the free energy. The implementation of this method requires typical transformations in derivation of the correlation free energy and is presented in the  Appendix. The resulting smeared expression for the free energy formally resembles that of standard DH theory,

where $κ̃2(q)=4πl¯BΓ^2(q)∑γσγϕγ$ is the now wavenumber-dependent inverse Debye length. The smearing function $Γ^(q)$ appearing above is defined as $Γ^(q)=exp(−12q2a2)$⁠, a Gaussian distribution in Fourier space with characteristic width a. The point-charge model is recovered when the smearing width vanishes. This width is expected to be close to the size of ions, and its value can be fixed by comparing the predicted ion solvation energies to experimental values. Quantitative predictions from our model depend on the choice of this smearing length—larger values result in a progressive damping of the correlation free energy and thus the electrostatic correlation effect on ion binding. To simplify the discussion, in the following, we set the smearing width to the monomeric radius, $a=12l$⁠.

To demonstrate the effects of proper handling of the self-energy on charge association behavior, we focus on the same single-chain model system as above and derive an exchange chemical potential for ion binding from the smeared DH free energy. Recognizing that $κ̃2(q)=4πl¯BΓ^2(q)(2σϕA+2ϕS)$⁠, the resulting potential to be inserted into Eq. (6) is

As seen in Fig. 3, this modified theory produces an association fraction that increases with ion concentration and approaches a non-zero value that depends on the choice of binding energy in the concentrated regime. Additionally, a sizable increase in α is observed when ΔG_α = 0, indicating a favorable driving force for binding arising from the electrostatic correlation exchange potential. We discuss the physical origin of this change in behavior in Sec. III C.

The results of Figs. 2 and 3 demonstrate that including the self-energy of ions within the correlation free energy significantly alters the concentration dependence of ion binding. To further elucidate the origin of this contrasting behavior, we consider a third variant of the DH free energy that has been applied in previous treatments of charge association¹⁹ and that does not exhibit the same high concentration behavior seen in Fig. 2. The correlation free energy from this model, often referred to as an extended Debye-Hückel (EDH) free energy, is given by^19,45,58

While the extended model is equivalent to the standard DH model only in the limit of infinite dilution (⁠$κ¯→0$⁠), it contains the same negative prefactor as in the standard DH model [Eq. (7)]. The resulting exchange chemical potential is

A comparison of the three exchange potentials derived from variants of the DH free energy can be found in Fig. 4. The DH and the extended DH models result in a positive exchange chemical potential, while the opposite is true for the smeared DH model. This difference stems from how the correlation free energy is regularized; subtracting off the nonlinear self-energy term results in the sign change in both the free energy and the exchange potentials for the DH and EDH models. This proves to be extremely influential for the overall binding equilibrium, as the correlation exchange potentials contribute alongside the local binding energy to the effective association constant in our model.

As inferred from Eq. (6), the positive exchange potentials derived after improper self-energy removal always oppose ion association, while the negative smeared potential always favors binding. At high concentrations, the DH and EDH potentials become large in magnitude and can dominate the binding equilibrium, resulting in the diminished bound fraction as seen in Fig. 2. From these observations, we interpret $μαDH$ and $μαEDH$ as a driving force for ionic dissociation due to favorable electrostatic correlations in the absence of an ion self-energy, while $μαSDH$ represents a driving force for ion binding due to a reduced self-energy as oppositely charged complexation occurs.

This observation is supported by recent work on the diverging “self” and attractive “correlation” contributions to the RPA free energy.⁵⁹ The RPA result for the electrostatic free energy contains the energetic cost of assembling charge from an infinitely disbursed reference state into a finite system volume and distribution of charged species. In addition, it contains a correlation contribution due to fluctuation-induced attraction of the already-assembled charge distribution,⁵⁹ which gives rise to favorable cohesive interactions. When the charge distribution is constant, this self-contribution can be safely removed from the RPA free energy. However, as shown above, removal of this nonlinear self-energy contribution is thermodynamically unjustified and results in unphysical behavior when considering reversible binding. We thus conclude that the self-energy contribution to the RPA free energy is essential when considering reversible charge association. Consequently, the exchange potential derived from the smeared charge model physically represents the driving force for binding due to a favorable reduction in the correlation free energy, which is in turn due to alleviation of repulsive self-interactions of dissociated charged species. As the electrostatic self-energy is dependent on the spatial distribution of assembled charge, we expect chain architecture to strongly influence this contribution,^22,59 and we explore these effects further in Sec. IV.

These physical insights related to the ion self-energy hold true for both the standard and extended Debye-Hückel models. However, the reduced magnitude of the exchange potential derived from the extended Debye-Hückel model (Fig. 4), especially at high concentrations, indicates that this expression is a lesser contributor to the overall binding equilibrium. This subtle difference may explain how the behavior shown in Fig. 2 has been overlooked in previous studies.¹⁹ Though, because the association trends with both the classical and extended DH models can result in non-sensible high concentration binding behavior, it is apparent that a more systematic treatment of electrostatic correlations and regularization of the ion self-energy are required to model systems with an inhomogeneous interaction strength. Finally, we stress that since all free energy models discussed above exhibit negative curvature as a function of composition, their effects on thermodynamic miscibility are similar when a treatment of ion association equilibrium is absent. For instance, previously reported analyses of coacervation behavior based on the Voorn-Overbeek model are not affected by such effects.¹⁵ 

The above analysis has highlighted the subtle nature of the ion self-energy in the electrostatic correlation free energy for molecules with no internal structure. In reality, the connectivity of charges in polyelectrolyte chains strongly affects the length-scales and strength of charge correlations. Here we extend our description of charge association to chains where connectivity is accounted for. We present a general expression for the electrostatic correlation free energy for molecules with an arbitrary macromolecular structure. We study two limiting chain structures, fixed Gaussian coils and rigid-rods, and examine how chain length and compactness in these cases alter binding behavior. We then generalize the role of molecular compactness by considering molecules with an arbitrary fractal dimension and chains with a concentration dependent, adaptive persistence length. Finally, we consider how the addition of salt affects binding behavior in both the low and high polymer concentration regimes.

We consider the same system as in Sec. II A, consisting of polyelectrolytes and salt ions in an aqueous solution. As before, the bound ions are treated as effectively smeared along the chains, resulting in an effective linear charge density of σ = 1 − α. Consequently, direct correlations between bound ions and monomer units on chains do not contribute to the free energy.⁴² As such, the effect of ion association manifests in our model on the charge density σ and free cation fraction $ϕ+f$⁠. The  Appendix shows that the electrostatic correlation free energy density adopts the form

where the inverse screening length accounting for chain structure is given by

As in Eq. (11), we rely here on the smeared charge model for both salt and polymers to regularize the divergent high-q contributions to the correlation free energy, without removing a nonlinear self-interaction term. Additionally, we note that this expression for the correlation free energy only considers the effects of charge density fluctuations, and not mass density fluctuations which are present even in neutral polymer systems.²⁰ 

The form of the screening length is similar to the Debye length, except for the dependence on the chain structure factor g(q). For the salt species, it is simply unity, but for the polyelectrolytes it imparts a dependence on the chain structure and bead connectivity. For Gaussian coils, this is given by the Debye function $gD(x)=2x2(x−1+e−x)$ with $x=q2Rg2$ and $Rg2$ being the radius of gyration of the polymer.²⁰ For molecules with rod-like configurations, this is given by the Neugebauer function $gN(x)=2x∫0xdtsin(t)t−1−cos(x)x$ with x = qL and L being the linear length of the rod.²⁰ In reality, however, we expect the polyelectrolyte to adopt a structure somewhere between these two limiting cases (the details depending on concentration). For this reason, we focus our analysis in Secs. IV B and IV C on the effects of chain stiffness and compactness on charge association behavior.

Integration of the above correlation free energy into our association model is straightforward. The exchange chemical potential derived from Eq. (15) is

where g(q) is to be substituted for the appropriate structure factor for a given chain structure. Similar to the case of the smeared DH model, this potential contains a negative prefactor indicating a favorable driving force for ion association. However, unlike Eq. (12), the numerator in Eq. (17) contains two terms, one stemming from the α-dependence of the free counterion fraction and another from the charge density along the chains. This second term scales with N_A and the q-dependent structure factor, implying that the strength of this driving force for ion binding is sensitive to chain length and architecture. As it stands, Eq. (17) can be inserted into the law of mass action Eq. (6) to predict equilibrium values of α.

We consider the effects of molecular architecture, chain length, and backbone stiffness on charge association. Molecular compactness is found to be important for enhancing the degree of ion binding, particularly in the dilute regime. The effects of compactness are shown to arise in both the large scale structure of the chain, as characterized by the chain structure factor, as well as the local chain segment length or stiffness. An intriguing non-monotonic concentration dependence of the association degree is also uncovered in the semidilute regime, which is reminiscent of experimental results in similar systems.^28,30,32

To highlight the effects of electrostatic correlations, we examine salt-free solutions consisting of polyanions and counterions and set the binding energy ΔG_α to zero. We again set the monomer volume to $v$⁠. The degree of polymerization of the polyanion is calculated by N_A = V_A/$v$⁠, where V_A is the molecular volume. For coil-like molecules, the chain size is parameterized by its radius of gyration $Rg2=NAb2/6$⁠, with b being the statistical segment length. For rigid rods, the contour length is specified as L = N_Ab, with an effective segment length b scaling inversely with the cross-sectional area. In both cases, b characterizes the relative molecular compactness of the chain, and we take this parameter in calculations to be of order the reference length l.

The exchange chemical potentials for rod-like and coil-like molecules are plotted in Fig. 5. For both chain structures, their magnitudes decrease with polymer concentration, implying stronger driving force for ion binding in the dilute regime. Furthermore, they increase with chain length, reaching an asymptotic limit at infinite chain length, indicating a significantly more favorable tendency to bind counterions for longer chains. This chain length dependence derives from the factor N_A in the numerator of Eq. (17). Physically, this dependence relates to the energetic cost of assembling like-charged species along a connected chain. Longer chains face a greater self-energy penalty, which encourages the association of oppositely charged ions.

At any given chain length, a significantly larger increase in magnitude of μ_α is observed with decreasing polymer concentration for the more compact coil-like molecules compared to rigid rod-like chains, suggesting a greater propensity for ion binding, particularly in the dilute regime, for this structure. This stark difference in magnitude of the exchange potential is attributed to the strength of self-interactions,²² which are characterized by the choice of structure factor for the chain. Coil-like chains are more strongly correlated with their like-charged neighbors as compared to rigid rod-like chains. Hence, this more compact arrangement of charge results in a greater driving force for oppositely charged ions to bind.

The degrees of ion association for these two structures evaluated using the law of mass action are shown in Fig. 6. The effect of chain length is weak in the more concentrated regime, as is expected. Since the screening length decreases with concentration, beyond a certain concentration, it falls below the molecular size, making the latter an irrelevant parameter. However, in the dilute regime, a strong dependence on chain length is seen for both rod- and coil-like molecules. Longer chains result in greater limiting values of α in the highly dilute regime, reflecting the trends in the exchange potential observed in Fig. 5.

The behavior of rod- and coil-like molecules is also noticeably different. For short rods, α ultimately decays to zero at very low concentrations owing to the increased translational entropy of counterions in this regime. Yet larger N_A, the driving force for ion binding due to correlation effects is enhanced to the point that a plateau value of α ≈ 0.34 persists to highly dilute concentrations. A similar type of dilute behavior is observed for coil-like molecules, though the limiting plateau values of α are much higher and approaching unity (i.e., the fully bound state) at long chain lengths. This trend implies stronger counterion binding for more compact coil-like molecules than rod-like ones. A systematic study of molecular compactness, as characterized by the structure factor of the chain, is presented in Sec. IV C.

The binding behavior is also sensitive to the monomer compactness, which is quantified in our model by the segment length b. Plots of the degree of ion association as a function of segment length are shown in Fig. 7. In the case of rod-like molecules, smaller values of b lead to larger limiting plateau values of α in the dilute regime. Conversely, larger values of b lead to vanishing values of α in the dilute regime. Thus, decreasing b at fixed N_A leads to similar effects in α as increasing N_A at fixed b. One way to understand this dependence on b is to consider the charge per unit length of the rod. The linear charge density in our model is defined as the fraction of charged polyelectrolyte monomer units, which for a fixed N_A does not change with the segment length b. However, decreasing b reduces the length of the rod and increases the total valence per length. This leads to a more compact arrangement of charged monomers and enhances the tendency for ions to bind.

A similar trend is observed for coil-like molecules, but the limiting values of α in the dilute regime are generally larger than the rod-like case for identical b values. Moreover, for larger segment lengths, a large dip in the association fraction is observed in the semidilute regime (around ϕ_A ≃ 0.01) before increasing again in the concentrated regime. This dip is more pronounced for larger statistical segment lengths.

Before proceeding to discuss the origin of the dipping and flattening behaviors in the semidilute regime, we note that the qualitative dependence of α on segment length b remains unchanged even if the smearing width a is allowed to vary in conjunction with b. Additionally, by varying the volume of the polyelectrolyte monomers, it is found that bulkier monomers tend to yield lower values of α compared to their more compact counterparts.⁶⁰ 

Flattening and dipping behaviors in the semidilute regime are evident in both Figs. 6 and 7 and are most pronounced in Fig. 7(b) where the dip in α with the largest segment length spans multiple decades in concentration. This effect derives from a trade-off between the translational entropy of counterions, which opposes binding as the concentration decreases, and the correlation exchange potential μ_α, which favors binding. In the more concentrated regime, ion binding is promoted by the higher counterion concentrations. In the dilute regime, the two effects are comparable and the equilibrium fraction of bound ions is sensitive to the choice of structure factor, chain length, and backbone stiffness. A non-monotonic crossover in α occurs from a balance between these competing factors at intermediate concentrations.

The non-monotonic behavior observed above is reminiscent of experimental trends regarding the effective valence of polyelectrolytes, obtained using conductiometric and dielectric spectroscopy measurements.^28,30 Both non-monotonic trends^28,30 and a plateau³² in charge density have been reported at semidilute concentrations. Carefully calibrated comparisons would be needed to convincingly show that those observations are caused by the charge connectivity effects noted above. Nonetheless, our analysis provides a clue and suggests that the effects of molecular architecture need to be thoroughly considered.

It is natural to generalize the above analyses to molecules with arbitrary fractal dimension. A previous study²⁰ has shown how the structure factor appropriate for fractal molecules can be combined with Eq. (15) to evaluate the correlation free energy. The intramolecular correlation function in q-space for an arbitrary fractal molecule scales as g(q) ∼ (ql)^−d, where l = $v1/3$ is the molecular diameter. The derived exchange chemical potential used in the above analysis enables us to investigate the effects of fractal dimension on ion binding with appropriate substitution for the fractal structure factor.

Figure 8(a) shows that the binding behavior for fractal molecules falls between the rod- and coil-like limits. For fractal dimensions below that of a rod (d < 1), complete counterion dissociation in the dilute regime is observed. However, in the range of 1 ≤ d ≤ 2, the behavior interpolates those of the rod- and coil-like limits. This trend is mirrored in Fig. 8(b), where we plot α versus fractal dimension at a various polymer concentrations. A strong dependence on fractal dimension is clearly seen for d > 1 in the dilute regime. A much weaker dependence is observed in the concentrated regime, as the driving force from the electrostatic correlation energy is relatively small in this regime compared to the effects of ion entropy and the binding energy ΔG_α.

The binding behavior explored here reinforces the notion of molecular compactness highlighted above. The intramolecular correlation functions contain information regarding the probability of finding two monomers on the same chain neighboring one another, which is higher for molecules with higher fractal dimension. With this in mind, we may interpret the driving force for ion binding due to charge correlations as a propensity for molecules to reduce their electrostatic self-energy. More compact, higher fractal dimension molecules feature a greater number of like-charged neighbors surrounding them, and a comparatively greater electrostatic self-repulsion. Counterion association reduces the effective charge density of the polyelectrolyte, and thus its harsh self-interaction. This effect is decidedly more pronounced in the dilute regime due to reduced charge screening.

The results thus far have relied on fixed-structure Gaussian fluctuation free energies to describe the charge correlation effects. However, in realistic polyelectrolyte solutions, the chain structure adapts to changes in concentration by adopting more extended conformations in the dilute regime and more collapsed conformations at higher concentrations. Assuming a fixed coil-like structure has been shown to overestimate charge correlations at low concentrations and produce erroneous results for the dilute branch of coacervate phase diagrams.^22,57 The importance of chain expansion and collapse in counterion binding have also been highlighted in a variety of simulation⁶¹ and theoretical^18,40 studies.

We demonstrate in the following that this adaptive behavior does not change the qualitative picture discussed above. Previous studies on adaptive structures have normally assigned a concentration-dependence to the persistence length l_eff of the chain.^22,43 This has been approximated as the Debye screening length,⁴³ analogous to the classical Odijk-Skolnik-Fixman limit.^1,62 More recently, a variational approach was developed²² that determines this persistence length by minimizing an interpolated form of the single-chain free energy for worm-like chains. We adopt the approach of Ref. 22 here, which entails simultaneously minimizing this single-chain free energy and evaluating the law of mass action in our model to self-consistently determine l_eff and α. This allows for the chains to collapse or stiffen in response to changes in concentration, or the effective charge density as counterion binding occurs.

The binding behavior obtained from the use of the adaptive structure predictions of Ref. 22 is shown in Fig. 9. This result closely resembles that of the rod-like behavior shown in Fig. 6(b). The value of α decays to zero for short chains and approaches a finite constant value for longer ones. This suggests that the persistence length is sufficiently high in this dilute regime, as a result of poor screening. Thus, in the low concentration regime where correlation effects tend to dominate the binding equilibrium, the worm-like chains behave qualitatively as highly extended rods.

Our analysis has thus far focused on salt-free solutions. Adding salt increases the number of free ions available to bind yet also reduces the screening length. The latter has been shown to strongly affect the correlation exchange potential. To highlight the competition of these two effects, we plot the degree of ion association versus polymer concentration in the presence of varying amounts of added salt in Fig. 10. A rod-like structure is chosen to illustrate this trend. Three values of the segment length b were tested. Regardless of the addition of salt, the limiting plateau values of α decrease with increasing b, consistent with the results in Fig. 7.

In the concentrated regime, for all segment lengths examined, higher levels of salt addition yield greater degrees of ion association. In the dilute regime, the screening length is large, and the limiting association behavior depends on the value of b. When b is small, correlation effects dominate over entropic considerations. Accordingly, lesser amounts of added salt lead to larger limiting values of α due to less effective screening of electrostatic self-interactions. On the other hand, when b is large, the chains are significantly stiffened and the driving force for association due to correlation effects is weakened. Thus, decreasing salt concentration yields lower values of α across the entire concentration range. At intermediate values of b, a crossover behavior shown in Fig. 10(b) is observed.

As a final note, we emphasize that the results in the high polymer concentration regime are consistent with our treatment of counterion adsorption as a reversible binding process. Under these conditions, electrostatic interactions are significantly screened, as characterized by a weakened exchange chemical potential μ_α. The association equilibrium then results from the competition between the translational entropy of ions and local binding energy. The degree of association increases with salt addition because more species are available for binding.

We have examined how the electrostatic correlation free energy affects counterion binding equilibrium in polyelectrolyte solutions. The law of mass action method we applied to model this behavior is entirely generalizable. One potential extension is to consider both charge and mass density fluctuations in the correlation free energy, yielding an exchange potential for counterion binding dependent on both factors. Furthermore, ion-dipole interactions between bound charge-pairs and free ions, or solvation effects that are sensitive to monomeric polarizability, may play a significant role in ion binding. Within our model, studying these effects simply requires introducing additional contributions to the excess free energy which would enter the law of mass action as a secondary component of the exchange potential μ_α, through the application of Eq. (6).

When treating electrostatic correlations, we have found that the straightforward application of the Debye-Hückel free energy implies a reduction in the degree of ion binding at high concentrations, due to the improper removal of ion self-energy implicit within the model. For this reason, we stress that the application of the Debye-Hückel free energy in reversible binding models, or any model with a concentration-dependent interaction strength, is not appropriate. This applies to alternative approaches treating reversible ion binding, such as the recently developed transfer matrix formalism for polyelectrolyte complexation.^37,38

Incorporating the self-energy explicitly in the electrostatic free energy is necessary in such binding models, as is true for models with a spatially varying dielectric permittivity.⁴⁷ In both cases, the ion self-energy couples nonlinearly to the composition profile. One way to account for the self-energy of both ions⁴⁷ and polyions²² is through charge-smearing, which effectively assigns a finite size to each charge. Adopting this approach, we showed that the self-energy promotes ion binding in order to alleviate electrostatic self-interactions.

Ion binding is opposed by the translational entropy of small ions. In the dilute solutions, this is balanced by a drive for binding derived from the electrostatic correlation free energy. We find that in this regime, molecular compactness enhances the degree of binding and chains with more flexible architectures or greater degrees of polymerization contain greater fractions of bound ions which, we argue, results from the strong self-interaction among neighboring like-charges.

In concentrated solutions, long-range electrostatic interactions are screened and their effect on binding is negligible. The bound fraction increases with ion concentration and is solely governed by the competition between the ion entropy and the short-ranged binding free energy. This trend is apparent for all chain architectures studied, which includes fixed rod- and coil-like chains, molecules with arbitrary fractal dimension, and semiflexible chains with adaptive structures.

This work has focused entirely on single polyelectrolyte solutions in the presence of counterions and added salt. We have neglected the highly relevant coacervate forming systems, consisting of multiple polyelectrolyte species and salt. A subject of future work is to extend our analysis to solutions of oppositely charged polyelectrolytes, and examine how electrostatic correlations and chain structure effect both reversible charge association and phase equilibrium behavior.

S.F. and J.Q. are grateful for helpful communications with Kevin Shen and Professor Zhen-Gang Wang. S.F. is supported by the Terman Faculty Fund. J.Q. is supported by the 3M Non-Tenured Faculty Award and the Hellman Scholar Award. R.G.L. and A.S. acknowledge the support of the National Science Foundation under Grant No. DMR-1707640. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation (NSF).

Here we present a general formalism for deriving the electrostatic correlation free energy for polyelectrolyte solutions. The formalism applies to a solution of polyanions, counterions, and added salt in a homogeneous dielectric medium, i.e., water. We consider only monovalent charged species with the same molecular volume $v$⁠. The solution contains n_A polyanions of length N_A, n₊ cations, n₋ anions, and n_W solvent molecules. N_A+ monomer units on the chain harbor bound cations, resulting in $n+b=NA+$ bound and $n+f=n+−n+b$ free cations. The microscopic density of each of these species is given by

The term α ≡ N_A+/n_AN_A is the average fraction of monomer units with a bound ion, which effectively smears the charges along the polymer backbone. This is equivalent to assuming that bound ions are free to permute positions along the chain, which gives rise to the combinatorial contribution to the free energy in Eq. (4).

Following the previous field-theoretic work on simple electrolytes,⁴⁷ polyelectrolytes,^22,57 and neutral polymeric systems,^54,56 the density of each individual particle is spread across a normalized Gaussian, $Γ(r)=exp(−|r|2/2a2)(2πa2)3/2$⁠, where the distribution width a corresponds to the monomeric size. The smeared microscopic density for each species is defined as $ρ̃(r)=Γ(r)*ρ^(r)≡∫dr′Γ(r−r′)ρ^(r′)$⁠. We then define a smeared, microscopic charge density as

where the signed components originate from the signed valencies of each species. Introduction of the smeared charge density at this stage ultimately regularizes the UV-divergent correlation free energy obtained from this derivation. Finally, we impose an incompressibility constraint by inserting $δ1−∑iϕ^i$ into the partition function, where $ϕ^i=vρ^i$ is the microscopic volume fraction operator and the summation index i runs over all bound and unbound species. We transform this constraint into its exponential representation by writing

which introduces an auxiliary field variable η in the partition function.

The canonical partition function for our system can then be written as

$Z0$ contains the ideal gas partition function for unbound salt, polyelectrolyte, and solvent species, $ZC$ is the combinatorial partition function which leads to Eq. (4), and $Zad=e−n+bΔGα$ contains the energetic contribution of $n+b$ bound counterions along the polyelectrolyte chains resulting in the adsorption free energy, Eq. (5). The trace operation denotes integration over all particle coordinates, i.e., Tr ≡ ∏_γ ∏_js∫dr_γjs, where the index γ runs over all unbound species, j over each molecule of that species, and s over each monomeric unit of that particle.

The Hamiltonian appearing above is defined as

where βH₀ is an arbitrary bonded interaction between monomer units on the polyelectrolyte chain and

is the Coulomb energy with $G(r,r′)$ satisfying

The Coulomb interaction in the partition function can be decoupled from particle coordinates through a standard Hubbard-Stratonovich transformation^22,47,53 which introduces an effective electrostatic potential Ψ(r), yielding the following canonical partition function:

where $ZΨ=(detG)1/2$ is a normalization factor from the Gaussian functional integrals and the effective field-based Hamiltonian is

The single-particle/chain partition functions appearing above are defined as

where $Ψ̃(r)=Γ(r)*Ψ(r)$ represents a spatial convolution of the field and $∫Dr$ denotes integration over space for all monomers in the chain. The volume factors appearing in the above equations have been absorbed into the definition of $Z0$⁠. Note that no single-particle partition function appears for the bound ions, though a factor of α does appear in Q_A due to coupling of the coordinates between bound ions and the polyelectrolyte.

No approximations have been made to this point. To evaluate the partition function analytically, we adopt the standard Gaussian fluctuation approach and expand the partition function to the quadratic order in the field variables around their homogeneous saddle-point values. This results in expressions for the electrostatic fluctuation free energy that include both intra- and inter-molecular correlations between species with a given molecular structure, as seen in Eqs. (9) and (15).

As the focus of our work is on the effect of electrostatic correlations on charge binding, we only perform this expansion in Ψ and leave η at the mean-field level, which removes the functional integral over η and simply recovers the incompressibility constraint. Doing so ultimately yields a correlation free energy describing only the effects of charge density fluctuations, rather than the mass density fluctuations present even in neutral polymer systems.²⁰ Including fluctuations in both fields may quantitatively alter results from our model, though we note that treatment of fluctuations in a single field is common in studies of both neutral⁶³ and charged polymer systems^20,22 and avoids introduction of additional, arbitrary excluded volume parameters through the context of an implicit solvent model.^55,57,64,65

Performing these approximations, we obtain

where δΨ · [*] · δΨ ≡ ∫dr dr′ δΨ[*]δΨ, δΨ = Ψ − Ψ^* is the deviation from the mean-field electrostatic potential Ψ^*, and $Ω̃γ$ are the charge structure factors for salt and polyelectrolyte species to be defined shortly. The field integral in δΨ is Gaussian and can be completed to obtain a factor $det(G−1+∑γΩ̃γ)−1/2$ in the partition function, which when merged with $ZΨ$ yields the Gaussian fluctuation correction to the free energy,

where ∫_q denotes (2π)⁻³∫d³q. In a homogeneous dielectric medium, the Coulomb operator takes the familiar Fourier representation $G(q)=4πlB/q2$⁠. The charge structure factors are defined as

where g(q) is the intramolecular structure factor, which is unity for the salt species. In isotropic fluids, both $Ω̃$ and g depend only on the amplitude of the wavevector. We emphasize that through our formalism, the effects of charge binding enter solely as an effective charge density on the chain, 1 − α, and the free cation volume fraction $ϕ+f$⁠. Explicit correlations between counterions localized along the polymer chain⁴² and fluctuations in the bound fraction of ions are neglected.

The above expression for the correlation free energy, Eq. (A19), along with Eqs. (A20)–(A22), is amenable for numerical evaluation. With the structure factor of the polyelectrolyte specified, these expressions are applied directly in Sec. IV to model the effects of the chain structure. When discussing the correlation free energy in Sec. III through the smeared Debye-Hückel model [Eq. (11)], we make an additional approximation for the polyelectrolyte structure factor. In this, we treat the polyelectrolyte as a set of structureless point charges with charged fraction (1 − α)ϕ_A, resulting in a charge structure factor of $Ω̃A(q)=Γ^2(q)(1−α)ϕA/v$⁠. This term is identical to the cation (counterion) structure factor in the salt-free system and appears in our definition of Eq. (11). This adjustment is needed in order to map our model to the form of the classical DH free energy, which neglects the squared charge density dependence in the polyelectrolyte terms apparent in a more rigorous derivation of the electrostatic free energy.