Electrostatic correlation free energy for finite polymer chains

Abstract

Electrostatic correlation free energy (ECF) is the basis for modeling the thermodynamic behavior of polyelectrolyte solutions. In the past, it has mainly been estimated using the Edwards approximation, valid for infinite chains. Here, we show that the leading contribution due to finite molecular size is of order N⁻¹, regardless of the fractal dimension d, where N is proportional to molecular weight. This contribution is a local effect, originating from the missing correlations among connected charges near chain ends. In contrast, the contribution from the long-wavelength or infrared regime is weaker, of order N^−3/d ln N. Closed-form expressions for the free energy are provided for polyelectrolytes exhibiting either coil- or rod-like statistics, in the absence or presence of small ions. The consequence of the end effect is demonstrated by evaluating the phase diagram, surface tension, and molecular weight-driven partitioning.

---

1 Introduction

Solutions of oppositely-charged polyelectrolytes (PEs) exhibit a demixing transition known as polyelectrolyte complexation (PEC), owing to the electrostatic interactions among the charged groups.^1–3 PEC results in the coexistence between a supernatant that is nearly depleted of polymers, and a polymer-rich complex.^1–3 The complex may behave as a solid-like precipitate, or flow and coalesce like a liquid.^1–3 Such liquid-like complexes are often referred to as coacervates, and the demixing transition as coacervation.^1–3 PEC and the concentrated coacervate phases it produces have found various applications in, for example, drug delivery devices,^4,5 adhesives,^6,7 and food products.^8,9 Moreover, recent interest in PEC has been further spurred by its recognition as an important mechanism controlling the organization of biomolecules in intracellular environments.^10,11

The simplest form of solution free energy required for modeling PEC is obtained by combining the mixing entropy with contributions from electrostatic interactions. For typical PEC resulting from macrophase separation, the charge density vanishes everywhere, owing to global charge neutrality, so the electrostatic interactions cannot be handled at the mean-field level. A minimal treatment of charge density fluctuations is needed. The interaction of these charge density fluctuations gives rise to the electrostatic correlation free energy (ECF), which is intrinsically a non-mean-field effect. With few exceptions,¹² most literature focus on the charge-neutral system, for which ECF is the leading electrostatic contribution.

Different levels of treatment have been adopted to evaluate the ECF. In the Voorn–Overbeek (VO) model,¹³ the ECF is approximated by that of point-like ions, i.e., the Debye–Hückel correlation free energy for the fully ionized plasma.¹⁴ Despite its simplicity, the VO model reproduces several key features of PEC such as the effect of added salt and the degree of ionization.^1,2,15–17 The apparent success of the VO model is surprising because it essentially neglects charge connectivity by “chopping” polyelectrolytes into point-like species,^1–3,15 and because the Debye–Hückel correlation is expected to work only in the regime many orders of magnitude more dilute than typical PECs.^2,18

Following the landmark work of Borue–Erukhimovich,^19,20 many studies have attempted to account for the effect of polymer structure on ECF.^21–28 When applied to weak polyelectrolytes, i.e., those with low charge density, these works all rely on some version of the random phase approximation (RPA),^29–31 which treats the fluctuation of charge density at different wavelengths independently. These efforts parallel the earlier developments on composition fluctuation in neutral polymers,^32,33 and have been summarized in several recent reviews.^2,15,34,35 The effects of a wide range of molecular scale features have been explored, which include chain stiffness,^23,36 liquid-state packing,^37–40 counterion condensation,^{28,36,38,41–43} charge pattern,^44–51 and the formation of nematic phases and microstructures.^52–54

The RPA-based approaches assume that the density fluctuation is weak, which limits their application in dilute regime that is especially relevant for the description of the supernatant phase. Several developments have been made to handle the strongly fluctuating regimes, including variational approaches,^55,56 full field theoretical sampling,^46,57–62 scaling arguments,^44,63–66 and coarse-grained simulations on polyelectrolyte aggregation.^67–70

Except for recent work,⁷¹ the majority of these studies adopted a continuum description of the solvent. The static dielectric permittivity of the solvent is used to capture the effects of screening, and its dependence on temperature is employed to capture the effects of solvent reorganization entropy.⁶⁸ Once a picture of continuum dielectric media is adopted, the excess free energy of ionic solutions can be decomposed into two types of contributions: self-energy and ECF. The self-energy is closely related to the Born solvation free energy, and has been carefully examined for both small ions⁵⁵ and polyelectrolytes.⁵⁶

As far as phase separation in charge neutral solutions is concerned, ECF is the most important electrostatic contribution to free energy, except when the degree of ionization is variable.⁷² The VO model, despite its limitations, is still widely used in recent studies,^2,15,73–75 largely because it offers a compact description to the ECF based on the familiar picture of Debye–Hückel screening. In contrast, the more sophisticated approaches^{23,36,38,52,56,72,76} often require numerical calculations. With the simple intuition and molecular parameters buried in algebra, application in experimental studies is somewhat discouraged.

The Edwards approximation^32,33 simplifies the ECF significantly. When the structure factor of polyelectrolytes is replaced with that appropriate for infinite chains, the ECF in a charge-neutral solution of symmetric Gaussian chains is found to be given by^{19,20,23,28,37}

f^(coil)_Ed = c₀ϕ^3/4, (1)

where ϕ is the total volume fraction of the pair of polycations and polyanions. The expression for the prefactor c₀ will be given in Section 4.1, and depends on monomer size, charge density, and dielectric permittivity of solvents, but not on molecular weight (MW). Eqn (1) generalizes the Debye–Hückel correlation free energy in the VO model to flexible polyelectrolytes.⁷⁷ In case of the VO model, the ECF is proportional to k_BT per cube of the Debye length, where the Debye length is inversely proportional to the square root of ion concentration, k_B is the Boltzmann constant, and T the absolute temperature. In eqn (1), the ECF is k_BT per cube of the screening length λ, with λ ∝ ϕ^−1/4. In the presence of both polyions and small ions, a similarly compact expression has also been obtained,^19,20,23,63 which smoothly interpolates the behavior at low and high salt concentrations.

It should be stressed that the effect of charge connectivity is solely contained in the dependence on concentration ϕ. Different conformational statistics will lead to a different expression.²³ For instance, the polyions adopting rod-like conformations have the ECF²³

f^(rod)_Ed = −c_rϕ ln ϕ (2)

which underlies the plateau in osmotic coefficient observed for long polyeletrolytes.⁵⁶ Similar to c₀, the coefficient c_r is a combination of molecular parameters but is independent of MW.

The predictions of eqn (1) are consistent with the scaling arguments developed for PE solutions.^30,63,78 Accordingly, these approximate theories are most commonly used alongside scaling arguments to substantiate their conclusions and provide exact expressions for numerical prefactors. For example, previous works have used this approach to examine the dependence of the coacervate composition,^{23,44,63,79,80} interfacial tension,^63,80 and salt-resistance⁶³ on the charge density. In addition, a small number of works have employed eqn (1) and the generalized version containing salinity,^23,45,81 to obtain phase diagrams for salted polyelectrolytes that agree qualitatively with experimental results.^{1–3,15,23,45}

2 Why bother with finite MW

Eqn (1) is obtained by applying the Edwards approximation, which captures the correlations in the internal sections of polymer chains and thus is, strictly speaking, applicable only to “infinite” chains.^{19,20,23,37,52} For finite chains, we will show that the ECF has a correction of order N⁻¹, where N is proportional to MW and counts the number of repeating units. (In the remainder of this work, we use the terms MW and N interchangeably.) The correction applies to chains with different conformational statistics. For the most familiar example, Gaussian coils, the ECF has the form

Here c₁ is a combination of molecular parameters independent of molecular weight, similar to c₀. The dependence on N is comparable to the contribution of mixing entropy. In the regime of low concentration, the correction term actually has a stronger ϕ-dependence than the Edwards term, i.e., ϕ^1/4 ≫ ϕ^3/4, which competes with the weakening factor N⁻¹. Eqn (3) is the simplest form of the ECF for binary mixtures of coil-like polyions and solvent. More general expressions, valid for solutions with small molecule salt in addition to polyions with coil or rod-like conformations, are provided in Section 4.4.

The Edwards approximation is usually justified by noticing that, when the MW is sufficiently high, the ECF is dominated by the structural correlations in the middle portion of polymer chains. This argument holds for polymers in good solvent, where the composition fluctuation is an additional contribution to the mean-field level excluded-volume interaction. For symmetric polyelectrolyte solutions, due to charge neutrality, the leading-order contribution to ECF is itself a result of charge density fluctuation. In this case, the correction due to finite MW can be significant.

The accuracy of eqn (1) for the Edwards approximation and eqn (3) are illustrated in Fig. 1(a), which shows the pairs of binodal and spinodal curves for a solution of symmetric polycation and polyanion, calculated using the exact chain structure factor, the Edwards approximation eqn (1), and the finite-N corrected ECF eqn (3). The Edwards approximation clearly overestimates the magnitude of ECF, resulting in a wider coexistence window. Moreover, the critical point, i.e., the minimum value of charge density needed to induce the phase separation, is absent or vanishes. In contrast, the predictions from the finite-N corrected expression agree with the full theory nearly quantitatively.

Fig. 1 (a) Critical points, spinodal curves, and binodal curves for polyion solutions with N = 100 predicted by the different approximations for f^(coil)_corr: f_Ed and f_NEd. Binodal curves give the coexisting supernatant and coacervate phase polymer content ϕ for a given value of charge density σ. No critical points are shown for f_Ed, as no nonzero critical points are predicted. (b) Illustration of the missing correlations near the chain ends, responsible for the finite-N correction.

The Edwards approximation is valid when the screening length λ is smaller than the chain size . For systems undergoing PEC, the supernatant phase and the critical point fall within the regime of low concentration, where the screening length λ is large. When the value of λ exceeds chain size, the Edwards approximation breaks down. This can be resolved by retaining the full structure factor, which requires numerical evaluation of the ECF. In contrast, eqn (3) is accurate and is as compact as eqn (1), convenient to use in practice.

It should be noted that an N-dependent ECF has already appeared in the literature.^26,82 In these works, an approximate structure factor is designed to interpolate the behavior in both high-q (ultraviolet, UV), where the Edwards approximation applies, and the low-q regimes (infrared, IR), where the finite size of chains is important. The N-dependence predicted by these studies is a consequence of applying the specific interpolated structure factor, derived from the IR contribution. We show below that the IR contribution, when explicitly worked out, results in a contribution of order N^−3/2 ln N for Gaussian coils, which is weaker than the order N⁻¹ term contained in eqn (3). The N⁻¹ term is instead derived from the sub-leading behavior in UV regime. Therefore, the physics identified in the two sets of work should be differentiated. The order N⁻¹ term is stronger and is comparable to the mixing entropy, which is particularly relevant for discussing partitioning of a mixture of finite chains.

The derivation of eqn (3) and the analogous results for alternative chain conformational statistics, as well as those for systems with added salts, constitute the main results of our work. We shall show that the N-dependent correction in eqn (3) can be attributed to the missing correlations for charges near the chain ends. As sketched in Fig. 1(b), the charges near chain ends have fewer connected neighbors, which weakens the inter-charge correlation, an effect particularly relevant for long-ranged electrostatic interactions. Such correction is proportional to the concentration of chain ends, hence the factor N⁻¹ in eqn (3). The precise form of concentration dependence is derived from the conformational statistics, which counts the number of missing neighbors near chain ends.

The remaining sections are organized as follows. First, we introduce the model, define the relevant molecular parameters, and provide the general expressions for ECF in a multi-component mixture. Then, the detailed derivation for eqn (3) is presented and generalized to alternative chain conformational statistics, chain architecture such as ring polymers, and systems with added salts. The results section examines how the N-dependent correction affects salt resistance, interfacial tension, and MW-dependent chain partitioning between phases. Since the N-dependence is not specific to electrostatic interactions, we discuss additionally the relevance of the end effects when other forms of interactions are involved, which includes in particular the short-ranged excluded volume potential. The influence of the interaction range mainly shows up in the consideration of the missing correlations depicted in Fig. 1(b). The last section summarizes our main findings.

3 Model and free energy

We consider a liquid mixture of m species in a volume V at temperature T. Each species has n_i molecules, molecular volume V_i, and volume fraction ϕ_i = n_iV_i/V, where 1 ≤ i ≤ m. The mixture is treated as incompressible, and volume fraction is the same as composition. For convenience, a common reference volume v is used to calculate the number of monomers per molecule, N_i ≡ V_i/v. For polyions, N_i is also referred to as polymer length. It is important to note that the choice of v is arbitrary and does not affect the consideration of any physically measurable quantity.

The molecular architecture is described by the intramolecular correlation function, which can be factored into the form [capital Omega, Greek, tilde]_i(q) = ϕ_iN_ig_i(q)/v. The function g_i(q) is the single-molecule structure factor and q is the wavevector. In isotropic liquids, g_i and [capital Omega, Greek, tilde]_i only depend on the magnitude q = |q|. Small molecular species, such as solvents and ions, are modeled as point-like with constant structure factor g = 1. For polymers adopting Gaussian coil conformation, the structure factor is the Debye function g_D(x) = 2(x − 1 + e^−x)/x², where x = q²R_g² and R_g = N^1/2b is the radius of gyration with b being the statistical segment length.⁸³ For polymers with rod-like conformation, the structure factor is the Neugebauer function , where x = qL and L is the molecular length.⁸⁴ The Edwards approximation is recovered by setting either L or N in the structure factor to infinity, which will be discussed in Section 4.2.

The Hamiltonian for such mixtures contains contributions from intramolecular interactions H₀ and intermolecular interactions H_int, H = H₀ + H_int. The exact form of H₀ is left unspecified as its effects are captured by the appropriate choice of [capital Omega, Greek, tilde]_i for each species. The intermolecular interactions are assumed pairwise and are expressed as

Here, the notation is used to represent integration over Fourier space. The microscopic composition fields [small phi, Greek, circumflex]_i(q) are the Fourier modes of the instantaneous composition fields [small phi, Greek, circumflex]_i(r), . Here δ(r) is the Dirac delta function and the summation is performed over all monomers belonging to a given species i.

The interaction kernel Ũ_ij has two contributions, Ũ_ij = vB + U_ij. The first piece B is an effective bulk modulus which serves to suppress density fluctuation. The value of B is set to infinity after the expression for free energy is obtained, so that only composition fluctuations compatible with a constant density are included. The second piece U_ij encodes all other pairwise intermolecular interactions present in the system. In this work, we are primarily interested in the effect of Coulomb interactions. In Fourier representation, they are written as U_ij = 4πl_Bσ_iσ_j/q², where l_B = e²/(4πεk_BT) is the Bjerrum length.^23,29 The dielectric permittivity ε is assumed to be equal to that of the solvent species, typically water. Under these conditions, the Bjerrum length is l_B ≈ 7.6 Å at 300 K. The terms σ_i are charge densities, defined by σ_i = Z_i/N_i, where Z_i is the valency, i.e., the total charge carried by each molecule of species i. The charge density of neutral species is zero.

It has been previously shown^29,31,85 that the free energy for incompressible homogeneous liquid mixtures with pairwise interactions specified by eqn (4) can be expanded as

Here, the first term captures the translational entropy of mixing and is present in non-interacting systems. The second term includes contributions from intermolecular interactions at the mean-field level and is evaluated using the average values of the density fields 〈ϕ_i(r)〉 = ϕ_i which are independent of r for a homogenous system. The final term f_corr contains all corrections arising from composition fluctuations around the mean-field values ϕ_i and depends on the underlying molecular architecture and nature of interactions.

The first two terms in eqn (5) comprise the mean-field contributions to the free energy and are commonly applied to neutral polymer blends via the Flory–Huggins (FH) model, where U_ij is the FH parameter χ_ij.³¹ In this work, we focus on the Coulomb interaction. Upon substitution of the Coulomb kernel U_ij = 4πl_Bσ_iσ_j/q², the second term in eqn (5) vanishes, because is the average total charge density, which vanishes in homogeneous and charge-neutral systems. In such systems, the nontrivial contribution from electrostatics has to come from the microscopic fluctuation of charge density. The minimum model of free energy requires the consideration of fluctuation corrections contained in f_corr, which is the same as the electrostatic correlation free energy, i.e., ECF.

Several previous works^{15,19,20,23,24} have shown that the leading contribution to ECF is given by the random phase approximation (RPA) or Gaussian fluctuation theory. In terms of the notation introduced above, the RPA contribution to the free energy for arbitrary mixtures can be written^31,85

Here, the m × m matrices [capital Omega, Greek, tilde](q) and Ũ(q) encode the molecular structure and interactions. The entries of Ũ(q) are the interaction kernels Ũ_ij(q) given above. The diagonal entries of [capital Omega, Greek, tilde](q) are the intramolecular correlation functions [capital Omega, Greek, tilde]_i(q). The off-diagonal entries of [capital Omega, Greek, tilde] represent intermolecular cross correlations, which may be present, for instance, when block copolymers are considered. For incompressible systems dominated by electrostatic interactions, eqn (6) can be reduced to²³

where Ũ(q) has been replaced by s(q) and u. The vector s is a factored form of the Coulomb interaction potential with entries . The vector u = (1, 1, …, 1)^T contains m ones and arises from the incompressibility constraint.

The final term of eqn (7) that depends on u excludes the contribution from charge density fluctuations that are incompatible with the constant liquid density. Nearly all previous studies using the RPA to estimate ECF neglected the contributions from this coupling. This coupling can be dropped for systems with symmetric correlations among the positive and negative species, e.g. polycation and polyanion with the same structure factor and charge density or cation and anion with the same valency and radius, since u^T[capital Omega, Greek, tilde]s = 0. However, most practically relevant systems do not fall into this category, which includes for instance solutions of polyanions and their counterions. Since the focus of this work is the effect of finite MW, we shall restrict our analysis to symmetric systems, in which u^T[capital Omega, Greek, tilde]s = 0. The effects of such coupling in asymmetric systems will be addressed in an upcoming work, which clarifies the condition of charge symmetry and examines the contribution from excluded volume in the asymmetric systems.

Eqn (7) for ECF is valid for mixtures with arbitrary number of components and molecular architecture. To highlight the effects of finite MW, we first focus on symmetric solutions containing oppositely charged and otherwise identical polycations and polyanions, then study the effects of added salts in a later section. The salt-free system has three species (m = 3): polyanion, polycation, and solvent. For simplicity, polycations and polyanions are assumed to have the same degree of polymerization N, structure factor g(q), and composition ϕ_a = ϕ_c = ϕ/2. The charge densities are opposite, so that σ_c = −σ_a = σ. Due to the constraint of incompressibility, the composition of solvent is ϕ_w = 1 − ϕ. The structure factor and the charge density vector [capital Omega, Greek, tilde](q) and s(q) are given by

Substituting the definitions of eqn (8) into eqn (7) gives the familiar form of the RPA expression for ECF in symmetric solutions

The analogous expression with salt addition will be provided in Section 4.4.

4 Finite-MW correction to ECF

Eqn (9) serves as the basis for our discussion of finite MW. The Edwards approximation, which has been extensively used in the literature to evaluate ECF, amounts to replacing the structure factor g(q) with its limiting form in the high-q regime. In the following, we compare the results obtained from the Edwards approximation with eqn (9), and show that the difference revealed in Fig. 1 can be associated with the end effects for polyelectrolytes with finite MW. Such end effects are considered for coil-like and rod-like molecules, as well as for more general fractal conformational statistics. The effect of molecular topology is illustrated by considering ring polymers. A compact, generic form of the correction to the Edwards approximation is provided for these cases, and generalized to systems with added salt. The physical interpretation of the correction term is discussed and contrasted with the related results of past studies.^19,26,82

4.1 End effect on electrostatic correlation

We first consider polymers with coil-like conformations. The full ECF is obtained by replacing g(q) in eqn (9) with the Debye function mentioned above. The Edwards approximation amounts to approximating the Debye function with g_D(x) ≃ 2/x. Substituting it into eqn (9) leads to eqn (1),for the free energy density. Here is the screening length.^19,23 It measures the decay length of the screened potential around a test charge introduced to the solution, analogous to the Debye length for point-like plasma. The screening length can be written equivalently as . Since polymer volume V_p = Nv, radius of gyration squared R_g², valency Z = Nσ, and volume fraction ϕ are all measurable quantities, the value of λ is independent of the choice of v. Moreover, the dependence on molecular weight N cancels out, making λ an intrinsic length scale, which mainly depends on composition ϕ and Bjerrum length l_B. Eqn (10) is obtained under the assumption that the electrostatic correlation is dominated by that among the internal chain segments. It is valid for large values of N and has been employed in several recent analyses of experimental results.^80,86

The finite N correction is incorporated by including the sub-leading term in the Debye function: g_D(q) ≃ 2(x − 1)/x². Substituting this to eqn (9) gives

The integral is convergent in both low and high-q regimes. To seek the leading correction, we note that 1/(q²R_g²) scales as N⁻¹. Expanding the integrand perturbatively, keeping the terms up to order [scr O, script letter O](N⁻¹), and completing the integrals results in

f^(coil)_NEd = f^(coil)_Ed + f^(coil)_N (12)

which is eqn (3). The first part f_Ed is given by eqn (10). The second part f_N scales inversely with N, which can be attributed to the effects of chain ends.

The above expansion makes it clear that f_N is derived from the term −2/x² in the Debye function. Section S.1 shows, by decomposing the Debye function into monomer–monomer correlations, that this term originates from the missing correlation for monomers near the chain ends: compared to monomers in the middle of chain, they miss about half the connected neighboring beads, as highlighted in Fig. 1(b).

The expression for f_N can be understood as follows. In the limit of infinite N, the ECF is k_BT per correlation volume λ³. For finite N, this ECF is reduced by the concentration of “defects”, chain ends. Each polymer chain may spread itself into Nb²/λ² number of correlation domains. The inverse, λ²/(Nb²), gives the fraction of correlation domains containing chain ends. The product of this factor and λ⁻³ gives the reduction of ECF in eqn (10). The correction f_N is comparable to f_Ed in the dilute regime, when ϕ ≤ N⁻².

The N-dependent correction to f_Ed is derived from the correction to the Edwards approximation in the high-q or UV regime. In contrast, previous works^26,82 have considered the N-dependent structure factor by interpolating the Edwards approximation with the behavior of the Debye function in the low q or IR regime. In the IR regime, with q < R_g⁻¹, the Debye function has the form 1 − q²R_g²/3 or simply 1. The error introduced by replacing the full Debye function with the Edwards approximation is

The upper bound is set to R_g⁻¹ as the border of the IR regime. The dominating contribution to the integral scales as . Since R_g ∝ N^1/2, we find that

The dependence on N in the IR regime is indeed weaker than f_N ≃ N⁻¹, so we conclude that the leading N-correction to ECF sits in the UV regime, and is an end effect.

The effects of the finite-N correction are demonstrated in Fig. 1(a), which shows both the binodal and the spinodal curves for N = 100 calculated from eqn (9), (10), and (12). The procedure for generating the binodal curves is detailed in Section S.6, and is based on the balance of osmotic pressure and chemical potentials. It is clear that, with the end-correction added, eqn (12) agrees to the prediction of the full Debye function nearly perfectly. In particular, the lower critical point that cannot be captured by the Edwards approximation is restored, and the accuracy of the polymer concentration in the coacervate branch is remarkable.

The accuracy of the Edwards approximation, eqn (10), is expected to improve asymptotically as N increases. This trend is verified in Fig. 2. Although f_Ed predicts a vanishing supernatant composition regardless of the value of N, the agreement of the coacervate concentration with that of the full theory continues to improve as N increases. This figure and Fig. 1(a) show the simplest comparison between the Edwards approximation and the N-corrected expression. The more substantial comparisons, on the effects of salt, surface tension, and MW-dependent partitioning in multi-component mixtures will be presented in Section 5.

Fig. 2 Comparison of the effect of molecular weight on binodal points predicted by the Edwards approximation, f_Ed and the N-corrected Edwards approximation, f_NEd.

4.2 Dependence on conformational statistics

The end effect for Gaussian chains identified in the previous section is expected to work for other conformational statistics. We first consider the rodlike case, which has a unique logarithmic dependence on composition, then generalize it to statistics with fractal dimension 1 < d < 2. As noted in Section 3, the structure factor for rodlike objects is the Neugebauer function. It can be expanded in the high q regime as follows

Here L is the rod length. If the cross-section area is denoted A, the chain length is given by N = LA/v. For convenience, we introduce [small script l] ≡ v/A as the length of segment filling the volume v, whose value depends on the choice of reference volume and should not be confused with the Kuhn length. Then we also have N = L/[small script l]. Likewise, the charge density σ = Z/N is the number of elementary charges carried within one reference volume. Inserting the above high-q expansion into eqn (9) gives

Here the screening length is given by .

The first term in eqn (17) has been discussed in previous works.^23,52 The integral is divergent in the UV regime, which can be renormalized by subtracting a linear function of ϕ, i.e., q²ϕ ln[1 + (ϕq³λ³)⁻¹]. Such renormalization removes the self-energy of polyions and does not affect the thermodynamic properties. After regularization, the correlation free energy for rod-like chains is

Here, the leading N-independent term is the ECF for rodlike objects given by the Edwards approximation, identical to eqn (11) of ref. 23. The second term is the negative finite-MW correction, of order N⁻¹, analogous to eqn (13) for Gaussian chains. Following the treatment in the coil-like case, it can be shown that the contribution from the IR regime is of order N³ ln N, much weaker than the N⁻¹ dependence caused by end effects. So eqn (18) captures the leading correction due to finite MW.

The similarity between eqn (12) and (18) implies that the finite-N correction in the rod-like case can also be attributed to end effects. To show this cleanly, we consider a generic polymer with fractal dimension d (1 ≤ d ≤ 2), so that the rod-like and coil-like cases correspond to d = 1 and d = 2 respectively. Let the polymer size be R = N^1/da, where a is the microscopic length, e.g., b for coil-like or [small script l] for rod-like molecules. In the high-q regime, with q > R⁻¹, the structure factor can be written

The prefactor is synonymous to the definition of the fractal dimension d, which is the product between the probability of randomly selecting a monomer 1/N, and the number of connected monomers 1/(qa)^d. The correction term in parenthesis counts the missing correlation: 1/N is the fraction of ends and α/(qa)^d is the number of connected monomers that would be present on an infinite chain, as shown in Fig. 1(b). The parameter α is a numeral depending on the conformation statistics. Section S.3 compares the exact structure factors for rodlike chains, coils, and ring polymers with the corresponding predictions of eqn (19), showing excellent agreement for q ≫ R.

Substituting eqn (19) into eqn (9) and following the same line of analysis as for Gaussian coils or rod-like chains, outlined in Section S.4, we obtain the following ECF

This expression has the same form as eqn (12) and (18). is a generalized version of the screening lengths for Gaussian rodlike chains. The first term says that the ECF in a homogeneous polyion solution is k_BT per screening volume.²³ The second term accounts for the finite-N correction due to the missed correlations. The ratio (R/λ)^d gives the number of screening volumes that one chain spans, and its inverse is the faction of such volumes containing chain ends. The picture is general, and the conformational statistics enters mainly through counting the number of missing correlated segments near ends.

4.3 Finite but no ends: ring polymer

Not all polymers with finite MW have ends. In the following, we examine the effect of topology by evaluating the ECF for ring polymers adopting Gaussian statistics. The structure factor is given bywhere is the Dawson integral⁸⁷ and x ≡ q²Nb²/6. Following the treatment in the earlier sections, the last step keeps the leading terms in the high-q regime. Note that the radius of gyration for ring polymer with N monomers is Nb²/12. The definition for x is kept to facilitate the comparison with linear chains.

Substitution of the asymptotic form of the structure factor into eqn (9) gives the ECF including the order N⁻¹ correction, analogous to eqn (12),

The first term is the ECF for infinitely large N, which is identical that for linear chains as the chain topology is not discernible in this limit. The second term has the expected N⁻¹ dependence, but is positive, opposite to the case of linear chains.

The positive contribution to the ECF is due to the enhanced correlation of ring polymers. A ring polymer can visualized by forcing the ends of a reference linear chain with the same N to close, which results in a more compact molecular size. All monomers on a ring polymer are statistically identical. So we can choose a monomer labeled i = 1 as a reference. The correlation between this reference monomer and the other monomers, indexed i = 2, 3, …N, decays as i increases, reaches the minimum at i = N/2, then grows between i = N/2 and i = N. Such non-monotonic correlation is a long-range effect derived from ring topology, and is distinct from the behavior of linear chains, for which correlation decays away from the reference monomer. This difference is illustrated in Fig. 3. To demonstrate concretely that the positive correction is a direct result of non-monotonic correlations, Section S.2 provides a derivation of the approximate structure factor for ring polymers analogous to that provided for linear chains in Section S.1.

Fig. 3 Schematic comparison of the finite-N effects for ring and linear polymers: missing correlations near chain ends for linear polymers become enhanced correlations for ring polymers due to the constraint of end-to-end connectivity.

4.4 Addition of small ions

We extend the above treatment to solutions with added ions, and present the closed-form expressions for free energy including the contribution of finite MW. With the addition of salt, the correlation array and charge vector, [capital Omega, Greek, tilde](q) and s(q) introduced in Section 3, contain 5 components each. Denote the volume fractions of cations and anions by ψ₊ and ψ₋, and their valencies by z₊ and z₋. We need

Here, the structure factors of both cations and anions are set to the q-independent constant 1. The total polymer fraction is ϕ = ϕ_a + ϕ_c, and the total ion fraction is ψ = ψ₊ + ψ₋. The solvent volume fraction is ϕ_w = 1 − ϕ − ψ, constrained by incompressibility.

As noted earlier, we study symmetric mixtures. The polycation and polyanion have same charge density, structure factor, and equal composition, ϕ_a = ϕ_c = ϕ/2. By charge neutrality, we further require z₊ψ₊ = z₋ψ₋. For this set of symmetric parameters, the term coupling charge density and mass density in eqn (7) vanishes. The ECF for asymmetric mixtures will be studied in a future work. Under these conditions, by substituting eqn (23) into eqn (7), we obtain the following form of ECF

Due to charge neutrality, the contributions from cations and anions are combined to z₊z₋ψ = z₊²ψ₊ + z₋²ψ₋, which is proportional to the ionic strength. Eqn (24) can be evaluated by applying the N-corrected Edwards approximation, for both coil-like and rod-like cases.

For Gaussian coils, the ECF is written as the sum, f^(coil)_NEd = f^(coil)_Ed + f^(coil)_N, with

In the above, R_g² = Nb²/6 is the radius of gyration squared. The screening length from salt-free solutions is used as a reference which scales with polymer concentration as ϕ^−1/4. The inverse Debye length squared κ² ≡ 4πl_Bz₊z₋ψ/v is proportional to the salt composition ψ. The term f_Ed gives the ECF for infinitely long chains, and has been obtained in the past.^19,20,23 We stress that it is valid only when the coupling between charge density and mass density vanishes, therefore is limited to the symmetric mixtures. It is straightforward to verify that f^(coil)_Ed reduces to eqn (10) when ψ = 0 and to the Debye–Hückel correlation free energy −κ³/(12π) when ϕ = 0. The second term represents the combined contribution of finite-N end effects and screening of small ions. In the high-salt regime, this N-dependent correction decays as .

For rod-like molecules, the ECF can be decomposed similarly, f^(rod)_NEd = f^(rod)_Ed + f^(rod)_N, with

Here L is the molecular backbone length, and the reference length is the screening length of a salt-free solution of rod-like chains, introduced in Section 4.2, which scales with polymer concentration as ϕ^−1/3. The parameter ρ is the screening length with added salts, given by the unique positive solution to the algebraic equation ρ⁶/λ⁶ + κ²ρ² = 1, which implies that both κρ and ρ/λ are less than unity. An explicit expression for ρ is provided in Section S.5. The parameter θ is given by . It is easy to see that (ρ,θ) = (λ,2π/3) when ψ = 0, and that (ρ,θ) = (κ⁻¹,π/2) when ϕ = 0. In the limit ψ = 0, f^(rod)_Ed = −ln ϕ/(12π²λ³), which was obtained previously,²³ and the end-correction reduces to the form in eqn (18). In the limit ϕ = 0, f^(rod)_Ed reduces to −κ³/(12π) as expected, and the end effect scales as .

5 Results

The effects of the finite-N correction are demonstrated by applying the ECF to three problems investigated experimentally: the stability of the complex phase to addition of salt, the interfacial tension between the coacervate and supernatant phases, and the partitioning of finitie-MW polymers between coexisting phases.

5.1 Salt resistance

Many previous experimental^{43,73,79,86,88–92} and theoretical^{16,19,20,23,36,42,93} studies have estimated phase diagrams for solutions of oppositely charged polyelectrolytes in presence of monovalent salts. In Fig. 4(a), we compare the two phase coexistence diagrams for fixed MW N = 200, obtained using the full RPA calculation eqn (23) and (25) for the Edwards approximation, and the finite-N corrected expression eqn (26) respectively. The equilibrium compositions are obtained by equating chemical potentials and osmotic pressure, as detailed in Section S.6.

Fig. 4 (a) Critical, spinodal and binodal points predicted for a salted flexible polyelectrolyte system with N = 200 and σ = 0.12 using different approximations for f^(coil)_corr: the Edwards approximation, the N-corrected Edwards approximation, and the exact RPA expression. (b) Effect of molecular weight on predicted critical and binodal points using the Edwards approximation and the N-corrected Edwards approximation.

In all three cases, the biphasic window narrows as the salt concentration ψ is increased. In the limit ψ = 0, the coexisting compositions are the same as given by Fig. 1(a), which shows the nearly quantitative agreement between the finite-N corrected free energy and full calculation. Here we observe that the agreement between the two is similarly quantitative as a finite amount of salts are introduced. In particular, the locations of the critical points are nearly exactly captured by the free energy with finite-N correction. In contrast, the Edwards approximation substantially overestimates the salt resistance, resulting in a wider two-phase window and a maximum salt concentration nearly twice that of the exact value. The overestimation is consistent with our observation that the end-effects weaken the electrostatic correlations, thus lowering the driving force for phase separation.

The convergence of the prediction from the N-corrected ECF to that of the Edwards approximation is illustrated in Fig. 4(b). Increasing N widens the two phase window slightly and results in higher salt-resistance. This is partially related to the entropy cost of phase separation, which is captured by both free energy models, and is also related to the finite-N correction to the ECF, which is only captured by eqn (26). For the highest MW studied, N = 1000, the composition in the low salt regimes between the two models are nearly indistinguishable, except near the critical point. In the limit N → ∞, the two models should yield identical results, and the polymer composition at the critical point vanishes.

5.2 Surface tension

Several previous works^{3,16,19,20,80,94} have estimated the interfacial tension between the supernatant and coacervate phase. In particular, ref. 80 considered the variation of interfacial tension with N in salt-free solutions. A generic dimensionless free energy for polymers in implicit solvents is used as starting point,

Here, the first term is ideal gas entropy, and the second term f_∞(ϕ) is the excess free energy, assumed to be independent of N. By further assuming the solvent to be ideal and using the Edwards approximation to estimate the ECF, the authors set

f_∞(ϕ) = wϕ³ + c_∞^3/4ϕ^3/4 (30)

in which w is the third virial coefficient. The constant c_∞ contains the effect of l_B, σ, b, and all numerical prefactors as shown in eqn (10). Plugging eqn (29) into the Cahn–Hilliard formalism gives the surface tension of the formwhere γ_∞ corresponds to the interfacial tension in the limit N → ∞. The coefficient h depends on the form of free energy. Expressions for γ_∞ and h are provided in Section S.7.

The finite-N correction we identified in this work can be combined with the first term in eqn (29). Then, essentially the same steps can be taken to evaluate the surface tension. The details of this analysis are provided in Section S.7. The comparison of this result and the original one from the Edwards approximation is given in Fig. 5(a). With the end-correction, the surface tension is lower due to the weakened ECF. However, the difference between the two predictions is less than 5% for relevant values of N, and is likely negligible in most cases. Moreover, the scaling of both predictions with σ is identical, indicating that the relative contributions of entropy and electrostatics to the N-dependence of the interfacial tension are independent of charge density at this level.

Fig. 5 (a) Approximate surface tension between the supernatant and coacervate phase as a function of polyion size N predicted by the Edwards approximation and N-corrected Edwards approximation with b = 1, σ = 0.5, and w = 1/6. (b) Effect of molecular weight on the interfacial profile predicted by the Edwards approximation and the N-corrected Edwards approximation with b = 1 and w = 1/6.

The authors of ref. 80 also numerically estimated the interfacial profile between the supernatant and coacervate phase along with its leading-order dependence on N. We repeat this analysis and compare the resulting interfacial profiles in Fig. 5(b). The details of this analysis can be found in Section S.7. Overall, the interfacial profiles predicted by the Edwards approximation exhibit a weaker dependence on N than those that account for the finite-N correction to the ECF. The interface is asymmetric in both cases; however, accounting for the additional MW dependence increases the asymmetry slightly and broadens the interfacial density profile. Similar to the surface tension, the difference between the predicted profiles is within about 10% and rapidly shrinks as N is increased.

Overall, the N dependence of the ECF has only a modest effect on the interfacial tension and profile within this framework. It is worth noting, however, that several approximations have been made to arrive at the results in Fig. 5. Most relevant is the assumption of zero polymer content in the supernatant phase. This assumption is necessary in the original formulation, because the original ECF is too strong to allow for the finite concentration in the supernatant phase, as shown in Fig. 1(a).

5.3 Chain partition: electrostatic vs. entropic

Many synthetic polyelectrolytes have nonnegligible dispersity and, based on the binodal curves shown in Fig. 2, polymers with disparate MWs are expected to be partitioned differently. One would expect that a disperse mixture of polyelectrolytes undergoing complexation expels short polyions in favor of longer chains, resulting in a shift of the MW distribution. This possibility has been explored and studied in our recent work.¹⁷ Our study rationalized the experimental observations using the MW-dependence of translational entropy. The results obtained in the previous sections allow us to examine effects of the N correction on ECF, and compare it to that from the entropic arguments.

First, we note that within the Edwards approximation, such a scenario is not possible. The condition of zero supernatant composition requires that the bulk polymer content is entirely contained within the coacervate phase. Thus, although the absolute concentration of each component can vary from its bulk value due to solvent exchange, their distribution must remain fixed and equal to its bulk value. To demonstrate this, we use the results in Fig. 2 to approximate the expected distribution of polymer in each phase when solutions of two oppositely charged polyion pairs with different MWs are mixed. In principle, an accurate treatment of such a system requires a slight generalization of the approach from Section 4.2 to include additional polyion species with unique values of N_i in the expression for f_d, enabling generation of complete phase diagrams for these multicomponent systems. For simplicity, however, we rely on the simple mixing rule described below to illustrate the effects of N-correction to the ECF.

As a simple example, we consider a solution of two types of polyion pairs with N₁ = N and N₂ = rN (r ≥ 1), and assume that the supernatant and coacervate phase compositions can be estimated as

ϕ^(coac) ≃ xϕ^(coac)(rN) + (1 − x)ϕ^(coac)(N) (32)

ϕ^(sup) ≃ xϕ^(sup)(rN) + (1 − x)ϕ^(sup)(N) (33)

where ϕ^(coac)(N) and ϕ^(sup)(N) are the supernatant and coacervate binodal compositions for a single type of polyion with a given N. The quantity x ranges between 0 and 1, which controls the mixing ratio of the two polyion species. Each value of x produces a unique tie line from which the supernatant and coacervate compositions can be used to estimate the distribution of species 1 and 2 in each phase:

In this case, x^(coac) and x^(sup) correspond to the concentration of the higher-MW polymer, species 2 with N₂ = rN, normalized to the total polymer content.

To show concretely that nonzero supernatant concentrations are required to recover nontrivial partitioning behavior, we must also consider mass balance constraints. For each species, this takes the form of the lever rule given in eqn (S.17) where ν is the volume fraction of the coacervate phase and ϕ_i is the bulk composition of species i. Considering the Edwards approximation, for which ϕ^(sup)(N) = 0 regardless of N, combining the mass balance with eqn (34) shows that x^(coac) = ϕ_i/ϕ where ϕ is the total bulk polymer content. In other words, distribution of species in the coacervate is identical to the bulk solution if the supernatant concentration is zero.

For a given value of charge density σ, the binodal compositions predicted by the N-corrected Edwards approximation f_NEd can be substituted into eqn (34) to predict the expected distribution of polymer species in each phase. The results of this analysis are shown in Fig. 6 for σ = 0.06 and N = 100. In this case, we select a value of x = 0.5 to represent equal mixing of polyion species in the bulk. As shown, the predicted value of x^(coac) increases monotonically with r while x^(sup) decreases. At a ratio of r = 10, more than 90% of the coacervate phase polymer content is comprised of the longer polyion species; in contrast, the supernatant is nearly devoid of higher-MW polymer under the same conditions.

Fig. 6 Fraction of high-MW (N₂ = rN) polyion in the (a) supernatant (x^(sup)) and (b) coacervate (x^(coac)) phase predicted using f_NEd and f_Ed as a function of the MW ratio r = N₂/N₁ with N₁ = 100 and σ = 0.06. (c) Schematic illustration of the difference in predicted partitioning between the Edwards approximation with ϕ^(sup) = 0 and the N-corrected Edwards approximation with ϕ^(sup) > 0.

Although we emphasize that these results for f_NEd are not exact due to the use of a simple linear mixing rule in eqn (32) and (33), the overall trends of x^(coac) and x^(sup) with increasing MW mismatch agree qualitatively with expectations. In addition, the failure of the Edwards approximation to predict nontrivial partitioning is a direct consequence of mass balance constraints and thus it is not affected by this approximation. With this in mind, the stark disagreement between the predictions in Fig. 6 is a general feature of the Edwards approximation that would still be observed even in a more complete treatment of the multicomponent phase equilibrium.

6 Effects of interaction range

The end effect discussed in the previous sections is not limited to electrostatic interactions. Other interaction forms, such as excluded volume interactions, may be affected as well. It is thus necessary to study the circumstances under which the N-correction is relevant. To reveal the trend, we consider the molecules with the structure factor described by eqn (19) in the high-q regime and approaching the constant 1 in the low-q regime. Similarly, the interaction potential is assigned a generic form U_ij = u/(qa)^η, in which u and η control the strength and the range of the interaction, and a is the same microscopic scale in eqn (19). For electrostatic interaction, u is proportional to the Bjerrum length and η = 2; for excluded volume interaction, u is the excluded volume parameter and η = 0.

Focusing first on the UV regime, we substitute the structure factor and interaction potential into eqn (6), and obtain the correlation free energy

The coupling with the compression mode u is absent because the potential U_ij does not depend on species type. Following the previous analyses, the free energy is given by

Here the screening length is given by λ ≡ a(uϕ)^−1/(d+η). The interpretation of the two contributions are the same as the previous sections, with the end effect of order N⁻¹.

In the IR regime, extrapolating the Edwards approximation to q → 0 results in an error given by

The upper bound is the inverse of the molecular size R = N^1/da. The integrand is dominated by q² ln(N/q), so we have

Therefore, we find that the IR contribution is weaker than the end-effect captured by the contribution in the UV regime, irrespective of the fractal dimension and the range of interaction.

The N-effect has not been emphasized for the more important problem of excluded volume interaction because the mean-field contribution is dominant. Consider polymers in good solvent with solvent treated implicitly. The osmotic pressure is given by⁸³

The first two terms are the ideal gas contribution and the mean-field excluded volume interaction. The correlation contribution f_corr is given by eqn (36) with η = 0. For 1 ≤ d ≤ 2, the two correlation terms are of order ϕ^3/d and ϕ^3/d−1/N, respectively, both weaker than the mean-field level ϕ² term. In contrast, for charge neutral systems, the mean-field contribution of electrostatic interactions vanishes, and the ECF is the leading term. The result of the Edwards approximation necessarily has to compete with the N⁻¹ correction in the UV regime.

Conclusion

Our primary concern in this work is the finite-N correction to the ECF, electrostatic correlation free energy, which in the past has been mainly evaluated using the Edwards approximation. We showed that, for charge-neutral symmetric solutions, an additional correction of order N⁻¹ is needed for polymers of different conformation statistics and in presence or absence of small ions. This correction is attributed to the reduced correlation around chain ends. Closed-form expressions for the end-corrected free energy are provided for polymers adopting both Gaussian statistics, eqn (25) and (26), and rod-like statistics, eqn (27) and (28), which should facilitate the theoretical modeling of experimental results.

The finite-N correction we identified is derived from the behavior of the structure factor in the UV regime, which is distinct from that identified in the IR regime.^26,82 We showed that the IR regime contributes a free energy of order N^−3/2 ln N for Gaussian chains and of order N⁻³ ln N for rod-like polymers, which are both weaker than the N⁻¹ correction attributed to the missing correlation at chain ends. The analyses in these two regimes can be generalized to the other type of interaction potentials, such as excluded volume interactions. However, the end-correction is particularly important for charge-neutral system, because the mean-field contribution vanishes and the fluctuation effect becomes the leading contribution.

A few examples are discussed to examine the importance of the N correction. For binary systems containing no added salts, we showed that the Edwards approximation is too strong, so that the concentration of polymers in the supernatant is always zero. This is fixed by restoring the N correction to the ECF. Similarly, when salts are added, we showed that the N-corrected ECF agrees with the prediction of the full expression nearly exactly. The contributions of the finite-N correction to the interfacial tension and the partitioning of polyions with different MWs are also discussed. Given the simplicity of eqn (25)–(28), we recommend them to be used in place of the full RPA expression for the discussion of finite chains.

The scaling of the leading order correction to the ECF with N is the same as that of the translational entropy typically employed to study polymer solution behavior. As a consequence, it is difficult to validate the results of this work directly from the N-dependence of measurable quantities since the parametric dependence on N⁻¹ is identical to existing models.⁸⁰ Given the distinct concentration dependence of translational entropy and the ECF, however, combined measurements of both the concentration and N-dependence of quantities such as the osmotic pressure could provide a means to verify the key results of this work.

All the results in this work are obtained for symmetric systems, so that the term coupling the charge density and mass density in eqn (7) vanishes. The importance of this term, as applied to systems with conformational asymmetry, charge asymmetry, and excluded volumes will be presented in the future.

Conflicts of interest

There are no conflicts to declare.

Data availability

The algorithms for generating all numerical results from the provided equations are described in the main text and SI. Supplementary information: Decomposition of the structure factor calculation for linear coils (S.1) and ring polymers (S.2) to reveal the end effects, numerical comparison between the exact and the approximate structure factors (S.3), derivation of the ECF for molecules with general fractal dimension (S.4), salt-dependence of the screening length for rod-like chains (S.5), details of phase diagram calculations (S.6), and additional details on interfacial property calculations (S.7). See DOI: https://doi.org/10.1039/d5sm00633c

The code used to generate phase diagrams, interfacial properties, and partitioning behavior will be made available upon reasonable request.

Acknowledgements

J. Q. acknowledges support from Stanford University as the Keck Faculty Scholar. J. H. acknowledges support from the Stanford Bio-X Bowes Graduate Student Fellowship.

References

1. S. Srivastava and M. V. Tirrell, Advances in Chemical Physics, John Wiley & Sons, Ltd, 2016, pp. 499–544.
2. C. E. Sing and S. L. Perry, Soft Matter, 2020, 16, 2885–2914.
3. A. M. Rumyantsev, N. E. Jackson and J. J. de Pablo, Annu. Rev. Condens. Matter Phys., 2021, 12, 155–176.
4. M. Zhang, Q. Xiong, Y. Wang, Z. Zhang, W. Shen, L. Liu, Q. Wang and Q. Zhang, Polym. Chem., 2014, 5, 4670–4678.
5. K. A. Black, D. Priftis, S. L. Perry, J. Yip, W. Y. Byun and M. Tirrell, ACS Macro Lett., 2014, 3, 1088–1091.
6. R. J. Stewart, Appl. Microbiol. Biotechnol., 2011, 89, 27–33.
7. R. J. Stewart, C. S. Wang and H. Shao, Adv. Colloid Interface Sci., 2011, 167, 85–93.
8. S. L. Turgeon, C. Schmitt and C. Sanchez, Curr. Opin. Colloid Interface Sci., 2007, 12, 166–178.
9. C. Schmitt and S. L. Turgeon, Adv. Colloid Interface Sci., 2011, 167, 63–70.
10. S. F. Banani, H. O. Lee, A. A. Hyman and M. K. Rosen, Nat. Rev. Mol. Cell Biol., 2017, 18, 285–298.
11. S. Boeynaems, S. Alberti, N. L. Fawzi, T. Mittag, M. Polymenidou, F. Rousseau, J. Schymkowitz, J. Shorter, B. Wolozin, L. Van Den Bosch, P. Tompa and M. Fuxreiter, Trends Cell Biol., 2018, 28, 420–435.
12. I. I. Potemkin and V. V. Palyulin, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 041802.
13. J. T. G. Overbeek and M. J. Voorn, J. Cell. Comp. Physiol., 1957, 49, 7–26.
14. L. D. Landau and E. M. Lifshitz, Statistical Physics, Part 1, Butterworth-Heinemann, Amsterdam Heidelberg, 3rd edn, 1980, vol. 5.
15. C. E. Sing, Adv. Colloid Interface Sci., 2017, 239, 2–16.
16. J. Qin, D. Priftis, R. Farina, S. L. Perry, L. Leon, J. Whitmer, K. Hoffmann, M. Tirrell and J. J. de Pablo, ACS Macro Lett., 2014, 3, 565–568.
17. K. P. Barker, J. Horne, J. Lou, Y. Xia and J. Qin, Macromolecules, 2025, 58, 6249–6259.
18. P. Debye and E. Hückel, Phys. Z., 1923, 24, 185–206.
19. V. Y. Borue and I. Y. Erukhimovich, Macromolecules, 1988, 21, 3240–3249.
20. V. Y. Borue and I. Y. Erukhimovich, Macromolecules, 1990, 23, 3625–3632.
21. M. Castelnovo and J.-F. Joanny, Langmuir, 2000, 16, 7524–7532.
22. M. Castelnovo and J.-F. Joanny, Eur. Phys. J. E: Soft Matter Biol. Phys., 2001, 6, 377–386.
23. J. Qin and J. J. de Pablo, Macromolecules, 2016, 49, 8789–8800.
24. M. O. de la Cruz, L. Belloni, M. Delsanti, J. P. Dalbiez, O. Spalla and M. Drifford, J. Chem. Phys., 1995, 103, 5781–5791.
25. F. J. Solis and M. O. de la Cruz, EPJdirect, 2000, 2, 1–18.
26. K. A. Mahdi and M. Olvera de la Cruz, Macromolecules, 2000, 33, 7649–7654.
27. P. González-Mozuelos and M. Olvera De La Cruz, J. Chem. Phys., 2003, 118, 4684–4691.
28. A. V. Ermoshkin and M. Olvera de la Cruz, Macromolecules, 2003, 36, 7824–7832.
29. G. Fredrickson, The equilibrium theory of inhomogeneous polymers, Oxford University Press, Oxford, 2005.
30. P.-G. de Gennes, Scaling concepts in polymer physics, Cornell University Press, Ithaca (N.Y.) London, 1991.
31. D. C. Morse, Ann. Phys., 2006, 321, 2318–2389.
32. S. F. Edwards, Proc. Phys. Soc., London, 1965, 85, 613.
33. S. F. Edwards, Proc. Phys. Soc., London, 1966, 88, 265.
34. M. Muthukumar, Macromolecules, 2017, 50, 9528–9560.
35. M. Muthukumar, Physics of Charged Macromolecules: Synthetic and Biological Systems, Cambridge University Press, Cambridge, 2023.
36. S. Friedowitz and J. Qin, AIChE J., 2021, 67, e17426.
37. A. Kudlay and M. Olvera de la Cruz, J. Chem. Phys., 2004, 120, 404–412.
38. A. Kudlay, A. V. Ermoshkin and M. Olvera de la Cruz, Macromolecules, 2004, 37, 9231–9241.
39. S. L. Perry and C. E. Sing, Macromolecules, 2015, 48, 5040–5053.
40. C. E. Sing and J. Qin, Macromolecules, 2023, 56, 5941–5963.
41. A. Kundagrami and M. Muthukumar, J. Chem. Phys., 2008, 128, 244901.
42. A. Salehi and R. G. Larson, Macromolecules, 2016, 49, 9706–9719.
43. J. Lou, S. Friedowitz, J. Qin and Y. Xia, ACS Cent. Sci., 2019, 5, 549–557.
44. A. M. Rumyantsev, A. Johner, M. V. Tirrell and J. J. de Pablo, Macromolecules, 2022, 55, 6260–6274.
45. B. Yu, A. M. Rumyantsev, N. E. Jackson, H. Liang, J. M. Ting, S. Meng, M. V. Tirrell and J. J. D. Pablo, Mol. Syst. Des. Eng., 2021, 6, 790–804.
46. S. P. O. Danielsen, J. McCarty, J.-E. Shea, K. T. Delaney and G. H. Fredrickson, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 8224–8232.
47. L.-W. Chang, T. K. Lytle, M. Radhakrishna, J. J. Madinya, J. Vélez, C. E. Sing and S. L. Perry, Nat. Commun., 2017, 8, 1273.
48. T. K. Lytle and C. E. Sing, Mol. Syst. Des. Eng., 2018, 3, 183–196.
49. T. K. Lytle, L.-W. Chang, N. Markiewicz, S. L. Perry and C. E. Sing, ACS Cent. Sci., 2019, 5, 709–718.
50. C. E. Sing, J. Chem. Phys., 2020, 152, 024902.
51. J. J. Madinya, L.-W. Chang, S. L. Perry and C. E. Sing, Mol. Syst. Des. Eng., 2020, 5, 632–644.
52. R. Kumar, D. Audus and G. H. Fredrickson, J. Phys. Chem. B, 2010, 114, 9956–9976.
53. A. M. Rumyantsev and J. J. de Pablo, Macromolecules, 2019, 52, 5140–5156.
54. A. M. Rumyantsev and A. Johner, ACS Macro Lett., 2025, 14, 472–483.
55. Z.-G. Wang, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 021501.
56. K. Shen and Z.-G. Wang, J. Chem. Phys., 2017, 146, 084901.
57. Y. O. Popov, J. Lee and G. H. Fredrickson, J. Polym. Sci., Part B: Polym. Phys., 2007, 45, 3223–3230.
58. J. Lee, Y. O. Popov and G. H. Fredrickson, J. Chem. Phys., 2008, 128, 224908.
59. R. A. Riggleman, R. Kumar and G. H. Fredrickson, J. Chem. Phys., 2012, 136, 024903.
60. K. T. Delaney and G. H. Fredrickson, J. Chem. Phys., 2017, 146, 224902.
61. S. P. O. Danielsen, J. McCarty, J.-E. Shea, K. T. Delaney and G. H. Fredrickson, J. Chem. Phys., 2019, 151, 034904.
62. M. Nguyen, N. Sherck, K. Shen, C. E. R. Edwards, B. Yoo, S. Köhler, J. C. Speros, M. E. Helgeson, K. T. Delaney, M. S. Shell and G. H. Fredrickson, Macromolecules, 2022, 55, 9868–9879.
63. A. M. Rumyantsev, E. B. Zhulina and O. V. Borisov, Macromolecules, 2018, 51, 3788–3801.
64. A. V. Dobrynin, R. H. Colby and M. Rubinstein, Macromolecules, 1995, 28, 1859–1871.
65. A. V. Dobrynin and M. Rubinstein, Prog. Polym. Sci., 2005, 30, 1049–1118.
66. M. Rubinstein, Q. Liao and S. Panyukov, Macromolecules, 2018, 51, 9572–9588.
67. S. Chen, P. Zhang and Z.-G. Wang, Macromolecules, 2022, 55, 3898–3909.
68. S. Chen and Z.-G. Wang, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2209975119.
69. S. Chen and Z.-G. Wang, Phys. Rev. Lett., 2023, 131, 218201.
70. A. Glisman, S. Mantha, D. Yu, E. P. Wasserman, S. Backer and Z.-G. Wang, Macromolecules, 2024, 57, 1941–1949.
71. M. Beckinghausen and A. J. Spakowitz, Macromolecules, 2023, 56, 136–152.
72. S. Friedowitz, A. Salehi, R. G. Larson and J. Qin, J. Chem. Phys., 2018, 149, 163335.
73. E. Spruijt, A. H. Westphal, J. W. Borst, M. A. Cohen Stuart and J. van der Gucht, Macromolecules, 2010, 43, 6476–6484.
74. Y. Ma, R. J. S. Ivancic, J. Obrzut, D. J. Audus and V. M. Prabhu, Soft Matter, 2024, 20, 7512–7520.
75. W. Kang, X. Ma, C. Liu, S. Wang, Y. Zhou, C. Xue, Y. Xu and B. Li, Mater. Sci. Eng., R, 2024, 157, 100762.
76. K. Shen and Z.-G. Wang, Macromolecules, 2018, 51, 1706–1717.
77. P. B. Warren, J. Phys. II, 1997, 7, 343–361.
78. H. Schiessel and P. Pincus, Macromolecules, 1998, 31, 7953–7959.
79. A. E. Neitzel, Y. N. Fang, B. Yu, A. M. Rumyantsev, J. J. de Pablo and M. V. Tirrell, Macromolecules, 2021, 54, 6878–6890.
80. D. J. Audus, S. Ali, A. M. Rumyantsev, Y. Ma, J. J. de Pablo and V. M. Prabhu, Phys. Rev. Lett., 2021, 126, 237801.
81. L. Li, A. M. Rumyantsev, S. Srivastava, S. Meng, J. J. de Pablo and M. V. Tirrell, Macromolecules, 2021, 54, 105–114.
82. A. M. Rumyantsev and I. I. Potemkin, Phys. Chem. Chem. Phys., 2017, 19, 27580–27592.
83. M. Doi and S. F. Edwards, The Theory of Polymer Dynamics, Oxford University Press, Oxford, New York, 1988.
84. T. Neugebauer, Ann. Phys., 1943, 434, 509–533.
85. J. Qin and D. C. Morse, J. Chem. Phys., 2009, 130, 224902.
86. S. Friedowitz, J. Lou, K. P. Barker, K. Will, Y. Xia and J. Qin, Sci. Adv., 2021, 7, eabg8654.
87. W. Burchard and M. Schmidt, Polymer, 1980, 21, 745–749.
88. L. Li, S. Srivastava, M. Andreev, A. B. Marciel, J. J. de Pablo and M. V. Tirrell, Macromolecules, 2018, 51, 2988–2995.
89. H. Li, Y. Liu, F. Lan, M. Ghasemi and R. G. Larson, Macromolecules, 2023, 56, 7909–7920.
90. Q. Wang and J. B. Schlenoff, Macromolecules, 2014, 47, 3108–3116.
91. D. Iyer, V. M. S. Syed and S. Srivastava, J. Polym. Sci., 2021, 59, 2895–2904.
92. J. Horne, K. Barker, Y. Xia and J. Qin, Effects of divalent salts on polyelectrolyte coacervation, 2024, https://chemrxiv.org/engage/chemrxiv/article-details/6658d15691aefa6ce1667b1a.
93. M. Radhakrishna, K. Basu, Y. Liu, R. Shamsi, S. L. Perry and C. E. Sing, Macromolecules, 2017, 50, 3030–3037.
94. J. Wang, X. Chen, E.-Q. Chen and S. Yang, Macromolecules, 2024, 57, 9698–9710.

---