Numerical insights on ionic microgels: structure and swelling behaviour

Abstract

Recent progress has been made in the numerical modelling of neutral microgel particles with a realistic, disordered structure. In this work we extend this approach to the case of co-polymerised microgels where a thermoresponsive polymer is mixed with acidic groups. We compare the cases where counterions directly interact with microgel charges or are modelled implicitly through a Debye–Hückel description. We do so by performing extensive numerical simulations of single microgels across the volume phase transition (VPT) varying the temperature and the fraction of charged monomers. We find that the presence of charges considerably alters the microgel structure, quantified by the monomer density profiles and by the form factors of the microgels, particularly close to the VPT. We observe significant deviations between the implicit and explicit models, with the latter comparing more favourably to available experiments. In particular, we observe a shift of the VPT temperature to larger values as the amount of charged monomers increases. We also find that below the VPT the microgel–counterion complex is almost neutral, while it develops a net charge above the VPT. Interestingly, under these conditions the collapsed microgel still retains a large amount of counterions inside its structure. Since these interesting features cannot be captured by the implicit model, our results show that it is crucial to explicitly include the counterions in order to realistically model ionic thermoresponsive microgels.

---

1 Introduction

Microgels are colloidal scale polymeric networks that can be dispersed in a good solvent.¹ They have recently become a favourite model system^2–5 thanks to their intrinsic softness and to the possibility to respond to external stimuli with changes in size. Such a phenomenon is commonly referred to as Volume Phase Transition (VPT)¹ and it is controlled by the properties of the constituent polymers. The prototype example is given by poly(N-isopropyl-acrylamide), PNIPAM, a thermoresponsive polymer which gives microgels in water the ability to reversibly increase or reduce their size after a change of temperature around the so-called VPT temperature T_VPT ∼32 °C. Another interesting case can be realized using pH-responsive ionic polymers, made of weak acidic or weak alkaline monomers. The resulting ionic (or simply called charged) microgels are able to adjust their bare charge in response to a pH variation by releasing H⁺ or OH⁻ ions due to the dissociation of a fraction of monomers.⁶

Out of the many possibilities provided by modern-day synthesis methods, co-polymerised PNIPAM-co-PAAc microgels are of particular interest,^6–9 as they combine the thermoresponsive properties of PNIPAM with the pH-responsive features of polyacrylic acid (PAAc), stemming from the weak acidic nature of AAc monomers. Indeed, at low pH almost all AAc monomers are not dissociated because of the high concentration of H⁺, which favours the inverse recombination reaction that leads to an almost neutral network. On the other hand, for high pH values, most of the acidic monomers dissociate, generating a charge distribution throughout the particle volume. It is important to note that the fraction and the distribution of the charges within the network depend on the chosen experimental conditions, such as the packing fraction, the specific molecular interactions, the local counterions concentration and the electrostatic interactions between nearest charged monomers, which can be optionally mediated by the presence of salt.^1,9,10

The multiresponsive character of ionic microgels makes them highly versatile. They are indeed responsive also to external alternating electric fields, through which their mutual interactions (and hence their phase behaviour) can be tuned.^11,12 Their single-particle properties have been extensively investigated in experiments as a function of both temperature and pH.⁹ Microgels with different content of AAc obtained through several synthesis methods have been analysed in order to assess the effects of inhomogeneities in the distribution of crosslinkers and charged monomers.^7,8 The tunability of ionic microgels has also been exploited in several fields of research, from biology¹³ to laser technology.¹⁴ For instance, their dual responsiveness makes them highly suitable to be employed in the smart design of optical switch devices³ based on colloidal photonic crystals.

Understanding the effects of electrostatic interactions in ionic microgels could also shed light on the behaviour of other kinds of microgels. Indeed, even those constituted by PNIPAM only show interesting features determined by the presence of charges, particularly for high concentrations^15,16 or above the VPT temperature.^10,17 In addition, microgels consisting of two different interpenetrated networks, made of PNIPAM and PAAc respectively, have recently gathered a lot of attention because of their suitability to study the problem of fragility in structural glasses.^18,19

From the theoretical point of view, several investigations of the swelling of charged microgels, mostly relying on a mean-field treatment of the polymer network based on the Flory–Rehner theory,^20,21 have been reported. In these works, electrostatic effects and steric interactions due to the presence of counterions have been taken into account by approximated theories such as the Poisson–Boltzmann equation,^22,23 the Ornstein–Zernike integral equation²⁴ and density functional theory.²⁵ Also an effective interaction potential has been derived using linear response theory,²⁶ which made it possible to draw a phase diagram for ionic microgels.²⁷ On the numerical side, the use of coarse-grained models^12,28 has allowed to go beyond the mean-field framework and tackle the behavior of ionic microgels at all concentrations. However, in order to refine the highly coarse-grained models required to study the bulk properties of microgel suspensions, it is important to first correctly capture the single-particle behavior, a task that has been tackled only relatively recently^29,30 due to the high computational cost of numerical studies reproducing microgels at monomer-resolved level.

The inclusion of long-range electrostatic interactions on complex objects such as microgels is a challenging and numerically demanding task, particularly if counterions are explicitly considered. Therefore, in several cases, an implicit treatment of counterions, for example based on the Debye–Hückel theory, has been employed to make it possible to perform simulations of relatively large systems.^31,32 However, a few numerical investigations have also been carried out in the explicit presence of the counterions. A pioneering work reported coarse-grained simulations of polyelectrolyte gel networks,³³ while simulations of single nanogel particles have appeared only later on.^34–39 Several techniques have been devised to treat charged networks. Particularly, recent Monte Carlo simulations^40–42 have been carried out to provide a coarse-grained description of the dissociation reaction on a statistical basis. These studies concluded that all investigated macroscopic properties mostly depend on the number of charges, rather than on their distribution, in agreement with experimental observations.⁸ Notwithstanding this, all coarse-grained studies of ionic microgels have so far been performed with networks built out of ordered topologies, e.g., based on the diamond lattice, which cannot take into account the disordered nature of real polymer networks.²⁹

In order to go beyond mean-field and to account in a more realistic way for the effect of the network topology, in this work we perform extensive simulations of charged microgels modelled as disordered networks. We start by preparing neutral microgel configurations following previous works,^43,44 ensuring that the internal microgel structure reproduces the swelling behavior and form factors of experimental non-ionic microgels. Then, we add a quenched charge distribution, varying the fraction of charged monomers that are randomly distributed throughout the network. Since the probability that a monomer is charged is lower near crosslinkers,^40,41 we add the constraint that the latter are always neutral. To account for charge–charge interactions we perform two different kinds of simulations: (i) we rely on the Debye–Hückel model in which charged monomers interact implicitly through a two-body Yukawa potential and (ii) we explicitly include counterions as charged coarse-grained particles. We calculate the density profiles and form factors of the microgels for both approaches and average over different charge realizations. We simulate microgels in swollen conditions and across the volume phase transition by using a solvophobic interaction between the monomers that models the different quality of the solvent as temperature varies.^43,44

Our work is important to understand the effects that inhomogeneous topologies and charge distributions beyond mean-field can have on the single-particle behavior of ionic microgels, filling a gap in the current literature. In addition, we provide significant insights on the difference between neutral and charged microgels across the volume phase transition. Indeed, the competition between the electrostatic repulsion and the solvophobic attraction, which develops at intermediate temperatures in between the swollen and collapsed regimes, could be important for the arising of a distinct phenomenology in the presence of charges. Finally, our work can be considered as a starting point for future investigations at finite concentrations, shedding light on the deswelling behavior of ionic microgels, which takes place at concentrations below the overlap one.^7,28

2 Models and methods

2.1 Monomer interactions

To analyse the role of charges on the single-particle properties and on the swelling behaviour of microgels, we exploit a recently proposed numerical protocol⁴³ to generate disordered, heterogeneous microgels that are structurally similar to real neutral ones. We start by preparing fully connected spherical networks by confining patchy particles in a cavity. We always consider an additional designing force acting only on the crosslinkers, which provides the typical core–corona structure of realistic microgels.⁴⁴ Once the network is formed, we freeze the topology of the network and adopt a monomer-resolved approach.⁴⁵ The beads that make up the polymers interact via a steric repulsion, modeled with the Weeks–Chandler–Anderson (WCA) potential:where ε and σ are respectively the energy and length units. In addition, connected beads interact via the Finitely Extensible Nonlinear Elastic potential (FENE):where l₀ sets the maximum bond distance and k_F is a stiffness parameter influencing the rigidity of the bond and the equilibrium bond-distance. This potential ensures that no covalent bonds between the monomers can be broken during the course of the simulations. In all cases, we use k_F = 15 and l₀ = 1.5.

All monomers also interact with each other by means of an effective solvophobic potential, named V_α, which implicitly takes into account the monomer–solvent interactions:⁴⁶

with and .⁴⁶ This potential represents an effective attraction, modulated by the solvophobic parameter α, arising between thermoresponsive monomers at high temperatures. In other words, α plays the role of an effective temperature: α = 0 represents good solvent conditions, while with increasing α the quality of the solvent worsens, leading to the aggregation of beads and to the shrinking of microgel particles.^43,59

We complement this model by adding electrostatic interactions between charges that are randomly assigned to a fraction f of the microgel monomers. This choice aims to model the dissociation of weak electrolyte groups, usually giving rise to negatively charged microgels, such as when acrylic acid is used as a co-monomer in the synthesis process. The neutrality of the overall suspension imposes the presence of positively charged counterions, which balances the total charge of the microgel–counterion complex. In the simplest approach, the effect of charges can be taken into account by using the Debye–Hückel potential, which models the charge–charge interaction as a screened Coulomb (or Yukawa) potential acting between each pair of charged beads as:⁴⁷

where λ_B and λ_D are the Bjerrum and the Debye lengths, respectively. The former represents the distance at which two ions of valence z feel a repulsive energy exactly equal to k_BT, thus quantifying the relative intensity of the electrostatic forces, and it is defined as:where ε₀ and ε_r are the vacuum and relative dielectric constants and e is the elementary unit charge. The Debye length instead is the screening length, depending on both λ_B and on the density of counterions ρ_ci as:

λ_D = (4πλ_Bρ_ci)^−1/2. (6)

The Debye–Hückel approach can be used in principle only for symmetric electrolytes, i.e., when the valence of positive and negative ions is the same, as it is for the present case.⁴⁸ We work with reduced units, with σ, m, ε being the units of length, mass and energy, respectively. Within this unit system, the experimental Bjerrum length, that is λ_B ≈ 0.7 nm for monovalent ions in water at room temperature, translates into a reduced Bjerrum length λ_B* ≈ 1, assuming σ ∼ 1.0 nm comparable to the Kuhn length of both neutral NIPAM and charged AAc monomers. This can be considered as a lower estimate of σ, according to different types of measurements for linear PNIPAM chains.⁴⁹ Note that the use of a larger value of σ would significantly decrease λ_B, thus resulting in a very small effect of the Debye–Hückel repulsion as compared to the neutral case.

Although the Debye–Hückel model is suitable to implicitly treat the role played by counterions in homogeneous systems and in the effective interactions among colloidal particles in dispersions, it should be avoided when studying electrostatic ion–ion interactions within inhomogeneous weak-electrolyte systems such as charged polymeric particles. In particular, one of the drawbacks of using this approach is that, for weak polyelectrolytes, there is not a simple link between the pH and the dissociation fraction of the acidic monomers, which determines the value of λ_D.⁵⁰ Moreover, this model cannot take into account other relevant effects due to the presence of counterions, such as their osmotic pressure. In order to overcome these issues it is crucial to explicitly take into account the counterions and thereby to adopt an alternative model where all charged beads interact via the bare Coulomb potential, as:

For ion–ion interactions this term is complemented by a steric repulsion, modeled again with the WCA potential (eqn (1)). This second approach significantly increases the computational cost of the simulations, but at the same time it yields a realistic representation of the counterion distributions within the network, which is important to correctly describe the behavior of the microgels across the volume phase transition. This type of study calls for some preliminary investigations, that are described in detail in the ESI.† In particular, we analyzed the dependence of our results on the choice of the simulation box (see Section S1, ESI†), discovering that there is a critical size of the box below which the long-range electrostatic forces are not correctly taken into account. In addition, we explored the role of the counterions diameter σ_c on the microgel swelling behavior (see Section S2, ESI†), finding that the use of too large counterions yields unrealistic excluded volume effects in the collapsed state of the microgel. We thus fix σ_c = 0.1σ throughout the rest of the manuscript.

2.2 Numerical simulations

We perform Molecular Dynamics simulations of single microgels with N ∼ 42 000 monomers at fixed crosslinker concentration c = 5%. Microgels are assembled as in ref. 44 in a spherical cavity of radius R₀ = 50σ, yielding an internal structure of the microgels which compares very well with experimental ones obtained through radical polymerisation techniques at the same value of c. Once a fully connected network is assembled, we assign a charge ze = 1 to a fraction f of the monomers, that are randomly chosen throughout the microgel. We then maintain this charge distribution fixed throughout the simulation run and we average results over four different network topologies and three different realizations of the charge distribution.

We study microgels for three different charge fractions, f = 0.05, 0.20, 0.95, and for several values of the solvophobic parameter α across the VPT. The equations of motion of the system are integrated via the velocity-Verlet algorithm.⁵¹ The equilibration of the system is carried out in the canonical ensemble using the Nosè–Hoover chains thermostat for 1.6 × 10⁶ simulation timesteps, while a long production run in the microcanonical ensemble of ∼2 × 10⁶ steps is used to obtain equilibrium averages of the thermodynamic observables under investigation. We used a cut-off of R_cut = 5λ_D for the Debye–Hückel potential, whereas the long-range Coulomb interactions are computed with the particle–particle–particle–mesh method.⁵² For the latter type of simulations we used the LAMMPS package.⁵³

2.3 Main observables

To assess the microgel size, we calculate the radius of gyration, defined as:where r⃑_i and r⃑_CM are the positions of the i-th monomer and of the microgel's center of mass, respectively.

To gain a better knowledge of the inner structure of the microgel we calculate its radial density profile, defined as the average density at a fixed distance from the center of mass:

where the brackets 〈·〉 indicate ensemble averages. We also compute the density profile of charged monomers, labelled as ρ_CH(r), and that of counterions only, labelled as ρ_CI(r). By adding the two latter quantities, weighted by the respective charge, we obtain the net-charge density profile ρ_Q(r) = −ρ_CH(r) + ρ_CI(r), which provides information on the charge distribution throughout the volume of the particle.

The counterpart of the density profile in Fourier space is the form factor P(q), which can be readily obtained in neutron or X-ray scattering experiments of dilute microgel suspensions. In simulations P(q) can be directly calculated as:

We have computed the rotationally invariant quantity P(q) as an average of P(q⃑) over 300 vectors q⃑ randomly picked onto a spherical surface of radius q.

Usually, experimental and numerical data of P(q) for neutral microgels are described by the fuzzy sphere model,⁵⁴ which is able to account for particles with a homogeneously dense core and a fuzzy corona, wherein the density gradually decreases away from the center of mass. This results in a density profile , where Erfc(·) is the complementary error function, while R and σ_surf are related to the extension of the core and of the corona, respectively. However, it has recently been suggested by super-resolution microscopy measurements⁵⁵ that the assumption of a homogeneous core may not be too accurate, and that the density profile of microgels could be better approximated by the function , which includes a linear growth of the monomer density inside the core modulated by the parameter s. Our microgels have thus been assembled through a numerical protocol that is able to reproduce such features.⁴⁴ The additional linear term in the density profiles modifies the shape of the form factor in an extended fuzzy sphere model:

This functional form is usually added to a Lorentzian term which takes into account the inhomogeneities of the network at large q. However, such a term was often found to be unsatisfactory in comparison to available experiments, especially for hydrogels.⁵⁶ A step forward is represented by the modified Lorentzian proposed by Shibayama and Tanaka,⁵⁶ which relies on the assumption that the spatial correlations of the network decay according to r^D−d, where d is the system physical dimension and D is the fractal dimension of the correlated domains. For large q the form factor can thus be written as:with ξ being the length over which concentration fluctuations are spatially correlated.

3 Results and discussion

3.1 Swollen microgels

In this section we discuss the properties of microgels in good solvent conditions. In our model this corresponds to α = 0, i.e., to monomers that interact via the bead-spring model plus the charges contribution only. To quantify the latter, we analyze both the Debye–Hückel approach and the simulations in the presence of explicit counterions, carrying out a comparison between these approaches and the neutral case.

3.1.1 Debye–Hückel microgels. We start by reporting in Fig. 1 the microgel radius of gyration R_g for the Debye–Hückel model as a function of the screening length λ_D for three different values of f. Data are normalized with respect to the neutral microgel case, for which f = 0. For all considered values of f, the microgel size increases with λ_D. We observe a progressive increase of the microgel size as f increases, with the fully charged microgel, which corresponds to f = 0.95 since crosslinkers are not charged, displaying the strongest variation of R_g with respect to the uncharged case. The fully charged situation was also analyzed in ref. 31 and 32 for a diamond-like microgel and we find comparable variation of the microgel size to that reported in these works.

Fig. 1 Gyration radius R_g as a function of Debye length λ_D for different values of f, up to the fully charged case (f = 0.95). Data are normalized with respect to the neutral microgel (f = 0).

To visualize the effect of charges on the internal structure of the microgels, we report in Fig. 2 the density profiles and the form factors of the microgels for a representative value of λ_D and different values of f, from the neutral case up to the fully charged one. As expected, we find that a larger presence of charges has the effect to lower the density of monomers in the core region and consequently to increase it in the corona, as shown in Fig. 2(a). Such a variation of the density profile is barely noticeable for f = 0.05 and very moderate for f = 0.20. However, the fully charged case displays a considerably different profile, where the core density is about half of that in the neutral case and the corona extends to distances larger by about 50% with respect to the neutral case. Small oscillations at short distances r disappear when averaging over a larger number of realizations of network topologies.^43,57

Fig. 2 (a) Monomer density profiles and (b) form factors as a function of f at fixed λ_D = 1.5σ for the swollen (α = 0) microgel, from neutral (f = 0) to fully charged (f = 0.95) conditions. For a matter of completeness we also show the statistical error on density profiles, which is appreciably large only for small values of r, because the sampling region of space is small for those points; in the plots that follow we omit the error bars for the sake of clarity. In (b) data are shifted on the vertical axis by a factor of 3 with respect to each other to improve visualization.

Corresponding P(q) are reported in Fig. 2(b) showing again tiny changes from f = 0 to f = 0.20: the first peak slightly shifts to smaller wavevectors, reflecting the larger size of the microgel, but no additional peaks are observed. In addition, the slope of the curves at high q remains the same. The case f = 0.95 shows the same features, but amplified by the large number of charges. Interestingly, we can compare the results in Fig. 2, with those reported in ref. 31 for a fully charged diamond lattice network where charges are also modelled by a Debye–Hückel potential. In that work, regular oscillations in the density profiles were observed, due to the underlying presence of a regular mesh of the network, as also discussed previously for non-ionic microgels.²⁹ Such oscillations were further enhanced in the presence of charges, leading to unrealistic density profiles. Similarly, the form factors were found to display strong deviations from the fuzzy sphere model, displaying a minimum at intermediate q. Such features are totally absent in the disordered network model examined here, regardless of the amount of charge. These results confirm once more the importance of a correct modeling of the underlying network topology to treat single-particle microgel properties, also for charged microgels.

We notice that at high pH the average fraction of ions that dissociate from the microgels may be considerably lower than the ideal one for dilute suspensions of AAc, resulting in a larger average distance between charged monomers.⁵⁰ This poses concerns about the use of too large values of f, which would be unrealistic under these conditions. Indeed, if we look more carefully, we notice that P(q) for f = 0.95 displays a sort of kink for qσ ∼ 1. Looking at the snapshots of the corresponding microgel (not shown), evident holes appear in the structure with a size comparable to this length scale, suggesting that such high-charge conditions are probably far from realistic ones for standard co-polymerized microgels. For these reasons, in the following, we will consider only the f = 0.05 and f = 0.20 cases. We kept the case f = 0.95 in the foregoing analysis (i) to compare with previous simulation studies in which f = 1 was used³¹ and (ii) to appreciate the qualitative trends of the analysed physical quantities.

3.1.2 Microgels with explicit counterions. We extend our study to the explicit counterions (EC) model by focusing on two values of f = 0.05, 0.20. Fig. 3 reports the resulting density profiles comparing the explicit model results to the Debye–Hückel ones (DH) for different values of λ_D.

Fig. 3 Evolution of the monomer density profiles for the Debye–Hückel microgels (DH) with different λ_D (full lines) and for the model with explicit counterions (EC, solid circles) with (a) f = 0.05 and (b) f = 0.20, in the swollen state (α = 0). The net charge density profile ρ_Q(r) for the EC microgel (solid diamonds) is also reported (scale on the right axis). The density profiles of the corresponding neutral microgel (dashed lines) are shown for comparison. Inset: Same data in semi-log representation.

For f = 0.05 the two models yield similar results, probably due to the limited presence of charged monomers. However, for f = 0.20 the microgel with explicit counterions exhibits a more extended corona than the Debye–Hückel model for all investigated values of the Debye length (see inset of Fig. 3). Even a large increase of λ_D, which has a qualitatively similar effect to the increase of f (since we find fewer monomers in the core and a more extended corona), gives rise to results that do not superimpose onto the explicit counterions case, suggesting an intrinsic different structure of the microgels between the two models. In an attempt to set up an effective Debye–Hückel model that mimics the explicit one, we have calculated an effective screening length λ_D* from eqn (6) by substituting ρ_ci with the average density of counterions that is present inside the microgel with explicit counterions within a sphere of radius (2/3)R_g. Such a value roughly takes into account the whole extent of the core region. In this way, we obtain λ_D* ≃ 2.4σ for f = 0.05 and λ_D* ≃ 1.3σ for f = 0.2, respectively. The resulting density profiles of the λ_D*-microgels are reported in Fig. 3, showing also in this case a different behavior with respect to the explicit model. Deviations are larger in Fig. 3(b) for the higher fraction of charges considered, where the effective Debye–Hückel result is actually much more similar to the neutral case than to the explicit one.

The use of explicit counterions makes it possible to monitor the total charge density of the microgels ρ_Q(r), also shown in Fig. 3. For both values of f we find that the complex microgel–counterions is globally neutral at all length scales. Indeed, charge density profiles are much smaller with respect to the average inner densities of charged monomers both for f = 0.05 (ρ_Q ∼ 1.5 × 10⁻²σ⁻³) and f = 0.20 (ρ_Q ∼ 5.0 × 10⁻²σ⁻³). Thus, the counterions are able to freely diffuse throughout the microgels, even within the core, so that they fully counteract the electrostatic repulsion. The presence of the counterions inside the network thus contributes to the increase in the size of the microgel.

The behavior of the form factors is shown in Fig. 4. We start by discussing the results for f = 0.05 in Fig. 4(a), where only very minor changes to P(q) are observed and no shift of the first peak position is found. We find that all curves corresponding to the Debye–Hückel model are quite similar to the neutral case, independently of λ_D. The only noticeable difference is a weakening of secondary peaks in the presence of charges. The explicit model is the only one with a significantly smaller peak height and a different behavior at larger q, with some small residual oscillations and an apparently different slope at intermediate wavevectors.

Fig. 4 Form factors for the Debye–Hückel microgels (DH) by varying λ_D (solid lines) and for the model with explicit counterions (EC, solid circles) with (a) f = 0.05 and (b) f = 0.20 in the swollen state (α = 0). The form factors for the corresponding neutral microgel (dashed lines) are reported for comparison.

These features are amplified for f = 0.20, where now also a shift of the first peak position to smaller q values is observed. This is actually more evident for the implicit, rather than for the explicit model, which displays the smallest peak intensity. Again, secondary peaks are suppressed and now the appearance of a different slope for P(q) in the second peak region is more evident. Hence, we confirm that the Debye–Hückel model cannot be superimposed on the one with explicit counterions, even with the use of an effective Debye–Hückel model with λ_D = λ_D*.

The fact that the implicit Debye–Hückel model fails to reproduce the features observed in the explicit counterions case can be attributed to at least two reasons. First, the permeable and inhomogeneous structure of microgels as well as the presence of a rough interface among its inner part and the solvent generate uneven distributions of charges. These in turn lead to different screening conditions in different regions of the particle, that cannot be captured by the single lengthscale of the Debye–Hückel model. Second, the counterions have to balance the electrostatic attraction which drives them close to the charged monomers of the network, and the entropic gain that pushes them to leave the microgel, the latter being particularly strong for small-sized nanogels.⁵⁸ Under these conditions, it is not a priori trivial to assess the relative contributions to the swelling of the electrostatic interactions and of the counterions osmotic pressure, respectively. In addition, these considerations make such a kind of implicit treatment not readily applicable to the study of finite-concentration suspensions (beyond the scope of this paper), because of the complex dependence, in thermosensitive soft colloids, of the local counterions concentration on the effective packing fraction and on temperature.

3.2 Temperature-driven swelling of charged microgels with explicit counterions

In this section, we analyze in detail the deswelling behavior of the microgels with explicit counterions by adding the solvophobic potential V_α between monomers (eqn (3)) to mimic the increase of temperature in experiments.^8,9

3.2.1 Swelling curves and distribution of counterions. We start by showing the swelling curves of the microgels in Fig. 5, reporting the radius of gyration R_g as a function of the parameter α in the presence of explicit counterions for two different values of f. The behavior of the neutral microgel model is also reported for comparison. The first important observation is that the value of α at which the VPT occurs, i.e. α_VPT, defined as the position of the maximum of |dR_g/dα|, shifts from α_VPT ∼ 0.63 for neutral microgels,^44,59 to α_VPT ∼ 0.69 for f = 0.05 and up to α_VPT ∼ 0.82 for f = 0.20, as reported in the inset of Fig. 5. Using the α − T mapping validated against experiments for neutral PNIPAM microgels with c = 5% and hydrodynamic radius of ≈400 nm,⁴⁴ the shifts would correspond to an increase from T ≈ 32 °C for neutral microgels to T ≈34 °C for microgels with f = 0.05 and T ≈37.5 °C for f = 0.20, respectively. These specific values should be taken with care, since the α − T mapping has been validated for non-charged microgels only and may not hold in the ionic case. Regardless, the observed trend of the increase of T_VPT with increasing charge is in qualitative agreement with experiments.^9,60,61

Fig. 5 Swelling curves (radius of gyration R_gversus effective temperature α) for the microgels with explicit counterions (EC) with f = 0.05 and f = 0.20, as compared to the neutral microgel. Inset: swelling curves normalized to the value of R_g(α = 0) (left axis). Dashed lines report |d[R_g/R_g(α = 0)]/dα| (right axis), whose maximum corresponds to the VPT transition. These curves are arbitrarily shifted along the y-axis to improve visualization.

Additionally, we find that, while the microgel radius of gyration becomes larger with increasing charge for small values of α, for α ≳ α_VPT all microgels have the same size, indicating that the collapsed state does not depend on the presence of charges. This result has the interesting consequence that, upon rescaling R_g by its value at the maximally swollen state (α = 0), as shown in the inset of Fig. 5, the swelling ratio becomes larger as f increases. Since such a ratio has been previously adopted as a measure of the particle softness,^18,19 this suggests that more charged microgels are softer than less charged or neutral ones, in agreement with experimental findings,^9,60,62 at least when the ionization is not too high. Indeed, in these studies it has also been shown that at high pH values (corresponding to large charges fraction f) the complete collapse of the microgel is no longer observed, differently from the present numerical results. This may be due to the fact that in our model we neglect the interplay between the hydrophilic character of the co-polymer and its charge content, while the presence of charges at high T could alter the monomer–solvent interactions. Thus, a more accurate analysis based on a systematic comparison with experiments should be done in the future, in order to establish a correct mapping, for ionic microgels, among the real temperature T and the effective temperature α. This should also take into account the effect of salt, which plays a major role when pH is different from 7. In addition, some of the observed discrepancy could be due to the random charge distribution that we have considered, which could be not entirely realistic and should be compared with different choices. In this respect an experimental measurement of how the charges are distributed throughout the microgel, e.g., a quantification of the charge fraction located on the surface rather than in the interior,^63,64 would be a very valuable input to the simulations.

The fact that the collapsed state is the same in all investigated cases could be misleadingly taken as an indication that all (or most of the) counterions are expelled from the interior of the microgel. However, this turns out not to be the case, as it can be seen from the internal charge distributions of the microgel reported in Fig. 6. Specifically, the evolution of the charged monomers and counterions density profiles is separately shown for a few selected α values across the VPT in Fig. 6(a), which only contains results for f = 0.20. The behavior for f = 0.05 is qualitatively similar and thus not shown. We find that the profiles of the charged monomers and counterions closely follow each other at all studied values of α. This indicates a residual presence of counterions inside the microgels, which actually increases with α in order to balance the increase of monomer charge density in the collapsed core. The fact that the presence of counterions inside the microgels does not affect the size of the collapsed state also indirectly confirms that the choice of a small size for the counterions in our simulations is appropriate.

Fig. 6 (a) Density profiles for charged monomers (lines) and counterions (symbols) as a function of the distance from the center of mass of the microgel for f = 0.20 from the swollen (α = 0) to the collapsed (α = 1.20) states. Included are values just below (α = 0.74) and just above (α = 0.90) the VPT. All the curves are normalized to the average number of charged monomers 〈fN〉, calculated over all realizations of the network topology and of the charge distribution; (b) net charge density profile for the same values of α.

Looking at the profiles in Fig. 6(a) more closely, we find a small difference between counterions and charged monomers profiles upon increasing α and close to the surface of the microgels. This can be better visualized in Fig. 6(b), which reports the net-charge density profiles ρ_Q(r) (defined in Methods) at different α-values. We find that in the swollen state the charge density is statistically zero at all distances, except for the outer corona region, where it takes a tiny negative contribution. This is caused by the outermost counterions that are entropically driven to freely move around the simulation box, even far from the microgel. This situation persists below the VPT. However a significant change occurs close and above the VPT temperature. Indeed, under these conditions the microgel still maintains a rather neutral core, but in the corona the charge density abruptly increases, leading to the formation of a charged double layer. This trend is enhanced as α increases, signalling that there is a large charge imbalance at the surface of the microgels, where the counterions tend to accumulate. Such a phenomenon can be tentatively explained as follows. For small α, the structure of the microgel is swollen and counterions can be close to the charged monomers at any distance from the center of mass, still retaining a large freedom of moving inside the network. However, when α increases, the asymmetry among the interactions experienced by charged monomers and counterions come into play. On one hand, charged monomers interact with the additional solvophobic potential V_α which partially counteracts their electrostatic mutual repulsion. These contributions combined with the presence of the polymer network, which constraints their positions, induce the charged monomers density ρ_CH to steeply decay close to the surface of the microgel. On the other hand, counterions are able to gain translational entropy and at the same time to reduce their mutual repulsion by positioning themselves close to the surface in a more dispersed way. This implies a smoother decay of ρ_CI. The asymmetry between these two behaviors causes the formation of the above-mentioned double layer. Hence for α > α_VPT, the microgels acquire an effective charge and are surrounded by a small counterion cloud. This result highlights the non-trivial arrangements of counterions with respect to the microgel structure. It can also be of guidance in the treatment of ionic microgels at finite concentrations, particularly for large ones where similar crowding effects may take place, which may cause the re-organization of the counterions distribution within the microgel even under swollen conditions.^7,28

3.3 Comparison between the explicit and implicit models across the volume phase transition

3.3.1 Swelling curves and snapshots. In this section we compare the behavior of the microgels with explicit counterions with those modelled with a Debye–Hückel approach across the volume phase transition for f = 0.20. In order to make a meaningful comparison, data have been averaged over the same topologies and distributions of charged monomers.

We start by reporting the swelling curves of implicit and explicit models in Fig. 7. For Debye–Hückel simulations we have computed two different swelling curves. The first one is obtained using the effective Debye length λ_D*, assuming it to be constant for all values of α. For the second swelling curve we have calculated the effective Debye length for each value of α to take into account the change in counterions density. Since the latter increases as a consequence of shrinking, the resulting λ_D*(α) decreases upon increasing α. The two swelling curves are very similar to each other, with small differences only visible close to the VPT, indicating that the transition occurs slightly earlier for the varying λ_D*(α) with respect to the constant one. However, in both cases, α_VPT is found to be close to the neutral microgel result, and hence smaller than that of the explicit one. Interestingly, making λ_D*(α) to vary with α leads to R_g predictions that are even further away from the explicit counterions case than those observed with a constant λ_D, suggesting that such an approach is deeply flawed. Our findings demonstrate that the charged microgel with explicit counterions retains a much larger structure for all α ≲ α_VPT.

Fig. 7 (a) Swelling curves of microgels with implicit (DH) and explicit (EC) counterions for f = 0.20, as well as the corresponding one for neutral microgels. For DH, we report both results for fixed λ_D* = 1.30σ, that is the effective screening length calculated from the α = 0 microgel with explicit counterions, and for varying λ_D*(α), calculated for each value of α; (b) same as in (a) but with curves rescaled by the respective values of R_g(α_VPT) on the y-axis and α_VPT on the x-axis.

The bottom panel of Fig. 7 shows the same swelling curves rescaled along both axis with the respective values of α and R_g at the VPT, in order to analyze the shape of the swelling curve with respect to each other. We find that, for α < α_VPT, both curves relative to the implicit model coincide with that of neutral microgels, while the explicit model significantly differs. For α > α_VPT, on the contrary, neutral, EC and DH curves are all different and, even for the implicit model, we find that the shrinking ratio is slightly increased with respect to the neutral case, confirming that charged microgels are softer also when modelled with the Debye–Hückel approach.

In order to better understand the main differences between the different models as the solvophobicity increases, in Fig. 8 we report representative snapshots of the system across the VPT. Data for the microgel with explicit counterions (top row) are compared to the Debye–Hückel model (intermediate row) and to the neutral system (bottom row) at similar values of α/α_VPT. All snapshots refer to the very same underlying network topology, in order to clearly discriminate the effects of charges. In the swollen regime, the microgel conformations are comparable, but the increase in microgel size as we go from neutral to DH to EC model is evident. By contrast, in the fully collapsed regime all microgels look very similar to each other. The most dramatic differences between the three situations can be immediately visualized close to the VPT. Under these conditions, corresponding to the second and third columns of Fig. 8, in the presence of explicit counterions the microgel appears to be made of a core and of a rather inhomogeneous corona. In fact, the most external chains do not completely collapse even when α = α_VPT, as they form small clusters between themselves while, at the same time, remaining clearly distinct from the homogeneous dense core. It is only when α significantly exceeds α_VPT that they get slowly incorporated within the core. We stress that this behavior is completely different from the observation of globule-like domains during a quench from the swollen to the collapsed state that was reported in several simulation works on neutral microgels.^57,59,65 In that case, such domains were found only during a transient regime, quickly disappearing at long times. Here such behavior is found in the final equilibrium states at intermediate α-values, genuinely indicating the stable presence of such inhomogeneities in the microgel structure across the VPT. On the other hand, these features are absent (in the long-time regime) both for neutral microgels and for microgels with implicit charges, where the collapse of the microgel is clearly homogeneous across the VPT, independently of the value of λ_D. These results can be explained by the fact that, for implicit charges, the competition between the electrostatic repulsion and the solvophobic attraction just shifts the occurrence of the VPT to larger values of α, because a larger amount of attraction is needed to compensate the additional monomer–monomer repulsion. However, when counterions are explicitly included, they provide the system with additional degrees of freedom, thus being able to compensate the balance between attraction and repulsion even locally. This creates inhomogeneities in the charge distributions which significantly alter the microgels internal profiles, giving rise to a distinct core–corona pattern close to the VPT.

Fig. 8 Simulation snapshots for the same microgel topology with f = 0.20 for different models and effective temperatures. The top row shows the charged microgel with explicit counterions (EC) from left to right: in the swollen state (α = 0), just before and after the VPT (α = 0.74 and α = 0.90), and in the collapsed state (α = 1.20). The intermediate and bottom rows display corresponding states for the implicit model (DH) with λ_D* = 1.30σ (α = 0.30, 0.65, 0.74 and 1.20) and for the neutral microgel (α = 0, 0.56, 0.64 and 1.00), respectively. Green (blue) particles represent neutral (charged) beads. Explicit counterions are shown as smaller red spheres. All snapshots refer to equilibrium states, where the microgel radius of gyration fluctuates around a constant value.

We note that, as evident from Fig. 8, in our model the outer chains of the microgel are found in closed loops, rather than as dangling ends. This is a consequence of our assembly process. While we plan to address possible differences in the structure of microgels in the presence of dangling ends in future works, we expect that the observed features of the microgel collapse in presence of charges should not depend on this aspect, since both neutral and charged microgels are assembled in the same way.

3.3.2 Form factors. In order to better quantify the behavior observed in the snapshots, we report the form factors of the microgels in Fig. 9, again comparing explicit, implicit and neutral cases at different values of α across the VPT. We find evidence that the neutral and implicit cases are quite similar to each other, and both are compatible with the extended fuzzy sphere model, as shown in the ESI,† Fig. S4. Instead, microgels with explicit counterions display a very different behavior in many aspects. First of all, we find that the first peak of P(q) is much smaller in intensity than for the other two cases for the investigated values of α < α_VPT. Indeed, it tends to only shift in position without growing much in amplitude upon increasing α. However, focusing on intermediate q-values beyond the first peak, P(q) considerably increases in height, a feature that is absent for implicit and neutral microgels and that cannot be captured by a fuzzy-sphere-like model (see below). No secondary peaks are observed. In addition, the behavior of P(q) looks almost discontinuous at the VPT temperature, sharply increasing for α > α_VPT and, at the same time, developing additional peaks. As the microgel approaches the fully collapsed state, it becomes again possible to describe its form factor with the extended fuzzy sphere model.

Fig. 9 Form factors for (a) explicit counterions (EC), (b) Debye–Hückel (DH, λ_D = 1.30σ) and (c) neutral microgels across the VPT. Data for charged microgels refer to f = 0.20.

To better discuss the features of the form factors with explicit counterions, a zoom of the data is reported in Fig. 10. For α < α_VPT, where we cannot rely on a fuzzy-sphere-like model, P(q) displays two distinct behaviors after the first peak, both of which are compatible with power law dependences. The first regime occurs for 0.2 ≲ qσ ≲ 0.6, where P(q) ∼ q^−δ₁ with the exponent δ₁ being rather constant for α < α_VPT, i.e. δ₁ = 0.75 ± 0.05. These q-values correspond to length scales within the corona region of the microgel. At larger q the form factors exhibit a crossover to a second regime characterized by a different apparent power law. The position of the crossover, marked with vertical lines in Fig. 10, shifts from qσ ∼ 0.55 at α = 0 to qσ ∼ 0.65 at α = 0.74. For such second regime, a power law description of the data as P(q) ∼ q^−δ₂ gives an exponent δ₂ strongly dependent on α (from ∼1.2 at α = 0 up to ∼1.8 close to the VPT). The fact that a similar power-law dependence in the first q-regime seems to hold for swollen microgels up to the VPT suggests that the outer corona structure remains roughly constant for this range of temperatures. By contrast, the increase of the apparent exponent at larger q-values suggests that for smaller length scales the structure feels the effect of the underlying interactions, which modify the fractal properties of the network. However, at such large values of q, beyond 2σ⁻¹, the data suffer from finite-size effects (as it can be observed by the onset of a minimum, which precedes the occurrence of the model-dependent monomer–monomer peak at qσ ∼ 2π⁵⁷). We have thus limited our analysis here and in the following to the range qσ ≲ 2, in which we have attempted a few types of different fits, going beyond the power-law behavior which cannot be considered to be very reliable in such a limited range of q (changing by only a decade).

Fig. 10 Zoom of the form factors of Fig. 9a for microgels with explicit counterions (EC) for f = 0.20 and relative fits. Symbols are simulation data below (α = 0, 0.56, 0.76) and above (α = 0.90, 1.00, 1.40) the VPT. Below the VPT: full lines are fits with a modified Lorentzian (eqn (12)) for the q-range that goes from the first peak of P(q) up to the correspondingly coloured vertical line; dashed lines are fits with a different modified Lorentzian function in the range starting from the vertical line up to qσ ≲ 2. In all cases, for qσ ≳ 2 data are affected from finite-size effects. Above the VPT: lines are fits according to an extended fuzzy sphere model (eqn (11)) plus a modified Lorentzian (eqn (12)). Data sets for different α are arbitrarily shifted on the y-axis to improve readability.

Among the available models, we found that the modified Lorentzian defined in eqn (12) is able to separately describe both regimes for α < α_VPT, as shown in Fig. 10. Interestingly, the fractal exponents D₁ and D₂ extracted from the fits in the two regimes, reported in the ESI† (see Section S3), closely match the apparent power-law exponents described above. Thus, a roughly α-independent value of D₁ is found for small q, while a larger value of D₂ is obtained, which rapidly increases with α. These two parameters refer to the fractal dimensions of the correlated domains in the network over the corresponding ranges of length scales. They are coupled to two characteristic lengths, ξ₁ and ξ₂, which quantify the correlation lengths among such domains.⁵⁶ These lengths are both found to decrease with α, in agreement with expectations. Most importantly, we find in all studied cases that ξ₁ > ξ₂. This suggests that the behavior of the form factors in the swollen state and up to VPT is compatible with a network with different characteristic domains occurring in the corona and in the core region, respectively. Within the corona region (first regime), the correlation length is quite large, reflecting the few domains (visible in the snapshots of Fig. 8) that are quite far apart from each other. The fractal dimension of such environments is rather low and unaffected by changes in α, reflecting the fact that the corona remains clearly distinct from the core, up to the VPT and beyond. Instead, within the core region (second regime), the domains correlation length is much smaller and it rapidly decreases with α, while the fractal dimension is larger and increases with α, consistent with the shrinking of the core. The trend of these parameters in the second regime is consistent with that found for α > α_VPT, in which we use an extended fuzzy sphere model plus a modified Lorentzian. Here the core–corona distinction gets less and less pronounced, suggesting that we do not need two Lorenztian terms any longer to fit the data. Further discussion on the reliability of the fits and a comparison of the extracted fit parameters with the implicit and neutral microgel models is reported in the ESI.†

Interestingly, a qualitatively similar behaviour has been reported in an experimental study of PNIPAM microgels at high crosslinking concentration.⁶⁶ In this work the form factors at various temperatures have been measured finding that, at intermediate and high q-values, a similar two-step decay of P(q) was observed. The authors attribute this behavior to the presence of a core and a surrounding shell having different fractal dimensions and correlation length of the polymer mesh. However, the nature of inhomogeneities is different between these experiments and our simulations. Indeed, in the former case, they arise from a sensible difference in the concentration of crosslinkers among the core and the external shell, which also leads to a difference in the solvophobic properties of the network in the two regions. Instead, in our model they stem from the competition among attraction and electrostatic repulsion. Nevertheless, this proves that such inhomogeneities can be observed in scattering experiments, also with microgels size of the order of ∼1 μm. It would be desirable to assess this scenario also in ionic microgels in the near future. At present, to the best of our knowledge, experimental form factors are available for an extended range of wavevectors only in a few studies, which unfortunately do not analyse their variation as a function of the temperature.^7,67 In addition, other studies report P(q) only in a limited q-range^68,69 or for different types of microgels, obtained through microfluidic techniques⁸ and IPN microgels,⁷⁰ for which it is not possible to make a meaningful comparison with respect to our numerical model.

From all the evidence gathered in this part, we can conclude that, in our model, below the VPT the competition between the solvophobic attraction and the repulsive electrostatic interactions, screened by the counterions, gives rise to a rather inhomogeneous structure, characterized by two different regimes. This complicated behaviour cannot be interpreted with a fuzzy-sphere-like description, and a new type of model would be needed to describe form factors in the whole q-range, perhaps inspired by multi-shell models.^71,72 On the contrary, neutral microgels and those with implicit charges display a simpler behaviour, and a modified fuzzy-sphere model with a fractal Lorentzian is sufficient to describe the data. To assess these structural features it would be also interesting to compare the present numerical findings with experimental data in the near future.

3.3.3 Comparison at the same microgel size. Finally, we perform a comparison for microgels with the same R_g obtained with the three employed models (neutral, implicit and explicit) in order to compare differences arising in the structures when they are roughly of the same size. We thus select the values R_g ∼ 26σ, R_g ∼ 21σ and R_g ∼ 17σ for which the system is respectively below the VPT temperature, slightly above it and in the fully collapsed state. The monomer density profiles of the microgels under these conditions are reported in Fig. 11.

Fig. 11 (a–c) Monomer density profiles and (d–f) form factors for microgels with the same radius of gyration below (R_g ∼ 26σ) and slightly above (R_g ∼ 21σ) the VPT, and in the highly collapsed state (R_g ∼ 17σ). The corresponding values of α for the three cases are: (i) α = 0.74, 0.90 and 1.40 for the explicit counterions (EC) model with f = 0.20 (symbols and lines), (ii) α = 0.45, 0.74 and 1.40 for the Debye–Hückel model (DH) with λ_D = λ_D* and f = 0.20 (lines), (iii) α = 0.30,0.65 and 1.40 for neutral microgels (full symbols).

For small values of α, we find that the monomer density inside the core is larger in the explicit charged microgels than in neutral or implicit ones (Fig. 11(a and b)). This is different for what observed for the maximally swollen case (α = 0), shown in Fig. 2(a), where the monomer density in the core was much smaller for the microgels with explicit counterions. This was due to the fact that at α = 0, the size of the microgels was very different. When we compare the models at the same R_g instead, we see that the explicit counterion microgel has more monomers in the core and for large distances, while it is less dense in an intermediate range of distances. These features are maintained even when we cross the VPT, where the corona of the explicit case is much more extended than the neutral and implicit ones. Finally when the complete collapse is achieved, all the microgels have an identical density profile within numerical uncertainty (Fig. 11(c)). We confirm no differences occurring between neutral and Debye–Hückel microgels at comparable R_g in the investigated range of temperatures and of electrostatic parameters.

Similar plots for the form factors are reported in Fig. 11(d–f). Interestingly, despite the microgels having the same R_g, the explicit charged ones clearly show that the first peak of P(q) is shifted to larger q-values. This is because R_g tells us how broad the mass distribution is, and we see from the density profiles that the microgels with explicit counterions have a more extended corona. This is then compensated by a smaller and denser core to produce the same R_g of the other two models. We can infer that the first peak of the form factor is mainly affected by the extension of the core rather than by R_g. Indeed, comparing the three kinds of microgels at the same value of α, R_g is much greater for the explicit model, which thus needs a much higher intensity of the attractive force to reduce its volume. Because of the underlying inhomogeneities, this leads to a denser core, within which the screening is stronger and the attraction larger. In addition, we also confirm that the intermediate q behavior of P(q) is completely different for the microgel with explicit counterions, even slightly above the VPT. It is only for very large values of α (≳1.20) that the structure is the same for all types of microgels.

The present results further suggest that the Debye–Hückel model is always very different from the one with explicit charges and that there seems to be no way to reconcile the two approaches. Thus, one could ask whether there is some other way to modify the parameters of the Debye–Hückel model in order to resemble the features observed in the presence of counterions. To this aim, we would need to bypass the standard definition of the Debye length in eqn (6), which clearly overestimates the screening effects of the counterions present in the core. One possibility would be a phenomenological-like approach in which we consider the value of λ_D as the one yielding the same R_g of the microgel in the maximally swollen conditions (α = 0). From Fig. 1, we observe that this would be achieved with a much larger value of the Debye length with respect to the effective one, i.e., λ_D = 3.0σ. We have thus performed additional simulations (not shown) of the Debye–Hückel model for this value of λ_D as a function of α to try to assess whether in this case the implicit treatment of the charges can give rise to inhomogeneous effects such as with explicit counterions. However, we find that the microgel undergoes a microphase separation at large α and does not resemble at all the case with explicit counterions. We will thus address the case of microgels with very high charge fractions in future works, concluding for the present study that the Debye–Hückel model yields results that appear to be too similar to those of neutral microgels, pointing to the crucial role of counterions in a correct treatment of single-microgel properties.

4 Conclusions

In this work, we carried out extensive numerical investigations of charged microgels, focusing on the single-particle structure and swelling behavior across the volume phase transition. Extending a realistic assembly protocol that we recently put forward in ref. 43 and 44, we now additionally included the effects of charges in two different ways. On one hand, we employed an implicit model where screening effects of the counterions are included through a Debye–Hückel treatment where we varied both the amount of charged monomers and the Debye length. On the other hand, we performed simulations in the presence of explicit counterions, interacting with the charged monomers through a Coulomb repulsion. In both frameworks, we addressed the importance of having an underlying disordered network topology with the desired core–corona architecture, similar to that featured in microgels synthesised through the most common routes.

Our results are consistent with common expectations for the behavior of thermoresponsive microgels where charged co-monomers are included in the synthesis. In particular, we find that the size of the microgels in the swollen state increases with the fraction of charged monomers included in the network. Such an increase is also responsible for the occurrence of larger swelling ratios for more charged microgels, confirming the link between charge and softness.¹⁸ This is found for both explicit and implicit counterion modeling, thus being a robust feature of charged microgels. We also confirm that the VPT temperature shifts to larger values as the amount of charge increases, in agreement with experimental results.^8,9,60,61 However, some of our findings are less obvious than could be naively thought. First of all, while charged microgels are very different with respect to the neutral case below and at the VPT temperature, we find no difference among their collapsed structures, independently on the presence of charges and of the treatment of the counterions, suggesting that at sufficiently high temperatures they eventually reach a homogeneous spherical structure of the same size. Furthermore, specific considerations have been made possible by the use of explicit counterions. In particular, we find that the fully collapsed microgel in the explicit model is not at all free of counterions, which therefore are not expelled from the interior of the microgel upon deswelling. Instead, they are retained inside it in order to balance the increased charge density of the collapsed structure. Thus, counterions freely permeate and screen monomer charges at all temperatures, acting as neutralizers for the polymer network. Only close to the microgel surface we observe the onset of a non-neutral local charge, which manifests just above the VPT temperature in agreement with electrophoretic measurements.¹⁰ As a consequence, charged microgels behave as overall neutral objects up to around the VPT temperature, at least in the dilute regime. This property also reflects on the stability of finite concentrations suspensions. Indeed, it was reported⁷³ in low pH experiments (corresponding to f ≈ 0) close to the VPT that both PNIPAM and PNIPAM-co-AAc microgels tend to aggregate as a consequence of the change in the solvent quality. However, this does not happen for neutral and high pH (corresponding to f > 0), for which the suspension of ionic microgels is found to remain stable against aggregation.

We also compared the structure and the swelling of the microgel using the Debye–Hückel model, this being a much more convenient way to treat charges from the theoretical and numerical point of view. It turned out that a qualitative agreement between the explicit and the implicit approaches is unachievable, even using an effective Debye length that was calculated from the density of counterions computed in the explicit case. Our findings indicate that the Debye–Hückel approach is not able to reproduce many important effects that arise in the presence of charges, being mainly able to describe the average effect of screening of counterions onto charged beads over the polymer network. In particular, it fails to take into account the osmotic pressure of both inner counterions, acting in favour of the microgel swelling, and external ones, acting against the swelling. Remarkably, the structural features observed for the Debye–Hückel model are actually much more similar to those of neutral microgels than to the explicit counterions case. The most prominent difference can be noticed in the snapshots of Fig. 8, where the inhomogeneous core–corona structure is augmented by the presence of charges and counterions, a fact that is completely missing in the implicit representation. Strikingly, this reflects on the form factors of the microgels, which show a profile that is incompatible with the a fuzzy-sphere-like model, but rather display the onset of two distinct regimes, each of them compatible with a modified Lorentzian. Gaining a strong theoretical understanding of these intriguing findings will be the subject of future work. A potentially interesting perspective would be to combine our simulations with theoretical approaches^24,25 in order to provide some description of the data and perhaps to develop a modified Debye–Hückel approach, which could take into account the inhomogeneity of the microgel, assigning different values of λ_D to the core and to the corona, respectively. However, how to determine these values a priori (i.e., without estimating them with explicit simulations) remains an open question.

The understanding of the single-particle properties of co-polymerised microgels can be considered as a first step toward a better understanding of interpenetrated network microgels (IPN), wherein PNIPAM and PAAc are organized into two independent, interpenetrated networks, so that the responsiveness to temperature and to pH can be decoupled.⁷⁴ This particular kind of microgels has recently gathered a lot of interest because of their intriguing fragility behavior: as the amount of charges increase, these systems exhibit features of strong glass-formers, a rather unique example in soft matter.^18,19 A recent work⁷⁵ has put forward the idea that this behavior directly stems from charge effects, which also increase the softness of the particles, as confirmed in the present work. However, other effects could play an important role under these conditions, in particular the ability of microgels to deform, compress and interpenetrate. The understanding of these issues is very timely, having been recently address in super-resolution microscopy experiments combined with rheological measurements⁷⁶ as well as in simulations of a simple elastic model.⁷⁷ It would thus be very interesting to address further these issues in future works for both IPN and regular microgels.

Finally, our aim will be to transfer the knowledge from single-particle properties to many-body systems by developing appropriate coarse-grained effective potentials, still retaining the essential ingredients of the microgels, in order to be able to address their structural and dynamical behavior at various concentrations. Hence, by calculating the effective potential between two charged microgels we could validate and refine the effective approaches carried out in recent works on the assembly properties of charged microgels in bulk.¹²

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank J. Ruiz-Franco for valuable discussions. We acknowledge support from the European Research Council (ERC Consolidator Grant 681597, MIMIC).

References

1. A. Fernandez-Nieves, H. Wyss, J. Mattsson and D. A. Weitz, Microgel suspensions: fundamentals and applications, John Wiley & Sons, 2011.
2. P. J. Yunker, K. Chen, M. D. Gratale, M. A. Lohr, T. Still and A. Yodh, Rep. Prog. Phys., 2014, 77, 056601.
3. L. A. Lyon and A. Fernandez-Nieves, Annu. Rev. Phys. Chem., 2012, 63, 25–43.
4. J. Brijitta and P. Schurtenberger, Curr. Opin. Colloid Interface Sci., 2019, 40, 87–103.
5. M. Karg, A. Pich, T. Hellweg, T. Hoare, L. A. Lyon, J. J. Crassous, D. Suzuki, R. A. Gumerov, S. Schneider, I. I. Potemkin and W. Richtering, Langmuir, 2019, 35(19), 6231–6255.
6. R. Pelton and T. Hoare, Microgel Suspensions: Fundamentals and Applications, 2011, vol. 1, pp. 1–32.
7. S. Nöjd, P. Holmqvist, N. Boon, M. Obiols-Rabasa, P. S. Mohanty, R. Schweins and P. Schurtenberger, Soft Matter, 2018, 14, 4150–4159.
8. D. Rochette, B. Kent, A. Habicht and S. Seiffert, Colloid Polym. Sci., 2017, 295, 507–520.
9. D. Capriles-González, B. Sierra-Martìn, A. Fernández-Nieves and A. Fernández-Barbero, J. Phys. Chem. B, 2008, 112, 12195–12200.
10. D. Truzzolillo, S. Sennato, S. Sarti, S. Casciardi, C. Bazzoni and F. Bordi, Soft Matter, 2018, 14, 4110–4125.
11. S. Nöjd, P. S. Mohanty, P. Bagheri, A. Yethiraj and P. Schurtenberger, Soft Matter, 2013, 9, 9199–9207.
12. T. Colla, P. S. Mohanty, S. Nöjd, E. Bialik, A. Riede, P. Schurtenberger and C. N. Likos, ACS Nano, 2018, 12, 4321–4337.
13. G. M. Eichenbaum, P. F. Kiser, A. V. Dobrynin, S. A. Simon and D. Needham, Macromolecules, 1999, 32, 4867–4878.
14. C. E. Reese, A. V. Mikhonin, M. Kamenjicki, A. Tikhonov and S. A. Asher, J. Am. Chem. Soc., 2004, 126, 1493–1496.
15. A. Scotti, U. Gasser, E. S. Herman, M. Pelaez-Fernandez, J. Han, A. Menzel, L. A. Lyon and A. Fernández-Nieves, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 5576–5581.
16. U. Gasser, A. Scotti and A. Fernandez-Nieves, Phys. Rev. E, 2019, 99, 042602.
17. A. M. Howe, S. Desrousseaux, L. S. Lunel, J. Tavacoli, H. N. Yow and A. F. Routh, Adv. Colloid Interface Sci., 2009, 147, 124–131.
18. J. Mattsson, H. M. Wyss, A. Fernandez-Nieves, K. Miyazaki, Z. Hu, D. R. Reichman and D. A. Weitz, Nature, 2009, 462, 83.
19. V. Nigro, R. Angelini, M. Bertoldo, F. Bruni, M. A. Ricci and B. Ruzicka, Soft Matter, 2017, 13, 5185–5193.
20. P. G. De Gennes, Scaling concepts in polymer physics, Cornell University Press, 1979.
21. Y. Levin, A. Diehl, A. Fernández-Nieves and A. Fernández-Barbero, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 65, 036143.
22. A. R. Denton and Q. Tang, J. Chem. Phys., 2016, 145, 164901.
23. T. Colla, C. N. Likos and Y. Levin, J. Chem. Phys., 2014, 141, 234902.
24. A. Moncho-Jordá, J. Chem. Phys., 2013, 139, 064906.
25. A. Moncho-Jordá and J. Dzubiella, Phys. Chem. Chem. Phys., 2016, 18, 5372–5385.
26. A. Denton, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 011804.
27. C. N. Likos, Microgel Suspensions: Fundamentals and Applications, 2011, pp. 163–193.
28. T. J. Weyer and A. R. Denton, Soft Matter, 2018, 14, 4530–4540.
29. L. Rovigatti, N. Gnan, L. Tavagnacco, A. J. Moreno and E. Zaccarelli, Soft Matter, 2019, 15, 1108–1119.
30. A. Martìn-Molina and M. Quesada-Pérez, J. Mol. Liq., 2019, 280, 374–381.
31. H. Kobayashi and R. G. Winkler, Polymers, 2014, 6, 1602–1617.
32. H. Kobayashi, R. Halver, G. Sutmann and R. G. Winkler, Polymers, 2017, 9, 15.
33. S. Schneider and P. Linse, J. Phys. Chem. B, 2003, 107, 8030–8040.
34. G. C. Claudio, K. Kremer and C. Holm, J. Chem. Phys., 2009, 131, 094903.
35. M. Quesada-Pérez, J. Ramos, J. Forcada and A. Martìn-Molina, J. Chem. Phys., 2012, 136, 244903.
36. M. Quesada-Pérez and A. Martìn-Molina, Soft Matter, 2013, 9, 7086–7094.
37. P. K. Jha, J. W. Zwanikken, F. A. Detcheverry, J. J. De Pablo and M. O. De La Cruz, Soft Matter, 2011, 7, 5965–5975.
38. R. Schroeder, A. A. Rudov, L. A. Lyon, W. Richtering, A. Pich and I. I. Potemkin, Macromolecules, 2015, 48, 5914–5927.
39. L. G. Rizzi and Y. Levin, J. Chem. Phys., 2016, 144, 114903.
40. C. Hofzumahaus, P. Hebbeker and S. Schneider, Soft Matter, 2018, 14, 4087–4100.
41. D. Sean, J. Landsgesell and C. Holm, Gels, 2018, 4, 2.
42. M. Quesada-Pérez, J. A. Maroto-Centeno, A. Martìn-Molina and A. Moncho-Jordá, Phys. Rev. E, 2018, 97, 042608.
43. N. Gnan, L. Rovigatti, M. Bergman and E. Zaccarelli, Macromolecules, 2017, 50, 8777–8786.
44. A. Ninarello, J. J. Crassous, D. Paloli, F. Camerin, N. Gnan, L. Rovigatti, P. Schurtenberger and E. Zaccarelli, 2019, arXiv preprint arXiv:1901.11495.
45. G. S. Grest and K. Kremer, Phys. Rev. A: At., Mol., Opt. Phys., 1986, 33, 3628.
46. T. Soddemann, B. Dünweg and K. Kremer, Eur. Phys. J. E: Soft Matter Biol. Phys., 2001, 6, 409–419.
47. R. J. Hunter and L. R. White, Foundations of colloid science, Clarendon Press, 1987.
48. Y. Levin, Braz. J. Phys., 2004, 34, 1158–1176.
49. C. G. Lopez, A. Scotti, M. Brugnoni and W. Richtering, Macromol. Chem. Phys., 2019, 220, 1800421.
50. J. Landsgesell, L. Nová, O. Rud, F. Uhlk, D. Sean, P. Hebbeker, C. Holm and P. Košovan, Soft Matter, 2019, 15, 1155–1185.
51. M. Tuckerman, Statistical mechanics: theory and molecular simulation, Oxford University Press, 2010.
52. M. Deserno and C. Holm, J. Chem. Phys., 1998, 109, 7678–7693.
53. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
54. M. Stieger, W. Richtering, J. S. Pedersen and P. Lindner, J. Chem. Phys., 2004, 120, 6197–6206.
55. G. M. Conley, S. Nöjd, M. Braibanti, P. Schurtenberger and F. Scheffold, Colloids Surf., A, 2016, 499, 18–23.
56. M. Shibayama, T. Tanaka and C. C. Han, J. Chem. Phys., 1992, 97, 6829–6841.
57. F. Camerin, N. Gnan, L. Rovigatti and E. Zaccarelli, Sci. Rep., 2018, 8, 14426.
58. H. Kobayashi and R. G. Winkler, Sci. Rep., 2016, 6, 19836.
59. A. J. Moreno and F. Lo Verso, Soft Matter, 2018, 14, 7083–7096.
60. P. Holmqvist, P. Mohanty, G. Nägele, P. Schurtenberger and M. Heinen, Phys. Rev. Lett., 2012, 109, 048302.
61. X. Li, J. Zuo, Y. Guo, L. Cai, S. Tang and W. Yang, Polym. Int., 2007, 56, 968–975.
62. J. Pinheiro, L. Moura, R. Fokkink and J. Farinha, Langmuir, 2012, 28, 5802–5809.
63. C. Huang, H. Kobayashi, M. Moritaka and M. Okubo, Polym. Chem., 2017, 8, 6972–6980.
64. O. Virtanen, A. Mourran, P. Pinard and W. Richtering, Soft Matter, 2016, 12, 3919–3928.
65. S. Nikolov, A. Fernandez-Nieves and A. Alexeev, Appl. Math. Mech., 2018, 39, 47–62.
66. A. Fernandez-Barbero, A. Fernandez-Nieves, I. Grillo and E. Lopez-Cabarcos, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 66, 051803.
67. J. J. Lietor-Santos, B. Sierra-Martìn, U. Gasser and A. Fernández-Nieves, Soft Matter, 2011, 7, 6370–6374.
68. J. S. Hyatt, A. M. Douglas, C. Stanley, C. Do, T. H. Barker and A. Fernández-Nieves, Phys. Rev. E, 2017, 95, 012608.
69. G. Romeo, L. Imperiali, J.-W. Kim, A. Fernández-Nieves and D. A. Weitz, J. Chem. Phys., 2012, 136, 124905.
70. V. Nigro, R. Angelini, M. Bertoldo, F. Bruni, V. Castelvetro, M. A. Ricci, S. Rogers and B. Ruzicka, J. Chem. Phys., 2015, 143, 114904.
71. M. Bergman, PhD thesis, Lund University, 2019.
72. O. Peña-Rodrguez, P. P. González Pérez and U. Pal, Int. J. Spectrosc., 2011 DOI:10.1155/2011/583743.
73. N. Al-Manasir, K. Zhu, A.-L. Kjøniksen, K. D. Knudsen, G. Karlsson and B. Nystrom, J. Phys. Chem. B, 2009, 113, 11115–11123.
74. V. Nigro, R. Angelini, M. Bertoldo, V. Castelvetro, G. Ruocco and B. Ruzicka, J. Non-Cryst. Solids, 2015, 407, 361–366.
75. A.-M. Philippe, D. Truzzolillo, J. Galvan-Myoshi, P. Dieudonné-George, V. Trappe, L. Berthier and L. Cipelletti, Phys. Rev. E, 2018, 97, 040601.
76. G. M. Conley, C. Zhang, P. Aebischer, J. L. Harden and F. Scheffold, Nat. Commun., 2019, 10, 2436.
77. N. Gnan and E. Zaccarelli, Nat. Phys., 2019, 15, 683–688.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c9sm01253b

---