Effect of Ion–Ion Correlations on Polyelectrolyte Gel Collapse and

Jun 4, 2013 - Ion Correlation-Induced Phase Separation in Polyelectrolyte Blends. Charles E. Sing , Jos W. Zwanikken , and Monica Olvera de la Cruz. A...
0 downloads 10 Views 3MB Size
Article pubs.acs.org/Macromolecules

Effect of Ion−Ion Correlations on Polyelectrolyte Gel Collapse and Reentrant Swelling Charles E. Sing, Jos W. Zwanikken, and Monica Olvera de la Cruz* Department of Materials Science and Engineering, Northwestern University, Evanston, Illinois 60208, United States ABSTRACT: Polyelectrolyte gels are known to exhibit strong response to changes in salt concentration, an effect attributed to osmotic effects. Traditional theories are insufficient to describe all the properties observed in experiments; namely, reentrant swelling effects at high salt concentrations and complete polyelectrolyte collapse at high salt valencies. We incorporate liquid-state integral equation theories to account for the many-body correlations between finite-sized ions in a strongly charged polyelectrolyte gel. This represents a theory for such systems that goes well beyond well-known Debye−Huckel (DH) and Poisson− Boltzmann (PB) theories. We demonstrate that charge correlations give rise to enhanced gel deswelling, a first-order collapse to a highly correlated state at large valencies, and counterion size-dependent reentrant swelling. Such ion-specific effects provide routes for practical control of stimuli-responsive polymeric materials based on polyelectrolyte gels, and deeper understanding of biological systems.



INTRODUCTION Polyelectrolyte gels are ubiquitous in both biological1−8 and synthetic9−14 polymer science due to their physical response and chemical versatility. It is imperative, therefore, that the physics of polyelectrolyte gels is fully understood in both a conceptual and rigorous manner. There is a long history of study on polyelectrolyte gels and related systems, with pioneering work in the first part of the twentieth century.15 These physical descriptions consider gel swelling as primarily driven by counterion entropy, wherein the polymer phase will expand to increase the volume accessible to maximize the entropy.15−18 Powerful theoretical frameworks provide a basis to study structural characteristics such as the inhomogeneities inherent to gelated polyelectrolytes,18 as well as mechanical instabilities that arise due to the presence of poor-solvent conditions.18,19 These results are largely confirmed by experimental observation of structural characteristics,20 and recently an abundance of similar effects have been studied in great detail.20,21 Nevertheless, the Coulombic correlations between counterions in polyelectrolyte gels have been largely ignored except in the context of Poisson−Boltzmann (PB) or Debye−Huckel (DH) assumptions which provide the correct limiting behavior22−24 but completely neglect the microscopic ion structure and its relevant effects.25−27 Results from a number of experimental investigations have frequently shown these traditional theories to be insufficient, especially when considering multivalent ions or extremely dense ionic situations with a significant local structure.13,28,29 Many of these effects are reproduced in simulation investigations of polyelectrolyte gels, reinforcing the importance of correlations when considering multivalent ions.30−33 Some theoretical insights have enabled the description of these behaviors in analogous nongelated polyelectrolyte solutions, © 2013 American Chemical Society

through the judicious use of approximations governing the possible states of these systems and the consideration of the strong-coupling limit.33−35 Nevertheless, a first-principles description of polyelectrolytes and polyelectrolyte gels has so far eluded the field. Indeed, in many strong coupling theory applications, such as the description of adsorbed polyelectrolytes onto oppositely charged surfaces,36 it is assume that the counterions and salt ions are point-like particles, and this approximation is carried out even in rigorous implementation of the strong coupling theory.37−39 Therefore, the approach fails to account for many body effects that depend on the size of the salt ions and counterions that explain phenomena such as reentrant transitions (a re-expansion of collapsed chains as the concentration of salt increases) and possible charge inversion as the salt concentration increases, as demonstrated by strongcoupling models that phenomenologically include finite ion sizes,35,40 which have been corroborated by simulations in saltfree polyelectrolytes,41 and in polyelectrolyte gels in the presence of salt.33 Liquid-state integral equation theories represent one of the most powerful ways of self-consistently calculating the full partition function, as it does not rely on the first few terms of a series expansion as is typical in a perturbation theory.25,42 Rather, the diagrammatic representation of different correlation functions can be directly related to each other to the sole exclusion of a small set of ″bridge diagrams″.25,42 Therefore, this approach can readily account for an infinite series of Mayer cluster diagrams in a nearly complete fashion.25,42 Recently, an accurate integral equation was developed to describe electrolyte Received: February 20, 2013 Revised: May 14, 2013 Published: June 4, 2013 5053

dx.doi.org/10.1021/ma400372p | Macromolecules 2013, 46, 5053−5065

Macromolecules

Article

solutions which permits the incorporation of this type of theoretical framework into the understanding of polyelectrolyte gels in such a fashion that other theoretical results can be assessed or clarified, and new physical phenomena can possibly be revealed.25 This allows for the treatment of ion−ion correlations that we demonstrate give rise to reentrant swelling at high salt concentrations, suppressed swelling at low concentrations, and multivalent ion-driven collapse at intermediate concentrations. These behaviors are widely observed in experimental situations concerning polyelectrolyte solutions, and arise even in a “primitive model” of gels that neglects dispersive or hydrophobic solvent effects, local correlations due to chain connectivity, and dielectric effects. These approximations, as well as the difficulties of incorporating rigorously the gel mechanical response to swelling and compaction, hinder quantitative matching to experiment and simulation, and we discuss their ramifications and possibilities for incorporation into the current framework. However, even with these simplifications qualitative comparison can be made to analogous experiments.28,29,43

Figure 1. Schematic showing the system described in this investigation. Two phases are considered, a solution phase (right) and a gel phase (left). The solution phase contains a number concentration ρS of a symmetric salt with charge q± that dissociates into ions of radius ai. The gel phase, which can expand from a dry volume V0 to a volume V contains a cross-linked polyelectrolyte with nx cross-links and a fraction f of positively charged monomers. These positive charges have associated counterions that are identical to the negative charges in the solution. The counterions of the polymer and the ions from the salt in solution are free to move across the boundary between the phases, such that there is a concentration of positive (ρ+,G) and negative (ρ−,G) ions in the gel phase.



THEORY The overall formulation of our numerical theory involves the articulation of liquid-state integral equation theory within macroscopic thermodynamics in such a way that structure of the ions is manifested on a larger scale through statistical mechanical connections. Therefore, this hybrid approach elucidates properties missing in the current theories.15−18,23,44 We describe both the theory at each length scale and subsequently the self-consistent scheme we use to fully describe the system. Gel Thermodynamic Parameters. The macroscopic thermodynamic definition of our system relies on modifying the classical Donnan model, with the eventual incorporation of the liquid state theory through the calculation of the chemical potential and pressures from local correlations. The classical Donnan model is widely used as the fundamental concept upon which the current understanding of polyelectrolyte gels is based;15,28,29,43,44 namely the osmotic pressure of the counterions that are localized in the gel due to charge neutrality drives the gel to swell as a way of maximizing the translational entropy of the ions.15,28,29,43,44 This is counteracted by the elasticity of the gel, which is dependent on the physical structure of the polyelectrolyte network.18 Figure 1 demonstrates the essential components of our model, which will consider two phases (subscripts S and G as the solution and gel respectively) that are open to the transport of unbound salt ions across the interface. The solution phase S consists of only an aqueous (relative dielectric constant, εr = 80.1) solution of a q+:−q− salt where qi is the valency of species i (+ for the positive, − for the negative ion). Both ions are free to cross in and out of the system by leaving or entering the gel, however we assume that the solvent phase is large enough compared to the gel phase that it maintains a constant number concentration of ρS(#/nm3) = ρ+,S = ρ−,S. For this paper we restrict ourselves to the case of symmetric ions q+ = −q− = q± for simplicity, and anticipate the study of asymmetric systems in the near future since there are many predictions of novel effects in these situations. For both the S and G (gel) phases, the open nature of these two phases indicates the appropriate potential to use is the grand potential Ω̃S, where the tilde denotes that the energy is rendered

dimensionless by the thermal energy kBT. For the S phase, we write: ̃ Ω̃S = -̃IG , S + -CORR ,S −

∑ μĩ ,S ni ,S

(1)

i

where ni,S and μ̃i,S are the number and chemical potential of ⎡ n Λ ⎤ + 1⎥ is the ideal gas species i in S. -̃ IG , S = ∑i ni , S⎢ln iV,S ⎣ ⎦ S free energy in a solution volume of VS and de Broglie ̃ wavelength Λ. -CORR , S is the free energy associated with liquid state correlations in the solution phase, and will be calculated through the liquid state theory integral equations.42 The gel phase G consists of an aqueous solution of the same q± salt that is in the solution phase, since the boundaries between the two systems are open. Furthermore, there are charges that are due to the presence of the polymer. We choose to consider only a strongly disassociating polymer, such as poly(styrenesulfonate) or DNA, that has a constant fraction f of charged monomers. The charged monomers and their counterions are identical to the salt species in the solution phase S, however the charged monomers are restricted to the gel phase G. The concentrations of each species is denoted by ρ+,G, ρ−,G, and ρP,G where the charged monomers are given a positive charge and the electroneutrality condition ρ+,G + ρP,G = ρ−,G is maintained. There is a grand potential Ω̃G for the gel phase that is given by

( )

̃ Ω̃G = -̃ IG , G + -̃ EL , G + -CORR ,G −

∑ μĩ ,G ni ,G i

(2)

̃ where -̃IG , G and -CORR , G are calculated the same as in the solution phase, and at equilibrium across the S−G boundary μ̃i,G = μ̃ i,S = μ̃i. There is an extra term -̃ EL , G that corresponds to the free energy of the elastic network as it changes volume V from its dry volume V0: 5054

dx.doi.org/10.1021/ma400372p | Macromolecules 2013, 46, 5053−5065

Macromolecules -̃ EL , G

Article

2/3 ⎡ ⎛ ⎡ ⎛ V ⎞2/3⎤⎤ 1 V ⎞ ⎢ ⎢ = Cnx ⎜ ⎟ − ln 1 − ⎜ ⎟ ⎥⎥ ⎢ 2 ⎝ Vmax ⎠ ⎢⎣ V ⎝ ⎠ ⎥⎦⎥⎦ max ⎣

where X = πV′/(12Nx2) and Nx is the degree of polymerization between cross-links, and aM is the monomer radius. P̃CORR values are pressures and the μ̃ CORR values are the chemical potentials associated with the ion−ion correlations, and are calculated from liquid state integral equations. The general form of these equations is F1({ρi}, V′) = F2({ρ′i }), and for the solution of these equations to be self-consistent {ρi} = {ρi′}. We therefore iteratively solve these equations by choosing test values {ρi} and V′ and solving directly for F2. The equation F1({ρi}, V′) = F2 is solved for the arguments of F1 where F2 is held constant, and the resulting values of {ρi} and V′ are used to solve once again directly for F2. F2 represents the correlation values Δμ̃CORR and ΔP̃ CORR, which are determined through liquid state theory, so the above thermodynamic equations only provide one-half of this iterative scheme. Liquid State Theory for Polyelectrolyte Solutions. To calculate the values for Δμ̃CORR and ΔP̃ CORR, we use the liquid state theory with the hard core−hypernetted chain (HC− HNC) approach.25 The HC−HNC approach requires the selfconsistent solution of the Ornstern−Zernike equation:25,42

(3)

where nx is the number of cross-links in the gel, Vmax = 8nxN3a3 is the maximum volume of the gel with monomer (and ion) radius a(= a±) = 0.8nm and degree of polymerization between cross-links N, and C is a phenomenological parameter that allows us to tune the overall stiffness of the gel. The ion radius we use is larger than typical hydrated ions,45 but we ultimately will show how the observed effects extrapolate to smaller ion sizes. We consider a highly stiff gel (C = 340, N = 100) in order to probe relatively dense gelated systems. This form of the elastic free energy represents an approximate version of the inverse Langevin representation of polymer elastic response, in line with the work by Jha, et al.44 There is also a hard cutoff that prevents the gel from shrinking below V/V0 = V′ = 1.0, and represents the excluded volume constraints of the noncharged polymer components that are not explicitly described in our integral equation theory. In this work we ignore contributions to both dispersive interactions between the various components (traditionally represented through the Flory χ parameter) and dielectric self-energy effects. These are short-range effects that, while they could be incorporated readily into the theory, only serve to expand the parameter space needlessly as our primary interest is the multi body correlations in the electrostatic interactions. Indeed, such short-range considerations are crucial for understanding temperature dependence of gels, and correlations would likely enter the conceptual picture.46−48 Existing literature allows for the conceptual extrapolation of these effects from our current constant dielectric case without short-range attractions, and would affect the magnitude and location of the behaviors investigated in this article.18,48,49 The equilibrium solution of these equations requires the simultaneous minimization of the free energy along the independent coordinates of gel counterion concentration ρ−,G and gel swelling ratio V′, which corresponds to chemical and mechanical equilibrium, respectively. The appropriate equations are ∂Ω̃G ∂Ω̃ =− S ∂V ∂V

(4)

∂Ω̃G ∂Ω̃S = ∂n± ∂n±

(5)

ĥ ij = c ̂ij + ρk c ̂ikĥ kj

with the HC−HNC closure relation25,42 c ij = −β u ij + h ij − ln(h ij + 1) rij > dij h ij = −1.0

β u ij =

(6)

and the mechanical equilibrium results in the relationship ⎛ 3X −1/3 − X1/3 ⎞ C ⎟ ⎜ 3 24(NxaM ) ⎝ 1 − X2/3 ⎠

̃ ̃ ̃ ̃ = (PCORR , +, G + PCORR , −, S) − PCORR , +, S + PCORR , −, S ̃ = ΔPCORR

(9)

qiqjlB rij

(10)

where qi is the valency of species i and lB = e2/(4πkBTε) is the Bjerrum length in a medium of dielectric constant ε. We have written these equations in tensor form, since we consider two separate components (+ and −) for each phase. In G, + and P are identical and thus calculations are carried out with no distinction between the two components. It has been shown that this method accurately reproduces correlation functions calculated by simulation for ionic solutions, with deviations only at high densities and low temperatures that are well outside the regimes represented in this work.25,50 We make the assumption that the along-the-chain connectivity (P − P) does not play a primary role in the gel thermodynamics, an assertion we examine later in comparison to simulation results in the literature. We expect that the primary effect of connectivity would be conformational and therefore would renormalize the gel elasticity, however we put this into the model at a phenomenological level.51 Previous work on electrolytes that incorporate connectivity using PRISM (a related integral equation method) to look at solutions, albeit with very different motivations and assumptions, but these results make no mention of whether or not the presence of chain connectivity appreciatively alters the surrounding ion correlations in a semidilute gel phase.52 It is unclear to what extent such connectivity considerations alter the gel thermodynamics, though it would be reasonable to expect some amount of enhancement of the local correlation effects; however, future

⎡ρ ρ ⎤ +, G −, G ⎥ = μCORR , +, S + μCORR , −, S − μCORR , +, G ln⎢ ⎢⎣ ρS2 ⎥⎦

2ρS − ρ+, G − ρ−, G +

rij < dij

where the hij and cij are the indirect and direct correlation functions respectively between species i and j, a hat denotes the Fourier transform of the correlation function, rij is the distance between i and j, dij is the hard-core distance between i and j, and uij is the Coulombic potential between species i and j:

where ∂n± corresponds to the movement of an ion pair rather than a single ion such that electroneutrality is maintained. Chemical equilibrium results in the relationship:

+ μCORR , −, G = ΔμCORR

(8)

(7) 5055

dx.doi.org/10.1021/ma400372p | Macromolecules 2013, 46, 5053−5065

Macromolecules

Article

work will be needed to test the magnitude of such effects and subsequently the validity of our hypothesis in a rigorous fashion. Since this work uses the HNC approximation, we can directly calculate the excess chemical potential due to the ion− ion correlations calculated by the equation42 μCORR = ̃ ,i =

̃ ∂-CORR V ∂ni ρj

∑ j

∫ 2 0



h ij[h ij − c ij]d r − ρj

∫0



c ij dr (11) Figure 2. V′ versus ρS for a f = 0.05,0.1,0.2, and 0.5 for a monovalent q± = 1.0 salt solution. Lines demonstrate the result of the Donnan theory D that neglects ion correlations, hard-sphere excluded volume, and electrostatics except through electroneutrality constraints. We use liquid-state integral equation theory to incorporate all of these effects in our HC−HNC scheme, resulting in the connected symbols. HC− HNC results demonstrate both an enhanced deswelling at low ρS and a reentrant swelling behavior at large ρS.

This represents one of two connections of the integral equation theory to the macroscopic thermodynamic variables. The other is the excess pressure due to the ion−ion correlations:42 ̃ PCORR =

̃ ∂-CORR

(1 − V ′−1)∂V ρj ρi ∂ũ ij (hij + 1.0) dr = −∑ −1 ∂r i , j 6(1 − V ′ ) ρj ρi (hij(aij) + 1.0)πaij 3 + 6(1 − V ′−1)



experimental literature on polyelectrolyte solutions.54−57 Traces of V′ as a function of ρS are shown in the context of both the D case (lines) and the HC−HNC case (symbols), and these contrasting situations demonstrate the full effect of the correlations by completely turning on and off the liquid-state theory component of the calculation. We note that this comparison represents the difference between a completely noncorrelated system and a near-complete description of correlations, with Poisson−Boltzmann and related calculations (that often put correlations in via perturbative or ad-hoc scaling arguments) representing a middle ground between these two extremes. We comment later on the limitations of these approaches, while for now we focus on the complete presence or absence of correlations. The non-negligible contribution of these effects are striking under all values of f, and a few pronounced characteristics are apparent. Most importantly, at low ρS the D and HC−HNC converge to the same trend, reaffirming that the D model does indeed provide the correct limiting behavior. Even at small values of ρS, however, there is an enhanced deswelling upon increase of ρS once correlations are included. This leads to a significant plateau region of minimal swelling, and ultimately to a reentrant swelling behavior at large ρS that is not observed at all in the traditional D model. These fundamental features are present at all values of f, however f importantly changes the degree of swelling in the gel. This is due to the increase in bound counterions upon a concomitant increase in f that leads to an increased entropic driving force toward the swollen state. In many polymers of technological interest, such as poly(styrenesulfonate), the entropic driving force toward the swollen state is mitigated somewhat by a large value of the effective χ-parameter of the polymer due to the inherent hydrophobicity of the hydrocarbon chain. This effect is not included in our work, but significant work has been done in understanding its ramifications through both simulation and theory.18,48 Here we focus on the effects of hard core and electrostatic correlations, without considering hydrophocity, which we expect to investigate further in the future. We identify three aspects of these curves that we will investigate in detail, of which two are already apparent and a

(12)

Here aij is the radial distance between the centers of i and j, and the factor (1 − V′−1) is an implicit correction for the actual volume available to the ions due to the excluded volume of the polymer. The first term corresponds to the repulsion and attraction due to the electrostatic interactions, and the second term is the portion of the pressure due to the hard core repulsions. A key advantage of the liquid state theory is that it is possible to “turn on” or “turn off” any given portion of the potential in the liquid. It is in principle possible to incorporate short-range dispersive or hydrophobic interactions at this level through additions to the potential uĩ j, though use of alternate closure approximations may be necessary.42 This would complement appropriate changes to the phenomenological elastic free energy of the underlying gel that include Flory χparameter contributions.53 Throughout the remainder of this work, we refer to three different situations. The pure Donnan theory (ΔμCORR = 0, ΔP̃CORR = 0) is denoted as D, and represents the most simplistic theory that is widely used.15 The incorporation of hard core interactions (hij = −1 for rij < aij and cij = 0 for rij > aij) is denoted as HC. For D and HC, the electrostatic interactions only manifest themselves as an electroneutrality constraint that confines the polymer charges and their counterions to the gel phase, and as such these equations hold for all valencies q±. The full calculation for the electrostatic correlations, as described in this section, will be denoted by HC−HNC.



RESULTS AND DISCUSSION The aforementioned calculations, once solved for {ρi} and V′, demonstrate the essential features of a polymer gel as a function of the elasticity constant C, the valency q±, the bead radii ai, the charge fraction on the polymer f, and the salt concentration in the solution, ρS. In Figure 2 we demonstrate the effect of varying the last two parameters by plotting V′ versus ρS for a number of different values of f. We focus on changing the value of ρS, as this corresponds to a change in the chemical potential of the ionic species and is what is traditionally controlled in the 5056

dx.doi.org/10.1021/ma400372p | Macromolecules 2013, 46, 5053−5065

Macromolecules

Article

third that will manifest itself upon further investigation of the previous two; the enhanced deswelling in the HC−HNC as compared to the D situation, the reentrant swelling observed in HC−HNC, and a low-V′ first order phase transition from a swollen gel into a collapsed Wigner crystal state that will appear in multivalent systems.34,35 We consider each component in turn. Correlation-Enhanced Deswelling. In the low-ρS regime for any value of f, it is apparent in Figure 2 that the deswelling of the gel begins to occur more rapidly upon increasing ρS upon inclusion of ion−ion correlations. We investigate the origin of this by considering both the HC case and the HC−HNC case with a variety of valencies q±. The results of these different calculations are demonstrated for f = 0.10 in Figure 3, which

Figure 4. V′ versus ρS with f = 0.1, a+ = 0.8 nm, and for a monovalent q± = 1.0 salt solution. The gray line represents the D theory, which does not depend on the ion size a±. The colored points (shown in detail in the inset) correspond to the HC calculations, and the lines correspond to the theory derived from the results of the van der Waals equation of state for different values of the counterion radius a−. At ρS < 0.4, this theory gives excellent results. Deviations occur at higher densities due to higher order terms not included in the theory, but are included in the HC calculation. Deswelling behavior is increased slightly by an increase in the ion size.

where ρij and aij are the number of contacts and distance respectively between species i and j.58 This calculates the volume fraction of excluded volume in a dilute system, and as such in the correct limit it is possible to calculate the inputs to the thermodynamic equations as an alternative to the liquid state theory as58

Figure 3. V′ versus ρS with f = 0.1, a = 0.8 nm, and for a monovalent q± = 1.0 salt solution. The gray line represents the D theory, which is highly idealized. The blue curve represents the calculation with hardcore ions and without electrostatic correlations HC, and demonstrates both the previously observed enhanced deswelling and reentrant swelling. The inclusion of electrostatics (red curves) in the HC−HNC case only slightly strengthens these behaviors for the monovalent q± = 1.0 case, but strongly increases both deswelling and reentrant swelling effects in the divalent q± = 2.0 case.

ΔμCORR = ̃

and

bS bG − 1 − bS /2 1 − bG /2

(14)

π 6

(15)

58

̃ ΔPCORR =

∑ ρi ρj aij3 ij

These approximations are crude, but effective for describing the fundamental physics demonstrated in the more rigorous liquid state theory. This is demonstrated in Figure 4, which plots both the van der Waals and HC results in the regime before they diverge and matching is indeed observed in the low-ρS limit. The two thermodynamic connections have effects that counteract each other; the change in chemical potential Δμ̃ CORR is such that ions from solution are prevented from entering the gel due to the large decrease in translational entropy due to the excluded volume of the existing counterions, leading to an enhanced deswelling. The hard-core repulsions, however, also apply a larger pressure in dense systems that tends to drive the gel toward a more swollen state. In this regime the chemical potential term is more important, however the pressure term is non-negligible. Figure 5 demonstrates the incorporation of electrostatic interactions to the correlations calculated using HC−HNC, as well as the results of HC and D for comparison. HC−HNC calculations are performed with q± = 1.0, 2.0, and 3.0 to illustrate the effect of using multivalent ions or (alternatively) decreasing the background dielectric constant and strengthening the Coulombic interactions. Especially at larger valencies, there is a marked enhancement of gel collapse due to the presence of these interactions, and we can in analogous fashion to the HC case investigate existing theoretical predictions and their limitations when compared to many-body liquid state

plots V′ versus ρS for D, HC, and HC−HNC for both q± = 1.0 and q± = 2.0. The enhanced deswelling is different in all cases, with both the hard core interactions and the electrostatic interactions contributing to deviations from the D case in the same direction. We can interrogate each individually, by changing the ion radius in the HC case and the valency in the HC−HNC case and comparing them to existing theory. In Figure 4, we plot a number of different counterion sizes (keeping the polymer and co-ions the same size) as a function of the number concentration of ions (ρS) in the solution. There is a distinct difference between the D and HC situations, and likewise there is a small difference in the traces as the counterion sizes are adjusted. Both of these differences can be accounted for using existing theoretical descriptions of the hard-core interactions that serve to approximate the values for Δμ̃CORR and ΔP̃ CORR that we have calculated in a near-exact fashion using liquid state theory. The current conceptual description of a hard-sphere fluid at low densities corresponds to the results of the well-known van der Waals equation of state, which accounts for the two-body repulsions but no higher-order (e.g., three or more body) effects.58 We define an excluded volume parameter bα as 2π bα = ∑ ρ aij ,α 3 3 ij ij , α (13) 5057

dx.doi.org/10.1021/ma400372p | Macromolecules 2013, 46, 5053−5065

Macromolecules

Article

predicted correction to the pressure is carried out in the same way, with ΔP̃ CORR → ΔP̃CORR,HC + ΔP̃ CORR,ES, and the DH prediction for the pressure:59 ̃ ΔPCORR , ES =

κG 3 ⎡ 1 −2 ⎢1 + 2κGa − 24π ⎣ 1 + 2κGa

⎤ κ3⎡ 1 −2 ln(1 + 2κGa)⎥ − S ⎢1 + 2κSa − 1 + 2κSa ⎦ 24π ⎣ ⎤ ln(1 + 2κSa)⎥ ⎦ Figure 5. V′ versus ρS with f = 0.1, and a± = 0.8nm. Different valencies q± = 1.0,2.0, and 3.0 are shown, with points representing calculations using HC−HNC. Debye−Huckel theory (lines - based on eqs 16 and 17) is shown as well, demonstrating matching with the HC−HNC calculation at the low-ρS limit. HC calculation is also shown, but is largely hidden behind the monovalent points in this regime (black line). Monovalent ions deviate very little from D−H predictions, but divalent and trivalent ions quickly deviate due to correlation effects at all but very low ρS. q± = 3.0 ions also demonstrate a pronounced collapse transition at ρS ≈ 0.007 that is due to a first order transition to a Wigner crystal.

The results of these calculations are plotted in Figure 5 along with the calculations using HC−HNC. The monovalent q± = 1.0 case essentially follows the HC behavior, and the deviations from this behavior match almost exactly the DH theory. This suggests that the DH (and therefore the PB) approximation is largely accurate for the case of monovalent ions so long as hard core corrections are included. These effects are required if reentrant swelling is to be observed, however, since hard core interactions exhibit this feature only if correlations are calculated with full correlations using integral equation theory (as per Figure 3). The same cannot be said for q± = 2.0 and q± = 3.0, which match in the low-ρS limit but quickly move away from the DH calculations. Furthermore, nowhere in the regime tested do they match the HC results, demonstrating that the electrostatics are important at almost all ρS when multivalent ions are used. The significance of these plots is that the rapid departure of these rigorous integral equation calculations from DH theory demonstrates that ion−ion correlations are nonnegligible in multivalent or low-ε solvents where it is now preferable to use integral-equation theories to obtain the full set of many-body effects in a rigorous fashion. Figure 5 demonstrates another feature that is noteworthy; there is a discrete drop from a swollen state to a fully collapsed state for trivalent q± = 3.0 ions at ρS ≈ 0.007. This feature is physical, and we consider it to be a manifestation of approaching a “Wigner crystal” phase that can be interrogated in a facile manner through examination of both thermodynamic parameters and structural information gleaned from the local liquid-state theory.34,35

theory. To accomplish this, we separate the hard-core interaction out such that the errors inherent in the theoretical description of HC just presented have no bearing on the accuracy of the electrostatic correlations to ensure that the errors due to the following theory are presented with full clarity. Our calculation therefore maps Δμ̃ CORR → Δμ̃ CORR,HC + Δμ̃CORR,ES such that the first term Δμ̃ CORR,HC is calculated via integral equation theory like the HC case and the second term Δμ̃CORR,ES is given by theory. For the theoretical term, we use the ubiquitous Debye−Huckel (DH) linearized Poisson− Boltzmann approach, which predicts a correction to the chemical potential:59 ΔμCORR = ̃ , ES

κSq±2lB (1 + 2κSa)



κGq±2lB (1 + 2κGa)

(17)

(16)

where we have introduced the inverse screening length κi = (8πq±2ρilB)1/2 and assume a uniform bead radius a. The

Figure 6. V′ versus ρS with a± = 0.8 nm and q± = 3.0. Charge fraction f is slowly changed from f = 0.05 to f = 0.5, demonstrating both the increase in swelling due to the larger osmotic pressure of an increasing number of counterions in the gel as well as an increase in the ρS at which the gel fully collapses to a Wigner crystal. This occurs at the abrupt drop from a slightly swollen state V′ > 1.0 to the collapsed state at V′ = 1.0. Diagrams on the side demonstrate the conceptual picture of a transition from a dilute, gas-like ion phase to a highly correlated, liquid-like Wigner crystal. 5058

dx.doi.org/10.1021/ma400372p | Macromolecules 2013, 46, 5053−5065

Macromolecules

Article

Figure 7. (a) Red curve shows V′ versus ρS for a± = 0.8 nm, f = 0.05, and q± = 2.0 (HC−HNC). Gray line corresponds to D case, which contrasts greatly from the red curve. Full collapse is seen at ρS ≈ 0.004, and subsequently reentrant swelling is observed. Arrows correspond to the values of ρS plotted in part b, which plots the potential of mean force Ω̃G,MF as a function of V′. At ρS < 0.004 there are two minima corresponding to the swollen (V′ > 1.0) and collapsed (V′ = 1.0) states, however at ρS > 0.004 the swollen state becomes unstable and the system spontaneously evolves to V′ = 1.0. At large ρS, the minimum once again increases to V′ > 1.0 during reentrant swelling.

Figure 8. (a) Potentials of mean force Ω̃MF(V′) as the valency of the ion q± is adjusted, for f = 0.05 and ρS = 0.0015. As q± is increased, the minimum of Ω̃MF moves toward lower V′ as increased electrostatic interactions deswell the gel. Simultaneously, the Wigner crystal phase at V = 1.0 becomes more favorable until the swollen state becomes unstable at q± ≈ 2.2. This is captured in part b, which plots both the position of the minimum value of Ω̃MF as a function of q± that represent the swollen state as well as the local maximum value of Ω̃MF that represents the barrier between the swollen minimum and the minimum at V′ = 1.0. When these two meet at q± ≈ 2.2, the swollen state is spontaneously driven to collapse.

ρS there is a single, broad potential well located at the value indicated by the plot in Figure 7a. As ρS is increased, the location of this minimum in the traces shown in Figure 7b shifts quickly to the left toward V′ → 1, which is likewise apparent in Figure 7a. At ρS ≈ 0.004, however, the second minimum at ρS = 1.0 deepens well below the swollen minimum. This second minimum is not apparent in the trajectory shown in Figure 7a until this point, due to its higher energy and the minimization scheme used to determine the equilibrium state. At this point ρS ≈ 0.004, the swollen minimum at V′ > 1.0 becomes unstable and the state of the system is driven completely into the collapsed state at V′ = 1.0. This minimum is sharp rather than broad due to the nature of the V′ ≥ 1.0 constraint imposed on the calculation that represents the effect of the excluded volume of the polymer itself, which is not rendered in this particular calculation. Despite this collapse, in a process that will be demonstrated in the next section to be roughly independent of the presence of the collapse, a broad minimum once again appears and moves to the right with a further increase in ρS representing reentrant swelling behavior. The double potential well demonstrated in Figure 7b is characteristic of a classical first-order phase transition, and predicts a complete and abrupt collapse of the gel at a transition salt concentration ρS*. A finite-sized gel/solvent system would similarly demonstrate coexistence upon constraining the

Multivalent Ion-Driven Gel Collapse. In Figure 6, we elaborate on the observation in Figure 5 by plotting a number of values of f for trivalent q± = 3.0 systems. The discrete transition to a fully collapsed state is observed at low values of f, and occurs at roughly the same value of V′ for all conditions. In this situation, it is at a rather low value of V′ ≈ 2.0 however we expect the actual location of this transition to be increased upon the incorporation of dielectric or solvent effects. An identical transition occurs at q± = 2.0 as well, which is demonstrated in Figure 7a for f = 0.05, however small values of f are required to observe the collapse. We can examine this transition by calculating the energy landscape that gives rise to this phenomenon. This is accomplished by using the pressure calculation that is used to maintain mechanical equilibrium in the unconstrained gel phase, ∂Ω̃G/∂V, to derive a potential of mean force Ω̃G,MF(V′) at a given ρS as the value of V′ is controlled: Ω̃G , MF (V ′) = Ω̃(V ′ref ) +

∫V

⎞ ⎜ G ⎟ dV ′1 ′ref ⎝ ∂V ′ ⎠V 1 ′ V ′ ⎛ ∂Ω̃

(18)

Figure 7a indicates a number of points along the trajectory of the q± = 2.0, f = 0.05 trace of V′ versus ρS with colored arrows. It is at the values of ρS indicated by these arrows that the potentials of mean force Ω̃G,MF are plotted in Figure 7b. At low 5059

dx.doi.org/10.1021/ma400372p | Macromolecules 2013, 46, 5053−5065

Macromolecules

Article

Figure 9. hij(r) versus r for q± = 2.0 and f = 0.05, where h11 refers to the like-charge correlation function and h12 corresponds to opposite charges. Two concentrations, offset by 0.3 are shown; ρ1 = 0.004 is immediately preceding the collapse transition and ρ2 = 0.007 is in the collapsed state. Solution correlation functions are in darker colors, gel correlation functions in brighter colors. The inset highlights the oscillatory nature of hij,ρ2 in the gel, which demonstrates that the gel is driven toward a Wigner crystal state. This is conceptually demonstrated to the side for both the swollen (ρ1 = 0.004) state and the collapsed (ρ2 = 0.007) state (colors in the diagram correspond to similarly colored traces). In the swollen case, both the gel and solution correlations do not oscillate and the central ion (in this representation, positive) induces a local preference for negative ions. In the collapsed case, however, the gel phase oscillates which we indicate by the ordering of positive-minus-positive that is represented by the oscillatory hij,ρ2 function. It is this ordering that is characteristic of the collapse transition.

Figure 10. (a) hij(r) versus r for f = 0.05 and q± = 1.0. No oscillatory behavior is seen in the correlations, upon comparing like-charge correlation functions (gray, h11) to opposite-charge correlation functions (purple, h12), however hard-sphere correlations are observed. From darkest to lightest, ρ = 0.001, 0.020, and 0.04. (b) Same plot as in part a, only with q± = 2.0. Upon increasing valency, oscillatory h(r) is seen at large values of ρS.

volume V′ at a small but nonzero V′ near one of these transitions. We can further interrogate the electrostatic nature of this transition by manipulating the valency of the ions, and choose a value for ρS,n = 0.0015 to observe the evolution of the landscape Ω̃G,MF(V′). The results are shown in Figure 8a, which demonstrates similar features to the behavior of Ω̃G,MF with ρS. At low valencies, a broad minimum is observed that is representative of the swollen state. As q± is increased incrementally, this minimum moves toward lower values of V′ and ultimately a second minimum is observed at V′ = 1.0 to represent the collapsed state. At q± ≈ 2.2, the swollen minimum becomes unstable and the system is always driven to the collapsed state. The observation that there is a first order transition to the collapsed state due to either q± or ρS indicates that the evolution to this state relies simultaneously on the energy of the collapsed gel (via the charge valency) and the osmotic driving force to swell the gel (via ρS). This draws analogy to the calculations performed for polyelectrolyte solutions by Solis and Olvera de la Cruz that treat the relative energies between the collapsed and swollen chain conformations using two-state approximations where the polymer is treated as either a collapsed Wigner crystal or an extended rod.34,35 We no longer

require the division of conformational space into these extreme states, and treat a system that is topologically different (gel versus free chain) but the essential physics is similar; there is a competition between the ion entropy associated with the expanded state and the electrostatic energy associated with the collapsed state that cannot be captured with traditional Poisson−Boltzmann theories. Crucial to the theories of Solis and Olvera de la Cruz is the assumption of a “Wigner crystal” state that provides a description of the behavior of multivalent ions in dense solutions and a requisite for the understanding of the chemical potential due to ion correlations.34,35 The use of liquid state theories in HC−HNC allows the direct calculation of these types of correlations, and subsequently demonstrates the approach to a highly correlated Wigner crystal state in the collapsed gel. This can be visualized by plotting the correlation function hij as a function of r. Figure 7 demonstrates a characteristic plot of V′ versus ρS for q± = 2.0 and f = 0.05, which demonstrates a collapse transition. From along this trajectory, we plot in Figure 9 the hij associated with the collapse transition (both for the solution and for the gel). Following the discrete transition observed in the thermodynamic variables, there is a corresponding discrete change in the 5060

dx.doi.org/10.1021/ma400372p | Macromolecules 2013, 46, 5053−5065

Macromolecules

Article

Figure 11. (a) V′ versus ρS for f = 0.05 and a+ = 0.8nm, in the HC calculation that neglects electrostatic interactions. a− is varied from 0.3−1.2nm, which is demonstrated to strongly effect the value of ρS where reentrant swelling occurs. (b) V′ versus VexρS − CHC [Δμ̃ IG − Δμ̃0.8,IG] (CHC = 0.25), based on eq 19, which describes the underlying physics governing reentrant swelling. Δμ̃ 0.8,IG serves as a reference such that the case of a− = 0.8 is not shifted by the ideal gas chemical potential term. All traces fall onto a universal curve upon rescaling via eq 19. This rescaling is conceptually illustrated in the bottom graphs whose colored frames correspond to the similarly colored arrows in part b. At low salt, the osmotic pressure inside the globule leads to large swelling due to the low osmotic pressure outside the globule. At medium salt, the number of ions outside are similar to the number of mobile (i.e., the nonorange ions on the gel diagram) and there is no driving force to swell. At high salt concentrations, the number of ions outside exceeds the mobile ions inside the gel and there is a chemical potential difference that drives these ions into the gel such that it undergoes reentrant swelling.

local structure. A distinct increase in the size of the oscillatory portion of the correlation function is associated with the formation of structure on larger length scales than the neighboring ions such that there is a clear excess of similarly charged ions two ion-diameters away from the central ion. These bridging correlations are characteristic of a transition toward a Wigner crystal phase that has crystalline-like order on small length scales, though full Wigner crystal behavior is prevented due to presence of the neutral polymer enforcing the V′ ≥ 1.0 constraint. This is observed also as the solution gets more dense for multivalent ions, which is observed regardless of the state of the gel phase. This is shown in Figure 10a, which plots hij as a function of r for a number of solution concentrations ρS. These clusters in solution manifest once again as oscillatory correlation functions, which become significant and result in a large cohesive energy at large ρS. This oscillation does not occur significantly in monovalent solutions, as demonstrated in Figure 9b, which plots the analogous correlation functions h for q± = 1.0. While we have restricted ourselves to symmetric valency salts, it is expected that asymmetric valencies and sizes would enrich this behavior further. We note that these features are indeed experimentally observed in polyelectrolyte gels, namely in the work by Horkay et al.28,29,32,43 While these investigations typically involve ion exchange (monovalent counterions are replaced with di- or trivalent counterions) that is not something investigated in this work, the physics of the observed discrete transitions can nevertheless be attributed to similar first-order correlationbased effects; the Coulombic attractions between charged species are being significantly increased upon introduction of highly valent counterions. These works attribute the discrete

and sudden collapse to concentration-dependent contributions to the χ-parameter, which is theoretically known to induce such chain collapse.18 Such theoretical models corroborate this statement, however leave open the question of the molecular interpretation of the phenomenological form used for the χparameter. Subsequent simulation also reproduced these features, and postulated the possibility of electrostatic effects leading to collapse.33 In our work we have assumed that the dispersive interactions implied by a χ-parameter are not rendered, and demonstrate that this collapse can indeed be observed solely through electrostatic correlations. While there are similarities between these two collapse transitions, namely the first order character of the collapse, there are some quantitative differences upon comparison of our calculations in Figure 5 to the work of Horkay et al.28,29,43 with regards to the value of V′ at which the transition is observed. This value of V′ is much larger in experiments. We note that some of the effects neglected in our work, namely correlations due to connectivity and dielectric self-energy effects, may lead to this disparity. This suggests further work into these aspects of the calculation will be needed to obtain complete quantitative reproduction of experimental data. Reentrant Swelling in Polyelectrolyte Gels. At large values of ρS for HC and HC−HNC the polymer gels begin to swell in a process known as reentrant swelling. This transition is observed in the literature,13 and recent work on traditional D models have shown that it is recovered upon addition of a phenomenological excluded volume term to the expression for the gel free energy.44 We no longer require such a phenomenological excluded volume term, though we still retain such a term for the polymer, due to the inclusion of excluded volume in both HC and HNC in a more rigorous 5061

dx.doi.org/10.1021/ma400372p | Macromolecules 2013, 46, 5053−5065

Macromolecules

Article

Figure 12. (a) V′ versus VexρS for f = 0.05 and a± = 0.8 nm. The HC−HNC case is used, with q± varied from 0.0−3.0. Reentrant swelling is shifted to the right upon increasing the electrostatic interactions through the increase in q±. (b) Using the prediction in eq 21, we plot V′ versus VexρS − CESq±2 for the same curves (CES ≈ 0.009). All of the traces lie along the same universal reentrant swelling curve, suggesting that eq 21 captures the essential physics to describe the effect of valency on reentrant swelling. This effect is shown conceptually in the bottom figures, whose colors correspond to the limiting colors in parts a and b. At low valencies, Coulombic attractions are small, but at high valencies the Coulombic attractions (symbolically represented with purple arrows) are large and counteract the driving force for these ions to move into the globule. A much larger difference in mobile ion concentrations must be realized before this attraction can be overcome and reentrant swelling is realized.

smaller ions. Normalizing for these effects, we aim to collapse all of the traces onto a single universal reentrant swelling curve. The argument for such universal behavior is that the chemical potential in a gel that is constrained to be slightly swollen (V′ ≈ 1−2) will undergo a transition from Δμ̃ < 0 to Δμ̃ > 0 as ρS is increased, and this represents the universal reentrant swelling point. The chemical potentials of each phase can be written as the two contributions, an ideal portion μIG ∼ ln (ρ+ρ−) and a hard core portion μHC ∼ ∑iρiai3 ≈ Vex(ρ− + ρ+). This results in the relationship for the critical reswelling concentration ρS*:

fashion. This allows us to probe more deeply the physics concerning the reentrant swelling and the effect that charge plays in manipulating this behavior. Reentrant swelling only appears upon the introduction of excluded volume effects, so we initially examine only the HC case to avoid effects due to electrostatics. It is possible to change the magnitude of these effects by adjusting the radius of the ions. We demonstrate this in Figure 11a, where the counterion size a− was changed from 0.3 to 1.2 nm. The location of the reentrant swelling changes drastically, with larger ions rapidly inducing reentrant swelling while small ions demonstrate this effect only at very large values of ρS and exhibit much slower reswelling as ρS increases. This suggests that the system would demonstrate the appropriate limit of ρS,reswelling → ∞ as ai → 0. Conceptually, the inclusion of a V′ > 1 constraint on the volume of the gel combined with the finite volume of the ions restricts the ρ+ and ρ− from continuously increasing as quickly as ρS; while in the D case the overall ρ in the gel is always slightly larger than outside, the HC case allows for 2ρS > ρ− + ρ+. This ultimately results in an increase in the external chemical potential in the solution relative to a gel that drives ions into the gel phase such that the gel must swell to accommodate the extra ions. The argument that the reentrant swelling is driven by excluded volume suggests that it is necessary to renormalize the concentration to be a volume percent ρS → ρSVex where Vex ∼ 4π(a−3 + a+3)/3 represents the excluded volume. The difference in the ideal gas portion of the chemical potential must be likewise taken to account, as the number difference at comparable ρSVex drives the ionic species into the gel; the ideal portion of the osmotic pressure is much smaller if there are a small number of large ions than if there were large numbers of

[ΔμIG + ΔμHC ]P = Δμ* = 0 ⎡ ⎛ ρ ρ ⎞⎤ ⎢V (ρ* − ρ ) ∼ ln⎜ −, G +, G ⎟⎥ G ⎜ *2 ⎟⎥ ⎢ ex S ⎝ ρS ⎠⎦ ⎣ P

(19)

where the subscript P indicates that this should be evaluated at the value of V′ of the plateau, and the value Δμ* is the change in chemical potential at the transition point due to excluded volume and ideal gas contributions. Δμ* = 0 for the HC case, but this will need to be revised upon including electrostatics. Renormalization based on this scaling argument is shown in Figure 11b using a proportionality constant CHC = 0.25 for all traces, which indeed demonstrates the collapse of the reentrant swelling using eq 19. We note that these arguments qualitatively describe the trends seen in Figure 2 as well; the larger values of f correspond to a larger concentration of ions within the gel, leading to a larger concentration of ions outside the gel needed to overcome the gel chemical potential. This in turn leads reentrant swelling at larger values of ρS for larger f, which is indeed observed. 5062

dx.doi.org/10.1021/ma400372p | Macromolecules 2013, 46, 5053−5065

Macromolecules

Article

experimental results on similar systems. The literature on polymer gels, both biological and synthetic, qualitatively demonstrates many of the effects described in this article. In biological gels, chromatin systems are known to exist in a slightly swollen state at physiological state that can swell upon both increasing and decreasing the monovalent salt concentration.3 This describes the reentrant and osmotic swelling regimes, respectively. Furthermore, the addition of divalent salts is known to hypercondense (completely deswell) chromatin in a fashion similar to the collapse seen in Figure 3.3 Experimental observation of both this same collapse transition and reentrant swelling is observed in synthetic gels, as well.13,28,29,43 The phase behavior of nongelated polyelectrolytes in solution provide classical results that allow for such comparisons; in particular, the observations of full collapse and reentrant swelling are observed in a great number of systems, leading to the longstanding observation that such effects are universal in polyelectrolytes.54−57,60,61 A universal and ubiquitous phase diagram describes a large number of systems by plotting ρS versus polymer concentration ϕ for both synthetic and biological polyelectrolytes.35,54 At low values of ρS, there is a swollen region that expands with polymer concentration ϕ and is characterized by a maximization of counterion entropy through solution formation.35,54 Solution formation in this case is analogous to the swelling of a polymer gel to maximize counterion entropy.35,54 Swelling and/or dissolution grow stronger with an increase in the number of counterions in the polymer phase, which corresponds to an increase in ϕ for the polyelectrolyte solution and f for the polyelectrolyte gel. The reentrant swelling is largely independent of ϕ for the polyelectrolyte solutions, and we would expect the same independence for f for a polyelectrolyte gel when the gel reaches the collapsed state due to the hard-core-driven nature of this transition.35,54 We can demonstrate that the same type of graph can indeed be reproduced for a gel, and plot ρS versus f in Figure 13. The general features are essentially the same as in polyelectrolyte solutions, namely a f-dependent initial collapse transition and a f-independent reentrant swelling transition. We emphasize that many of the comparisons to experimental work that we just described involve somewhat different systems; practical differences involving ion exchange, ion binding, and dispersive or hydrophobicity effects are not incorporated into our model that is essentially the so-called “primitive model” of polyelectrolyte gels. We do not presume to capture every physical nuance of these systems, nevertheless it is indeed significant that qualitative agreement occurs solely upon introduction of ion correlation effects in a rigorous fashion, as it suggests that successful theories of these more complicated systems will require a similar level of description of ion correlations. We can provide a more direct comparison to systems that are described in simulation work on the same “primitive model”. In that spirit, we plot in Figure 14 the swelling V′ of our model versus the renormalized Bjerrum length (Γ = lB/lB,0, where lB,0 is the Bjerrum length of water, is changed by decreasing the dielectric constant ε in the expression for uij in eq 10 for both solution and gel phases). This is plotted for the two limiting cases seen in Yin et al.32 (q± = 1.0 and q+ = 1.0/q− = 2.0, f = 1.0), which is simulation work that likewise considers the “primitive model” for polyelectrolyte gels. For convergence of the integral equations, we use ρS = 0.001 to approximate the salt-free situation in the simulations in

This approach can be extended to include the effect of electrostatics through a similar scaling approach. Figure 12a shows the effect of turning on the electrostatic interactions in the HC−HNC scheme. As the valency q± is increased, there is a large impact on the location of the reentrant swelling away from the original point where ΔμIG + ΔμHC = Δμ* = 0. The value Δμ* was defined to facilitate the introduction of electrostatic effects, since we can expand around the transition solution concentration ρ*S Δμ*(ρS ) = Δμ*(ρS*) + K1(ρS − ρS*) + + ···

K2 (ρ − ρS*)2 2 S (20)

where the first term Δμ*(ρ*S ) = 0 due to the previous definition of the variable, and if we remain close to ρ*S we can truncate this expansion to the first term. Ki are phenomenological coefficients whose values will be unnecessary to know in the scaling analysis. With the inclusion of the electrostatics in the HC−HNC case, we require another term in the balance of the chemical potential to capture the difference due to the electrostatics. We write this chemical potential as μES ∼ Mq2± /a, where M is the Madelung constant for a dense phase that is dependent on the structure of the ion ordering but typically of order unity.34 The structure is different between the gel and the solution phase, meaning that there is a difference ΔM between the phases. We thus write the criteria for the reentrant swelling concentration ρ*S : Δμ*(ρS* − ρS*,0 ) + ΔμES ⎡ ΔM ⎤ 2 q ≈ K1(ρS* − ρS*,0 ) + ⎢ ⎣ a ⎥⎦ ± =0 ΔM 2 CES 2 q = q ρS* − ρS*,0 = K1a ± Vex ±

(21)

where ρS,0 * is the transition solution concentration for the HC case and the bracketed term in the second line is representative of the different structure between the solution and gel phases, and is a relatively small value and assumed to not be a function of q±. CES is a proportionality constant that carries structural information through the incorporation of the ΔM between the two phases. We shift the curves in Figure 12b such that they all fall along a universal reentrant swelling curve based on the predicted rescaling of ρS → VexρS − CESq±2. These results demonstrate the underlying physics governing the reentrant swelling behavior, namely that the hard-core repulsion of the ions constrains the volume such that the concentration of ions in the slightly swollen or collapsed gel can be surpassed by the solution concentration. This leads to a reversal in the chemical potential that contains contributions due to hard-core repulsions, translational entropy, and electrostatic attractions. Importantly, these reversals are largely independent of the chain ionization f and the gel elasticity C, which will be of interest in our qualitative comparisons to the experimental literature in the next section.



QUALITATIVE COMPARISON TO EXPERIMENTAL AND SIMULATION LITERATURE While we leave a quantitative comparison with experimental literature for future work, we can compare our results to existing experiments on gels and draw analogies to well-known 5063

dx.doi.org/10.1021/ma400372p | Macromolecules 2013, 46, 5053−5065

Macromolecules

Article

calculation of ion correlations. Despite these differences, the collapse of the gel at the appropriate values of Γ importantly establishes that we capture most of the essential elements of that work without requiring an “effective χ-parameter”, and establishes charge correlations as of primary importance in the understanding of polyelectrolyte gels.



CONCLUSIONS We have incorporated integral equation liquid state theory into a macroscopic thermodynamic description of a polyelectrolyte gel using a self-consistent scheme that allows the volume and concentration of the species within the gel to evolve to the equilibrium gel structure. This permits the realization of both system observables and molecular-level details in a rigorous fashion that accounts for many-body interactions that are important for the description of ion−ion correlations in a nearcomplete fashion. The primary assumption made in our work consist of the traditional HNC closure approximation, which involves neglecting what are known as “bridge diagrams” that are almost negligible in the thermodynamic limits described here. There are also assumptions made in neglecting connectivity effects on ion−ion correlations in the gel phase, dielectric effects, and dispersive forces, as well as in the simplistic model for the gel elasticity. Many of these effects, now that this model has been clarified in their absence, can be included into this formalism in a facile manner to provide an increasingly accurate description of a polymer gel. Nevertheless, most of the experimentally observed effects are readily apparent in this model; we have demonstrated an enhanced deswelling due to the addition of salt that can readily undergo a first-order collapse transition due to multivalent ions upon inclusion of correlation effects. This collapse has been observed in numerous simulation and experimental studies on gels28,29,43 and polyelectrolyte solutions,3,56 and we qualitatively match the trends observed in these systems. Furthermore, we observe a reentrant swelling transition that likewise matches the qualitative behavior of experimental results of polyelectrolyte gels and solutions,13,56,60 and can describe reentrant swelling in terms of a chemical potential reversal such that reentrant swelling curves for different ion sizes and valencies can be collapsed onto a single universal curve. This work represents a theoretical description of the response of a strongly charged polyelectrolyte gel due to the chemical control of an external salt solution. We expect that this approach, which can also be applied to asymmetric electrolyte solutions, could be extended in the near future to include such effects as varying dielectric constants and dispersive interactions involving either ions or gel thermodynamics, and weakly charged polymers that respond to pH. Furthermore, such theories may enable the quantitative matching of theory to experimental systems of interest that are both biological and synthetic in nature.

Figure 13. Log−log phase diagram of ρS versus f, which is analogous to the ρS versus ϕ (polymer concentration) plots ubiquitous in experiments of polyelectrolyte solutions. The collapsed and swollen regimes are indicated, corresponding to soluble and phase separated regimes in experimental literature. The red points correspond to the reentrant swelling transition, and the black points correspond to the collapse transition. Lines are not fits, and simply meant to guide the eye. Schematics of the various phases are shown; the collapsed state is a nearly dry gel with correlated ions, the low-salt swollen state is an expanded gel due to the osmotic pressure of the fixed ions and their counterions, and the high-salt reentrantly swollen state is expanded to incorporate the salt ions that are driven into the gel from the surrounding solution.

Figure 14. Swelling behavior (V′ versus Γ) for two gels with f = 1.0. The black curve corresponds to q+ = 1.0 and q− = 2.0, while the red curve corresponds to the symmetric case q± = 1.0. The values of the transition Γ for both cases correspond to the limiting values in Figure 3 in Yin, et al.,32 demonstrating that our model reproduces the qualitative behaviors observed in simulations that include connectivity constraints and explicit polymer elasticity. We note that differences in gel elasticity models between the two systems results in different limiting values as Γ → 0, and differences in the apparent sharpness of the transition may be due to the lack of connectivity in the ion correlation calculation.



Yin et al.32 Besides this approximation and the different model for the chain elasticity, which is expected to decrease the overall swelling in the limit of low Γ, we observe qualitative matching between the two situations as Γ is changed. Important features are reproduced, such as chain collapse for q± = 1.0 around Γ ≈ 7−8 and a transition for q+ = q−/2.0 = 1.0 around Γ ≈ 2. There is a correspondence, roughly, in the values of Γ at which these transitions occur. There are some apparent differences between Figure 14 and Yin et al.32 in the sharpness of the transitions, which may be due to the absence of connectivity in the

AUTHOR INFORMATION

Corresponding Author

*E-mail: (M.O.d.l.C) [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors acknowledge support from NSF Grant Number DMR-0907781. C.E.S. thanks the Northwestern International 5064

dx.doi.org/10.1021/ma400372p | Macromolecules 2013, 46, 5053−5065

Macromolecules

Article

(34) Solis, F. J.; Olvera de la Cruz, M. J. Chem. Phys. 2000, 112, 2030−2035. (35) Solis, F. J.; Olvera de la Cruz, M. Eur. Phys. J. E. 2001, 4, 143− 152. (36) Dobrynin, A. V.; Deshkovski, A.; Rubinstein, M. Phys. Rev. Lett. 2000, 84, 3101−3104. (37) Naji, A.; Netz, R. R. Eur. Phys. J. E 2004, 13, 43−59. (38) Naji, A.; Jungblut, S.; Moreira, A. G.; Netz, R. R. Physica A 2005, 352, 131−170. (39) Rouzina, I.; Bloomfield, V. A. J. Phys. Chem. 1996, 100, 9977− 9989. (40) Solis, F. J. J. Chem. Phys. 2002, 117, 9009−9015. (41) Hsiao, P. Y.; Luijten, E. Phys. Rev. Lett. 2006, 97, 148301. (42) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids, 3rd ed.; Elsevier: Boston, MA, 2006. (43) Horkay, F.; Basser, P. J.; Hecht, A. M.; Geissler, E. Macromol. Biosci. 2002, 2, 207−213. (44) Jha, P. K.; Zwanikken, J. W.; Olvera de la Cruz, M. Soft Matter 2012, 8, 9519. (45) Nightingale, E. R. J. Phys. Chem. 1959, 63, 1381−1387. (46) Quesada-Pérez, M.; Ibarra-Armenta, J. G.; Martín-Molina, A. J. Chem. Phys. 2011, 135, 094109. (47) Capriles-González, Sierra-Martín, B.; Fernández-Nieves, A.; Fernández-Barbero, A. J. Phys. Chem. B 2008, 112, 12195−12200. (48) Jha, P. K.; Zwanikken, J. W.; Detcheverry, F. A.; de Pablo, J. J.; Olvera de la Cruz, M. Soft Matter 2011, 7, 5965. (49) Wang, Z. G. Phys. Rev. E 2010, 81, 021501. (50) Quesada-Pérez, M.; Martín-Molina, A.; Hidalgo-Á lvarez, R. J. Chem. Phys. 2004, 121, 8618. (51) Dobrynin, A. V.; Rubinstein, M. Prog. Polym. Sci. 2005, 30, 1049−1118. (52) Yethiraj, A. J. Chem. Phys. 1998, 108, 1184−1192. (53) We emphasize that care must be taken in including and interpreting a χ-parameter within the context of a charged system. It is practical and commonplace to use χ to represent all the interactions in the system, including charges, however, such a practice lacks a sound theoretical basis. It is well-known in the theory literature that χparameters cannot be derived from charged systems due to the divergence of the spatial integral of 1/r. Therefore, χ-parameters should only be strictly used to represent short-range interactions. (54) Olvera de la Cruz, M.; Belloni, L.; Delsanti, M.; Dalbiez, J. P.; Spalla, O.; Drifford, M. J. Chem. Phys. 1995, 103, 5781−5791. (55) Raspaud, E.; Olvera de la Cruz, M.; Sikorav, J. L.; Livolant, F. Biophys. J. 1998, 74, 381−393. (56) Raspaud, E.; Chaperon, I.; Leforestier, A.; Livolant, F. Biophys. J. 1999, 77, 1547−1555. (57) Jia, P.; Zhao, J. J. Chem. Phys. 2009, 131, 231103. (58) Kardar, M. Statistical Physics of Particles; Cambridge University Press: Cambridge, U.K., 2007. (59) McQuarrie, D. A. Statistical Mechanics; University Science Books: Sausalito, CA, 2000. (60) Zhang, F.; et al. Proteins: Structure 2010, 78, 3450−3457. (61) Widom, J.; Baldwin, R. L. J. Mol. Biol. 1980, 144, 431−453.

Institute for Nanotechnology for an International Institute for Nanotechnology Postdoctoral Fellowship. The computational cluster and J.W.Z. are funded by the Office of the Director of Defense Research and Engineering (DDR&E) and the Air Force Office of Scientific Research (AFOSR) under Award no. FA9550-10-1-0167.



REFERENCES

(1) Tasaki, I.; Byrne, P. M. Biopolymers 2004, 32, 1019−1023. (2) Tasaki, I. Jpn. J. Physiol. 1999, 49, 125−138. (3) Poirier, M. G.; Monhait, T.; Marko, J. F. J. Cell. Biochem. 2002, 85, 422−434. (4) Lieleg, O.; Ribbeck, K. Trends Cell Biol. 2011, 21, 543−551. (5) Button, B.; et al. Science 2012, 337, 937−941. (6) Marudova, M.; MacDougall, A. J.; Ring, S. G. Carb. Res. 2004, 339, 1933−1939. (7) Fels, J.; Orlov, S. N.; Grygorczyk, R. Biophys. J. 2009, 96, 4276− 4285. (8) Guaqueta, C.; Sanders, L. K.; Wong, G. C. L.; Luijten, E. Biohpys. J. 2006, 90, 4630−4638. (9) Kang, Y.; Walish, J. J.; Gorishnyy, T.; Thomas, E. L. Nat. Mater. 2007, 6, 957−960. (10) Walish, J. J.; Fan, Y.; Centrone, A.; Thomas, E. L. Macromol. Rapid Commun. 2012, 33, 1504−1509. (11) Richtering, W.; Pich, A. Soft Matter 2012, 8, 11423−11430. (12) Sidorenko, A.; Krupenkin, T.; Taylor, A.; Fratzl, P.; Aizenberg, J. Science 2007, 315, 487. (13) Henderson, K. J.; Zhou, T. C.; Otim, K. J.; Shull, K. R. Macromolecules 2010, 43, 6193−6201. (14) Schmidt, J.; Cha, C.; Kong, H. Soft Matter 2010, 6, 3930−3938. (15) Flory, P. J. Principles of Polymer Chemistry; Cornell University Press: Ithaca, NY, 1953. (16) Tanaka, T.; Fillmore, D.; Sun, S. T.; Nishio, I.; Swislow, G.; Shah, A. Phys. Rev. Lett. 1980, 45, 1636−1639. (17) Quesada-Perez, M.; Maroto-Centeno, J. A.; Forcada, J.; HigaldoAlvarez, R. Soft Matter 2011, 7, 10536. (18) Onuki, A. Phase Transition Dynamics; Cambridge University Press: Cambridge, U.K.,2002. (19) Jha, P. K.; Solis, F. J.; de Pablo, J. J.; Olvera de la Cruz, M. Macromolecules 2009, 42, 6284−6289. (20) Matsuo, E.; Tanaka, T. Nature 1992, 358, 482−485. (21) Dervaux, J.; Amar, M. B. Annu. Rev. Condens. Matter Phys. 2012, 3, 311−332. (22) Hua, J.; Mitra, M. K.; Muthukumar, M. J. Chem. Phys. 2012, 136, 134901. (23) Longo, G. S.; Olvera de la Cruz, M.; Szleifer, I. Macromolecules 2011, 44, 147−158. (24) Jha, P. K.; Zwanikken, J. W.; de Pablo, J. J.; Olvera de la Cruz, M. Curr. Opin. Solid State Matter 2011, 15, 271−276. (25) Zwanikken, J. W.; Jha, P. K.; Olvera de la Cruz, M. J. Chem. Phys. 2011, 135, 064106. (26) Zwanikken, J. W.; Olvera de la Cruz, M. Proc. Natl. Acad. Sci. U.S.A. 2013, 110, 5301−5308. (27) Zwanikken, J. W. Moderately Coupled Charged Fluids near Dielectric Interfaces and in Confinement. In Electrostatics of Soft and Disordered Matter; Dean, D. S.; Dobnikar, J.; Naji, A.; Podgornik, R., Eds.; Pan Stanford: Singapore, 2013. (28) Horkay, F.; Tasaki, I.; Basser, P. J. Biomacromolecules 2000, 1, 84−90. (29) Horkay, F.; Tasaki, I.; Basser, P. J. Biomacromolecules 2001, 2, 195−199. (30) Schneider, S.; Linse, P. Macromolecules 2004, 37, 3850−3856. (31) Yin, D. W.; Yan, Q.; de Pablo, J. J. J. Chem. Phys. 2005, 123, 174909. (32) Yin, D. W.; Horkay, F.; Douglas, J. F.; de Pablo, J. J. J. Chem. Phys. 2008, 129, 154902. (33) Yin, D. W.; Olvera de la Cruz, M.; de Pablo, J. J. J. Chem. Phys. 2009, 131, 194907. 5065

dx.doi.org/10.1021/ma400372p | Macromolecules 2013, 46, 5053−5065