Modeling deswelling, thermodynamics, structure, and dynamics in ionic microgel suspensions
MModeling deswelling, thermodynamics, structure, and dynamics in ionic microgel suspensions
Mariano E. Brito, Alan R. Denton, and Gerhard Nägele Institute of Complex Systems, ICS-3, Forschungszentrum Jülich GmbH, 52425 Jülich, Germany Department of Physics, North Dakota State University, Fargo, ND 58108-6050 USA (Dated: October 3, 2019)Ionic microgel particles in a good solvent swell to an equilibrium size determined by a balance of electrostaticand elastic forces. When crowded, ionic microgels deswell owing to a redistribution of microions inside andoutside the particles. The concentration-dependent deswelling affects the interactions between the microgels,and consequently the suspension properties. We present a comprehensive theoretical study of crowding effectson thermodynamic, structural, and dynamic properties of weakly cross-linked ionic microgels in a good solvent.The microgels are modeled as microion- and solvent-permeable colloidal spheres with fixed charge uniformlydistributed over the polymer gel backbone, whose elastic and solvent-interaction free energies are describedusing Flory-Rehner theory. Two mean-field methods for calculating the crowding-dependent microgel radius areinvestigated, and combined with calculations of the net microgel charge characterizing the electrostatic part of aneffective microgel pair potential, with charge renormalization accounted for. Using this effective pair potential,thermodynamic and static suspension properties are calculated including the osmotic pressure and microgel pairdistribution function. The latter is used in our calculations of dynamic suspension properties, where we accountfor hydrodynamic interactions. Results for diffusion and rheological properties are presented over ranges ofmicrogel concentration and charge. We show that deswelling mildly enhances self- and collective diffusionand the osmotic pressure, lowers the suspension viscosity, and significantly shifts the suspension crystallizationpoint to higher concentrations. The paper presents a bottom-up approach to efficiently computing suspensionproperties of crowded ionic microgels using single-particle characteristics.
PACS numbers:
I. INTRODUCTION
Mesoscopically sized microgel particles are of both funda-mental and technological interest owing to their strong sensi-tivity to external parameters, including temperature (solventquality), pH, ionic strength, and concentration, which con-trol their equilibrium size and soft interactions [1, 2]. Be-cause of their profound environmental adaptability, and capac-ity to partially interpenetrate and host small molecular species,microgels have a multitude of possible applications, e.g.,as drug-delivery vehicles, functionalized colloids, switchablemembrane filters, and tunable microreactors [3].An important subgroup is formed by ionic microgels, whichare typically globular particles consisting of cross-linkedpolyelectrolyte chains. When dispersed in a polar solvent un-der good solvent conditions, an ionic microgel particle be-comes charged due to the dissociation of counterions fromionizable groups on the polymer backbone. For weak cross-linking, the soft particle swells to an equilibrium size that canbe substantially larger than that in the dry state. It has beenshown experimentally [4–6] and theoretically [5–11] that theparticles deswell with increasing concentration, which is mostpronounced for low background ionic strength and smallerconcentrations well below the threshold where the microgelsstart to overlap. This behavior distinguishes ionic microgelsfrom nonionic ones, since for the latter deswelling (or inter-penetration and facetting) is observed only at very high con-centrations [12, 13].In theory-simulation studies by Denton et al. [9, 10],where the interplay of elastic and electrostatic influencesis accounted for using a coarse-grained model, microgeldeswelling with increasing concentration is explained by a re- distribution of counterions. These ions increasingly permeatethe microgel, to maintain balance of the electro-elastic pres-sure inside and outside the particles. The shrinkage of micro-gels with increasing concentration is thus accompanied by adecreasing net microgel charge and narrowing of the particlesize distribution around the equilibrium mean radius. At lowbackground ionic strength (i.e., low salt content), the meancounterion concentration outside a microgel is much smallerthan inside, with an accordingly strong sensitivity of the meanmicrogel size to concentration changes. With increasing saltcontent, the strong inside-outside counterion concentrationgradient, and the concomitant ion pressure gradient, are flat-tened out, reducing the sensitivity of the microgel size to con-centration variations.Owing to the counterion-induced deswelling, sufficientlysoft ionic microgels can penetrate apertures considerably nar-rower than their dilute-concentration size, at concentrationsbelow particle overlap [14]. This ability has potentially im-portant applications for drug delivery, microfluidics and fil-tration. For example in pressure-driven membrane filtrationused to concentrate and purify microgels, deswelling can leadto an unwarranted enhanced clogging of membrane pores. Onthe other hand, deswelling reduces the formation of a fluidconcentration-polarization layer and a solid filter cake layer ofparticles accumulated on the membrane surface, which causeadditional (effective) hydraulic resistance. The formation ofthese inhomogeneous layers is determined by concentration-dependent transport properties of crowded microgel suspen-sions in conjunction with osmotic pressure effects, namely bythe collective diffusion coefficient and the suspension viscos-ity [15–17].The filtration example illustrates the demand for studyingdiffusion and rheological transport properties of soft ionic mi- a r X i v : . [ c ond - m a t . s o f t ] O c t crogel suspensions in general, and the effects of counterion-regulated deswelling in particular. In this paper, we present acomprehensive theoretical exploration of dynamic and equi-librium microstructural properties of fluid-phase suspensionsof ionic microgels in the swollen state. Being of interestin their own right, microstructural properties such as the ra-dial distribution function (rdf) and static structure factor arealso required as input in the calculation of dynamic sus-pension properties, including generalized sedimentation andcollective diffusion coefficients and the high-frequency andzero-frequency viscosities. Following earlier work by Den-ton et al. [9, 10, 18, 19], we model the ionic microgelsin a coarse-grained way as microion- and solvent-permeablemonodisperse elastic colloidal spheres, with the charged sitesof the cross-linked polymer gel backbone described by a uni-form charge distribution. This description is reasonable, un-der the proviso that the cross-linker density does not varystrongly along the particle radius. We describe the elastic andsolvent-interaction free energy contributions of a microgel us-ing Flory-Rehner theory [20–22] for uniform cross-linker dis-tribution. For calculating the electrostatic semi-grand free en-ergy contribution of microgels in a concentrated suspension,in Donnan equilibrium with a 1:1 strong electrolyte reser-voir, we use two different mean-field methods, namely, thespherical Poisson-Boltzmann cell model (PBCM) approach ofDenton and Tang [9] and a first-order thermodynamic per-turbation theory (TPT) method of Weyer and Denton [10],based on a multi-center linear-response approach. The equi-librium microgel radius is obtained from minimizing the to-tal suspension free energy, equivalent to enforcing the bal-ance of total pressure inside and outside a particle. In com-bination with an effective electrostatic pair potential expres-sion for ionic microgels derived from the multi-center ap-proach [18, 23, 24], we determine the pressure and osmoticcompressibility of microgel suspensions as well as the mi-crogel pair distribution function and static structure factorusing the hypernetted chain (HNC) and thermodynamicallyself-consistent Rogers-Young (RY) integral-equation methods[25]. For non-overlapping particles, the effective pair po-tential is of a screened-Coulomb form, akin to the potentialfor ion-impermeable charge-stabilized colloidal particles, butwith a coupling strength that decreases with increasing con-centration and ionic strength of the suspension. The net mi-crogel charge likewise decreases with increasing concentra-tion. For overlapping microgels, the effective electrostaticpotential remains finite, and is augmented in our model bya soft Hertz potential accounting for elastic repulsion at mod-est overlap [26, 27]. Van der Waals attraction between theweakly cross-linked microgels can be neglected due to theirhigh solvent content.The microgel pair distribution function is used as input toour calculations of dynamic suspension properties. Semi-analytic methods are used to calculate dynamic properties,whose good performance has been established, by compari-son with elaborate dynamic computer simulations, for a vari-ety of colloidal model systems describing globular proteins,impermeable charge-stabilized colloids, and non-ionic spher-ical microgels. These methods account for the salient hydro- dynamic particle interactions mediated by intervening solventflow. In our assessment of deswelling effects, the results ob-tained for various static and dynamic suspension propertiesare compared with the ones for a (fictitious) reference suspen-sion of constant-sized microgels.The paper is structured as follows: In Sec. II, we reviewthe derivation of the effective one-component model of ionicmicrogels by integrating out the microion and particle-internalpolymer degrees of freedom, resulting in a state-dependent ef-fective pair potential and volume pressure contribution. Sec-tion III gives the essentials of the linear TPT and nonlinearPBCM methods used for calculating the equilibrium micro-gel radius as a function of concentration, backbone charge,and reservoir ionic strength, and further describes how thenet microgel charge and the electrostatic screening constantcharacterizing the pair potential are obtained. The methodsused for calculating structural and thermodynamic propertiesare discussed in Sec. IV. The diffusion and rheological prop-erties explored in this work are summarized in Sec. V, to-gether with the analytic methods for their calculation withinthe one-component model framework. In Sec. VI, resultsare presented for static properties, including the microgelswelling ratio, pair distribution function, and osmotic pres-sure, as well as for dynamic properties, including the hydro-dynamic function, collective diffusion coefficient, and low-and high-frequency suspension viscosities. Finally, Sec. VIIsummarizes our conclusions. II. EFFECTIVE ONE-COMPONENT MODEL
We describe here the employed microgel model and theessential steps of tracing out the microion and polymer-backbone monomeric degrees of freedom, leading to aneffective one-component suspension description of pseudo-microgels interacting via a state-dependent effective pair po-tential [9, 10, 12]. This potential determines the equilibriummicrostructure of the microgels and, in conjunction with astructure-independent volume grand free energy contribution,also the (osmotic) thermodynamic properties of the whole sus-pension, including its phase behavior. For conciseness, we usethe compact notation of [10], to which we refer for further de-tails, focusing here on the physical aspects.The distributions of backbone polyelectrolyte monomers,cross-linkers, and backbone charges of ionic microgels de-pend on the synthesis method. For simplicity, we assume hereuniform distributions [10]. The charged microgel backbonepolymers and cross-linkers coexist with polymer-releasedcounterions and salt ions dissolved in the solvent. For tem-peratures T > T cr higher than the lower critical solution tem-perature (LCST) T cr of the corresponding polymer solution,the microgels are collapsed into a dry state, characterized by adry radius a < a , where a is the microgel equilibrium radiusof the swollen microgel at a temperature lower than the LCST.The swollen microgel radius depends, in addition to tempera-ture and solvent quality, on the elastic properties of the back-bone network and the backbone charge, and furthermore onthe microgel concentration and background (reservoir) ionicstrength. Two methods used in this paper to calculate theswelling ratio α = a / a are described in the next section(Sec. III). Assuming that a single microgel consists of a uni-form polymer network with N mon monomers, the dry microgelradius a is well approximated by a ≈ ( N mon / φ rcp ) / a mon ,where φ rcp = .
64 is the volume fraction for random close-packing of spherical monomers and a mon is the monomer ra-dius. It is assumed here that random close-packing is theunstressed polymer backbone structure in the collapsed state[8, 20–22].We consider in this paper a monodisperse microgel sus-pension formed by N spherical microgels, each of negativebackbone charge − Ze , with e the proton charge, dispersed ina volume V of water at room temperature T . The suspen-sion is assumed to be in (Donnan) osmotic equilibrium witha 1:1 strong electrolyte reservoir of ion concentration 2 n res ,via an ideal membrane permeable to the microions and sol-vent only. The counterions dissociated from the backbonesare likewise taken as monodisperse. The microgel concentra-tion (number density) n = N / V determines the volume frac-tion φ = π a n / φ = π a n / Z > Z shouldbe viewed as net backbone valence, already accounting forthe possibility of Manning counterion condensation on poly-mer sites. Global electroneutrality implies ZN = (cid:104) N + (cid:105)−(cid:104) N − (cid:105) ,where N s = (cid:104) N − (cid:105) is the equilibrium number of monodispersecoions in the system, equal to the number N s of salt ion pairs,and (cid:104) N + (cid:105) is the equilibrium number of monovalent coun-terions. The concentration (number density) n s = N s / V ofsalt ion pairs in the suspension is determined by the equal-ity, µ ± = µ res , of the microion chemical potentials of cationsand anions, µ ± , in the suspension and the microion chemicalpotential, µ res = k B T ln (cid:0) Λ n res (cid:1) , in the reservoir, assumingequal thermal de Broglie wavelength Λ for all microions. InDonnan equilibrium, the salt pair concentration n s in the sus-pension is determined by the given reservoir salt pair concen-tration (number density) n res . A closed suspension of givensalt content can be straightforwardly mapped to an equivalentDonnan equilibrium system using an accordingly selected saltconcentration n res ≥ n s .Our starting point in deriving the one-component modelof pseudo-microgels is a semi-grand canonical description ofuniform-backbone spherical microgels with the solvent de-grees of freedom already integrated out. This amounts to de-scribing the solvent statically as a dielectric continuum of di-electric constant ε and Bjerrum length λ B = e / ( ε k B T ) , anddynamically as a Newtonian solvent of shear viscosity η . Inthis McMillan-Mayer implicit solvent picture, the semi-grandcanonical partition function of the suspension reads Ξ = (cid:104)(cid:104)(cid:104) e − β ( K + U m + U mm + U m µ + U µµ ) (cid:105) p (cid:105) µ (cid:105) m . (1)Here β = / ( k B T ) , K is the total kinetic energy of all poly-meric and ionic suspension constituents, and the angular brackets denote canonical traces over polymer (p) and center-of-mass microgel (m) coordinates, and grand-canonical tracesover the microion ( µ ) coordinates. The polymer coordinatesare particle-internal degrees of freedom associated with themotion of segments and associated fixed charges constitut-ing the cross-linked polymer chains. In the Boltzmann fac-tor, U m is the single-microgel energy, comprising both poly-meric and electrostatic self energies, U mm incorporates poly-meric and electrostatic energies of interaction between the mi-crogels, and U m µ and U µµ account, respectively, for micro-gel–microion and microion–microion interactions.Performing the trace over polymer coordinates, one obtains Ξ = e − β ( U e + F p ) (cid:104)(cid:104) e − β ( K m , µ + U mm + U m µ + U µµ ) (cid:105) µ (cid:105) m (2)where U e is the sum of the electrostatic self energies of the N microgels, which for uniformly distributed backbone chargesis U e ( a ) = N ∑ i = u e ( a ) = N (cid:18) Z e ε a (cid:19) , (3)with the equilibrium radius a of swollen microgels. Further-more, K m , µ is the translational kinetic energy associated withthe center-of-mass microgel (m) and microion ( µ ) coordi-nates.The free energy associated with the non-electrostatic poly-meric degrees of freedom of the N microgels is F p = N ∑ i = f p ( a ) . (4)We use Flory-Rehner theory [20–22] to approximate the poly-mer free energy per microgel, f p ( a ) , for a particle networkwith uniformly distributed cross-linkers that is divided into N ch chains, i.e., β f p ( a ) = N mon [( α − ) ln ( − α − ) + χ ( − α − )] ++ N ch ( α − ln α − ) , (5)where χ is the Flory solvency parameter, α the swelling ratio,and N mon the total number of polymer monomers in a micro-gel. The first term on the right-hand side is the ideal mixingentropy of microgel monomers and solvent molecules. Thesecond term accounts for polymer-solvent interactions in amean-field approximation, by neglecting interparticle corre-lations. The last term accounts for the elastic free energy forisotropic stretching of the microgel network, with the poly-mers treated as Gaussian coils. As argued in [10, 12], the ap-proximations employed here for the microgel backbone self-energies are reasonable for loosely cross-linked, uniformlystructured microgels.Tracing out in a second step the microion degrees of free-dom for fixed configuration of microgels leads to the expres-sion [10] Ξ = (cid:104) e − β H eff (cid:105) m , (6)with the effective Hamiltonian of pseudo-microgels, H eff = K m + U e + F p + E V ( n ) + U eff ( n ) , (7)where K m accounts for the translational kinetic energy of themicrogels, E V ( n ) is the microgel configuration-independentvolume energy, and U eff ( n ) is the configuration-dependenteffective N -particle interaction energy of pseudo-microgels.The latter, which incorporates electrostatic screening by thetraced-out microions, consists of the bare interaction energy, U mm , comprising the concentration-independent Coulomband elastic inter-microgel interactions, and a concentration-and temperature-dependent contribution, related to the freeenergy of microions in the presence of the microgels.With F = − k B T ln Ξ denoting the semi-grand suspensionfree energy, the pressure, p , of the multi-component sus-pension, consisting of polymer networks with charged sitesand microions, is then determined by the generalized one-component virial equation [19], p = − (cid:18) ∂ F ∂ V (cid:19) res = p V + p se + nk B T − V (cid:10) N ∑ i = r i · ∂ U eff ∂ r i (cid:11) eff − (cid:10) ∂ U eff ∂ V (cid:11) eff , (8)invoking an extra volume derivative term due to the concentra-tion dependence of U eff ( X ; n ) . Here, X = { r , · · · , r N } are thecenter-of-mass positions of microgels and p V = − ∂ E V / ∂ V isthe pressure contribution of the volume energy. There is ananother pressure contribution, p se = n ∂ [ u e ( a ) + f p ( a )] / ∂ n ,originating from the electrostatic and polymeric self-energies(se) per particle, owing to their implicit concentration depen-dence via the equilibrium radius a ( n ) . This contribution is ab-sent for incompressible particles. The angular brackets (cid:104)· · · (cid:105) eff denote the canonical average with respect to the equilibriumdistribution function, P eq ( X ) ∝ exp [ − β U eff ( X )] , of pseudo-microgels, not to be confused with the canonical microgeltrace (cid:104)· · · (cid:105) m over microgel center positions and momenta. Thevolume derivative of F in Eq. (8) is for fixed reservoir ionchemical potentials and, hence, fixed n res . The generalizedvirial equation does not suffer from ambiguities introducedwhen state-dependent pair potentials are combined in an adhoc manner with the compressibility and virial equation ofstate expressions for one-component simple liquids [28, 29].In Donnan equilibrium, the reduced microgel osmotic com-pressibility can be expressed via the Kirkwood-Buff (KB) re-lation [30, 31], k B T (cid:18) ∂ n ∂ p (cid:19) res = + n (cid:90) d r [ g ( r ; n ) − ] = S ( q → n ) , (9)solely in terms of the solvent-averaged microgel radial dis-tribution function g ( r ; n ) = g mm ( r ; n ) , which in turn is solelydetermined by the effective interaction energy U eff . In con-trast, p is not determined by U eff ( n ) alone, since it has alsoa pressure contribution arising from the volume energy. TheKB relation follows from the isothermal differential Gibbs-Duhem relation in Donnan equilibrium, d p = n d µ , where µ is the microgel chemical potential, in conjunction with therelation β S ( n ) = ( ∂ ln n / ∂ µ ) T , µ ± for the zero-wavenumberstructure factor of the microgels confined to the suspension.The radial distribution function is basically the inverse Fouriertransform of the microgel structure factor [32], S ( q ; n ) = + π n (cid:90) ∞ dr r [ g ( r ; n ) − ] sin ( qr ) qr , (10)determined in static scattering experiments as a function ofthe scattering wavenumber q . The zero-wavenumber limit of S ( q ; n ) is proportional to the osmotic compressibility.The volume and concentration derivatives, taken in Eqs. (8)and (9) respectively, are for fixed reservoir properties, i.e.,fixed n res and T , so that ∂ p / ∂ n = ∂ π os / ∂ n , where π os = p − k B T n res (11)is the osmotic pressure of the suspension, measured relativeto the reservoir pressure p res . Non-ideality contributions tothe reservoir pressure are negligible for the considered reser-voir ionic strengths of monovalent electrolyte ions. Noticethat both the generalized virial equation and the KB relationsare valid also for a non-pairwise additive U eff .For weak particle overlap, it is reasonable to assume pair-wise additive elastic forces which for swollen volume frac-tions φ (cid:46) β u H ( r ) = (cid:40) ε H (cid:0) − r a (cid:1) / , r ≤ a , r > a , (12)where r is the center-to-center separation of two particles.The Hertz soft particle radius is identified with the equilib-rium (swollen) radius a . The reduced interaction (softness)parameter ε H is determined by the single-particle elastic mod-uli, independent of temperature and particle volume [24], andit scales linearly with N ch . The pairwise-additive bare micro-gel interaction energy is thus U mm ( X ) = N ∑ i < j [ u H ( r i j ) + u C ( r i j )] , (13)where u C ( r i j ) is the Coulomb interaction energy between mi-crogels i and j at center-to-center distance r i j , modeled as theelectrostatic energy of two uniform spherical charge clouds,each of charge − Ze and radius a .To obtain the effective potential energy, U eff ( X ; n ) , of N pseudo-microgels, the concentration-dependent free energycontribution due to the traced out microions needs to be cal-culated. Assuming weak perturbation of the microion distri-bution by the uniform microgel backbone charges, this can bedone using the linear-response approximation method of Den-ton [18, 23, 24]. This method invokes a random phase approx-imation for the static response functions of a reference plasmaof pointlike-assumed microions, resulting in a linear superpo-sition of isotropic coion and counterion concentration profiles(orbitals), n ± ( | r − r i | ) , centered at the respective microgel po-sitions r i . The effective N -microgel interaction energy hasthen only two-body contributions, so that U eff ( X ; n ) = N ∑ i < j [ u H ( r i j ) + u eff ( r i j ; n )] . (14)The effective electrostatic pair potential in the linear-responseapproximation is of different functional form for overlappingand non-overlapping microgels, i.e., u eff ( r ; n ) = (cid:40) u Y ( r ; n ) , r > au ov ( r ; n ) + u H ( r ) , r ≤ a . (15)For non-overlapping pairs, the effective electrostatic pair po-tential has the functional form of a screened-Coulomb poten-tial [18, 19, 23], β u Y ( r ; n ) = λ B [ Z net ( n )] (cid:18) e κ a + κ a (cid:19) e − κ r r , (16)with net microgel valence Z net = Z ( + κ a )( κ a ) e κ a (cid:20) cosh ( κ a ) − sinh ( κ a ) κ a (cid:21) , (17)depending on the product, κ a , of the Debye screening con-stant κ and the swollen radius a . The net microgel valence isobtained from Z net ( n ) = Z − π (cid:90) a [ n + ( r ) − n − ( r )] r dr , (18)using the linear-response prediction for the equilibrium coun-terion and coion concentration orbitals n ± ( r ) . The Debyescreening constant has the form κ = πλ B ( Zn + n s ) , (19)with dependence on the microgel concentration n and salt pairconcentration n s in the suspension. The somewhat lengthyexpression for the repulsive effective electrostatic potentialof overlapping microgels, u ov ( r ; n ) , which depends on Z and κ a , is given elsewhere [18, 19, 23]. As an ultrasoft poten-tial, u ov ( r ; n ) is bounded with zero slope (no repulsion) forfull overlap of two spherical microgels and connects smoothlywith u Y ( r ; n ) at r = a , where the first derivatives are equal.Since u H ( r ) and its first derivative are zero at contact distance,also u eff ( r ) crosses over smoothly at r = a .Owing to the pairwise additivity of U eff ( n ) in the invokedlinear-response electrostatic plus Flory-Hertz elasticity ap-proximations, the generalized virial equation for the suspen-sion pressure reduces to [25] β p = n − π n (cid:90) ∞ drr g ( r ) ∂ β u eff ( r ) ∂ r + π n (cid:90) ∞ drr g ( r ) ∂ β u eff ( r ) ∂ n + β p V + β p se . (20)From the explicit expression for the volume energy E V ( n ) in the linear-response approximation, the corresponding pres- sure contribution is [19] β p V = n (cid:18) ∂ β ε V ∂ n (cid:19) res = Zn + n s + Z λ B a n (cid:20) − κ +
94 ˜ κ −
154 ˜ κ + (cid:18)
32 ˜ κ +
214 ˜ κ +
152 ˜ κ +
154 ˜ κ (cid:19) e − κ (cid:21) , (21)where ˜ κ = κ a , and ε V = E V / N is the volume energy per parti-cle. Note that the (reduced) kinetic pressure of the microions, Zn + n s , is included in p V .In calculations of the osmotic compressibility, one can takeadvantage of a theorem by Henderson [33], asserting that fora one-component system with only pairwise interactions, foreach considered thermodynamic state (concentration n ) thereis a one-to-one correspondence between g ( r ; n ) and the under-lying pair potential, up to an irrelevant additive constant forthe latter. As thoroughly discussed in [29, 31], at given con-centration n and temperature, the osmotic compressibility canthus be obtained also from the concentration derivative of thesuspension pressure, p OCM , for a fictitious system with state-independent pair potential u ( r ) = u eff ( r ; n , T ) . Explicitly, (cid:18) ∂ p ∂ n (cid:19) res = (cid:18) ∂ p OCM ∂ n (cid:19) u eff , (22)where p OCM is the one-component model (OCM) pressure ofthe fictitious system, given by the right-hand side of Eq. (20)without volume pressure p V and without the integral invokingthe concentration derivative of u eff ( r ; n ) . The concentrationderivative on the right-hand side of Eq. (22) is taken for fixed u eff , by discarding any concentration dependence of the effec-tive pair potential, which amounts to keeping κ a and Z net fixedto their values at the considered concentration. A consistencytest of the approximations used in calculating u eff ( r ; n ) , theequilibrium radius a , and g ( r ; n ) follows from comparing nu-merical values for the reduced osmotic compressibility iden-tified by S ( q → n ) with values obtained from Eq. (22). III. MICROION-INDUCED DESWELLING
So far, we have treated the microgel radius as a given quan-tity. However, as pointed out in the introduction, it is actuallya state-dependent thermodynamic variable whose equilibriummean value, a , for temperatures T < T cr in the swollen state, isdetermined from minimization of the semi-grand free energy F ( a t ) of the suspension with respect to trial radius values a t .The necessary condition for determining a is thus ∂ F ∂ a t (cid:12)(cid:12)(cid:12)(cid:12) N , Z , res = a t = a . In addition to the dry radius a and the electrostaticparameters Z and n res , the elasticity-related Flory-Rehner andHertz potential parameters, χ , N m , N ch , and ε H , are kept con-stant in taking the size derivative. In this way, a is determinedas a function of the control parameters n , Z , and n res for fixedtemperature and microgel elastic properties.The minimization of F with respect to a t is equivalent to themechanical requirement that the intrinsic pressure difference[8, 9], Π ( a t ) = Π g ( a t ) + Π e ( a t ) , (24)between the interior and exterior of a single microgel is zeroat thermodynamic equilibrium, where a t = a . The radius-dependent gel-elasticity pressure contribution, Π g , due to sol-vency, elasticity, and mixing entropy of individual microgelnetworks and the Hertz elastic pair interactions, is given inthe present microgel model by Π g ( a t ) = − ∂∂ v t (cid:104) f p ( a t ) + n (cid:10) u H ( r ; a t ) (cid:11) eff (cid:105) , (25)where v t = π a / Π ( a t ) is Π e ( a t ) = − ∂∂ v t (cid:104) u e ( a t ) + ε V ( a t ) + n (cid:10) u eleff ( r ; a t ) (cid:11) eff (cid:105) , (26)where u e ( a t ) is the electrostatic self energy of the uniformbackbone microgel charge, and u eleff ( r ; a t ) is the effective elec-trostatic pair potential [Eq. (15)]. For conditions where over-lap distances are very unlikely, the Hertz potential energy doesnot contribute to Π ( a t ) and the canonical average (cid:104)· · · (cid:105) eff overthe center positions of pseudo-microgels of radius a t is deter-mined alone by the Flory-Rehner and electrostatic parameters.The equilibrium radius is determined by the competition be-tween Π g , which is negative for a t sufficiently larger than a favoring deswelling, and the positive-valued Π e ( a t ) favoringswelling. Physically, the microion distribution in the micro-gel interior and the self-repulsion of the charged sites of thepolymer backbone network generate an outward electrostaticpressure that swells the macroion. This swelling is limitedby the inward elastic restoring forces due to the cross-linkedpolymer gel. In equilibrium, the balance between these op-posing pressures determines the microgel size.The microgel surface plays here the role of a mobile semi-permeable membrane, permeable to microions and solventonly where the outer and inner pressures balance to Π ( a ) = π os = p − p res , acting across a (mentally pic-tured) fixed semi-permeable membrane separating the suspen-sion from the microion reservoir.In the following, we describe and contrast two methodsused for calculating the state-dependent equilibrium radius a as a function of microgel concentration, backbone valence,and reservoir salt concentration. The first method makes di-rect use of Eq. (23) and of the one-component multi-centerpicture of pseudo-macroions interacting electrostatically bythe linear-response effective pair potential u eff ( r ) , using athermodynamic perturbation theory (TPT) approximation forthe semi-grand free energy. The second method invokes aspherical cell model (CM) approximation for the semi-grand free energy of a single macroion with nonlinear Poisson-Boltzmann (PB) distributions of microions, referred to ac-cordingly as the PBCM method. The two methods differ inthe manner in which they treat inter-microgel electrostatic in-teractions and correlations. A. Thermodynamic Perturbation Theory
In the thermodynamic perturbation theory (TPT) method,the equilibrium radius a is obtained by minimizing the semi-grand free energy per microgel [10], F ( a t , n ) N = u e ( a t ) + ε V ( a t ) + f p ( a t ) + f ex ( a t , n ) , (27)with respect to trial radius values a t . We have disregarded herethe kinetic (ideal gas) free energy contribution, ln (cid:0) Λ n (cid:1) − F / N where Λ m denotes the thermal de Broglie wavelengthof microgels, since it is independent of a t . The excess semi-grand free energy per microgel, f ex ( a t , n ) , is due to the ef-fective interactions between the pseudo-microgels. Provided U eff ( X ; n ) is pairwise additive, f ex ( a t , n ) is exactly given bythe charging-process ( λ -integration) expression [25, 29] f ex ( a t , n ) = n (cid:90) d r [ u H ( r ; a t ) + u eff ( r ; a t , n )] × (cid:90) d λ g λ ( r ; a t , n ) , (28)irrespective of whether the pair potential is state-dependent ornot. Here, g λ ( r ; a t , n ) is the rdf corresponding to the pair po-tential λ [ u H ( r ) + u eff ( r )] at charging fraction λ , which rangesfrom g λ ( r ) = λ = λ =
1. In principle, the above two-step integral expressioncan be used in Eq. (27) to determine a by minimization of F ( a t , n ) . Moreover, it provides another route to determine thesuspension pressure p and the osmotic compressibility fromthe first and second volume derivatives of F ( a , n ) .To avoid the cumbersome double integration involving thecalculation of a large number of rdfs for different values of λ ,we approximate f ex ( a t , n ) instead using the first-order pertur-bation expression [25] given by the right-hand side of f ex ( a t , n ) ≤ min ( d ) (cid:110) f EHS ( d , n )+ π n (cid:90) ∞ d drr g EHS ( r ; d , n ) u eff ( r ; a t , n ) (cid:111) , (29)which invokes a reference system of effective hard spheres(EHS) of diameter d , rdf g EHS ( r ; d , n ) , and free energy perparticle f EHS ( d , n ) . The EHS free energy, f EHS , is accuratelydescribed by the analytic Carnahan-Starling free energy ex-pression, and the EHS rdf by the semi-analytic Percus-Yevickresult [34] with Verlet-Weis correction [25]. The above per-turbation expression provides an upper bound to the actualexcess free energy f ex ( a t , n ) for all values of the effective di-ameter d , as follows from the Gibbs-Bogoliubov inequality[25]. The equilibrium radius a results from the (double) min-imization of F ( a t , n ) / N in Eq. (27) with respect to a t , aftersubstitution of the right-hand side of Eq. (29) for the excessfree energy minimized with respect to d >
0. For u eff and ε V ,we use the analytic linear-response expressions of Denton etal. [18, 19], and for f p ( a t ) the Flory-Rehner expression givenin Eq. (5).In the TPT, the suspension pressure can be computed from f ( a , n ) = F ( a , n ) / N using the thermodynamic relation, p = n (cid:18) ∂ f ( a , n ) ∂ n (cid:19) res , (30)where the concentration dependence of a ( n ) must be ac-counted for, giving rise, in particular, to the extra pressurecontribution p se . In taking the concentration derivative, theelectroneutrality condition n s = (cid:104) N + (cid:105) / V − nZ must be main-tained for given Z . The suspension salt pair concentration n s ,which affects κ ( n , n s ) , and hence the range of the effectivepair potential in the TPT expression for f ex ( a , n ) in Eq. (29),is determined, in turn, from equating the microion chemicalpotentials in suspension and reservoir, using (cid:104) N − (cid:105) = N s , ac-cording to ∂∂ n s (cid:104) n ( ε V ( a ) + f ex ( a , n )) (cid:105) n = k B T ln (cid:0) Λ n res (cid:1) . (31)The TPT method was successfully tested in earlier worksfor deswelling ionic microgels [10], incompressible ionic mi-crogels [19], and impermeable charged colloids [35]. Themethod self-consistently incorporates effective microgel pairinteractions for low to moderately high Z , where linear-response theory can be used. B. Poisson-Boltzmann Cell Model
The PBCM applies to suspensions of ionic microgels,where on average around each microgel there is a region voidof others [36]. This condition requires sufficiently strong andlong-ranged electrostatic repulsion between the microgels andconcentrations small enough that particle overlap is unlikely.In this case, a Wigner-Seitz (WS) cell tessellation can be used,with each WS cell subsequently approximated by an overallelectroneutral spherical cell of radius R = ( / π ) / n − / ,containing a single spherical microgel of radius a t at its center.In Donnan equilibrium, the cell is in osmotic contact with a1:1 strong electrolyte reservoir of salt pair concentration 2 n res .In the PBCM, the radially symmetric concentration profiles n ± ( r ) of the pointlike monovalent microions dissolved in astructureless dielectric solvent of Bjerrum length λ B are de-scribed in a mean-field way by the Boltzmann distributions, n ± ( r ) = n res e ∓ Φ ( r ) , where Φ ( r ) = ψ ( r ) e / k B T is the reducedform of the total electrostatic potential ψ ( r ) due to all chargesin the cell. As in the TPT method, polarization and imagecharge effects are disregarded, which can be justified by thehigh solvent content of weakly cross-linked, swollen micro-gels. While the cell model focuses on only a single microgel,with the semi-grand suspension free energy being N times thatof the cell, the presence of other microgels is implicitly ac-counted for through the cell radius R and the associated (trial)volume fraction φ t = ( a t / R ) . Assuming, as in the TPT, a uniform backbone charge distri-bution inside each microgel, the electrostatic potential in thecell region 0 < r < R is obtained from solving the nonlinearPoisson-Boltzmann (PB) equation, Φ (cid:48)(cid:48) ( r ) + r Φ (cid:48) ( r ) = κ sinh Φ ( r ) + Z λ B a , < r ≤ a t κ sinh Φ ( r ) , a t < r ≤ R , (32)where κ = πλ B n res is the square of the reservoir Debyescreening constant. The solution for Φ ( r ) is uniquely deter-mined by the boundary conditions, Φ (cid:48) ( ) = = Φ (cid:48) ( R ) , on theelectric field at the cell center and edge, and by the continuityconditions, Φ ( a − t ) = Φ ( a + t ) and Φ (cid:48) ( a − t ) = Φ (cid:48) ( a + t ) , at the mi-crogel surface. Once Φ ( r ) , and hence the microion concen-tration distributions, are determined by numerically solvingEq. (32) for given boundary conditions and microgel trial ra-dius a t , the intrinsic osmotic pressure Π ( a t ) in the PBCM fol-lows from Eqs. (25) and (26) taken for (cid:104) u H (cid:105) eff = = (cid:104) u eff (cid:105) eff ,and for ε V replaced by (cid:104) u m µ ( r ) (cid:105) µ , i.e., by the electrostatic in-teraction energy between the uniform central microgel back-bone charge and pointlike microions, weighted by the mi-croion number density profiles V R n ± ( r ) and averaged over thecell volume V R = ( π / ) R . Considering the variation of theelectrostatic component of the free energy with respect to themicrogel radius leads to an exact statistical mechanical rela-tion for the electrostatic pressure [9]: β Π e ( a t ) v t = Z λ B a t (cid:18) Z − (cid:104) N + (cid:105) + (cid:104) N − (cid:105) + (cid:104) r (cid:105) + − (cid:104) r (cid:105) − a (cid:19) , (33)where (cid:104) N ± (cid:105) = π (cid:90) a t n ± ( r ) r dr (34)and (cid:104) r (cid:105) ± = π (cid:90) a t n ± ( r ) r dr (35)are the mean numbers of internal microions and the secondmoments of the interior microion number density profiles, re-spectively.Using Eq. (5), the polymer gel contribution to the intrinsicosmotic pressure for trial radius a t is [9] β Π g ( a t ) v t = − N m [ α ln ( − α − ) + χα − + ] − N ch ( α − / ) . (36)According to Eq. (24), the equilibrium microgel radius is ob-tained from setting the sum of the intrinsic pressure contribu-tions in Eqs. (33) and (36) equal to zero.Once the microgel equilibrium radius is determined, thepressure in the cell model due to the mobile microions, p µ ,follows from the contact theorem [37] β p µ = n + ( R ; a ) + n − ( R ; a ) , (37)i.e., from the microion concentrations at the cell edge, wherethe electric field vanishes due to overall electroneutrality. Inaddition to the kinetic (ideal gas) microgel pressure nk B T ,there are microgel-correlation-induced pressure contributionsto the suspension pressure, p , not accounted for in the cellmodel. Therefore, p µ can differ significantly from p , ex-cept for relatively low reservoir salt concentrations, where thedominant number of backbone-released counterions ( ZN (cid:29) N s ) contribute most to p [31]. As shown elsewhere [38], in thiscounterion-dominated regime, where Zn (cid:29) n res , the dom-inant contribution to p in Eq. (20) stems from the volumeenergy-related pressure p V . Akin to the cell model pressure p µ , the pressure contribution p V arises from the microions inthe presence of fixed microgels. The positive-valued OCMpressure on the right-hand side of Eq. (20) is nearly compen-sated at low salinity by the negative-valued pressure contri-bution from the concentration derivative of u eff ( r ; n ) . Whilethis compensation is observed for non-permeable charge-stabilized colloids [38], it likely holds also for ionic micro-gels. One should not infer from this compensation, how-ever, the practical identity of p µ and p V in the counterion-dominated concentration region, since the underlying models,i.e., spherical cell versus multi-center model, and the respec-tively invoked approximations (linear-response theory versusPB aproximation) in the pressure calculations are distinctlydifferent.In the cell model, the net microgel valence Z net is calculatedby means of Eq. (18) using the microion number density pro-files n ± ( r ; a ) , and the suspension salt pair concentration n s ,by integrating the coion (anion) profile over the cell volumeaccording to n s = π V R (cid:90) R n − ( r ; a ) r dr . (38)While TPT is self-consistently linked to the effective pairpotential u eff ( r ) in Eq. (15), characterized by Z net and κ a forgiven backbone valence Z , such a direct link does not exist inthe single-microgel PBCM, which does, however, incorporatea nonlinear electrostatic response of the microions that is ne-glected in the TPT. However, an ad hoc link between PBCMand the linear-response u eff ( r ) is readily established, for given Z , by identifying Z net and κ a in the no-overlap Yukawa po-tential in Eq. (16) with the PBCM-calculated values Z ∗ net and κ ∗ a ∗ , respectively, where ( κ ∗ ) = πλ B ( nZ ∗ + n ∗ s ) (39)and the asterisk labels PBCM-calculated properties. An ap-parent backbone valence Z ∗ is defined here as a function of Z ∗ net and κ ∗ a ∗ by Z ∗ net = Z ∗ ( + κ ∗ a ∗ )( κ ∗ a ∗ ) e κ ∗ a ∗ (cid:20) cosh ( κ ∗ a ∗ ) − sinh ( κ ∗ a ∗ ) κ ∗ a ∗ (cid:21) , (40)which when used in the expression for the overlap elec-trostatic potential u ov ( r ) , according to the substitution { Z , a , κ } → { Z ∗ , a ∗ , κ ∗ } , maintains the continuity of the effec-tive potential and its first derivative at r = a . Substitution of Z ∗ into Eq. (39) gives an implicit equation for κ ∗ , which canbe solved iteratively. For lower backbone valences Z ≤ Z ∗ is close to Z , so that the latter can be used instead as inputin Eq. (39).Most results presented here are for ionic microgel sys-tems with electrostatic coupling strengths Γ el ≡ Z net λ B / a (cid:46) Z and κ , whichcan be obtained, e.g., from linearization of the potential Φ ( r ) in the cell model with respect to its value at the cell bound-ary [8], or with respect to the cell volume averaged potentialvalue. While microgel charge-renormalization is not in thefocus of this paper, in the framework of PBCM we use it toassess the concentration shift at a fluid-solid freezing transi-tion caused by the deswelling of strongly charged microgels. IV. THERMODYNAMICS AND STRUCTURE
Once the mean radius a is determined for given system pa-rameters n , Z , n res , χ , N mon , N ch and ε H , we are in the positionto calculate thermodynamic, structural, and dynamic proper-ties of the one-component suspension of pseudo-microgels in-teracting via the effective pair potential in Eq. (15). As weshow below in the Results section (Sec. VI), the TPT andPBCM predictions for a are quantitatively different, as re-flected in the calculated static and dynamic properties.Our methods for calculating dynamic properties of the mi-crogel suspension require the static structure factor, S ( q ) ,of microgels and the associated radial distribution function, g ( r ) , as the only input. Since u eff ( r ; n ) is purely repulsive,we can use the thermodynamically self-consistent Rogers-Young (RY) integral-equation scheme [25] for calculatingthese structural properties. This hybrid scheme, which usesa closure mixing function interpolating between the hyper-netted chain (HNC) and Percus-Yevick (PY) integral-equationschemes [25], is known from comparisons with computer sim-ulation data to make accurate structural predictions for a vari-ety of repulsive interaction potentials, including the screened-Coulomb potential [39, 40] used to model non-overlappingionic microgels. The mixing parameter α in the RY mixingfunction is determined self-consistently from enforcing equal-ity of the microgel osmotic compressibility obtained from theone-component compressibility and virial equation of states,respectively, i.e., from demanding1 S ( q = α ) = β (cid:18) ∂ p OCM ( α ) ∂ n (cid:19) u eff (41)in accord with Eqs. (9) and (22).With α determined self-consistently, the pressure p can becalculated in the TPT using the thermodynamic relation inEq. (30). Alternatively, the pressure can be calculated fromEq. (20) using the RY- g ( r ) as input, in conjunction with thevolume energy-related pressure contribution p V in Eq. (21).Differences in the predictions for p by the two routes reflectthe accuracies of the approximations going into the TPT andRY methods. V. DYNAMIC PROPERTIESA. General Theory
The employed methods for calculating dynamic proper-ties of microgel suspensions are based on the one-componentmodel of pseudo-microgels interacting by the state-dependenteffective pair potential in Eq. (15). With regard to dynamicproperties, different colloidal time regimes need to be distin-guished [32, 41].We focus mainly on the colloidal short-time regime, char-acterized by correlation times t for which τ B (cid:28) t (cid:28) τ I holds,i.e., for times t well separated from the long-time regimewhere t (cid:29) τ I . Here, τ B = M / ( πη a h ) , with M the parti-cle (microgel) mass and a h the hydrodynamic particle radius,is the particle momentum relaxation time characterizing thetime range where momentum changes (i.e. inertia) matters.Moreover, τ I = a / d is the structural relaxation time, where d = k B T / ( πη a h ) is the Stokes-Einstein-Sutherland trans-lational free diffusion coefficient of a spherical colloidal par-ticle. Moreover, η is the shear viscosity of the suspendinglow-molecular-weight Newtonian solvent (i.e., water). Ow-ing to the low hydrodynamic permeability of (ionic) micro-gels [42], we identify for simplicity the hydrodynamic radius a h of the microgels with the equilibrium radius a determinedin the TPT and PBCM, respectively.During times t (cid:28) τ I , over which particle displacements byBrownian motion are minuscule compared to the particle ra-dius, short-time dynamic properties are influenced solely bythe inter-microgel hydrodynamic interactions (HIs), which arequasi-instantaneously transmitted by intervening solvent-flowperturbations. Short-time transport properties can thus becalculated as genuine equilibrium averages of configuration-dependent hydrodynamic mobilities. The non-dynamic inter-actions embodied in u eff ( r ; n ) are only indirectly influentialthrough their effect on the equilibrium microstructure encodedin g ( r ) and S ( q ) . Long-time transport properties, such as thezero-frequency, steady-shear suspension viscosity η and thelong-time self-diffusion coefficient d l , with the latter coeffi-cient being proportional to the long-time slope of the parti-cle mean-squared displacement, are influenced additionallyby u eff ( r ; n ) via non-instantaneous caging (i.e., memory) ef-fects, whose description requires, in general, more elaboratecalculations.The short-time diffusion of microgels is commonly probedexperimentally by measuring the q -dependent dynamic struc-ture factor S ( q , t ) using dynamic light scattering. At shorttimes, S ( q , t ) decays exponentially according to [32, 40, 43] S ( q , t (cid:28) τ I ) = S ( q ) exp {− q D ( q ) t } , (42)where D ( q ) is the wavenumber-dependent short-time diffu-sion function characterizing the decay of concentration fluc- tuations of wavelength 2 π / q . A statistical-mechanical expres-sion for D ( q ) follows from the generalized Smoluchowski dif-fusion equation of interacting Brownian particles in the formof the ratio [32, 40, 43], D ( q ) = d H ( q ) S ( q ) , (43)where H ( q ) is the so-called hydrodynamic function given bythe equilibrium average [32], H ( q ) = (cid:42) N µ q N ∑ l , j = q · µ l j ( X ) · q e i q · ( r l − r j ) (cid:43) eff , (44)over the positional configurations X of the microgels.Here, k B T µ = d and µ l j ( X ) are the translational N -spheremobility tensors linearly relating the hydrodynamic force ona sphere j to the instant velocity change of sphere l causedby the solvent-transmitted HIs. These tensors depend on theinstantaneous configuration, X , of the N microgel centersthrough the specified hydrodynamic surface boundary con-ditions. The positive-valued function H ( q ) is a measure ofthe influence of HIs on short-time diffusion over the lengthscale ∼ / q . In the (hypothetical) case of hydrodynamicallynon-interacting particles, H ( q ) ≡
1, independent of q and theparticle concentration. Deviations of H ( q ) from the infinitedilution value of one thus hallmark the influence of HIs.According to H ( q ) = d s d + H d ( q ) , (45)the hydrodynamic function is the sum of a self-part equal tothe short-time self-diffusion coefficient d s (in units of d ),quantifying the initial slope of the particle mean-square dis-placement, and a wavenumber-dependent distinct part, H d ( q ) ,accounting for hydrodynamic cross correlations between themicrogels. The latter part decays to zero at large q . For large qa (cid:29)
1, the hydrodynamic function becomes thus equal to d s / d , while for small wavenumbers qa (cid:28) K ( n ) = H ( q → n ) ofa homogeneous suspension subjected to a weak (gravitational)force field. The associated short-time collective diffusion co-efficient, d c ( n ) = d ( n ) K ( n ) S ( q → n ) = d ( n ) K ( n ) k B T ( ∂ n / ∂ p ) res , (46)is even for a concentrated suspension only slightly larger (by afew percent) than the long-time collective diffusion coefficientappearing in the macroscopic Fickean constitutive law, whichlinearly relates the particles current to the concentration gra-dient [40]. This behavior should be distinguished from self-diffusion, where d l ≈ . × d s right at the fluid-crystal freez-ing transition point of a three-dimensional colloidal system[44, 45].An important feature distinguishing (ionic) mi-crogels from impermeable solid particles is that d ( n ) = k B T / ( πη a ( n )) = d dry0 / α ( n ) depends on con-centration. Here, d dry0 is the Stokes-Einstein diffusion0coefficient of collapsed (dry) microgels, and α ( n ) = a ( n ) / a the swelling ratio at concentration n . In our calculationsof diffusion and rheological properties, we identify thehydrodynamic microgel radius for simplicity with the ther-modynamic mean particle radius a ( n ) as obtained by theTPT/PBCM methods. While on first sight this appears to bea severe approximation owing to the solvent permeabilityof weakly cross-linked (ionic) microgels, calculations showthat the hydrodynamic penetration depth related to the Darcypermeability of microgels is actually quite small so thatsolvent-permeability effects can be disregarded, as they playa noticeable role only at high concentrations [42].A non-diffusional, rheological short-time property char-acterizing the microgel suspension as a whole is the high-frequency viscosity, η ∞ , for low shear rates. This property lin-early relates the average deviatoric suspension shear stress tothe applied rate of strain in a low-amplitude, oscillatory shearexperiment at frequencies ω (cid:29) / τ I , where shear-inducedperturbations of the microstructure away from the equilibriumspherical symmetry are negligible. Experimentally, η ∞ can bedetermined using a torsional rheometer operated at high fre-quencies and low amplitudes. The high-frequency viscosityis a purely hydrodynamic property, whose statistical physicsexpression is given, owing to isotropy, by (see, e.g., [46]) η ∞ = η + lim q → ∑ α , β = (cid:68) V N ∑ l , j = µ ddl j , αββα ( X ) e i q · ( r l − r j ) (cid:69) eff , (47)where µ ddl j , αββα are the Cartesian components of the fourth-rank dipole-dipole hydrodynamic tensor µ ddl j relating the sym-metric hydrodynamic force dipole moment tensor of microgelsphere l to the rate of strain tensor evaluated at the center of asphere j . The zero-wavenumber limit is taken after the ensem-ble averaging over a macroscopic system, guaranteeing in thisway convergence of the integrals following from the averag-ing over the spatially slowly decaying hydrodynamic tensors[47].As an important colloidal long-time property, we computealso the low shear rate, zero-frequency viscosity η > η ∞ , mea-sured in a suspension subjected to steady-state weak shearflow. The viscosity η is the sum [48], η = η ∞ + ∆ η , (48)of η ∞ and a shear stress relaxation contribution denoted ∆ η .The latter contribution is related to the additional dissipationin the suspension originating from stress relaxations of theshear-perturbed next-neighbor particle cages formed aroundeach microgel, and it is influenced both by direct and hydro-dynamic interactions. The viscosity part ∆ η can be calculatedbased on an exact Green-Kubo relation for the time integralof the equilibrium stress time auto-correlation function whereHIs are included [48].In the employed one-component model of ionic pseudo-microgels, electro-kinetic effects due to a non-instantaneousdynamic response of the microion clouds formed inside andoutside the microgels are disregarded. These effects tend tolower d c and d l , and to increase η , but in general by only small amounts. Electrokinetic effects on diffusion and rheology areof secondary importance, in particular, when non-dilute sus-pensions are considered and when the microions are smallcompared to the microgels, which is commonly the case. B. Methods of Calculation
For the calculation of H ( q ) , we use the well-established an-alytic BM-PA scheme [49]. This scheme is a hybrid of thesecond-order Beenakker-Mazur method (BM), used here forthe wavenumber-dependent distinct part H d ( q ) , and the hydro-dynamic pairwise-additivity approximation (PA) used for the q -independent self part d s / d . The BM-PA scheme combinesthe advantages of the BM and PA methods. It requires themicrogel S ( q ) and g ( r ) as its only input, for which the RY re-sults based on u eff ( r ) , and the TPT/PBCM results for a ( n ) andhence for d ( n ) , are used. The overall good accuracy of theBM-PA scheme was assessed by the comparison with elabo-rate dynamic simulation results, where many-particles HIs areaccounted for, and with experimental H ( q ) data, for a varietyof colloidal model systems, including solvent-permeable hardspheres (non-ionic microgels), charge-stabilized rigid parti-cles, and globular proteins exhibiting short-range attractionand long-range repulsion [39, 40, 49, 50]. For details aboutthe employed BM-PA method, we refer to [49].For the here considered low-salinity microgel suspensions,which show counterion-induced deswelling, we calculate thehigh-frequency viscosity η ∞ using a modified Beenakker-Mazur mean-field method described in [49]. In this semi-analytic method invoking one-dimensional integrals only,many-particles HIs are approximately accounted for, but lu-brication is disregarded. Lubrication is irrelevant, however,for solvent-permeable microgels. Just as the BM-PA schemefor H ( q ) , the modified BM method for η ∞ has S ( q ) as its onlyinput. The modified BM expression for η ∞ is [49] η ∞ η = + φ ( + φ ) − λ + λ + λ , (49)with so-called zeroth and second-order BM viscosity contri-butions, λ ( φ ) and λ ( φ ) , respectively, whose explicit formsare given in [49]. As explained in detail in this reference,the invoked modification of the standard BM expression for η ∞ / η is the subtraction of the structure-independent BM part1 / λ , and its replacement by the structure-independent pair-wise additive viscosity contribution 1 + . φ ( + φ ) , whichis known to give the dominant contribution at low salinityand small volume fractions. The modified BM expression isin very good agreement with Stokesian Dynamics simulationdata for the high-frequency viscosity of low-salinity charge-stabilized suspensions, even up to the freezing transition con-centration.As noted above, the calculation of the shear relaxation con-tribution ∆ η to the zero-frequency viscosity, η = η ∞ + ∆ η , ismore demanding since it is explicitly influenced by direct andhydrodynamic interactions. Starting from an exact but formalGreen-Kubo relation for ∆ η , mode-coupling theory (MCT)1integro-differential equations with HIs included have been de-rived for its approximate calculation, whose numerical solu-tion is quite involved. We use therefore a simplified MCT the-ory expression for ∆ η , constituting the first-iteration step inthe self-consistent numerical solution of the MCT equations.This simplified MCT expression is [48] ∆ ηη = π (cid:90) ∞ dy y ( S (cid:48) ( y )) S ( y ) H ( y ) , (50)where y = qa and S (cid:48) ( y ) = dS ( y ) / dy . HIs enter here onlythrough the dynamic structure factor S ( q , t ) , which in turn isapproximated by its short-time form given by the right-handside of Eq. (42) involving H ( q ) . Since for correlated particles S ( q , t ) decays more slowly than exponentially at longer times, ∆ η is somewhat underestimated by Eq. (50), as compared tothe fully self-consistent MCT viscosity solution. This under-estimation becomes more pronounced at higher φ . VI. RESULTS
To analyze the influence of counterion-induced deswellingon thermodynamic, structural, and dynamic properties ofionic microgel suspensions, and to make contact with a recentstudy by Weyer et. al. [10], in which TPT results for the meanmicrogel radius were compared against computer simulationsfor salt-free systems, we use the following system parame-ters, corresponding to aqueous suspensions at lower salinity:solvent Bjerrum length λ B = .
714 nm (i.e., water at temper-ature T = K ), backbone valences Z = , a =
10 nm, monomer number per mi-crogel N mon = × , polymer chain number per microgel N ch = χ = .
5, and Hertz softnessparameter ε H = . × . For the 1:1 electrolyte reservoirconcentration, we use c res = µ M, if not stated otherwise,so that n res = c res N A , where N A is the Avogradro number. Val-ues of the dry volume fraction φ = ( π / ) n a in the rangefrom 2 × − − × − are considered.Note here that φ ∝ n has the meaning of a dimensionlessmicrogel concentration. The Debye screening length, 1 / κ , inEq. (19) attains values from 40 − . φ varied from 0 . − .
05. It is noteworthythat, for most of the considered suspensions, κ is determinedby the mean concentration, Zn , of monovalent counterions re-leased from the microgel polymer backbone, which is signifi-cantly higher than the salt pair concentration n s .The electrostatic repulsion between the microgels is quan-tified by the reduced electrostatic coupling strength Γ el ≡ Z net λ B / a , which in the present study is in the range of 1-9, comparable to values for typical ionic microgel systems[4, 5, 51]. The electrostatic repulsion between the micro-gels is here strong enough that configurations of microgelsthat are in contact or overlapping are very unlikely, such that g ( r ≤ a ) ≈
0. On the other hand, nonlinear screening effectsare in most cases weak enough that the linear TPT method canbe used for determining a , in addition to the PBCM method. A. Equilibrium Radius Predictions
In the following, we analyze TPT and PBCM predic-tions for the concentration-dependent microgel swelling ratio α ( φ ) = a ( φ ) / a , the suspension salt concentration n s ( φ ) ,and the electrostatic coupling strength Γ el ( φ ) , where φ ∝ n is the non-dimensional microgel concentration. r / a [ n + ( r )- n - ( r ) ] a [ n + (r)-n - (r) ] a r / a n _ ( r ) a f = 0.03 f = 0.03 f = 0.005 f = 0.005 f = 0.005 Z = 500 a = 10 nm c res = n res / N A = 100 m M FIG. 1: PBCM predictions for radial profile of microion charge den-sity ρ el ( r ) in units of 1 / ( a e ) , versus radial distance r from center ofcell (units of dry radius a ) for microgel valence Z = a =
10 nm, reservoir salt concentration c res = µ M, and dry vol-ume fractions φ = .
03 and 0 .
005 (red and black solid curves). In-set: Reduced coion number density n − ( r ) a (dashed-dotted curves).All curves terminate at the cell radius r = R = a / φ / . Vertical linesegments indicate equilibrium swelling ratio, α = a / a , computedfrom zero balance of intra-particle pressure contributions in Eqs. (33)and (36). The physical mechanism leading to counterion-induceddeswelling with increasing concentration can be reasoned onthe basis of Fig. 1, showing PBCM results, at two differentconcentrations for the radial dependence of the (reduced) totalmicroion charge density, ρ el ( r ) = [ n + ( r ) − n − ( r )] e , and of thecoion concentration n − ( r ) (displayed in the inset) inside andoutside of a negatively charged microgel centered at r = φ , the counterions con-stitute the dominant microion species where n + ( r ) (cid:29) n − ( r ) ,and hence ρ el ( r ) ≈ n + ( r ) e holds for the total microgel chargeconcentration inside the cell up to its boundary at radius R = a / φ / , where the curves in Fig. 1 terminate.One clearly notices that both the equilibrium microgel ra-dius a , marked by the vertical line segments in the figure,and the cell radius R decrease with increasing concentration.The counterion concentration profile rises with increasing sys-tem concentration, while the coion concentration profile falls.With increasing concentration, the volume exterior to the mi-crogels is reduced, making it less favorable (entropically) forcounterions to reside outside the oppositely charged micro-gel backbone region. Consequently, a fraction of the outsidecounterions permeates into the backbone region, thereby low-ering the expansive intrinsic PBCM pressure contribution Π e Π g is established at a smaller equilibrium radius. Theenhanced counterion permeation of microgels with increas-ing concentration is reflected in the lowering of the net mi-crogel valence Z net , defined in Eq. (18), which for the back-bone valence Z =
500 is given by Z net =
223 at φ = . Z net =
154 at φ = .
03. The counterion-induceddeswelling becomes weaker with increasing salt concentra-tion, which causes a flattening of the microion concentrationprofiles across the microgel surface.In the PBCM method, the mean salt concentration n s in thesuspension is obtained by integrating the coion concentrationprofile over the cell volume according to Eq. (38). In the TPTmethod, n s is computed using the equality of the chemicalpotentials of the microions in the suspension and reservoir.In Donnan equilibrium, n s is a state-dependent quantity. TheTPT and PBCM predictions for the concentration dependenceof n s are depicted in Fig. 2 (red and black curves, respec-tively), for reservoir microion concentration c res = µ Mand backbone valences Z = n s with increasing φ , and hence with in-creasing number of backbone-released counterions, is due toan increasing expulsion of salt ion pairs into the reservoir, ne-cessitated to maintain global electroneutrality in the suspen-sion.At high dilution, φ →
0, where the concentration of saltcounterions greatly exceeds the concentration of backbone-released counterions, the exact limit n s → n res is recoveredby both methods. For the moderately high valences Z = Z =
200 considered here, the TPT and PBCM curves for n s ( φ ) in Fig. 2 lie close to each other, but with a slightlystronger salt expulsion predicted in the PBCM. Pronounceddifferences are observed for the high valence Z =
500 andintermediate φ , where the concentration, Zn , of backbone-released counterions is comparable to the salt-counterion con-centration. While the PBCM predicts a decreasing n s withincreasing Z , in accord with physical expectation, this trend isreversed for φ (cid:46) .
05 by the TPT curve for Z = φ . The PBCM accounts for nonlinear elec-trostatic effects, but not for inter-microgel correlations, whichthe TPT accounts for on a linear level.The inset of Fig. 2 displays TPT and PBCM predictions forthe Debye screening constant κ in Eq. (19), which on the scaleof the inset are practically equal. In dimensionless form, thescreening constant is ( κ a ) = ( κ c a ) + ( κ s a ) = φ Z λ B a + πλ B n s a . (51)The first term on the right-hand side is the contribution bythe backbone-released counterions (subscript c). The secondterm, proportional to n s , is the salt-ion contribution (subscripts). This splitting of κ into released-counterion and salt-ion contributions allows to identify the counterion-dominatedregime by the condition κ c (cid:29) κ s and the salt-dominatedregime by κ c (cid:28) κ s . At very low microgel concentrations, f n s a f k a n res a N mon =2e5, N ch =100, c =0.5, B/k B T =1.5e4 c res =100 m M, a =10 nm, Z = 100, 200, 500 FIG. 2: Reduced suspension salt concentration n s a versus micro-gel concentration φ . Inset: Reduced Debye screening constant κ a .TPT predictions are in red, and PBCM predictions in black for back-bone valence Z =
100 (dotted), 200 (solid), and 500 (dash-dotted).Reservoir salt concentration is c res = n res / N A = µ M. f Z ne t l B / a Z = 100 Z = 200 Z = 500 FIG. 3: Electrostatic coupling parameter Γ el = Z net λ B / a versus mi-crogel concentration φ , where Z net is the net microgel valence. Sys-tem parameters, colors, and linetypes are same as in Fig. 2. i.e., in the salt-dominated regime where κ ≈ κ s , the TPT andPBCM predictions for κ differ due to differing values for n s . However, these differences are not visible in the inset.At higher concentrations in the counterion-dominated regimewhere κ ≈ κ c , both methods predict practically the same κ ≈ κ c , determined by Z and φ . With increasing backbonevalence, κ c increases while κ s decreases, owing to increasedsalt expulsion. The total screening constant κ increases mono-tonically with increasing concentration, more steeply so forhigher Z .Figure 3 shows the electrostatic coupling strength Γ el asa function of φ . Notice that Γ el depends on the equilib-rium radius a and net microgel valence Z net , both of whichare monotonically decreasing with increasing φ . The de-crease of Z net due to inside-permeated counterions is morepronounced than the decrease of a with increasing concentra-3tion, which explains the monotonic decrease of Γ el . The over-all behavior of the coupling strength as function of concentra-tion and backbone valence is similar in the TPT and PBCM,but the TPT consistently predicts a stronger electrostatic cou-pling than PBCM. The greatest differences are visible for lowconcentrations and for the highest considered backbone va-lence Z =
500 where Γ el >
5, such that nonlinear electrostaticeffects, not accounted for in the linear TPT, come into play[38, 52, 53]. We stress here that, in contrast to suspensions ofimpermeable, solid particles, a reduction in the concentrationof permeable, compressible particles results in a strengtheningof the electrostatic interparticle repulsion. f s w e lli ng r a t i o a f f FIG. 4: Swelling ratio α = a / a versus reduced concentration φ for backbone valence Z =
500 (dashed-dotted), 200 (solid), and 100(dotted) at c res = µ M. Inset: Swollen microgel volume fraction φ = φ α versus φ . The straight dashed line in the inset depicts φ = φ α (cid:0) φ ∗ (cid:1) for a fixed TPT microgel radius taken at φ ∗ = . × − and backbone valence Z = Figure 4 depicts the concentration dependence of the equi-librium microgel swelling ratio, α = a / a , for three differ-ent backbone valences. At a given φ , the swelling ratio in-creases with increasing Z , owing to an enhanced electrostaticrepulsion between the Z monovalently charged backbone sitesfor a constant reservoir salt concentration c res = µ M.Deswelling in the counterion-dominated regime displayed inthe figure is most pronounced at smaller φ , and a decreaseshere more strongly for higher Z . For Z = α with increasing φ .An important quantity characterizing the swollen microgelsis the volume fraction φ = φ α , whose concentration depen-dence is shown in the inset for Z = φ increases sublinearly with increasing φ . Differences betweenthe TPT and PBCM predictions for φ are small except forsmall concentrations, where nonlinear electrostatic couplingis significant.To assess quantitatively the effect of deswelling on struc-tural and dynamic properties, it is useful to compare findingsfor the actual suspension of deswelling microgels with thosefor a fictitious reference suspension of non-swelling particles. We select the system parameters of the reference system to bethe same as in the actual one, except for the microgel radius a ref , which is fixed to the equilibrium value of the deswellingsystem at the lowest considered concentration, φ ref0 , wherenonlinear screening by the microions can still be disregarded.Explicitly, we set a ref = a ( φ ref0 ) using φ ref0 = . c res = µ M, and backbone valencesrestricted to values Z ≤ α predicted by the twomethods, compared with the respective constant value α ( φ = . ) (dashed horizontal lines) for the reference system.Note that the reference system microgel radius is differentfor the two methods, namely, a ref ≈ . a ref ≈ . Zn = n s andhence κ c = κ s . At very small concentrations where κ c < κ s , α changes only little with concentration. f s w e lli ng r a t i o a desw., PBCMref., PBCMdesw., TPTref., TPT 0.001 0.01 0.05 f s w e lli ng r a t i o a nZ = 2 n s FIG. 5: Predictions of TPT and PBCM for swelling ratio α = a / a versus microgel concentration φ compared with corresponding ref-erence system fixed value (dashed lines). Inset: Swelling ratio α forlow concentrations where salt-dominated regime is resolved. Systemparameters: Z =
200 and c res = µ M. It was shown in [36, 52–56] that nonlinear electrostatic cou-pling, which comes into play for Γ el (cid:38)
5, can be incorpo-rated into linear Yukawa-type effective pair potentials, suchas in Eq. (16), by using renormalized values of the particle(backbone) charge and of the Debye screening constant. Dif-ferent renormalization schemes were developed for this pur-pose for charge-stabilized suspensions of impermeable parti-cles [36, 52–56], but considerably less applications to ionicmicrogels were reported so far [5, 8, 51, 57, 58].In the remainder of this paper, we study fluid-phase suspen-sions mostly for conditions Γ el < Z net , this condition lim-its us to concentrations φ > φ ref0 = .
005 in the counterion-dominated regime where deswelling is most pronounced.4
B. Potential Parameters and Pressure Contributions
Having introduced the reference system of constant-sizemicrogels, we analyze next the parameters characterizingthe effective pair potential of deswelling microgels, in com-parison with the reference system values. For the consid-ered system parameters, the likelihood of particle overlap issmall. The microgel interaction is thus determined by thenon-overlapping (Yukawa) part of the effective pair potential, u Y ( r ; n ) , in Eq. (16). The Yukawa potential, which is charac-terized by Z net and κ , can be expressed in the form β u Y ( r ; n ) = a A Y e − κ r r , (52)where A Y = β u Y ( a ; n ) exp ( κ a ) is an interaction strengthparameter. f Z ne t desw., PBCMref., PBCMdesw., TPTref., TPT f k a desw., PBCMref., PBCMdesw., TPTref., TPT FIG. 6: Concentration dependence of net microgel valence Z net . In-set: Reduced Debye screening constant versus φ . Predictions ofTPT and PBCM for deswelling systems (solid red and black lines, re-spectively) are compared with constant-size reference system predic-tions (dashed lines). System parameters: Z =
200 and c res = µ M. In Fig. 6, the concentration dependence of the net microgelvalence Z net and the Debye screening constant κ of deswellingparticles are compared with the reference system predictions.While Z net decreases with increasing concentration, κ mono-tonically increases. This trend can be attributed to the asso-ciated increase in the number of counterions inside the mi-crogels. Deswelling slightly increases Z net , but has almost noeffect on κ , which in the counterion-dominated regime is de-termined solely by Z and φ , independent of n s [see Eq. (19)].At low φ , the Z net curves merge with those of the referencesystem, since a ref becomes at φ = .
005 equal to the radius a of the deswelling microgels. Deswelling enlarges the volumeavailable to the microions outside the microgels by a factor V ( φ − φ ref ) , where φ ref = φ ( a ref / a ) is the volume fractionof the reference system. The resulting gain in entropy forcounterions leaving the deswelling microgels is nearly com-pensated by a greater work required to expel these ions, asthe backbone charge density of opposite sign is increased bya factor ( a ref / a ) . The net effect is an only slightly increased f A Y [ x ] desw., PBCMref., PBCMdesw., TPTref., TPT f b u Y ( a ) (a) (b) FIG. 7: Influence of deswelling on interaction strength parameters (a) A Y = β u Y ( a ; n ) exp ( κ a ) and (b) β u Y ( a ; n ) of effective Yukawapair potential for non-overlapping microgels, as predicted by TPT.System parameters: Z =
200 and c res = µ M. Z net for the deswelling microgel system. Both TPT and PBCMpredict such a slight enhancement of Z net at higher φ , but withconsistently higher values for TPT.According to Figs. 7(a) and (b), A Y grows with increasingconcentration, while β u Y ( a ) decreases. The order relation Z net ( φ ) ≥ Z refnet ( φ ) is valid, which implies the order relation β u Y ( a ) ≥ β u refY ( a ) for the effective potential at contact dis-tance 2 a . The opposite order A Y ≤ A refY holds for the inter-action parameter A Y in Eq. (52). To understand these rela-tions, recall with Eq. (17) that A Y is proportional, in addi-tion to Z , to a geometric factor depending on κ a , and thisfactor is higher for the reference system. Fig. 7(b) quanti-fies the aforementioned peculiarity of ionic microgel systemsthat, with decreasing concentration, the electrostatic couplingstrength measured at contact distance is increased. f p a b ( p m - p r e s ) Z = 100 Z = 200 Z = 500 FIG. 8: PBCM prediction for microion osmotic pressure p µ − p res (in reduced units) versus φ for backbone valences Z as indicated andreservoir pressure p res = n res k B T . For Z = p µ ,calculated using the contact theorem in Eq. (37). As expected,for given φ , p µ grows rapidly with increasing backbone va-lence. For Z = φ ref > φ , of thereference system.It is instructive to compare the PB cell model pressure p µ with the total suspension pressure p from TPT and the vol-ume energy-derived contribution p V . Such a comparison isshown in Fig. 9 for a system with Z =
200 and c res = µ M.All pressures are measured relative to the reservoir osmoticpressure p res . Here, p µ is calculated according to Eq. (37)using the PBCM microion concentrations n ± ( R ) at the celledge, p according to Eq. (30) with TPT input for f ( a , n ) [Eq. (29)], and p V according to the linear-response expres-sion in Eq. (21). Also shown is the kinetic microion pressure p kin = ( Zn + n s ) k B T , with n s calculated from TPT. In thedilute limit ( φ → p res , andthe system salt concentration n s tends to n res . The dry volumefraction at which Zn = n s is φ ≈ . p to p V , inter-microgel correlationcontributions to p are significant for φ (cid:38) .
04 where p be-comes distinctly higher than p V . This comparison shows fur-ther that the pressure contribution p se , generated by the con-centration dependence of a ( n ) in the single-particle energies u e ( a ) and f p ( a ) in Eq. (27), is negligible at lower concentra-tions. The PBCM pressure p µ exceeds p V for non-zero con-centrations and is overall close to p , except at high φ . At thisrelatively low salt concentration, the kinetic microion pres-sure difference p kin − p res (dotted curve) is practically equalto the reduced ideal gas pressure of counterions, Znk B T (or3 Z φ in reduced units), up to a small negative correction pro-portional to 2 ( n s − n res ) , owing to the salt expulsion (Donnan)effect (cf. Fig. 2). While in the concentration range of Fig. 9the counterions contribute most strongly to the suspension os-motic pressure, due to the electrostatic attraction of the fixedbackbone charge they behave distinctly non-ideal, which is re-flected in the non-constant, radially decaying counterion con-centration profile n + ( r ) (see Fig. 1). This is why p kin is higherthan p . C. Structural Properties and Charge Renormalization
As explained in Sec. IV, using u eff ( r ; n ) with associated val-ues for the equilibrium radius a , net valence Z net , and screen-ing constant κ , one can calculate the microgel radial dis-tribution function (rdf) g ( r ) and static structure factor S ( q ) characterizing pair correlations in real and Fourier space,respectively. For this purpose, we use the Rogers-Young(RY) integral-equation scheme. This thermodynamically self-consistent scheme is known, from comparisons with computersimulation results, to be very accurate for fluids of charge-stabilized particles interacting via a repulsive Yukawa-type f p a b ( p - p r e s ) p (TPT) p V (TPT) p kin (TPT) p m (PBCM) FIG. 9: Reduced pressure p of microgel suspension from TPT[Eqs. (29) and (30)], volume energy contribution p V [Eq. (21)], andPBCM pressure p µ [Eq. (37)] versus microgel concentration φ . Allpressures are relative to reservoir pressure p res . System parametersare Z =
200 and c res = µ M. Also shown is the TPT prediction forthe kinetic (ideal gas) pressure, p kin = Zn + n s , where n s is systemsalt density (nearly identical to PBCM prediction). potential below concentrations where the suspension crystal-lizes [40, 59, 60].To illustrate the high accuracy of the RY method for mi-crogel particles interacting via the pair potential in Eq. (15),in Fig. 10 the RY results for g ( r ) and S ( q ) are compared withMonte-Carlo (MC) simulation data obtained using the methodin [10] for a salt-free suspension with Z = Z = r / s g ( r ) MC sim.RYHNC 0 5 10 15 q s S ( q ) FIG. 10: Results of RY, HNC, and MC for radial distribution func-tion, g ( r ) , and static structure factor, S ( q ) (inset), of a salt-free sus-pension with φ = . Z = α = . κ a = . α is computed using TPT. Length unit is microgeldiameter σ = a . S ( q m ) and the osmotic com-pressibility factor S ( ) (inset) are displayed in Fig. 11. Theequilibrium radius a in the potential u eff ( r ; n ) , on which theRY calculations are based, is determined here using the TPTand PBCM. The contact value, g ( σ ) , of the associated rdf re-mains small in the considered concentration range, showingthat the no-overlap potential, u Y ( r ) , essentially determinesthe microstructure of the microgels. The overlap potential, u ov ( r ) + u H ( r ) in Eq. (15), comes into play only at high con-centrations. From comparison with the reference system peakheight predictions, one notices that S ( q m ) is reduced whendeswelling is accounted for, though only slightly, since for Z =
200 the decrease in the microgel radius relative to the ref-erence value remains small even at higher concentrations (cf.Fig. 5).For a given concentration, the TPT predicts a more struc-tured system than the PBCM, as reflected by the higher val-ues of S ( q m ) . This difference originates from the higher netcharge and larger microgel radius in the TPT, as discussed al-ready in relation to Figs. 3 and 4. f S ( q m ) desw., PBCMref., PBCMdesw., TPTref., TPT f S ( ) x FIG. 11: RY structure factor peak height, S ( q m ) , and osmotic com-pressibility factor, S ( q → ) (inset), versus φ for Z =
200 and c res = µ M. Results are presented for deswelling microgels withradius a computed in TPT and PBCM and compared with corre-sponding results for reference system (dashed curves). Other param-eters as in Fig. 5. In discussing the Kirkwood-Buff relation [Eq. (9)], wenoted that, for a monodisperse suspension in osmotic equi-librium with a salt reservoir, S ( ) = S ( q → ) equals the os-motic compressibility factor. Rogers-Young predictions forthe concentration dependence of S ( ) (inset of Fig. 11) showthat deswelling slightly increases the osmotic compressibility.The increase of S ( ) predicted by both methods follows fromthe reduced volume fraction of deswelling particles, which islower by the factor ( a / a ref ) than that of the reference system.In its effect on S ( ) , this reduction in volume fraction over-compensates the small deswelling-induced increase of Z net (see Fig. 6). The PBCM yields distinctly higher compressibil-ities than the TPT, since it predicts smaller equilibrium radiiand net charges. The peak value of the structure factor, S ( q m ) , can be used asan indicator for the proximity of a fluid suspension to a freez-ing transition. The frequently cited empirical Hansen-Verletcriterion, S ( q m ) = .
85, applies only to the freezing of a hard-sphere fluid. It does not apply to suspensions with longer-range, soft inter-particle repulsion. As shown in detail in[60], a somewhat higher freezing indicator value, S ( q m ) = . κ n − / (cid:46) g ( r m ) ≈ . r m [60].To illustrate how the freezing transition concentration isdetermined using the citerion S ( q m ) = .
1, we consider astrongly charged microgel suspension with Z =
500 and c res = µ M, for which g ( a ) ≈ et al. [8] in linearizing the PB equation [Eq. (32)] aroundthe nonlinear potential value Φ R = Φ ( R ) at the cell boundary.This procedure leads to a linear PB equation, ∆Φ l ( r ) = κ [ Φ l ( r ) − Φ R + γ R ] + λ B Z ren a Θ ( a − r ) , (53)for r ≤ R , where κ = κ cosh ( Φ R ) , and γ R = tanh ( Φ R ) ,with the latter quantity being negative due to the negativebackbone charge. Here, Z ren is the yet unknown renormalizedbackbone valence, and Θ ( r ) is the unit step function.The unique solution for the linearized potential, Φ l ( r ) , in-side and outside the microgel sphere can be obtained analyti-cally using the boundary conditions Φ l ( R ) = Φ R , Φ (cid:48) l ( R ) = Φ (cid:48) l ( ) =
0, in conjunction with the continuity of Φ l ( r ) andits first derivative at r = a . These five conditions determine Z ren , together with the four integration constants arising fromthe integration of the linearized PB equation [Eq. (53)] insideand outside a microgel sphere. Input parameters are here Φ R and the equilibrium radius a , determined independently. Notethat the linearized potential Φ l ( r ) gives rise to the same poten-tial and electric field values at the cell boundary as the nonlin-ear PBCM potential. We refrain from quoting the somewhatlengthy analytic expressions for Φ l ( r ) and Z ren given in [8].As discussed in [8] (see also [61]), due to the monotonic in-crease of the microgel radius with increasing bare backbonevalence Z , the renormalized valence Z ren ≤ Z does not reach asaturation value beyond the linear regime, as it does for non-permeable rigid colloids [62, 63]. Instead, Z ren grows mono-tonically with increasing Z , showing only a slight indicationof a plateau behavior in the regime of intermediately high Z and low suspension salinity.In Donnan equilibrium, the renormalized net microgel7charge number, Z rennet , is obtained in the PBCM as Z rennet = − a Φ (cid:48) l ( a ) λ B = tanh ( Φ R ) κ eff λ B (cid:104) κ eff ( a − R ) cosh ( κ eff ( a − R ))+ ( κ aR − ) sinh ( κ eff ( a − R )) (cid:105) , (54)which follows alternatively from Eq. (18), wherein Z and n ± ( r ) on the right-hand side are replaced, respectively, by Z ren and the linearized microion profiles n l , ± ( r ) = n res e ∓ Φ R [ ∓ ( Φ l ( r ) − Φ R )] , (55)whose values at the cell boundary match the non-linearizedones.To implicitly account for nonlinear effects, we use Z rennet given by Eq. (54) as the input for the net valence in the linear-response, no-overlap Yukawa potential u Y ( r ) in Eq. (16). Inaddition, κ eff might be identified as the renormalized input forthe screening constant in u Y ( r ) , as done based on the origi-nal description by Alexander et al. [62] for the case of non-permeable rigid spheres [56, 63]. However, to recover thePBCM screening constant in Eq. (39) in the limiting case oflow backbone charges, where charge renormalization is notoperative, we determine the renormalized screening constant, κ ren , to be substituted into u eff ( r ) , in a manner that maintainsthe smoothness of the effective potential at r = a for the non-linear case. Explicitly, we determine κ ren as ( κ ren ) = πλ B (cid:0) nZ renapp + n rens (cid:1) , (56)where n rens = π V R (cid:90) R n l,- ( r ) r dr . (57)The apparent renormalized backbone valence, Z renapp , is definedby Eq. (40) using the substitutions Z ∗ → Z renapp , κ ∗ → κ ren , and Z ∗ net → Z rennet . Notice that κ ren is given here only implicitly, sothat an iteration procedure with starting seed κ eff is used forits calculation. For κ eff , we obtain the expression ( κ eff ) = πλ B − γ R (cid:20) nZ ren + n rens + γ R (cid:21) , (58)identical to the one for non-permeable charged colloids [63].The screening constants κ eff , κ ren , and κ mutually differ, ex-cept in the limit Z →
0, where γ R → { n s , n rens } → n res , inwhich case all three quantities then equal the reservoir screen-ing constant, κ res .Rogers-Young results for the concentration dependence of S ( q m ) of strongly charged microgels are displayed in Fig. 12.The effective pair potential parameters are determined hereusing the PBCM charge-renormalization method describedabove. The inset shows the size ratio a ( φ ) / a ref with a ref = a ( φ = . ) , as predicted by the nonlinear PBCM method.The RY values for g ( a ) are practically zero (i.e. g ( a ) < . φ , so that S ( q m ) = . f S ( q m ) desw., PBCMref., TPT f a » » FIG. 12: RY peak height, S ( q m ) , for strongly repelling microgelswith Z =
500 and c res = µ M. Results are presented for deswelling(solid curve) and reference system constant-size microgels (dashedcurve), based on the charge-renormalized PBCM. The dotted hori-zontal line marks the freezing criterion S ( q m ) = .
1. Inset: Swellingratio α = a / a in charge-renormalized PBCM. strongly charged microgels, deswelling significantly increasesthe freezing transition concentration by about 16 %, corre-sponding to a 4 % decrease in the swelling ratio α (see inset).For deswelling microgels, the freezing transition concentra-tion determined by S ( q m ) = . φ ≈ .
01. The RY rdf peakheight is here g ( r m ) = .
7, which is close to the freezing tran-sition value 2 . κ coll n − / (cid:46) κ coll by κ ren , where ( κ ren ) = πλ B nZ renapp for the present counterion-dominated microgel system, we ob-tain κ coll n − / ≈ .
3, consistent with a fluid-bcc freezing tran-sition quite close to the fluid-bcc-fcc triple point [60].
D. Diffusion and Rheological Properties
We explore next dynamic properties of ionic microgelsuspensions, using the one-component model of pseudo-microgels interacting via u eff ( r ; n ) . The deswelling ratio α ( φ ) , net valence Z net ( φ ) , and Debye screening constant κ ( φ ) in this model are determined using the TPT and PBCMmethods. As explained in Subsec. V B, the employed methodsfor calculating dynamic properties depend on u eff ( r ; n ) onlyimplicitly via the radial distribution function g ( r ) and staticstructure factor S ( q ) . On taking into account that solvent per-meability effects are very small for non-overlapping ionic mi-crogels [4, 42], we identify the hydrodynamic microgel radius a H with the equilibrium radius a ( φ ) .Figure 13(a) displays our results for the positive definite hy-drodynamic function H ( q ) at concentration φ = . S ( q ) as the only input. This inputis calculated using the RY scheme, which gives somewhat8different results in the TPT and PBCM, respectively, owingto their different predictions for a and Z net . For example, at φ = . φ = .
067 in the PBCM and φ = . S ( q ) cause less pronounced dif-ferences in H ( q ) , since the latter depends on S ( q ) only in aglobal (functional) way [42, 49]. The differences in H ( q ) aregreatest at the peak, which is located at practically the samewavenumber q m as the principal peak of S ( q ) . There are pro-nounced undulations in H ( q ) due to strong HIs between themicrogels. In the absence of HIs, H ( q ) = q . The peak height H ( q m ) exceeds unity for φ = . φ , where the hard coreof the colloidal particles is masked by the strong and long-range electrostatic repulsion [40, 64]. q s H ( q ) f H ( ) desw., PBCMref., PBCMdesw., TPTref., TPT f = 0.005 f = 0.02 f = 0.005TPT referenceTPT deswellingTPT deswellingPBCM deswelling f = 0.02 (b)(a) FIG. 13: (a) BM-PA results for the hydrodynamic function H ( q ) asfunction of reduced wavenumber, q σ , for two concentrations φ asindicated. The lower one, φ ≈ . φ = . H ( q ) of deswelling particles with TPT calculated radius is com-pared with that of the reference system. (b) Concentration depen-dence of sedimentation coefficient K = H ( q → ) . Solid curves:deswelling particles in the TPT (red) and PBCM (black). Systemparameters: Z =
200 and c res = µ M. Also displayed in Fig. 13(a) is the hydrodynamic function(with TPT input for a ) of a more concentrated suspension at φ = .
02, which corresponds in the TPT to the volume frac- f H ( q m ) f D ( q m ) / d ( n ) desw., TPTdesw., PBCM (a) (b) FIG. 14: BM-PA results for (a) hydrodynamic function peak height, H ( q m ) , and (b) reduced cage diffusion coefficient, D ( q m ) / d ( n ) ,as functions of volume fraction φ = φ α ( n ) using TPT (solid redcurve) and PBCM (sold black curve) for α ( n ) and compared withhard-sphere results (dashed curves). System parameters: Z = c res = µ M. tion φ = .
26. The principal peak of H ( q ) is here signifi-cantly below one. The reduced short-time self-diffusion co-efficient, d s / d ( a ) = H ( q σ (cid:29) ) , is accordingly significantlylower than its value for φ = . H ( q ) ’s of deswelling and referencemicrogels (solid and dashed curves, respectively, in Figs. 13and 14) are small, and basically due to the higher volumefraction of the reference system. This is also the reason forthe slight downshift of the reference-system H ( q ) relative tothe one of the deswelling system. The microgel H ( q ) bears aqualitative similarity to the one of colloidal hard spheres (hs)at the same volume fraction φ = .
26, in particular regard-ing its peak value and location. The hydrodynamic functionof hard spheres, H hs ( q ) , is likewise characterized by a peakheight below one, and the peak is located at q m σ ≈ π . TPTbased explicit values are H ( q m ) = .
81 (0 .
65) for the peakheight, d s / d = .
52 (0 .
51) for the short-time self-diffusioncoefficient, and K = H ( q → ) = .
11 (0 .
18) for the sedimen-tation coefficient, where given in brackets are the respectivevalues for colloidal hard spheres, obtained using the analyticexpressions [42, 43] H hs ( q m ) = − φ / φ cp = − . φ d hss / d = − . φ (cid:0) + . φ − . φ (cid:1) K hs = − . φ × (cid:0) − . φ + . φ − . φ + . φ (cid:1) g hs ( r m = a + ) = − . φ ( − φ ) , (59)9which are accurate for volume fractions up to the hard-spherefreezing transition value φ = . H hs ( q m ) with increasing volume fraction, whichholds to high accuracy for the complete liquid-phase concen-tration range. In the above expression for H hs ( q m ) , φ cp = π / (cid:16) √ (cid:17) ≈ .
74 is the highest possible volume fraction, at-tained for monodisperse hard spheres in close-packed fcc andhcp crystalline structures.Eq. (59) quotes also the accurate Carnahan-Starling expres-sion for the height, g hs ( a + ) , of the principal peak of the hard-sphere rdf, located at the contact distance r m = a + . TheBM-PA values of H ( q ) at q = q m and in the q → ∞ limitare somewhat higher than the corresponding hard-sphere val-ues. There are also differences between the microgel g ( r ) andthe hard-sphere g hs ( r ) (not shown here). The microgel rdffor φ = .
02 has the peak height g ( r m ) = .
50 at pair dis-tance r m = . σ , whereas g hs ( σ + ) = .
15. The differencesfrom the hard-sphere values are due to the electrostatic repul-sion between the microgels, which is here of shorter range1 / κ = . a . The pair potential contact value, β u Y ( a ) ≈ k B T (see Fig. 7(b)), reflected in a nearly zero probability, g ( a ) < − , of finding two microgels in contact.Owing to HIs, two microparticles in contact sediment fasterthan at larger separations. This underlies the fact that the sed-imentation coefficient K for a homogeneous ionic microgelsuspension is lower than the one for hard spheres at the same φ . The monotonic decline of K = V sed / V with increasingconcentration is shown in Fig. 13(b). Owing to stronger sol-vent backflow, the sedimentation velocity, V sed ( φ ) , is lowerin a more concentrated suspension than in a less concentratedone. The maximal sedimentation velocity, V sed ( φ = ) = V ,is thus attained at infinite dilution, where K =
1. Since themajor effect of deswelling is to lower φ , K is higher fordeswelling microgels than for the constant-size reference par-ticles, which explains the slightly higher values of K in thePBCM, since a PBCM < a TPT .As seen in Fig. 14(a), the H ( q m ) of ionic microgels hasa non-monotonic volume fraction dependence. Starting froma value of one at infinite dilution, with increasing φ , H ( q m ) increases towards its maximal value ∼ . φ ≈ .
07 cor-responding to φ ≈ . φ (cid:38) .
2. This behaviorshould be contrasted with the strictly linear decrease of H ( q m ) for hard spheres (curved, dashed line on the lin-log scale).In contrast to the non-monotonic H ( q m ) , both K and d s / d (latter not shown here) decrease monotonically with increas-ing φ . Furthermore, unlike the swollen radius a , which de-creases with increasing φ , the reduced Debye screening con-stant κ a increases monotonically from κ a ≈ .
24 at φ = . κ a ≈ . φ = . H ( q m ; φ ) > H hs ( q m ; φ ) , d s ( φ ) > d hss ( φ ) , and K ( φ ) < K hs ( φ ) were pre-viously demonstrated [65, 66]. These relations hold also forionic microgels provided particle overlap is very unlikely, i.e.,provided g ( a ) ≈
0. For conditions not encountered in this paper, where over-lap of microgels is likely and their softness matters, such asfor low Z or high salt content, the expected effect on H ( q ) isa flattening of its oscillations at larger q , possibly to an ex-tent that H ( q m ) ≈ d s / d . Moreover, particle softness tends toenhance K , while H ( q m ) is lowered. This behavior of H ( q ) is observed indeed in a model system of particles interact-ing by the Hertz potential [42]. Softness effects in non-ionicand weakly charged microgel systems will be the subject of aforthcoming study.The short-time diffusion function, D ( q ) , measured in unitsof d ( n ) , is determined according to Eq. (46) by the ratio ofthe hydrodynamic factor H ( q ) and S ( q ) , the latter being inde-pendent of HIs. The principal minimum of D ( q ) is located,for repulsive interactions, at practically the same wavenum-ber q m at which S ( q ) and H ( q ) attain their respective max-ima, with S ( q m ) being in general distinctly higher than H ( q m ) .The so-called cage diffusion coefficient, D ( q m ) , quantifiesthe slow relaxation of concentration fluctuations of a wave-length 2 π / q m comparable with the diameter of the dynamiccage formed around each particle by its neighbors. For hardspheres, D ( q m ) / d decreases monotonically with increasing φ , which reflects a dynamical stiffening of the next-neighborcage. The cage diffusion coefficient of hard spheres is quan-titatively described, within 2% of accuracy up to the freezingvolume fraction, by the polynomial D hs ( q m ) d = − φ − . φ + φ , (60)according to which D hs ( q m ) follows closely a linear declinewith slope − φ ∼ .
3. At freezing,where H hs ( q m ) ≈ .
33 and S hs ( q m ) ≈ . D hs ( q m ) ≈ . × d .The cage diffusion coefficient of ionic microgels is plottedin Fig. 14(b) as a function of φ , where D ( q m ) is normalizedby the concentration-dependent single-microgel diffusion co-efficient d ( n ) , allowing direct comparison with the reducedcage diffusion coefficient of hard spheres [Eq. (60)]. Un-like H ( q m ) , the reduced cage diffusion coefficient monoton-ically decreases with increasing φ . The only remnant of thepeak in H ( q m ) is a shallow inflection point in D ( q m ) / d ( n ) at φ ≈ .
07. Owing to the electrostatic repulsion, the next-neighbor cage of microgels is more structured than that ofhard spheres at the same φ , reflected in an accordingly higherstructure factor peak and lower cage diffusion coefficient. Thedistinctly higher values of S ( q m ) in the TPT, in comparison tothe PBCM, lead to lower values of D ( q m ) / d ( n ) in the TPT,which explains the reverse order in the curves of H ( q m ) and D ( q m ) / d ( n ) in Figs. 14(a) and (b), respectively.While D ( q ) is minimal at q m , it attains its maximum at q = d c = D ( q → ) .The maximum reflects the fast relaxation of long-wavelengthconcentration fluctuations by a collective diffusive motion ofparticles. In this context, recall that H ( q ) and S ( q ) are bothminimal at q =
0. At a given φ , however, S ( ) appear-ing in the denominator of d c = d ( n ) K / S ( ) is clearly below0 f d c / d r y f nZ = 2 n s no HIs FIG. 15: BM-PA results for the reduced collective diffusion coef-ficient of deswelling microgels d c / d dry0 versus φ (solid curves) forswollen radius a calculated using TPT (red) and PBCM (black). Dot-ted curves are results without HIs where H ( ) =
1. Vertical line seg-ments mark concentrations at which κ c = κ s [cf. Eq. (51)]. Systemparameters: Z =
200 and c res = µ M. Inset: Comparison withreference system results (dashed lines) for concentrations exceedingpeak position value φ ≈ . K = H ( ) , as noticed from Figs. 11 and 13(b), with a conse-quentially high value of d c .Figure 15 displays the d c for deswelling ionic microgels,obtained using the BM-PA method with respective TPT andPBCM input for a . To uncover its genuine concentration de-pendence, d c is divided, in lieu of d ( n ) , by the concentrationindependent single particle diffusion coefficient, d dry0 , of col-lapsed microgels, implying that d c / d dry0 → a / a for φ → d c is observed, with a pronounced maximum of d c at φ ≈ . H ( q m ) isgreatest. The non-monotonic concentration dependence of d c is explained on noting first that K and S ( ) are both monoton-ically decreasing with increasing φ . At low φ , the decreaseof S ( ) with increasing φ is stronger than that of K , givingrise to a growing d c . At higher concentrations, the slowinginfluence of HIs on K is strong enough that the increase of d c is turned into a monotonic decline. To show explicitly that themaximum of d c , and its decline at higher φ , are due to HIs,results for d c without HIs are included in the figure for com-parison. Without HIs, K = d c ( φ ) without HIs are monotonicallyincreasing, and they converge to the ones with HIs at very lowconcentrations only. It is further noticed that the higher val-ues of d c in the TPT are due to the lower osmotic compress-ibility values predicted by this method, and this even though d ( n ) ∝ / a in the TPT is lower than in the PBCM. Quiteinterestingly, the concentration in Fig. 15 where the numberof backbone-released counterions equals the number of saltcounterions marks an inflection point, where the shape of thecurve of d c ( φ ) changes from convex to concave.The influence of deswelling on d c at higher φ is assessed in the inset of Fig. 15, in comparison with the reference sys-tem predictions (dashed lines). Deswelling slightly enhancescollective diffusion, as predicted by both the TPT and thePBCM. This enhancement can be attributed to weaker HIs be-tween deswelling microgels, with a corresponding increase in K overcompensating the increase in S ( ) .A short discussion is in order regarding the BM-PA schemeresults for H ( q ) at wavenumbers q (cid:28) q m , where its accuracyis known to worsen with increasing concentration, up to a de-gree where non-physical negative values for K are predicted[49, 65]. This is mainly due to the self-diffusion contribu-tion to H ( q ) = H d ( q ) + d s / d , which in the hybrid schemeis calculated using the pairwise additivity (PA) approxima-tion. The PA method fully accounts for two-body HIs but ne-glects three-body and higher-order contributions. These com-plicated higher-order contributions account for the reductionin the strength of the HIs between two particles, due to a hy-drodynamic shielding by intervening particles. The disregardof this hydrodynamic shielding effect by the PA scheme leadsat higher concentrations to an underestimation of d s . The lat-ter contributes to H ( q ) most significantly at q = H d ( ) , is negative. For this reason, we showBM-PA results for K = H ( ) and d c ∝ K for concentrationsup to φ = .
02 only where the small- q BM-PA predictionsare trustworthy. f h¥ / h desw., PBCMref., PBCMdesw., TPTref., TPT FIG. 16: Modified BM theory results for the reduced high-frequencyviscosity, η ∞ / η , as function of φ , for interaction parameters and a calculated using TPT and PBCM. Solid curves are for deswellingmicrogels, while dashed curves are for the reference system. Dottedcurve is the prediction by Eq. (61) using φ = φ α ( φ ) , with α ( φ ) calculated in the PBCM. System parameters: Z =
200 and c res = µ M. Having discussed (short-time) diffusion properties of ionicmicrogel suspensions, we finally consider rheological prop-erties, namely the high-frequency (short-time) viscosity η ∞ and the zero-frequency viscosity η introduced in Eqs. (47)and (48), respectively. Our analysis is limited here to weaklysheared suspensions, where nonlinear phenomena such asshear thinning and the buildup of normal stress differencesare negligible. Just as for the diffusion properties, we iden-tify the hydrodynamic particle radius with a . As described in1Subsec. V B, η ∞ is calculated using the modified Beenakker-Mazur (BM) expression in Eq. (49). The shear stress re-laxation contribution ∆ η in η = η ∞ + ∆ η is calculated us-ing the simplified mode-coupling theory (MCT) expressionin Eq. (50). The only input to these methods is S ( q ) , whichis calculated in RY approximation based on u eff ( r ; n ) , with a obtained in the TPT and PBCM, respectively. HIs are incor-porated into the simplified MCT expression via H ( q ) , deter-mined using the BM-PA method.Figure 16 presents results for η ∞ (in units of the solvent vis-cosity η ) as a function of φ . With increasing concentration, η ∞ grows gradually, to a value at φ = .
03 only three timeshigher than the solvent viscosity. Such a modest growth withincreasing concentration is a characteristic feature of η ∞ . Inaddition, η ∞ is known to be rather insensitive to the form ofthe pair potential [39, 49] and hence to changes in the equilib-rium radius a , as reflected in the nearly coincident curves for η ∞ with a obtained from the TPT and PBCM methods. At agiven concentration, the reference microgel suspension has ahigher volume fraction than the deswelling microgels system,which explains the mildly higher viscosity values.For comparison, we show the prediction for η ∞ from thepolynomial expression, η ∞ η ≈ + φ ( + φ ) + . φ , (61)derived in [39]. This expression is a good viscosity approx-imation for dilute suspensions of strongly repelling charge-stabilized spheres with prevailing two-body HIs and low val-ues of S ( ) . As shown in Fig. 16, Eq. (61) is in qualitative ac-cord with the modified BM results, but underestimates η ∞ athigher concentrations. Note that Eq. (61), though not a virialexpansion to third order in φ , reduces to the linear Einsteinviscosity formula for very low volume fractions where the par-ticles are uncorrelated, and thus ∆ η =
0. For the hypotheticalcase of vanishing HIs, the particles remain uncorrelated onshort time scales for all fluid-phase volume fractions. In thiscase, η ∞ / η = + [ η ] φ holds for all φ , with [ η ] = / η / η , of deswellingmicrogels is plotted in Fig. 17 (solid lines). The pronouncedincrease of η at higher φ is mainly due to the shear stressrelaxation part ∆ η . The latter is more sensitive to changes inthe pair potential than η ∞ , as reflected in visibly higher val-ues of η , for φ (cid:38) .
03, when the TPT radius input is used.The higher volume fractions of the reference system in com-parison to the system of deswelling microgels imply a lowerzero-frequency viscosity for the deswelling particles, visiblein the inset at higher concentrations.
VII. CONCLUSIONS
We have presented a comprehensive theoretical study ofthe influence of concentration on deswelling, thermodynamic,structural, and dynamic properties of suspensions of weaklycross-linked, ionic microgels dispersed in a good solvent and f h / h f h / h desw., PBCMref., PBCMdesw., TPTref., TPT FIG. 17: Reduced zero-frequency viscosity, η / η , versus φ for sys-tem parameters Z =
200 and c res = µ M. The viscosity contribu-tion η ∞ is calculated using modified BM theory, and the shear stressrelaxation contribution ∆ η using simplified MCT. Inset: Comparisonwith reference system viscosity (dashed curves). in osmotic equilibrium with an electrolyte reservoir. To cal-culate microion density profiles, single-particle and bulk os-motic pressures, and state-dependent, equilibrium microgelswelling ratios, we implemented two mean-field methods andassessed their respective pros and cons. We consistently com-bined these methods – a thermodynamic perturbation theoryand a Poisson-Boltzmann spherical cell model – with calcu-lations of the net microion valence Z net and Debye screen-ing constant κ , characterizing the electrostatic part of the ef-fective one-component microgel pair potential derived fromlinear-response theory. On the basis of the effective one-component model of microion-dressed microgels, we deter-mined static structural properties, including S ( q ) and g ( r ) , byMonte-Carlo simulation and the self-consistent Rogers-Youngintegral-equation method, and used these properties as inputto the calculation of dynamic suspension properties, with thesalient hydrodynamic interactions included.At salt concentrations high enough that salt ions contributesignificantly to electrostatic screening, the microion distribu-tion inside and outside the microgels is relatively uniform andcounterion-induced deswelling is consequently weak. There-fore, our study focused on the counterion-dominated regime,with salt and microgel concentrations low enough, and micro-gel valences high enough, that deswelling is pronounced evenwithout significant particle overlap.The TPT method neglects nonlinear electrostatic effects,but accounts for inter-microgel correlations. In contrast,the PBCM method accounts for nonlinear screening by mo-bile microions, but neglects inter-microgel correlations, ex-cept for the remnant concentration dependence of the cell ra-dius. Unlike impermeable surface-charged colloidal particles,ionic microgels are characterized by electrostatic interactionswhose strength, as measured at mutual contact, increases withdecreasing microgel concentration. This property restricts theapplicability of the TPT method to non-vanishing microgelconcentrations.2While both methods predict the same trends for the effec-tive microgel pair potential, there are quantitative differencesin the swelling ratio, net valence, and the potential value atcontact, whose values are in general higher in the TPT thanin the PBCM. In the counterion-dominated regime, the range1 / κ of the electrostatic repulsion is equal in both methods.The greatest differences in the pair potential parameters occurat very low concentrations and high backbone valences, whichcan be partially attributed to the linear-response approxima-tion inherent in the TPT. The relative variation in the microgelradius with changing concentration is less pronounced in thePBCM, in which nonlinear response confines the counterionsmore strongly to the microgel interior.Differences in predictions of the TPT and PBCM methodsare more pronounced for static (thermodynamic and struc-tural) properties than for dynamic properties, which can be ex-plained by the fact that dynamic properties depend only glob-ally (i.e., functionally) on S ( q ) . The only exception is the col-lective diffusion coefficient d c which is directly proportionalto the inverse of the static compressibility factor 1 / S ( ) .Owing to the dominance of the electrostatic interactions inthe considered microgel systems, their dynamic behavior re-sembles that of charged-stabilized suspensions of imperme-able solid particles. In particular, the peak, H ( q m ) , of the hy-drodynamic function has a non-monotonic concentration de-pendence, with a maximum higher than one at an intermedi-ate concentration value, reflected in a concomitant inflectionpoint of the cage diffusion coefficient. The collective diffusioncoefficient, d c , behaves likewise non-monotonically and hasits maximum at the same concentration as H ( q m ) . This maxi-mum was shown to arise from the slowing effect of HIs, whichbecomes more influential with increasing concentration. Theelectric repulsion between the microgels distinctly enhancesthe zero-frequency viscosity at higher concentrations, as com-pared to suspensions of uncharged particles.The comparison with corresponding results for the refer-ence system of constant-sized microgels revealed that the ma-jor influence of deswelling on structural and dynamic proper-ties is via the reduced volume fraction φ , which grows onlysublinearly with increasing concentration n ∝ φ . The effectof counterion-induced deswelling on structural and dynamicproperties is overall quite weak, for valences where nonlin-ear electrostatic contributions to the microion distributions arenegligible, and changes of α = a ( φ ) / a with concentrationare accordingly small.At higher concentrations, deswelling slightly enhances S ( ) and the hydrodynamic function H ( q ) for all wavenumbers q .Deswelling reduces the zero- and high-frequency viscosityand slightly enhances collective diffusion. From the behavior of d c = d H ( ) / S ( ) , one notices that the deswelling-inducedenhancement of d ∝ / a is nearly counterbalanced by the ac-companying de-enhancement of H ( ) / S ( ) .The most pronounced effect of deswelling is to shift thefreezing (crystallization) transition to higher concentrationvalues, as we have determined from an empirical freezingrule for the static structure factor peak height. This concen-tration shift is more pronounced for strongly charged micro-gels, in which case the nonlinear PBCM method can still beused to determine the swelling ratio α . To determine the con-centration shift, however, the PBCM must be combined witha charge renormalization procedure to determine renormal-ized values of the microgel net valence and screening con-stant from a linearized Poisson-Boltzmann equation in the cellmodel. The renormalized parameters are used in the linear-response pair potential [Eq. (16)], where they summarily ac-count for the enhanced accumulation of counterions insideand close to the spherical backbone region. We illustratedsuch a renormalization procedure using a linearization of thenonlinear PB equation around the potential at the cell bound-ary. While such a linearization is most commonly used inrenormalization schemes applied to non-permeable and per-meable colloidal particles, it is not the only choice. Thereare sound reasons to use instead a linearization around themean (i.e., cell-volume-averaged) electrostatic potential value[63, 67]. A proper assessment of the pros and cons of dif-ferent charge renormalization schemes, formulated also in theframework of TPT, was not in the scope of the present work,but is the subject of a forthcoming paper [38].Finally, in the presented generic study, we considered onlyuniformly cross-linked microgels, modeled by a uniformlydistributed backbone charge. In future work, extensions tononuniformly charged microgels can be explored where, e.g.,the backbone charge is concentrated near the particle periph-ery. Such a nonuniform charge distribution is expected to sig-nificantly affect deswelling and the strength of the effectivepair potential, and consequently also structural and dynamicproperties. Acknowledgments
M.B. and G.N. thank J. Riest (Viega Holding GmbH & Co.KG, Attendorn, Germany) and G.-W. Park (FZ Jülich, Ger-many) for many helpful discussions. This work was underappropriation of funds from the Deutsche Forschungsgemein-schaft (SFB 985, project B6). [1] Fernández-Nieves, A., Wyss, H., Mattsson, J., and Weitz, D. A.,editors,
Microgel Suspensions: Fundamentals and Applica-tions , Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim,2011.[2] Lyon, L. A. and Fernández-Nieves, A., Annu. Rev. Phys. Chem. (2012) 2.1. [3] Plamper, F. A. and Richtering, W., Acc. Chem. Res. (2017)131.[4] Holmqvist, P., Mohanty, P. S., Nägele, G., Schurtenberger, P.,and Heinen, M., Physical Review Letters (2012) 1.[5] Nöjd, S. et al., Soft Matter (2018) 4150.[6] Gasser, U., Scotti, A., and Fernandez-Nieves, A., Phys. Rev. E (2019) 042602.[7] Romeo, G., Imperiali, L., Kim, J.-W., Fernández-Nieves, A.,and Weitz, D. A., J. Chem. Phys. (2012) 124905.[8] Colla, T., Likos, C. N., and Levin, Y., J. Chem. Phys. (2014) 234902.[9] Denton, A. R. and Tang, Q., J. Chem. Phys. (2016) 164901.[10] Weyer, T. J. and Denton, A. R., Soft Matter (2018) 4530.[11] Hofzumahaus, C., Hebbeker, P., and Schneider, S., Soft Matter (2018) 4087.[12] Urich, M. and Denton, A. R., Soft Matter (2016) 9086.[13] Mohanty, P. S. et al., Scientific Reports (2017) 1487.[14] Nir, O., Trieu, T., Bannwarth, S., and Wessling, M., Soft Matter (2016) 6512.[15] Roa, R., Zholkovskiy, E. K., and Nägele, G., Soft Matter (2015) 4106.[16] Roa, R. et al., Soft Matter (2016) 4638.[17] Park, G. W. and Nägele, G., to be submitted (2019).[18] Denton, A. R., Phys. Rev. E (2003) 011804.[19] Hedrick, M. M., Chung, J. K., and Denton, A. R., J. Chem.Phys. (2015) 034904.[20] Flory, P. J. and Rehner, J., J. Chem. Phys. (1943) 512.[21] Flory, P. J. and Rehner, J., J. Chem. Phys. (1943) 521.[22] Flory, P. J., Principles of Polymer Chemistry , Cornell Univer-sity Press, 1953.[23] Gottwald, D., Likos, C. N., Kahl, G., and Löwen, H., J. Appl.Cryst. (2005) 074903.[24] Riest, J., Mohanty, P. S., Schurtenberger, P., and Likos, C. N.,Z. Phys. Chem. (2012) 711.[25] Hansen, J.-P. and McDonald, I. R.,
Theory of Simple Liquids ,Elsevier, 2013.[26] Landau, L. D. and Lifschitz, E. M.,
Theory of Elasticity , Else-vier, 1986.[27] Rovigatti, L., Gnan, N., Ninarello, A., and Zaccarelli, E.,Macromolecules (2019) 4895.[28] Louis, A. A., J. Phys.: Condens. Matter (2002) 9187.[29] Hoffmann, N., Likos, C. N., and Löwen, H., J. Chem. Phys. (2004) 7009.[30] Kirkwood, J. G. and Buff, F. P., J. Chem. Phys. (1951) 774.[31] Dobnikar, J., Castaneda-Priego, R., Von Grünberg, H. H., andTrizac, E., New Journal of Physics (2006).[32] Nägele, G., Phys. Rep. (1996) 215.[33] Henderson, R. L., Phys. Lett. A (1974) 197.[34] Henderson, D., Condensed Matter Physics (2009) 127.[35] Denton, A. R., Phys. Rev. E (2006) 041407.[36] Colla, T. E., Levin, Y., and Trizac, E., J. Chem. Phys. (2009) 074115.[37] Wennerstrom, H., Jonsson, B., and Linse, P., J. Chem. Phys. (1982) 4665.[38] Brito, M. E., Riest, J., Denton, A. R., and Nägele, G., Criticalassessment of methods for calculating effective interactions andpressure in charge-stabilized dispersions, to be submitted.[39] Banchio, A. J. and Nägele, G., J. Chem. Phys. (2008)104903. [40] Banchio, A. J., Heinen, M., Holmqvist, P., and Nägele, G., J.Chem. Phys. (2018) 134902.[41] Nägele, G., Colloidal hydrodynamics, in Physics of Com-plex Colloids , edited by Bechinger, C., Sciortino, F., and Ziherl,P., Proceedings of the International School of Physics "EnricoFermi", 2013.[42] Riest, J., Eckert, T., Richtering, W., and Nägele, G., Soft Matter (2015) 2821.[43] Pamvouxoglou, A., Bogri, P., Nägele, G., Ohno, K., and Pe-tekidis, G., J. Chem. Phys. (2019) 024901.[44] Löwen, H., Palberg, T., and Simon, R., Phys. Rev. Lett. (1993) 1557.[45] Nägele, G., Kollmann, M., Pesché, R., and Banchio, A. J.,Molecular Physics (2002) 2921.[46] Abade, G. C., Cichocki, B., Ekiel-Jezewska, M. L., Nägele, G.,and Wajnryb, E., J. Chem. Phys. (2010) 084906.[47] Szymczak, P. and Cichocki, B., J. Stat. Mech. (2008)P01025.[48] Nägele, G. and Bergenholtz, J., J. Chem. Phys. (1998)9893.[49] Heinen, M., Banchio, A. J., and Nägele, G., J. Chem. Phys. (2011) 154504.[50] Das, S. et al., Soft Matter (2018) 92.[51] Braibanti, M., Haro-Pérez, C., Quesada-Pérez, M., Rojas-Ochoa, L. F., and Trappe, V., Phys. Rev. E (2016) 1.[52] Denton, A. R., J. Phys.: Condens. Matter (2008) 494230.[53] Denton, A. R., J. Phys.: Condens. Matter (2010) 364108.[54] Trizac, E., Bocquet, L., Aubouy, M., and Von Grünberg, H. H.,Langmuir (2003) 4027.[55] Pianegonda, S., Trizac, E., and Levin, Y., J. Chem. Phys. (2007) 014702.[56] Boon, N., Guerrero-García, G. I., van Roij, R., and Olvera dela Cruz, M., Proceedings of the National Academy of Sciences (2015) 9242.[57] Baulin, V. A. and Trizac, E., Soft Matter (2012) 6755.[58] Aguirre-Manzo, L. A. et al., Phys. Rev. E (2019) 032602.[59] Gapinski, J., Nägele, G., and Patkowski, A., J. Chem. Phys. (2012) 024507.[60] Gapinski, J., Nägele, G., and Patkowski, A., J. Chem. Phys. (2014) 124505.[61] Levin, Y., Diehl, A., Fernández-Nieves, A., and Fernández-Barbero, A., Phys. Rev. E (2002) 036143.[62] Alexander, S. et al., J. Chem. Phys. (1984) 5776.[63] Trizac, E., Bocquet, L., Aubouy, M., and Von Grünberg, H. H.,Langmuir (2003) 4027.[64] Westermeier, F. et al., J. Chem. Phys. (2012) 114504.[65] Heinen, M., Holmqvist, P., Banchio, A. J., and Nägele, G., J.Appl. Cryst. (2010) 970.[66] Gapinski, J., Patkowski, A., and Nägele, G., J. Chem. Phys. (2010) 054510.[67] Deserno, M. and von Grünberg, H.-H., Phys. Rev. E66