In Silico Synthesis of Microgel Particles.

Microgels are colloidal-scale particles individually made of cross-linked polymer networks that can swell and deswell in response to external stimuli, such as changes to temperature or pH. Despite a large amount of experimental activities on microgels, a proper theoretical description based on individual particle properties is still missing due to the complexity of the particles. To go one step further, here we propose a novel methodology to assemble realistic microgel particles in silico. We exploit the self-assembly of a binary mixture composed of tetravalent (cross-linkers) and bivalent (monomer beads) patchy particles under spherical confinement in order to produce fully bonded networks. The resulting structure is then used to generate the initial microgel configuration, which is subsequently simulated with a bead-spring model complemented by a temperature-induced hydrophobic attraction. To validate our assembly protocol, we focus on a small microgel test case and show that we can reproduce the experimental swelling curve by appropriately tuning the confining sphere radius, something that would not be possible with less sophisticated assembly methodologies, e.g., in the case of networks generated from an underlying crystal structure. We further investigate the structure (in reciprocal and real space) and the swelling curves of microgels as a function of temperature, finding that our results are well described by the widely used fuzzy sphere model. This is a first step toward a realistic modeling of microgel particles, which will pave the way for a careful assessment of their elastic properties and effective interactions.

Introduction

Soft colloids, combining properties of hard-sphere colloids and polymers, offer the interesting possibility to tailor their macroscopic behavior and flow at the molecular level.[1] Hard sphere colloids have served for decades as a reference model to shed light on many physics problems, such as the structure of atomic liquids,[2] crystal nucleation,[3] and the glass transition.[4−6] Recently, however, soft colloids have become even more popular in the scientific community.[7] Typically, soft particles have an internal polymeric architecture, which allows them to reach the paste regime,[8] where particles are in very dense, squeezed (or jammed) states. Under such conditions, the polymeric chains experience a substantial interpenetration,[9] and the internal structure of particles themselves changes. Since the single-particle elasticity alters the properties of the resulting macroscopic material, establishing the crucial link between microscopic properties and macroscopic response requires the knowledge of the effective interactions among particles. Yet, while synthesis and experiments have produced an immense database of soft polymeric colloids,[1] theoretical efforts still lag behind. The most dramatic example of such a dichotomy is offered by microgels, which are soft particles with tunable swelling properties and an extremely broad range of materials applications. Microgels are nano- or microsized particles made by cross-linked polymer networks, which can respond to external stimuli, such as changes to temperature T or pH, by swelling and deswelling (on a short time scale as compared to macrogels). The important experimental advantage of using microgels as a model system is that, thanks to the aforementioned responsivity, the size of the particles can be finely controlled. Such a sensitivity is usually exploited to carefully tune the sample volume fraction without changing the particle number density. This property allows to explore classical physics problems from new points of view: for example, by using confocal microscopy, it is possible to follow the melting of colloidal crystals by tuning the temperature, allowing to investigate fluidization events at grain boundaries[10] or to study the nature of the fluid–hexatic–crystal transition in 2D.[11] Similarly, a slow annealing induces a reversible assembly of crystals even in the presence of defects,[12] while size changes can be used to smartly vary the lattice spacing and to obtain color-tunable crystals.[13]

Among microgels, largely studied are those based on PNIPAM cross-linked networks, which are thermoresponsive and display a volume phase transition (VPT) at T ∼ 32 °C, from a swollen state at low T to a compact state at high T. The first synthesis of PNIPAM microgels in 1986[14] was followed by hundreds of publications, including several reviews on synthesis,[15,16] applications,[17,18] and physical aspects[7,19] as well as several books.[20−22] Despite microgels being one of the most experimentally studied soft matter systems, most theoretical efforts have adopted so far an approach that neglects the polymeric nature of the particles, even though at high density this aspect becomes overwhelmingly important. As underlined in a recent review by Lyon and Fernandez-Nieves,[19] the polymer/colloid duality of microgels, which is at the core of their technological and fundamental relevance, has largely been overlooked. Thus, the description of microgel phase behavior and dynamics has mostly been based on coarse-grained effective interactions[23] as simple as the Hertzian model for elastic spheres.[24−26] Only in the case of ionic microgels an effective Hamiltonian has been derived from first-principles, thanks to the fact that electrostatics is the dominant contribution.[27,28]

The purpose of this work is to build up a flexible numerical protocol able to design individual microgel particles in silico with properties comparable to the experimental ones. It is worth stressing that we do not seek to realistically reproduce the kinetics of the chemical synthesis, but rather to deliver a final product with characteristics as close as possible to the synthesized particles. In this work we will thus build the microgel particles and compare their swelling behavior with experimental results. To do so, we specifically focus on the case of very small microgels, or nanogels, whose diameter in the swollen regime is approximately ∼50 nm. This value appears to be close to the smallest possible size for which stable microgel particles can be synthesized using surfactant-based methods.[29] We choose to work with such small microgels because in this case we are able to reproduce the network in a monomer-resolved way by using the classic bead–spring model for polymers. In this respect, each bead has a size comparable to the Kuhn length of the PNIPAM chains,[30] thus avoiding complications due to coarse-graining.

The strategy we devise makes use of a binary mixture of patchy particles under spherical confinement to self-assemble fully bonded disordered networks. These networks are then topologically constrained and simulated using classical polymer models. By contrast, previous numerical efforts in microgel modeling have focused on unrealistic networks formed by chains of the same length, often connecting cross-linkers placed on crystalline lattice sites.[31−36] A notable exception is provided in refs (37 and 38), where closed polymer networks are constructed by directly joining randomly distributed cross-linkers with polymer chains.

In addition to their disordered nature, the networks we generate possess the additional advantage of having their swelling properties dependent on the extent of the confinement employed during the initial assembly. We find an excellent agreement between the measured swelling curve for small microgels and our numerical results upon rather strong confinement. At the same time, we show that the experimental behavior cannot be reproduced using a crystal-based network microgel of the same size. We further investigate the internal structure of the particles and successfully compare the form factors with the widely used fuzzy sphere model,[21] thus achieving a complete characterization of the particles. Our study is thus a first step toward a realistic microgel modeling which fully incorporates and takes into account the polymeric nature of the particles.

Materials and Methods

Microgel Synthesis and Experimental Characterization

We specifically focus on the modellization of PNIPAM microgels synthesized via precipitation polymerization through the combination of 0.1929 g of sodium dodecyl sulfate (SDS, Duchefa Biochemie), 1.471 g of NIPAM, and 0.0647 g of BIS in 96.29 g of water. NIPAM was recrystallized in hexane, and all other chemicals were used as received. The mixture was heated to 70 °C and bubbled with argon, and 0.0539 g of KPS in 2.0145 g of water was added to start the reaction. The reaction was then left for 6 h under an argon atmosphere. The particle suspensions were cleaned by three centrifugation and redispersion series. The cross-linker concentration is thus 3.2 mol % with respect to NIPAM.

The characterization of the particles was carried out through dynamic light scattering (DLS) in pseudo-2D cross-correlation mode with laser wavelength λ = 660 nm (LS instruments, Switzerland). DLS measurements were performed over a range of temperature in steps of 2 K to yield an accurate swelling curve. The hydrodynamic radius was extracted using a first-order cumulant analysis averaged over an angular range of 60°–100° every 10°. The microgels have a hydrodynamic radius RH ∼ 24 nm under the most swollen investigated conditions (T = 282 K).

A sketch of the procedure employed in this work to synthesize in silico microgels. The monomer and cross-link molecules (leftmost panel) are initially modeled as patchy particles with two and four patches, respectively. This is shown in the central panel, where bivalent patchy particles (mimicking NIPAM monomers) are colored in blue, tetravalent ones, acting as BIS cross-linkers, are red, and patches are represented in gray. After the network has assembled, the resulting configuration is used as a starting point for simulations performed with a bead–spring model (rightmost panel).

Design, Assembly of the Network, and Swelling Simulations

To build the initial microgel configuration, we consider a mixture of bivalent and tetravalent patchy particles under spherical confinement. Each patchy particle is modeled as a sphere whose surface is decorated by either four patches of species A arranged on a tetrahedron or two patches of species B placed on the poles. A patch μ on particle i is identified by the unit vector . The two-body interaction potential between particles i and j readswhere r is the vector connecting i and j, r is its length, {p} is the set of patches of particle k, and rμν is the distance between patch μ on particle i and patch ν on particle j. The two contributions VWCA and Vpatchy encode a short-range repulsion and a short-range attraction, respectively. The repulsion is modeled with a Weeks–Chandler–Andersen potential[39]where σ is the particle diameter, which is taken as the unit of length, and ϵ controls the energy scale. The patch–patch interaction readswhere σ sets the position of the attractive well (of depth ϵμν) and rc is chosen by imposing Vpatchy(rc) = 0. Here we set σ = 0.4 and rc = 1.5σ. In addition, we set ϵAB = ϵBB = ϵ and ϵAA = 0, so that only bonding between BB and AB patches is possible. The two-body interaction given in eq is complemented by a three-body potential acting on triplets of close patches.[40] This term provides an efficient bond-swapping mechanism that makes it possible to easily equilibrate the system at extremely low temperatures, while, at the same time, retaining the single-bond-per-patch condition. Additional details can be found in ref (40).

We focus on systems of NA tetravalent particles, acting as cross-linkers, and NB bivalent particles, mimicking the NIPAM monomers. We set NA + NB = 42 000 and NA = 0.032(NA + NB) = 1344 to impose the same cross-linker concentration as in experiments. Each monomer corresponds to a bead of size compatible with the Kuhn length[30] of PNIPAM polymers,[41,42] i.e., σ ∼ 1 nm. Thus, the total number of monomers employed for the assembly has been estimated by imposing that the volume of the microgel particles is the same in experiments and simulations. We perform molecular dynamics (MD) simulations at fixed temperature T* = kBT/ϵ = 0.05, where kB is the Boltzmann constant. Thanks to such a low temperature, the system tends to maximize the number of bonds. In addition, owing to the bond-swapping mechanism, the system is able to continuously restructure itself, until the large majority of possible bonds are formed.

To ensure that the resulting network retains the spherical shape of a microgel colloid, we confine particles in a sphere of radius Z, with no periodic boundary conditions applied. The simulations are left running until the great majority of the particles (95–99%, depending on the system parameters) have self-assembled into a single network-like cluster. The assembled microgels contain a negligible fraction of dangling ends, since nearly all the bivalent particles belonging to the largest cluster are involved in two bonds. Indeed, the ratio between the number of formed bonds and the total number of possible bonds is always larger than 0.998. This assembly stage takes from a few hours to a few days on K80 NVIDIA GPUs, depending on cross-linker concentration and confinement. In particular, the aggregation rate is slower the fewer the cross-linkers or the larger the confinement. At the end of the procedure, we remove all those particles that are not in the largest cluster and use the resulting structure (typically comprising N ∼ 41 000 particles) as the initial configuration for the bead–spring model simulations.

Once the network is formed, we preserve its topology by replacing the directional interactions among the particles with a standard set of interactions designed to reproduce the behaviour of polymers.[43] To mimic the interactions between polymers, we adopt the classic bead–spring model in which bonded monomers i and j interact through the sum of a WCA potential (the same as eq ) and a FENE potential VFENE:[43−45]where kF = 15 is the spring constant and R0 = 1.5 is the maximum extension of the bond. In this way the reversible bonds formed by patchy particles are replaced by permanent bonds between connected monomers which preserve the initial network topology at all times. Non-bonded monomers interact only through the WCA potential in eq .

In order to grasp the swelling behavior of PNIPAM microgels with temperature, we need to take into account the tunable quality of the solvent. We do so in an implicit way by adding a term Vα to the interactions between all monomer pairs. The strength of this term is set by the parameter α, which controls the solvophobicity of the monomers and plays the role of a temperature.[45,46] This term readswhere we set γ = π(2.25 – 21/3)−1 and β = 2π – 2.25γ.[45] In the α = 0 limit, the system is in good solvent conditions, and no additional attractive interactions are at work. We perform MD simulations of a single microgel particle at constant temperature T* = 1.0 using the Nosè–Hoover thermostat and a leapfrog integration scheme with a time step δt = 0.001, with the same units as described above. We increase α until a complete collapse of the particles is observed.

We note that in the standard chemical synthesis[47] the polymerization kinetics takes place at temperatures ∼70 °C to facilitate aggregation. It is well-known that under these bad solvent conditions the cross-linkers react faster than the NIPAM monomers, soon forming a rather homogeneous branched core, which later incorporates additional chains and smaller aggregates. Other components, such as surfactants, are often added to obtain particles with small sizes.[47] Other parameters, such as the initiator concentration,[48] may influence the synthesis and the final density profile in ways that are not fully understood yet.[49,50] Therefore, we stress that in our approach we do not aim to reproduce the experimental protocol for the synthesis in its different aspects, but rather to adjust the several involved parameters, such as the confinement but also the patchy potential between different species, in order to achieve network realizations that are as similar as possible to the experimental particles in their basic properties, such as swelling behavior and elastic properties. The latter will be addressed in future work.

Finally, we compare the structural properties of our microgels with those obtained from a microgel generated out of an ordered network based on a diamond lattice, as recently done in refs (32−34 and 36). In this approach the sites of the lattice represent the cross-linkers, and each two sites are connected by a polymer chain. The spherical shape of the microgel is then obtained by cutting the network along a spherical surface. By construction, all chains have the same number of monomers. Consequently, the size of the microgel cannot be tuned by maintaining the same total number of particles. In the following we will employ a diamond-lattice-based microgel with N = 42 000 monomers and c = 3.2%.

Results

Swollen Regime

We generate microgel configurations for different values of the confinement radius in the range 30σ < Z < 70σ. For each value, we replicate the assembly protocol a few times (2–4 independent realizations are already enough to suppress any significant numerical noise) in order to average properties over independent conformations. After a fully bonded network is obtained, the patchy interactions are replaced by the bead–spring ones (with α = 0) until equilibration is obtained.

Figure shows typical equilibrated snapshots of the microgels for different choices of cross-linker concentration c and Z. Immediately we notice that the qualitative effect of reducing the number of cross-linkers is the same as that of increasing the confinement radius. In both cases the microgel tends to be larger and less compact, even though the properties of the gel networks should be quantitatively different. However, the reduction of the confining volume has a more dramatic effect on the size of the microgels than the increase of the number of cross-linkers for the range of parameters considered here. In addition, as the confining radius gets larger and larger, the conformation of the microgels becomes more irregular. Indeed, even though it retains on average a spherical shape, there are more and more portions of long chains that stick out of the corona before looping back toward the microgel core. In the framework of the simplified core–corona model often invoked to describe microgel behavior,[51] it is expected that the corona properties, particularly the elasticity and the effective interactions, will be largely affected by these outer chains. Indeed, recent work by Boon and Schurtenberger[52] has shown that these need to be taken into account to fully describe the experimental density profiles in the corona. Our approach is thus able to incorporate this effect.[53]

Figure 2

Equilibrated simulation snapshots of microgels generated at (a) c = 1.4%, (b) c = 3.2%, and (c) c = 5.0% and three different confinements: from left to right, Z = 30σ, 50σ and 70σ.

Regarding the distribution of chain sizes inside the network, Figure shows N, the number of chains of size l, which is defined as segments formed by l bonded bivalent particles, as a function of confinement Z. Quite strikingly, the distribution is found to be always exponential and independent of the confinement. This is a consequence of the equilibrium nature of the self-assembly process of the constituent patchy particles. Indeed, for bivalent particles it was found that N is always exponential, provided that the system undergoes an equilibrium polymerization.[54] Here, we show that this result holds even in the presence of cross-linkers and is not affected by the presence of the confinement.

Figure 3

Chain size distribution for c = 3.2% microgels at different confinement radius Z.

From now on, we focus on a fixed cross-linker concentration c = 3.2%, which is the same used in the experiments reported in this work. First we calculate the monomer density profile ρ(r) as a function of the distance r from the center of the microgel. This is shown in Figure a for several values of the confining radius. We see that for all values of Z we have a well-defined core, characterized by a homogeneous, flat profile, followed by a corona where the density decays to zero with a slope that depends on Z. Therefore, the choice of the confinement allows us to tune the extent of the core and the corona profile. It is thus a useful parameter that can be adjusted to compare with experiments.

Figure 4

(a) Monomer density profiles for c = 3.2% microgels at different confinement radius Z as a function of the distance r from the microgel center of mass. Here symbols are simulations results; lines are obtained from the fits of the form factors using the fuzzy sphere model and applying eq . The inset shows the same data rescaled by the polymer density of the inner core, which changes with Z, on a log–linear scale. The shaded area shows the approximate extent of the region where the curves obtained through eq well-reproduce the numerical data. The value of r at which each data set intersects the dashed gray line defined by ρ(r)/ρ0 ≡ 10–3 provides an estimate of the hydrodynamic radius, RH. (b) Form factors P(q) for microgels generated with different confinements. Symbols are simulation data, and lines are fits to eq . The (dark red) thick lines in (a) and in the bottom-rightmost panel of (b) refer to a microgel generated by a diamond-lattice configuration.

We also calculate the density profile of a diamond-lattice-based microgel, for which it is possible to obtain only one density profile once the number of cross-linkers is fixed. This is shown in Figure a for a diamond-lattice microgel with N ≈ 42 000 and c = 3.2%. It is interesting to note that the density profile of the diamond-lattice microgel is characterized by oscillations of period ≈8σ inside the core, a signature of the ordered polymer network.

Except for very recent super-resolution microscopy measurements,[55] typically the density profiles in experiments are indirectly obtained from the form factors P(q), which can be measured by neutron or X-ray scattering as a function of the wavenumber q. In general, microgel form factors are well-described by the fuzzy sphere model,[56] which assumes that the microgel particle has a homogeneous density profile in the core, being then surrounded by a corona of decreasing density. A Lorentzian function is also incorporated in the model to account for the network fluctuations due to the presence of static and dynamic inhomogeneities in the polymer network. Thus, the fit expression for the form factors iswhich depends on several adjustable fit parameters: the particle radius R′, the smearing parameter σsurf which corresponds to about half the thickness of the fuzzy shell, the average correlation length of the PNIPAM network ξ, and the amplitude of the long wavelength contribution of the network fluctuations to the intensity I(0).[21] In simulations we directly calculate the form factors aswhere the angular brackets indicate an average over different configurations. The results are reported in Figure b, showing a consistent decrease of the amplitude and number of oscillations as Z increases, signaling that the internal structure becomes increasingly soft. Together with the numerical data, we also show the corresponding fits by eq , which provide clear evidence that the fuzzy sphere model fully captures the shape of the form factors for all values of the confining radius. We also report in the bottom-rightmost panel of Figure b the form factor of the diamond-lattice microgel. We find that the latter well overlaps at low and high q with the disordered microgel one generated for Z = 70σ but displays a different slope at intermediate q values, being almost flat until qσ ≈ 0.8, where a small peak appears. The latter comes from the cross-linkers, which are spatially correlated due to the underlying ordered structure of the network; indeed, the length scale associated with the position of the small peak roughly corresponds to the width of the oscillations observed in the density profile. The flat behavior of P(q) at slightly higher q values is a consequence of such geometric correlations. Thus, the diamond-lattice microgel is able to reproduce the fuzzy sphere shape for small wave vectors and to retain the self-similarity of the chains at higher q values. However, it crucially fails to describe the network at intermediate scales due to its crystalline order.

According to the fuzzy sphere model, the normalized radial polymer density profile ρ(r)/ρ0, where ρ0 is the polymer density of the inner core, is:[56]where Rc = R′ – 2σsurf is the radius of the constant-density part of the particle, R′ is usually referred as the core radius, and R″ = R′ + 2σsurf is the total radius including the fuzzy shell. The latter is the estimate of particle radius that is commonly obtained from small-angle scattering experiments. From fitting the form factors, we extract the density profiles predicted by the model in eq and plot them as lines in Figure a together with those calculated from simulations (symbols). Quite remarkably, we find an excellent agreement between the two. At very small distances, a small residual noise is still present in the numerical data, and averaging over different realizations is crucial to correct the r → 0 behavior. We thus find confirmation that the fuzzy sphere model fully describes the microgel internal structure, providing an accurate estimate of the monomer density profiles.

In experiments the particle size obtained from the fits of the density profiles is commonly found to be systematically lower than the one estimated from dynamic light scattering measurements, which is usually referred to as the hydrodynamic radius RH. This systematic underestimation can be understood by looking at the density profiles close to the tail of the distributions. Indeed, as shown in the inset of Figure a, there is a residual tail, approximately exponential at large distances, which is not captured by the fuzzy sphere model. We adopt a qualitative threshold for determining the hydrodynamic radius: we define RH as the distance r at which the normalized density profile becomes of the order of 10–3. The values for RH extracted in this way are plotted together with R″ and Rc in Figure . As expected, all radii increase monotonically with Z.

Figure 5

Microgel radii obtained from fits to P(q), R′, and R″ (see eq ) and directly extracted from simulations, RH, for different values of the confinement Z.

Temperature Behavior

After having discussed the properties of the microgel in the swollen regime, we analyze the temperature behavior by setting α > 0, thus in the presence of the solvophobic term Vα. After equilibration and for each value of the solvophobic parameter α, we compute the swelling curves by calculating the ratio between the radius of gyration of the particles for a given α, R(α), and its value at α = 0, that is, in the maximally swollen case, RgMAX. The radius of gyration is directly calculated from simulations as Rg = [1/N∑(r – rCM)2]1/2 being rCM and r the position of the center of mass of the microgel and of the ith bead, respectively. We have verified that Rg and R″ are proportional to each other for all studied values of Z and α. In order to compare experimental and numerical data, a relation between the parameter α and temperature must be established. To this aim, we scale the numerical data to match the volume phase transition temperature of experimental data, which, for PNIPAM microgels, is known to occur at T ≈ 307 K. Figure shows the resulting normalized radius of gyration Rg/RgMAX as a function of temperature T, together with experimental results. The comparison demonstrates that the solvophobic potential correctly captures the thermoresponsive behavior of microgels, allowing us to establish, for the present case, the linear relationship T = 280.23 K + 38.33 Kα. Figure reports several swelling curves at different confining radii Z, showing that the change in Rg increases by more than a factor of 2 if the microgel is generated in a sphere with radius Z = 70σ compared to the Z = 30σ case. This is a consequence of the looser structure of the microgel generated in the larger confining sphere which, consequently, has a larger available free volume in the swollen state compared to the one obtained with a smaller Z, which is more compact (see Figure ). Interestingly, the swelling curves for all confinements can be well fitted by the Flory–Rehner theory[57−59] which provides the equation of state of the microgel particle when the net osmotic pressure of the microgel is zero:Here ϕ0 and R0 are the particle volume fraction and the radius of gyration in the collapsed state, respectively, while Ngel is the average chain size, defined as the average number of monomers in a chain connecting two cross-linkers. In eq , χ(T,ϕ) is a series expansion in ϕ of the Flory parameter, which also depends on the volume phase transition temperature TVPT.[21] Inverting eq , an expression for the radius of gyration Rg as a function of T is obtained, which can be used as the fitting function.[60]

Figure 6

Normalized radius of gyration Rg/RgMAX as a function of α (top axis) and temperature (bottom axis). Open symbols are numerical results obtained for microgels with N ∼ 41 000 monomer beads and a concentration of cross-linkers = 3.2% for different values of the confining radius Z. Filled symbols are results from multiangle DLS experiments; here the radius of gyration is normalized by its value at the lowest measured temperature, T = 282 K. Solid lines are fits of the numerical data using the Flory–Rehner theory. The dashed line refers to the case of a microgel generated by a diamond-lattice configuration.

In general, the fitting function depends on all the free parameters listed above. We fix TVPT to 307 K as well as R0 ≃ 16.8σ and Ngel = 14, using the Z-independent values of these quantities calculated in simulations. With these choices, the resulting curves, represented as full lines in Figure , are able to fully describe the numerical data.

The microgels obtained for different confining radii display quite a distinct topology among themselves. Even though Ngel is only weakly dependent on Z, a small confinement gives rise to a dense, possibly interwined, network, while large confinements allow particles to aggregate in a less intricate network due to the larger volume at disposal. We find that only the most strongly confined microgels, Z = 30σ, obey the same small variation in size observed in the experimental data for small microgels. Thus, Z represents an useful control parameter that can be tuned to obtain microgels with different topologies and density profiles, while maintaining the same percentage of cross-linkers. This feature is lost for instance for the case of a microgel with N ≃ 41 000 beads and c = 3.2% generated from a diamond lattice; the resulting curve is also reported in Figure (dashed line). It is evident that these data, despite referring to a microgel having the same number of monomer beads and the same percentage of cross-linkers, cannot describe the experimental swelling curve. Indeed, the microgel obtained from the diamond lattice displays a swelling curve that is too steep compared to the small-microgel experimental data examined in this work. We further note that since the topology of the microgel cannot be modified in the lattice synthesis approach, this is expected to also have important consequences on the mechanical properties at the single-particle level.

The successful comparison of the experimental swelling curve for small microgels with c = 3.2% and that of our numerically generated ones for Z = 30σ provides us with an estimate of the bead size σ. Indeed, if we compare the experimental hydrodynamic radius at the lowest T, RH = 24 nm, with the numerical counterpart estimated for α = 0, i.e., RH ≈37σ, we obtain that our bead size is σ ∼ 0.65 nm. If R″(α = 0) is used instead, we would get a slightly larger value (σ ≈ 0.72 nm). Both values are consistent with the expected Kuhn length of PNIPAM polymers.[41,42] Thus, the choice to work with small microgels allowed us to avoid coarse-graining and to quantitatively match the number of cross-linkers in the simulations. However, we need to point out two drawbacks. The first one is that the experimental microgels are so small that their effective interaction, which can be approximated by Hertzian model,[61] is very soft. Indeed, we estimated that at full overlap they do not exceed ∼100kBT.[62] Hence, when brought above the VPT temperature, van der Waals attraction quickly takes over, leading to aggregation into large clusters. This is the reason for the lack of experimental data points for T ≳ 305 K. Second, due to the small size of the microgels, SAXS measurements of the form factors have not been carried out, preventing a direct comparison with the numerical P(q). Both these aspects will be tackled in the future by addressing the case of larger microgels with hydrodynamic radius of the order of a few hundred nanometers. We note that in that case the monomer-resolved representation would require ≈107 monomers. Thus, to undertake an investigation of such large particles, future approaches will need to be complemented by some coarse-graining procedure and by the analysis of system size dependence.

We have shown in Figure that the microgel particles generated in a confining sphere of radius Z = 30σ are the best candidates to reproduce the swelling properties of small microgels. We now focus on this value of Z and investigate the behavior of P(q) as a function of T. The results are obtained after averaging over distinct microgel realizations. In Figure b, we report P(q) for different values of T. Similarly to what found in experimental form factors,[48,63−65] we find that, on increasing T, the peaks of P(q) move to higher values of qσ, signaling a decrease of the size of the microgel. In addition, the number of peaks increases with T, a signature of the fact that the microgel becomes more and more compact. Snapshots of the corresponding microgels are shown in Figure a, providing a visual confirmation of the occurrence of the VPT. The form factors are well described by the fuzzy sphere model of eq at all T, mimicking the excellent agreement exhibited for the purely swollen case. The density profiles extracted from the fuzzy sphere fits of P(q) eq are shown in Figure c and are compared with those calculated numerically. Once again, a very good agreement between calculated and estimated density profiles is observed, confirming that the fuzzy sphere model is adequate to describe the microgel structure at all temperatures.

Figure 7

(a) Snapshots of a microgel particle exhibiting the typical volume phase transition from swollen to compact. (b) Numerical form factors, averaged over four different realizations, for microgels with N ∼ 41 000 monomers and c = 3.2% of cross-linkers generated in a sphere of radius Z = 30σ (symbols) at different values of the solvophobic parameter α, corresponding to the reported temperatures as obtained from the swelling curve comparison to experiments in Figure . Solid lines are fits of the curves using the fuzzy sphere model of eq . (c) Averaged density profiles obtained from MD simulations of a microgel with N ∼ 41 000, c = 3.2%, and Z = 30σ (symbols) and from the fit of the form factors using the fuzzy sphere model (solid lines).

Conclusions

In this work we have introduced a novel approach to synthesize microgels in silico. We exploit the self-assembly properties of patchy particles to reproduce the polymerization of NIPAM monomers (two-patch particles) and BIS cross-linkers (four-patch particles). The two are fixed in concentrations by the molar ratio used in experiments. We specifically synthesized and measured the swelling curve of unusually small microgels, with hydrodynamic radius of the order of 24 nm in the swollen regime, in order to avoid complications due to coarse-graining, that will be tackled in the future. By using an elegant method for achieving equilibration at very large attractive strengths,[40] we are able to build a fully bonded network. While the kinetics of the network formation is clearly different from the experimental one, our aim is to obtain a particle which resembles the ones produced in laboratory as much as possible possible. We find that an important control parameter of our numerical synthesis is the volume in which particles initially form, and to this aim we confine the patchy particles in spheres of different radii. The confining radius Z turns out to crucially control the network organization. Indeed, a strong confinement gives rise to a very intricate network with lots of entanglements, which is able to swell and deswell much less than a microgel generated in a larger volume. Once the network is formed, we preserve its topology by substituting the patchy model, useful to obtain the initial assembly, with a standard bead–spring polymer model. We further introduce the temperature dependence by adding a solvophobic term in the potential.[45]

Our computer-generated microgels are shown to share the main characteristics of their experimental counterparts. Their internal structure is characterized by a homogeneous core followed by a soft corona. The extension of the corona can be controlled in the preparation of the initial network by changing the confining radius. The microgel form factors are also shown to closely follow the fuzzy sphere model, widely used to describe the experimentally measured form factors. In addition, we have performed a consistency check for the density profiles, comparing the results of the direct calculations with those resulting from the fuzzy sphere assumption, finding perfect agreement for all investigated Z.

We next addressed the thermoresponsive properties of the microgels by comparing the swelling curve to the experimentally measured one for our test-case small microgels. We found that the experimental data can be well described only by microgels generated under strong confinement. This is likely due to the fact that the size of the experimental microgels is close to the lower limit imposed by the standard synthetic protocols. Thus, we expect that future comparisons with larger microgels will require less extreme confinements. Indeed, while the interpretation of the confinement is not straightforward in terms of the actual chemical synthesis, it surely correlates with the size. For the investigated experimental microgel we are able to reproduce the full temperature behavior and the occurrence of the volume phase transition, whose location is used to establish a relationship between the model solvophobic parameter α and the temperature. However, a deeper investigation will be required to assess how the chosen parameters in the hydrophobic interaction affect the swelling behavior of the microgel particles. Finally, the internal structure of the microgels at different T is well-reproduced by the fuzzy sphere model, once again confirming the potential of our model to faithfully describe experimental data. The numerical protocol we developed was also shown to perform better than the ordered lattice approach used so far, which is not able to quantitatively reproduce the experimental swelling curve.

Unfortunately a direct comparison with experimental form factors for the investigated small microgels was not possible due to the difficulty of the measurements. We will address this key issue in the near future by devising a coarse-graining strategy to tackle the assembly of microgels one order of magnitude larger in size (and thus roughly three orders of magnitude bigger in the number of particles). A careful design of the initial patchy particles and the conditions under which they will self-assemble will be crucial to achieve this goal. The present work is thus the first necessary step toward the more challenging task of designing a fully fledged realistic microgel in the computer.

Our promising new approach to microgel computer simulations will also be crucial to evaluate effective interactions and to go beyond the widely used, but oversimplified, Hertzian model.[26] It will also allow us to carefully evaluate the role of dangling chains and dangling ends in determining the single-particle elastic properties as well as interparticle interactions. In addition, our in silico synthesis protocol will make it possible to elucidate the important role of chain entanglements in the swelling behavior as well as in the mechanical properties of microgels. Indeed, we have shown that their abundance in our particles is likely controlled by the radius of the confining sphere used during the generation of the network. This is a topic that has recently gained a lot of interest, and a variety of methods[66,67] are being developed to evaluate entanglement effects, for which our model will provide an interesting test case. Finally, starting from the neutral (at least in the swollen case) PNIPAM microgels case, we will extend the model by adding charges and, in a more ambitious project, to study interpenetrated networks microgels,[68−70] which have been shown to exhibit an intriguing fragile-to-strong transition,[68] a very uncommon feature in colloidal systems.