Numerical insights on ionic microgels: structure and swelling behaviour

TL;DR

This study uses numerical simulations to analyze how charges and counterions influence the structure and swelling behavior of ionic microgels, highlighting the importance of explicit counterion modeling for accuracy.

Contribution

It extends numerical microgel modeling to ionic systems with explicit counterions, revealing significant structural and transition behavior differences from implicit models.

Findings

Charges significantly alter microgel structure near the VPT.

Explicit counterion modeling aligns better with experimental data.

VPT temperature shifts with increased charged monomers.

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\"uckel 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.

Equations

$$V_{\text{WCA}}(r)=\{\begin{array}{ll}4ϵ\left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right]+ϵ& \text{if}r\le 2^{1/6}\sigma \\ 0& \text{if}r>2^{1/6}\sigma \end{array}$$

$$V_{\text{FENE}}(r)=-ϵk_F{l_0}^2\log \left[1-{\left(\frac{r}{l_0\sigma }\right)}^2\right],r<l_0\sigma$$

$$V_{\alpha }(r)=\{\begin{array}{ll}-ϵ\alpha & \text{if}r\le 2^{1/6}\sigma \\ \frac{1}{2}\alpha ϵ\left\{\cos \left[\gamma {\left(\frac{r}{\sigma }\right)}^2+\beta \right]-1\right\}& \text{if}{2}^{1/6}\sigma <r\le l_0\sigma \\ 0& \text{if}r>l_0\sigma \end{array}$$

$$V_{DH}(r)$$

$${\lambda }_B=\frac{z^2{e}^2}{4\pi {ϵ}_0{ϵ}_r{k}_{B}T}$$

$${\lambda }_D={\left(4\pi {\lambda }_B{\rho }_{ci}\right)}^{-1/2}.$$

$$V_{\mathrm{coul}}(r)$$

$$R_g={\left(\frac{{\sum }_{i=1}^N({\vec{r}}_i-{\vec{r}}_{CM}{)}^2}{N}\right)}^{1/2}$$

$$\rho (r)=⟨\frac{{\sum }_{i_{=1}}^N\delta (|{\vec{r}}_i-{\vec{r}}_{CM}|-r)}{N}⟩.$$

$$P(\vec{q})=⟨\frac{1}{N}\sum _{i=1}^N\sum _{j=1}^N\exp \left[i\vec{q}\cdot ({\vec{r}}_i-{\vec{r}}_j)\right]⟩,$$

$$P(q)$$

$$P(q)\propto \frac{1}{{\left[1+\frac{D+1}{3}{\xi }^2{q}^2\right]}^{D/2}}$$

Full text

1 Introduction

Microgels are colloidal-scale polymer networks which have recently become a favourite model system1, 2, 3, 4 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) 5 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 the ability to reversibly increase or reduce their size after a change of temperature around the so-called VPT temperature $T_{\rm VPT}\sim 32^{\circ}$ 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 \chH+ or \chOH- ions due to the dissociation of a fraction of monomers 6.

Out of the many possibilities provided by modern-day synthesis methods, co-polymerised PNIPAM-co-PAAc microgels are of particular interest 6, 7, 8, 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 \chH^+, 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 distribution of these 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 salt5, 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 tuned11, 12. Their single-particle properties have been extensively investigated in experiments as a function of both temperature and pH9. 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 monomers8, 7. The tunability of ionic microgels has also been exploited in several fields of research, from biology13 to laser technology14. For instance, their dual responsiveness makes them highly suitable to be employed in the smart design of switch optical devices 2 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 (often considered as neutral microgels) show interesting features, particularly above the VPT temperature, which suggest the presence of charge effects15, 10. 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 glasses16, 17.

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 theory18, 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 equation19, the Ornstein-Zernike integral equation20 and density functional theory21. Also an effective interaction potential has been derived using linear response theory22, which made it possible to draw a phase diagram for ionic microgels23. On the numerical side, the use of coarse-grained models12, 24 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 recently25, 26 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 27, 28. 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 29, while simulations of single nanogel particles have appeared only later on 30, 31, 32, 33, 34. Several techniques have been devised to treat charged networks. Particularly, recent Monte Carlo simulations35, 36, 37 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 observations8. 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 networks25.

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 works38, 39, 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 crosslinkers35, 36, 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 varies38, 39.

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, differently from neutral ones, already takes place at concentrations below the overlap one7, 24.

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 protocol38 to generate disordered, heterogeneous microgels that are structurally similar to real neutral ones. We start by preparing fully connected spherically shaped networks by confining patchy particles in a cavity with a designing force on the crosslinkers, which provides the typical core-corona structure of realistic microgels39. Once the network is formed, we freeze the topology of the network and adopt a monomer-resolved approach40. The beads that make up the polymers interact via a steric repusion, modeled with the Weeks-Chandler-Anderson (WCA) potential:

[TABLE]

where $\epsilon$ and $\sigma$ are respectively the energy and length units. In addition, connected beads interact via the finitely extensible nonlinear elastic potential (FENE):

[TABLE]

where $l_{0}$ 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_{0}=1.5$.

All monomers also interact with each other by means of an effective solvophobic potential, named $V_{\alpha}$, which implicitly takes into account the monomer-solvent interactions41:

[TABLE]

with $\gamma=\pi\left(\frac{9}{4}-2^{1/3}\right)^{-1}$ and $\beta=2\pi-\frac{9}{4}\gamma$41. This potential represents an effective attraction, modulated by the solvophobic parameter $\alpha$, arising between thermo-responsive monomers at high temperatures. In other words, $\alpha$ plays the role of an effective temperature: $\alpha=0$ represents good solvent conditions, while with increasing $\alpha$ the quality of the solvent worsens, leading to the aggregation of beads and to the shrinking of microgel particles.

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 42:

[TABLE]

where $\lambda_{\scriptscriptstyle B}$ and $\lambda_{\scriptscriptstyle 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_{\scriptscriptstyle B}T$, thus quantifying the relative intensity of the electrostatic forces, and it is defined as:

[TABLE]

where $\epsilon_{0}$ and $\epsilon_{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 $\lambda_{\scriptscriptstyle B}$ and on the density of counterions $\rho_{ci}$ as:

[TABLE]

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 43 since the multiple dissociation of single polyelectrolyte monomers is highly unlikely. We work with reduced units, with $\sigma$, $m$, $\epsilon$ being the units of length, mass and energy, respectively. Within this unit system, the experimental Bjerrum length, that is $\lambda_{\scriptscriptstyle B}\approx 0.7nm$ for monovalent ions in water at room temperature, translates into a reduced Bjerrum length $\lambda_{\scriptscriptstyle B}^{*}\approx 1$, assuming $\sigma\sim 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 $\sigma$, according to different types of measurements for linear PNIPAM chains44. Note that the use of a larger value of $\sigma$ would significantly decrease $\lambda_{\scriptscriptstyle 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 we cannot easily compare the gyration radius as a function of $\lambda_{\scriptscriptstyle D}$ with the experimentally measured diameter of the microgel varying the pH of the suspension or the salt concentration. Indeed, for weak poly-electrolytes, there is not a simple link between the pH and the dissociation fraction of the acidic monomers, which determines the value of $\lambda_{\scriptscriptstyle D}$ 45. 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:

[TABLE]

For ion-ion interactions this term is complemented by a steric repulsion, modeled again with the WCA potential (Eq. 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 Supporting Information. In particular, we analyzed the dependence of our results on the choice of the simulation box (see section S1), 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 $\sigma_{c}$ on the microgel swelling behavior (see section S2), finding that the use of too large counterions yields unrealistic excluded volume effects in the collapsed state of the microgel. We thus fix $\sigma_{c}=0.1\sigma$ throughout the rest of the manuscript.

2.2 Numerical simulations

We perform Molecular Dynamics simulations of single microgels with $N\sim 42000$ monomers at fixed crosslinker concentration $c=5\%$. Microgels are assembled as in Ref. 39 in a spherical cavity of radius $R_{0}=50\sigma$, 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 randomly assign a charge $ze=1$ to a fraction $f$ of the monomers, maintaining this charge distribution fixed throughout the simulation run. We average results over four different topologies and three different charge distributions.

We study microgels for three different charge fractions, $f=0.05$, $0.20$, $0.95$, and for several values of the solvophobic parameter $\alpha$ across the VPT. The equations of motion of the system are integrated via the velocity-Verlet algorithm 46. The equilibration of the system is carried out in the canonical ensemble using the Nosè-Hoover chains thermostat for $1.6\times 10^{6}$ simulation timesteps, while a long production run in the microcanonical ensemble of $\sim 2\times 10^{6}$ steps is used to obtain equilibrium averages of the thermodynamic observables under investigation. We used a cut-off of $R_{\scriptscriptstyle cut}=5\lambda_{\scriptscriptstyle D}$ for the Debye-Hückel potential, whereas the long-range Coulomb interactions are computed with the particle-particle-particle-mesh method 47. For the latter type of simulations we used the LAMMPS package48.

2.3 Main Observables

To assess the microgel size, we calculate the radius of gyration, defined as:

[TABLE]

where $\vec{r}_{i}$ and $\vec{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 density profile, defined as the average density at a fixed distance from the center of mass:

[TABLE]

We also compute the density profile of charged monomers, labelled as $\rho_{CH}(r)$, and that of counterions only, labelled as $\rho_{CI}(r)$. By adding the two latter quantities, weighted by the respective charge, we obtain the net-charge density profile $\rho_{Q}(r)=-\rho_{CH}(r)+\rho_{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:

[TABLE]

where the brackets $\langle\cdot\rangle$ indicate ensemble averages and $r_{i}$ is the position of the $i$-th monomer. We have computed the rotationally invariant quantity $P(q)$ as an average of $P(\vec{q})$ over 300 vectors $\vec{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 model49, 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 $\rho(r)\propto\operatorname{Erfc}\left(\frac{r-R}{\sqrt{2}\sigma_{\rm surf}}\right)$, where $\operatorname{Erfc}(\cdot)$ is the complementary error function, while $R$ and $\sigma_{\rm surf}$ are related to the extension of the core and of the corona, respectively. However, it has recently been shown by super-resolution microscopy measurements 50 that the assumption of a homogeneous core is not accurate, and that the density profile of microgels is better approximated by the function $\rho(r)\propto\operatorname{Erfc}\left(\frac{r-R}{\sqrt{2}\sigma_{\rm surf}}\right)(1-sr)$, 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 39. The additional linear term in the density profiles modifies the shape of the form factor in an extended fuzzy sphere model:

[TABLE]

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 hydrogels51. A step forward is represented by the modified Lorentzian proposed by Shibayama and Tanaka51, 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:

[TABLE]

with $\xi$ 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 $\alpha=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 $\lambda_{\scriptscriptstyle 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 $\lambda_{\scriptscriptstyle 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 Refs.27, 28, 35, 36, 37 for a diamond-like microgel and we find comparable variation of the microgel size to that reported in these works.

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 $\lambda_{\scriptscriptstyle 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 topologies38, 52.

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 larger number of charges. Interestingly, we can compare the results in Fig. 2, with those reported in Ref.27 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 microgels25. 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 monomers45. 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\sigma\sim 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.

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 $\lambda_{\scriptscriptstyle D}$.

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 and a less dense core than the Debye-Hückel model for all investigated values of the Debye length. Even a large increase of $\lambda_{\scriptscriptstyle 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 $\lambda_{\scriptscriptstyle D}^{*}$ from Eq. 6 by substituting $\rho_{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 $\lambda_{\scriptscriptstyle D}^{*}\simeq 2.4\sigma$ for $f=0.05$ and $\lambda_{\scriptscriptstyle D}^{*}\simeq 1.3\sigma$ for $f=0.2$, respectively. The resulting density profiles of the $\lambda_{\scriptscriptstyle 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 $\rho_{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$ ($\rho_{Q}\sim 1.5\cdot 10^{-2}\sigma^{-3}$) and $f=0.20$ ($\rho_{Q}\sim 5.0\cdot 10^{-2}\sigma^{-3}$). Thus, the counterions are able to freely diffuse throughout the microgels, even within the core, so that they fully counteract the electrostatic repulsion. Tthe presence of the counterions inside the network thus contributes to the increase of 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 $\lambda_{\scriptscriptstyle 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.

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 $\lambda_{\scriptscriptstyle D}=\lambda_{\scriptscriptstyle 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 53. 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_{\alpha}$ between monomers (Eq. 3) to mimic the increase of temperature in experiments9, 8.

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 $\alpha$ 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 $\alpha$ at which the VPT occurs, i.e. $\alpha_{VPT}$, defined as the position of the maximum of $|dR_{g}/d\alpha|$, shifts from $\alpha_{VPT}\sim 0.63$ found in neutral microgels 54, 39, to $\alpha_{VPT}\sim 0.69$ for $f=0.05$ and up to $\alpha_{VPT}\sim 0.82$ for $f=0.20$, as reported in the inset of Fig. 5. Using the $\alpha-T$ mapping validated against experiments for neutral PNIPAM microgels with $c=5\%$ crosslinkers and hydrodynamic radius of $\approx 400$nm39, the shifts would correspond to an increase from $\sim 32^{\circ}$C for neutral microgels to $\approx 34^{\circ}$C for microgels with $f=0.05$ and $\approx 37.5^{\circ}$C for $f=0.20$, respectively. These specific values should be taken with care, since the $\alpha-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 experiments9, 55.

Additionally, we find that, while the microgel radius of gyration becomes larger with increasing charge for small values of $\alpha$, for $\alpha\gtrsim\alpha_{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 ($\alpha=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 softness16, 17, this suggests that more charged microgels are softer than less charged or neutral ones, in agreement with experimental findings.

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 $\alpha$ 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 not shown. We find that the profiles of the charged monomers and counterions closely follow each other at all studied values of $\alpha$. This indicates a residual presence of counterions inside the microgels, which actually increases with $\alpha$ 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.

Looking at the profiles in Fig. 6(a) a bit more closely, we find a small difference between counterions and charged monomers profiles upon increasing $\alpha$ and close to the surface of the microgels. This can be better visualized in Fig. 6(b), which reports the net-charge density profiles $\rho_{Q}(r)$ (defined in Methods) at different $\alpha$-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 $\alpha$ 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 $\alpha$, 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 $\alpha$ 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_{\alpha}$ 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 $\rho_{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 $\rho_{CI}$. The asymmetry between these two behaviors causes the formation of the above-mentioned double layer. Hence for $\alpha>\alpha_{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 conditions7, 24.

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 $\lambda_{\scriptscriptstyle D}^{*}$, assuming it to be constant for all values of $\alpha$. For the second swelling curve we have calculated the effective Debye length for each value of $\alpha$ to take into account the change in counterions density. Since the latter increases as a consequence of shrinking, the resulting $\lambda_{\scriptscriptstyle D}^{*}(\alpha)$ decreases upon increasing $\alpha$. 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 $\lambda_{\scriptscriptstyle D}^{*}(\alpha)$ case with respect to the constant one. However, in both cases, $\alpha_{VPT}$ is found to be close to the neutral microgel result, and hence smaller than that of the explicit case. Interestingly, the $\lambda_{\scriptscriptstyle D}^{*}(\alpha)$ case leads to $R_{g}$ predictions that are even further away from the explicit counterions case than those observed with a constant $\lambda_{\scriptscriptstyle 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 $\alpha\lesssim\alpha_{VPT}$.

The bottom panel of Fig. 7 shows the same swelling curves rescaled along both axis with the respective values of $\alpha$ 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 $\alpha<\alpha_{VPT}$, both curves relative to the implicit model coincide with that of neutral microgels, while the explicit model significantly differs. For $\alpha>\alpha_{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 $\alpha/\alpha_{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 $\alpha=\alpha_{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 $\alpha$ significantly exceeds $\alpha_{VPT}$ that they get slowly incorporated within the core. This behavior is completely absent 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 $\lambda_{\scriptscriptstyle D}$. These findings 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 $\alpha$, 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.

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 $\alpha$ 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 Supporting Information, Fig. S3. 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 $\alpha<\alpha_{VPT}$. Indeed, it tends to only shift in position without growing much in amplitude upon increasing $\alpha$. 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 $\alpha>\alpha_{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.

To better discuss the features of the form factors with explicit counterions, a zoom of the data is reported in Fig. 10. For $\alpha<\alpha_{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\lesssim q\sigma\lesssim 0.6$, where $P(q)\sim q^{-\delta_{1}}$ with the exponent $\delta_{1}$ being rather constant for $\alpha<\alpha_{VPT}$, i.e. $\delta_{1}=0.75\pm 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\sigma\sim 0.55$ at $\alpha=0$ to $q\sigma\sim 0.65$ at $\alpha=0.74$. For such second regime, a power law description of the data as $P(q)\sim q^{-\delta_{2}}$ gives an exponent $\delta_{2}$ strongly dependent on $\alpha$ (from $\sim 1.2$ at $\alpha=0$ up to $\sim 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\sigma^{-1}$, 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\sigma\sim 2\pi$52). We have thus limited our analysis here and in the following to the range $q\sigma\lesssim 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).

Among the available models, we found that the modified Lorentzian defined in Eq. 12 is able to separately describe both regimes for $\alpha<\alpha_{VPT}$, as shown in Fig. 10. Interestingly, the fractal exponents $D_{1}$ and $D_{2}$ extracted from the fits in the two regimes, reported in the Supporting Information (see Section  S3), closely match the apparent power-law exponents described above. Thus, a roughly $\alpha$-independent value of $D_{1}$ is found for small $q$, while a larger value of $D_{2}$ is obtained, which rapidly increases with $\alpha$. 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, $\xi_{1}$ and $\xi_{2}$, which quantify the correlation lengths among such domains51. These lengths are both found to decrease with $\alpha$, in agreement with expectations. Most importantly, we find in all studied cases that $\xi_{1}>\xi_{2}$. 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 $\alpha$, 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 $\alpha$, while the fractal dimension is larger and increases with $\alpha$, consistent with the shrinking of the core. The trend of these parameters in the second regime is consistent with that found for $\alpha>\alpha_{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 Supporting Information.

It is interesting to compare our findings to the few available small-angle neutron scattering measurements that we are aware of, namely results for large PNIPAM-co-PAAc microgels8 and for IPN microgels56. In both cases, the probed range of wavevectors is limited to a portion of the interior of the network, which lies well within the core region. Unfortunately, our simulated microgels are too small to make a direct comparison over a sizeable range of wave-vectors. However, in qualitative terms, both sets of experiments report an increase of the signal with $T$ for smaller values of $q$, while at larger values of the wavevectors a crossing of the data is found, so that the signal actually decreases in $T$. This seems to agree well with our numerical findings. Moreover, experimentally no peaks are detected in the probed $q$-range, and the data seems to be compatible with a power-law behavior. It remains a challenge for the future to compare our numerical form factors to light or x-ray scattering measurements that would be able to probe the corona region to confirm the lack of peaks and the power-law-like behaviors that we observe.

From all the evidence gathered in this part, we can conclude that 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 models57, 58. 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.

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 of roughly the same size. We thus select the values $R_{g}\sim 26\sigma$, $R_{g}\sim 21\sigma$ and $R_{g}\sim 17\sigma$ 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.

For small values of $\alpha$, 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,b)). This is different for what observed for the maximally swollen case ($\alpha=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 $\alpha=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 at all 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 $\alpha$, $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 $\alpha$ ($\gtrsim 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 Eq. 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 $\lambda_{\scriptscriptstyle D}$ as the one yielding the same $R_{g}$ of the microgel in the maximally swollen conditions ($\alpha=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. $\lambda_{\scriptscriptstyle D}=3.0\sigma$. We have thus performed additional simulations (not shown) of the Debye-Hückel model for this value of $\lambda_{\scriptscriptstyle D}$ as a function of $\alpha$ 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 $\alpha$ 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 and perspectives

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 Refs.38, 39, 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 in 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 softness16. 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 results9, 8, 55. 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 measurements10. As a consequence, charged microgels behave as overall neutral objects up to around the VPT temperature, at least in the dilute regime.

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 obtained within 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 approaches20, 21 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 $\lambda_{\scriptscriptstyle 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 decoupled59. 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 16, 17. A recent work60 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. It would thus be very interesting to address the behavior of IPN microgels in future works.

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 bulk12.

Acknowledgments

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

Supporting Information

Assessment of the size effects by varying the side of the simulation box, preliminary analysis on the choice of the counterion size, discussion on the form factors fits and corresponding parameters for the models analyzed in the main text.