Criticality and Connectivity in Macromolecular Charge Complexation

Nov 4, 2016 - We examine the role of molecular connectivity and architecture on the complexation of ionic macromolecules (polyelectrolytes) of finite ...
0 downloads 8 Views 1MB Size
Article pubs.acs.org/Macromolecules

Criticality and Connectivity in Macromolecular Charge Complexation Jian Qin*,† and Juan J. de Pablo‡,§ †

Department of Chemical Engineering, Stanford University, Stanford, California 94305, United States Institute for Molecular Engineering, University of Chicago, Chicago, Illinois 60637, United States § Argonne National Laboratory, Argonne, Illinois 70439, United States ‡

ABSTRACT: We examine the role of molecular connectivity and architecture on the complexation of ionic macromolecules (polyelectrolytes) of finite size. A unified framework is developed and applied to evaluate the electrostatic correlation free energy for point-like, rod-like, and coil-like molecules. That framework is generalized to molecules of variable fractal dimensions, including dendrimers. Analytical expressions for the free energy, correlation length, and osmotic pressure are derived, thereby enabling consideration of the effects of charge connectivity, fractal dimension, and backbone stiffness on the complexation behavior of a wide range of polyelectrolytes. Results are presented for regions in the immediate vicinity of the critical region and far from it. A transparent and explicit expression for the coexistence curve is derived in order to facilitate analysis of experimentally observed phase diagrams.



interact through a Debye−Hückel energy term.5,9 Several efforts have been presented in attempts to incorporate the effects of connectivity.10−14 Most of them have focused on the limit of infinitely long chains, and the details have been summarized in a recent review.15 Such field-based approaches can be collectively grouped under the umbrella of the so-called “Gaussian theory”, which describes charge density fluctuations by relying on the random phase approximation (RPA).16 Higher order corrections to RPA-type theories have also been considered,17 but they have only been applied to investigate their effects on the spinodal curve. More elaborate developments have considered the possibility of microphase separation for relatively complex polymer architectures,18,19 some have coupled a field-base approach to liquid-state theory,8 and others have investigated the effects of charge discreteness.20 Elegant field theoretic sampling approaches based on complex Langevin simulations have also proven to be valuable21−23 but have the disadvantage of being highly computationally demanding. In this work, we focus on the effects of charge connectivity and molecular compactness as characterized by the fractal dimension. We first present a concise theory for thermodynamics of polyelectrolyte coacervation solutions without small ions and then examine how the original predictions of the VO model are altered by charge connectivity. We include in our treatment the full intramolecular structural correlation between charges along the polymers and present closed-form free energy expressions for chains of finite molecular weight. We then present a detailed analysis of critical behavior in polyelectrolyte

INTRODUCTION Polyelectrolytes of opposite charge in solution exhibit strong attractive forces that can lead to precipitation of solid aggregates or to the formation of a second liquid phase that is generally referred to as a “complex coacervate”.1 There is considerable interest in understanding polyelectrolyte complexation, which arises in a wide range of technologies in the food, pharmaceutical, medical, and coatings industries, to name a few.2,3 Polyelectrolyte complexation is driven by electrostatic interactions; under suitable conditions of salinity and pH, a coacervate phase typically contains 20−30% of polymer. The coexisting supernatant phase is nearly depleted of polymer.4 The formation of a coacervate has been generally understood to result from a balance between electrostatic, cohesive interactions, and the entropy of mixing. Coacervation has primarily been discussed within the context of the so-called Voorn−Overbeek (VO) model,5 which has been able to rationalize several features of the experimentally observed phase behavior, including the effects of macromolecular charge density and salinity. The model accounts for contributions to the free energy in the form of Debye−Hückel correlations6 and the entropy of mixing of all species present in solution. Its strength resides in its simplicity and the ease with which it can be adapted to describe relatively complex molecular systems. It does, however, neglect many potentially relevant physics that are thought to contribute to the thermodynamics of ionic polymers, such as charge regulation7 and the local ionic liquid structure.8 These deficiencies have been long recognized, and it is only recently that they have begun to be addressed. A particularly severe approximation within the VO theory is the absence of connectivity between charges along the polymer backbone. All charges are assumed to be isolated and to only © XXXX American Chemical Society

Received: September 27, 2016 Revised: October 20, 2016

A

DOI: 10.1021/acs.macromol.6b02113 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

ρ̂ i (r). The latter are uniquely specified by molecular configuration and are operationally defined by ϕ̂ i(r) = v∑α∈iδ(rα − r) and ρ̂i(r) = v∑α∈iσiδ(rα − r), where δ is the Dirac delta function and the summation runs over all monomers belonging to the ith molecule. To convert from mass density to charge density, we introduce the charge density σi on molecule i, which is defined by Zi/Ni, the ratio between molecular valence and size. The charge density for noncharged molecules, such as solvents, is set to 0. The first kernel Bi,j = vB suppresses the overall density fluctuations, with B being an effective bulk compressibility. The second kernel Ũ i,j(q) = 4πlB/ q2 represents the Coulomb interaction, where lB = e2/(4πϵkBT) is the Bjerrum length.27 In this work, we assume that the dielectric permittivity of the polyelectrolyte solution is the same as that of water. Thus, the Bjerrum length is close to 7 Å at room temperature. Finally, we note that by replacing the charge density with the mass density using ρi(r) = σϕi(r), B and Ũ can be combined, and only one interaction kernel U with Ui,j = vB + 4πlBσiσj/q2 need be evaluated. The free energy can be deduced from the Hamiltonian and partition function by perturbation theory. It has been shown26,28 that for homogeneous multicomponent liquids with negligible compressibility or large B values the free energy can be expanded as (eqs 47 and 59 in ref 26)

complexation and discuss how it is altered by the inclusion of charge connectivity and by molecular compactness. Finally, we generalize the electrostatic correlation free energy results to solutions with an arbitrary number of ionic species. We use these to derive a closed-form expression for the free energy of systems with added ions and examine the effects of added salt when connectivity is included. Particular emphasis is placed on the subtleties of the renormalization associated with the high-q behavior of molecular conformations, which has not gained sufficient attention in past work.



COMPLEXATION FREE ENERGY Our analysis is focused on ionic macromolecules dissolved in aqueous solution. Macromolecular architecture is described in terms of the intramolecular correlation function or the structure factor. The electrostatic cohesive energy of these solutions is of order kBT/λ3, where λ is the correlation length for charge density fluctuations. We present explicit results for the solution free energy and correlation length in what follows and show that the VO model is just a special case of a larger family of models. Including charge connectivity affects the nature of complexation in a significant, qualitative manner, and the result correlates strongly with the fractal dimension24 of the underlying macromolecules. We consider a liquid mixture composed of m components. Each component i (1 ≤ i ≤ m) occupies a volume fraction ϕi. Its internal structure is described by the intramolecular ϕN correlation function Ω̃i(q) ≡ i i g (q). Here v is a reference v

f≡ =

i

2 x

x

(∫

0

dt

sin t t



1 − cos x x

),

∑∫

q

i,j

1 + 2

i,j

q

ln(ϕi) +

v 2

∑∫ i,j

q

⟨ϕi(q)⟩Ui , j(q)⟨ϕj(− q)⟩ + fcorr

Here the first term represents the mixing entropy. The second term is the mean field interaction term, which is evaluated from average density fields ⟨ϕi(q)⟩ = ϕi and is independent of q value for homogeneous mixtures. The third term, fcorr, represents all fluctuation corrections and depends on molecular architecture (or structure factor). Explicit expressions can be found in the literature.26 Keeping the first two terms in eq 2 yields the mean field free energy and generalizes the Flory−Huggins lattice theory to a fluid in a continuum. To reveal the latter fact, we replace Ui,j by the Flory−Huggins χ-parameter χi,j. However, this level of description cannot be applied to ionic systems that are dominated by Coulomb interactions. Since the total charge density vanishes for charge neutral homogeneous systems, the interaction term that depends on lB vanishes, and the free energy only contains the mixing entropy term. This leads to one point we would like to stress: the theory for ionic systems cannot be of a mean field type, and fluctuations fcorr must be included. The lowest order fluctuation contribution to fcorr is the Gaussian theory. It has been worked out previously (see e.g., ref 28), and in our notation it is given by

in which x ≡ qL

fcorr =

ϕi(̂ q)Bi , j (q)ϕj(̂ −q)

∑∫

Ni

(2)

and L is the length of the molecular rod. The Hamiltonian for such a mixture can be symbolically written as H = H0 + Hint. Here H0 refers to contributions from intramolecular correlations, which may be taken as the Gaussian energy for coil-like molecules, and is left unspecified since its effects can be represented by the intramolecular correlation function Ω̃i.26 For liquids with pairwise interactions, including dispersion and Coulomb interactions, the interaction energy within a reference volume v may be written in terms of the microscopic mass density fields ϕ̂ i and charge density fields ρ̂i as Hint 1 = kBT 2

∑ i

volume, and Ni is the number of repeat units of the ith molecule, defined by the ratio between the molecular volume Vi and v. The function g is the structure factor, and q is the wave vector. In isotropic liquids, both Ω̃ and g depend only on the amplitude of the wave vector. For point-like solvents, the structure factor g(q) is unity. For ideal Gaussian chains, the structure function is the Debye 2 function g (q) = gD(x) ≡ 2 (x − 1 + e−x), in which x ≡ q2Rg2 x and Rg is the radius of gyration of the Gaussian coil.16 For stiff chains, which adopt a rod-like configuration, the structure function is the Neugebauer function, 2 5 g (q) = gN(x) ≡

Fv kBTV ϕi

v 2

∫q ln det(I + Ω̃(q)U(q))

(3)

Here the intramolecular correlation matrix Ω̃ has dimension m × m. Its diagonal entries are Ω̃i(q), and the off-diagonal terms represent cross correlations. The interaction matrix U has been defined previously, and I is the m × m identity matrix. Although eq 3 can be applied to a general multicomponent fluid, in this section we focus on ternary mixtures containing a solvent, a polyanion, and a polycation. The solvent is typically

ρî (q)Uĩ , j(q)ρĵ ( −q) (1)

For Fourier integrals the convention ∫ q = (2π)−3∫ d3q is adopted. The density fields ϕ̂ i(q) and ρ̂i(q) are Fourier modes of the instantaneous mass and charge density fields ϕ̂ i(r) and B

DOI: 10.1021/acs.macromol.6b02113 Macromolecules XXXX, XXX, XXX−XXX

Macromolecules



EXAMPLES 1. Point-like Objects. These correspond to the usual ionic solutions, for which we have Ω̃2 = ϕ/(2v) and σ = 1. The Coulomb correlation free energy for such systems reads

water and is treated as a structureless medium (g = 1). The polyanion and polycation are symmetric, in the sense that they have identical molecular architectures and carry equal amounts of charge, albeit of opposite sign. To ensure charge neutrality, the volume fractions of polycation and polyanion are identical and are denoted by ϕ/2. The solvent volume fraction is therefore given by 1 − ϕ. The degree of polymerization of the polyions is the same and denoted by N = Vpolyion/v. The choice of reference volume does not affect the solution thermodynamics. For convenience, we use the volume of a water molecule 30 Å3 as the reference. The polyionic valences are Z = Z+ = −Z−, and the charge densities are σ = σ+ = −σ−, respectively. Applying these identities simplifies U to ⎛B B ⎞ B ⎜ ⎟ U(q) = v⎜ B B + ϵ B − ϵ ⎟ ⎜ ⎟ ⎝ B B − ϵ B + ϵ⎠

fcorr =

∫0



⎛ κ2 ⎞ dq q2 ln⎜1 + 2 ⎟ q ⎠ ⎝

(4)

where the dimensionless quantity ϵ ≡ 4πlB/(vq ) measures the strength of Coulomb interactions. Similarly, the diagonal entries of the intramolecular correlation matrix can be denoted by Ω̃11 = Ω̃1 and Ω̃22 = Ω̃33 = Ω̃2 = ϕNg(q)/(2v), respectively. Substituting the above simplification into eq 3 yields the desired result. In this work, we are primarily interested in incompressible fluids. To reach this limit, we keep all B-dependent terms in the intermediate steps and then set them to infinity in the final expression. After a few algebraic manipulations, we find

(point) f corr =

v 4π 2

∫0



⎡ ⎛ κ 2σ ⎞ κ 2σ ⎤ dq q2⎢ln⎜1 + 2 ⎟ − 2 ⎥ ⎢⎣ ⎝ q ⎠ q ⎥⎦

(4π )1/2 ⎛ lB3 ⎞ ⎜ ⎟ 3 ⎝ v ⎠

ϕ3/2 (8)

This is precisely the Debye−Hückel free energy that was originally built into the Voorn−Overbeek model:5,9 the correlation free energy for a plasma of disconnected ions is −vκ3/(12π) or −kBT per Debye volume. To recover the VO expression, one replaces the plasma fraction ϕ by the value converted from the polyion, i.e., σϕpolyion. Equations 2 and 8 show explicitly that the Debye−Hückel correlation free energy does not come from a mean field treatment. It is in fact the lowest order (Gaussian) term for fluctuation correlation contributions to the mixture’s free energy. A minimal treatment of the free energy inevitably involves composition fluctuations, at least at the order of the Debye length κ−1. 2. Rod-like Objects. Such objects can be treated ϕN analogously. We replace Ω̃2 by 2v gN(qL) and find

(5)

In eq 5, we have only kept terms up to the linear order in the modulus B. The other terms are of lower order in powers of B. The modulus B that approaches infinity drops out from the other terms after taking the logarithm in eq 3. It does not depend on the composition ϕ; it is therefore thermodynamically inconsequential and can be neglected. The correlation free energy then becomes

∫q [ln(vΩ̃1 + 2vΩ̃2) + ln(1 + 2vΩ̃2ϵσ 2)]

(7)

1/2

=−

det(I + Ω̃(q)U(q)) ≃ (1 + 2ϵσ 2v Ω̃ 2)(Ω̃1 + 2Ω̃ 2)vB

v 2

v 4π 2

Here the term κ2 ≡ 4πϕlB/v is proportional to the ionic strength and κ−1 is the Debye length. The integral is divergent in the high-q regime, which is associated with the short distance divergence of the Coulomb 1/r potential for r → 0. The divergence can be regulated by subtracting from the integrand the self-energy contribution,20,30,31 which physically amounts to choosing the electrostatic energy for a plasma as the reference. We therefore write

2

fcorr =

Article

(6) (rod) = f corr

The first term does not depend on Coulombic interactions and arises from osmotic compositional fluctuations, which Edwards has investigated in the study of neutral polymer solutions28,29 and which is neglected in several related fluctuation field theories.13,17 The original VO model neglected this term, presumably because it was constructed by coupling the mixing entropy with the Debye−Hückel correlation free energy in an ad hoc fashion. It therefore lacked several details that are revealed by our more systematic treatment. In this work, we focus on how connectivity changes the predictions of the original Voorn−Overbeek formulation, and we therefore also neglect the contributions arising from osmotic compositional fluctuations. Its effects and a generalization to more complicated systems will be considered in a future publication. The second term in eq 6 corresponds to the electrostatic correlation free energy. Replacing the term Ω̃2 by the proper structure function generates the lowest-order complexation free energy for the corresponding system. In what follows, we consider three examples that are frequently studied: point-like, rod-like, and coil-like molecules. The first example is closely related to the original VO theory.

v 4π 2

∫0



⎛ ⎞ κ 2σ 2 L dq q2 ln⎜1 + 2 gN(qL)⎟ q p ⎝ ⎠

(9)

Here gN is the structure factor for rod-like molecule first calculated by Neugebauer.25 The effective packing length p is defined by p ≡ v/A, where A is the intersection area of the rod. So the effective degree of polymerization N can be calculated by N = LA/v = L/p, L being the rod length. The low-q behavior of gN is 1 − q2L2/36 = 1 − q2Rg2/3 with Rg2 = L2/12. The high-q behavior of gN is π/(qL). Similar to the case of point-like objects, the free energy expression eq 9 exhibits an ultraviolet divergence. It can be regulated by subtracting from eq 9 a linear function of composition ϕ (rod) f corr =

v 4π 2

∫0



⎡ ⎛ ⎞ κ 2σ 2 L dq q2⎢ln⎜1 + 2 gN(qL)⎟ ⎢⎣ ⎝ q p ⎠

⎛ 4π 2lBσ 2 ⎞⎤ ⎟⎥ − ϕ ln⎜1 + vpq3 ⎠⎥⎦ ⎝ C

(10) DOI: 10.1021/acs.macromol.6b02113 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules This regulation or renormalization procedure is validated by noting that the linear terms of ϕ are inconsequential for discussion of thermodynamic properties. Using the above expression to analyze the complexation behavior of rod-like molecules is somewhat tedious. Since the results are dominated by contributions from high-q modes, we replace the structure function by its high-q asymptotic form, gN ≃ π/(qL), and find (rod) f corr ≃

v 4π 2

∫0



(coil) f corr ≃

=

(11)



Equation 11 can also be written as −σ2Aκ2 ln(ϕ)/(12π) in terms of the Debye length. We further note that this form can be assimilated into the mixing entropy term but has a negative prefactor that favors complexation. The molecular weight dependence drops completely since we kept only contributions from high-q modes. 3. Coil-like Objects. This model is more appropriate for polyelectrolyte solutions. We substitute the intramolecular ϕN correlation function Ω̃ 2 = g (q2R g 2) and find v 4π 2

∫0



⎛ ⎞ κ 2Nσ 2 2 2 ⎟ dq q2 ln⎜1 + g ( q R ) g q2 D ⎝ ⎠



⎛ ⎞ 2κ 2Nσ 2 dq q2 ln⎜⎜1 + 4 2 ⎟⎟ q Rg ⎠ ⎝

3 1/4 23/2 ⎛ lB v ⎞ ⎟ σ 3/2ϕ3/4 ⎜ 6 1/4 b (3π ) ⎝ ⎠

(13)

EXTENSION TO FRACTAL MOLECULES

The examples discussed in the previous section can be generalized to molecules of arbitrary fractal dimension d. The point-like, rod-like, and coil-like molecules can be regarded as molecules with fractal dimension d = 0, 1, and 2, respectively. In this section, we show how the generic free energy expression, eq 3, can be applied to arrive at explicit results for fractal molecules of variable dimensions. This practically may be applied to dendrimers with varying generation indices and fractal dimensions.32 We begin by noting that the main input to eq 3 is the intramolecular correlation function, which is related to the pair distribution function g(r). At the scaling level, the pair distribution function for a molecule with fractal dimension d behaves as g(r) ≃ rd−3. This translates into an intramolecular ϕ 1 correlation function in q-space of the form Ω̃ 2(q) = 2v d .

2v D

(coil) = f corr

∫0

which agrees with the expression obtained previously for infinite chains.11 When expressed in terms of the Debye length, the above expression may also be written as (bp2κ3)1/2/(31/4π), where p ≡ v/b2 is the packing length. As expected, the dependence on molecular weight drops out. The concentration dependence becomes ϕ3/4 instead of the form −ϕ3/2 that appears in the original VO model.10 Note that both forms serve to stabilize polyelectrolyte solutions. A quantitative analysis is presented later in this work.

⎡ ⎛ 4π 2lBσ 2 ⎞ dq q2⎢ln⎜1 + ϕ⎟ ⎢⎣ ⎝ vpq3 ⎠

⎛ lBσ 2 4π 2lBσ 2 ⎞⎤ ⎥ ⎟ − ϕ ln⎜1 + = − ϕ ln(ϕ) 3p vpq3 ⎠⎥⎦ ⎝

v 4π 2

(12)

In the low-q regime, the integral behaves as −1/(2π2)∫ 0dq q2 ln(q) < 0. In the high-q regime, the integral behaves as 3κ2Nσ2/ (π2Rg2)∫ ∞dq q−2 and is convergent. This result should be contrasted with past statements in the literature,19,23 which proposed that the electrostatic correlation free energy is divergent in the high-q regime. In past work, a Gaussian smearing function was introduced for the charges, and the width of the smearing function became an adjustable parameter in the model. We show here that this smearing procedure is not unavoidable in order to regularize the theory. Physically, this is because that spreading the charges along the Gaussian chains continuously reduces the strength of self-interactions. In the case of point-like objects, only one length scale arises, namely the Debye screening length κ−1. For rod-like and coillike objects, an additional length scale, the radius of gyration of the polymer Rg, is also needed. The Debye length κ−1 characterizes the screening length. The radius of gyration Rg characterizes the range of charge connectivity. The two terms collectively dictate the strength of correlation contributions to the free energy. A quantitative assessment of such contributions is made in the next section in order to illustrate how charge connectivity affects phase behavior. In doing so, we also discuss the consequence of chain stiffness, as reflected in the value of Rg for fixed N, on complexation. In practice, similar to the case of rod-like molecules, eq 12 is not particularly convenient to use because of the integral over the Debye function. However, since the lower-q behavior of the integrand is suppressed greatly by the factor “dq q2”, we may focus on the high-q regime. In this regime, the Debye function for a Gaussian chain is well approximated by 2/(q2Rg2), the Edwards function.16 Substituting it into eq 12 yields

(qa)

Here a is a microscopic length scale corresponding to the reference volume v. For flexible polymers with statistical segmental length b, the Edwards expression gives d = 2 and a = b/ 12 . For rod-like molecules, it was shown in the previous section that d = 1 and a = p/π = v/(πA). For pointlike objects, d = 0 and no microscopic length a is needed. Then the correlation free energy, eq 3, can be written in a succinct manner as (d) f corr =

v 4π 2

∫0



⎡ ⎛ 4πlBσ 2 ϕ ⎞ ⎟⎟ dq q2⎢ln⎜⎜1 + ⎢⎣ ⎝ q2v (qa)d ⎠

⎛ 4πlBσ 2 1 ⎞⎤ ⎟⎟⎥ − ϕ ln⎜⎜1 + q2v (qa)d ⎠⎥⎦ ⎝

(14)

Here, as was done previously for the case of point-like and rodlike objects, we have subtracted a linear term in ϕ to regulate the free energy. Equation 14 can in fact be solved to yield the following analytical expression (d) f corr

csc =

( d3+π 2 )

12π

v (ϕ3/(d + 2) − ϕ) λ03

(15)

where we introduced an auxiliary length λ0 ≡ (vad/ (4πlBσ2))1/(d+2). Neglecting inconsequential linear term, we can further write D

DOI: 10.1021/acs.macromol.6b02113 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 1. Left: spinodal (dashed) and binodal (solid) coexistence curves for symmetric polyelectrolyte solutions with N = 100. Computed from the Voorn−Overbeek model and from eq 12 for coil-like polymers using a packing length p = lB. Right: comparison between the full theory, eq 12, and the approximation based on the Edwards kernel, eq 13.

Figure 2. Molecular weight dependence of the critical composition and charge density for coil-like molecules. Calculated for b2/(6lB2) = 1 using the water molecular volume as the reference v. The Bjerrum length lB is calculated by using the dielectric permittivity of water at room temperature.

(d) f corr

csc =

ϕ* ≃ (v/a3)(a/Rg)d−3, which depends on molecular weight. Above this concentration, the solution may be microscopically regarded as a fishnet with characteristic mesh size decaying according to ξ = ϕ1/(d−3). In that limit the free energy no longer depends on the size of the molecules. A qualitative conclusion emerging from eq 16 is that complexation is more likely to occur for molecules with higher fractal dimension d. Therefore, we expect that in going from point-like to rod-like to coil-like macromolecules, the chances for complexation increase.

( d3+π 2 ) v

12π

λ3

(16)

( d3+π 2 )/(12π )

Equation 16 may be conceived as csc

per

−1/(d+2)

. correlation volume, with a correlation length λ ≡ λ0 ϕ The correlation length is the same as the Debye length for point-like molecules but may depend on other monomer-level length scales for molecules having a more complicated architecture. This free energy density is negative for d ≤ 1 (point-like and rod-like) and is positive for d > 1. For the special case of d = 1, we can replace the diverging term csc



CONNECTIVITY EFFECTS AND CRITICAL BEHAVIOR Charge connectivity adds entropic contributions to the free energy. We examine its effects by comparing coexistence curves in the charge density−composition domain calculated using the VO model and using eq 12. The results are shown in the left panel of Figure 1. The spinodal curves are obtained by setting to zero the second derivatives of the free energy with respect to composition. The binodal curves are calculated by equating the chemical potentials of all components. In both cases, a minimum charge density σc is required to induce phase separation. Above σc, the two-phase coexistence window opens up; the left side and the right side correspond to the supernatant and complex coacervate phase, respectively. The main difference between the two cases is that the minimum charge density is much lower for the case with charge connectivity, nearly one-third of the value predicted by the VO model for this particular molecular weight. Furthermore, the two-phase window is considerably wider. Both observations suggest that charge connectivity enhances the ability to form complex coacervates. To show that the difference revealed by Figure 1 does not depend on the particular choice of molecular weight, we analyzed the scaling of the critical (minimum) charge density σc and the critical composition ϕc with molecular weight, for fractal molecules of all fractal dimensions in the Appendix. We

( d3+π 2 ) by −ln(ϕ)/π, which is obtained by taking the limit of

eq 15 as d → 1. For coil-like object, the scaling of correlation length with composition noted previously11,17 agrees with our general results. The concentration dependence of the correlation length λ can be understood as follows. Because of connectivity, within the range λ surrounding any test charge, one can always find σλd other charges that belong to the same molecule. Each of these charges interacts with λ3ϕσ other charges inside the correlation volume λ3. The interaction strength is governed by the Bjerrum length and scales as lB/λ. Setting the total Coulomb energythe multiplication of these three termsto be 1 (or kBT), we find λ dσ(λ 3ϕσ )

lB =1 λ

(17) −1/(d+2)

which leads to the scaling λ ∝ λ0ϕ that was derived in the preceding paragraphs. This correlation length represents the length over which the electrostatic potential is screened. Equation 16 represents a key result of this work, which reduces to eqs 8, 11, and 13 when proper substitutions are made for a. We expect the above approximation to be valid in the semidilute regime, i.e., above the overlapping concentration E

DOI: 10.1021/acs.macromol.6b02113 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 3. Terminal exponent for molecular weight dependence of critical volume fraction ϕc and charge density σc for different fractal dimensions. The apparent jump in exponent of ϕc at d = 1 is compensated by the asymptotically large prefactor J3(sc) in eq 39.



found that the critical composition scales with N−1 for both the VO model and for coil-like molecules, whereas the critical charge density σc scales as N−1/3 for the VO model and as N−2/3 for coil-like molecules. The scaling for the VO model has been reported previously,9 and the scaling for the coil-like molecules is numerically verified in Figure 2. More generally, we find that the critical behavior depends sensitively on the fractal dimension d: for 3/2 < d < 3, we have ϕc ∝ N−1 and σc ∝ N−2/3; for 1 < d ≤ 3/2, we have ϕc ∝ N−(3−d)/d and σc ∝ N−(2d−1)/(2d); for 0 < d ≤ 1, we have ϕc ∝ N−1 and σc ∝ N−(2d+1)/6. These results are summarized in Figure 3. Note that these behaviors are expected for asymptotically large molecular weights. The absolute value of the exponent for σc decreases with d for d < 3/2, suggesting that it is easier to form complexes for molecules of higher fractal dimension. But for d > 3/2, which includes the case of coil-like molecules as an example, the exponent saturates to the value 2/3, so we do not expect significant enhancements to complexation. Although we have presented the scaling exponents for dimension d ranging between 0 and 3, it is worth noting that the relevant physical dimension of macromolecules ranges between 1 and 3. This is natural in that a unidimensional rod is the most extended shape that is possible for a macromolecule. A value between 0 and 1 would have resulted in the so-called Cantor dust.24 In experiments, it is possible to vary molecular fractal dimension systematically.32 In this work, particular attention is paid to the case of d = 0 because it relates to the VO model. In the Appendix, we also show that the critical charge density σc does not depend on the stiffness of the molecular backbone when the fractal dimension is greater than 3/2. It does affect the width of the coexistence window, as we show in the next section. Our key results related to the critical scaling for rodlike and point-like molecules are summarized in Table 1.

OFF-CRITICAL COEXISTENCE CURVE The phase behavior may also be investigated by examining the osmotic pressure. When the coacervate phase is in equilibrium with a supernatant phase, the osmotic pressures of the two phases must be equal. From the free energy expressions eqs 2 and 16, we may evaluate the osmotic pressure by using Πv = ϕf′(ϕ) − f(ϕ). We find 3π ⎛ ⎞ ϕ ⎜ d − 1 csc d + 2 v ⎟ 3/(d + 2) Πv = − ϕ − ln(1 − ϕ) − ⎜ ϕ ⎜ d + 2 12π N λ 0 3 ⎟⎟ ⎝ ⎠ (18)

( )

The first three terms arise from the mixing entropy. The last term arises from the electrostatic correlation free energy. Note that the singularity in the csc term scales as 1/(d − 1) and is compensated by the factor d − 1 in the numerator. For any dimension d other than 1, the last term is nonanalytic. So a virial expansion for the osmotic pressure is not possible, and the electrostatic contribution cannot be absorbed into an effective Flory−Huggins χ parameter. This view should be contrasted with the recent notion of an effective χ presented in the literature.33 This behavior differs from that of polymers in a poor solvent, for which a virial expansion is possible and the osmotic pressure can be written as ϕ 1 Πv = N + 2 (1 − 2χ )ϕ2 + ..., from which a well-defined θ condition χ = 1/2 can be identified. One can appreciate from our analysis that the physics of polyelectrolyte complexation is distinct from that of neutral polymers, and analogies between these two types of materials and their thermodynamic behavior should be regarded with caution. It is possible to derive a convenient result for the coacervate branch of the binodal curve by relying on the osmotic pressure. Equation 19 can be used to evaluate the osmotic pressure for both the coacervate and the supernatant phase. Since the supernatant is nearly depleted of polymers (ϕ = 0), its osmotic pressure Πsupernatant vanishes. Equating eq 18 for the complex with Πsupernatant = 0 leads to the following expression for the coacervate branch of the binodal

Table 1. Scaling Behaviors of Free Energy, Characteristic Microscopic Length a, Correlation Length λ, Critical Volume Fraction ϕc, and Critical Charge Density σc for Complexation of Molecules with Different Fractal Dimension da model

d

free energy

a

λ

ϕc

σc

point-like rod-like coil-like

0 1 2

−ϕ3/2 −ϕ ln(ϕ) ϕ3/4

− v/(πA) b/√12

ϕ−1/2 ϕ−1/3 ϕ−1/4

N−1 N−1/2 N−1

N−1/6 N−1/2 N−2/3

3π ⎛ λ03 ⎜ d − 1 csc d + 2 (binodal) = ⎜ ⎜ d + 2 12π v ⎝

( ) ⎞⎟

×

⎟⎟ ⎠

⎛ϕ ⎞−1 − ϕ − ln(1 − ϕ)⎟ ϕ3/(d + 2) ⎝N ⎠



(19)

Note that λ0 can be converted to charge density σ. This expression works for arbitrary fractal objects and may be used in practice to estimate the polymer concentration in the

a

The off-critical scaling is deduced from eq 16. Note that the VO model is not exactly identical to the case of a point-like model (see text on point-like objects). F

DOI: 10.1021/acs.macromol.6b02113 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Density fluctuations in incompressible fluids are penalized with a high free energy cost. Composition fluctuations are mainly responsible for the correlation energy. To split the contributions from these two terms, we write the m × m matrix U as a summation U(0) + U(1), for which m is the number of components in the system. The term U(0) suppresses the overall density fluctuation and is written as U(0) = Bv uuT, where u = (1, 1, ..., 1) and B is the compressible modulus introduced earlier, which will eventually be set to infinity. The term U(1) represents any additional pairwise interactions, such as Coulombic forces. The interaction kernel in eq 3 can be written as

coexisting coacervate phase. A few representative curves for different fractal dimensions are shown in Figure 4. The results show a clear trend of wider binodal windows for increasing fractal dimension.

det(I + Ω̃U(0) + Ω̃U(1)) = det(I + Ω̃U(0)) det(I + (I + Ω̃U(0))−1Ω̃U(1))

Figure 4. Approximate binodal curves obtained by balancing osmotic pressure, for fractal dimensions d = 2, 3/2, and 1, calculated using eq 19 using N = 100 and a2 = v/(12lB) and using water volume as the reference. The vertical dashed line indicates the supernatant, with vanishing osmotic pressure. The critical regime is shaded. The curve labeled d = 2 is the same as the approximate binodal in Figure 1b.

(20)

The first term captures the effects of composition fluctuations, which arise even in noninteracting (U(1) = 0) fluids. The second term depends on the nature of the interactions and will be explicitly treated later. The composition fluctuation term can be simplified considerably, given the dyadic structure of U(0). Applying the matrix determinant lemma,34 det(A + abT) = (1 + bTA−1a)det(A), we find

The validity of the approximation involved in eq 19 can be appreciated in the right panel of Figure 1 for coil-like molecules. The approximate theory works extremely well as long as the critical regime is avoided. This is mainly because, near the critical point, the correlation length for composition fluctuations becomes much larger than the size of the chains; in that limit, the molecular structural correlation on length scales of the size of the chains (or Rg) becomes relevant, and the Edwards approximation breaks down. Two useful observations can be made from the approximate binodal curve. First, eq 19 exhibits a weak dependence on the molecular weight N. This is due to the fact that in the offcritical region polymer molecules interpenetrate, and the free energy is governed by the structure of a mesh with size much smaller than the polymer radius of gyration, leading to the predicted weakening of the dependence on molecular size. Second, eq 19 implies that the binodal window is narrower for rigid molecules. The stiffness of a molecule is measured by the value of a. Since λ0 is proportional to ad/(d+1)/σ2/(d+2), a fixed value of ϕ results in a greater value of σ for larger a. Conversely, a greater value of σ is needed for the coacervate to coexist with the supernatant, effectively leading to a wider two-phase window. These two observations can also be understood by examining the correlation length, which does not depend on molecular weight and therefore leads to a weak N dependence for the binodal curves. The correlation length grows with molecular stiffness, reducing the electrostatic correlation free energy and requiring a greater charge density to induce phase separation.

det(I + Ω̃U(0)) = 1 + Bu TΩ̃u ≃ Bv u TΩ̃u

(21)

In the last step, we used the fact that B approaches infinity. The strength of the composition fluctuations is controlled by the summation of all entries in the matrix Ω̃. The interaction term depends on the form of the inverse (I + ΩU(0))−1, which can be completed by using the Sherman− Morrison formula,35 (A + abT)−1 = A−1 − ing this result, we find (I + Ω̃U(0))−1 = I −

A −1abTA −1 . 1 + bTAa

BΩ̃uu T Ω̃uu T ≃ − I 1 + Bu TΩ̃u u TΩ̃u

Substitut-

(22)

for incompressible mixtures. Combining the above with the term depending on U(1), we arrive at det(I + (I + Ω̃U(0))−1Ω̃U(1)) ⎛ ⎞ ⎛ Ω̃uu T ⎞ = det⎜⎜I + ⎜I − T ⎟Ω̃U(1)⎟⎟ ⎝ u Ω̃u ⎠ ⎝ ⎠

(23)

We note that both eqs 23 and 21 are valid for a general form of the intramolecular correlations and the pairwise interactions. In the special case when the interaction potential is dominated by electrostatics, the matrix U(1) can be written as a bilinear form, U(1)(q) = ssT, with entries si(q) given by (4πlB)1/2σi/q, where σi is the charge density on the ith species. The value of σi for neutral species is 0. Using this form, eq 23 can be written as



INCOMPRESSIBLE MULTICOMPONENT FLUID The previous sections focus on symmetric polyelectrolyte solutions consisting of a solvent and the polyelectrolytes, in the absence of ions. To analyze the effects of added ions, we start from the generic expression, eq 3, and derive an electrostatic correlation free energy that is analogous to eq 6 but that is completely general and applicable to incompressible mixtures with an arbitrary number of ionic species and with arbitrary molecular architecture, including multiblock polyelectrolytes. In the next section, we apply these results to investigate how added salts alter the complexation behavior.

⎛ ⎛ (u TΩ̃s)2 Ω̃uu TΩ̃s ⎞ T⎞ T ̃ ⎟ det⎜⎜I + ⎜Ω̃s − s 1 s s = + Ω − ⎟ ⎟ ⎝ u TΩ̃u ⎠ ⎠ u TΩ̃u ⎝ (24)

Here we have used the matrix determinant lemma again and used the fact that Ω̃ is symmetric. For symmetric electrolytes, such as solutions of cations and anions that differ only in the sign of the charge, and of polycations and polyanions that differ only in the sign of the charge, it is apparent that uT Ω̃ s G

DOI: 10.1021/acs.macromol.6b02113 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 5. Ternary phase diagram showing the effects of ion addition with connectivity included. The statistical segmental length is b = (v/lB)1/2. The charge density σ is 0.15 for models with connectivity and is 0.20 for the VO modelthe two-phase window in the VO model vanishes for σ = 0.15. Left: two-phase window widens greatly when charge connectivity effects are included. Right: the supernatant phase is nearly fully depleted of polymers; the polymer content in coacervates is well approximated by the results for infinite chains. The greater the N value, the closer the polymer contents to N = ∞. Note that for most polymers 100 repeating units is sufficiently close to N = ∞: the volume of even the smallest monomer, ethylene, is 1.5 times that of a water molecule, resulting in a value of N well above 100. The salt concentration in the shaded regimes is high, which is susceptible to another thermodynamic instability associated with ions salting out, and is relevant for this work since the type of theory developed here is not applicable in the high-salt regime.



EFFECTS OF ADDED IONS We consider a mixture of solvents, polycations, polyanions, cations, and anions. The polycations and polyanions are symmetric, similarly to the systems considered in the previous sections; cations and anions are also symmetric, each of valence 1. The volume fractions of polycation, polyanion, cation, and anion are ϕ/2, ϕ/2, ψ/2, and ψ/2, respectively. By applying eq 25 to this system, we arrive at the following free energy expression v ln(1 + 2v Ω̃ 2ϵσ 2 + 2v Ω̃3ϵ) fcorr = 2 q (26)

vanishes. In these cases, the electrostatic correlation free energy is governed by 1 + sT Ω̃ s. The Gaussian free energy with both osmotic and electrostatic terms can therefore be succinctly written as fcorr =

v 2





∫q ⎢⎢⎣ln(vu TΩ̃u) + ln⎜⎝1 + sTΩ̃s −

(u TΩ̃s)2 ⎞⎤ ⎟⎥ u TΩ̃u ⎠⎥⎦ (25)

where the nonessential B-dependent term has been discarded. The above expression may be applied to block copolyelectrolytes as well. For instance, substituting the intramolecular structure factors for symmetric “ABA + CBC” triblock copolymers to the second term of eq 25, one recovers exactly the key results, eqs 18 and 19 of ref 19. A salient feature of eq 25 is that at the Gaussian level the free energy contributions from the composition fluctuations (the first term) and from the interaction (the second term) are additive. The composition-fluctuations contribution originates from osmotic fluctuations and is already present in neutral solutions, when the interactions are completely turned off. This term has been identified before in neutral polymeric systems29,36 but is rarely discussed in theories of polyelectrolyte solutions. Reference 22 considered this term in the context of an implicit solvent model, and their results are rigorously valid for dilute solutions, since composition-independent excluded volume parameters were included. The effect of these osmotic fluctuations is to regularize the virial coefficients of the free energy densities and has been discussed in previous work on neutral solutions.28,29,36 When needed, for the sake of capturing experimental or simulation results, these effects can be included by introducing the renormalized χ parameter13 or higher order virial coefficients. In this work, for simplicity, we shall focus on the contribution of the electrostatic correlation free energy. Finally, we note that the above expression does not depend on the compression modulus B. All terms depending on the modulus have been either canceled or dropped (ln B) since they only involve a constant shift of the free energy density. This way of “renormalizing” the theory eliminates the unwanted dependence on short-wavelength parameters such as B, which has plagued field theories in the past (e.g., ref 37, eqs 50−53, and discussions therein).



Here, as before, Ω̃2 is the structure factor of the polymers, and ψ Ω̃3 = 2v is that of the ions. On the basis of the result of Figure 1b, which shows that that Edwards approximation works well above the critical charge density, and focusing on the Gaussian chain model, we replace Ω2 by 12/(2vq2Nb2) and find from eq 26 fcorr =

v 2

⎡ ⎛





∫q ⎢⎢⎣ln⎜⎝1 + ϵψ + ϵσ 2ϕ q122b2 ⎟⎠ − ϵψ ⎥⎥⎦

(27)

Here, as in previous sections, we have regularized the free energy by subtracting a term linear in the ion concentration. The above integral can be evaluated analytically to yield 3 1/2 2 ⎛ πl ⎞ fcorr = sgn(Φ − 1) ⎜ B ⎟ (ψ 3 − 3cψϕ + 2c 3/2ϕ3/2)1/2 3⎝ v ⎠

(28)

where we have introduced Φ ≡

σ 2ϕ 3v ψ 2 πlBb2

and c ≡ σ2(3v)/

(πlBb2). Note that the free energy changes sign at Φ = 1. In the low-salt regime, the above expression reduces to the corresponding expression for the binary system, eq 13; in the high-salt regime, the free energy has the exact form of the VO model, eq 8. Alternatively, it can be shown that the free energy can be factored and rewritten as 2(πl 3 / v)1/2

B fcorr = (σ ′ϕ1/2 − ψ )(2σ ′ϕ1/2 + ψ )1/2 , where σ′ is 3 the charge density times a factor representing chain stiffness, 3v defined by σ ′ = σ( πl b2 )1/2 . This expression reduces to eq 47 of

B

ref 10 when proper substitutions are made for the structure factors of connected ions. H

DOI: 10.1021/acs.macromol.6b02113 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 6. Left: s-dependence of J2 integrals for rod-like (d = 1) and coil-like (d = 2) molecules. The low-s slopes are 3/2 in both cases. The high-s slopes are 3/4 for coils, 1 for rods, and 3/(d + 2) for fractal molecules. Right: ratios between J2 and J3 for rod-like and coil-like molecules, demonstrating that the asymptotic exponents of J2 and J3 are identical. Dashed lines indicate asymptotic scalings calculated using prefactors c2, c3, e2, and e3 provided in the text.

of the critical charge density and the polymer composition are found to be solely determined by the fractal dimension of the complexing molecules. Different scaling regimes were identified and, for 3/2 < d < 3, the scaling components were found to remain constant and insensitive to molecular stiffness. We anticipate these results to be useful for actual design of complexing molecules. It should be noted, however, that the Gaussian approach employed in this work breaks down for the supernatant or dilute branch of the coexistence curve. The resulting critical exponents only provide a qualitative guide to the effects of charge connectivity. The validity of this approximation needs to be examined by using a nonperturbative (non-Gaussian) analysis or using well-calibrated simulation results. The fractal dimension affects the complexation behavior considerably, and its effects can be calculated through eq 19. At a qualitative level, such effects can be inferred from the scaling behavior of the critical charge density with molecular weight. We find that more compact objects, with greater d values, have a greater chance to form coacervates. Therefore, we expect that flexible polyelectrolytes are more favorable “coacervate formers” than rod-like molecules. This prediction could be tested in future experiments by designing dendrimers32 with variable generation number. We also derived in this work transparent expressions for the electrostatic correlation free energy of incompressible mixtures including an arbitrary number of ionic species and having arbitrary molecular architectures specified by the structure factor, eq 25. Importantly, this result shows that contributions to the free energy arising from osmotic compositional fluctuations and from electrostatic correlations are additive and are both well-defined in the incompressible limit. The electrostatic correlation free energy interpolates the regimes encompassed in the range between the absence of small ions and the presence of small ions. When this expression is specialized to ternary systems of solvents, polyelectrolytes, and ions, we find that the charge connectivity enhances considerably the stability window of ternary systems (Figure 5). Several important physics, such as charge regulation,7 counterion condensation,38 secondary molecular interactions,39 charge correlation,8 ion pairing,40 excluded volume interaction,13,19 and the possibility of a nematic phase transition for rod-like polymers,14 have been neglected in our theory. Our field-type theory could be further improved by incorporating the fact that charges are distributed discretely, rather than smoothly,20 by adopting a variational approach, rather than a perturbative method,31 and by investigating ternary systems

For systems with added ions, eq 28 smoothly interpolates the free energy expressions for solutions of point-like electrolytes and polyelectrolytes. To examine the effects of charge connectivity in the presence of added ions, in Figure 5 we show binodal coacervation phase diagrams calculated using eq 28 at charge density σ = 0.15 and using the VO model (no charge connectivity) at σ = 0.20. Note that a charge density of 0.20 is used for the VO model because the VO model with a value 0.15 would have predicted a vanishing two-phase window. Charge connectivity leads to three main differences from the VO model. First, the two-phase window widens considerably. Second, the coacervate phase is more stable with respect to the addition of salts. Third, the polymer content in the supernatant phase is much lower, essentially zero. The right panel of Figure 5 also shows that as in the VO model9 the effect of molecular weight N is small. For relevant molecular weights (N ≫100), the polymer content in coacervates is nearly identical to that corresponding to infinite molecular weight chains. This can be used to simplify calculations of phase behavior.



SUMMARY A closed-form expression, eq 16, has been presented for the electrostatic correlation free energy of polyelectrolyte molecules of arbitrary fractal dimension. That expression essentially states that the correlation free energy is about 1 kBT per correlation volume. The theory presented here is insensitive to the short wavelength cutoffin particular, it is insensitive for coil-like moleculesand can be readily used to examine and interpret experimentally observed coacervation behaviors. In the off-critical regime, a convenient analytical expression, eq 19, was also derived to determine the polymer concentration in coacervates phases in equilibrium with a supernatant phase. That expression should facilitate interpretation of experimental measurements, where the norm has been to fit phase diagrams using approximate versions of the VO model, which lacks correlation effects. Our results show that the phase behavior is largely dictated by the correlation length, which grows with chain stiffness and shrinks with the Bjerrum length and polymer concentration. A smaller correlation length leads to higher electrostatic cohesive energy and a wider coexistence region. We stress that the theory developed here is invariant with respect to the choice of reference volume v. Any choice of v can be made for eq 19 and other expressions provided in this work. Using the volume of the solvent as a reference, however, is particularly convenient. In the critical regime, compositional fluctuations at all wavenumbers are relevant. The molecular weight dependencies I

DOI: 10.1021/acs.macromol.6b02113 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules with added salts.9 These all represent possible future directions. We hope that the expressions derived here will motivate more systematic, quantitative experimental studies of phase diagrams, single-chain structure, and collective scattering, which would provide the ultimate test of the assumptions invoked in our theory. Similarly, additional computer simulations41−43 of polyelectrolyte structure and phase behavior would be particularly helpful in that regard.

These values determine the prefactors but not the exponents in the N-dependence of critical properties. Numerical tests of these scaling behaviors are provided in Figure 6. The ratio J2(s)/J3(s) approaches 4 in the small-s regime and approaches a d-dependent constant that is greater than e2/e3 = 1.43 for d = 3 in the large-s regime. To find the critical point, we require that the second and third derivatives of the total free energy vanish. Merging the contributions from mixing entropy and electrostatic correlation free energy, we obtain



APPENDIX. SCALING PROPERTIES NEAR THE CRITICAL POINT We consider the free energy of fractal objects of dimension d. Their intramolecular correlation function is denoted by ϕN Ω̃(q) = 2v g (qR g). Here g(qRg) is a generalized Debye function with the following properties: in the low-q regime, g(qRg) = 1 − (qRg)2/3; in the high-q regime, g(qRg) ≃ (qRg)−d. The radius of gyration can be denoted by Rg = N1/da, with a being a monomeric length scale. Equation 14 mainly captures the high-q behavior of the Debye function. The low-q properties are needed in order to analyze the critical asymptotics. The second- and third-order differentials of the electrostatic correlation free energy can be written in the following compact form: (d) ∂ 2f corr 2

=−

∂ϕ

(d) ∂ 3f corr

∂ϕ3

1 1 v N −3/ dϕ−2J2 (s) + = 2 3 Nϕ 1−ϕ 4π a ϕ 1 v N −3/ dϕ−2J3(s) − = Nϕ (1 − ϕ)2 2π 2a3

Solution of the above conditions yields the critical values of ϕc and sc. Although the behavior of J2(sc) and J3(sc) is not known at this stage, their ratio is bounded between 1 and 4 as shown in Figure 6. So we introduce hc ≡ J2(sc)/(2J3(sc)) and take the ratio of the two equations in eq 32. The result can be simplified and reads 1 + (hc − 1)ϕc hc − 1 = Nϕc (1 − ϕc)2

ϕc ≃ (hc − 1)N −1

s≡

∫0



⎛ g (s1/2Q ) ⎞m ⎟⎟ , dQ Q 2⎜⎜ 2 1/2 ⎝ Q + g (s Q ) ⎠

4πlBa 2 1 + 2/ d 2 ϕσ N v

Here, the nonessential constant 1/(2π2) has been discarded. The magnitude of the term on the right-hand side clearly depends on the value of d. We discuss, in the following, the possible scenarios separately. Regime a: d > 3/2. The term N3/d−2 becomes vanishingly small for sufficiently large N, so we have J3(sc) ≪ 1, which implies a vanishingly small sc. Figure 6 shows that in this regime J3(s) behaves as sc3/2 ∝(lBa2/v)3/2N3/d(hc− 1)3/2σ3. So we further have

(30)

Note that Q in eq 30 represents qRg/s1/2. We transformed the original expression to this form so that the effects of all the controlling parameters can be considered collectively. Before moving on to the discussion of the spinodal and the critical point, we note that Jm(s) has the following asymptotic properties. For s → 0, g(s1/2Q) approaches 1 and Jm(s) behaves 2 as cms3/2, where the prefactor cm is given by ∫ ∞ 0 dQ Q (1 + 2 −m 1/2 −d/2 −d Q ) . For s → ∞, g(s Q) approaches βds Q and Jm(s) behaves as ems3/(2+d), where βd is a geometric constant corresponding to the fractal object and the prefactor em is 2 2+d −m given by βd3/(2+d)∫ ∞ 0 dQ Q (1 + Q ) . The parameter βd is 2 for coil-like and π for rod-like molecules. To get the last expression for em, we have rescaled Q by sd/(4+2d). It can be shown that c2 = π/4, c3 = π/16, and

⎞1/6 ⎛ v ⎟ N −2/3 σc ≃ ⎜ 3 ⎝ lB (hc − 1) ⎠

1 3/(2 + d) ⎛ 3(d + 1) ⎞ ⎛ d + 5 ⎞ ⎟ ⎟Γ⎜ β Γ⎜ ⎝ d + 2 ⎠ ⎝d + 2⎠ 6 d

(36)

The critical charge density scales with N−2/3. Another important point is that the dependence on backbone stiffness a drops completely. Recall that in order to obtain the solution ϕc, we have required that hc > 1. This is indeed satisfied for sc ≪ 1, so the results are self-consistent. Moreover, Figure 6 shows that hc → 2 in this regime, so we can neglect the hc dependence and write ϕc ≃N−1 and σc ≃N−2/3. Regime b: d ≤ 1. The term N3/d−2 becomes asymptotically large, so we have J3(sc) ≫ 1 and sc ≫ 1. But in this regime J3(sc)

⎛ 3π ⎞ d − 1 ⎟ e 2 = πβd 3/(2 + d) csc⎜ ⎝d + 2⎠d + 2 e3 =

(34)

The critical composition scales with N−1, irrespective of the fractal dimension (however, see eq. 39 below). Further, it is evident that the solution exists only for hc > 1, which constrains the range of sc. Substituting eq 34 back to eq 32, we infer that both sides of the two equations need to be of order 1, and the solution exists only if we satisfy v J (sc) ≃ (hc − 1)N3/ d − 2 (35) a3 3

(29)

Here we introduced the auxiliary functions J2(s) and J3(s), whose argument s depends on both volume fraction ϕ and charge density σ. They both are derived from a more general auxiliary Jm(s) function defined by Jm (s) ≡ s 3/2

(33)

Since ϕc ≪ 1 and hc are bounded, the right-hand is of order unity. So we have

v N −3/ dϕ−2J2 (s) 4π 2a3

v N −3/ dϕ−3J3(s) = 2 3 2π a

(32)

(31) J

DOI: 10.1021/acs.macromol.6b02113 Macromolecules XXXX, XXX, XXX−XXX

Macromolecules



scales as sc3/(2+d) ∝ (N2/d (hc − 1)σc2)3/(2+d). Comparing this with eq 35, we get (d − 1)/3

σc ≃ (hc − 1)

N

−(1 + 2d)/6

(37)

(38)

Alternatively, it can be written ϕc ≃ N −(3/ d − 1)J3(sc)

(39)

The value of sc approaches 0 in regime a and ∞ in regime b. What happens in regime c is that sc approaches a fixed constant. We have numerically verified this for d = 3/2 and for several values of d sampled between 1 and 3/2. Further, the value of sc increases with decreasing d and becomes ∞ at d = 1. In case of d = 3/2, J3(sc) is finite. So eq 39 represents the actual scaling, i.e., ϕc ≃ N−1, which is compatible with the scaling in regime a. In case of d → 1, the apparent scaling N−2 is compensated by the diverging factor J3(sc) ∝ sc. The scaling of σc follows from the fact that sc is fixed and from the scaling of ϕc. From sc ∝N1+2/dϕcσc2, we deduce that ⎛ s ⎞1/2 σc ≃ ⎜⎜ c ⎟⎟ N −(2d − 1)/(2d) ⎝ J3(sc) ⎠

(40)

Note that since J3(sc) ≃ sc for d → 1, the above expression represents a well-defined scaling rule, σc ≃ N−(2d−1)/(2d). The exponent becomes −3/2 for d = 3/2 and −1/2 for d = 1, smoothly interpolating the values in regimes a and b.



REFERENCES

(1) de Jong, H. G. B.; Kruyt, H. R. Coacervation: partial miscibility in colloid systems. Proc. Koninkl. Med. Akad. Wetershap. 1929, 32, 849− 856. (2) van der Gucht, J.; Spruijt, E.; Lemmers, M.; Cohen Stuart, M. A. Polyelectrolyte complexes: bulk phases and colloidal systems. J. Colloid Interface Sci. 2011, 361, 407−422. (3) Wang, Q.; Schlenoff, J. B. The polyelectrolyte complex/ coacervate continuum. Macromolecules 2014, 47, 3108−3116. (4) Spruijt, E.; Westphal, A. H.; Borst, J. W.; Cohen Stuart, M. A.; van der Gucht, J. Binodal compositions of polyelectrolyte complexes. Macromolecules 2010, 43, 6476−6484. (5) Overbeek, J.; Voorn, M. Phase separation in polyelectrolyte solutions: theory of complex coacervation. J. Cell. Comp. Physiol. 1957, 49, 7. (6) Fowler, R.; Guggenheim, E. A. Statistical Thermodynamics, 2nd ed.; Cambridge University Press: 1949. (7) Muthukumar, M.; Hua, J.; Kundagrami, A. Charge regularization in phase separating polyelectrolyte solutions. J. Chem. Phys. 2010, 132, 084901. (8) Perry, S. L.; Sing, C. E. PRISM-based theory of complex coacervation: excluded volume versus chain correlation. Macromolecules 2015, 48, 5040−5053. (9) Qin, J.; Priftis, D.; Farina, R.; Perry, S. L.; Leon, L.; Whitmer, J.; Hoffmann, K.; Tirrell, M.; de Pablo, J. J. Interfacial tension of polyelectrolyte complex coacervate phases. ACS Macro Lett. 2014, 3, 565. (10) Borue, V. Y.; Erukhimovich, I. Y. A statistical sheory of weakly charged polyelectrolytes: fluctuations, equation of state, and microphase separation. Macromolecules 1988, 21, 3240. (11) Borue, V. Y.; Erukhimovich, I. Y. A statistical theory of globular polyelectrolyte complexes. Macromolecules 1990, 23, 3625−3632. (12) Mahdi, K. A.; Olvera de la Cruz, M. Phase diagram of salt-free polyelectrolyte semidilute solutions. Macromolecules 2000, 33, 7649− 7654. (13) Kudlay, A.; Ermoshkin, A. V.; Olvera de la Cruz, M. Complexation of oppositely charged polyelectrolytes: effect of ion pair formation. Macromolecules 2004, 37, 9231. (14) Kumar, R.; Audus, D.; Fredrickson, G. H. Phase separation in symmetric mixtures of oppositely charged rodlike polyelectrolytes. J. Phys. Chem. B 2010, 114, 9956−9976. (15) Sing, C. E. Development of the modern theory of polymeric complex coacervation. Adv. Colloid Interface Sci. 2016, DOI: 10.1016/ j.cis.2016.04.004. (16) de Gennes, P. G. Scaling Concepts in Polymer Physics; Cornell University Press: 1979. (17) Castelnovo, M.; Joanny, J. F. Complexation between oppositely charged polyelectrolytes: beyond the random phase approximation. Eur. Phys. J. E: Soft Matter Biol. Phys. 2001, 6, 377. (18) Borue, V. Y.; Erukhimovich, I. Y. Weak crystallization and structural phase transitions in weakly charged polyelectrolyte systems. Sov. Phys. JETP 1991, 72, 751. (19) Audus, D. J.; Gopez, J. D.; Krogstad, D. V.; Lynd, N. A.; Kramer, E. J.; Hawker, C. J.; Fredrickson, G. H. Phase behavior of electrostatically complexed polyelectrolyte gels using an embedded fluctuation model. Soft Matter 2015, 11, 1214−1225. (20) Potemkin, I. I.; Palyulin, V. V. Complexation of oppositely charged polyelectrolytes: effects of discrete charge distribution along the chain. Phys. Rev. E 2010, 81, 041802. (21) Popov, Y. O.; Lee, J.; Fredrickson, G. H. Field-theoretic simulations of polyelectrolyte complexation. J. Polym. Sci., Part B: Polym. Phys. 2007, 45, 3223. (22) Lee, J.; Popov, Y. O.; Fredrickson, G. H. Complex coacervation: a field theoretic simulation study of polyelectrolyte complexation. J. Chem. Phys. 2008, 128, 224908. (23) Riggleman, R. A.; Kumar, R.; Fredrickson, G. H. Investigation of the interfacial tension of complex coacervates using field-theoretic simulations. J. Chem. Phys. 2012, 136, 024903.

The critical charge density scales as N−(1+2d)/6. The selfconsistency check is ensured by noting that hc > 1 for sc ≫ 1. Regime c: 1 < d ≤ 3/2. For the particular case of d = 3/2, eq 35 implies that J3(sc) and sc reach a constant. For the cases 1 < d < 3/2, eq 35 implies that J3(sc) and sc become asymptotically large with N, analogous to the case of d < 1. However, the argument breaks down because Figure 6 shows that the ratio hc = J2(sc)/2J3(sc) is less than 1 for sc ∝ ∞ and for d > 1, which contradicts the presumption of hc > 1 behind eq 34. To determine the correct scaling relation, we examine eq 32 more closely. We first notice that the J3- and J2-dependent terms are of the same order: the ratio J2/J3 is bounded and the prefactors are identical. Then we note that J3 cannot be of the same order as or lower order than ϕ/(1 − ϕ)2. Had it been the case, J2 would be of lower order than 1/(1 − ϕ), and the first equation could not be balanced. So we conclude that J3 is of higher order than ϕ/(1 − ϕ)2. To balance the second equation, we require that while neglecting the nonessential prefactors 1 ≃ N −3/ dϕc−2J3(sc) Nϕc

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail [email protected] (J.Q.). Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work is supported by the Department of Energy, Basic Energy Sciences, Division of Materials Research and Engineering. K

DOI: 10.1021/acs.macromol.6b02113 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules (24) Mandelbrot, B. B. The Fractal Geometry of Nature; W.H. Freeman and Company: 1983. (25) Neugebauer, T. Berechnung der Lichtzerstreuung von Fadenkettenlösungen. Ann. Phys. 1943, 434, 509. (26) Morse, D. C. Diagrammatic analysis of correlations in polymer fluids: cluster diagrams via Edwards’ field theory. Ann. Phys. 2006, 321, 2318−2389. (27) Fredrickson, G. H. The Equilibrium Theory of Inhomogeneous Polymers; Oxford University Press: 2006. (28) Qin, J.; Morse, D. C. Renormalized one-loop theory of correlations in polymer blends. J. Chem. Phys. 2009, 130, 224902. (29) Olvera de la Cruz, M.; Edwards, S. F.; Sanchez, I. C. Concentration fluctuations in polymer blend thermodynamics. J. Chem. Phys. 1988, 89, 1704. (30) Landau, L. D. Statistical Physics, Part I, 3rd ed.; ButterworthHeinemann: 1980; section 78. (31) Wang, Z.-G. Fluctuation in electrolyte solutions: the self energy. Phys. Rev. E 2010, 81, 021501. (32) Rathgeber, S.; Gast, A. P.; Hedrick, J. L. Structural properties of star-like dendrimers in solution. Appl. Phys. A: Mater. Sci. Process. 2002, 74, S396−398. (33) Spruijt, E.; Sprakel, J.; Cohen Stuart, M. A.; van der Gucht, J. Interfacial tension between a complex coacervate phase and its coexisting aqueous phase. Soft Matter 2010, 6, 172. (34) Harville, D. A. Matrix Algebra from a Statistician’s Perspective; Springer-Verlag: 1997. (35) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes; Cambridge University Press: 2007. (36) Fredrickson, G. H.; Liu, A. J.; Bates, F. S. Entropic corrections to the Flory-Huggins Theory of polymer blends-architectural and conformational effects. Macromolecules 1994, 27, 2503−2511. (37) Delaney, K. T.; Fredrickson, G. H. Recent developments in fully fluctuating field-theoretic simulations of polymer melts and solutions. J. Phys. Chem. B 2016, 120, 7615−7634. (38) Manning, G. S. Limiting laws and counterion condensation in polyelectrolyte solutions I. colligative properties. J. Chem. Phys. 1969, 51, 924. (39) Perry, S. L.; Leon, L.; Hoffmann, K. Q.; Kade, M. J.; Priftis, D.; Black, K. A.; Wong, D.; Klein, R. A.; Pierce, C. F., III; Margossian, K. O.; Whitmer, J. K.; Qin, J.; de Pablo, J. J.; Tirrell, M. Chirality-selected phase behaviour in ionic polypeptide complexes. Nat. Commun. 2015, 6, 6052. (40) Salehi, A.; Desai, P. S.; Li, J.; Steele, C. A.; Larson, R. G. Relationship between polyelectrolyte bulk complexation and kinetics of their layer-by-layer assembly. Macromolecules 2015, 48, 400−409. (41) Ou, Z.; Muthukumar, M. Entropy and enthalpy of polyelectrolyte complexation: Langevin dynamics simulations. J. Chem. Phys. 2006, 124, 154902. (42) Hoda, N.; Larson, R. G. Explicit- and implicit-solvent molecular dynamics simulations of complex formation between polycations and polyanions. Macromolecules 2009, 42, 8851−8863. (43) Hoffmann, K. Q.; Perry, S. L.; Leon, L.; Priftis, D.; Tirrell, M.; de Pablo, J. J. Chirality in charge-driven polypeptide complexation. Soft Matter 2015, 11, 1525−1538.

L

DOI: 10.1021/acs.macromol.6b02113 Macromolecules XXXX, XXX, XXX−XXX