Numerical modelling of non-ionic microgels: an overview

Microgels are complex macromolecules. These colloid-sized polymer networks possess internal degrees of freedom and, depending on the polymer(s) they are made of, can acquire a responsiveness to variations of the environment (temperature, pH, salt concentration, etc.). Besides being valuable for many practical applications, microgels are also extremely important to tackle fundamental physics problems. As a result, these last years have seen a rapid development of protocols for the synthesis of microgels, and more and more research has been devoted to the investigation of their bulk properties. However, from a numerical standpoint the picture is more fragmented, as the inherently multi-scale nature of microgels, whose bulk behaviour crucially depends on the microscopic details, cannot be handled at a single level of coarse-graining. Here we present an overview of the methods and models that have been proposed to describe non-ionic microgels at different length-scales, from the atomistic to the single-particle level. We especially focus on monomer-resolved models, as these have the right level of details to capture the most important properties of microgels, responsiveness and softness. We suggest that these microscopic descriptions, if realistic enough, can be employed as starting points to develop the more coarse-grained representations required to investigate the behaviour of bulk suspensions.

that happen at the atomic and molecular length-and time-scales but whose effects extend well beyond them. For example, it is known that the volume-phase transition (VPT) exhibited by PNIPAM microgels is connected to the good-to-bad solvent transition that single PNIPAM chains undergo in water at the lower critical solution temperature (LCST). Below the LCST the polymeric chains are in an extended conformation, while above this temperature they collapse into a globule state due to the complex interplay between PNIPAM-PNIPAM and PNIPAM-water interactions, which are temperature dependent, as schematically shown in Figure 1(a). This phenomenon can be (and has been) carefully investigated with all-atom simulations. In particular, the high resolution of the atomistic models have helped evaluating the molecular parameters that can selectively tune the LCST value and therefore the corresponding VPT temperature (VPTT), with direct impact on the technological applications of microgels. Indeed, atomistic simulations showed that the LCST increases by decreasing the degree of polymerization, thereby allowing to correlate the dependence of the transition temperature on the degree of polymerization to the effects of the chain length on the accessible conformations [13]. Stereochemistry is another important molecular factor that influences the LCST value and that can be better understood with all-atom models. Simulation studies carried out on PNIPAM oligomers with different isotactic content revealed that the LCST value of PNIPAM chains with a high meso diad content, the isotactic PNIPAM, is lower than that of atactic oligomers because the diad composition affects the size, conformation and water affinity of the polymeric chain [14]. Below the LCST the isotactic stereoisomer prefers conformations with a lower radius of gyration and shows higher hydrophobicity. The effect on the LCST of another form of isomerism, such as the structural isomerism, can also be tackled with an atomistic approach [15]. Interestingly, these detailed simulations demonstrate that polymeric chains with the same atomic compositions but different interconnections do not show the same LCST value, since hydrophobic interactions are affected by the spatial arrangement of the functional groups.
Atomistic simulations are also invaluable tools to gain insights into the molecular mechanism that drives the coilto-globule transition [16,17]. The chain transformation from an extended to a globule state was shown to occur with a relatively small loss of hydration water molecules and through a cooperative process that can be ascribable to the breaking of the hydrogen bonding network formed by water molecules in the proximity of hydrophobic groups. PNIPAM chains appear largely hydrated even above the LCST, and the coil-to-globule transition takes place with a significant rearrangement of the hydration water structure.
Atomistic modelling can be extended beyond single linear chains, as more complex molecular architectures can be represented and studied with all-atom simulations. This is the case, for example, of systems of crosslinked PNIPAM oligomers [18], polymeric membranes [19] or even portions of microgels [20]. As shown in Figure 1(b), it is now possible to simulate a cubic section (of linear size 5 nm) of a microgel particle. These simulations have allowed to gather information on the temperature-dependence of its dynamical properties, finding a quantitative agreement with experimental results [20].
Another example of the application of atomistic modelling is provided by the phenomenon of co-non-solvency, which occurs when organic solvents such as short chain alcohols are added to aqueous solutions of PNIPAM. Even though water and alcohol are individually good solvents for PNIPAM, a mixture of the two cosolvents induces a collapse of the FIG. 2: Snapshots of microgels generated with different protocols. (a) a diamond-lattice microgel, (b) a microgel assembled by pre-formed chains [41], (c)-(d) microgels built from mixtures of patchy colloids assembled under two different concentration conditions [42].
polymer for intermediate mixing ratios [21,22]. In the context of microgels, this counterintuitive phenomenon can be used to realize microgels that swell upon increasing the temperature above the VPTT [23]. For example, in the case of methanol-water mixtures it was shown that adding a small amount of methanol to a water solution of PNIPAM promotes a deswelling of the microgel, as PNIPAM and methanol molecules experience a favourable interaction that drives the collapse. By contrast, when methanol is added in excess the microgel re-swells in force of a favourable entropic contribution [24]. Recently, detailed models have been used to quantitatively understand the origin of the co-non-solvency experienced by polymers dispersed in water/alcohol mixtures in both good and bad solvents [25], providing an additional tool to control the responsiveness of microgels.

III. MONOMER-RESOLVED MODELS
In order to explore longer time-and larger length-scales, it is common to simulate complicated systems by using coarse-grained representations, which map groups of atoms or molecules onto single interaction sites [26]. In the context of polymeric systems, the most prominent properties of polymers are known to be scale-invariant (at least in the limit of high polymerisation degree) [27], and hence rather insensitive to the microscopic (atomistic) details. This feature has been exploited to devise techniques that allow to systematically coarse-grain polymeric systems [28][29][30], and chains are often modelled as collection of beads of size σ m connected by springs. In this representation, the size of the single beads is taken to be comparable with the Kuhn length of the real chain [31] (∼ 1 nm [27,32]). For polymers in good-solvent conditions, the actual force-field used in MD simulations is usually the Kremer-Grest set of interactions [33] which is the de-facto gold standard. In this model the connectivity between bonded neighbours is provided by a finite-extension non linear elastic (FENE) spring, whereas the steric repulsion between all beads (bonded and non-bonded) is modelled through a Weeks-Chandler Andersen (WCA) interaction [34]; the parametrisation of these potentials makes sure that, under ordinary conditions of temperature and concentration, chains do not cross each other, so that the overall topology of the system is preserved. The Kremer-Grest model can be augmented by additional terms that can describe more specific cases, such as charged or semiflexible polymers [35][36][37]. The great majority of the numerical work done on polymer-based objects deals with systems made of chains or aggregates with simple topologies such as rings [38], star polymers [39], or dendrimers [40] since, from the numerical point of view, the disordered nature of polymer networks poses significant challenges, in terms of both modelling complexity and computer time. Several strategies for the generation of suitable network topologies have been put forward, which are discussed and compared in depth in the next few paragraphs. A visual overview of some of the models discussed is presented in Fig 2. The description of the different preparation protocols is completed by results on the swelling transition and to the swelling/deswelling kinetics of in silico microgels.
A. Protocols for the numerical design of coarse-grained microgel particles

Microgel formation from a crystalline lattice
Early attempts of modelling particles made by a crosslinked polymer network rely on placing the crosslinkers on a crystalline lattice (usually diamond) and connecting them with chains having the same size. The network is then cut-out from a sphere to obtain the shape of the particle. Following the same approach it is possible to generate both standard and hollow microgels, which have been simulated to study the uptake and release of neutral species [5] or their behaviour at the liquid-liquid interface [3]. Recently, the diamond lattice procedure has also been employed to generate microgels made of interpenetrated polymer networks [43]. In addition, the possibility of having chains of the same length represents an advantage when comparing numerical simulations with theoretical approaches. For instance, one can exploit scaling arguments based on the Flory theory for polymers or can apply the blob tension model for evaluating the polymer stretching for different crosslinker concentrations as for instance done for ionic microgels [44,45].
While this approach is suitable to showcase the potentialities of microgels for applications and for a direct test of the theory, it is far from being a realistic model. Indeed, lattice-based topologies suffer from several drawbacks: (i) all chains have the same length, which is quite unlikely to happen in real microgels; (ii) the distribution of crosslinkers is uniform by construction, while in the polymerization process employed to synthesize standard neutral microgels it is known that crosslinkers react faster than monomers, giving rise to an inhomogeneous distribution, with more crosslinkers being in the core than in the outer corona [1, 7]; (iii) the corona is obtained from the spherical cut of the crystalline lattice and its extension can be controlled only by changing the chain length; (iv) the model has no loop-like defects and few or no dangling ends [46], which are fundamental for a comparison with real microgels, as they contribute to the hydrodynamic radius R h and might also play a role in the rheological properties of microgel suspensions at moderate concentrations [47]; (v) polymer chains do not entangle, which means that the elastic properties and the permeability of the diamond microgel depend on the length of the chains but not on the degree of entanglement of the polymers. Finally, the underlying crystalline structure of the diamond-lattice-based microgel affects the numerical density profiles and the form factors [42], making a comparison with the experimental data difficult.

Microgels from randomly distributed crosslinkers
A step forward with respect to the previous protocol can be taken by randomly distributing crosslinkers within a cubic simulation box; close-by crosslinkers are then connected by polymer chains, for example by choosing a given cutoff distance [48]. This allows to generate non-ordered networks made of polymers chains that are slightly polydisperse. As for the crystalline-lattice-based microgels, the spherical shape is obtained from cutting out a sphere from the cubic simulation box. The main advantage of this method is that the crosslinker distribution can be controlled easily, even though the connectivity among crosslinkers is not completely satisfied and thus cannot be fully controlled. We further note that this protocol makes it possible to easily generate core-shell [49] or even hollow microgels [50], since the idea behind the assembly of the particle is similar to that employed for crystalline-lattice-based methods.

Microgels from the self-assembly of a gel network
Taking inspiration from the synthesis process of PNIPAM microgels, a recent numerical protocol has been developed to design coarse-grained microgel particles. The approach is based on the self-assembly of patchy particles, i.e. hard-sphere particles decorated with attractive sites, which have shown to form gel networks at low and moderate densities [51]. In Ref. [42] bivalent and tetravalent patchy particles are used to mimic, respectively, monomers and crosslinkers. Inter-and intra-species bonds are allowed except for crosslinkers particles that cannot form bonds among themselves. Although the interactions among monomers and crosslinkers resemble those occurring in the synthesis of real microgels, the dynamic processes that bring to the network formation in the two cases have little in common. Indeed, while the polymerization mechanism in real microgels is an off-equilibrium process, the network formed by patchy particles occurs in equilibrium, and is facilitated by a "swapping" mechanism [52] which allows the system to easily equilibrate even at low temperatures, where the fully-bonded-network condition can be almost accessed. The self-assembly process allows to generate a disordered network where the length of the polymer chains has an exponential distribution which can be predicted following a heuristic argument based on the Flory theory in the fully bonded limit [53]. Instead of cutting out a spherical region from a bulk homogeneous network, a spherical confinement (mimicking confinement in a droplet) is employed during the self-assembly process. Such external field acts as an extra parameter which can be tuned to influence the topology of the network: by varying the radius of the spherical confinement at fixed crosslinkers concentration it is possible to generate microgel particles with different density, topology and degree of entanglement. The resulting microgels can range from compact to floppy with several dangling ends, experiencing very different swelling behaviors [42]. As a result, this protocol makes it possible to investigate the role of topology on the dynamics of swelling in a unique way. As for the case of crystalline-latticebased microgels, the resulting network possesses an homogeneous distribution of crosslinkers, being the network an equilibrated binary mixture. However, the corona spontaneously arises from the interfacial region formed by the system due to the presence of a confining field; with this approach the width of the corona can be controlled by the thermodynamic properties of the network, i.e. by temperature and density.

Microgels from the assembly of functionalized chains
Differently from the standard synthesis process based on precipitation polymerization, another route to synthesize microgels is through microfluidics fabrication using droplets of macromolecular precursor chains that are later photocrosslinked [54]. Inspired by this technique, a new numerical protocol for assembling microgel particles has been recently devised, exploiting the self-assembly of pre-formed chains which are functionalised with reactive groups [41] and placed in a spherical confinement. In the method presented in Ref. [41], a fraction f of reactive groups are placed randomly on the polymer chain, with the constraint that consecutive reactive sites are not allowed in the backbone sequence. During the dynamics reactive sites form permanent bonds; this gives rise to a fast stage in which the majority of reactive groups bond together, followed by a second slow stage in which non-bonded reactive sites seek for other reactive groups until full crosslinking is achieved. The latter stage is sped up by selecting randomly two non reacted sites and by applying an attractive external field between the two that allows them to get in contact and form a bond. Differently from crystalline-lattice-based microgels, this procedure allows to design microgels with entangled polymer chains of different sizes, together with a conformational polydispersity which is fundamental for investigating the role of topology in the swelling dynamics of the particle. In addition, the number of the precursor chains is independent on the number of crosslinkers, which allows to prepare microgels with different densities at fixed crosslinker concentrations. Finally the presence of a spherical confinement provides a spherical shape to the assembled-network without the need to cut it out from the bulk of the polymer network. A similar procedure was employed in Ref. [55] where the reactive sites are not chosen at the beginning of the simulation, but at the end. Namely, a number of polymer chains are equilibrated within a confined network and then crosslinkered by selecting randomly monomers among those separated by a maximum distance. Unlike in the procedure of Ref. [41], only the crosslinking of sites belonging to different chains is allowed, which artificially suppresses the formation of intrachain loops. If the final crosslinker concentration is smaller than the one desired, the maximum bonding distance is increased and the procedure is repeated. Although similar to the technique described previously in this section, this strategy gives rise to non-compact microgels even when a high crosslinkers concentration is chosen. The form factor of microgels generated with different numerical protocols across the volume-phase transition. The two microgels assembled with the protocol of Ref. [42] have been generated within two different spherical confinements, here indicated by the overall density of the mixture ρi.

Comparison of the assembly protocols: topology and form factors in the swollen state
We conclude the overview of the protocols used to build different topologies by comparing some of the microscopic architectures generated. Figure 2 shows representative snapshots of the diamond microgel and of microgels generated with the methods of Refs. [41] and [42]. Figure 3 shows the form factor P (q) of microgels composed of N ≈ 21000 monomers and a crosslinker concentration of ≈ 1.2%. The form factor of the diamond-lattice microgel, which we use here as a reference, displays a peak at a position that roughly corresponds, in real space, to the size of the particle. At larger wavevectors (0.2 ≤ qσ m ≤ 0.4), the form factor exhibits a weak dependence on q that is due to the underlining ordered structure of the network and can therefore be considered spurious [42]. The other curves refer to disordered topologies, either generated with the method by Moreno and Lo Verso [41], or through the assembly of patchy particles binary mixtures [42]. In all cases we observe a peak or a shoulder around qσ m ≈ 0.15 that is linked to the size of the microgels and, for qσ m ≥ 0.3, very similar decays, reflecting the self-avoiding character of the strands [27]. Contrarily to the diamond case, the form factors of each disordered topology are compatible with the fuzzy-sphere model [1]. The difference between the different topologies is concentrated in the intermediate q-region: microgels generated at higher densities display a more structured P (q). The lack of well-resolved peaks in the form factor of the microgel assembled with functionalized chains [41] signals the larger heterogeneity of the network compared to the case of microgels generated by assembling patchy mixtures.

B. Swelling and solvent effects
After assembling the network, the next relevant issue to address is reproducing the swelling/deswelling transition of microgels. To this aim, the solvent effects must be taken into account, and this can be done at various levels of coarse-graining. Apart from the atomistic route discussed extensively in Section II, even for a single microgel the simulations in the presence of a coarse-grained solvent can be computationally expensive. Thus, it is preferable to use implicit solvent models to reproduce those features of the volume phase transition that are independent of the actual presence of the solvent, such as the thermodynamic and geometrical properties across the VPT, and then to resort to explicit solvent models to tackle specific problems for which the presence of the solvent is absolutely needed, such as for example solvent expulsion, kinetic aspects and interfacial problems.
We now start to discuss the so-called implicit solvent models, where an effective 'solvophobic' potential between the monomers which takes into account the affinity between polymer and solvent is introduced. The solvophobicity is normally modulated by a control parameter which plays the role of the temperature (or of the external parameter controlling the swelling). From the practical standpoint, the effective potential is a monomer-monomer interaction that is zero under good solvent (or maximally swollen) conditions, and becomes very attractive under bad solvent (or collapsed) conditions. To this purpose, the use of a simple Lennard-Jones potential may give rise to unphysical non-monotonic behavior of the microgel size [56] with increasing quench depth due to the relative contribution of the attraction and the repulsion when the potential depth increases. Instead, the potential initially proposed by Soddemann and coworkers [57] was found to well reproduce the swelling behavior for microgels assembled in different ways [41,42]. In this model, the polymer-solvent interaction is modulated by the solvophobic parameter α which plays the role of the temperature in experiments. We further note that the specific form of this solvophobic interaction is not expected to play a major role in the swelling behavior, at least from the qualitative point of view, and thus different choices could be adopted, as for example the potential used to model star polymers [58] or telechelic star polymers [59] in solvents of different quality. The relative change of gyration radius of microgels generated with different numerical protocols across the volumephase transition. The two microgels assembled with the protocol of Ref. [42] have been generated within two different spherical confinements, here indicated by the overall density of the mixture ρi. Here α is the parameter of the interaction that controls the quality of the solvent, playing the role of the temperature in real PNIPAM microgels.
Interestingly, Fig. 4 shows that the swelling curve is not very sensitive to the inner topology of the network. Indeed, microgels generated in different ways at approximately the same crosslinker concentration, including also the diamond lattice microgel, display very similar swelling properties and even the same VPTT [41]. Of course, with an underlying ordered lattice, the only way to tune the topology is to vary the chain length, affecting the crosslinker concentration. Instead, using disordered assemblies, the variation of the confining volume can be used to significantly alter the swelling properties [42,53]. The incorporation of an explicit solvent is the next step of description. Of course, this cannot be done at the atomistic level of accuracy, but still sometimes simple potentials like Lennard-Jones or modifications there-of have been adopted. These have the main disadvantages that excluded volume of the solvent molecules can be sometimes overestimated [60]. Thus, it is much better to rely on a coarse-grained solvent representation, where groups of solvent molecules are treated as soft beads. This is precisely the aim of the Dissipative Particle Dynamics (DPD) technique, which has the advantage of correctly reproducing hydrodynamics at long times by imposing locally the conservation of momentum [61]. In addition, the DPD method has been mapped to polymer-solvent interactions and provides a way to directly relate the parameters of the involved soft potentials to the Flory-Huggins solvency parameter [62]. However, in order to do so, it is necessary to use such soft potentials among all species involved, including monomer-monomer interactions. This may give rise to unphysical crossing between polymer chains and care must be taken when adopting this method [43].
DPD has been used in a number of studies [63][64][65] performed with a regular diamond network. It was also used by Nikolov and coworkers for a topology obtained by using randomly distributing crosslinkers [48]. However, in order to compare the effect of the explicit solvent with the widely used implicit ones, a one-to-one correspondence must be established. This was the aim of a recent work [60] where identical microgel configurations, interacting through the usual Kremer-Grest force field, were compared in implicit and explicit solvent conditions for both swelling curves and form factors. It was shown that a DPD treatment of the solvent gives a faithful representation of the implicit solvophobic potential in all aspects, putting the basis for a systematic use of the explicit solvent to investigate interfacial aspects of microgels. For example, an interesting aspect to model is the flattening of the soft colloids, and particularly microgels with their inhomogeneous core-corona structure, at a liquid-liquid interface [66], that is relevant for applications as emulsion stabilizers [67].
Finally, some studies have also adopted the more accurate Multi-Particle-Collision-Dynamics to treat the solvent [45,56], but only for the diamond-lattice topology. These simulations have shown that the monomer dynamics under swollen conditions agree with the predictions of the Zimm model [27,68]. However, as the microgel shrinks by decreasing the solvent quality, the dynamics progressively deviates from this theoretical model and, in the fully collapsed state, hydrodynamic interactions are screened out while the dynamics approaches the predictions of the Rouse model expected for polymer melts [27,68].

C. Kinetics of swelling and deswelling
Monomer-resolved simulations make it possible to study in detail the evolution of the microgel internal structure under changing the solvent conditions -a feature that is not easily accessible in experiments -and to unravel the effect of the network microstructure on the kinetics of swelling and deswelling. The time evolution of the microgel radius of gyration during its collapse, R g (t), was analyzed in Refs. [41,49,56,60] for diamond and disordered microgels. Ref. [48] also analyzed the swelling, finding consistency with Tanaka's theory [48]. Some general trends were observed for all models of microgels. In particular, higher degrees of crosslinking [48], higher regularities of the network [41] and deeper quenches (to poorer solvent conditions) [41] all result in faster and less stretched decays of R g (t). Very interestingly, the shape of R g (t) is apparently independent of the solvent model [60].
As the microgel collapses when it is driven beyond the VPT point, the monomers start to form local globules that progressively merge into interconnected larger domains, until the whole structure is finally joined into a single dense spherical globule. This phenomenon is known as 'coarsening' and is universally observed in phase separating systems [69][70][71]. Direct visual inspection of simulation snapshots at intermediate times of the coarsening reveals rather different conformations depending on the topology of the network. Both regular and disordered networks with relatively high degree of crosslinking are approximately spherical at all times, from the initial swollen to the final collapsed state [41,48]. Instead, disordered ones with low degree of crosslinking display at intermediate times irregular conformations, with significant asphericity [41,49] and large globulated protrusions [41,60]. Two illustrative examples of these conformations are presented in Fig. 5(a)-(b).
The homogeneous or heterogeneous character of the collapse of the outer shell is also reflected in the monomer density profiles calculated with respect to the center of mass. In regular networks or in dense disordered ones, the initial stage of the coarsening is characterized by a strong monomer aggregation in the outer shell, while the core of the microgel remains 'hollow' [56]. This effect is much less pronounced in the heterogeneous collapse of disordered low-density microgels (see main panel of Fig. 5(c)), for which a flat density profile in the core is quickly reached. The coarsening kinetics of the microgel deswelling has been characterized in Ref. [41] by measuring the length of the growing domains. A clear difference was found between the regular diamond networks and the disordered microgels constructed by crosslinking of functionalized chains. The latter show a power-law time dependence for the growing domain length, which has been suggested [41] to be an intermediate scenario between liquid-gas phase separation [72] and collapse of linear chains [73]. Though a similar power-law domain growth is observed at early times in the diamond networks, an accelerated growth is found at the late stage of the coarsening [41]. The observed master functions for the domain growth in both kinds of microgels are independent of the depth of the quench (i.e., of the solvent quality parameter) [41]. A similar result has been found in Ref. [49]. Remarkably, in analogy with general observations for critical phenomena, a scaling relation between dynamic correlators and the growing domain size has been found [41], with the scaling function being independent of the network microstructure. Finally, it is worth mentioning that, although no quantitative analysis has been reported, the conformations presented as simulation snapshots in Ref. [48] display no significant coarsening in swelling microgels [48], suggesting a much more homogeneous character of the expansion of the monomers starting from the collapsed globular state.

IV. FURTHER COARSE-GRAINING
Microgel suspensions have become a model system in fundamental physics, allowing to shed light on diverse phenomena such as jamming [74,75] and glass [76] transitions, charge effects [77,78], and more [1]. In this context, the complex internal architecture and the resulting responsiveness of the single microgels are crucial ingredients that can be harnessed to steer the collective behaviour of the system. However, a numerical description of a bulk system that make use of monomer-resolved models, which include these features by construction, is out of reach. Indeed, such a detailed description would require the total number of degrees of freedom to be so large that accessing the lengthand time-scales that are characteristic of the phenomena of interest would be impossible with modern-day numerical resources. The obvious solution is to employ much simpler models that contain the minimal number of ingredients required to observe the desired bulk behaviour.
A simple approximation is to describe microgels as spheres that interact through a soft potential, which usually FIG. 6: The effective interaction between two microgels with 5% crosslinker concentration, as computed with monomer-resolved numerical calculations (points) and as estimated by using the simple Hertizan form of Eq. (1) (line) [82]. The snapshots on the left and on the right show two representative conformations at large and small separations r, as indicated by the circles.
takes the form where U 0 is a prefactor linked to the overall softness of the interaction, σ is an effective particle diameter and n is commonly set to 2 (harmonic potential) or 5/2 (Hertzian potential). The latter choice can be theoretically justified by leveraging the classical elasticity theory (CET) [79]. Indeed, in the CET framework two elastic spheres in contact experience an Hertzian effective repulsion. The CET assumes that the two objects are homogeneous and is strictly valid only in the small deformation regime, i.e. when the separation between the two particles is not too small. In early works, microgels were compared to hard sphere behavior, particularly for the dependence of the zero-shear viscosity on packing fraction φ and found that up to about φ ≈ 0.5, no significant differences were observed [80]. However, above this threshold, while the hard sphere viscosity would diverge close to φ ≈ 0.58, data measured for microgels show a clear deviation. The rheological data of elastic moduli were found to obey a power law increase that would be compatible with a soft sphere potential, particularly of an inverse power law form with exponent between 9 and 12. However, later on, experimental evidence based on microscopy measurements in dilute conditions lended support to an effective Hertzian potential [74]. These findings were confirmed by quantitative comparisons based on confocal microscopy experiments in the fluid phase and simulations which showed that, treating the microgels as elastic Hertzian spheres, a good description of the radial distribution functions across the whole fluid concentration region is obtained [81]. In order to provide a numerical test of the CET assumptions and of the range of validity of the Hertzian model, the calculation of the effective potential between two microgels is required. Such a study was done for two diamondlattice-based small microgels by Ahualli and coworkers [83], who reported that the resulting effective potential was not compatible with the Hertzian model, but rather with a generalized form of it, where the exponent n in Eq. 1 varies, being best described by the value 3.5. However, this description is purely phenomenological and not based on any elasticity theory concept. Going one step forward and really verifying the Hertzian model would require the calculation of the elastic moduli of individual microgels, entering in the prefactor of Eq. 1 for n = 5/2, namely U 0 = (2Y σ 3 )/(15(1− ν 2 )) where Y is the Young modulus and ν the Poisson's ratio. Using microgels with a disordered topology, generated from the self-assembly of a gel network with the patchy mixture method described in Section III A 3, recent simulations provided the numerical evaluation of all elastic moduli as well as of the microgel-microgel effective interaction. This work, differently from [83], confirmed that the effective potential is well described by the Hertzian model at small deformations or, equivalently, for repulsion energies of few K B T 's [82]. Above these limits, strong deviations are found, as shown in Fig. 6, where it is evident that the Hertzian picture breaks down when particles come into close contact, becoming anisotropic. These results have also been quantified in a packing fraction range of validity of the Hertzian model, by means of bulk simulations with the Hertzian model and with the numerically calculated potential, yielding good agreement between the two up to packing fractions φ ∼ 0.8−1.0, depending on the crosslinker concentration [82]. These findings confirm and extend earlier results in the fluid regime, but cast relevant doubts about the use of the Hertzian interaction to describe dense microgel suspensions that undergo glass or jamming transitions, a procedure that is largely employed in the literature [74,84,85].
At high packing fractions, not only the Hertzian picture breaks down, but also (and probably most severely) the two-body approximation. Thus many-body effects must be taken into account. This is currently a rapidly evolving field of research activity. Among the early attempts to go beyond the Hertzian model, with a simple but effective modification is the numerical work of Urich and Denton [86], introducing the possibility of changing the particle size isotropically. This approach has also been extended to charged microgels [87], showing that the model is able to capture the swelling behavior as a function of concentration observed in experiments [88]. A step forward would be to also allow shape changes, as done in an early work related to crystallization aspects only [89]. Several recent works have also accounted for osmotic compression and deswelling to compare with experimental systems within simple models. In particular, van der Scheer and coworkers provide a phenomenological description of the internal equation-of-state of soft particles to predict the behavior of the collective relaxation [90], while de Aguiar et al employ a modification of the Flory-Rehner theory which explicitly includes the stretching of the chains [91]. All these efforts represent concrete steps forward with respect to simple pair-wise models, but additional work will be needed in the near future to develop a more microscopic model which takes into account the internal degrees of freedom of realistic soft particles, thus being able to describe interpenetration [92], deformation and faceting [75,91]. To consider these aspects in more refined micromechanical models [93] will be an important issue for a more quantitative prediction of rheological properties of dense and jammed states.

V. PERSPECTIVES
There are several directions towards which future work aiming at providing a numerical description of microgels should proceed. Regarding atomistic simulations, there are a few aspects that could be explored in order to better understand the structural and dynamical properties of polymer-solvent interactions. In particular, more complex topologies, such as polymeric networks of different degree of disorder, could be exploited to investigate the molecular behaviour across the VPT and to evaluate the dominant interactions throughout the process. Another possibility, that takes advantage of the explicit description of the solvent at the atomistic level, is to use different architectures to better understand the interplay between water and polymer and to characterize to what extent the solvent affects the microgel properties. In this context, a very promising use of these atomistic approaches would be to quantify the interactions taking place polymeric oligomers at the molecular level and incorporate them in more coarse-grained models.
Concerning the monomer-resolved approaches that have been proposed so far to produce disordered networks, we also foresee several strategies for improving the current understanding. Starting from the protocol based on patchy particles self-assembly, it is important to stress that, while this was not designed in order to reproduce realistically any experimental synthesis process, it aims to achieve a fine control of the resulting assembly product. To achieve this goal, it appears utterly important to be able to control the distribution of crosslinkers inside the network, in order to produce more and more realistic network topologies which would closely resemble experimental measurements, e.g. by recent super-resolution microscopy [94]. In parallel, such a control would open up the possibility to tailor the resulting density profiles at will. Since the assembly of the patchy mixture occurs in equilibrium, such a finer degree of control might be obtained by either changing the force field or by adding external fields. Another possible extension of the protocol presented in Ref. [42] would be to change the numerical procedure so as to explicitly reflect the experimental synthesis at the microscopic level. For example, since polymerization is a one-directional process, in a realistic numerical synthesis the chains would start growing from a small fraction of one-patch particles, which would play the role of the initiators of the polymerization process. The free monomers, or the end monomers of a growing chain, would then only react with the end monomers of other growing chains. Likewise, branching would only occur at the ends of the growing chains. Another ingredient that might be incorporated in the model is the difference between the reaction rates of chain growth and branching. Thus, the fraction of initiators and the rates for chain growth and branching would act as additional control parameters, which can be adjusted to control the topology of the resulting network.
About the numerical protocol based on crosslinking of pre-existing chains [41], it already reflects the experimental synthesis from a microfluidic approach in a reasonably realistic fashion. However, possible improvements of this method should take into account the directionality of the interactions between the reactive groups, thereby penalizing the formation of intra-chain small loops. In addition, a question to address in the future is how the mechanisms implemented to accelerate the synthesis, such as random crosslinking of the last unreacted groups in [41] but also bond swapping in [42], may affect the resulting microstructure. For instance, random crosslinking of the last groups makes the non-trivial assumption that the energy barriers that impede the formation of the few remaining unbound pairs can be overcome by waiting long enough. On the other hand, bond swapping may prevent the freezing of entropically unfavourable local structures that would emerge in a purely irreversible process, making the system less heterogeneous. Understanding the effect of these mechanisms on the final network topology would further enrich the range of possibilities to generate different microscopic architectures.
The ultimate aim of an effective multi-scale modeling approach is to transfer the knowledge obtained at a smaller scale to the next level of description. In this respect, a very promising step forward would be to use the results obtained through accurate monomer-resolved models as a guidance to develop more coarse-grained models that include some internal degrees of freedom and thus naturally include many-body interactions. A few examples of this kind have been recently proposed and could be potentially promising for microgels. Among these, we recall the Voronoi model widely used to describe biological tissues, which encodes the particle elasticity in the Hamiltonian [95] and the liquid drop model, which has been used to calculate the phase diagram of polymeric particles at large compressions [96]. Finally, explicit models that are able to capture particle deformation and shrinking will be crucial to tackle by simulations the problem of the glass transition at high enough particle volume fractions. A recent study put forward a first model in this direction [97].
As a final note, we stress that, although the present review has focused on non-ionic microgels, many of the issues that we have discussed will be also relevant to accurately model ionic microgels. However, for highly charged systems, the electrostatic force usually dominates over the other contributions, and thus the microscopic details become somewhat less important for the swelling transition [77]. In this context, important questions linked to the counter-ion distribution [98] and to the chemical equilibria of ions in weak polyelectrolyte nanogels [99] have been the subject of recent work. Interestingly, since PNIPAM microgels are also weakly charged [100][101][102], with the effect of such a charge showing up close to the VPTT [103], some of these results may be broadly relevant. Understanding the interplay between the electrostatic interactions and the onset of the swelling transition will be crucial to develop coarse-grained models that can be used across (and beyond) the VPT.

VI. CONCLUSIONS
Here we have presented an overview of the numerical methods and models that have been used to investigate the behaviour of non-ionic microgels at many different length-scales, from the molecular level up to much more coarse-grained descriptions. We have put particular emphasis on the modelling of PNIPAM (thermoresponsive) microgels, which are increasingly being used as model systems to investigate fundamental problems in condensedmatter physics, but many of the results and methods reported here can be extended to other types of microgels. We have highlighted the inherent multi-scale nature of microgels. Indeed, part of the behaviour of individual microgels, such as the temperature at which thermoresponsive microgels deswell, can be directly traced back to the properties of the polymeric repeating units they are made of. However, other fundamental quantities, such as the swelling ratio or the single-particle elastic moduli, depend on the mesoscopic architecture of the network. The bulk (macroscopic) behaviour, in turn, is controlled by all these properties as well as by external parameters such as temperature, pH and salt concentration. Modelling microgels is thus a multi-faceted challenge that cannot be addressed with a single tool or technique.
We have shown that in the last years there has been a flourishing of numerical studies on microgels. However, most of the effort has been devoted towards building realistic microgels to understand the single-particle rather than the bulk properties. In the meantime, the fast development of synthesis, imaging and scattering techniques have made it possible to experimentally probe dense suspensions of microgels to shed light on important open issues such as the jamming and glass transitions. Time is ripe now for the numerical community to catch up and use the knowledge gained by investigating the single-particle properties as a guidance for developing models which are simple enough to allow for bulk simulations but still take into account some of the details of the inner structure of microgels, as prefigured in the Perspectives section. A first challenging task for such a model would be to provide a microscopic explanation of the fragile-to-strong transition observed in suspensions of microgels of increasing softness [76,104].