# Dynamics of Nanoparticles in Polydisperse Polymer Networks: from Free Diffusion to Hopping

Valerio Sorichetti, Virginie Hugouvieux, Walter Kob

► To cite this version:

Valerio Sorichetti, Virginie Hugouvieux, Walter Kob. Dynamics of Nanoparticles in Polydisperse Polymer Networks: from Free Diffusion to Hopping. Macromolecules, 2021, 54 (18), pp.8575-8589. 10.1021/acs.macromol.1c01394 . hal-03358744

HAL Id: hal-03358744

https://hal.inrae.fr/hal-03358744v1

Submitted on 11 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamics of nanoparticles in polydisperse polymer networks: From free diffusion to hopping

Valerio Sorichetti,*,†,‡,¶ Virginie Hugouvieux,*,¶ and Walter Kob*,§

†Laboratoire de Physique Théorique et Modèles Statistiques (LPTMS), CNRS, Université Paris-Saclay, F-91405 Orsay, France

‡Laboratoire Charles Coulomb, Univ. Montpellier, CNRS, F-34095, Montpellier, France

¶IATE, Univ. Montpellier, INRAE, Institut Agro, F-34060, Montpellier, France

§Laboratoire Charles Coulomb (L2C), Univ. Montpellier, CNRS, F-34095, Montpellier, France

E-mail: valerio.sorichetti@universite-paris-saclay.fr; virginie.hugouvieux@inrae.fr; walter.kob@umontpellier.fr

## Abstract

Using molecular dynamics simulations we study the static and dynamic properties of spherical nanoparticles (NPs) embedded in a disordered and polydisperse polymer network. Purely repulsive (RNP) as well as weakly attractive (ANP) polymer-NP interactions are considered. It is found that for both types of particles the NP dynamics at intermediate and at long times is controlled by the confinement parameter $C = \sigma_N/\lambda$, where $\sigma_N$ is the NP diameter and $\lambda$ is the dynamic localization length of the crosslinks. Three dynamical regimes are identified: i) For weak confinement ($C \lesssim 1$) the NPs can freely diffuse through the mesh; ii) For strong confinement ($C \gtrsim 1$) NPs proceed by means of activated hopping; iii) For extreme confinement ($C \gtrsim 3$) the mean squared displacement shows on intermediate time scales a quasi-plateau since
the NPs are trapped by the mesh for very long times. Escaping from this local cage
is a process that depends strongly on the local environment, thus giving rise to an
extremely heterogeneous relaxation dynamics. The simulation data are compared with
the two main theories for the diffusion process of NPs in gels. Both theories give a
very good description of the C-dependence of the NP diffusion constant, but fail to
reproduce the heterogeneous dynamics at intermediate time scales.

# 1 Introduction

When a nanoparticle (NP) is embedded in a polymer network, its dynamics can slow down dramatically.¹,² Understanding what factors govern this slowing down is of primary importance in many fields, such as material science (e.g. with application to thin films³⁻⁵ and polymer-based sensors³,⁶,⁷), biophysics,⁸⁻¹⁶ and medicine, in particular for applications to drug delivery.⁶,¹⁷,¹⁸ While in recent years the diffusion of NPs in polymer solutions and melts has been the subject of numerous theoretical¹⁹⁻²³ and simulation studies,²⁴⁻³⁵ only few investigations have dealt with the problem of NP diffusion in permanently crosslinked networks, despite its importance in many applications.³⁶⁻⁴⁴ In some of the earliest simulation studies, this network was simply modeled as an array of fixed obstacles,³⁶ which is clearly a far cry from a physically realistic description. Other authors have included connectivity and flexibility in the network model, but most of them considered only regular structures, in which the crosslinks are placed on the vertices of a regular lattice and connected either by chain segments³⁷,⁴⁵ or directly by springs.³⁸⁻⁴⁰,⁴² In the latter case, since there is no actual strand connecting the crosslinks, strand dynamics and entanglement effects are not accounted for. Moreover, real-life networks such as hydrogels,⁴⁶ vulcanized rubbers,⁴⁷ or networks produced by electron irradiation,⁴⁸ are often disordered and polydisperse, with a continuous distribution of strand lengths, properties that lead to an additional complexity in the dynamics of the NP. Recently, a small number of simulation studies adopting more realistic models for the network has been published,⁴¹,⁴³,⁴⁴ however none of these studies
has, to the best of our knowledge, taken explicitly into account the effect of disorder and polydispersity.

Also analytical studies dealing with NP dynamics in permanently crosslinked networks are rather scarce. 1,2 In entangled polymer liquids, the relaxation dynamics of particles of size larger than the tube diameter$^{49}$ *d* can proceed through the release of the entanglements (constraint release$^{22}$), which happens on time scales of the order of the disengagement (repetition) time, $\tau_d \propto N^{3.4}$, 49 with *N* the degree of polymerization of the chains. In polymer solids containing irreversible crosslinks, like dry networks (e.g. rubbers) and gels, the constraint release mechanism is completely turned off and the only available process for the motion of large NPs is the one of hopping, i.e., activated motion triggered by local fluctuations of the entanglement/crosslink mesh. 1,2 Dell and Schweizer$^1$ have developed a theory of hopping based on a combination of a nonlinear Langevin equation and PRISM$^{50}$ theory, 51 showing that the quantity which controls the NP dynamics is the so-called *confinement parameter*, i.e. the ratio $C = \sigma_N/d$ between the NP diameter $\sigma_N$ and the effective tube diameter *d*, resulting from both crosslinks and entanglements. If the effect of the entanglements can be neglected one has $d \approx \xi$, where $\xi$ is the average mesh size of the network. 2 The conclusion that *C* is the parameter controlling the NP dynamics has also been reached by Cai *et al*. using scaling theory. 2 We note, however, that the two approaches predict qualitatively different behaviors for the relevant dynamical quantities, such as the NP diffusion coefficient as a function of *C*, see below for details. The importance of the confinement parameter for the description of the diffusion of NPs in polymer networks has also been confirmed in experiments$^{14,15,52-55}$ and simulations.$^{42-44}$ However, even recent simulations$^{43,44}$ have not explored the strong confinement regime, $C \gtrsim 3$, due to the extremely slow dynamics that characterizes it and hence the dynamics of the NP in this range of parameters is at present not known.

In this work, we present a simulation study of NP diffusion in polymer networks which are both disordered and polydisperse. In particular, we probe for the first time in simulations
the strong confinement regime, considering confinement parameters up to $C \simeq 4$, for which NP motion is dramatically slowed down. The remaining part of the paper is organized as follows: In Sec. 2, the model and the simulation method are presented. In Sec. 3, we analyze the structural properties of the network and of the NPs for different NP diameters at low NP concentration. In Sec. 4, we present the NP dynamics and in particular their diffusion coefficient and van Hove function, showing that for $C > 1$, NP diffusion proceeds through hopping motion. We conclude with a summary in Sec. 5.

## 2 Model and simulation method

We performed molecular dynamics (MD) simulations of spherical NPs embedded in a polydisperse, disordered and permanently crosslinked polymer network. The network is generated following the procedure described by Gnan *et al.*,⁵⁶ initially developed for the simulation of microgels,⁵⁷,⁵⁸ which is based on the self-assembly of particles with limited valence (“patchy” particles): $N_m = 4 \times 10^5$ particles (monomers) are placed in a cubic box of volume $V$ with periodic boundary conditions, and thus the monomer density is $\rho_m = N_m/V$. Of these particles, $N_{cl} = cN_m$, with $c = 0.1$, play the role of crosslinks, in that they can form three bonds, whereas the $N_m - N_c$ others can only form two bonds (bivalent particles). These bivalent particles have patches on the opposite site of the particles, while the crosslinks have three patches forming 120 degrees with each other. Two crosslinks are forbidden to bind to each other, whereas a bivalent particle can bind to any other particle. A NVT MD simulation is started and it is stopped when at least 99% of all the possible bonds are formed and a percolating network is generated. At this point, all the monomers which are not part of the percolating network (at most 0.4% of all the monomers in the systems considered here) are removed. (Although for the sake of computer time we chose to stop the reaction before reaching the fully-bonded ground state of the system, reaching this state is in principle possible by making a greater computational effort.) The system obtained from this procedure
contains only a small number of dangling ends, i.e., it is an almost perfect network. We find that the number of monomers forming dangling ends is less than 2.5% of the total. However, since a dangling end of length $n$ has a relaxation time which growth exponentially with $n$,⁵⁹,⁶⁰ we remove recursively all the dangling ends in the system, so that at the end a fully-bonded network is obtained. We note that in this procedure the value of $N_m$ slightly decreases as well as the fraction $c$ of crosslinks since a crosslink becomes a bivalent particle whenever a dangling end attached to this crosslink is removed; however, given the small number of dangling ends, these changes can be considered to be negligible.

Previous studies have found that the networks generated using this procedure have a chain length distribution $p(n)$ which decays exponentially in $n$, Ref.,⁶¹ and is given by the Flory-Stockmayer expression,⁶²,⁶³

$$p(n) = \frac{N_m - N_c}{\langle n \rangle^2} \left(1 - \frac{1}{\langle n \rangle}\right)^{n-1}, \quad (1)$$

where $n$ is the chain length, i.e., the number of beads between two crosslinks, see Supplementary Material, Sec. S2.1. We note that $p(n)$ is independent of the monomer density $\rho_m$ in the $\rho_m$-range here considered, as it is to be expected given the equilibrium nature of the assembly procedure.⁵⁶,⁶¹ The mean chain length $N_x$ is obtained as $N_x \equiv \langle n \rangle = 2(c^{-1} - 1)/3$,⁶⁴ which for $c = 0.1$ gives $N_x = 6$ (the actual value is slightly larger, since the removal of the dangling ends makes that $c$ decreases slightly, as discussed above).

Once the dangling ends are removed, the topology is frozen and the interaction potential between the particles is changed from that of Gnan et al.⁵⁶ to the Kremer-Grest potential,⁶⁵ in which all monomers interact via a Weeks-Chandler-Andersen (WCA) potential,⁶⁶
$$ U_{mm}(r) = \begin{cases} 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} + \frac{1}{4} \right] & r \le 2^{1/6}\sigma \\ 0 & \text{otherwise.} \end{cases} \qquad (2) $$

In addition, bonded monomers interact via a finite extensible nonlinear elastic (FENE) potential,

$$ U_{\text{bond}}(r) = -\frac{kr_0^2}{2} \ln \left[ 1 - \left( \frac{r}{r_0} \right)^2 \right], \qquad (3) $$

where $k = 30\epsilon/\sigma^2$ and $r_0 = 1.5\sigma$. With this choice of parameters the minimum of the potential is at $r_b = 0.961$. The combined effect of the FENE and the WCA potentials prevents the chains from crossing each other at the thermodynamic conditions considered here.⁶⁵ In the following, all quantities are given in Lennard-Jones (LJ) reduced units. The units of energy, length and mass are thus, respectively, $\epsilon$, $\sigma$ and $m$, where $\epsilon$, and $\sigma$ are defined by Eq. (2) and $m$ is the mass of a monomer. The units of temperature and time are, respectively, $T^* = \epsilon/k_B$ and $\tau^* = \sqrt{m\sigma^2/\epsilon}$, where $k_B$ is Boltzmann's constant, which we set equal to 1.

Once the network is generated, we embed in it spherical nanoparticles and perform MD simulations to anneal the system. We consider NPs of diameter $\sigma_N$ ranging from $\sigma_N = 1$ to $\sigma_N = 10$ (see Supplementary Material, Sec. S1). For a given set of parameters we use a single realization of each NP-enriched network, since the system is large enough to be self-averaging. Note that the network is the same for each value of $\sigma_N$. The interaction between monomers and NPs and between two NPs is given by an expanded Lennard-Jones potential:³³,⁴³,⁶⁷,⁶⁸
<figure><img src="image_2.png" /><figcaption>Figure 1: Left: Snapshot of the system for $\rho_{m0} = 0.290$, $\sigma_N = 6.0$ for the repulsive nanoparticles. NPs, crosslinks and bivalent monomers are represented, respectively, by red, blue and green spheres. For clarity, bivalent monomers are shown as transparent. Right: Close-up of a NP with the surrounding polymer mesh. Shown are all monomers within a distance $3\sigma_N/2$ from the NP center.</figcaption></figure>

$$ U_{N\alpha}(r) = \begin{cases} 4\epsilon \left[ \left( \frac{\sigma}{r - \Delta_{N\alpha}} \right)^{12} - \left( \frac{\sigma}{r - \Delta_{N\alpha}} \right)^{6} \right] + E_{N\alpha} & r \le \Delta_{N\alpha} + r_{N\alpha}^c \\ 0 & \text{otherwise,} \end{cases} \quad (4) $$

with $\alpha \in \{N, m\}$ and where for the NP-monomer interaction $\Delta_{Nm} = (\sigma_N + \sigma)/2 - \sigma = (\sigma_N - \sigma)/2$ and for the NP-NP interaction $\Delta_{NN} = \sigma_N - \sigma$. The addition of $E_{N\alpha}$ ensures that at the cutoff distance $\Delta_{N\alpha} + r_{N\alpha}^c$ the potential $U(\Delta_{N\alpha} + r_{N\alpha}^c)$ is continuous. The NP-NP interaction is purely repulsive ($r_{NN}^c = 2^{1/6}$), whereas for the NP-monomer interaction we consider both repulsive ($r_{Nm}^c = 2^{1/6}$) and attractive ($r_{Nm}^c = 2.5$) interactions. In the following, we will refer to attractive NPs as ANP and to repulsive NPs as RNP. We assume that the NPs have the same mass density as the monomers, $\rho_{mass} = 6m/\pi\sigma^3$, and therefore the mass of the NPs is $m_N = m(\sigma_N/\sigma)^3$. In Fig. 1 we show a snapshot of the system for $\sigma_N = 6$ and density $\rho_{m0} = 0.290$ (defined below) for the case of RNPs.

All the simulations are carried out using the LAMMPS software.⁶⁹,⁷⁰ The simulation box is cubic and periodic boundary conditions are applied in all directions. Initially, the centers of the NPs are inserted at random positions in the network. The possible overlaps are then
removed using the “fast push-off” method,⁷¹ that increases the NP diameter from 0 to 1; subsequently the NP size is gradually increased until the diameter reaches the value $\sigma_N$. After the NPs have reached the desired size, we perform an NPT run with Nosé-Hoover chains⁷² and allow the system to reach pressure $P = 0$ at temperature $T = 1.0$. During this NPT run, the box sides are coupled so that they fluctuate together, i.e. $L_x = L_y = L_z$. Once the system has adjusted to $P = 0$, we perform another short NPT run, during which the mean volume $V$ is measured. Subsequently we switch to the NVT ensemble by fixing the system’s volume to $V$ and perform an equilibration run before starting production at constant volume, i.e. the dynamics of the particles is not perturbed by unphysical volume fluctuations due to a barostat. Since the Nosé-Hoover thermostat does not produce a realistic dynamics,⁷² we switch to a Langevin thermostat, so that the force experienced by particle $i$ (monomer or NP) is given by⁷³

$$m_i \ddot{\mathbf{r}}_i = -\nabla_i U(\{\mathbf{r}_k\}) - m_i \gamma_i \dot{\mathbf{r}}_i + \zeta(t). \quad (5)$$

Here $\mathbf{r}_i$ is the position vector, $m_i$ the mass, and $U(\{\mathbf{r}_k\})$ is the total interaction potential acting on the particle, with $\{\mathbf{r}_k\}$ representing the set of coordinates of all the particles in the system. The second term on the right-hand side of Eq. (5) represents viscous friction, with $\gamma_i$ the friction coefficient. The term $\zeta$ is a stochastic force which represents the collisions with solvent molecules, and satisfies $\langle \zeta(t) \rangle = 0$ and $\langle \zeta_{\alpha}(t) \zeta_{\beta}(t') \rangle = 2m_i \gamma_i k_B T \delta_{\alpha,\beta} \delta(t-t')$, with $\zeta_{\alpha}$ its spatial components. The monomer friction coefficient is $\gamma_m = 0.1$, whereas the NP friction coefficient is $\gamma_N = \gamma_m(m\sigma_N/m_N\sigma) = \gamma_m(\sigma/\sigma_N)^2$, so that an isolated monomer and an isolated NP experience the same solvent viscosity.⁶⁷ Although with this thermostat hydrodynamic interactions are neglected, we expect that due to the slow dynamics these interactions are not relevant. The velocity Verlet algorithm is employed to integrate the equations of motion, and the integration time step is $\delta t = 0.006$. The duration of the equilibration run is
between $3 \times 10^4$ and $1.5 \times 10^6$ time units, depending on the system considered. We note that the relaxation time of the NPs is expected to increase at least exponentially in $\sigma_N/d$, where $d$ is an effective tube diameter resulting from both crosslinks and entanglements.¹,² For this reason, not all systems presented here have reached equilibrium, as it will be discussed below. Nevertheless, we decided to include also the results for these systems since they are affected only weakly by aging effects and hence are still instructive. The duration of the production runs is between $6 \times 10^4$ and $3 \times 10^6$, depending on the NP diameter and on the network density.

In order to avoid that the addition of the NPs influences significantly the structural and dynamical properties of the network and to study single-NP dynamics, the number $N_N$ of NPs embedded in the network is kept small. The value of $N_N$ is chosen in such a way that the total NP volume, $V_N \equiv \pi N_N \sigma_N^3 / 6$, is 2% of the mean volume of the neat system $V_0$, i.e., the NP volume fraction of the unrelaxed system is

$$ \phi_{N0} \equiv \frac{\pi \sigma_N^3 N_N}{6 V_0} = 0.02 \quad . \qquad (6) $$

To make sure that we are indeed probing the dilute limit, we have simulated also some systems with $\phi_{N0} = 0.005$, but, unless explicitly specified, the shown results are for $\phi_{N0} = 0.02$. No significant difference in the dynamical properties is found between these two values, confirming that the dilute limit is already reached at $\phi_{N0} = 0.02$. We consider three different starting networks, with monomer densities $\rho_{m0} = N_m/V_0 = 0.190, 0.290$ and $0.375$ in the neat state. We note that in general $V_0$ is different from the final volume $V$ (and therefore $\rho_0 \neq \rho$ and $\phi_{N0} \neq \phi_N$), since the addition of NPs can cause swelling or shrinkage (depending on $N_N$ and $\sigma_N$) of the network, as also observed in non-crosslinked nanocomposites.⁶⁷ This point is discussed in more detail in the Supplementary Material, Sec. S2.2, where one can also find the details of the simulated systems (Sec. S1).
## 3 Static properties

In this section we discuss the structural properties of the system, notably the radial distribution function, the static structure factor, and the pore size distribution.

### 3.1 Radial distribution function and static structure factor

<figure><img src="image_3.png" /><figcaption>Figure 2: Monomer-NP (a) and NP-NP (b) radial distribution functions for ANPs and RNPs embedded in the $\rho_{m0} = 0.290$ network.</figcaption></figure>

The simplest way to characterize the structure of the system is by means of the radial distribution function (RDF). In Fig. 2a we show for the systems with $\rho_{m0} = 0.290$ the monomer-NP radial distribution function, defined as$^{74}$

$$g_{mN}(r) \equiv \frac{N_m + N_N}{4\pi N_m N_N (\rho_m + \rho_N) r^2} \sum_{k=1}^{N_m} \sum_{j=1}^{N_N} \langle \delta(|\mathbf{r} + \mathbf{r}_k - \mathbf{r}_j|) \rangle, \quad (7)$$

where $\rho_m$ and $\rho_N$ are, respectively, the monomer and NP number density. Note that on the horizontal axis we plot $r - \sigma_N/2$ in order to remove the trivial $\sigma_N$ dependence of the curves. One recognizes that for repulsive NPs (RNPs) as well as for attractive NPs (ANPs), the data collapse approximately onto a master curve. For the ANPs, $g_{mN}(r)$ displays a main peak at $r \simeq (\sigma_N + 1)/2 = r_1$, which corresponds to the first shell of monomers touching the
NP, followed by a smaller peak at $r \approx (\sigma_N + 3)/2 = r_2$, corresponding to the second shell. Since this is not a dense system, the presence of these peaks demonstrates that, as expected, monomer-NP contacts are favored. One can see that the height of the main peak increases slightly if $\sigma_N$ is decreased: This can be explained by the fact that the total surface of the NP increases with decreasing $\sigma_N$, making that more monomers can touch this (energetically favorable) surface, thus increasing the height of the nearest neighbor peak. Furthermore this attraction has the effect that neighboring strands are pulled closer to each other, making that on overall the system shrinks (see Supplementary Material, Sec. S2.2).

For repulsive NPs, $g_{mN}(r)$ is significantly smaller than 1 for $r_1 < r < r_2$, in agreement with the fact that monomer-NP contacts are unfavorable. A consequence of this correlation hole is that the effective radius of repulsive NPs is larger than that of ANPs with the same $\sigma_N$, as we will discuss below, and hence the mobility of the ANP is higher than the one of the RNP, see next section. For the other two networks ($\rho_{m0} = 0.190$ and 0.375) we find qualitatively similar results (see Supplementary Material, Sec. S2.3).

In Fig. 2b we show for the same systems the NP-NP radial distribution function, which is obtained by setting in Eq. (7) $m = N$. For the RNPs, the RDF displays a pronounced single peak at the contact distance $r_c \equiv \sigma_N + (2^{1/6} - 1) = \sigma_N + 0.122$.⁶⁷ From the rather large value of $g_{NN}(r_c)$ and the absence of further peaks one can deduce that the RNPs tend to form clusters, likely by filling the largest holes in the mesh. The height of this peak increases with increasing $\sigma_N$, which is simply due to the fact that the NP density is lower, at constant NP volume fraction, for larger NPs (see also Supplementary Material, Sec. S2.3). For the ANPs, one finds at $r = \sigma_N + 0.122$ a weak shoulder and a small peak at $r = \sigma_N + (2^{7/6} - 1) = \sigma_N + 1.245$, which corresponds to a configuration in which two neighboring NPs are separated by a polymer strand.⁶⁷ This is a consequence of the fact that the NPs are well dispersed, and that each NP is surrounded by a layer of polymers. Such local structures are expected, since it is known that weak attractive polymer-NP interactions lead to a good NP dispersion in nanocomposites.⁴³,⁷⁵⁻⁷⁸ For a more detailed discussion of
$g_{NN}(r)$, we refer to the Supplementary Material, Sec. S2.3.

<figure><img src="image_4.png" /><figcaption>Figure 3: (a) Monomer-monomer structure factor for $\rho_{m0} = 0.290$ and different values of $\sigma_N$. (b) NP-NP structure factor as a function of the rescaled wavevector $q\sigma_N/2\pi$, for $\rho_{m0} = 0.290$ and different values of $\sigma_N$. The curves for RNPs are shifted up by 0.2 for clarity.</figcaption></figure>

Although the volume fraction of NPs added to the network is less than 2% in all the systems we considered, the presence of the NPs can still modify the structure of the network. In order to quantify this effect, we study two quantities: The monomer-monomer structure factor $S_{mm}(q)^{79}$ and the pore size distribution $P(r)^{80-82}$

The structure factor is defined as$^{79}$

$$S(\mathbf{q}) = \frac{1}{N} \sum_{k,j=1}^{N} \langle \exp[-i\mathbf{q} \cdot (\mathbf{r}_k - \mathbf{r}_j)] \rangle, \quad (8)$$

where $\mathbf{q}$ is the wavevector. Since our configurations are isotropic, we will consider the spherically averaged structure factor $S(q)$, with $q = |\mathbf{q}|$.

In Fig. 3a we show $S_{mm}(q)$ for $\rho_{m0} = 0.290$, for both ANPs and RNPs. The addition of ANPs leaves $S_{mm}(q)$ basically unchanged with respect to the neat system (thick black line*), except for $\sigma_N = 1$, for which we observe a shoulder at around $q \simeq 3$, i.e., $r \simeq 2\pi/q \simeq 2$, which results from configurations with two monomers separated by a NP. We note that for

<footnote>*The structure factor for the neat system has been computed for a system of $5 \times 10^4$ instead of $4 \times 10^5$ particles.</footnote>
small $q$, the structure factor shows a power-law, $S_{mm}(q) \propto q^{-\alpha}$, indicating that the network has a fractal nature.⁸³ For the case of RNPs, the presence of the NPs influences the structure factor much more strongly in that $S_{mm}(q)$ increases noticeably for $q \lesssim 1$, i.e., $r \gtrsim 6$. This is a consequence of the fact that the addition of RNPs causes a swelling of the network on intermediate length scales, as we will also discuss below when analyzing the pore size distribution (see also Supplementary Material, Sec. S2.2). This interpretation is also in agreement with the fact that the modification of $S_{mm}(q)$ is independent of NP size. Note that we do not consider RNPs with $\sigma_N = 1$, since they induce an excessive swelling of the network, as discussed in the Supplementary Material, Sec. S2.2. For additional details on $S_{mm}(q)$, see also Sec. S2.4 in the Supplementary Material.

In Fig. 3b we show the NP-NP structure factor $S_{NN}(q)$ for $\rho_{m0} = 0.290$ as a function of the rescaled wavevector $q\sigma_N/2\pi$. (The corresponding data for the other densities are shown in the Supplementary Material Sec. S2.5.) Studying this quantity is useful to determine whether or not the NPs are distributed homogeneously in the system. For the ANPs, with the exception of $\sigma_N = 1$, the data fall on a master curve. The curve is basically flat for $q\sigma_N/2\pi > 1/2$, $S_{NN}(q) \simeq 1$, like that of a gas, signaling that the NPs are homogeneously dispersed and aggregation is essentially absent. The only exception is $\sigma_N = 1$, for which a peak of modest height appears at $q \to 0$, signaling a weak NP aggregation in which several NPs (occasionally) fill the holes in the mesh. For the RNPs, on the other hand, the presence of weak aggregation is clear for all the systems, in that all the curves fall on a master curve with a main peak at $q \to 0$ and a smaller one at $q\sigma_N/2\pi = 1$, which results from direct contacts between NPs. The main peak at small $q$ is likely due to the fact that these NPs fill the largest cavities of the network. However, for $\sigma_N=5$ we find that the main peak is no longer at the smallest accessible wave-vector, likely because there are not sufficiently large cavities to accommodate a substantial number of large NPs. These trends are consistent with what we have concluded from the radial distribution function $g_{NN}(r)$ (Fig. 2b). For additional details on $S_{NN}(q)$, see also Sec. S2.5 in the Supplementary Material.
## 3.2 Pore size distribution

<figure><img src="image_5.png" /><figcaption>Figure 4: Distribution of pore diameters, $\hat{P}(2r)$, for $\rho_{m0} = 0.290$ and for different NP diameters. Dashed black line: $\hat{P}_0(2r)$ (neat system distribution). (a): ANP (b): RNP.</figcaption></figure>

A further quantity that is useful for the characterization of the network structure is the pore size distribution (PSD) $P(r)$.⁸² The PSD is calculated by randomly sampling points in the void space and finding the radius $r$ of the largest sphere that can be inserted without touching a polymer strand.⁸⁰⁻⁸² In Fig. 4 we show $P(r)$ for $\rho_{m0} = 0.290$, while the data for the other densities are presented in the Supplementary Material, Sec. S2.6. (To facilitate the comparison with the NP diameter $\sigma_N$, we actually show the distribution of the pore diameters, $\hat{P}(2r)$, a distribution that is trivially related to $P(r)$.) For the sake of comparison we have also included in the graph the PSD of the neat system, $P_0(r)$ (dashed black line). For the ANPs (Fig. 4a) and $\sigma_N > 1$, $P(r)$ is basically identical to $P_0(r)$, with the exception of a sharp peak appearing at $2r = \sigma_N$, which results from the NPs creating “cavities” in the network. For $\sigma_N = 1$, $P(r)$ is slightly shifted to smaller $r$ with respect to $P_0(r)$, with the main peak shifting to smaller $r$: The network contracts as a consequence of the attractive monomer-NP interaction, as already observed in non-crosslinked nanocomposites.⁶⁷ For the RNPs (Fig. 4b), $P(r)$ shifts to the right at small $r$, rationalizing the swelling of the network, and a broad peak appears centered at $2r \approx \sigma_N + \delta$ i.e., shifted with respect to the peak observed for ANPs. For $\rho_{m0} = 0.290$, one finds $\delta \approx 0.5$, but this value decreases with increasing network density (see Sec. S2.6 in the Supplementary Material). This shifting of
the peak is related to the fact that the RNPs have a larger effective diameter than the ANPs,
in agreement what discussed for $g_{mN}(r)$ (Fig. 2a). The peak is significantly broader than
the one seen for the ANPs for the same value of $\sigma_N$: This is likely due to the fact that for
ANPs, the network strands are attracted to the surface of the NPs, giving thus rise to a
well defined distance, whereas the RNPs repel the strands and locally deform the network.
Since the amplitude of this deformation will depend on the local properties of the mesh, the
resulting distribution of the pore size is broad.

The position of the main peak of $P(r)$, $r_{\max}$ (most probable pore radius), can be taken
as a measure of the mean mesh size $\xi$ of the system, i.e., $\xi \equiv r_{\max}$. In a previous study,⁸² it
was found that for purely polymeric systems $r_{\max}$ as well as the mean $\langle r \rangle$ can be taken as
reliable estimates of $\xi$. However, for the networks studied in the present work, $r_{\max}$ provides
a more meaningful value for the mean mesh size, since $P(r)$ can have a double peak structure
because of the presence of the NP. We note that $\xi$ has a weak dependence on the NP diameter
$\sigma_N$ and on the character of the monomer-NP interaction. However, this dependence can be
neglected for all but the smallest value of $\sigma_N$, and hence $\xi \approx \xi_0$, where $\xi_0$ is the mean mesh
size of the neat system.

# 4 Dynamic properties

## 4.1 Theoretical background

The dynamics of NPs in permanently crosslinked networks has been studied among oth-
ers by Dell and Schweizer¹ using the nonlinear Langevin equation (NLE) theory,⁵¹ and
independently by Cai, Panyukov and Rubinstein² using scaling theory. Dell and Schweizer
considered a polymer network characterized by an effective tube diameter $d$ resulting from
both the crosslinks and the entanglements, with $d \approx b(N_e^{\text{eff}})^{1/2} \dagger$ where $b$ is the Kuhn length

<footnote>†Here and in the following, we will use $\approx$ to signify equivalence apart from a dimensionless constant of order 1, and $\simeq$ to signify a numerical approximation.</footnote>
and $N_e^{\text{eff}}$ is the effective entanglement length $^{\ddagger}$. Making some simplifying assumptions regarding the structure of the polymeric matrix,¹ and introducing the confinement parameter $C \equiv \sigma_N/d$, one finds that the onset of localization happens at $C_c$ slightly larger than 1, as intuition suggests (the exact value of $C_c$ depends on the density and compressibility of the matrix). For $C$ slightly above $C_c$, the mean hopping time, i.e., the mean time between two hopping events, increases approximately exponentially with $C$, whereas for stronger confinements this increase becomes stronger than exponential, making that for large $C$ the NPs can be considered as indefinitely trapped on experimental time scales. The authors concluded therefore that hopping can be observed experimentally only in the weak confinement regime.

Cai et al.² considered a generic network with both crosslinks and entanglements, containing non-sticky NPs. Depending on the ratio between the tube diameter resulting from topological entanglements, $d_e \approx bN_e^{1/2}$ ($N_e$ = entanglement length), and the one resulting from crosslinks, $d_x \approx bN_x^{1/2}$ ($d_x$ being equivalent to the geometrical mesh size $\xi$, $N_x$ = mean strand length), they identified two dynamical regimes: the entangled regime ($N_e < N_x$, or equivalently $d_e < d_x$), where entanglements dominate, and the unentangled regime ($N_e > N_x$), where crosslinks dominate. Since all our systems are unentangled, see Sec. 4.2, we will only consider the latter regime for which one has $d \simeq d_x$. In such unentangled networks, NPs of diameter $\sigma_N < d$ can freely diffuse through the mesh, whereas NPs of diameter $\sigma_N > d$ can only diffuse via hopping processes. The resulting diffusion coefficient is predicted to have the following $C$ dependence:

$$D_N \propto \frac{b^2 \exp(-C^2)}{\tau_x C}, \qquad (9)$$

where $\tau_x \propto N_x^2$ is the Rouse relaxation time of a network strand and $C \equiv \sigma_N/d$.

<footnote>$^{\ddagger}$The effective entanglement length, $N_e^{\text{eff}}$, can be obtained experimentally from the plateau modulus $G_0$, since $G_0 \approx \rho k_B T / N_e^{\text{eff}}$</footnote>
## 4.2 Mean-squared displacement and diffusion coefficient

<figure><img src="image_6.png" /><figcaption>Figure 5: MSD, $\langle r_N^2(t) \rangle$, of the ANPs (a) and the RNPs (b) and of the crosslinks, $\langle r_{cl}^2(t) \rangle$, for $\rho_{m0} = 0.375$ and for different NP diameters. Dashed curves represent systems which have not reached the diffusive regime (see discussion in the text). Dashed horizontal line: $(2\xi_0)^2$ (see Tab. 1) Labels give the values of the confinement parameter C, Eq. (12).</figcaption></figure>

In order to test the validity of the theoretical predictions summarized in Sec. 4.1, we have determined the NP diffusion coefficient from the mean-squared displacement (MSD). In Fig. 5, we show the time dependence of the MSD,

$$ \langle r_N^2(t) \rangle \equiv \langle |\mathbf{r}_N(t) - \mathbf{r}_N(0)|^2 \rangle, \qquad (10) $$

for ANPs (a) and RNPs (b), for the densest network ($\rho_{m0} = 0.375$), and compare it with the MSD of the crosslinks, $\langle r_{cl}^2(t) \rangle$. (The corresponding curves for the other densities are in the Supplementary Material, Sec. S3.2.) Dashed lines denote systems for which the NPs have not reached the diffusive regime within the duration of the simulation, i.e., $\langle r_N^2(t) \rangle = 6D_N t$, where $D_N$ is the diffusion coefficient of the NPs. We stress that the MSD of the crosslinks is shown for all values of $\sigma_N$ but since the curves superimpose very well, only one curve is visible. We can therefore conclude that this quantity is basically independent of $\sigma_N$, confirming that the NP volume fraction we consider is sufficiently low to represent the dilute-NP limit.⁴³ Due to the fixed network topology, the MSD of the crosslinks quickly reaches a plateau,
<table><caption>Table 1: Properties of the neat polymer networks. $\rho_{m0}$: monomer density. $\lambda$: crosslink localization length, Eq. (11). $\xi_0$: mean mesh size in the neat state ($2\xi_0$ = most probable pore diameter). $d$: Estimate for the effective tube diameter, $d \approx \lambda(N_e/N_x)^{1/4} = 1.94 \times \lambda$ ($N_e = 85, N_x = 6$).</caption><thead><tr><th>ρ<sub>m0</sub></th><th>λ</th><th>2ξ<sub>0</sub></th><th>d</th></tr></thead><tbody><tr><td>0.190</td><td>3.12</td><td>3.06</td><td>6.05</td></tr><tr><td>0.290</td><td>2.01</td><td>2.06</td><td>3.90</td></tr><tr><td>0.375</td><td>1.51</td><td>1.58</td><td>2.93</td></tr></tbody></table>

$\langle r_{\text{cl}}^2(t) \rangle = \lambda^2 = \text{const.}$, the height of which defines the (squared) localization length,⁸⁴ i.e. the mean amplitude of the fluctuations of the crosslinks around their equilibrium position:

$$ \lambda \equiv \left[ \lim_{t \to \infty} \langle r_{\text{cl}}^2(t) \rangle \right]^{1/2}. \qquad (11) $$

For $\rho_{m0} = 0.375$, we find $\lambda = 1.51$, a value that is close to $2\xi_0$ (see horizontal dashed line in the graph), the most probable pore diameter $r_{\text{max}}$ in the neat state (see Fig. 4). Since for the other values of $\rho_{m0}$ we find the same relation, see Tab. 1, we can conclude that the amplitude of the oscillations of the crosslinks around their equilibrium positions are basically equal to the mean mesh size of the network, a result that is certainly very reasonable. We also note that $\lambda$ is proportional to the effective tube diameter $d$ of the system, since $d \approx \lambda(N_e/N_x)^{1/4}$.⁴³,⁸⁵ Assuming $N_e = 85$⁸⁶ and since $N_x = 6$, we obtain $d \approx 1.94 \times \lambda$. We also remark in passing that $N_e \gg N_x$, i.e., all our systems are in the unentangled regime.

Since previous studies have shown that it is useful to discuss the dynamics as a function of the so-called confinement parameter¹,²,⁵² we define

$$ C \equiv \frac{\sigma_N}{\lambda}. \qquad (12) $$
Note that this definition is slightly different from the one of previous studies in that we use
here $\lambda$ as the length scale for the normalization instead of the tube diameter $d$. This choice
is motivated by the fact that $\lambda$ can be measured with high precision from the simulations
while $d$ is a length scale that rests on a theoretical concept that allows to define $d$ only up
to a numerical prefactors. As discussed above, however, for all simulated systems we have
$\lambda \approx d/1.94$, a relation which allows to express our results also in terms of $d$. We also note
that $\lambda \simeq 2\xi_0$, where $2\xi_0$ is the mesh size of the neat system (Tab. 1).

The values of *C* are reported in the labels of the curves in Figs. 5a and b. One can see
that for *C* ~ 1, the MSD becomes diffusive, i.e., ⟨r<sup>2</sup><sub>N</sub>(t)⟩ ~ t, immediately after the ballistic
regime, irrespective of the type of NP. More interesting is the dynamics of the NPs with
*C* ~ 1, in that their MSD shows between the ballistic regime and the diffusive behavior
at long times, a subdiffusive regime, ⟨r<sup>2</sup><sub>N</sub>(t)⟩ ~ t<sup>β</sup>, with an exponent *β* < 1 that depends
on *C*. This subdiffusive regime begins when ⟨r<sup>2</sup><sub>N</sub>(t)⟩ ~ d<sup>2</sup>σ/σ<sub>N</sub> ~ 3.76 × λ<sup>2</sup>σ/σ<sub>N</sub><sup>1.2</sup> (see
Supplementary Material, Sec. S3.1 and discussion below). Note that for the largest value
of *C*, the systems do not reach the diffusive regime neither during equilibration nor during
production, since the simulation time required would be prohibitively long (dashed lines).
Although strictly speaking these systems have not fully equilibrated, we have not observed
any noticeable time dependence in any of the structural and thermodynamic quantities (i.e.,
no aging is observed), and have therefore decided to include them in the analysis as well.
Moreover, as discussed in Sec. 4.1, for strong confinements we expect the typical hopping time
of the NPs to diverge at least exponentially in *C*. This implies that, even in experiments, a
slight increase of *C* will prevent the system from reaching equilibrium. Therefore, we expect
that our data for the extreme confinement regime is still useful for the comparison with
experimental data.

The MSD from the simulations can be compared with the theoretical predictions by Cai
et al.,² who for the case of large particles in an unentangled network predict four regimes in
the dynamics: After the initial ballistic regime, one finds a subdiffusive regime ⟨r<sup>2</sup><sub>N</sub>(t)⟩ ∝ t<sup>β</sup>
with $\beta = 1/2$, which lasts up to the Rouse relaxation time of a network strand, $\tau_x \approx \tau_0 N_x^2$, where $\tau_0 \approx m\gamma_m\sigma^2/k_B T$ is the monomer relaxation time $^§$. For $t > \tau_x$, the NP is trapped by the mesh, and the MSD shows a plateau $\langle r_N^2(t) \rangle \approx d^2\sigma/\sigma_N$. For times larger than the hopping time, i.e. $t > \tau_x N_x^{1/2} \exp(-\sigma_N^2/d^2)$, the NP finally escapes from the local trap formed by the mesh and starts to diffuse via hopping motion. Comparing these predictions with our simulation data shows that we do not observe the intermediate subdiffusive $t^{1/2}$ regime, which is likely due to the fact that the average length of the network strands, $N_x = 6$, is too small to be in the scaling regime. For additional details on the MSD of the NPs, we refer to Sec. S3.2 of the Supplementary Material.

<figure><img src="image_7.png" /><figcaption>Figure 6: Apparent subdiffusive exponent, Eq. (13), for systems of RNPs (a) and ANPs (b) for different $\rho_{m0}$ and similar values of C.</figcaption></figure>

In order to study the dependence of the subdiffusive dynamics on $\sigma_N$, we consider the apparent subdiffusive exponent $\beta(t)$, defined as

$$ \beta(t) \equiv \frac{d \ln[\langle r_N^2(t) \rangle]}{d \ln(t)}. \qquad (13) $$

In Fig. 6, we show $\beta(t)$ for systems with different $\rho_{m0}$ and similar values of C, for the RNPs (a) and the ANPs (b). The short-time ballistic regime ($\beta(t) = 2$) is followed by a sharp decrease of $\beta(t)$ when the MSD of the NPs leaves the ballistic regime. For $C \le 1$ the exponent

<footnote>$^§$We recall that $m$ is the monomer mass, $\sigma$ the diameter, and $\gamma_m$ the friction coefficient.</footnote>
becomes unity, thus indicating the diffusive motion. If the confinement parameter is larger
than 1.0 the exponent $\beta$ drops to a value below unity and subsequently increases slowly with
time. (Note that $\beta(t)$ is independent of the density $\rho_{m0}$ but decreases if $C$ is increased, which
is further evidence that $C$ is the relevant parameter for the dynamics.) In the subdiffusive
regime the time dependence of $\beta$ is compatible with a logarithmic dependence, a result that
indicates that in this time window the dynamics of the NP is very heterogeneous. This result
can be rationalized by the fact that the time at which a given particle starts to leave the cage
formed by the mesh depends on the details of this cage (number of strands, their length,...)
giving rise to a very broad distribution of local relaxation times. Note that if this type of
disorder is absent it can be expected that the escape times from this cage is not broadly
distributed, making that for a long time the NP will stay confined and the MSD shows
a plateau. Hence we conclude that the plateau that is predicted in Cai et al.² is washed
out due to the intrinsic disorder of the gel. For more details on the apparent subdiffusive
exponent, we refer to the Supplementary Material, Sec. S3.1.

<figure><img src="image_8.png" /><figcaption>Figure 7: NP diffusion coefficient as a function of the confinement parameter $C$, Eq. (12). Unless specified, $\phi_{N0} = 0.02$. Solid lines: ANPs. Dashed lines: RNPs. Open symbols: Upper bounds for $D_N$.</figcaption></figure>

As discussed above, the NPs are eventually able to escape from the mesh and start to
diffuse. The dynamics in the diffusive regime can be characterized by the diffusion coefficient
of the NPs, which can be obtained from the MSD via the Einstein relation:⁷⁹
$$D_N = \lim_{t \to \infty} \frac{1}{6t} \langle r_N^2(t) \rangle. \qquad (14)$$

In Fig. 7, we show the NP diffusion coefficient $D_N$ as a function of the confinement parameter $C$ for all the systems studied here. For the systems which have not reached the diffusive regime, we report an upper bound for $D_N$, represented by open symbols. One sees that for the ANPs as well as the RNPs the data falls nicely on a master curve when plotted as a function of $C$, which indicates that the confinement parameter is the relevant quantity determining also the long-time diffusive dynamics. We also note that a relatively small increase of $C$ causes a dramatic decrease of $D_N$. For example, for the ANPs in the $\rho_{m0} = 0.375$ system, $D_N$ decreases from $2.23 \times 10^{-1}$ to less than $2.95 \times 10^{-6}$ (upper bound) when going from $C = 0.664$ to $C = 3.98$, i.e., a reduction of five orders of magnitude.

It is interesting to observe that for $C < 1$ (weak confinement), the diffusivity of ANPs is smaller than that of RNPs, while for $C > 1$ the reverse is true. This can be rationalized as follows: In the weak confinement regime the NPs can freely diffuse through the mesh, and do not have to wait for the relaxation of the local mesh. Thus in this regime, an attractive monomer-NP interaction will causes an increase of the effective friction felt by the NP, and hence the $D_N$ is smaller for ANPs. For $C > 1$ (strong and extreme confinement), these short-range effects lose their importance since they are dominated by the confinement effects induced by the mesh on the NP. In this regime it is thus the effective diameter of the NP that is the relevant quantity and, as we will discuss below, this diameter is smaller for the ANPs than the one for the RNPs, thus explaining why at large $C$ the diffusion of the RNPs is slower than the one of the ANPs.

In order to verify that the NP volume fraction ($\phi_N \simeq \phi_{N0} = 0.02$) is sufficiently low to be in the single-NP (dilute) limit, we also include in Fig. 7 data for $\rho_{m0} = 0.190$ and NP volume fraction $\phi_N \simeq \phi_{N0} = 0.005$. One can see that the difference between the higher and
lower volume fraction is negligible, confirming the validity of our choice.

<figure><img src="image_9.png" /><figcaption>Figure 8: Same data as in Fig. 7, rescaled according to Eq. (15) (a) and Eq. (16) (b). $C_{eff} \equiv BC$ (a) and $B'C$ (b) is the effective confinement parameter. (a): RNPs: $A = 0.265$, $B = 0.928$. ANPs: $A = 0.206$, $B = 0.809$. (b): RNPs: $A' = 0.970$, $B' = 1.56$, $\delta = 1.49$. ANPs: $A' = 0.989$, $B' = 1.35$, $\delta = 1.53$. Dashed lines: exponential with slope $-1$.</figcaption></figure>

As discussed above, all our systems are in the unentangled regime, i.e. $N_x < N_e$: According to the theory of Cai et al.² we thus expect $D_N$ to be given by (see Eq. (9)):

$$D_N = A \frac{\exp[-(BC)^2]}{BC} \equiv A \frac{\exp(-C_{eff}^2)}{C_{eff}}. \qquad (15)$$

where $A, B$ are positive constants and $C_{eff} \equiv BC$ is the effective confinement parameter. Fitting the data in Fig. 7 with the expression (15), we find for the RNPs $A = 0.265$ and $B = 0.928$, so that $C_{eff} \simeq C$, confirming that for the RNP $\sigma_N/\lambda$ represents indeed the relevant confinement parameter (Eqs. (9) and (12)). For the ANPs, we restrict the fit to $C > 2$, in order to probe the intermediate-strong confinement regime, for which hopping is expected to take place. We find $A = 0.206$ and $B = 0.809$, so that once again $C_{eff} \simeq C$. We note, however, that the ANPs behave as if they experience an effective confinement parameter $C_{eff} = BC = B\sigma_N/\lambda$, i.e., their effective size is smaller, as discussed for $g_{Nm}(r)$ (Sec. 3.1). Note that we find for both types of particles a value of $B$ which is close to unity, indicating that $C = \sigma_N/\lambda$ is the relevant parameter for a quantitative prediction of $D_N$.
In Fig. 8a, we show the same data as in Fig. 7, rescaled according to Eq. (15) and using the fit parameters mentioned above for the RNPs and the ANPs. One can see that the data fall very well on a master curve, indicating that the expression by Cai et al. for the diffusion constant does give a good description of the data. The only exceptions to this are the ANPs at small values of C in that one can note some small but systematic deviation from the theoretical curve. This disagreement is likely due to the fact that for the ANPs the mechanism of diffusion at small C is dominated by the small-scale friction felt by the NPs, as discussed above.

Finally, we compare our data with the theoretical prediction by Dell and Schweizer.¹ For $C \gtrsim C_c$, diffusion proceeds via hopping (see Sec. 4.3 below), and the hopping diffusion coefficient can be estimated as $D_N \approx \Delta_h^2/\tau_h$, where $\Delta_h^2$ is the mean jump length and $\tau_h$ the mean hopping time.¹ Since the dependence of $\Delta_h$ on C is weak for small values of $C - C_c \simeq C - 1$, to a first approximation one has $D_N \propto \tau_h^{-1} \propto \exp(-F_B/k_B T)$, where $F_B$ is the free-energy barrier associated with the hopping process. This free energy barrier is found to increase with C as a power law for small values of $C - C_c$, so that $D_N$ is predicted to be given by

$$D_N = A' \exp[-(B'C)^\delta] = A' \exp[-C_{\text{eff}}^\delta], \quad (16)$$

with $A'$, $B'$, and $\delta$ positive constants. Fitting the data in Fig. 7 with Eq. (16), we find for the RNPs $A' = 0.970$, $B' = 1.56$, and $\delta = 1.49$, whereas for the ANPs (restricting also in this case the fit to $C > 2$, see discussion above) we find $A' = 0.989$, $B' = 1.35$, and $\delta = 1.53$. Thus the value of the exponent $\delta$ we find is independent of the type of particles considered and close to the one from the theory, i.e., $1.7 \pm 0.1$, valid for value of C that are not too large.¹ Furthermore we note that when using Eq. (16), one finds (again) that the effective diameter of the ANPs is slightly smaller than that of the RNPs.
In Fig. 8b, we report the data from Fig. 7, but now rescaled according to Eq. (16). We observe an excellent agreement with the theoretical prediction, in that the data collapse on the same master curve. Also in this case, there are small deviations from the master curve at small *C* for the ANPs, for the reasons discussed above. We conclude this analysis by observing that, since both Eq. (15) and Eq. (16) are in very good agreement with our data, we cannot, at present, conclude that one of the two theoretical descriptions is significantly better than the other. In order to draw any further conclusion, longer simulations are needed, so that more precise estimates of *D*<sub>N</sub> in the *C* &gt; 2.5 regime (strong/extreme confinement) can be obtained.

Before concluding this discussion, we also mention that the confinement parameter is not only able to describe the dynamics of the relaxation dynamics of the NP but also, e.g., the amount of shrinking/expansion of the network due to the NP (see Sec. S2.2 of the Supplementary Material for more details). We thus conclude that *C* is not only a relevant parameter for the dynamics, but also for certain static properties of the NP-gel system.

## 4.3 Analysis of the hopping dynamics: van Hove function

<figure><img src="image_10.png" /><figcaption>Figure 9: (a): Squared displacement of three RNPs denoted A,B, and C (solid lines) and MSD of the crosslinks (dashed line) for $\rho_{m0} = 0.190$ and $\sigma_N = 8.0$ ($C = 2.56$). (b): Trajectories of the NPs A,B and C, created by taking snapshots of the NP position at regular time intervals of length $\Delta t = 6 \times 10^3$ from the time $t = 0$ to $t = 6 \times 10^5$. To facilitate visualization, the NPs are represented as having half their real diameter.</figcaption></figure>
To provide a better understanding of the hopping dynamics, we compare in Fig. 9a the squared displacement (SD) $r_N^2(t)$ of three RNPs (labeled A, B and C) with the MSD of the crosslinks. One can see that particle C is trapped by the mesh, in that it displays, after the brief ballistic regime, a completely flat SD, with a localization length approximately equal to the localization length of the crosslinks, $r_{cl}^2(t) \simeq \lambda$. Also the A and B particles are initially trapped by the mesh, but eventually manage to escape. The SD of the NPs consists in a series of abrupt jumps, each of which is followed by a flat region. This type of interrupted progression is one of the typical signatures of a hopping process. In order to clarify the nature of the motion, we show in Fig. 9b snapshots of the particles taken at regular time intervals of duration $\Delta t = 6 \times 10^3$, from the time $t=0$ to $t=6 \times 10^5$. These 100 snapshots clearly reveal the existence of hopping motion, in that the positions of the mobile particles A and B at different times form well-separated clusters. The same qualitative behavior is observed, under similar conditions, for the ANPs.

<figure><img src="image_11.png" /><figcaption>Figure 10: Self part of the van Hove function, Eq. (17), for $\rho_{m0} = 0.375$, $\sigma_N = 3.5$ for the ANPs (a) and the RNPs (b). $\tau_D \equiv \sigma_N^2/6D_N$ is the NP diffusion time. Dashed lines: Gaussian approximation, Eq. (18).</figcaption></figure>

In order to characterize the hopping dynamics in a quantitative manner, we probe the self part of the van Hove function, which is defined as:⁷⁴
$$G_{s,N}(\mathbf{r}, t) = \langle \delta[\mathbf{r} - (\mathbf{r}_{N,i}(t) - \mathbf{r}_{N,i}(0))] \rangle . \qquad (17)$$

Since our configurations are isotropic, we consider the spherically-averaged van Hove function, $G_{s,N}(r, t)$. From $G_{s,N}$ one can define the distribution of the displacements, $\Delta_N(r, t)dr = 4\pi r^2 G_{s,N}(r, t)dr$, which represents the probability to find a particle at time $t$ at distance between $r$ and $r+dr$ from its original position. At short times the motion is ballistic and hence $G_{s,N}(r, t)$ is a Gaussian$^{79}$ with variance $2\langle r_N^2(t) \rangle/3$:

$$G_{s,N}(r,t) = \left( \frac{3}{2\pi \langle r_N^2(t) \rangle} \right)^{3/2} \exp \left( -\frac{3r^2}{2\langle r_N^2(t) \rangle} \right). \qquad (18)$$

At long times the NP motion becomes uncorrelated and Gaussian behavior is recovered again.$^{79}$ In Fig. 10, we show $\Delta_N(r, t)$ for $\rho_{m0} = 0.375$, $\sigma_N = 3.5$ for the ANPs (a) and the RNPs (b). Let us introduce the diffusion time,

$$\tau_D \equiv \frac{\sigma_N^2}{6D_N}, \qquad (19)$$

which represents the time it takes a NP with diffusion coefficient $D_N$ to move over a distance of the order of its diameter. At short times ($t \ll \tau_D$), both for ANPs and RNPs, $G_{s,N}(r, t)$ is well described by Eq. (18) (dashed curves in Fig. 10), and $\Delta_N(r, t)$ displays a single peak at $r = [2\langle r_N^2(t) \rangle/3]^{1/2}$, which progressively moves towards larger $r$ values. When the position of the peak approaches the crosslink localization length $\lambda$, the NP interacts strongly with the mesh and the dynamics is sub-diffusive, see Fig. 5. For longer times $\Delta_N(r, t)$ maintains the peak at $r \simeq \lambda$, while at the same time developing a weak shoulder at $r \simeq 2\lambda$. This double peaked distribution is the clear signal of heterogeneous hopping dynamics:$^{87}$ While
some particles are still trapped in their initial cage (primary peak at $r \simeq \lambda$), other particles have escaped from the cage and started a diffusive motion. The presence of the second peak indicates that the escaping NPs perform “jumps” with a typical length comparable to the mesh size, as predicted by Cai et al.² In simulations of NP-charged networks with a regular structure, multiple peaks can be observed at r-values which are multiples of the mesh size.⁴⁴ In contrast to this, we do not find in our system such multiple peaks, likely because our network is disordered. For $t > \tau_D$, the primary peak eventually disappears completely as all the particles have escaped from their initial cage.

<figure><img src="image_12.png" /><figcaption>Figure 11: (a): Non-Gaussian parameter for the ANPs with different values of $\sigma_N$ at $\rho_{m0} = 0.375$, shown as a function of time (a) and of the rescaled time $t/\tau_D$, with $\tau_D \equiv \sigma_N^2/6D_N$ (b). Labels in (a) give the value of the confinement parameter C. Dashed curves represent systems which have not reached the diffusive regime.</figcaption></figure>

The heterogeneity of the relaxation dynamics can be characterized by probing the deviations of $G_s(r, t)$ from a Gaussian, which can be quantified using the non-Gaussian parameter,⁸⁷

$$ \alpha_2(t) = \frac{3\langle r^4(t) \rangle}{5\langle r^2(t) \rangle^2} - 1. \qquad (20) $$

Thus if $\alpha_2 = 0$, the dynamics is Gaussian, whereas a value which is appreciably different from 0 signals strong heterogeneities in the dynamics. In Fig. 11 we report $\alpha_2(t)$ for $\rho_{m0} = 0.375$
for the ANPs (dashed lines denote systems where the diffusive regime was not reached). At short times $\alpha_2$ is very small, i.e., the dynamics of the NPs is Gaussian because of the ballistic motion. This regime is followed by an increase of $\alpha_2$, which corresponds to the onset of the subdiffusive regime in the MSD. For $C \lesssim 2$, this increase is very modest, and $\alpha_2$ remains smaller than 1. This result is in agreement with the fact that for small C the MSD crosses over from the ballistic regime directly to the diffusive regime, see Fig. 5. However, for strong confinement, $C \gtrsim 3$, $\alpha_2$ reaches values significantly larger than 1. These values are comparable to those obtained from experimental studies of polymer networks containing NPs$^{52}$ or in simulations of deeply supercooled liquids,$^{87}$ i.e., systems which show pronounced dynamical heterogeneities. In our system, as in experimental ones,$^{52}$ this pronounced non-Gaussian dynamics is directly related to the structural heterogeneity of the network on the length scale of the NP size: Some NPs diffuse freely through the gaps in the network, whereas other NPs end up trapped in the mesh for a long time. The time at which $\alpha_2(t)$ reaches its maximum also increases strongly with C, i.e., the time needed to make one hopping movement depends as expected strongly on the confinement. In Fig. 11b we plot $\alpha_2$ as a function of the rescaled time $t/\tau_D$, where $\tau_D$ is defined in Eq. (19). We observe that for $C \lesssim 3$ the descending parts of all the curves fall on the same master curve, confirming that it is indeed the time scale $\tau_D$ which controls the long-time dynamics. We also note that for small C-values, $\alpha_2(t)$ reaches its maximum for $t \simeq \tau_D$, as also observed in polymer solutions containing NPs.$^{67}$ For extreme confinement, one cannot conclude whether or not the Gaussian dynamics is recovered at long times, since we are unable to reach the diffusive regime. From this graph one also recognizes that the width of the peak in $\alpha_2(t)$ increases quickly with increasing C, a further indication that the dynamical heterogeneity of the system is increasing. We also mention that the time dependence of $\alpha_2$ for the RNPs display the same qualitative behavior as the one for the ANPs, and the same is true when different values of $\rho_{m0}$ are considered (see Supplementary Material, Sec. S3.3).

To probe the collective dynamics of the NPs we investigate the distinct part of the van
<figure><img src="image_13.png" /><figcaption>Figure 12: Distinct part of the van Hove function $G_{d,N}(r,t)$ for $\rho_{m0} = 0.375$, $\sigma_N = 2.0$ (a, b) and $\rho_{m0} = 0.375$, $\sigma_N = 4.0$ (c, d). (a, c) show the RNP and (b, d) show the ANP. $\tau_D \equiv \sigma_N^2/6D_N$ is the NP diffusion time. Inset of (a): $r|G_d(r,t) - 1|$ on a logarithmic scale.</figcaption></figure>
Hove function which is given by<sup>79</sup>

$$G_{d,N}(\mathbf{r}, t) = \frac{1}{N_N} \sum_{\substack{k=1 \\ j \neq k}}^{N_N} \langle \delta[\mathbf{r} - (\mathbf{r}_{N,j}(t) - \mathbf{r}_{N,k}(0))] \rangle. \quad (21)$$

For our isotropic system the quantity $4\pi r^2 \rho_N G_{d,N}(r,t)dr$ is thus proportional to the probability to find at time $t$ a particle at a distance between $r$ and $r+dr$ from the position occupied by another particle at time 0. For $t=0$ one recovers the NP-NP radial distribution function $g_{NN}(r)$. In Fig. 12a we show $G_{d,N}(r,t)/\rho_N$ for a system with RNP with $\rho_{m0} = 0.375$ and $\sigma_N = 2$ ($C = 1.99$). One observes that with increasing $t/\tau_D$ the correlation hole at $r < \sigma_N$ is slowly filled and the contact peak at $r/\sigma_N \simeq 1$ is washed out. If the NPs were able to completely explore the system's volume at $t \gg \tau_D$, we would expect $G_{d,N}(r,t)/\rho_N \simeq 1$, however this is not what we find, as shown in the inset of Fig. 12a. The slope on the right hand side of the contact peak becomes independent of time for all values of $t/\tau_D$, even though the NPs reach the diffusive regime (Fig. 5b) and the dynamics is basically Gaussian (see Supplementary Material). At the same time we find at long times a pronounced peak at $r \approx 0$, i.e., a probability that is enhanced with respect to the one for an ideal gas. This behavior is due to the permanent nature of the network: A given region of space that at $t=0$ is occupied by a NP is likely to have a NP also at infinitely long times, since the local network structure makes that this spot has a higher probability than average or, put otherwise, the frozen in disorder of the network makes that the stationary distribution of the NP is non-uniform in space. Qualitatively the same behavior is observed for the ANP, Fig. 12b. Although in this case there is no pronounced nearest neighbor peak, see discussion in the context of Fig. 2, we find a marked peak at small distances even at long times, i.e., the signature that certain regions in the network are highly preferential sites for the NP. This is thus evidence that the NP spend a substantial time in these sites and when they leave them they move quickly to another favorable site, i.e., that there is hopping dynamics.<sup>83,87</sup>
For stronger confinement, hopping motion becomes even more relevant, and the peak at $r = 0$ is dramatically more pronounced, as shown in Figs. 12c-d, where we report $G_d(r, t)/\rho_N$ for $\rho_{m0} = 0.375$ and $\sigma_N = 4$ ($C = 2.66$) for RNPs (c) and ANPs (d) (note the logarithmic scale for the ordinate). Although the statistics is worse than for the case of smaller values of $C$ and hence does not allow us to study accurately the long-range behavior of $G_d(r, t)$, we clearly observe the appearance of a peak at $r \simeq 0$. At long times, this peak grows significantly, signaling that NP motion is dominated by hopping dynamics.

# 5 Summary and conclusions

We have carried out molecular dynamics simulations of a disordered and polydisperse polymer network in which we have embedded nanoparticles at a fixed volume fraction $\phi_N \simeq 0.02$. Three different networks were studied, all with trivalent crosslinks, with densities in the neat state $\rho_{m0} = 0.190, 0.290$, and $0.370$. The size of the NP ranged between 1 and 10, allowing to probe weak as well as extreme confinement by the network mesh, and we considered NP-polymer interactions that were either attractive or repulsive.

The static structure factor indicates that the ANPs are well dispersed in the network, whereas RNPs show a weak tendency to cluster. From the analysis of the pore size distribution of the gel, $P(r)$, one can conclude that the NP locally deform the mesh and cause the appearance of a secondary peak in $P(r)$, the location of which is directly related to the size of the NPs.

The analysis the MSD of the NP shows that the confinement parameter $C = \sigma_N/\lambda$, where $\lambda$ is the localization length of the crosslinks, is the relevant parameter that determines the dynamics, in agreement with theoretical predictions.¹,² Three dynamic regimes can be identified: Weak confinement ($C \lesssim 1$), strong confinement ($1 \lesssim C \lesssim 2$), and extreme confinement ($C \gtrsim 3$). In the weak confinement regime, the NPs can freely diffuse through the mesh. In the strong confinement regime, the MSD of the NP displays a subdiffusive transient,
$\langle r_N^2(t) \rangle \propto t^\beta$ with $\beta < 1$, before eventually recovering diffusive behavior. In the extreme confinement regime, a marked subdiffusive regime appears in the MSD on intermediate time scales during which the NP undergo a very heterogeneous dynamics. This heterogeneity is directly related to the frozen in disorder of the network and is an ingredient which, to the best of our knowledge, has so far not been properly been taken into account in the theoretical approaches to describe these systems.

In the strong and extreme confinement regimes, NP motion can basically only proceed through the mechanism of activated hopping,¹,² *i.e.*, by waiting for a thermal fluctuation of the mesh that allows a NP to jump to a nearby cage. We find that the exponent $\beta(t)$ and the long-time NP diffusion coefficient, $D_N$, are controlled by $C$, but depend on the type of NP. The diffusion coefficient decreases dramatically with increasing $C$, displaying in one case a drop of five orders of magnitude upon an increase of $C$ by a factor 6. The $C$-dependence of $D_N$ can be described very well with the theoretical predictions by Cai *et al.*² but also with the ones by Dell and Schweizer¹ and hence we conclude that with the present set of data it is not possible to decide which theory is more reliable. More quantitative calculations using these two approaches should therefore be done in the future, including, if possible, the disordered nature of the network.

Finally we study the details of the hopping dynamics by analyzing the van Hove function of the NPs and the non-Gaussian parameter. The time and space dependence of the self part of the van Hove function confirms that, in the strongly confined regime, the NPs move indeed through activated hopping in that the function shows several peaks from which one can infer that the jump length is close to the mean mesh size, in agreement with the predictions of Cai *et al.*² and with a recent simulation study.⁴⁴ At the same time the motion is strongly non-Gaussian with values of the non-Gaussian parameter comparable to those found in deeply supercooled liquids.⁸⁷ The prevalence of hopping motion is confirmed by the analysis of the distinct part of the van Hove function

In conclusion, we have studied for the first time in simulations the diffusion of NPs
in disordered and polydisperse networks, probing dynamic regimes up to extreme confine-
ments of the NPs by the mesh. This work represents the natural extension of recent efforts
in simulations of crosslinked nanocomposites. 42-44 Understanding the dynamics of NPs in
crosslinked networks still poses a formidable challenge both for theoretical approaches and
for simulations, in the latter case due to the enormous relaxation times of these systems.
Although at present it is not possible to fully equilibrate such systems in the extreme con-
finement regime, many questions remain whose answer is within reach, and which will help
to understand better real-life systems. Examples are the interplay between crosslinks and
entanglements in polydisperse systems, the role of chemical *versus* physical crosslinks, the
effect of NP shape etc., questions that remain to be clarified. Some of these questions can
probably be addressed by equilibrating the systems using clever Monte Carlo algorithms,
such as parallel tempering. Moreover, a more detailed analysis is required in order to com-
pare qualitatively different theoretical predictions. 1,2 We are therefore convinced that the
study of such systems will remain a fruitful and challenging topic of research in the future.

Acknowledgement

We thank K. Schweizer and M. Lenz for useful discussions. The analysis of the pore size
distribution has been performed with the open-source *baggianalysis* software by L. Rovi-
gatti (DOI: 10.5281/zenodo.4588503). This work has been supported by LabEx NUMEV
(ANR-10-LABX-20) funded by the "Investissements d'Avenir" French Government program,
managed by the French National Research Agency (ANR).
# SUPPLEMENTARY MATERIAL

In this Supplementary Material we give more details on the properties of the systems, such as the swelling behavior, the structure factor, pore size distribution, and the mean squared displacement.

## S1 Properties of the simulated systems

Here we list some properties of the systems we simulated. Tables S1 and S2 are for the repulsive and attractive NP, respectively.

## S2 Structure

### S2.1 Chain length distribution

<figure><img src="image_14.png" /><figcaption>Figure S1: Rescaled chain length distribution $p(n)/p(1)$ for the simulated networks. Solid line: Theoretical prediction given by Eq. (1).</figcaption></figure>

In the main text we have mentioned that the distribution of the chain length, $p(n)$, is given by an exponential. In Fig. S1 we show $p(n)$, divided by $p(1) = (N_m - N_c)/(\langle n \rangle)^2$, for the three simulated networks. We observe that $p(n)$ is independent of $\rho_{m0}$ and decreases indeed exponentially,⁶⁴ i.e. the functional form given by the Flory-Stockmayer formula,⁶²,⁶³ Eq. (1) of the main text.
<table><caption>Table S1: Properties of the simulated systems (RNPs). $\sigma_N$: NP diameter. $N_N$: Number of NPs. $\rho_{m0} = M/V_0$: monomer number density in the neat system (independent of $\sigma_N$). $\phi_m = \pi\rho_N\sigma^3/6$: monomer volume fraction. $\phi_N = \pi\rho_N\sigma_N^3/6$: NP volume fraction. $C \equiv \sigma_N/\lambda$: Confinement parameter ($\lambda$ = crosslink localization length). $D_N$: NP diffusion coefficient. $\tau_D \equiv \sigma_N^2/6D_N$: Diffusion time. Values marked with * are upper and lower bounds, for $D_N$ and $\tau_D$ respectively, for systems which have not equilibrated.</caption><thead><tr><th>σ<sub>N</sub></th><th>N<sub>N</sub></th><th>ρ<sub>m0</sub></th><th>φ<sub>m</sub></th><th>φ<sub>N</sub></th><th>C</th><th>D<sub>N</sub></th><th>τ<sub>D</sub></th></tr></thead><tbody><tr><td>2.0</td><td>9755</td><td>0.1902</td><td>0.0823</td><td>0.0165</td><td>0.641</td><td>3.758 × 10<sup>-1</sup></td><td>1.774 × 10<sup>0</sup></td></tr><tr><td>3.0</td><td>2890</td><td>0.1902</td><td>0.0900</td><td>0.0181</td><td>0.962</td><td>1.331 × 10<sup>-1</sup></td><td>1.127 × 10<sup>1</sup></td></tr><tr><td>4.0</td><td>1219</td><td>0.1902</td><td>0.0928</td><td>0.0186</td><td>1.282</td><td>4.862 × 10<sup>-2</sup></td><td>5.484 × 10<sup>1</sup></td></tr><tr><td>5.0</td><td>624</td><td>0.1902</td><td>0.0943</td><td>0.0189</td><td>1.603</td><td>1.634 × 10<sup>-2</sup></td><td>2.550 × 10<sup>2</sup></td></tr><tr><td>6.0</td><td>361</td><td>0.1902</td><td>0.0951</td><td>0.0191</td><td>1.923</td><td>4.230 × 10<sup>-3</sup></td><td>1.419 × 10<sup>3</sup></td></tr><tr><td>7.0</td><td>228</td><td>0.1902</td><td>0.0956</td><td>0.0192</td><td>2.244</td><td>9.627 × 10<sup>-4</sup></td><td>8.483 × 10<sup>3</sup></td></tr><tr><td>8.0</td><td>152</td><td>0.1902</td><td>0.0958</td><td>0.0192</td><td>2.564</td><td>*2.480 × 10<sup>-4</sup></td><td>*4.301 × 10<sup>4</sup></td></tr><tr><td>2.0</td><td>6538</td><td>0.2898</td><td>0.1376</td><td>0.0181</td><td>0.997</td><td>1.444 × 10<sup>-1</sup></td><td>4.615 × 10<sup>0</sup></td></tr><tr><td>3.0</td><td>1937</td><td>0.2898</td><td>0.1428</td><td>0.0188</td><td>1.496</td><td>3.222 × 10<sup>-2</sup></td><td>4.655 × 10<sup>1</sup></td></tr><tr><td>3.5</td><td>1220</td><td>0.2898</td><td>0.1442</td><td>0.0190</td><td>1.745</td><td>1.343 × 10<sup>-2</sup></td><td>1.520 × 10<sup>2</sup></td></tr><tr><td>4.0</td><td>817</td><td>0.2898</td><td>0.1448</td><td>0.0191</td><td>1.994</td><td>4.995 × 10<sup>-3</sup></td><td>5.339 × 10<sup>2</sup></td></tr><tr><td>4.5</td><td>574</td><td>0.2898</td><td>0.1456</td><td>0.0192</td><td>2.244</td><td>1.596 × 10<sup>-3</sup></td><td>2.115 × 10<sup>3</sup></td></tr><tr><td>5.0</td><td>418</td><td>0.2898</td><td>0.1459</td><td>0.0192</td><td>2.493</td><td>*5.169 × 10<sup>-4</sup></td><td>*8.062 × 10<sup>3</sup></td></tr><tr><td>6.0</td><td>242</td><td>0.2898</td><td>0.1468</td><td>0.0193</td><td>2.991</td><td>*3.887 × 10<sup>-5</sup></td><td>*1.544 × 10<sup>5</sup></td></tr><tr><td>2.0</td><td>5028</td><td>0.3751</td><td>0.1833</td><td>0.0187</td><td>1.327</td><td>5.989 × 10<sup>-2</sup></td><td>1.113 × 10<sup>1</sup></td></tr><tr><td>3.0</td><td>1490</td><td>0.3751</td><td>0.1878</td><td>0.0191</td><td>1.991</td><td>5.832 × 10<sup>-3</sup></td><td>2.572 × 10<sup>2</sup></td></tr><tr><td>3.5</td><td>938</td><td>0.3751</td><td>0.1885</td><td>0.0192</td><td>2.323</td><td>1.287 × 10<sup>-3</sup></td><td>1.586 × 10<sup>3</sup></td></tr><tr><td>4.0</td><td>629</td><td>0.3751</td><td>0.1892</td><td>0.0193</td><td>2.655</td><td>*3.148 × 10<sup>-4</sup></td><td>*8.471 × 10<sup>3</sup></td></tr><tr><td>4.5</td><td>441</td><td>0.3751</td><td>0.1895</td><td>0.0193</td><td>2.987</td><td>*4.429 × 10<sup>-5</sup></td><td>*7.620 × 10<sup>4</sup></td></tr><tr><td>5.0</td><td>322</td><td>0.3751</td><td>0.1897</td><td>0.0193</td><td>3.319</td><td>*1.004 × 10<sup>-5</sup></td><td>*4.151 × 10<sup>5</sup></td></tr></tbody></table>
<table><caption>Table S2: Properties of the simulated systems (ANPs). $\sigma_N$: NP diameter. $N_N$: Number of NPs. $\rho_{m0} = M/V_0$: monomer number density in the neat system (independent of $\sigma_N$). $\phi_m = \pi\rho_N\sigma^3/6$: monomer volume fraction. $\phi_N = \pi\rho_N\sigma_N^3/6$: NP volume fraction. $C \equiv \sigma_N/\lambda$: Confinement parameter ($\lambda$ = crosslink localization length). $D_N$: NP diffusion coefficient. $\tau_D \equiv \sigma_N^2/6D_N$: Diffusion time. Values marked with * are upper and lower bounds, for $D_N$ and $\tau_D$ respectively, for systems which have not equilibrated.</caption><thead><tr><th>σ<sub>N</sub></th><th>N<sub>N</sub></th><th>ρ<sub>m0</sub></th><th>φ<sub>m</sub></th><th>φ<sub>N</sub></th><th>C</th><th>D<sub>N</sub></th><th>τ<sub>D</sub></th></tr></thead><tbody><tr><td>1.0</td><td>78038</td><td>0.1902</td><td>0.1054</td><td>0.0212</td><td>0.321</td><td>5.086 × 10<sup>-1</sup></td><td>3.277 × 10<sup>-1</sup></td></tr><tr><td>2.0</td><td>9755</td><td>0.1902</td><td>0.1027</td><td>0.0206</td><td>0.641</td><td>1.810 × 10<sup>-1</sup></td><td>3.682 × 10<sup>0</sup></td></tr><tr><td>3.0</td><td>2890</td><td>0.1902</td><td>0.1005</td><td>0.0202</td><td>0.962</td><td>9.593 × 10<sup>-2</sup></td><td>1.564 × 10<sup>1</sup></td></tr><tr><td>4.0</td><td>1219</td><td>0.1902</td><td>0.0995</td><td>0.0200</td><td>1.282</td><td>5.340 × 10<sup>-2</sup></td><td>4.994 × 10<sup>1</sup></td></tr><tr><td>5.0</td><td>624</td><td>0.1902</td><td>0.0991</td><td>0.0199</td><td>1.603</td><td>2.832 × 10<sup>-2</sup></td><td>1.471 × 10<sup>2</sup></td></tr><tr><td>6.0</td><td>361</td><td>0.1902</td><td>0.0988</td><td>0.0199</td><td>1.923</td><td>1.339 × 10<sup>-2</sup></td><td>4.482 × 10<sup>2</sup></td></tr><tr><td>7.0</td><td>228</td><td>0.1902</td><td>0.0987</td><td>0.0198</td><td>2.244</td><td>5.149 × 10<sup>-3</sup></td><td>1.586 × 10<sup>3</sup></td></tr><tr><td>8.0</td><td>152</td><td>0.1902</td><td>0.0985</td><td>0.0198</td><td>2.564</td><td>1.520 × 10<sup>-3</sup></td><td>7.017 × 10<sup>3</sup></td></tr><tr><td>9.0</td><td>107</td><td>0.1902</td><td>0.0983</td><td>0.0197</td><td>2.885</td><td>*3.423 × 10<sup>-4</sup></td><td>*3.943 × 10<sup>4</sup></td></tr><tr><td>10.0</td><td>78</td><td>0.1902</td><td>0.0981</td><td>0.0197</td><td>3.205</td><td>*1.112 × 10<sup>-4</sup></td><td>*1.499 × 10<sup>5</sup></td></tr><tr><td>1.0</td><td>52300</td><td>0.2898</td><td>0.1635</td><td>0.0216</td><td>0.499</td><td>3.060 × 10<sup>-1</sup></td><td>5.446 × 10<sup>-1</sup></td></tr><tr><td>2.0</td><td>6538</td><td>0.2898</td><td>0.1525</td><td>0.0201</td><td>0.997</td><td>9.891 × 10<sup>-2</sup></td><td>6.740 × 10<sup>0</sup></td></tr><tr><td>2.5</td><td>3347</td><td>0.2898</td><td>0.1515</td><td>0.0200</td><td>1.246</td><td>6.042 × 10<sup>-2</sup></td><td>1.724 × 10<sup>1</sup></td></tr><tr><td>3.0</td><td>1937</td><td>0.2898</td><td>0.1507</td><td>0.0199</td><td>1.496</td><td>3.703 × 10<sup>-2</sup></td><td>4.051 × 10<sup>1</sup></td></tr><tr><td>3.5</td><td>1220</td><td>0.2898</td><td>0.1505</td><td>0.0198</td><td>1.745</td><td>2.100 × 10<sup>-2</sup></td><td>9.721 × 10<sup>1</sup></td></tr><tr><td>4.0</td><td>817</td><td>0.2898</td><td>0.1502</td><td>0.0198</td><td>1.994</td><td>1.096 × 10<sup>-2</sup></td><td>2.433 × 10<sup>2</sup></td></tr><tr><td>4.5</td><td>574</td><td>0.2898</td><td>0.1499</td><td>0.0198</td><td>2.244</td><td>4.880 × 10<sup>-3</sup></td><td>6.916 × 10<sup>2</sup></td></tr><tr><td>5.0</td><td>418</td><td>0.2898</td><td>0.1498</td><td>0.0197</td><td>2.493</td><td>1.836 × 10<sup>-3</sup></td><td>2.269 × 10<sup>3</sup></td></tr><tr><td>1.0</td><td>40225</td><td>0.3751</td><td>0.2040</td><td>0.0208</td><td>0.664</td><td>2.234 × 10<sup>-1</sup></td><td>7.462 × 10<sup>-1</sup></td></tr><tr><td>2.0</td><td>5028</td><td>0.3751</td><td>0.1953</td><td>0.0199</td><td>1.327</td><td>5.412 × 10<sup>-2</sup></td><td>1.232 × 10<sup>1</sup></td></tr><tr><td>2.5</td><td>2574</td><td>0.3751</td><td>0.1940</td><td>0.0198</td><td>1.659</td><td>2.537 × 10<sup>-2</sup></td><td>4.105 × 10<sup>1</sup></td></tr><tr><td>3.0</td><td>1490</td><td>0.3751</td><td>0.1940</td><td>0.0198</td><td>1.991</td><td>1.059 × 10<sup>-2</sup></td><td>1.417 × 10<sup>2</sup></td></tr><tr><td>3.5</td><td>938</td><td>0.3751</td><td>0.1936</td><td>0.0197</td><td>2.323</td><td>3.132 × 10<sup>-3</sup></td><td>6.519 × 10<sup>2</sup></td></tr><tr><td>4.0</td><td>629</td><td>0.3751</td><td>0.1935</td><td>0.0197</td><td>2.655</td><td>*8.280 × 10<sup>-4</sup></td><td>*3.221 × 10<sup>3</sup></td></tr><tr><td>5.0</td><td>322</td><td>0.3751</td><td>0.1930</td><td>0.0196</td><td>3.319</td><td>*2.680 × 10<sup>-5</sup></td><td>*1.555 × 10<sup>5</sup></td></tr><tr><td>6.0</td><td>186</td><td>0.3751</td><td>0.1924</td><td>0.0196</td><td>3.982</td><td>*2.954 × 10<sup>-6</sup></td><td>*2.031 × 10<sup>6</sup></td></tr></tbody></table>
The mean chain length can be thus be calculated as$^{64}$ $N_x \equiv \langle n \rangle = 2(c^{-1} - 1)/3$ and for $c = 0.1$ one obtains $N_x = 6$.

## S2.2 Swelling and shrinking of the network upon the addition of NPs

<figure><img src="image_15.png" /></figure>

Figure S2: (a): Mean mesh size (most probable pore radius) $\xi \equiv r_{max}$ as a function of $\sigma_N/2\xi_0$, where $\xi_0$ is the mean mesh size of the neat system. The value of $\xi_0$ is 1.53, 1.03 and 0.79 for $\rho_{m0} = 0.190, 0.290$ and $0.375$, respectively. Filled symbols: RNP. Open symbols: ANP. Solid lines: $\xi_0$. (b): Ratio $\eta \equiv (V-V_0)/V_N$ as a function of the confinement parameter $C$. Dashed line: $\eta = \alpha + \beta C^\gamma$, with $\alpha = 2.16, \beta = -2.57$, and $\gamma = -0.61$. Solid line: same as in (c). (c): Representation of Eq. (S2). Solid line: Linear fit, $\eta^{1/3} = 0.93 + 0.80 \times C^{-1}$.

The addition of the NP to the neat system makes that the density of the latter changes. In Fig. S2a we show for all simulated systems with $\phi_{N0} = 0.02$ the mean mesh size $\xi \equiv r_{max}$ as a function of the NP diameter $\sigma_N$ divided by $2\xi_0$, with $\xi_0$ the mesh size of the neat
system (represented by the horizontal solid lines). Naively one could expect that if a volume fraction $\phi_{N0} = 0.02$ of NPs is added to the network, the mesh size would increase by a factor of $(1 + \phi_{N0})^{1/3} = 1.0067 \simeq 1$. For $\sigma_N/2\xi_0 > 1$, we observe indeed an increase of $\xi$ for both the ANPs and the RNPs, although somewhat larger than this value. However, for values of $\sigma_N/2\xi_0 < 1$, a reduction of $\xi$ is observed for the ANPs, i.e. the network is shrinking, whereas for the RNP a strong increase of $\xi$ is found (network swelling). Furthermore we find that for the ANPs, the maximum swelling is reached for $\sigma_N \simeq 2\xi_0$. These observations can be understood intuitively as follows: Small NPs are sufficiently small so that their excluded volume does not perturb the local mesh. However, they have, for a fixed NP volume fraction ($\phi_{N0} = 0.02$), a much larger total surface than large NPs, $S_N = \pi\sigma_N^2 N_N = 6\phi_{N0}/\sigma_N V_0$, and therefore give rise to a significant change in the total energy of the system (here $V_0$ is the volume of neat system). For the case of the ANP, the NP-monomer interaction gives rise to an effective attraction between the network strands, and thus to a shrinkage of the system, while for the RNP the effective interaction is repulsive, leading to an expansion of the network.⁶⁷

Large NPs also perturb the local mesh, as one can see from the secondary peak in Fig. 4, but have a limited effect on the network as a whole, since their total surface is much smaller at fixed volume fraction, as one can conclude from the fact that the main peak in Figs. 4a and b remains basically unperturbed, i.e. $\xi \simeq \xi_0$.

In Fig. S2b we report the quantity

$$ \eta \equiv (V - V_0)/V_N \qquad (S1) $$

as a function of the confinement parameter $C$, where $V$ is the volume of the filled system, and $V_N \equiv \pi\sigma_N^3 N_N/6$ is the total NP volume. The choice of this quantity is motivated by the fact that, if one approximates the network as an incompressible medium that expands
uniformly due to the excluded volume of the NPs, one would have $\eta \simeq 1$ (we recall that the NP volume fraction with respect to $V_0$, $\phi_{N0} \equiv V_N/V_0$, is constant, $\phi_{N0} = 0.02$). The graph demonstrates that for both types of NP, $\eta$ depends only on $C$, and it is independent of the network's density. Hence one can conclude that the confinement parameter controls the swelling of the network. This result is highly non-trivial in that it hints that $C$ is not only the relevant parameter for the dynamics, but also for the static properties of the system. If this conclusion is indeed true, it justifies the use of theoretical approaches that use the static properties of the system in order to predict its dynamics.

Fig. S2b also shows that the swelling/shrinking induced by the NPs is, in absolute value, significantly different from the simple estimate $V_N$ (horizontal line $\eta = 1$ in Fig. S2b) for intermediate and small values of $C$. For the RNPs, the behavior of $\eta(C)$ can be interpreted as follows: We have seen how the RNPs have a larger effective diameter than the ANPs. This is also seen in the behavior of the pore size distribution, since the presence of the NPs creates pores of diameter $\simeq \sigma_N + \delta$ (Fig. 4 in the main text and Sec. S2.6). The expansion caused by the RNPs can therefore estimated to be

$$V = V_0 + V_N \left(1 + \frac{\delta}{\sigma_N}\right)^3 \rightarrow \eta^{1/3} = \left(1 + \frac{\delta}{\lambda}C^{-1}\right). \qquad (S2)$$

In Fig. S2c, we plot $\eta^{1/3}$ as a function of $C^{-1}$ for the RNPs: The data are well described by a straight line, as expected from Eq. (S2); In particular, we find with good accuracy $\eta^{1/3} = 0.93 + 0.80 \times C^{-1}$, which implies $\delta = 0.80 \times \lambda$, which is in qualitative agreement with the pore size distribution data, since one finds that $\delta$ decreases (like $\lambda$) with increasing network density (Sec. S2.6). We note that the intercept is slightly smaller than unity, which is not surprising since the network is not incompressible.

For the ANPs, one can fit $\eta(C)$ by a (phenomenological) power-law, see solid line in Fig. S2b, the origin of which is currently unknown (the values of the parameters are given in
the caption of Fig. S2). For values of C corresponding to extreme confinement we find that $\eta$ is 1.4 for the RNP and very close to 1.0 for the ANP. We can however expect that these limits might depend on the details of the architecture of the gel and hence are not universal.

## S2.3 NP radial distribution function

In the main text we have discussed the radial distribution function for different values of $\sigma_N$ at a given density. For the sake of completeness we show in Fig. S3 $g_{NN}(r)$ for the other densities as well as further values of $\sigma_N$.

For the RNPs, the RDF displays a peak at the contact distance $r_c \equiv \sigma_N + (2^{1/6} - 1) = \sigma_N + 0.122$. Since the height of this peak is rather large, one can deduce that the probability to find two NPs in contact with each other is significantly higher than in a NP fluid with the same density $\rho_N$. The decay of $g_{NN}$ on the right hand side of the peak is basically exponential, as shown in the inset of Fig. S3e, where we report $|g_{NN}(r) - 1|$ as a function of $r - \sigma_N$ for $\rho_{m0} = 0.375$ in a semi-logarithmic representation. We note that for $\sigma_N > 2$ the shape of the peak, i.e. the decay length, depends only very weakly on the NP diameter, as evidenced by the fact that data for different values of $\sigma_N$ fall almost on the same master curve. This result is thus in harmony with the discussion in the main text that RNP have the tendency to form small clusters and that the probability to find such a cluster with a given size decreases exponentially with its size.

These figures also demonstrate that the height of the peak of the RDF increases with increasing $\sigma_N$. This trend is simply due to the fact that $g(r)$ is divided by the monomer density $\rho_m$, and $\rho_m$ decreases with increasing $\sigma_N$ since $\phi_{N0}$ is kept fixed. In order to better compare the peak heights, it is therefore useful to consider the coordination number, defined as$^{79}$

$$n(r) = 4\pi\rho_m \int_0^r g(r')r'^2dr'. \qquad (S3)$$
<figure><img src="image_16.png" /><figcaption>Figure S3: NP radial distribution function for the various simulated systems (see labels). Inset of Fig. (e): $|g_{NN}(r) - 1|$ as a function of $r - \sigma_N$. Dashed line: Exponential with slope -2.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure S4: Coordination number of the RNP (Eq. (S3)), divided by the ideal gas value (Eq. (S4)) for $\rho_{m0} = 0.290$ and different values of $\sigma_N$.</figcaption></figure>

For the interpretation of the $r$-dependence of $n(r)$ it is useful to divide it by the same quantity calculated for an ideal gas:

$$n_{\text{ig}}(r) \equiv \frac{4}{3}\pi\rho_m r^3 . \qquad (\text{S4})$$

Thus the ratio $n(r)/n_{\text{ig}}(r)$ gives us the coordination number of the NPs normalized by the one for an ideal gas. In Fig. S4, we show this ratio for the case of the RNP at $\rho_{m0} = 0.290$. One can see that for distances $1 < r/\sigma_N \lesssim 3$, this ratio is basically independent of $\sigma_N$ for $\sigma_N = 3.0, 3.5, 4.0$ and $5.0$, while it is smaller for $\sigma_N = 2.0$ and $\sigma_N = 6.0$. These two exceptions can be understood by arguing that small NP will be able to integrate into the matrix without much problem, i.e. the particles are well dispersed. On the other hand very large NP will have difficulty to find holes in the mesh that are sufficiently large to host several NP, leading to a reduction of the coordination number. (We can, however, not exclude the possibility that the systems with very large NP are not fully equilibrated.) Overall, the absolute value of $n(r)/n_{\text{ig}}(r)$ remains modest, signaling that the NP clustering is not severe, i.e., only small transient clusters are formed.

For the ANPs, panels d-f, the main peak of the RDF shifts to $r = \sigma_N + (2^{7/6} - 1) =$
$\sigma_N+1.245$, which corresponds to a configuration in which two neighboring NPs are separated by a polymer strand.⁶⁷ This is a clear signal of the fact that the NPs are well dispersed, since each NP is surrounded by a layer of polymers. Although a shoulder at the contact distance $r = \sigma_N + 0.122$ is still present, it reaches a height comparable with that of the main peak only for $\sigma_N = 1$: This is a signal of the fact that for $\sigma_N = 1$ even ANPs have a weak tendency to form clusters, which likely happens when several of them occupy the same hole in the mesh.

## S2.4 Monomer structure factor

<figure><img src="image_18.png" /><figcaption>Figure S5: Monomer structure factor. (a): $\rho_{m0} = 0.190$, RNP (b): $\rho_{m0} = 0.190$, ANP (c): $\rho_{m0} = 0.375$, RNP (d): $\rho_{m0} = 0.375$, ANP. Dashed line: Neat system ($N_m = 5 \times 10^4$).</figcaption></figure>

To study the structural properties of the system on large length scales it is useful to consider the static structure factor, as already defined in the main text. In Fig. S5 we show
the monomer structure factor $S_{mm}(q)$ for $\rho_{m0} = 0.190$ (a-b) and 0.375 (c-d), both for ANPs and RNPs. The data are compared with $S_{mm}(q)$ for a neat network with $N_m = 5 \times 10^4$ monomers. As also observed in the main text for $\rho_{m0} = 0.290$, one can see that $S_{mm}(q)$ does not change much upon the addition of the ANPs, except if $\sigma_N$ is very small. This is consistent with the fact that ANPs are well dispersed, thus inducing only moderate deformation on the local mesh. The RNPs, on the other hand, induce a noticeable deformation of the mesh, causing an increase of $S_{mm}(q)$ at low $q$ which signals the appearance of larger holes and the swelling of the network, in agreement with the results shown in Fig. S2b. This is due to the fact that RNPs partially aggregate to form small clusters, which deform the local mesh in a more substantial manner.

## S2.5 NP structure factor

In Fig. S6 we present the NP structure factor $S_{NN}(q)$ as a function of the rescaled wavevector $q\sigma_N/2\pi$ for $\rho_{m0} = 0.190$ and 0.375, for both ANPs and RNPs. The $q$-dependence of $S_{NN}(q)$ is qualitatively the same as the one which was observed for $\rho_{m0} = 0.290$ in the main text: For the RNPs, the data basically fall on a master curve, with a small peak at small $q$ which is due to partial NP clustering (as also observed from the RDF, see Sec. S2.3). In contrast to this one finds that for the higher density $\rho_{m0} = 0.375$ the height of this peak depends significantly on $\sigma_N$, and the peak height increases with decreasing $\sigma_N$. This $\sigma_N$-dependence can be understood by realizing that with increasing density the typical size of the holes in the mesh decreases and hence the NP will have an increasing hard time to find holes that can host several NP. This has the consequence that the clustering of the NP is suppressed, i.e. they become dispersed well in the gel, and hence the peak at small $q$ decreases.

Also for the ANPs the data fall on a master curve; in this case, the absence of a peak at $q \to 0$ confirms that the NPs are well dispersed. The only exception is the system with $\sigma_N = 1$, for which clustering is observed, in agreement with what was observed from the RDF in Sec. S2.3.
<figure><img src="image_19.png" /><figcaption>Figure S6: NP-NP structure factor as a function of the rescaled wavevector $q\sigma_N/2\pi$ for the RNPs (a,c) and the ANPs (b,d), for different values of $\sigma_N$ and $\rho_{m0}$.</figcaption></figure>
S2.6 Pore size distribution

<figure><img src="image_20.png" /><figcaption>Figure S7: Distribution of pore diameters, $\hat{P}(2r)$ for the RNPs (a,c) and the ANPs (b,d) for different values of $\sigma_N$ and $\rho_{m0}$.</figcaption></figure>

In Fig. S7 we report the pore size distribution (PSD) of the pore diameters, $\hat{P}(2r)$, for the RNPs (a,c) and the ANPs (b,d) for different values of $\sigma_N$ and $\rho_{m0} = 0.190$ (a,b) and 0.375 (c,d). We observe a behavior which is qualitatively similar to the one shown in the main text for $\rho_{m0} = 0.290$ (Fig. 4): The NPs create cavities in the network, resulting in peaks in the PSD located at $2r = \sigma_N$ for the ANPs and at $2r \approx \sigma_N + \delta$ for the RNPs, where the value of $\delta$ decreases slightly with increasing network density. For the RNPs, the peak is significantly broader than for the ANPs with the same diameter, due to the fact that the RNPs repel the strands and deform locally the network. We refer to the main text for a more detailed discussion of these features. Note that the distribution becomes more narrow if the density is increased, a dependence that indicates that on average denser systems have
smaller pores, in agreement with naive expectations.

# S3 Dynamics

## S3.1 Subdiffusive regime: onset and apparent exponent

<figure><img src="image_21.png" /><figcaption>Figure S8: MSD of the ANPs (a) and the RNPs (b), $\langle r_N^2(t) \rangle$, multiplied by the NP diameter $\sigma_N$, for $\rho_{m0} = 0.375$ and for different $\sigma_N$ values. Dashed horizontal line: $d^2\sigma$, with $d$ the tube diameter (Tab. 1). Dashed curves represent systems which have not reached the diffusive regime (see discussion in the text). Labels give the value of the confinement parameter $C$, Eq. (12).</figcaption></figure>

In Fig. S8, we present the MSD of the NPs, $\langle r_N^2(t) \rangle$, multiplied by the NP diameter $\sigma_N$ for the RNPs (a) and the ANPs (b), for $\rho_{m0} = 0.375$ and for different NP diameters. The theory of Cai et al.² predicts that $\langle r_N^2(t) \rangle$ reaches a plateau of height $d^2\sigma/\sigma_N$, with $\sigma = 1$ the monomer diameter and $d \simeq 1.94\xi(\rho_{m0})$ the tube diameter, i.e. the product plotted in the graph should have a plateau that is independent of $\sigma_N$. For our systems, the plateau is not fully developed even for the largest values of $C$, but one notes that the onset of the subdiffusive regime (which develops into a plateau for large values of $C$) occurs at $\langle r_N^2(t) \rangle \sigma_N \simeq d^2$ for all the simulated systems, in agreement with the theoretical prediction.

In Fig. S9 we show the apparent subdiffusive exponent $\beta(t)$ of the NPs, Eq. (13) in the main text, for the simulated systems. The time dependence of these curves are qualitatively similar to the one shown in the main text, Fig. 6: Initially, $\beta = 2$ (ballistic regime), then
<figure><img src="image_22.png" /><figcaption>Figure S9: Apparent subdiffusive exponent $\beta(t)$ of the NPs, Eq. (13) in the main text, for the simulated systems (see labels). Dashed curves represent systems which have not reached the diffusive regime. (a)-(c): RNP, (d)-(f): ANP.</figcaption></figure>
there is a sharp decrease as the NP enters in the subdiffusive regime, followed by a very slow transition to the diffusive regime $\beta = 1$. So these results show that the MSD of the NPs does not really show a plateau at intermediate times, which would correspond to $\beta = 0$, since the structural heterogeneity of the network makes that NPs start to leave their cage on time scales that are extremely broadly distributed.

## S3.2 NP mean-squared displacement

<figure><img src="image_23.png" /><figcaption>Figure S10: MSD of the RNPs (a,c), and the ANPs (b,d), $\langle r_N^2(t) \rangle$, and of the crosslinks, $\langle r_{cl}^2(t) \rangle$, for different NP diameters and $\rho_{m0}$. Dashed curves represent systems which have not reached the diffusive regime (see discussion in text). Labels give the value of the confinement parameter C, Eq. (12).</figcaption></figure>

In Fig. S10, we display the MSD of the RNPs (a,c), and the ANPs (b,d), $\langle r_N^2(t) \rangle$, and of the crosslinks, $\langle r_{cl}^2(t) \rangle$, for different NP diameters and $\rho_{m0} = 0.190$ (a,b) and 0.375 (c,d). These results are qualitatively the same as those reported in the main text for $\rho_{m0} = 0.375$:
Small NPs slip through the mesh, and go directly from the ballistic to the diffusive regime, whereas larger NPs are transiently trapped in the mesh, showing therefore a subdiffusive regime at intermediate times. At longer times, diffusive behavior is recovered as the NP is able to escape its local cage and diffuse *via* hopping motion. We refer to the main text for a more detailed discussion.

### S3.3 Non-Gaussian parameter

<figure><img src="image_24.png" /><figcaption>Figure S11: Non-Gaussian parameter for ANPs with different values of $\sigma_N$ and $\rho_{m0}$ (see labels), shown as a function of time (a,c) and of the rescaled time $t/\tau_D$, with $\tau_D \equiv \sigma_N^2/6D_N$ (b,d). Dashed curves represent systems which have not reached the diffusive regime.</figcaption></figure>

The time dependence of the non-Gaussian parameter $\alpha_2(t)$ is shown in Fig. S11 (ANP) and Fig. S12 (RNP). This quantity is shown as a function of time (left side of the figures) and as a function of the rescaled time $t/\tau_D$, with $\tau_D \equiv \sigma_N^2/6D_N$. For the ANPs, the results are qualitatively the same as those reported in the main text for the ANPs for
<figure><img src="image_25.png" /><figcaption>Figure S12: Non-Gaussian parameter for RNPs with different values of $\sigma_N$ and $\rho_{m0}$ (see labels), shown as a function of time (a,c,e) and of the rescaled time $t/\tau_D$, with $\tau_D \equiv \sigma_N^2/6D_N$ (b,d,f). Dashed curves represent systems which have not reached the diffusive regime.</figcaption></figure>
$\rho_{m0} = 0.375$, and thus we refer to the main text for a detailed discussion. Also for the RNPs the behavior of $\alpha_2$ is qualitatively similar to that of the same quantity for ANPs, with the difference that the curves are somewhat broader. This is consistent with the observation that the time needed to reach the diffusive regime is larger for RNPs, as one can see from Fig. S10. Furthermore one notices that the height of the maximum in $\alpha_2$ increases quickly with increasing density (keeping $C$ constant), which shows that the dynamical heterogeneities become more pronounced.

## References

(1) Dell, Z. E.; Schweizer, K. S. Theory of localization and activated hopping of nanoparticles in cross-linked networks and entangled polymer melts. *Macromolecules* **2014**, *47*, 405–414.

(2) Cai, L.-H.; Panyukov, S.; Rubinstein, M. Hopping diffusion of nanoparticles in polymer matrices. *Macromolecules* **2015**, *48*, 847–862.

(3) Huang, C. J.; Dostalek, J.; Knoll, W. Long range surface plasmon and hydrogel optical waveguide field-enhanced fluorescence biosensor with 3D hydrogel binding matrix: On the role of diffusion mass transfer. *Biosensors and Bioelectronics* **2010**, *26*, 1425–1431.

(4) Flier, B. M.; Baier, M. C.; Huber, J.; Müllen, K.; Mecking, S.; Zumbusch, A.; Wöll, D. Heterogeneous diffusion in thin polymer films as observed by high-temperature single-molecule fluorescence microscopy. *Journal of the American Chemical Society* **2012**, *134*, 480–488.

(5) Bhattacharya, S.; Sharma, D. K.; Saurabh, S.; De, S.; Sain, A.; Nandi, A.; Chowdhury, A. Plasticization of poly (vinylpyrrolidone) thin films under ambient humidity: Insight from single-molecule tracer diffusion dynamics. *The Journal of Physical Chemistry B* **2013**, *117*, 7771–7782.
(6) Riedinger, A.; Pernia Leal, M.; Deka, S. R.; George, C.; Franchini, I. R.; Falqui, A.; Cingolani, R.; Pellegrino, T. "Nanohybrids" based on pH-responsive hydrogels and inorganic nanoparticles for drug delivery and sensor applications. *Nano letters* **2011**, *11*, 3136-3141.

(7) Zhai, D.; Liu, B.; Shi, Y.; Pan, L.; Wang, Y.; Li, W.; Zhang, R.; Yu, G. Highly sensitive glucose sensor based on Pt nanoparticle/polyaniline hydrogel heterostructures. *ACS nano* **2013**, *7*, 3540–3546.

(8) Amblard, F.; Maggs, A. C.; Yurke, B.; Pargellis, A. N.; Leibler, S. Subdiffusion and anomalous local viscoelasticity in actin networks. *Physical Review Letters* **1996**, *77*, 4470.

(9) Wong, I.; Gardel, M.; Reichman, D.; Weeks, E. R.; Valentine, M.; Bausch, A.; Weitz, D. A. Anomalous diffusion probes microstructure dynamics of entangled F-actin networks. *Physical Review Letters* **2004**, *92*, 178101.

(10) Fritsch, C. C.; Langowski, J. Anomalous diffusion in the interphase cell nucleus: the effect of spatial correlations of chromatin. *The Journal of Chemical Physics* **2010**, *133*, 07B602.

(11) Stylianopoulos, T.; Poh, M.-Z.; Insin, N.; Bawendi, M. G.; Fukumura, D.; Munn, L. L.; Jain, R. K. Diffusion of particles in the extracellular matrix: the effect of repulsive electrostatic interactions. *Biophysical journal* **2010**, *99*, 1342–1349.

(12) Peulen, T.-O.; Wilkinson, K. J. Diffusion of nanoparticles in a biofilm. *Environmental science & technology* **2011**, *45*, 3367–3373.

(13) Yu, M.; Xu, L.; Tian, F.; Su, Q.; Zheng, N.; Yang, Y.; Wang, J.; Wang, A.; Zhu, C.; Guo, S. Rapid transport of deformation-tuned nanoparticles across biological hydrogels and cellular barriers. *Nature communications* **2018**, *9*, 1–11.
(14) Cherstvy, A. G.; Thapa, S.; Wagner, C. E.; Metzler, R. Non-Gaussian, non-ergodic, and non-Fickian diffusion of tracers in mucin hydrogels. *Soft Matter* **2019**, *15*, 2526–2551.

(15) Burla, F.; Sentjabrskaja, T.; Pletikapic, G.; Van Beugen, J.; Koenderink, G. H. Particle diffusion in extracellular hydrogels. *Soft Matter* **2020**, *16*, 1366–1376.

(16) Debets, V. E.; Janssen, L. M.; Šarić, A. Characterising the Diffusion of Biological Nanoparticles on Fluid and Elastic Membranes. *bioRxiv* **2020**,

(17) Cho, K.; Wang, X. U.; Nie, S.; Shin, D. M.; Others, Therapeutic nanoparticles for drug delivery in cancer. *Clinical Cancer Research* **2008**, *14*, 1310–1316.

(18) Ward, M. A.; Georgiou, T. K. Thermoresponsive polymers for biomedical applications. *Polymers* **2011**, *3*, 1215–1242.

(19) Cai, L.-H.; Panyukov, S.; Rubinstein, M. Mobility of nonsticky nanoparticles in polymer liquids. *Macromolecules* **2011**, *44*, 7853–7863.

(20) Egorov, S. A. Anomalous nanoparticle diffusion in polymer solutions and melts: A mode-coupling theory study. *The Journal of Chemical Physics* **2011**, *134*, 84903.

(21) Yamamoto, U.; Schweizer, K. S. Theory of nanoparticle diffusion in unentangled and entangled polymer melts. *The Journal of Chemical Physics* **2011**, *135*, 224902.

(22) Yamamoto, U.; Schweizer, K. S. Microscopic theory of the long-time diffusivity and intermediate-time anomalous transport of a nanoparticle in polymer melts. *Macromolecules* **2014**, *48*, 152–163.

(23) Dong, Y.; Feng, X.; Zhao, N.; Hou, Z. Diffusion of nanoparticles in semidilute polymer solutions: A mode-coupling theory study. *The Journal of Chemical Physics* **2015**, *143*.

(24) Bedrov, D.; Smith, G. D.; Smith, J. S. Matrix-induced nanoparticle interactions in a polymer melt: A molecular dynamics simulation study. *The Journal of Chemical Physics* **2003**, *119*, 10438–10447.
(25) Liu, J.; Cao, D.; Zhang, L. Molecular dynamics study on nanoparticle diffusion in polymer melts: a test of the Stokes-Einstein law. *Journal of Physical Chemistry C* **2008**, *112*, 6653-6661.

(26) Kalathi, J. T.; Yamamoto, U.; Schweizer, K. S.; Grest, G. S.; Kumar, S. K. Nanoparticle Diffusion in Polymer Nanocomposites. *Physical Review Letters* **2014**, *108301*, 1-5.

(27) Patti, A. Molecular dynamics of spherical nanoparticles in dense polymer melts. *The Journal of Physical Chemistry B* **2014**, *118*, 3731-3742.

(28) Li, Y.; Kröger, M.; Liu, W. K. Dynamic structure of unentangled polymer chains in the vicinity of non-attractive nanoparticles. *Soft Matter* **2014**, *10*, 1723.

(29) Kalathi, J. T.; Kumar, S. K.; Rubinstein, M.; Grest, G. S. Rouse mode analysis of chain relaxation in polymer nanocomposites. *Soft Matter* **2015**, *11*, 4123-4132.

(30) Volgin, I. V.; Larin, S. V.; Abad, E.; Lyulin, S. V. Molecular Dynamics Simulations of Fullerene Diffusion in Polymer Melts. *Macromolecules* **2017**, *50*, 2207-2218.

(31) Karatrantos, A.; Composto, R. J.; Winey, K. I.; Clarke, N. Polymer and spherical nanoparticle diffusion in nanocomposites. *The Journal of Chemical Physics* **2017**, *146*, 203331.

(32) Chen, R.; Poling-Skutvik, R.; Nikoubashman, A.; Howard, M. P.; Conrad, J. C.; Palmer, J. C. Coupling of nanoparticle dynamics to polymer center-of-mass motion in semidilute polymer solutions. *Macromolecules* **2018**, *51*, 1865-1872.

(33) Chen, A.; Zhao, N.; Hou, Z. The effect of hydrodynamic interactions on nanoparticle diffusion in polymer solutions: a multiparticle collision dynamics study. *Soft Matter* **2017**, *13*, 8625-8635.

(34) Yamamoto, U.; Carrillo, J.-M. Y.; Bocharova, V.; Sokolov, A. P.; Sumpter, B. G.;
Schweizer, K. S. Theory and simulation of attractive nanoparticle transport in polymer melts. *Macromolecules* **2018**, *51*, 2258–2267.

(35) Du, Y.; Jiang, H.; Hou, Z. Study of active Brownian particle diffusion in polymer solutions. *Soft Matter* **2019**, *15*, 2020–2031.

(36) Netz, P. A.; Dorfmüller, T. Computer simulation studies of diffusion in gels: Model structures. *The Journal of Chemical Physics* **1997**, *107*, 9221–9233.

(37) Sonnenburg, J.; Gao, J.; Weiner, J. Molecular dynamics simulations of gas diffusion through polymer networks. *Macromolecules* **1990**, *23*, 4653–4657.

(38) Licinio, P.; Teixeira, A. Anomalous diffusion of ideal polymer networks. *Physical Review E* **1997**, *56*, 631.

(39) Zhou, H.; Chen, S. B. Brownian dynamics simulation of tracer diffusion in a cross-linked network. *Physical Review E* **2009**, *79*, 021801.

(40) Godec, A.; Bauer, M.; Metzler, R. Collective dynamics effect transient subdiffusion of inert tracers in flexible gel networks. *New Journal of Physics* **2014**, *16*, 092002.

(41) Kamerlin, N.; Elvingson, C. Tracer diffusion in a polymer gel: simulations of static and dynamic 3D networks using spherical boundary conditions. *Journal of Physics: Condensed Matter* **2016**, *28*, 475101.

(42) Kumar, P.; Theeyancheri, L.; Chaki, S.; Chakrabarti, R. Transport of probe particles in a polymer network: effects of probe size, network rigidity and probe–polymer interaction. *Soft matter* **2019**, *15*, 8992–9002.

(43) Chen, Y.; Ma, R.; Qian, X.; Zhang, R.; Huang, X.; Xu, H.; Zhou, M.; Liu, J. Nanoparticle Mobility within Permanently Cross-Linked Polymer Networks. *Macromolecules* **2020**,
(44) Cho, H. W.; Kim, H.; Sung, B. J.; Kim, J. S. Tracer Diffusion in Tightly-Meshed Homogeneous Polymer Networks: A Brownian Dynamics Simulation Study. *Polymers* **2020**, *12*, 2067.

(45) Xu, Z.; Dai, X.; Bu, X.; Yang, Y.; Zhang, X.; Man, X.; Zhang, X.; Doi, M.; Yan, L.-T. Enhanced Heterogeneous Diffusion of Nanoparticles in Semiflexible Networks. *ACS nano* **2021**, *15*, 4608–4616.

(46) Peppas, N. A.; Moynihan, H. J.; Lucht, L. M. The structure of highly crosslinked poly (2-hydroxyethyl methacrylate) hydrogels. *Journal of biomedical materials research* **1985**, *19*, 397–411.

(47) Gehman, S. Network chain distribution and strength of vulcanizates. *Rubber chemistry and technology* **1969**, *42*, 659–665.

(48) Falcao, A.; Pedersen, J. S.; Mortensen, K. Structure of randomly crosslinked poly (dimethylsiloxane) networks produced by electron irradiation. *Macromolecules* **1993**, *26*, 5350–5364.

(49) Rubinstein, M.; Colby, R. H. *Polymer physics*; Oxford University Press New York, 2003.

(50) Schweizer, K. S.; Curro, J. G. Integral equation theories of the structure, thermodynamics, and phase transitions of polymer fluids. *Advances in Chemical Physics* **1997**, *98*, 1–142.

(51) Schweizer, K. S.; Saltzman, E. J. Entropic barriers, activated hopping, and the glass transition in colloidal suspensions. *The Journal of Chemical Physics* **2003**, *119*, 1181–1196.

(52) Parrish, E.; Caporizzo, M. A.; Composto, R. J. Network confinement and heterogeneity
slows nanoparticle diffusion in polymer gels. *The Journal of Chemical Physics* **2017**, *146*, 203318.

(53) Parrish, E.; Seeger, S. C.; Composto, R. J. Temperature-Dependent Nanoparticle Dynamics in Poly (N-isopropylacrylamide) Gels. *Macromolecules* **2018**, *51*, 3597–3607.

(54) Anderson, S. J.; Matsuda, C.; Garamella, J.; Peddireddy, K. R.; Robertson-Anderson, R. M.; McGorty, R. Filament Rigidity Vies with Mesh Size in Determining Anomalous Diffusion in Cytoskeleton. *Biomacromolecules* **2019**, *20*, 4380–4388.

(55) Wang, Y.; Li, Z.; Ouyang, J.; Karniadakis, G. E. Controlled release of entrapped nanoparticles from thermoresponsive hydrogels with tunable network characteristics. *Soft Matter* **2020**, *16*, 4756–4766.

(56) Gnan, N.; Rovigatti, L.; Bergman, M.; Zaccarelli, E. In silico synthesis of microgel particles. *Macromolecules* **2017**, *50*, 8777–8786.

(57) Rovigatti, L.; Gnan, N.; Tavagnacco, L.; Moreno, A. J.; Zaccarelli, E. Numerical modelling of non-ionic microgels: an overview. *Soft Matter* **2019**, *15*, 1108–1119.

(58) Ninarello, A.; Crassous, J. J.; Paloli, D.; Camerin, F.; Gnan, N.; Rovigatti, L.; Schurtenberger, P.; Zaccarelli, E. Modeling Microgels with a Controlled Structure across the Volume Phase Transition. *Macromolecules* **2019**, *52*, 7584–7592.

(59) Curro, J. G.; Pincus, P. A theoretical basis for viscoelastic relaxation of elastomers in the long-time limit. *Macromolecules* **1983**, *16*, 559–562.

(60) Duering, E. R.; Kremer, K.; Grest, G. S. Structure and relaxation of end-linked polymer networks. *The Journal of Chemical Physics* **1994**, *101*, 8169–8192.

(61) Sorichetti, V.; Ninarello, A.; Ruiz-Franco, J. M.; Hugouvieux, V.; Kob, W.; Zaccarelli, E.; Rovigatti, L. Effect of Chain Polydispersity on the Elasticity of Disordered Polymer Networks. *Macromolecules* **2021**, *54*, 3769–3779.
(62) Flory, P. J. *Principles of polymer chemistry*; Cornell University Press, 1953.

(63) Stockmayer, W. H. Theory of molecular size distribution and gel formation in branched-chain polymers. *The Journal of Chemical Physics* **1943**, *11*, 45-55.

(64) Rovigatti, L.; Gnan, N.; Zaccarelli, E. Internal structure and swelling behaviour of in silico microgel particles. *Journal of Physics: Condensed Matter* **2017**, *30*, 044001.

(65) Kremer, K.; Grest, G. S. Dynamics of entangled linear polymer melts: A molecular-dynamics simulation. *The Journal of Chemical Physics* **1990**, *92*, 5057-5086.

(66) Weeks, J. D.; Chandler, D.; Andersen, H. C. Role of repulsive forces in determining the equilibrium structure of simple liquids. *The Journal of Chemical Physics* **1971**, *54*, 5237-5247.

(67) Sorichetti, V.; Hugouvieux, V.; Kob, W. Structure and dynamics of a polymer-nanoparticle composite: Effect of nanoparticle size and volume fraction. *Macromolecules* **2018**, *51*, 5375-5391.

(68) Chen, R.; Poling-Skutvik, R.; Nikoubashman, A.; Howard, M. P.; Conrad, J. C.; Palmer, J. C. Coupling of nanoparticle dynamics to polymer center-of-mass motion in semidilute polymer solutions. *Macromolecules* **2018**, *51*, 1865-1872.

(69) LAMMPS website. https://lammps.sandia.gov, [accessed 4 October 2019].

(70) Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. *Journal of Computational Physics* **1995**, *117*, 1-19.

(71) Auhl, R.; Everaers, R.; Grest, G. S.; Kremer, K.; Plimpton, S. J. Equilibration of long chain polymer melts in computer simulations. *The Journal of Chemical Physics* **2003**, *119*, 12718-12728.

(72) Tuckerman, M. *Statistical mechanics: theory and molecular simulation*; Oxford university press, 2010.
(73) Schneider, T.; Stoll, E. Molecular-dynamics study of a three-dimensional one-component model for distortive phase transitions. *Physical Review B* **1978**, *17*, 1302.

(74) Binder, K.; Kob, W. *Glassy materials and disordered solids: An introduction to their statistical mechanics*; World scientific, 2011.

(75) Hooper, J. B.; Schweizer, K. S. Theory of Phase Separation in Polymer Nanocomposites. *Macromolecules* **2006**, 5133–5142.

(76) Liu, J.; Gao, Y.; Cao, D.; Zhang, L.; Guo, Z. Nanoparticle dispersion and aggregation in polymer nanocomposites: Insights from molecular dynamics simulation. *Langmuir* **2011**, *27*, 7926–7933.

(77) Meng, D.; Kumar, S. K.; Cheng, S.; Grest, G. S. Simulating the miscibility of nanoparticles and polymer melts. *Soft Matter* **2013**, *9*, 5417–5427.

(78) Karatrantos, A.; Clarke, N.; Composto, R. J.; Winey, K. I. Polymer conformations in polymer nanocomposites containing spherical nanoparticles. *Soft Matter* **2015**, *11*, 382–388.

(79) Hansen, J.-P.; McDonald, I. R. *Theory of simple liquids*; Elsevier, 1990.

(80) Gelb, L. D.; Gubbins, K. E. Pore size distributions in porous glasses: a computer simulation study. *Langmuir* **1999**, *15*, 305–308.

(81) Bhattacharya, S.; Gubbins, K. E. Fast method for computing pore size distributions of model materials. *Langmuir* **2006**, *22*, 7726–7731.

(82) Sorichetti, V.; Hugouvieux, V.; Kob, W. Determining the Mesh Size of Polymer Solutions via the Pore Size Distribution. *Macromolecules* **2020**, *53*, 2568–2581.

(83) Roldán-Vargas, S.; Rovigatti, L.; Sciortino, F. Connectivity, dynamics, and structure in a tetrahedral network liquid. *Soft Matter* **2017**, *13*, 514–530.
(84) Zaccarelli, E.; Buldyrev, S.; La Nave, E.; Moreno, A.; Saika-Voivod, I.; Sciortino, F.; Tartaglia, P. Model for reversible colloidal gelation. *Physical Review Letters* **2005**, *94*, 218301.

(85) Hsu, H.-P.; Kremer, K. Static and dynamic properties of large polymer melts in equilibrium. *The Journal of Chemical Physics* **2016**, *144*, 154907.

(86) Hoy, R. S.; Foteinopoulou, K.; Kröger, M. Topological analysis of polymeric melts: Chain-length effects and fast-converging estimators for entanglement length. *Physical Review E* **2009**, *80*, 31803.

(87) Kob, W.; Andersen, H. C. Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture. I: The van Hove correlation function. *Physical Review E* **1995**, *51*, 4626.