Charge affinity and solvent effects in numerical simulations of ionic microgels

Ionic microgel particles are intriguing systems in which the properties of thermo-responsive polymeric colloids are enriched by the presence of charged groups. In order to rationalize their properties and predict the behaviour of microgel suspensions, it is necessary to develop a coarse-graining strategy that starts from the accurate modelling of single particles. Here, we provide a numerical advancement of a recently-introduced model for charged co-polymerized microgels by improving the treatment of ionic groups in the polymer network. We investigate the thermoresponsive properties of the particles, in particular their swelling behaviour and structure, finding that, when charged groups are considered to be hydrophilic at all temperatures, highly charged microgels do not achieve a fully collapsed state, in favorable comparison to experiments. In addition, we explicitly include the solvent in the description and put forward a mapping between the solvophobic potential in the absence of the solvent and the monomer-solvent interactions in its presence, which is found to work very accurately for any charge fraction of the microgel. Our work paves the way for comparing single-particle properties and swelling behaviour of ionic microgels to experiments and to tackle the study of these charged soft particles at a liquid-liquid interface.


Introduction
Soft matter is a very active branch of condensed matter physics, which comprises, among other systems, colloidal suspensions, whose constituent particles can greatly vary in shape, softness and function. Soft matter encompasses not only synthetic particles, but also constituents of many biological systems, such as proteins, viruses and even cells, whose size ranges between the nano and the micrometer scale. A peculiar aspect of soft matter systems is the great variety of amorphous states they can form, including glasses [1,2] and gels [3,4,5]. Indeed, a large amount of work in this field is devoted to the study of these non-ergodic states which may form due to different kind of interactions, such as steric, hydrophobic or electrostatic ones, both of attractive and repulsive nature.
Sometimes, a single colloidal particle is already quite a complex object whose behaviour at the collective level is strongly connected to the microscopic features of the particle itself. This situation is typical of soft colloids, i.e. deformable particles with internal degrees of freedom strongly influencing their mutual interactions, which makes them already intrinsically multi-scale. For these systems a theoretical approach is quite challenging even at the single-particle level, thus it is convenient to rely on the development of suitable coarse-grained models [6] that allow to greatly reduce the system complexity, still capturing the important ingredients to be retained for a correct description of the collective behavior. This strategy is very profitable for the case of microgel particles [7] that, combining together the properties of colloids and polymers, can be viewed as a prototype example of soft particles [8,9]. A microgel is a microscale gel whose internal polymeric network controls its peculiar properties. By varying the constituent monomers, microgels can be made responsive to temperature, pH or to external forces [7]. For their intriguing properties, they are employed in a wide variety of applications, ranging from biomedical purposes [10,11] to paper restoration [12].
In order to be able to predict the behaviour of dense microgel suspensions and the formation of arrested states, it is important to properly take into account the internal degrees of freedom of the particles, by modelling their effective interactions in such a way that the resulting object can still shrink, deform and interpenetrate [13,14,15,16]. Hence, an accurate modelling of a single microgel is a necessary pre-requisite for a correct description of bulk suspensions. To validate numerical models at the single-particle level, there are a number of different experiments we can refer to. One of the most straightforward is the measurement of the effective size of the microgels via dynamic light scattering experiments. Upon varying the controlling parameter of the dispersion, the so-called swelling curves can be determined. For instance, microgels synthesized by employing a thermoresponsive polymer, such as Poly(N-isopropyl-acrylamide) (PNIPAM), undergo a so-called Volume Phase Transition (VPT) [7] at a temperature T VPT ≈ 32 • C from a swollen to a collapsed state. In addition, form factors can be measured by small-angle scattering experiments of dilute microgel suspensions, either using neutrons [17], x-rays [18] or even visible light for large enough microgels [19]. This observable directly provides information on the inner structure of the microgels and shows that microgels prepared via precipitation polymerization [20] can be modelled as effective fuzzy spheres [17], where a rather homogeneous core is surrounded by a fluffy corona, giving rise to what is usually called a core-corona structure. A more complex situation arises when ionic groups are added to the synthesis to make microgels responsive also to external electric fields [21,22] and to pH variations [23]. A case study of such these copolymerized microgels is the one made of PNIPAM and polyacrylic acid (PAAc) [20,24,25,26], that is pHresponsive due to the the weak acidic nature of AAc monomers.
An increasing amount of work in the last years has focused on modelling single-particle behaviour both of neutral [27,28,29] and ionic microgels [30,31,32,33]. For the latter case, we have recently shown [33] that it is important to take into account both the disordered nature of the network, as opposed to the diamond lattice structure [29], and an explicit treatment of charges and counterions. Indeed, meanfield approaches completely neglect the complex, heterogeneous nature of the charge distribution within the microgel.
In this work, we extend our previous effort by going one step further towards a realistic numerical treatment of ionic co-polymerized microgels.
In Ref. [33], we modelled a single microgel particle such that all of its monomers, including charged ones, experienced a solvophobic attraction on increasing temperature.
Here, instead, charged monomers experience Coulomb and steric interactions only. This is expected to be more realistic, since charged or polar groups always remain hydrophilic irrespectively of temperature, thus having a distinct behaviour with respect to all other monomers. We examine the consequences of this difference on the microgel swelling behaviour as well as on its structure and charge distributions across the VPT. In the second part of the manuscript, we consider the presence of an explicit solvent, to examine its effects on the structural properties of the microgel. In this way, we aim to make our model suitable for situations where solvent effects become important. In particular, this will enable us to study the effect of charges for microgels adsorbed at liquid-liquid interfaces, similarly to what we recently put forward for neutral microgels [34,35].

Methods
The coarse-grained microgels used in this work are prepared as in Refs. [27] starting from N patchy particles of diameter σ, which sets the unit of length, confined in a spherical cavity. A fraction c = 0.05 of these particles has four patches on their surface to mimic the typical crosslinker connectivity in a chemical synthesis, while all the others have two patches to represent monomers in a polymer chain. During the assembly, an additional force is employed on the crosslinking particles to increase their concentration in the core of the microgel in agreement with experimental features [18]. Once a fully-bonded configuration is reached (when the fraction of formed bonds is greater than 0.995), a permanent topology is obtained by enforcing the Kremer-Grest bead-spring model [36]. In this way, all particles experience a steric repulsion via the Weeks-Chandler-Anderson (WCA) potential, where sets the energy scale and r is the distance between the centers of two beads. Additionally, bonded particles interact via the Finitely Extensible Nonlinear Elastic potential (FENE), with kF = 15 which determines the stiffness of the bond and R0 = 1.5 which determines the maximum bond distance. The resulting microgel is thus constituted by a disordered polymer network with a core-corona structure and form factors across the VPT that closely resemble experimental ones for microgels synthesized via the precipitation polymerization procedure [18]. Once the microgels are prepared, we add electrostatic interactions to mimic co-polymerized polymer networks with charged groups. To this aim, we randomly assign a negative charge −e * to a given fraction f of microgel monomers, to mimic the acrylic acid dissociation in water, where e * = √ 4πε0εrσ is the reduced unit charge (which roughly corresponds to the elementary charge e, considering ≈ kBT at room temperature and σ as the polymer's Kuhn length), and ε0 and εr are the vacuum and relative dielectric constants. Accordingly, we insert in the simulation box an equal number of positively charged counterions with charge e * to ensure the neutrality of the system. Interactions among charged beads are given by the reduced Coulomb potential V coul (r) = qiqjσ e * 2 r , where qi and qj are the charges of counterions or charged monomers. We adopt the particle-particle-particle-mesh method [37] as a long-range solver for the Coulomb interactions. As discussed in a previous contribution [33], the size of the counterions is set to 0.1σ to facilitate their diffusion within the microgel network and to avoid spurious excluded volume effects. They interact with all other species simply via the WCA potential. The swelling behaviour of a thermoresponsive microgel can be reproduced in molecular dynamics simulations either by means of an implicit solvent, namely by adding a potential that mimics the effect of the temperature on the polymer, or by explicitly adding coarse-grained solvent particles within the box. In the first case, we employ a solvophobic potential with γ = π 9 4 − 2 1/3 −1 and β = 2π − 9 4 γ [38]. This potential introduces an effective attraction among polymer beads, modulated by the parameter α, whose increase corresponds to an increase in the temperature of the dispersion. For α = 0 no attraction is present, which corresponds to fully swollen, i.e. low-temperature, conditions. For neutral microgels, the VPT is encountered at α ≈ 0.65, while a full collapse is usually reached for α 1.2.
In the first part of this work, we analyze two different models, based on a different treatment of the interactions of charged ions on the microgels: in Model I all monomers experience a total interaction potential where Vα (Eq. 4) is added to the Kremer-Grest interactions, as previously done in Ref. [33]; in Model II only neutral monomers experience this total interaction potential, while the charged monomers do not interact with Vα, i.e. α = 0 for them in all cases. This second situation is equivalent to leaving unaltered the behaviour of charged groups of the microgel as the solvent conditions change, so that they always retain a good affinity for the solvent (α = 0). A similar treatment is also adopted for the counterions, for which α = 0 for both Model I and Model II.
In the second part of this work, we explicitly consider the presence of the solvent in driving the Volume Phase Transition. Solvent particles are modelled within the Dissipative Particle Dynamics (DPD) framework in order to avoid spurious effects which may arise from the use of a standard Lennard-Jones potential due to the excessive excluded volume of the solvent [39]. In the DPD scheme, two particles i and j experience a force where F C ij is a conservative repulsive force, with w(rij) = 1 − rij/rc for rij < rc and 0 elsewhere, F D ij and F R ij are a dissipative and a random contribution of the DPD, respectively; γ is a friction coefficient, θ is a Gaussian random variable with average 0 and unit variance, and ∆t is the integration time-step. We set rc = 1.75σ and γ = 2.0 in all the simulations. Here ai,j quantifies the repulsion between two particles i and j, which effectively allows the tuning of the monomer-solvent (m,s) and solvent-solvent (s,s) interactions. Following our previous work [39], we fix as,s = 25.0 while we vary am,s ≡ a between 5.0 and 16.0, that is the range where the collapse of a neutral microgel takes place. The reduced DPD density is set to ρsr 3 c = 3.9 (with ρs the actual number density of solvent beads). With this choice of parameters, we previously showed that this model reproduces the swelling behaviour and structural properties of a neutral microgel particle, in quantitative agreement with the implicit solvent model that was explicitly tested against experiments [18]. To compare the explicit solvent model with the implicit one, we only consider Model II, where charged monomers always retain a high affinity for the solvent. We will show later that, in the explicit treatment, this translates to having charged monomers-solvent interactions (ch,s) always set to a ch,s = 0. All other interactions are identical to the implicit solvent model. Simulations are performed with LAMMPS [40]. The equations of motion are integrated with the velocity-Verlet algorithm. All particles have unit mass m, the integration time-step is ∆t = 0.002 mσ 2 / and the reduced temperature T * = kBT / is set to 1.0 by means of a Nosè-Hoover thermostat for implicit solvent simulations or via the DPD thermostat for explicit solvent ones. In the former case the number of monomers in the microgels is fixed to N ≈ 42000, while for the latter case we limit the analysis to N ≈ 5000 due to the very large computational times owing to the presence of the solvent.
From equilibrated trajectories, we directly calculate the form factor of the microgel as, where rij is the distance between monomers i and j, while the angular brackets indicate an average over different configurations and over different orientations of the wavevector q (for each q we consider 300 distinct directions randomly chosen on a sphere of radius q). Also, we determine the radius of gyration Rg of the microgels as a measure of their swelling degree. This is calculated as where rCM is the position of the center of mass of the microgel. For each swelling curve, representing Rg as a function of the effective temperature α (implicit solvent) or a (explicit solvent), we define an effective VPT temperature, either αVPT or aVPT , as the inflection point of a cubic spline interpolating the simulation points. Finally, the radial density profile for all monomers is defined as By restricting the sum in Eq. 10 to only charged monomers or to counterions, we also calculate ρCH (r) and ρCI (r), that are the radial density profiles of charged microgel monomers and of counterions, respectively. In addition, we define the net charge density profile as ρQ(r) = −ρCH (r) + ρCI (r). (11) 3. Results and Discussion Here we focus on different swelling stages of the microgels upon increasing α. We notice immediately that a large amount of inhomogeneities persists in Model II at large α, in contrast with the behavior of Model I where a full collapse is achieved. This can be better quantified by the swelling curves, reporting the variation of the radius of gyration R g versus the effective temperature α, that are shown in Fig. 2 for different values of the charge fraction f . For both models we observe that the increase of f shifts the transition towards larger effective temperatures, but important differences arise at large α, as displayed in the snapshots. In Model I, where charged beads experience Coulomb as well as solvophobic interactions, the VPT occurs at all studied f , as shown in Fig. 2(a). Using the α-temperature mapping established in Ref. [18] through a comparison to experiments, the VPT temperature observed for f = 0.2 microgels would correspond to T ≈ 38 • C. However, experiments on ionic microgels, for which the amount of charges was systematically varied [41,26], have shown that even for values of f smaller than 0.2, the microgel does not collapse below 40 • C. As hypothesized in our earlier work [33], Model I neglects the interplay between the hydrophilic character of the co-polymer and its charge content. However, charged or polar groups, such as AAc groups, are known to remain hydrophilic even at high temperatures [42], which would increase the stability of the microgel in the swollen state with increasing f . We thus incorporate such a feature in Model II by removing solvophobic interactions for charged microgel beads. The resulting swelling curves, shown in Fig. 2(b), clearly demonstrate that for f = 0.20 the VPT is not encountered within the investigated solvophobicity range, up to values of α that would correspond to temperatures above 50 • C, in qualitative agreement with experimental observations [41,26,43].

Structural properties
It is now important to compare the two models from the structural point of view, to check whether major differences arise. We start the analysis by looking at the form factors which, in our previous work on Model I [33], were shown to exhibit novel features with respect to neutral microgels. In particular, we found evidence that for α < α VPT , the standard fuzzy-sphere-like model was not able to describe the numerical form factors, Instead, the emergence of two distinct power-law behaviours was found immediately after the first peak, at intermediate and high q values, respectively [33]. This was attributed to the presence of charges in the inhomogeneous structure of the microgel, which gives rise to different features for core and corona regions, each being characterized by a different domain size.
It is now crucial to verify whether such distinctive behaviour also persists when the interactions among charged beads are modelled more realistically.  neutral case (f = 0), at different values of α. For f = 0.032, the amount of charges in the microgel is still too low to observe differences between the two models and the neutral microgel. Also, at large α, the form factor is that of a collapsed microgel in all cases, as expected from the swelling curves in Fig. 2. For f = 0.2 and low enough α, the behaviour of the two models is again very similar, with the form factors of ionic microgels showing a first peak that is systematically smaller with respect to that of the  neutral case. At intermediate α, we find that two power-law-like behaviours are compatible with both sets of data for charged microgels, while the neutral case does not show such a feature. This finding, already elaborated in Ref. [33], appears to be a distinctive feature of our numerical model of ionic microgels and is the result of the combination of a random charge distribution within a disordered, heterogeneous network topology with the explicit treatment of ions and counterions. Such a distinctive feature was tentatively attributed to the different degree of swelling of the corona and of the core, but still awaits a direct experimental confirmation. However, hints of a similar two-step decay for P (q) were reported in Ref. [44] and would certainly deserve further investigation in future experiments. On the other hand, major differences between the two charged models arise for large values of α. Indeed, in Model I the microgel approaches and crosses the VPT leading to a fully collapsed state, while in Model II it remains in a quasi-swollen configuration for all studied α. Consequently, for high α values, the form factor does not resemble that of a homogeneous sphere, with only a second peak becoming evident, as opposed to the neutral case where many sharp peaks emerge. We notice that Model I fully coincides with the neutral case for very large α, even for f = 0.2.
In Fig. 4, we compare the monomers density profiles for the two models as a function of α. These data further indicate that, for Model II, the microgel does not achieve a collapsed state, as also visible from the behaviour of the profiles of charged monomers and of counterions, respectively. These are reported in the bottom panels of Fig. 4, showing that, for both models, the counterions are always found to be very close to the charged monomers, in order to neutralize the overall charge of the microgel. However, all profiles remain much more extended for Model II as compared to Model I, for all α. We stress that the comparison is performed for microgels with different affinity of the charged monomers for the solvent at the same α, which corresponds to very different swelling conditions, as evident from Fig. 2. Additional information can be extracted by comparing the two cases for a similar value of R g , as reported in Fig. 5. Also in this case, we find that Model II displays a more slowly decaying radial profile, albeit having a very similar mass distribution with respect to Model I, which is due to the presence of more stretched external dangling chains. Similar results also apply to ions and counterions profiles, that are shown in the inset of Fig. 5: even at the same R g , there is a surplus of charges at the surface in the case where the affinity of charged monomers for the solvent does not change with the effective temperature (Model II ). Overall, these findings confirm an enhanced stabilization of the swollen configuration operated by the charged groups of the microgel, hindering the tendency of the remaining (neutral) monomers to collapse.
To complete the structural analysis of the two models, it is instructive to consider the net charge density profile inside the microgels, that is reported in Fig. 6 for both f = 0.032 and f = 0.2. We confirm that, for both models, the net charge of the core region is roughly zero. However, it was shown in Ref. [33] that in the collapsed configuration a charged double layer arises at the surface of microgels, signalling the onset of a charge imbalance that grows with α. This feature, that is clearly visible in the behaviour of Model I at high α for all values of f , is also present for Model II for the low charge case (f = 0.032). However, the double peak in the net charge distribution is smeared out for f = 0.2, due to the fact that, up to the largest explored values of α, the microgel does not fully collapse. In this way, it maintains a low concentration of charged beads, that is always roughly balanced by counterions, resulting in a rather uniform charge profile. Instead, the peaks at the surface appear when the microgel collapses: this is indeed the case for both models at low charge fraction and even for large f when charged monomers are assigned a solvophobic behaviour (Model I ).
We conclude from this analysis that the hydrophilicity of the charged monomers at all effective temperatures enhances the tendency of the microgel to remain swollen, even when most of the monomers experience a very large solvophobic attraction. Thanks to the charge neutralization operated by counterions, the microgel remains very stable in a rather swollen configuration up to very large α, avoiding collapse for large enough values of f . This scenario agrees well with experimental observations, where the suppression of the VPT [26,43,42] is found when the concentration of charged hydrophilic groups in the polymer network is large enough. These considerations imply that Model I should not be used to describe microgels with high charge content. Indeed, its identical treatment of the solvophilic character of both neutral and charged monomers leads the particle to collapse at extremely high α. Incidentally, we report that this was observed also for unrealistic values of f up to 0.4 (not shown), in evident contrast with experiments. We will thus rely on Model II in the future to correctly incorporate charge effects in modelling microgels in a realistic fashion.

Solvent effects
We now go one step further in modelling ionic microgels, by explicitly adding the solvent to the simulations. This is a necessary prerequisite to tackle phenomena that cannot be described with an implicit solvent, e.g. situations in which hydrodynamics or surface tension effects at a liquid-liquid interface [34] play a fundamental role. In this subsection, we compare results for swelling behaviour and structural properties of the microgels for implicit and explicit solvent simulations. In particular, we restrict our discussion to Model II, having established this to be more in line with experimental observations. Since simulations with an explicit solvent require a much higher computational effort, we limit the following discussion to microgels with N ≈ 5000.

Swelling curves and explicit-implicit (a-α) mapping
We start by reporting the swelling curves of charged microgels, stressing the point that they have been obtained by fixing the value of a ch,s , which tunes the solvophilic properties of charged beads and counterions. We find that setting a ch,s = 0, while a m,s ≡ a varies, the explicit model is essentially equivalent to the implicit one. This means that it is possible to find a relation that links every implicit system with a certain value of the solvophobic attraction α to an explicit one with solvophobic parameter a that shows the same structure and swelling properties. In order to establish such a a-α mapping, we explored two different routes. The first one, referred to as linear mapping in the following, is based on the assumption that the dependence of a on α is linear, as previously adopted for neutral microgels [39]. In this way, the mapping relation is obtained through a horizontal rescaling of the relative swelling curves R imp g (α)/R imp g (α = 0) and R exp g (a)/R exp g (a = 0) for the neutral implicit and explicit microgels onto each other. Specifically, given two points for each curve, (a 1 , a 2 ) and (α 1 , α 2 ), the rescaled x-coordinate is calculated using the following relationship: where x = 0.5(x 1 + x 2 ) and ∆x = x 1 − x 2 with x = a, α. The second mapping a num (α), referred to as numerical mapping, has been obtained by numerically inverting the equation after spline fitting the two swelling curves. We report both mapping relations in Fig. 7(a), finding that they fall onto each other for almost the entire range of investigated solvophobic parameters in the two models, confirming the overall correctness of the assumption of linearity. However, we find some differences in the region α > 1.0 (a > 15). Having established the mapping for neutral microgels, we now use it to directly remap also the results for ionic microgels for all studied f without any further adjustments.

Swelling behaviour
The normalized swelling curves with varying charge fraction f , comparing implicit and explicit solvent, are reported in Fig. 7(be). Data from implicit simulations are mapped via both methods described above. For the neutral case, the presence of the solvent does not affect the swelling behaviour, as shown in Fig. 7(b), where no appreciable differences are found between linear and numerical mapping even at high α. Using the same relations for comparing charged microgels in explicit and implicit solvent, we find that, remarkably, the same swelling behaviour works for all charge contents. The swelling curves are virtually identical, which ensures that the inclusion of the solvent does not alter the microgel behavior in temperature even in the presence of charges. Small deviations, as expected from Fig. 7(a), appear only at large α values, being more pronounced for high charge content. This confirms the robustness of the DPD model which, as already discussed in Ref. [39], does not induce spurious effects, e.g. due to excluded volume, even in the collapsed state. An important result of this work is that, even in the presence of an explicit solvent, the microgel at high f does not fully collapse at large α, being entirely equivalent to implicit Model II and compatible with experimental findings.

Structural properties
In this subsection, we will show that the implicit and explicit solvent treatments with the newly established numerical mapping (Eq. 13) lead to identical structural features of the microgels. Small differences arise when using the linear mapping (Eq. 12) at high f and large values of α. We show in Fig. 8 the monomer (top panels), ion and counterion (middle panels) and charge (bottom panels) density profiles only for the f = 0.1 case, since similar results are also found for the other studied charge fractions. Reported data for different values of monomer-solvent interactions show an overall similarity between implicit and explicit solvent descriptions at all swelling conditions. Small deviations arise only for f = 0.2 for states with the highest values of a or α, when using the linear mapping: as we can observe from the rightmost panels of Fig. 8, the linear mapping fails to associate implicit and explicit states in the most collapsed state, where a visible difference arises between the profiles.
The distribution of ions and counterions within the microgel is an observable that should be more sensitive to the presence of the solvent. However, quite remarkably, also in this case, we find excellent agreement between the two models, as shown in the middle panels of Fig. 8. In particular, the emergence of a clear double-peak structure in the ion distribution is found in both models for large α (implicit) and a (explicit), signalling an accumulation of ions at the exterior surface of the microgels. This can be understood from the fact that ions, remaining always hydrophilic, never completely collapse onto the core of the particle. Thus, the appearance of a peak at distances corresponding to the outer region of the microgel is the result of an attempt of ions to maximize their contact with solvent. This is preceded by a minimum, which indicates a region where ions are  This feature is the echo of the minimum that arises in the net charge density distributions, already anticipated for the large microgel treated with the implicit model in Fig. 6. Importantly, a minimum also occurs in ρ Q (r) for smaller microgels, as shown in the bottom panels of Fig. 8, for the most collapsed conditions. Here a charged double layer is clearly present, with an excess of positive charges inside the microgel corona due to the increased amount of counterions in this region. At the same time, a negative charge surplus is found at the surface of the microgel, since charged ions preferably remain in contact with solvent particles. The net charge distribution is also identical for explicit and implicit solvent when using the numerical mapping, with again very small differences arising for the linear mapping at large α.
It is important to notice that, although a double layer was also observed with the implicit solvent in Ref. [33] (equivalent to Model I ), the two distributions (the one in Fig. 8 of the current manuscript and that reported in Fig. 6 of Ref. [33]) have opposite signs. Indeed in Ref. [33] the superposition of electrostatic and solvophobic effects led to an accumulation of counterions at the microgel surface, with the onset of a seemingly Donnan equilibrium [45]. Notwithstanding the different origin of the double layer, both models demonstrate that an almost perfect neutrality is achieved within the core of the microgel, and it is only at the surface that inhomogeneous distributions appear. Besides, the reduced size of the microgels studied with the explicit solvent facilitates the onset of peaks due to the increased surface-to-volume ratio of the microgels. A more precise assessment of size effects and a careful comparison to experiments will be the subject of future works.
Finally, the explicit solvent model allows us to quantify the amount of solvent that is located inside the microgel as temperature increases. This is illustrated in Fig. 9, where the solvent density profile ρ s (r) is reported for different values of a and all investigated charge fractions. These plots confirm the reduced tendency to collapse of charged microgels which retain a large amount of solvent within the network structure. No inhomogeneities within the microgel are in general observed. At large f and α some oscillations arise which may be due to reduced statistics. Finally, this study confirms that even at temperatures above the VPT there is quite a residual amount of solvent within the microgel, that is significantly enhanced by increasing the charge. These findings are in line with expectations [43,46], that are thus confirmed by our simulations.

Conclusions
In this work we report an extensive numerical study of single microgel particles, a prototype of soft colloids that is of great interest for the colloidal community, particularly for the formation of arrested states with tunable rheological properties [9], including glasses [47,48] and gels [49]. The use of different polymers within the microgel network allows to exploit responsiveness to different control parameters, such as temperature and pH, giving rise also to unusual responses in the fragility of the system [47,50,16]. In order to be able to model dense suspensions of these soft particles, we can rely on two possible strategies. On one hand, we can exploit highly coarsegrained models, such as the Hertzian one, which completely neglect the polymeric degrees of freedom of the particles and thus cannot reproduce the complex phenomenology observed in experiments in the gel or glassy regimes, such as shrinking, faceting and interpenetration [51,15]. On the other hand, we can try to model a single microgel in a realistic way, aiming to reproduce its structural properties and, from this, to build effective interactions which retain the polymeric features of the single particle.
Adopting the second strategy, the aim of this work is to improve the current numerical modelling of single ionic microgels with randomly distributed charged groups, aiming to describe PNIPAM-co-PAAc microgels across the Volume Phase Transition. In particular, we assess two different ways to model the interactions of the charged monomers belonging to the polymer network, either considering or not a solvophobic attraction that mimics their hydrophilic/hydrophobic interactions. We find that, as long as the charged groups maintain the same affinity for the solvent, the tendency of the microgel to remain in swollen conditions is enhanced even at high effective temperatures. Thus, for a charge fraction of f = 0.2 we find no evidence of the collapse of the microgel within the investigated range of our simulations, in agreement with experimental observations that are currently available [41,26,43,42]. This result is different from the case where charged beads also attract each other like neutral monomers upon increasing temperature, which undergoes a Volume Phase Transition to a fully collapsed state [33]. Despite this fundamental difference, the structural properties of the microgels treated with both models are rather similar, especially at low and intermediate temperatures. For instance, we confirm that the peculiar power-law regimes observed in the form factors are independent of the chosen model.
Having established the most appropriate modelling for charged monomers, we then performed another necessary step in the modelling of ionic microgels, namely to explicitly consider the presence of the solvent, which may affect the rearrangement of the charges during the swelling-deswelling transition. To this aim, we build on previous results showing that for neutral microgels a description with an explicit solvent can be directly and quantitatively superimposed to the implicit modelling by using a DPD representation of the solvent, leaving unchanged the treatment of the polymer network with a bead-spring model. In this way, the solvophobic potential in Eq. 4, modulated by the parameter α, is replaced by the DPD repulsive interactions between monomers and solvent. The latter is varied through a change of the parameter a controlling the repulsion between non-charged beads, while the interaction between charged monomers always retains a solvophilic nature.
We have thus carried out a careful comparison between explicit and implicit solvent treatments, finding quantitative agreement between the two. Interestingly, the relation among a and α established by the comparison of neutral microgels can be used also to compare charged microgels, even with large values of f (some deviations occur only at f = 0.2 and large α, a values), where the same correspondence between implicit and explicit solvent states is retrieved. We showed that a linear mapping between the two control parameters of the interactions in the implicit and explicit case is sufficient to obtain a very good agreement between the two descriptions.
From our analysis of the internal structure of the microgels across the VPT, we found that counterions have a rather similar distribution within the microgel core, effectively neutralizing the internal charge at small distances, but being in excess close to the surface. This gives rise to a charged double layer for large values of a and α. Interestingly, such peaks in the charge density distributions are swapped with respect to the case of Model I, where ions do not experience a tendency to remain at the surface, since they are also treated as solvophobic. These detailed predictions will have to be compared to experiments on ionic microgels as a function of charge fraction, pH and T , in order to establish the limit of validity of our model and to further improve it, towards a more realistic description of experimental microgels.
In perspective, this work paves the way to study realistic charged microgels in diffusing conditions, such as in electrophoresis and thermophoresis experiments [52], or at liquid-liquid interfaces and to calculate their effective interactions, similarly to what has been done for neutral microgels [34,35]. In this way, we will be able to determine the conditions under which electrostatic effects play a dominant role over elastic ones. Another important line of research will be the assessment of the role of the network topology: examples of interesting cases whose properties could be investigated are microgels consisting of two interpenetrated networks [50,53] or ultra-low crosslinked [54,55] and hollow [56,57] microgels. Finally, we hope that our theoretical efforts will stimulate further experimental activity on charged microgels to verify the predicted behaviour so that it will be possible to tackle the investigation of dense suspensions in the near future.