Microgels at interfaces behave as 2D elastic particles featuring reentrant melting

Soft colloids are increasingly used as model systems to address fundamental issues such as crystallisation and the glass and jamming transitions. Among the available classes of soft colloids, microgels are emerging as the gold standard. Since their great internal complexity makes their theoretical characterization very hard, microgels are commonly modelled, at least in the small-deformation regime, within the simple framework of linear elasticity theory. Here we show that there exist conditions where its range of validity can be greatly extended, providing strong numerical evidence that microgels adsorbed at an interface follow the two-dimensional Hertzian theory, and hence behave like 2D elastic particles, up to very large deformations. We are also able to estimate the Young's modulus of the individual particles and, appropriately comparing with its counterpart in bulk conditions, we demonstrate a significant stiffening of the polymer network at the interface. Finally, by analyzing dynamical properties, we predict multiple reentrant melting phenomena in experimentally accessible conditions: by a continuous increase of particle density, microgels first arrest and then re-fluidify due to the high penetrability of their extended coronas. The present work establishes microgels at interfaces as a new model system for fundamental investigations, paving the way for the experimental synthesis and research on unique high-density liquid-like states. In addition, these results can guide the development of novel assembly and patterning strategies on surfaces and the design of novel materials with desired interfacial behavior.


I. INTRODUCTION
Mesoscopic assemblies of colloids and nanoparticles display features that depend critically on the microscopic details of the building blocks, e.g. composition, size and shape, as well as on the specific macroscopic physical conditions such as the thermodynamic control parameters. It is by carefully choosing and tuning these variables that one can induce the formation of different structures and explore various states, such as liquid-like fluids, glasses or crystals [1]. At the core of this collective behaviour is the interparticle interaction, which ultimately dictates the phase behaviour and dynamics of the assembly: information at the multi-particle level is thus fundamental to determine the properties of the material. If a system is made of rigid building blocks that interact only through excluded volume interactions, it can be approximately mapped to a hard-sphere system and its behaviour can be investigated through packing models. These have been used for a long time to successfully answer fundamental questions in physics and material science whenever simple constituent units are involved [2].
However, in certain cases, the complexity that resides at the microscopic level cannot be described in these terms. This is especially true for soft polymeric colloids that possess internal degrees of freedom endowing them with elasticity and deformability. Among the available library of soft deformable particles, microgels, colloidal-sized crosslinked polymer network, are one of the finest illustrations of this concept [3][4][5]. Their structure is determined by the chemical synthesis conditions that, in the common procedure of precipitation polymerization [6], lead to the formation of spherical particles made of a compact core and a fluffy external corona [7]. Although microgels are often considered as simple elastic particles that can be modeled with the classical elasticity theory, the presence of multiple length scales in their internal architecture makes their effective interactions in bulk more complex than a pure Hertzian model, and calls for more refined treatments that range from a phenomenological multi-Hertzian model [8] to descriptions that depend on the concentration regime [9,10].
The intrinsic softness of microgels and, in general, of soft deformable objects, is fully revealed, and can be taken advantage of, at interfaces, which can be used to fulfil different purposes. In fact, if the interfacial tension is large enough, it is possible to coat an interface with nano-or microsized particles which then remains adsorbed and can form stable monolayers. This concept can be used, for instance, to stabilize emulsions and biological membranes [11][12][13][14] or to study fundamental self-assembly phenomena in 2D such as crystal [15] or quasi-crystal formation [16,17]. Moreover, by changing in situ the single-particle properties and the local environment, it is possible, in principle, to finely control the stability and the structure of the whole monolayer/emulsion [18,19]. Noteworthy, in this respect, is the possibility to realize complex patterns whose application as etching masks in nanolithography can lead, for instance, to the fabrication of nanowire arrays [20,21].
When adsorbed at interfaces, microgels flatten out and adopt the so-called "fried-egg" shape, making them very different from their bulk counterparts [22][23][24]. In particular, at an interface, the polymer network tends to minimize as much as possible the surface tension between the two liquids by taking a stretched configuration already in the dilute regime. So far, the effects that an interfacial confinement induces on such particles have been mainly limited to the study of their characteristic shape, and hence to their structural arrangement at oil-water and air-water interfaces. These aspects have been widely investigated experimentally [25,26] and, more recently, also numerically [24,27,28].
Instead, little is known on the collective behavior of such particles, mostly by means of indirect experimental feedback [29][30][31][32][33][34], which has neither allowed to extract a functional form for the interaction potential nor to properly understand the role of the surface tension. In the same way, a true characterization of the elasticity of the polymer network within the interfacial plane is still missing, being limited both experimentally and numerically by subpar techniques and models. Clearly, a simple transfer of results from bulk to interface would be highly inappropriate, due to the dramatic change of conditions between the two cases.
This lack of understanding hampers the progress towards further applications, since an established fundamental knowledge of the basic constituents would make it possible to a priori design and guide the assembly of innovative materials and nanostructures. From a theoretical standpoint, it also prevents the adoption of microgels at interfaces as model systems for the study of open questions in fundamental science [35,36]. In this sense, it is crucial to provide a microscopic understanding of such system.
In this work, we address this problem by calculating both the effective interactions between two microgels at a liquid-liquid interface and their individual elastic properties, using this knowledge to predict their multi-particle response at high densities. Our approach relies on state-of-the-art modeling of singleparticle microgels that was shown to quantitatively capture the internal topology of laboratory microgels both in the bulk [37] and at the interface [24]. In the latter case an appropriate framework, that explicitly takes into account the effects of the surface tension between the two liquids, has been developed in order to correctly describe the deformation of the microgels [24]. Despite the complex arrangement of the polymer network at the interface and the intrinsic presence of a core-corona structure, the calculation of the effective interactions for different crosslinker concentrations c reveals a remarkable agreement with the 2D expression of the Hertzian model for elastic disks for all investigated distances. This is clearly different from what was found for the same system in bulk conditions [38] and establishes the validity of the two-dimensional Hertzian model for microgels at interfaces up to large compression regimes. The Young's modulus determined from the effective potential is also directly compared to explicit calculations based on elasticity theory for small and intermediate deformations. Thanks to this method we are able to achieve a full characterization of the elastic response of the microgels in the two-dimensional interfacial plane and we can thus establish a sound comparison to the three-dimensional bulk case. Notably, our results show that the elastic moduli, once converted to their three-dimensional counterparts, are roughly one order of magnitude larger at the interface than for the same microgels in the bulk. This highlights the key role of the interfacial tension in stiffening the microgels due to the stretching of their coronas.
Having determined how such complex particles interact with each other, we are finally able to carry out our study also at the collective level by investigating the dynamical phase behavior of an ensemble of these 2D effective elastic disks. When softness and elasticity are taken into account in the interparticle interaction a rich behavior is, in general, expected [39][40][41]. In particular, we find the presence of multiple reentrant melting phenomena, where a glass is melted simply by an increase in particle concentration. Despite similar findings have long been predicted for simple soft models [42], here for the first time such scenario is found for microscopically-motivated effective interactions and, most importantly, for potential parameters that can be realized in experiments.
The extensive analysis of microgels at the interface presented here, together with the notion that the Hertzian model can be used up to very large deformation energies, sheds light on how single microgel properties and collective response are coupled, and demonstrates that this system is an incredibly promising model for probing the collective behaviour of 2D elastic particles up to large densities. Their unusual dynamical features are of wide interest for the preparation of colloidal monolayers with non-monotonic viscoelastic properties that could be used for a variety of different applications.

II. THE MODEL
Microgel monomers interact through the Kremer-Grest bead-spring model [43,44] in explicit solvent. The polymer network has a disordered topology [45], with density profiles and form factors comparable to experimental ones thanks to an inhomogeneous distribution of crosslinkers within the network [37]. The solvent is treated within the Dissipative Particle Dynamics framework [46] and, for simulations at the interface, their mutual interactions are tuned to mimic a water/hexane interface. Under these conditions the microgel spontaneously adopts the characteristic fried-egg shape when placed close to the interfacial plane. Its structural characterization is in good agreement with experiments both in terms of flattening on the interfacial plane and of protrusion on the preferred water side [24]. In the present work, for the study of Figure 1. Microgels interacting at the interface. Top and side simulation snapshots of two microgels with c = 5% at the water-oil interface at a representative distance r ≈ 40σm. The effective diameter of the microgel is σ eff . Solvent particles are not shown for clarity. the microgel-microgel effective interactions, we simulate particles with N ≈ 5000 monomers of diameter σ m , that defines the unit of length, and crosslinker molar fraction c = 3%, 5% and 10%. Due to the exceptional computational cost to carry out the simulations, we limit our study to a single microgel topology for each studied value of c, checking the consistency of the results with a second topology for c = 5% (see Supplemental Material). The elastic properties of smaller microgels with 2000 and 3000 monomers are then analyzed for further considerations on their collective behavior. Additional details on the interface modeling and simulations are provided in the Methods section and in Refs. [24,37,46]. The typical conformation taken by two interacting microgels at the interface is reported in the simulation snapshots of Fig. 1.

III.1. Effective interaction potential
The two-body effective potential V eff (r) between the microgels at a water/hexane interface is evaluated by means of extensive simulations exploiting the Umbrella Sampling technique [47][48][49], as also explained in the Methods section, and it is shown, rescaled by β = 1/k B T , in Fig. 2  where r is the distance between the centers of mass of the microgels at the interface, σ eff quantifies the extension of the microgel on the interfacial plane and Y is the Young's modulus of the individual particle. The agreement between the numerical results and the theoretical fits is remarkable for all probed distances and all studied values of c. Therefore, it clearly emerges from these findings that two microgels confined at an interface effectively behave as 2D elastic objects, further confirming the soft repulsive nature of their mutual interactions. Further discussion on the functional form of the potential can be found in the Supplemental Material. Experimentally, small microgels -those whose diameter does not exceed few hundreds nanometersare the best candidates to interact in this way, since they do not experience long-range attractions due to capillary effects [32]. Indeed, the latter have been widely reported [52][53][54] and found to be relevant only for microgels large enough to induce a local deformation of the water/oil interface [32]. By contrast, our solvent modeling is aimed essentially at reproducing the surface tension and the microgels solubility, both of which have a direct influence on the conformation of the particle. We can thus directly probe the elastic interactions between the microgels without the interference of attractive capillary forces.
These outcomes also evidence the presence of a single characteristic length in the potential up to a center-to-center distance as small as the interaction radius of the microgel (∼ σ eff /2) for the case c = 5%, which we have probed up to a repulsion of ≈ 500k B T . The observed behavior is strikingly different from the corresponding bulk one, where the Hertzian potential was found to be valid only up to a few k B T s [38]. Indeed, in bulk the distinction between core and corona imposes to consider different kind of interactions, depending on the probed distances [8], that would describe different inner regions of the particle with changing elastic properties. Instead, at an interface, the microgel behaves as if the polymer network were more homogeneous and uniform: the 2D Hertzian behavior persists even in the core region, which here corresponds to r 30σ m for c = 5%, as also reported in Fig. 2(a).
By fitting the calculated potential with Eq. 1, we can obtain the effective diameter σ eff of the flattened microgel and its Young's modulus Y . Interestingly, the latter quantity can be also directly estimated from the fit of the calculated potentials, at odds with the corresponding 3D case where two elastic parameters, namely Y and the Poisson's ratio ν, are contained in the Hertzian prefactor [38]. The resulting fit parameters are shown in Fig. 2(b-c). In particular, the effective diameter is found to be very close, at all c, to the microgel extension σ ext , that can be estimated . Elastic moduli at the interface and in the bulk. Bulk modulus K, shear modulus G, Poisson's ratio ν and Young's modulus Y for the same microgel topology at the interface (full symbols) and in bulk (empty symbols) with explicit solvent as a function of c. In the last row, the theoretical results for Y are also compared to the ones obtained from the effective potential fits with the Hertzian model (Eq. 1), also reported in Fig. 2(c). K, G and Y are in units of kBT /σ 3 m to appropriately compare bulk and interface moduli, where the latter are divided by the thickness of the shell at the interface (≈ σm); ν is dimensionless. Error bars estimated from the fits of P (J) for K and P (I) for G (see Methods) are propagated in the calculation of ν and Y .
by taking opposite edges of the microgel on the interfacial plane [24], as displayed in Fig. 2(b). The slight underestimation of σ ext as compared to σ eff is associated to the fact that effective interaction calculations are also sensitive to the outer dangling chains. This information is partially lost by averaging over the distance of all opposite sites on the plane of the interface. As expected, the extension of the microgel at the interface decreases as a function of c in agreement with experiments [24], since softer microgels deform more strongly, and hence spread more at the interface. The corresponding values of Y are reported in Fig. 2(c), showing that higher crosslinking leads to stiffer networks, following expectations and in agreement with findings for microgels in bulk [38].

III.2. Elasticity theory calculations
The estimate of the Young's modulus extracted from the fit can be compared to the one obtained through the use of elasticity theory in 2D.
In this framework, one can evaluate the area and shape fluctuations of the microgel on the plane of the interface, writing the elastic energy U as a function of the two strain invariants of the strain tensor [55,56]. In this case, we write U according to the phenomenological Mooney-Rivlin theory, that is known to be valid also beyond the linear elastic regime. Within this theoretical approach, we can calculate all the elastic moduli of the microgel from equilibrium simulations, as previously done for microgels in bulk [38]. In particular, the moduli refer to the two-dimensional projection of the microgel on the interface, assuming that these are dominated by the corona fluctuations. In order to compare with the corresponding bulk properties, we also perform a similar procedure in 3D for the same microgel topologies in the presence of explicit solvent (see Ref. [38] for the implicit solvent treatment). However, bulk and interfacial moduli are naturally given in different units, so we adopt the plane-stress approximation for the 2D moduli. In this way, we assume that the stress normal to the interface is zero [57,58], a legitimate assumption for two-dimensional objects. Hence, we consider the small thickness of the microgels normal to the interface to be roughly comparable to the monomer size (≈ σ m ), and divide the obtained 2D moduli by this length. We can finally convert them into the correspoinding 3D moduli for very thin three-dimensional objects using the relations reported by Torquato [59] for plane-stress conditions. More details on these calculations are provided in the Methods section and in the Supplemental Material.
The resulting elastic moduli are reported in Fig. 3 as a function of the crosslinker concentration both for microgels at interfaces (left panels) and in bulk (right panels). Overall, we observe a monotonic increase of K, G and Y as a function of c, while ν remains nearly constant. These trends are preserved both in the bulk and at the interface. We stress that our two independent estimates of the Young's modulus, namely the one provided by the Mooney-Rivlin theory and that obtained by the 2D Hertzian fitting, also reported in Fig. 3, are consistent with each other. We highlight in this way how the single-particle properties of a microgel at an interface are fully reflected in the multiparticle behavior. The most striking result of this analysis is the fact that all three moduli at the interface are significantly larger, approximately by one order of magnitude, than their respective bulk counterparts. As for the Poisson's ratio, even though we find similar values in both cases, it should be noted that its upper limit in 2D is 1.0 while in 3D is 0.5 [60].
These findings provide a robust evidence of the reduced flexibility of the microgels at a liquid-liquid interface, an issue that up to now has been either extracted from indirect results or sometimes attributed  to charge effects [61,62]. Instead, we directly prove that it is entirely attributable to the presence of the interface, where microgels assume a much more stretched configuration in contrast to their standard arrangement in bulk. Under these conditions, microgels are much more resistant to deformation. Indeed, the corona is completely extended and confined at the interface with the polymer chains being much less responsive to external forces than in bulk, while still minimizing the surface tension. We further note that no available experimental results have so far reported the lateral elastic response of the microgels on the interfacial plane but rather the perpendicular one over a solid substrate [63]. The former is supposed to be the relevant one for the formation of thin microgel layers or for pattern formation on surfaces [20].

III.3. Multi-particle dynamical response
The level of coarse-graining adopted up to now has allowed us to describe how the properties of single constituents affect their mutual interactions. Now we go one step further by investigating the multi-particle behavior, i.e. the condition where many microgel particles interact on the interfacial plane. To shed light on this aspect, we simulate a system of particles whose interaction potential is the one we extracted previously, that is the 2D Hertzian potential. In this way, by further coarse-graining our system, we are able to assess for instance the dynamical response of microgels that are adsorbed on the interfacial plane.
The research on the phase behavior of soft colloids has recently gained much interest: being the archetype potential to describe interactions among elastic particles, the Hertzian phase diagram has been studied both in three [42,64] and in two dimensions [65,66]. In the latter case, however, the inves-tigations that have been carried out were limited to a change in the value of the exponent of the well-known 3D Hertzian without considering that a variation in the dimensionality of the problem implies a change in the functional form itself. Indeed, the logarithmic correction arising in Eq. 1 cannot be properly captured by a simple variation in the Hertzian exponent.
We perform Langevin Dynamics 2D simulations of Hertzian particles for different area fractions φ and varying the strength of the 2D Hertzian A = πY σ 2 eff /(2 ln 2), which corresponds to the r → 0 limit in Eq. 1. In order to have access to the dynamical response, we avoid crystallization by introducing polydispersity in the system (see Methods). Figure 4(a) reports the self-diffusion coefficients D extracted from the long-time behavior of the mean-squared displacements of the effective microgels for different values of φ and A. We consider the system to fall out-ofequilibrium on the simulation timescale when D decreases by roughly three orders of magnitude with respect to its low-density value. Thus, we consider the system to be arrested for D 2 × 10 −5 .
Importantly, we reveal the onset of two clear reentrant melting phenomena where the diffusivity, at first, decreases leading to the formation of a glassy system and then it grows again. This takes place primarily for φ 1.5 with a local maximum emerging at φ ∼ 1.9. For higher densities, after a further slowdown, the system re-fluidifies again acquiring a finite diffusion coefficient. Interestingly, at the new local maximum appearing for φ ∼ 2.5, the value of D is even larger than for the previous one. We stress that similar reentrant features in the dynamics have long been predicted in the three dimensional version of the Hertzian potential [42] and in extensive simulations of monomer-resolved single-chain nanoparticles [67], but they have never been found in experiments of soft [68] and ultrasoft colloids [69].
Previous works have shown that one can estimate the locus of the glass transition by monitoring the socalled iso-diffusivity (iso-D) lines [70][71][72], along which D remains constant. Importantly, it has been shown that the iso-D lines always maintain, for not too large values of the probed D, the same shape as the ultimate line of arrest. Thus, by extrapolating to the D → 0 limit, it is possible to locate the glassy region of a system. By taking a set of different isodiffusivity lines in Fig. 4(a), we draw the corresponding fluid-glass state diagram for the 2D Hertzian model, shown in Fig. 4(b). We notice that for the present system a fluid-like region persists at high densities for βA 1100.
As compared to the cases in which a generalized version of the Hertzian potential has been employed [17,66], our treatment allows to trace back the conditions for which this behavior can actually be observed experimentally. Indeed, the functional form of the potential we employ is physically meaningful, being connected to the Young's modulus and to the effective diameter of the single-particle microgels under investigation. We thus consider microgels whose combined spreading and elastic properties at the interface would fall into this range. To do this, we need to focus on microgels with reduced size. In fact, a reduction of the particle diameter can strongly affect the value of the Hertzian strength, which depends quadratically on it. Hence, we perform additional simulations of singleparticle microgels made of 2000 and 3000 monomers, besides those with ≈ 5000 monomers. In order to avoid long computational times for the calculation of the effective interactions, we directly determine the Hertzian strengths via elasticity theory calculations and by measuring σ ext for single particles with different sizes and crosslinker concentrations at the interface.
We report the estimated strengths as a function of c in Fig. 5 and find that soft and small microgels have an Hertzian strength that falls in the range where a reentrant behavior of the dynamics is present, according to the phase diagram in Fig. 4(b). We also confirm that the value of the Young's modulus does not show a strong size dependence, especially for c = 3% and 5% (see Supplemental Material), in qualitative agreement to experimental findings on microgels of different sizes [73][74][75]. Hence, from this analysis, we conclude that highly crosslinked microgels will always have glassy dynamics at the interface, independently on their size. Instead, low-crosslinked microgels whose Young's modulus at the interface is around 0.1 − 0.3k B T /σ 2 m and whose extended size is between ≈ 35−50σ m are expected to show a reentrant dynamics. Using the mapping performed in previous comparisons with experiments [24,37], this would correspond to microgels with hydrodynamic diameters in bulk 200 nm, well within the experimental range, where also capillary effects are less relevant. Therefore, adsorbed microgels constitute a realistic model system to experimentally investigate the presence of a reentrant dynamics, long postulated in the 2k 3k 5k realm of soft colloids. It is also instructive to think where this regime can be observed in terms of compression isotherms to which experiments typically refer to. From the present calculations, we estimate that the value of the area fraction is reduced by about a half as compared to the corona-corona contact at low densities. Even though these are not too high compressions [32], a number of critical issues may emerge and these are ultimately linked to the real-time visualization of the microgels at the interface. Currently in fact most of the studies are performed ex-situ by means of AFM on silica wafer or similar techniques. In these cases, despite the structural information on the microgel layer being reliable, no information on the dynamics of such particles can be observed. In this regard, we believe that our predictions will stimulate experimental work to confirm the predicted dynamical behavior for microgels at interfaces.

IV. CONCLUSIONS AND PERSPECTIVES
In summary, in this work we have provided the first numerical estimate of the two-body effective interaction potential of microgel particles adsorbed at an interface. Its simple functional form reveals that such particles interact like effective elastic disks with a Young's modulus that increases with the crosslinker concentration. Notably, the values of the elastic moduli at the interface, after appropriate rescaling, are found to be roughly one order of magnitude higher than the one measured in bulk, as also confirmed by elasticity theory calculation of single microgel particles. This can be tentatively attributed to the dominant effects of the interfacial tension, which controls the response of the polymer network to an external stress, making it much stiffer and resistant to deformation with respect to the same network in good solvent conditions. This result has profound consequences on the properties of a generic interfacial assembly of soft colloids, not limited to microgel particles. Indeed, we expect that the reduced mobility of the polymer chains and their enhanced stiffness should be taken into account in the development of novel materials that rely on deformable constituents of any kind. As demonstrated by our results, this effect should be expected at interfacial conditions with large surface tensions, independently of the presence of intrinsic charges in the material or in the fluids. Furthermore, our study clearly demonstrates that the knowledge that is gained on the bulk properties of soft colloids cannot be transferred to the interface, which should be considered as a completely separate case.
From a more fundamental perspective, we have numerically shown the existence of soft colloidal particles that behave as Hertzian disks even at short separations, beyond the small deformation regime. The extensive analysis of the multi-particle dynamics has further evidenced the emergence of reentrant melting transitions, where the system behaves as an ergodic fluid up to very large densities, well above individual particles contact, sometimes loosely called jamming. Experimentally, small (nano-sized) soft microgels appear to be the ideal candidates to verify our theoretical predictions, as indicated by the values of the Young's modulus and of the interfacial extension at which the reentrance is observed. In addition, small colloids are the least likely to experience capillary attractions at the interface, and hence will behave more similarly to the ones we have simulated.
It will be important in the future to extend this study to crowded configurations to investigate the validity of the present results at considerably high packing fractions where faceting phenomena can be foreseen. In the latter case, many-body effects should be taken into account for more refined treatments, where the microgel-microgel interaction may not be pairwise additive. Our approach provides, in this respect, a first step towards a comprehensive description of microgel interactions at a microscopic level. Similar considerations should be extended to microgels in bulk conditions, for which high-density states still require appropriate theoretical assessment. The analysis can be further broadened to other microgel topologies that have recently gained considerable attention, such as hollow [76], ultra-low-crosslinked [23] or anisotropic ones [77].
All in all, our study opens the way for the investigation of microgels at the interface as a simple realization of 2D elastic particles. We expect that the evidence reported here will have important consequences on the study of two-dimensional elastic objects at the fundamental level [78,79] and for the clever design of composite materials [80][81][82].
V. METHODS

V.1. Modeling and simulation details
Microgels are assembled in a disordered fashion starting from an ensemble of two and four-folded patchy particles in a spherical cavity. Bivalent and tetravalent particles mimic respectively NIPAM (N-Isopropylacrylamide) monomers and BIS (N,N'methylenebisacrylamide) crosslinkers in a chemical synthesis [45]. The topology of the polymer network, whose monomers have diameter σ m , is then fixed by means of a classical bead-spring model [43] that amounts to the Weeks-Chandler-Andersen (WCA) potential for non-bonded monomers, and a sum of the WCA and the Finitely Extensible Nonlinear Elastic (FENE) potentials for bonded ones: with k F = 15 a dimensionless spring constant and R 0 = 1.5 the maximum extension of the bond. To reproduce in a realistic fashion the form factors and the density profiles of realistic microgels, we add a designing force on the crosslinkers as described in Ref. [37]. In the present work, we study microgels with ≈ 5000 monomers and three different crosslinkers concentrations c, namely 3%, 5% and 10%. The radius of the confining sphere into which microgels are assembled is set to 25σ m . We also evaluate the elasticity of smaller microgels with ≈ 2000 and 3000 particles assembled in the same way, maintaining the same internal monomer density. According to previous works [24,46], the solvent is modeled within the Dissipative Particle Dynamics (DPD) framework [83]. The interactions are described by three forces, conservative F C ij , dissipative F D ij and random F R ij , of form: where r ij = r i − r j , with r i the position of particle i, r ij = | r ij |,r ij = r ij /r ij , r c the cutoff radius, v ij = v i − v j with v i the velocity of particle i, a is the maximum repulsion between two particles, θ ij is a Gaussian random number with zero mean and unit variance and ξ is the friction coefficient [83]. To ensure that Boltzmann equilibrium is reached, w D (r ij ) = [w R (r ij )] 2 and σ 2 R = 2ξk B T with k B the Boltzmann constant and T the temperature.
In order to reproduce a water/hexane (w/h) interface, we choose a ww = a hh = 8.8, a hw = 31.1 whereas, for the monomer-solvent interactions a mw = 4.5 and a mh = 5.0. The cut-off radius is always chosen to be r c = 1.9σ m and the reduced solvent density ρ DPD = 4.5. Depending on the microgel size, up to ≈ 750000 particles are inserted in the simulation box. To analyze the elastic properties of microgels in bulk, we also run bulk simulations with explicit solvent. In there, the solvent-solvent parameters are not varied with respect to interfacial simulations, while a ms = 1.0. A more detailed discussion on how these parameters have been determined can be found in Ref. [24].
The phase behavior of the 2D Hertzian potential is assessed by means of molecular dynamics simulations in two dimensions with 5000 particles interacting via Eq. 1. To avoid crystallization, we set the polydispersity to p = 0.2. We analyze a range of area fractions φ = π 4 σ 2 m ρ, with ρ the number density, that goes from 0.8 to 2.8 for a Hertzian strength, defined as A = πY σ 2 eff /(2 ln 2), that goes from 220 to 1150k B T . To determine the glass region in the 2D phase diagram, we run simulations for ∼ 2 × 10 7 timesteps and we calculate the mean-squared displacement ∆r 2 of the particles, extracting the long-time self-diffusion coefficients D: where t is the simulation time. Since we are only interested in providing a state diagram assessment, we do not perform an extensive characterization of the glassy dynamics of the system and we just monitor the onset of non-ergodicity within the timescale of our simulations. We attribute this condition to state points where we find D 2.5 × 10 −5 , roughly three orders of magnitude lower than the corresponding lowdensity value. Under these conditions, the system has become so slow that aging is present within our simulation time window.
All simulations are carried out using the LAMMPS simulation package [84]. The equations of motion are integrated with a velocity-Verlet algorithm. The reduced temperature T * = k B T / is always set to 1.0 via the DPD thermostat in explicit solvent simulations [46] and via the Langevin thermostat for simulations in two dimensions. Length, mass, energy and time are given in units of σ m , m, and mσ 2 m / , respectively. DPD repulsion parameters a are in units of /σ m .

V.2. Calculation of the effective interaction potential
The two-body effective potential between the microgels at the interface is evaluated by means of the Umbrella Sampling technique in explicit solvent [47][48][49]. This method allows to uniformly sample all distances between the centers of mass of the microgels by adding a harmonic potential between them. For each sampled window i, we evaluate the probability distribution P (r, ∆ i ) of finding the microgels' centers of mass at distance r given the equilibrium length of the spring ∆ i . The final probability for the entire range of explored distances is obtained by first removing the contribution of the bias potential and by subsequently merging P (r, ∆ i ) into P (r) for all the windows via a least-square method. Finally, the potential of mean force V eff is retrieved knowing that where C is such that V eff (r → ∞) = 0. The major drawback of studying the effective interactions with an explicit solvent model is the computational cost required to carry out the simulations: about two months on about 80 CPUs are needed to investigate a range from 20 to 30∆ i .

V.3. Assessment of the elastic moduli
The elastic energy U of a two-dimensional object can be written as a function of the invariants J and I of the strain tensor as U (J, I) = U 0 + W (J) + W (I) (9) where U 0 is the energy of a reference configuration that is taken as the average ellipse adopted by equilibrium configurations of the microgels at the interface. Its semi-axes s 1 and s 2 are obtained in turn by the gyration tensor built via the two-dimensional convex hull on the plane of the interface. We approximate W with the corresponding potentials of mean force with X = J, I, P (X) the respective probability distribution and D an arbitrary constant. These potentials can then be fitted to appropriate functions with γ = 2 when X = J and γ = 1 when X = I, to obtain M J and M I . In the Supplemental Material, we report, as an example for c = 10%, the simulation outcomes and their relative fit both for the microgel at the interface and in bulk. The elastic moduli are then readily obtained as with A = πs 1 s 2 . Y and ν only depend on K and G as [60] Similar expressions can be derived for the 3D case and can be found, for instance, in Ref. [38] (see also below).
The particular choice of W as a function of J and I depends on the specific elastic model employed. Here, we considered the Mooney-Rivlin model for which the elastic energy reads [55,56] U (J, We have further checked that the obtained results do not crucially depend on the specific form of the employed W . To this aim, we also employed the linear elastic model (Hookean) [55,85] and the Saint-Venant-Kirchhoff model [86,87], finding results for the moduli, particularly the Young's modulus, that are very close to the ones presented in the main text. They display the same increase with respect to the bulk model and a similar monotonic increase with c.
To convert the 2D moduli into 3D ones, we consider that for two-dimensional objects the stress normal to the interfacial plane is zero. Under these conditions, there exist relations to convert 2D moduli into 3D ones, by assuming that the 2D object has a given (small) thickness h: where X (3) indicates the converted 3D moduli (in units of k B T /σ 3 m ) from the 2D results X (2) (in units of k B T /σ 2 m ) with X = G, Y, K. Also, we have that [59] ν (3) = ν (2) .
In our case, we consider h to be roughly equal to the monomer size, σ m , as in the outer shells chains do not pile up, but remain confined to the interfacial plane, providing the dominant contribution to the elastic response of the microgels.
Due to the high computational cost of the calculations of the effective potential, we limit the analysis in the main text to a single microgel topology. However, we checked that the presented results are robust in their main conclusions by analyzing a different microgel topology for the case c = 5%. We report in Fig. S1 the interaction potential for such a different topology, limiting the analysis, again due to the high computational cost, to small compression regimes. Also in this case, we confirm the validity of the fitting procedure with the bidimensional Hertzian potential (Eq. 2 in the main text). The extracted fitting parameters are σ eff = 54.5σ m and Y = 0.26k B T /σ 2 m , in good qualitative agreement with the ones of the microgel presented in the main text.

FURTHER CONSIDERATIONS ON THE EFFECTIVE POTENTIAL
It is interesting to compare the effective potential between microgels at interfaces to a gaussian functional form, that was used to describe brush-coated spherical nanoparticles in bulk [88,89] and at an interface [90]. Such a simple model can be written as where b, d, e are fitting parameters and r is the distance between the centers of mass of the particles. From a structural perspective, the conformation that microgels retain at interfaces may resemble the one of such particles, given the presence of extended, flexible polymer chains surrounding a more compact core. For polymer brushes, there exists a scaling theory for the fitting parameters b and d [89]. Surprisingly, the functional form in Eq. 21 was found to describe the calculated interactions for these systems quite well. Nonetheless, there should not be any physical reason for this framework to be applicable to our system. In addition, the gaussian fit does not account for any deformation of the polymer at the interface, being developed for 3D bulk systems. We report in Fig. S2 our calculated potentials and the corresponding fits with Eq. 21 and with the 2D Hertzian model described in the main text. We find that the latter agrees much better with data also at large distances between the microgels, while the gaussian form fails in this regime. Although this is the region in which data are most affected by statistical noise being the probed energy of the order of a few k B T s, the gaussian fit would give rise to a potential which tends to zero at distances that are clearly noncompatible with the dimensions of the microgel particles analyzed here (see also the comparison with σ ext   in Fig. 2 of the main text). While we could think of operating the gaussian fit in a reduced region of distances, i.e. only for short ones, it is important to stress that the parameters that we would extract from such fits cannot be related to any physical feature of our system. For the sake of completeness, the fitting parameters for the gaussian functional form are reported in Table I, where it is evident that in the case of parameter e, we cannot even identify a clear trend as a function of the crosslinker concentration. We thus conclude that the 2D Hertzian description presented in the main text is the most appropriate to treat the effective interactions between microgel particles at an interface.

ELASTIC MODULI
As also explained in the Methods section of the main text, the elastic energy of the microgel can be written as a function of a reference configuration energy and of W (X) with X = J, I, being J and I the invariants of the strain tensor [55,56]. W can be approximated with the potentials of mean force W (X) = −k B T ln P (X) + C with P (X) the respective probability distribution and C an arbitrary constant. This can be fitted to functions of form f (X; M X , X 0 , γ, C) = M X (X − X 0 ) γ + C with γ = 2 when X = J and γ = 1 when X = I, to obtain M J and M I . Similar considerations apply to the 3D case, for microgels in bulk [38]. In Fig. S3, we report, as an example for c = 10%, the simulation outcomes and their relative fit both for the microgel at the interface and in bulk. From these, all the elastic moduli are readily obtained. In Fig. S4, we show the dependence of the Young modulus as a function of the size of the microgel. Except for the smallest microgels with c = 10%, we observe only a slight dependence on the size of the microgel at fixed crosslinker concentration. This is in qualitative agreement with experimental findings [75].

MULTI-PARTICLE DYNAMICAL RESPONSE
We show in Fig. S5, for a representative case, the mean-squared displacement of particles interacting via