Matthew
Urich
and
Alan R.
Denton
*
Department of Physics, North Dakota State University, Fargo, ND 58108-6050, USA. E-mail: alan.denton@ndsu.edu
First published on 17th October 2016
Microgels are soft colloidal particles that, when dispersed in a solvent, swell and deswell in response to changes in environmental conditions, such as temperature, concentration, and pH. Using Monte Carlo simulation, we model bulk suspensions of microgels that interact via Hertzian elastic interparticle forces and can expand or contract via trial moves that allow particles to change size in accordance with the Flory–Rehner free energy of cross-linked polymer gels. We monitor the influence of particle compressibility, size fluctuations, and concentration on bulk structural and thermal properties by computing particle swelling ratios, radial distribution functions, static structure factors, osmotic pressures, and freezing densities. For microgels in the nanoscale size range, particle compressibility and associated size fluctuations suppress crystallization, shifting the freezing transition to a higher density than for the hard-sphere fluid. As densities increase beyond close packing, microgels progressively deswell, while their intrinsic size distribution grows increasingly polydisperse.
Over the past two decades, numerous experimental and modeling studies have characterized the elastic properties of single microgel particles9–24 and the equilibrium and dynamical behavior of bulk suspensions.25–36 Connections between single-particle properties, such as swelling ratio, and bulk properties, such as osmotic pressure, thermodynamic phase behavior, pair structure, and viscosity, have been probed experimentally by static and dynamic light scattering, small-angle neutron scattering, confocal microscopy, and osmometry.37–45 While the swelling/deswelling behavior of microgels has been extensively studied, the full implications of elasticity and compressibility of these soft colloids for bulk suspension properties remain only partially understood.
Previous modeling studies have examined the influence of elastic interparticle interactions on structure and phase behavior. In an extensive simulation survey, Pamiès et al.46 mapped out the thermodynamic phase diagram of a model of elastic, but incompressible, spheres interacting via a Hertz pair potential.47 For a related model of ionic microgels, one of us and coworkers48 recently applied molecular dynamics simulation and thermodynamic perturbation theory to compute the osmotic pressure and static structure factor of a suspension of Hertzian spheres that interact also electrostatically, via an effective Yukawa pair potential.49 At the single-particle level, numerous studies have established that the classic Flory–Rehner theory of swelling of cross-linked polymer networks,50 though originally developed for macroscopic gels, provides a reasonable description also of the elastic properties of microgel particles.3,4,12–22,33,34,51 However, the implications of particle compressibility and intrinsic size polydispersity for bulk properties of microgel suspensions have not been fully explored. The purpose of this paper is to analyze the combined influences of particle elasticity, compressibility, and associated size fluctuations on bulk thermal and structural properties of microgel suspensions.
The remainder of the paper is organized as follows. In Section 2, we describe a model of microgels as compressible spheres, whose swelling is governed by a single-particle free energy derived from the Flory–Rehner theory of cross-linked polymer networks, and whose interparticle interactions are represented by a Hertz elastic pair potential. In Section 3, we outline Monte Carlo simulation methods by which we modeled bulk suspensions of compressible, size-fluctuating microgels. Section 4 presents numerical results for thermodynamic properties, specifically osmotic pressure and liquid–solid phase behavior, and structural properties, including particle volume fraction, radial distribution function, and static structure factor. Finally, in Section 5, we summarize and conclude.
For the uniform-sphere model, with particle size swelling ratio α ≡ a/a0, the Flory–Rehner theory of polymer networks combines mixing entropy, polymer–solvent interactions, and elastic free energy to predict the total Helmholtz free energy of the network:50
(1) |
P0(α) ∝ exp[−βF(α)]. | (2) |
(3) |
(4) |
(5) |
The elastic properties of bulk, water-swollen hydrogels have been measured using scanning force microscopy.59 For poly(N-isopropylacrylamide) (PNIPAM) hydrogel with 0.25 mol% cross-linker, Young's modulus was determined as 0.33 kPa in the fully swollen state at 10 °C and 13.9 kPa in the collapsed state at 40 °C. In contrast, poly(acrylamide) (PA) hydrogels with cross-linker between 1 and 5 mol% have Young's moduli measured in the range from 90 to 465 kPa at 20 °C. Typical Poisson ratios for these gels are ν ≃ 0.5. In Section 4, we appeal to such measurements to select realistic pair potential parameters for our simulations.
(6) |
To guide analysis of phase behavior, we computed structural properties that probe interparticle pair correlations and thermodynamic properties that govern phase stability. First, we determined the radial distribution function g(r) by standard means, histogramming into radial bins the radial separation r of particle pairs in each configuration and averaging over configurations. Second, we computed the static structure factor S(q), proportional to the Fourier transform of g(r) and to the intensity of scattered radiation at scattered wave vector magnitude q. For a system of N particles, the orientationally averaged two-particle static structure factor is defined as
(7) |
(8) |
Fig. 2 Typical snapshot from an MC simulation of a suspension of compressible, spherical microgel particles in a cubic box with periodic boundary conditions. |
In the low-density (dilute) regime, the initial lattice structure rapidly fell apart, as particles freely diffused through the system, implying a stable fluid phase. In the high-density (concentrated) regime, strongly interacting particles remained in the fcc lattice structure, consistent with a stable or metastable crystalline phase. At intermediate densities, relative stability of fluid and solid phases depended on the interparticle interactions. We emphasize that we did not attempt to determine thermodynamic phase stability, which would require more extensive simulations and thermodynamic integration to compute total free energies and perform coexistence analyses. Nor did we investigate the glass transition,69 for which purpose molecular dynamics simulation would be better suited. Nevertheless, our results suggest that compressibility of microgels may significantly alter the phase diagram of strictly Hertzian spheres.46
To explore dependence of bulk properties on particle size and pair interactions, we studied two systems distinguished by widely different particle sizes and interactions. The first system is a suspension of relatively small microgels (nanogels) of dry radius a0 = 10 nm, comprising Nm = 2 × 105 monomers – consistent with monomers of radius 0.3 nm – and Nch = 200 chains, corresponding to 0.05 mol% of cross-linker (assuming two chains per cross-linker). Focusing on polymers in a good solvent, e.g., PNIPAM or PAA in water, we took χ = 0. For these loosely cross-linked particles, the Flory–Rehner free energy in the dilute limit [eqn (1)] attains a minimum at swollen radius as = 35.055 nm (α = 3.5055). For particles of this size with Young's modulus Y = 100 kPa and Poisson ratio ν = 0.5,59eqn (5) yields a reduced Hertz pair interaction amplitude ε* ≃ 1.5 × 103. To explore sensitivity to particle elasticity, we also considered microgels with ε* ≃ 1.5 × 104.
The second system studied contains much larger particles, of dry radius a0 = 100 nm, comprising Nm = 2 × 108 monomers, more closely matching many experimentally studied materials. To again model loosely cross-linked particles in a good solvent, we chose Nch = 2 × 105 (0.05 mol% cross-linker) and χ = 0. The Flory–Rehner free energy [eqn (1)] is now a minimum at swollen radius as = 350.55 nm (dilute limit). Young's moduli in the range Y = 50–1000 kPa correspond to reduced Hertz pair interaction amplitudes ε* ≃ 106–107. We note that the densities and temperatures (ε* values) studied here correspond to thermodynamic states that fall well within the range of fcc crystal stability in the known phase diagram of Hertzian spheres.46 Indeed, we never observed structural reassembly from fcc into bcc or any other crystalline structure.
Normalized probability distributions for the equilibrium swelling ratio are shown in Fig. 3 and 4 over a range of concentrations. Below close packing (ζ < 0.74), steric (Hertzian) interparticle interactions are sufficiently weak that the particles are not significantly compressed, as reflected by the distributions differing negligibly from the intrinsic dilute-limit distribution [eqn (2)]. With increasing ζ, as close packing is approached and exceeded, the distributions not only shift toward smaller α, as particles become compressed, but also broaden, i.e., become more polydisperse. The smaller particles (Fig. 3) have a polydispersity (ratio of standard deviation to mean) that increases from roughly 2% to 4% as ζ increases from 0 to 2. In contrast, the larger particles (Fig. 4) have negligible polydispersity when uncrowded, but fluctuate significantly in size for ζ ≥ 1. For both small and large particles, with increasing Hertz pair potential amplitude, the degree of compression and the breadth of polydispersity increase. The reason why the polydispersity distribution broadens with increasing concentration is that compression forces particles into a range of sizes in which the curvature of the Flory–Rehner free energy [eqn (1)] with respect to α is weaker than for uncompressed particles. Measurements of particle size in suspensions of microgels whose dry radii have a static polydisperse70 or bidisperse71 distribution show that swollen particles have an equilibrium polydispersity that decreases with increasing particle compression. Conversely, our study demonstrates that microgels whose dry radius is monodisperse have swollen polydispersity – associated purely with fluctuations in swelling ratio – that increases with particle compression.
Fig. 4 Normalized probability distribution P(α) of swelling ratio α of compressible microgels of dry radius a0 = 100 nm, composed of Nm = 2 × 108 monomers with Nch = 2 × 105 chains in a solvent with Flory solvency parameter χ = 0 at generalized volume fractions ζ = 1, 1.5, 2 (right to left). The particles interact via a Hertz pair potential with reduced amplitude (a) ε* = 106 and (b) ε* = 107, where the ζ = 1 distribution (off scale to right) is extremely narrow and sharply peaked. Note the change in scale from Fig. 3. |
Particle compressibility is illustrated further in Fig. 5, which plots the volume fraction ϕ vs. generalized volume fraction ζ. Below close packing, where particle compression is negligible, ϕ = ζ. With increasing density, however, interparticle interactions become stronger and so energetically costly, relative to the free energy cost of compression, that the particles more readily contract to minimize interactions. As a result, the volume fraction is lower than the generalized volume fraction (ϕ < ζ) for ζ > 0.74, the density at which nearest-neighbor particles in the fcc crystal would begin to overlap. For these systems, significant particle compression occurs only in the solid phase. Moreover, the higher the Hertz pair potential amplitude – for given values of Nm, Nch, and χ in the Flory–Rehner free energy – the more easily the particles yield to compression.
To explore the evolution of structure as a function of density and interparticle interactions, and to aid our diagnosis of the liquid–solid phase transition, we computed radial distribution functions and static structure factors, as described in Section 3. Fig. 6 shows our results for g(r) over a range of densities that span the freezing transition. The emergence of distinct peaks signals the onset of crystallization at ζ ≃ 0.58 for a0 = 10 nm and ζ ≃ 0.49 for a0 = 100 nm (black curves in Fig. 6). For both particle sizes, the freezing density is insensitive to ε* over the range of values considered. Suspensions of larger microgel particles thus crystallize at the same volume fraction as the hard-sphere fluid, while suspensions of smaller nanogel particles remain fluid up to significantly higher volume fractions. This interesting difference in phase stability can be attributed to the relatively softer Hertz pair repulsion between the smaller particles.
Fig. 6 Radial distribution functions g(r) vs. radial distance r (in units of uncrowded swollen radius as) of microgel suspensions, of generalized volume fractions ζ straddling freezing transition (black curves), with same particle parameters as in Fig. 5 and Hertz pair potential reduced amplitudes (a) ε* = 1.5 × 103, (b) ε* = 1.5 × 104, and (c) ε* = 107. Red and blue curves correspond to fluid and solid phases, respectively. |
Fig. 7 shows corresponding results for S(q), computed using eqn (7) by averaging over the same configurations used to compute g(r) in Fig. 6. With increasing density, the peaks progressively sharpen and shift toward larger q. At the densities at which distinct peaks in g(r) indicate crystallization, the main peak of S(q) attains a height in the range S(qmax) = 2.8–3. This freezing threshold is consistent with the Hansen–Verlet criterion for fluids interacting via a Lennard-Jones pair potential,72 according to which S(qmax) ≃ 2.85 at freezing. In passing, we note that a different freezing criterion proposed specifically for ultrasoft pair-potential fluids73 does not apply here, since the Hertz pair potential, although bounded, has amplitudes in our system that far exceed the thermal energy kBT. As a consistency check, we also computed the Lindemann parameter L, defined as the ratio of the root-mean-square displacement of particles from their equilibrium lattice sites to the lattice constant, and confirmed that L < 0.15 and remains constant in the solid phase.
Fig. 7 Static structure factors S(q) vs. scattered wave vector magnitude q [from eqn (7)] corresponding to radial distribution functions in Fig. 6. Freezing occurs when main peak height S(qmax) > 2.8 (black curves). Red and blue curves correspond to fluid and solid phases, respectively. |
Finally, we turn to the osmotic pressure, which we computed from our simulations using eqn (8). Fig. 8 presents plots of osmotic pressure vs. generalized volume fraction (equation of state) for the two microgel systems studied. In the dilute regime, the osmotic pressure is very close to that of the hard-sphere fluid, which is consistent with the swelling ratio results (Fig. 5). At volume fractions above 20%, deviations emerge and steadily grow with increasing volume fraction. The compressible microgels have osmotic pressures consistently lower than for the hard-sphere fluid. At volume fractions approaching close packing, the microgel osmotic pressure rapidly rises, yet remains finite beyond close packing, in sharp contrast to the hard-sphere fluid, whose pressure diverges at close-packing. Furthermore, as is clear from Fig. 8a, the osmotic pressure rises more gradually with concentration the smaller and softer the particles. This equilibrium trend may have implications for rheological properties, such as shear viscosity. In fact, recent measurements74 of the zero-shear viscosity of aqueous suspensions of monodisperse, highly branched, phytoglycogen nanoparticles75 are consistent with significant compression of the particles.
Fig. 8 Equation of state: reduced osmotic pressure βπas3vs. generalized volume fraction ζ [from eqn (8)] for suspensions of compressible microgels with χ = 0 and (a) a0 = 10 nm, as = 35.055 nm, Nm = 2 × 105, Nch = 200, and (b) a0 = 100 nm, as = 350.55 nm, Nm = 2 × 108, Nch = 2 × 105, for various values of Hertz pair potential reduced amplitude ε*. Insets display broader ranges of ζ and π. Pressure of hard-sphere fluid (ε* → ∞) is shown for reference. With increasing ε*, osmotic pressure becomes more hard-sphere-like. |
At the apparent freezing density of the microgel suspensions, as diagnosed by the main peak height of the static structure factor, the osmotic pressure plateaus, consistent with a thermodynamic phase transition between fluid and solid. Thus, the freezing behavior of the compressible microgels modeled here seem to be accurately described by the Hansen–Verlet criterion. Our methods are not suited, however, to identifying a glass transition, which could preempt crystallization.
Finally, to probe the dependence of system properties on solvent quality, we have repeated several of the calculations for χ = 0.5, modeling a theta solvent. As χ increases, the particles become progressively compressed, as would be expected with worsening solvent quality, and exert lower osmotic pressure. However, with varying concentration, the systems display similar qualitative trends in polydispersity, volume fraction, structure, and osmotic pressure.
Radial distribution functions, static structure factors, and osmotic pressures all display behavior similar to that of a hard-sphere fluid at low densities – volume fractions below about 0.3 – but reveal the particles' intrinsically soft nature at higher densities. Swelling ratios confirm that, with increasing density above close packing, the particles become progressively compressed and polydisperse, compared with their dilute (uncrowded) states. While suspensions of relatively large, stiff particles (microgels of dry radius 100 nm) crystallize at the same volume fraction as a fluid of hard spheres, internal degrees of freedom associated with particle compressibility and size polydispersity enable suspensions of smaller, softer particles (nanogels of dry radius 10 nm) to remain in a stable fluid phase up to densities significantly beyond freezing of the hard-sphere fluid. For both nanogels and microgels, the freezing transition is accurately predicted by the Hansen–Verlet criterion.
Our calculations of equilibrium swelling ratios and bulk osmotic pressures are qualitatively consistent with observations of significant compression of deformable particles, albeit ionic microgels, only at volume fractions approaching and exceeding close packing.17,24 Similarly, our finding that particle softness lowers the osmotic pressure is consistent with experimental measurements of reduced zero-shear viscosity in suspensions of soft nanoparticles.74 Furthermore, our calculations may guide the choice of system parameters in future experiments and molecular-scale simulations and could be extended to map out equilibrium phase diagrams.
The coarse-grained model of nonionic, spherical microgels studied here may be refined by incorporating a more realistic swelling free energy56 that goes beyond the limitations of the Flory–Rehner theory noted in Section 2.1, and may be further developed in several directions. For example, the model may be extended to ionic microgels, whose properties emerge from an interplay between elastic and electrostatic interparticle interactions. Such an extension would enable analyzing the coupled influences of fluctuating particle size and counterion distribution on swelling, structure, and phase behavior of concentrated suspensions of soft, charged colloids. By incorporating shape fluctuations, our simulation method also may be generalized to model dynamically deformable particles, as in recent studies of elastic colloids,60–62 and extended to study the influence of added hard crowders, as in our previous studies of polymer–nanoparticle mixtures.64–66 Finally, rheological properties, such as diffusion coefficients, viscosity, elastic moduli, jamming and glass transition densities69 could be explored by wedding our model of soft, compressible, size-fluctuating colloids to appropriate simulation methods, such as Brownian dynamics, lattice Boltzmann, or molecular dynamics of micromechanical models.61,62
This journal is © The Royal Society of Chemistry 2016 |