# Structure and Dynamics of a Polymer-Nanoparticle Composite: Effect of Nanoparticle Size and Volume Fraction

Valerio Sorichetti, Virginie Hugouvieux, Walter Kob

► To cite this version:

Valerio Sorichetti, Virginie Hugouvieux, Walter Kob. Structure and Dynamics of a Polymer-Nanoparticle Composite: Effect of Nanoparticle Size and Volume Fraction. Macromolecules, 2018, 51 (14), pp.5375-5391. 10.1021/acs.macromol.8b00840. hal-01870804

HAL Id: hal-01870804

https://hal.science/hal-01870804v1

Submitted on 10 Feb 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
# Structure and dynamics of a polymer-nanoparticle composite: Effect of nanoparticle size and volume fraction

Valerio Sorichetti

Laboratoire Charles Coulomb (L2C), Univ Montpellier,
CNRS, Montpellier, France and SPO, Univ Montpellier,
INRA, Montpellier SupAgro, Montpellier, France

Virginie Hugouvieux

SPO, Univ Montpellier, INRA, Montpellier SupAgro, Montpellier, France

Walter Kob

Laboratoire Charles Coulomb (L2C), Univ Montpellier, CNRS, Montpellier, France

(Dated: September 18, 2018)

We use molecular dynamics simulations to study a semidilute, unentangled polymer solution containing well dispersed, weakly attractive nanoparticles (NP) of size ($\sigma_N$) smaller than the polymer radius of gyration $R_g$. We find that if $\sigma_N$ is larger than the monomer size the polymers swell, while smaller NPs cause chain contraction. The diffusion coefficient of polymer chains ($D_p$) and NPs ($D_N$) decreases if the volume fraction $\phi_N$ is increased. The decrease of $D_p$ can be well described in terms of a dynamic confinement parameter, while $D_N$ shows a more complex dependence on $\sigma_N$, which results from an interplay between energetic and entropic effects. When $\phi_N$ exceeds a $\sigma_N$-dependent value, the NPs are no longer well dispersed and $D_N$ and $D_p$ increase if $\phi_N$ is increased.

## I. INTRODUCTION

Understanding the motion of nanoparticles (NP) and macromolecules in complex fluids, such as polymer solutions and melts, is a problem of broad importance, with applications to many different fields. In material science, understanding how NPs move in a polymer matrix is fundamental for the production of nanocomposites with mechanical, thermal, optical or electrical properties superior to those of pure polymeric materials [1–3]. In biophysics, the dynamics of macromolecules in the cytoplasmic environment can have a strong influence on cellular functions, such as enzymatic reactions and self-assembly of cellular structures [4–6]. Also in medicine there is a growing interest in the topic, with the objective to develop new and more efficient forms of NP-mediated drug delivery [7, 8], a practice which is already in use for cancer treatment [9, 10].

In the past years, a lot of attention has been dedicated to the study of the motion of polymers and NPs in polymer solutions and melts, using theoretical [11–15] and experimental [16–25] approaches, as well as computer simulations [26–32]. Polymer-NP mixtures represent a tough challenge for theoretical physics mainly because of the large number of different length scales present: the NP diameter $\sigma_N$, the monomer diameter $\sigma$, the Kuhn length $\ell_K$, the radius of gyration $R_g$, and, in the case of concentrated solutions and melts, the mesh size $\xi$ and the diameter of the Edwards tube a [33, 34]. The behavior of these systems strongly depends on the length and time scale at which they are probed, and in certain conditions it is possible to observe interesting dynamical phenomena, like anomalous diffusion [11, 24, 25, 35, 36] or the breakdown of the Stokes-Einstein relation [11, 16–

18, 27, 28, 37, 38]. Also the interaction between the different components, which depend on the microscopic details, can have a great impact on the system's structure and dynamics [13, 29, 39–42]. Understanding how all these factors affect the static and dynamic properties of the NPs in polymer solutions and melts is thus crucial for practical applications.

When studying polymer-NP mixtures, two main regimes can be identified depending on the NP diameter $\sigma_N$: the “colloid limit”, where the polymers are much smaller than the NPs ($2R_g/\sigma_N \ll 1$) and the “protein limit” or “nanoparticle limit” [43], where the size of the polymers is larger or comparable to that of the NPs ($2R_g/\sigma_N \gtrsim 1$). The colloid limit has been studied extensively and it is nowadays well understood in terms of effective depletion pair potentials [44, 45]. The protein limit, on the other hand, is much more problematic, since an accurate description in terms of effective pair potentials is not possible [46, 47]. In the present work, we will focus on the protein limit, using molecular dynamics simulations of a coarse-grained model.

With few exceptions [27, 30, 32], most of the previous simulation studies of polymer-NP mixtures have focused on the dilute NP regime, in which the NPs can be assumed not to interact with each other and the properties of the polymer solution are expected to be unchanged by the presence of the NPs. Thus, the purpose of the present work is to study the diffusion of polymers and NPs in an unentangled, semidilute polymer solution in a wide range of NP volume fractions and NP diameters, up to values where the interaction between NPs cannot be neglected.

The paper is organized as follows: In Section II, the model and the simulation method and details are presented. In Section III we discuss the structural prop-
erties of the polymer-NP solution for different NP sizes and volume fractions, with a special focus on the structure of single polymer chains. In Section IV, we study the dynamical properties of the system in the presence of good NP dispersion, and in particular the diffusion coefficient of the centers of mass of the chains and of the NPs. Finally, in Section V we investigate the behavior of the system at high NP volume fraction, where the NP dispersion becomes progressively poorer until large polymer-free regions are formed. We conclude with a summary in Section VI.

## II. MODEL AND SIMULATION METHOD

We performed NVT molecular dynamics simulations of a system of $N_p = 500$ polymer chains of length (degree of polymerization) $N = 100$ and a variable number $N_N$ of nanoparticles of different diameters $\sigma_N$. To simulate the polymer chains, we used the bead-spring model of Kremer and Grest [48]. All monomers interact via a Weeks-Chandler-Andersen (WCA) potential [49],

$$ 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} \quad (1) $$

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

$$ U_{\text{bond}}(r) = -\frac{1}{2}kr_0^2 \ln[1 - (r/r_0)^2], \quad (2) $$

where $k = 30\epsilon/\sigma^2$ and $r_0 = 1.5\sigma$. With this choice of parameters the bond length at the minimum of the potential is $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 [48].

In the following, all quantities are given in Lennard-Jones (LJ) reduced units. The units of energy, length and mass are respectively $\epsilon$, $\sigma$ and $m$, where $\epsilon$, and $\sigma$ are defined by Eq. (1) and $m$ is the mass of a monomer. The units of temperature, pressure, volume fraction and time are respectively $[T] = \epsilon/k_B$, $[P] = \epsilon\sigma^{-3}$, $[\phi] = \sigma^{-3}$ and $[t] = \sqrt{m\sigma^2/\epsilon}$.

For the interaction potentials involving the NPs, we use an “expanded Lennard-Jones” (expanded LJ) potential, which is a LJ potential shifted to the right by a quantity $\Delta_{ij}$: Thus, as opposed to the standard LJ potential, in the expanded LJ potential the “softness” (slope) of the potential does not change when the NP size varies, as one can see in Fig. 1. Since experiments have shown that the thickness of the interfacial region surrounding a NP in a polymer matrix changes only weakly with the size of the NP [50], the expanded LJ is a better choice than the standard LJ potential when simulating polymer-NP mixtures [26, 27, 51].

<figure><img src="image_2.png" /><figcaption>FIG. 1. Monomer-NP potential $U_{Nm}(r)$ (expanded LJ) compared to a standard LJ potential $U_{LJ}(r) = 4[(\sigma_{Nm}/r)^{12} - (\sigma_{Nm}/r)^{6}]$, where $\sigma_{Nm} = (\sigma_N + \sigma)/2$. The expanded LJ potential is cut and shifted at $2.5\sigma + (\sigma_N - \sigma)/2$ and the LJ potential is cut and shifted at $2.5\sigma_{Nm}$. For $\sigma_N = \sigma$, the two potentials coincide.</figcaption></figure>

The interaction between monomers and NPs and between two NPs has thus the following general form:

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

where $\Delta_{Nm} = (\sigma_N + \sigma)/2 - \sigma = (\sigma_N - \sigma)/2$ and $\Delta_{NN} = \sigma_N - \sigma$. The quantity $\delta_{ij}$ is such that $U_{ij}(r_{ij}^c + \Delta_{ij}) = 0$. The cutoff distances are $r_{Nm}^c = 2.5$ for the NP-monomer interaction and $r_{NN}^c = 2^{1/6}$ for the NP-NP interaction. The interaction between monomers and NPs is therefore attractive, while the interaction between NPs is purely repulsive. A moderate attractive interaction between polymers and NPs is required in order to prevent aggregation (and eventually phase separation) of the NPs [40–42, 52].

In this study, we consider NP diameters $\sigma_N = 1, 2, 3, 4, 5$, and $7\sigma$. 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$.

We define the NP volume fraction as $\phi_N = \pi\sigma_N^3 N_N/6V$, where $V$ is the total volume of the simulation box; the monomer volume fraction $\phi_m$ is defined in an analogous way. In our simulations, $\phi_m$ is larger than the overlap volume fraction [33], which can be estimated from the polymer’s radius of gyration at infinite dilution (see below) and for the pure polymer system has the value $\phi_m^* = 2.98 \cdot 10^{-2}$. Moreover, since the entanglement length $N_e$ [33] for this model is $N_e \approx 85$ at $\rho_m = 0.85$ [53] and since $N_e$ scales approximately as $\rho_m^{-2}$ [54–56], we are always in the unentangled regime [57].

All the simulations were carried out using the LAMMPS software [58]. The simulation box is cubic
and periodic boundary conditions are applied in all directions. The initial configurations are prepared by randomly placing the polymers and the NPs in the box; initially, the NPs have diameter equal to that of the monomers ($\sigma_N = \sigma$) and overlaps between particles are allowed. The overlaps are then removed by using a soft potential whose strength is increased over a short amount of time ("fast push-off" method [59]). After the overlaps are removed, the diameter of the NPs is gradually increased until the desired value is reached, and finally the system is allowed to adjust its density until we reach pressure $P = 0.1$ at temperature $T = 1.0$. In the pure polymer systems, these parameters correspond to a monomer volume fraction $\phi_m = 0.147$ (monomer density $\rho_m = 0.280$). Finally, we switch to the NVT ensemble and perform an equilibration run before starting the production run. During the NVT simulations, the pressure fluctuations are always less than 14%.

The length of both the equilibration and the production runs is $10^8 \delta t = 3 \cdot 10^5$, where $\delta t = 3 \cdot 10^{-3}$ is the integration time step. In all cases, we verified that during the equilibration runs the NPs (polymers) diffused on average over a distance equal to several times their diameter (radius of gyration), and that their motion became diffusive (see below).

Both during the NPT and the NVT runs, the temperature is kept fixed by means of a Langevin thermostat, so that the force experienced by a particle $i$ (monomer or NP) is

$$m_i \ddot{\mathbf{r}}_i = -\nabla_i U(\{\mathbf{r}_k\}) - m_i \Gamma_i \dot{\mathbf{r}}_i + \sqrt{2m_i \Gamma_i k_B T} \zeta(t), \quad (4)$$

where $\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 of the right side of Eq. (S1) represents viscous damping, with $\Gamma_i$ the friction coefficient, and the last term is a random, uncorrelated force representing the collisions with solvent particles. The Langevin thermostat acts therefore as an implicit solvent, in which every particle interacts independently with the solvent "molecules", but hydrodynamic interactions between solute particles are not accounted for. We note that it has been recently pointed out that hydrodynamic interactions can affect the long-time dynamics of NPs in a polymer solution even at high monomer volume fractions [60], an observation which warrants further investigation.

The damping constant for the monomers is $\Gamma_m = 0.1$, while that of the NPs is chosen by imposing that the viscosity of the pure solvent calculated via the Stokes formula, $\eta_s = \Gamma_i m_i / 3\pi\sigma_i$, is the same for both free monomers and NPs. Therefore we have $\Gamma_N = \Gamma_m \sigma_N / m_N$. For a more detailed discussion, see Sec. SI in the S.I.

Additional details about the simulations can be found in Tab. S2 in the S.I.

<figure><img src="image_3.png" /><figcaption>FIG. 2. Snapshots of systems containing NPs of diameter $\sigma_N = 4$ at different volume fractions (yellow spheres: monomers, blue spheres: NPs). In systems (a) and (b), the dispersion of NPs is good, while in (c) and (d) it is poor. In system (d) the formation of large polymer-free regions is evident. The length of the simulation box edges are respectively 57.45 (a), 59.15 (b), 67.13 (c), and 93.57 (d).</figcaption></figure>

### III. STRUCTURE

#### A. Nanoparticles

To give a feeling of what the simulated system looks like, we show in Fig. 2 some snapshots for $\sigma_N = 4$ and different values of the NP volume fraction $\phi_N$. We can see how the NP dispersion, which is initially good (Figs. 2a-b), becomes progressively poorer as $\phi_N$ is increased (Fig. 2c), until eventually large polymer-free regions are formed (Fig. 2d). In order to characterize the structure of the systems when $\sigma_N$ and $\phi_N$ are varied, we start by analyzing some basic quantities, such as the radial distribution function $g(r)$ and the structure factor $S(q)$.

The radial distribution function $g(r)$ can be obtained from the pair correlation function $g(\mathbf{r})$ by performing a spherical average [61]. We recall that the pair correlation function of a system of $M$ particles with number density $\rho$ is defined as [62]

$$\rho g(\mathbf{r}) = \frac{1}{M} \sum_{\substack{k=1 \\ j \neq k}}^{M} \langle \delta(\mathbf{r} + \mathbf{r}_k - \mathbf{r}_j) \rangle, \quad (5)$$

where $\langle \cdot \rangle$ denotes the thermodynamic average.
<figure><img src="image_4.png" /><figcaption>FIG. 3. (color online) NP-NP radial distribution function at different NP volume fractions $\phi_N$ for $\sigma_N = 2$ (a) and $\sigma_N = 4$ (b). The contact peak at $r_c = \sigma_N + 0.122$ corresponds to a configuration in which two NPs are touching, while the secondary peak at $r_s = \sigma_N + 1.245$ represents a configuration in which two NPs are separated by a polymer strand (a schematic representation of the two cases is shown in (a)): At low and intermediate $\phi_N$, this last configuration is favored.</figcaption></figure>

Figure 3 shows the NP-NP radial distribution function $g_{NN}(r)$ for $\sigma_N = 2$ and 4 and different values of the NP volume fraction $\phi_N$. For low values of $\phi_N$, $g_{NN}$ shows a peak at $r_s = \sigma_N + (2^{1/6} - 1) = \sigma_N + 1.245$, which corresponds to twice the distance at the minimum of the monomer-NP potential. This indicates that the NPs are well dispersed in the polymer solution and configurations in which two neighboring NPs are separated by a polymer strand are favored (this kind of configuration is schematically represented in Fig. 3a). We call this peak secondary peak.

When $\phi_N$ increases, another peak appears at $r_c = \sigma_N + (2^{1/6} - 1) = \sigma_N + 0.122$, which corresponds to the cutoff of the NP-NP potential and represents a configuration in which two NPs are touching; we therefore call it contact peak. Eventually, the contact peak becomes higher than the secondary peak, an evidence of the formation of large polymer-free regions (Fig. 2d).

The structure factor is defined as [62]

<figure><img src="image_5.png" /><figcaption>FIG. 4. (color online) NP-NP structure factor at different NP volume fractions $\phi_N$ for $\sigma_N = 2$ (a) and $\sigma_N = 4$ (b).</figcaption></figure>

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

where $\mathbf{q}$ is the wavevector. Since our configurations are isotropic, also in this case we will consider the spherically averaged structure factor $S(q)$ [61].

Because $S(q) - 1$ is the Fourier transform of $\rho g(r)$ [62], we can in principle find in the NP-NP structure factor $S_{NN}(q)$ the same information that we find in $g_{NN}(r)$. If the position of the main peak of $g(r)$ is $r = r_0$, the main peak of $S(q)$ will be at $q_0 \approx 2\pi/r_0$, although the precise value of $q_0$ depends on temperature and density [63]. Hence, we expect to find the main peak of $S_{NN}(q)$ at $q \approx 2\pi/r_s$ at low NP volume fraction and at $q \approx 2\pi/r_c$ at high NP volume fraction, as we indeed observe in Figs. 4a-b. We also notice that, while in the $g(r)$ we can clearly distinguish two peaks at intermediate values of $\phi_N$ (Fig. 3), in the $S(q)$ their contributions interfere with each other and result in a single peak that is shifted towards higher wavevectors as $\phi_N$ is increased. Therefore, interpretation of $S_{NN}(q)$ might not always be straightforward.

In what follows, we will mainly consider those systems in which the NPs are well dispersed in the polymer so-
lution (Fig. 2a-b). As a qualitative criterion, we define a system with good NP dispersion as one where the secondary peak of $g_{NN}(r)$ is higher than or comparable to the contact peak. It should be noticed that the maximum volume fraction that we can reach while keeping a good NP dispersion depends on the NP diameter $\sigma_N$. To see this, we consider the interparticle distance $h$ [19–22], which represents the average spacing between the surfaces of neighboring nanoparticles. In the literature, the following expression for $h$ has often been used [19–22]:

$$h^{\text{th.}} = \sigma_N \left[ \left( \frac{\phi_N^M}{\phi_N} \right)^{1/3} - 1 \right], \quad (7)$$

where $\phi_N^M$ represents the maximum NP volume fraction, at which $h = 0$. The value $\phi_N^M = 2/\pi \approx 0.637$, corresponding to the (ill-defined [64]) random close packing, is often employed [19–22]. However, this definition presents some issues (Sec. AI in the Appendix), and therefore we have chosen to measure $h$ directly from the data using the pore-size distribution [65]. This approach is similar to the one used by Li et al. in Ref. 30, with the difference that they used an Euclidean distance map. For the details on how $h$ can be extracted from the pore-size distribution, see Sec. AI in the Appendix.

In Fig. 5, we show the interparticle distance calculated from the pore-size distribution, $h$, versus the NP volume fraction: filled (open) symbols represent systems with a good (poor) NP dispersion (according to the above defined criterion). We also report for comparison the “theoretical” interparticle distance $h^{\text{th.}}$, Eq. (7), with $\phi_N^M = 0.637$ (continuous lines); we note that the two quantities are very similar, with $h^{\text{th.}}$ being on average slightly larger than $h$. As we can see, NP dispersion starts to become poor when $h \approx 1$, i.e., when the average distance between the surface of neighboring NPs becomes comparable with the monomer size, in qualitative agreement with the snapshots shown in Fig. 2c-d.

## B. Polymers

The radial distribution function of the monomers, $g_{mm}(r)$, is strongly dominated by the short-distance signals coming from the chain bonds (see Fig. S2 in S.I.), and therefore it is not easy to extract from it information about the medium and long range distribution of the monomers. Hence we focus our attention on the monomer-monomer structure factor, $S_{mm}(q)$, shown in Fig. 6.

Figure 6a, shows $S_{mm}(q)$ for $\sigma_N = 4$. At $\phi_N = 0$ (pure polymer solution), there is a small peak at $q^* \approx 1.4$ (inset of Fig. 6a), which in real space corresponds to a distance $r^* = 2\pi/q^* \approx 4.5$. This peak reveals the presence of a typical length scale in the NP-free system, which can be interpreted as the average size of the holes in the polymer matrix [66]. The main peak of $S_{mm}(q)$

<figure><img src="image_6.png" /><figcaption>FIG. 5. Interparticle distance as a function of the NP volume fraction $\phi_N$. Filled symbols: systems with good NP dispersion. Open symbols: systems with poor NP dispersion. Continuous lines: “theoretical” interparticle distance $h^{\text{th.}}$, calculated from Eq. (7), with $\phi_N^M = 0.637$.</figcaption></figure>

is at $q_0 = 7.8 \approx 2\pi/r_b$, where $r_b = 0.96$ is the average monomer-monomer bond length.

For $\phi_N > 0$, the spatial arrangement of the NPs starts to be visible as a modulation in $S_{mm}(q)$, with a main peak appearing approximately at the same wavevector as the main peak of $S_{NN}(q)$, as we can see from Fig. 6b, where $S_{mm}(q)$ is compared to $S_{NN}(q)$. At even higher NP volume fraction, a signal starts to appear at $q = 0$, due to the fact that the polymers are getting far from each other (see Fig. 2d). If $\phi_N$ was increased even more, eventually the monomer volume fraction would become smaller than the overlap volume fraction (dilute regime) and $S_{mm}(0)$ would saturate to $N$ [33].

Another important quantity to characterize the structure of the polymer mesh is the correlation length or mesh size $\xi$, which for the pure polymer solution ($\phi_N = 0$) can be estimated via scaling considerations [33, 34]:

$$\xi = \begin{cases} (R_g^*/\sqrt{3}) & \rho_m < \rho_m^* \\ (R_g^*/\sqrt{3}) (\rho_m/\rho_m^*)^{-\nu/(3\nu-1)} & \rho_m > \rho_m^* \end{cases} \quad (8)$$

where $R_g^*$ is the radius of gyration of an isolated chain, $\rho_m^*$ is the overlap monomer concentration and $\nu = 0.588$ is the Flory exponent [33] [67].

In this work, $R_g^* = 7.48$ and consequently, defining $\rho_m^* = N(4\pi R_g^*/3)^{-1}$ (other definitions are possible [34]), we get $\rho_m^* = 5.70 \cdot 10^{-2}$. Using these values we obtain from Eq. (8) $\xi = 1.27$ for the pure polymer solution ($\rho_m = 0.280$).

The structure of the individual polymer chains can be characterized by the function $p(r) = 4\pi r^2 \rho g_1(r)/(N-1)$, where $\rho g_1(r)$ is obtained by applying Eq. (5) to a single polymer chain (and, as usual, taking the spherical average). The quantity $p(r)dr$ represents the probability to find a monomer
<figure><img src="image_7.png" /><figcaption>FIG. 6. (color online) (a) Monomer-monomer structure factor $S_{mm}(q)$ (log-lin scale) at different NP volume fractions for $\sigma_N = 4$. Inset: $S_{mm}(q)$ in the pure polymer solution (linear scale) (b) Comparison between the structure factors of monomers and NPs for $\sigma_N = 4$. For clarity, the curves representing $S_{mm}(q)$ have been shifted up by a factor of 3.</figcaption></figure>

belonging to the same chain at distance between $r$ and $r + dr$ from a given monomer.

For a Gaussian chain, $p$ has the following expression [34]:

$$p_G(x) = \frac{8Nx}{\sqrt{\pi}(N-1)R_g} \left[ \frac{\sqrt{\pi}}{2}(1+2x^2)\text{erfc}(x) - xe^{-x^2} \right], \quad (9)$$

where $x = r/2R_g$ and erfc($x$) is the complementary error function. This probability density peaks at $r \approx 0.74R_g$ [34].

In Fig. 7a, we show $p(r)$ for different values of $\sigma_N$ and $\phi_N$, along with $p_G$ for a Gaussian chain, Eq. (9), with $R_{g0} = R_g(\phi_N = 0)$ (dashed line). We observe that, for small $\phi_N$, $p_G$ provides a good approximation of $p$ at intermediate and large $r$ (at small $r$, $p(r)$ is dominated by excluded volume interaction between nearest neighbors). For all values of $\phi_N$, $p(r)$ shows a very high peak at $r \approx r_b = 0.96$, corresponding to the first nearest neighbor, and a smaller peak at $r \approx 2r_b = 1.92$, corresponding to the second nearest neighbor. For values of

<figure><img src="image_8.png" /><figcaption>FIG. 7. (color online) (a) The function $p(r)$ for different values of the NP diameter $\sigma_N$ and volume fraction $\phi_N$. The long black arrow indicates the direction of increasing $\phi_N$. The small arrows indicate the “bump” at $r \approx \sigma_N + 1$. Dashed lines: Gaussian chain approximation $p(x) = p_G(x)$ with $R_g = R_g(\phi_N = 0)$. For the sake of clarity, for $\sigma_N \le 5$ every set of curves has been shifted upwards with respect to the previous set by 0.035. (b) The ratio $p(r)/p_0(r)$ for $\sigma_N = 7$ and different values of $\phi_N$. The small arrows help locating the two “bumps” at $r \approx \sigma_N + 1$ and $2(\sigma_N + 1)$. Cartoon: schematic representation of the polymer structure in the presence of the NPs.</figcaption></figure>

$r$ larger than $2r_b$, this signal gets washed out, and ultimately $p(r)$ decays to zero. When $\phi_N$ increases, we observe two effects: The third nearest neighbor peak becomes more pronounced and the curve becomes broader. This indicates that the presence of the NPs stretches the chains, causing them to become locally more ordered. We also note that $p(r)$ shows a modulation of wavelength $\approx \sigma_N + 1$, the first peak of which is clearly visible in Fig. 7a as a “bump” at $r \approx \sigma_N + 1$ (colored arrows). The presence of this modulation can be better appreciated by plotting the ratio $p(r)/p_0(r)$, where $p_0(r) = p(\phi_N = 0, r)$. In Fig. 7b, we report $p(r)/p_0(r)$ for $\sigma_N = 7$: As we can see, the effect of the NPs is to produce a “hole” in the range $0 \lesssim r \lesssim \sigma_N + 1$, but also to stretch the chain, increasing $p(r)$ significantly at larger distances. The modulation is clearly visible, with two bumps appearing at
<figure><img src="image_9.png" /><figcaption>FIG. 8. Reduced radius of gyration of the polymers as a function of the NP volume fraction $\phi_N$ (a) and of interparticle distance $h$ (b). Inset: $R_g/R_{g0}$ as a function of the inverse interparticle distance $1/h$. Continuous lines: $R_g = R_{g0}(1 + 0.0762 \cdot h^{-2.26})$.</figcaption></figure>

$r \approx \sigma_N + 1$ and $2(\sigma_N + 1)$ (small arrows).

Chain swelling in the presence of NPs has already been predicted theoretically [51] and observed in both simulations [52] and experiments [68–70]. In particular, Karantzos et al. [52] have shown that polymer chains are unperturbed by the presence of repulsive NPs, while attractive NPs cause the polymer chains to be stretched and flattened when $2R_g > \sigma_N$ (which is always the case for the systems that we considered). Using the SC/PRISM theory, Frischknecht et al. [51] reached the same conclusions. For a recent review discussing the influence of NPs on polymer size and local structure in simulations, see Ref. 71.

In order to quantify the expansion of the chains, we measure the radius of gyration $R_g$, defined as [33]

$$R_g^2 = \frac{1}{N} \sum_{i=1}^{N} \langle (\mathbf{r}_i - \mathbf{R}_{\text{CM}})^2 \rangle, \quad (10)$$

where $\mathbf{R}_{\text{CM}}$ is the position of the center of mass of the polymer. In the pure polymer solution, we have $R_{g0} = R_g(\phi_N = 0) = 6.28 \pm 0.02$. In Fig. 8a we present

the reduced radius of gyration $R_g/R_{g0}$ as a function of NP volume fraction for different values of $\sigma_N$. With the exception of $\sigma_N = 1$, there is a modest but clear increase of $R_g/R_{g0}$ with increasing NP volume fraction. We also notice that at fixed NP volume fraction the increase is stronger for smaller NPs, which suggests that chain expansion is mainly controlled by the NP excluded volume, which, at fixed $\phi_N$, is larger for smaller NPs [72].

Fig. 8a also shows that for $\sigma_N = 1$, $R_g$ decreases with increasing NP volume fraction. The reason is that NPs of this size have the largest surface-to-volume ratio, making the monomer-NP interaction (which scales approximately with the NP surface) very relevant. The consequence is that while in this range of $\phi_N$ the effect of the excluded volume is small, the effect of the interaction is large: Small NPs produce an effective attractive interaction between the monomers, which results in a decrease of $R_g$ and of the overall monomer volume fraction $\phi_m$ (we recall that all the simulations were performed at the same average pressure $P = 0.1$; see also Sec. V). We can therefore say that in this range of $\phi_N$, the NPs of size $\sigma_N = 1$ act like a poor solvent, promoting chain contraction.

In Fig. 8b we plot the reduced radius of gyration as a function of the interparticle distance $h$. For $\sigma_N \ge 3$, the data fall on a master curve, which can be approximated by the empirical expression $R_g = R_{g0}(1 + 0.0762 \cdot h^{-2.26})$ (continuous line in Fig. 8b). This confirms that in this range of NP size chain expansion is a geometrical effect, dominated by excluded volume: the NPs force the chains to take less tortuous paths, therefore increasing their effective size. The fact that larger particles have a locally “flatter” surface that could enhance chain expansion does not seem to play a role in this size range, as we can conclude from the fact that data for different $\sigma_N$ fall on the same master curve. For $\sigma_N = 1$ and 2 the data do not fall on the master curve, for the reasons explained above (high surface-to-volume ratio promotes chain contraction). To provide a better resolution for small values of $h$, in the inset of Fig. 8b we plot $R_g/R_{g0}$ as a function of $1/h$ [51].

We can summarize our results by saying that NPs of size $\sigma_N \ge 2$ act like a good solvent, swelling the polymer chains, while NPs of size $\sigma_N = 1$ act like a poor solvent, causing them to contract. We note that this effect is expected to depend on the strength of the monomer-NP interaction: With stronger interactions, chain contraction could be observed also for $\sigma_N > 1$. Further study is needed in order to clarify this point.

## IV. DYNAMICS

### A. Mean squared displacement

To characterize the dynamics of the system, we study the mean squared displacement (MSD) of the NPs and of the centers of mass (CM) of the chains. We recall that
<figure><img src="image_10.png" /><figcaption>FIG. 9. (color online) (a) MSD of the NPs divided by time. Dashed line: $\langle r^2(t) \rangle / t \propto t^{-1/2}$, predicted in Ref. 11 for the subdiffusive regime of intermediate size NPs. Dash-dotted line: $\langle r^2(t) \rangle / t \propto t$ (ballistic regime). The curves for $\sigma_N = 7$ have been shifted downwards by a factor 2 to facilitate the visualization. (b) MSD of the centers of mass of the chains divided by time. Dashed line: $\langle r^2(t) \rangle / t \propto t^{-0.13}$ (subdiffusive transient). The curves for $\sigma_N = 4$ have been shifted down by a factor 2 and those for $\sigma_N = 7$ by a factor 8 for clarity.</figcaption></figure>

the MSD of a system of M particles is defined as [62]:

$$ \langle r^2(t) \rangle = \frac{1}{M} \sum_{i=1}^{M} \langle [\mathbf{r}_i(t) - \mathbf{r}_i(0)]^2 \rangle. \quad (11) $$

In order to visualize more clearly the transition between the short-time ballistic regime, $\langle r^2(t) \rangle \propto t^2$, and the long time diffusive regime, $\langle r^2(t) \rangle \propto t$, we show in Fig. 9 the MSD divided by time $t$ for the NPs (Fig. 9a) and for the polymers (Fig. 9b). At low $\phi_N$, the motion of the NPs shows the same qualitative behavior for all the values of the NP diameter $\sigma_N$ (Fig. 9a): After the initial ballistic regime, the motion becomes almost immediately diffusive, with the exception of the system with $\sigma_N = 7$, which shows a weak subdiffusive transient, $\langle r^2(t) \rangle \propto t^\beta$ ($\beta < 1$), between these two regimes. A clear transient subdiffusive regime appears between the ballistic and diffusive regimes at intermediate and high values of $\phi_N$. The MSD

of the chains, on the other hand, shows a weak subdiffusive transient for all values of $\phi_N$ and $\sigma_N$, with an exponent $\beta \approx 0.87$ that is not much influenced by the value of $\phi_N$ (Fig. 9b). This transient, which is most likely due to non-Gaussian dynamics caused by intermolecular correlations [73], has been previously observed in experiments [73–75] and simulations [48, 73, 74, 76, 77] of polymer melts, where the measured subdiffusive exponent was $\beta \approx 0.8$. The fact that in our case the exponent is slightly larger than 0.8 is likely due to the fact that the density considered here is significantly smaller than that of a melt, $\rho \approx 0.85$. The different regimes (ballistic, subdiffusive, diffusive) and the transitions between them can also be studied systematically through the function $d \log \langle r^2(t) \rangle / d \log t$, which is a generalization of the subdiffusive exponent $\beta$. For a detailed analysis of this quantity, see Sec. SIII in the S.I.

Using scaling arguments, Cai et al. have formulated a theory for the diffusion of single nonsticky NPs in polymer liquids [11]. For an unentangled polymer mixture, they predicted that NPs of diameter $\sigma_N < \xi$, where $\xi$ is the mesh size, should always move diffusively, whereas the MSD of larger NPs, $\xi < \sigma_N < 2R_g$, should behave as follows:

$$ \langle r_N^2(t) \rangle \propto \begin{cases} t & t < \tau_\xi \\ t^{1/2} & \tau_\xi < t < \tau_{\sigma_N} \\ t & \tau_{\sigma_N} < t, \end{cases} \quad (12) $$

where $\tau_\xi \propto \eta_s \xi^3 / k_B T$ and $\tau_{\sigma_N} \propto \tau_\xi (\sigma_N / \xi)^4$ are, respectively, the relaxation times of polymer segments of size $\xi$ and $\sigma_N$ (here $\eta_s$ is the viscosity of the pure solvent). For NPs larger than the polymers, $\sigma_N > 2R_g$, $\tau_{\sigma_N}$ must be replaced by $\tau_{R_g} \propto \tau_\xi (R_g / \xi)^4$. The crossover from subdiffusive to diffusive motion for NPs in polymer solutions has also been observed in experiments [24, 25].

Since at low NP volume fraction in our system $\xi \approx 1.3$, according to the scaling theory of Cai et al. one expects the MSD of the NPs of diameter $\sigma_N > 2$ to show subdiffusive behavior with exponent $\beta = 1/2$ at small $\phi_N$. However, no such behavior is observed for any value of $\sigma_N$. For small NPs, this may be due to the fact that the time window in which the subdiffusive behavior is expected to be present, i.e., $\tau_\xi < t < \tau_{\sigma_N}$, is too small, since $\tau_{\sigma_N} \propto \tau_\xi (\sigma_N / \xi)^4$. For larger NPs, this time window regime should be large enough to observe subdiffusion, and indeed for $\sigma_N = 7$ we observe a very weak subdiffusive transient, but the exponent $\beta$ is close to 1. This is in agreement with previous simulations, which have also found that $\beta$ is not always equal to 1/2 in the subdiffusive regime, but rather gradually approaches this value as $\sigma_N$ is increased [78, 79].

In addition to the mean squared displacement, we have also studied the van Hove function [62] and the non-Gaussian parameter $\alpha_2(t)$ [80] of the NPs, finding that their dynamics is with a good approximation Gaussian, in agreement with experiments [24, 25] (see Sec. SIV in
the S.I.). These results indicate that the dynamics of the NPs is not heterogeneous.

## B. Polymer diffusion

In order to make a more quantitative characterization of the dynamical properties of the polymers and the NPs, we now focus on the self diffusion coefficient $D$ (which for simplicity we will refer to as “diffusion coefficient”), which can be obtained from the MSD, Eq. (11), through Einstein’s relation [62]:

$$D = \lim_{t \to \infty} \frac{\langle r^2(t) \rangle}{6t}. \qquad (13)$$

It is known that measurements of $D$ in systems with periodic boundary conditions suffer from finite-size effects because of long-ranged hydrodynamic interactions [81, 82]. Although an analytical expression for the correction to $D$ is available [81, 82], it is not evident whether it can be applied to the motion of polymer chains and NPs in a concentrated polymer solution. For the NPs, such an expression is most likely not adequate when, as in our case, the NP size is smaller than the polymer size [28]. Therefore, for consistency we choose not to apply any finite size correction to the measured diffusion coefficients.

In the pure polymer system ($\phi_N = 0$), the diffusion coefficient of the CM of the chains is $D_{p0} = (1.14 \pm 0.02) \cdot 10^{-2}$. In Fig. 10a we plot the reduced diffusion coefficient of the polymer chains $D_p/D_{p0}$ as a function of the NP volume fraction $\phi_N$. We can observe that $D_p/D_{p0}$ decreases with increasing NP volume fraction, with the decrease being stronger, at fixed $\phi_N$, for smaller NPs. The data can be fitted to the empirical functional form $D_p = D_{p0}[1 - (\phi_N/\phi_{N0})^\alpha]$, where $\alpha$ increases with NP size (the values of $\phi_{N0}$ and $\alpha$ for the different NP diameters are reported in Tab. S1 in the S.I.). We note that this functional form implies that $D_p$ becomes zero at $\phi_N = \phi_{N0}$, i.e., that the dynamics of the polymers is completely arrested. However, for larger values of $\phi_N$ the dependence of $D_p$ on $\phi_N$ changes (see Sec. V) and thus this dynamic transition is avoided; hence, the above functional form is valid only if $\phi_N/\phi_{N0}$ is small. By using this relation, we can interpolate between the data points and plot $D_p/D_{p0}$ as a function of the NP diameter $\sigma_N$ for different volume fractions (Fig. 10b) and we observe that $D_p/D_{p0}$ increases monotonically with $\sigma_N$ at fixed $\phi_N$.

There are two possible causes (or a combination of the two) that can lead to the slowing down of the chains with increasing NP volume fraction: the increase of the number of obstacles to polymer motion and the increase of polymer-NP interfacial area, which, since the interaction between polymers and NPs is attractive, can result in a reduced chain mobility. A predominance of the first effect would imply that the slowing down of the polymers

<figure><img src="image_11.png" /><figcaption>FIG. 10. (a) Reduced diffusion coefficient of the centers of mass (CM) of the chains $D_p/D_{p0}$ as a function of NP volume fraction, where $D_{p0} = D_p(\phi_N = 0)$. Continuous lines: empirical fits with the functional form $D_p = D_{p0}[1 - (\phi_N/\phi_{N0})^\alpha]$. (b) $D_p/D_{p0}$ as a function of NP diameter $\sigma_N$.</figcaption></figure>

is a mostly entropic effect, while if the second effect is the most important the dynamics of the polymers is mainly controlled by enthalpy.

Composto and coworkers [19–22] observed a similar slowing down of chain motion in a series of experimental studies on polymer nanocomposites containing large NPs ($\sigma_N \gtrsim 2R_g$). They found that the reduced diffusion coefficient of the polymers falls on a master curve when plotted versus a “confinement parameter”, defined as $h/2R_g$, where $h$ is the interparticle distance, which the authors computed using Eq. (7) with $\phi_N^M = 2/\pi$. Since the collapse of the data was independent of the strength of the polymer-NP interaction [21], the authors concluded that the slowing down of the polymers is entropic in origin, caused by the reduction of chain entropy as the chain passes through bottlenecks formed by neighboring NPs (entropic barrier model) [19]. An analogous reduction in polymer mobility due to the presence of NPs was observed by Li et al. [30] in molecular dynamics simulations of unentangled melt of short chains ($N = 35, \rho_m = 0.85$) containing repulsive NPs. The slowing down was weaker than that observed by Composto
<figure><img src="image_12.png" /><figcaption>FIG. 11. Reduced diffusion coefficient of the centers of mass (CM) of the chains as a function of the dynamic confinement parameter $h/\lambda_d$ for different values of the NP diameter $\sigma_N$. Continuous line: $D_p = D_{p0}[1 - \exp(-h/\lambda_d)]$. Inset: $D_p/D_{p0}$ as a function of $h/2R_g$. Continuous and dash-dotted lines: $D_p = D_{p0}[1 - \exp(-a \cdot h/2R_g)]$, with a respectively equal to 5.44 and to 3.22.</figcaption></figure>

and coworkers in Refs. 19, 20, an effect which the authors attributed to the absence of chain entanglements. Karatrantos et al. [32] also observed a monotonic decrease in the polymer diffusion coefficient with increasing NP volume fraction in molecular dynamics simulations of NPs in unentangled and weakly entangled melts, and attributed this phenomenon to the increase in the polymer-NP interfacial area. Desai et al. [39], on the other hand, have reported that the polymer diffusion coefficient in a simulated lightly entangled melt ($N = 80, \rho_m = 0.85$) containing repulsive/weakly attractive NPs initially increases with $\phi_N$, reaches a maximum around $\phi_N = 5\%$ and decreases for higher values. An enhancement of chain diffusivity at low $\phi_N$ has also been observed in simulations by Kalathi et al. [28], possibly because attractive monomer-monomer interactions were considered in their work. It is therefore clear that, despite the fact that some general trends can be identified, the dynamics of the polymers can depend strongly on the details of the simulated systems.

Following Composto and coworkers [19–22], we plot in the inset of Fig. 11 the reduced diffusion coefficient of the chains as a function of the confinement parameter $h/2R_g$. We recall that in our case $h$ is not defined by Eq. (7), but rather computed from the pore-size distribution (see Sec. AI in the Appendix). In addition to the $T = 1.0$ data, we also show the results from simulations at $T = 0.7$ with $\sigma_N = 3$.

The first observation is that, with the exception of $\sigma_N = 1$, all the $T = 1.0$ data fall on the same master curve, which is well approximated by the empirical expression $D_p = D_{p0}[1 - \exp(-ah/2R_g)]$, with $a = 5.44$ (continuous line in the inset of Fig. 11). The fact that also in our case $D_p/D_{p0}$ is only a function of the confinement parameter is rather surprising, since Composto

and coworkers mainly considered NPs of size comparable to that of the polymers or larger, which could be considered as basically immobile [19–22], whereas in our case $\sigma_N < 2R_g$ and the NPs diffuse faster than the chains in almost all the systems considered (see Fig. S7 in the S.I.).

We notice, however, two important differences: The first one is that while in our case the diffusion coefficient of the pure polymer solution ($D_p/D_{p0} = 1$) is recovered at $h/2R_g \approx 1$, in Refs. 19–22 it is recovered only at much higher values of the confinement parameter, $h/2R_g \approx 20$. Our finding is similar to what observed by Li et al. [30], who attributed the discrepancy between their data and those of Composto and coworkers to the absence of entanglement in their simulated system. The second difference is that the $T = 0.7$ data clearly do not fall on the same master curve. Since a decrease in temperature is approximately equivalent to an increase in the strength of the polymer-NP interaction, this result suggests that in our system the polymer-NP interaction plays a relevant role, in contrast with Ref. 21, where the authors concluded that the confinement parameter captures the polymer slowing down independently of the polymer-NP interactions. We propose in the following a possible solution to these discrepancies.

The confinement parameter $h/2R_g$ is a purely static quantity, which only depends on the spatial configuration of the polymers and the NPs in the system. However, there are several cases in condensed matter physics in which two systems with identical structure show a completely different dynamics: A well-known example is that of the glass transition, where a supercooled liquid shows structural properties identical to those of a liquid at higher temperature, but completely different dynamical properties [83–85].

It seems therefore more appropriate to introduce a dynamic confinement parameter $h/\lambda_d$, where $\lambda_d$ is a dynamic length scale which will in general depend on temperature, density and on the details of the simulated system. We have already seen that the $T = 1.0$ data are well approximated by the function $D_p = D_{p0}[1 - \exp(-ah/2R_g)]$, with $a = 5.44$ (continuous line in the inset of Fig. 11); the $T = 0.7$ data are well approximated by the same functional form, but with a different coefficient, $a = 3.22$. In light of what we discussed above, we make the hypothesis that the reduced diffusion coefficient of the polymers can be expressed as

$$D_p = D_{p0} \left[ 1 - \exp \left( - \frac{h}{\lambda_d(T)} \right) \right], \quad (14)$$

where we have explicitly reported the dependence of $\lambda_d$ on temperature. Since $R_g$ does not change more than 8% with respect to the pure polymer solution value $R_{g0}$ (see Fig. 8), we can estimate $\lambda_d$ as $\lambda_d = 2R_{g0}/a$: This gives $\lambda_d(1.0) = 2.31$ and $\lambda_d(0.7) = 3.90$. We show $D_p/D_{p0}$ as a function of $h/\lambda_d$ in Fig. 11: In this plot, the data for different temperatures fall on the same master curve, showing that the dynamic confinement parameter $h/\lambda_d$
is more successful than the “static” confinement parameter $h/2R_g$ in capturing the slowing down of the polymers (inset of Fig. 11). However, one question remains: What does the dynamic length scale $\lambda_d$ represent exactly, and why does it increase when temperature is decreased? Our answer is that $\lambda_d$ is a cooperativity length scale, i.e., it represents the typical length scale of the spatial rearrangement needed for a polymer segment to escape its local cage. Similarly to what happens in a supercooled liquid [83, 86], this cooperative length scale is expected to increase when $T$ is decreased. In our system, an important role could also be played by the attractive polymer-NP interactions, which become more relevant when $T$ is decreased and could reduce the mobility of polymer segments close to the polymer-NP interface. We also expect $\lambda_d$ to increase with monomer density, since a higher density naturally leads to a locally more constrained dynamics: This could explain why the data of Li et al. [30], who simulate NPs in a dense melt, are compatible with a larger cooperativity length scale. Another factor that is expected to play a major role is the stiffness of the chain, with stiffer chains expected to lead to a larger $\lambda_d$.

To sum up, we propose a modification of the confinement parameter theory of Composto and coworkers [19–22]: our hypothesis is that the dynamics of the polymers is controlled by a dynamic confinement parameter $h/\lambda_d$, where $\lambda_d$ is a cooperativity length scale which will depend in general on the thermodynamic parameters and on the details of the model. Further study is required to test the validity of this hypothesis, and to understand how $\lambda_d$ depends on the properties of the physical system.

### C. Single nanoparticle diffusion

The diffusion coefficient of a hard-sphere probe particle of diameter $\sigma_N$ in a continuum solvent with shear viscosity $\eta$ is given by the Stokes-Einstein equation [62]:

$$D_{N0} = \frac{k_B T}{f \pi \eta \sigma_N}, \qquad (15)$$

where $f$ is a number between 2 and 3 which depends on the choice of the hydrodynamic boundary conditions: $f=2$ for pure slip and $f=3$ for pure stick boundary conditions [87]. If the particle is not a perfect hard sphere, for example because its shape is not perfectly spherical or because there is adsorption of solvent molecules on its surface, $\sigma_N$ must be replaced with an effective hydrodynamic diameter $\sigma_h$ [88, 89]. It is well-known that Eq. (15) is inadequate to describe the motion of particles smaller than the polymer size in a polymer solution/melt [11, 16–18, 24, 25, 27, 28, 37, 38], because the continuum assumption breaks down when the size of the probe particle becomes comparable to the characteristic length scale of the solvent.

Cai et al. [11] have predicted three regimes for the diffusion of a NP of diameter $\sigma_N$ in an unentangled poly-

<figure><img src="image_13.png" /><figcaption>FIG. 12. Single NP diffusion coefficient multiplied by the effective hydrodynamic diameter of the NP $\sigma_h$ as a function of $\sigma_h/2R_g$ (blue circles), compared with the results from Ref. 27 (red diamonds). Filled symbols: $\sigma_h = \sigma_N$. Open symbols: $\sigma_h = \sigma_N + 1$. Continuous lines: slope of -2. The dashed vertical line marks the value $\xi/2R_g$ for our system.</figcaption></figure>

mer mixture: For $\sigma_N < \xi$, where $\xi$ is the mesh size, the NP diffusion coefficient should follow the Stokes-Einstein law: $D_{N0} \propto k_B T/\eta_s \sigma_N$, with $\eta_s$ the viscosity of the pure solvent (small size regime). If $\sigma_N > \xi$, the motion of the NPs becomes coupled to the segmental relaxation of the polymer mesh, so that $D_{N0} \propto k_B T/\eta_s \sigma_N^3$ (intermediate size regime). The relation $D_{N0} \propto \sigma_N^{-3/2}$ for intermediate size NPs was originally proposed by Wyart and de Gennes using scaling arguments [37], and was also predicted by Yamamoto and Schweizer using mode-coupling theory [13] and subsequently a self-consistent generalized Langevin equation approach [14]. This prediction has also been confirmed by simulations [27, 79].

For even larger diameters, the Stokes-Einstein law (15) is eventually recovered, but with the viscosity of the pure solvent replaced by the bulk viscosity of the solution $\eta$: $D_{N0} \propto k_B T/\eta \sigma_N$ (large size regime). For this second crossover size, different versions of mode-coupling theory predict different values depending on the level of approximation and on polymer density, all of them of the order of the polymer diameter: $\sigma_N = 2R_g$ (unentangled solutions and melts) [12], $\sigma_N = 3R_g$ (unentangled melts and unentangled concentrated solutions) [13] and $\sigma_N = 2R_g^*$ (semidilute solutions) [15], where $R_g^*$ is the radius of gyration of an isolated chain. The prediction that the crossover to Stokes-Einstein behavior should occur when the NP diameter is of the order of $2R_g$ has also been confirmed by experiments [17, 18, 24].

To test the validity of the Stokes-Einstein formula, Liu et al. [27] have used simulations to measure the single particle diffusion coefficient of NPs in simulations of a dense, unentangled melt ($N=60, \rho_m = 0.84$). The results of their simulations are shown in Fig. 12 (red diamonds). The authors argued that the effective hydrodynamic radius of the particle, $R_h = \sigma_h/2$, should have the
value $R_h = (\sigma_N + 1)/2$, which corresponds to the contact
distance between a NP and a monomer (the same argu-
ment can be found in Ref. 90). By fitting their data in
the size range $\sigma_N < 2R_g$ with a power law $D_{N0} \propto R_h^{-\gamma}$,
they found $\gamma \approx 3$ (open red diamonds in Fig. 12), and for
diameters $\sigma_N > 2R_g$ they recovered the Stokes-Einstein
relation. The results of Ref. 27 are therefore in agree-
ment with the prediction that $D_{N0} \propto \sigma_N^{-3}$ [11, 13, 37]
if one replaces the NP diameter $\sigma_N$ with the effective
hydrodynamic diameter $\sigma_h = \sigma_N + 1$. However, when
plotting $D_{N0}$ as a function of $\sigma_N$ instead of $R_h$ and fit-
ting with a power law $D_{N0} \propto \sigma_N^{-\gamma}$, one obtains instead
$\gamma \approx 2$ (filled red diamonds in Fig. 12). Hence one must
conclude that the value of the exponent $\gamma$ depends on
the exact definition of the effective NP diameter, which
makes the comparison of simulation data with theoretical
predictions a delicate matter, especially when the size of
the NP is of the same order of magnitude as the monomer
size (for large NPs, $\sigma_N + 1 \approx \sigma_N$).

In order to test these predictions, we have performed
additional simulations at low NP volume fraction ($\phi_N <$
0.015) for $\sigma_N = 10, 12$, and 14. In Fig. 12 we show
$D_{N0}\sigma_h$ as a function of $\sigma_h/2R_g$, with $\sigma_h$ alternatively
defined as $\sigma_N$ and $\sigma_N + 1$ (blue circles). Also included
are the data from Ref. 27 (red diamonds). We can see
that $D_{N0}$ decreases continuously for $\sigma_h < 2R_g$, whereas
at $\sigma_h \approx 2R_g$ Stokes-Einstein behavior, $D_{N0}\sigma_h = \text{const.}$,
is recovered. Taking $\sigma_h = \sigma_N + 1$ we find in the range
$0.6 \lesssim \sigma_h/2R_g \lesssim 1$ a slope of approximately $-2$, which
agrees with the theoretical predictions. However, the
range in which we observe this slope is rather small and
hence we cannot claim that our data confirm the theory.
In particular, if we use the scaling estimate of Eq. (8) for
the mesh size, $\xi \approx 1.27$, we can see that there is no sign of
the transition from $D_{N0}\sigma_h = \text{const.}$ to $D_{N0}\sigma_h \propto \sigma_h^{-2}$ at
$\sigma_h \approx \xi$ (which corresponds to $\sigma_h/2R_g \approx 0.10$) predicted
in Ref. 11. However, a caveat is in order: We have veri-
fied that, as also reported in previous studies [27, 38, 39],
the diffusion coefficient of small NPs decreases when the
NP mass increases at fixed NP volume, i.e., when the
mass density is increased. The effect becomes progres-
sively weaker as $\sigma_N$ is increased, and at $\sigma_N = 7$ no mass
density dependence is observed. Nevertheless, this effect
should be taken into account when comparing the results
of simulations to those of experiments or to theoretical
predictions. For a detailed discussion, see Sec. AII in the
Appendix.

D. Nanoparticle diffusion

In the previous section, we have dealt with the mo-
tion of a single NP in the polymer solution, i.e., we have
considered the dilute NP limit: We will now discuss the
dynamics of NPs at higher NP volume fraction $\phi_N$.

Only few simulation studies have considered high NP
volume fractions. Liu et al. [27] have observed a reduc-
tion of the NP diffusion coefficient with increasing $\phi_N$,

<figure><img src="image_14.png" /><figcaption>FIG. 13. Reduced NP diffusion coefficient as a function of the NP volume fraction $\phi_N$ (a) and of NP diameter $\sigma_N$ (b).</figcaption></figure>

and attributed the phenomenon to polymer-mediated in-
teractions, i.e., to the formation of chain bridges be-
tween neighboring NPs that would hinder NP motion;
it is not clear, however, what the typical lifetime of such
bridges should be, and thus whether this explanation is
correct. Karatrantos et al. [32] have observed a similar
reduction in NP mobility and argued that it is due “to
both nanoparticle-polymer surface area and nanoparticle
volume fraction” [32], implying that pure geometry and
polymer-NP attraction both play a role. The importance
of polymer-NP interaction in NP dynamics is beyond
dispute: Patti [29] showed that the diffusion coefficient
of NPs in an unentangled melt decreases monotonically
when the strength of the polymer-NP interaction is in-
creased, with the decrease being stronger for smaller NPs.
A monotonic decrease of NP diffusivity with the strength
of the polymer-NP interaction was also observed by Liu
et al. [27]. We mention, however, that this trend can
be reversed (NP diffusivity *increasing* with increasing in-
teraction strength) in strongly entangled systems, where
the dynamics of the NPs is dominated by density fluctu-
ations on length scales of the order of the tube diameter
[13].

In Fig. 13a, we show the reduced diffusion coefficient
of the NPs, $D_N/D_{N0}$, where $D_{N0}$ is the diffusion coeffi-
cient of a single NP in the polymer solution, as a function of the NP volume fraction $\phi_N$; also shown are data for $\sigma_N = 3$ and $T = 0.7$. Similarly to the diffusion coefficient of the chains, $D_N$ decreases with increasing NP volume fraction. The first thing that one can notice is that the decrease of $D_N/D_{N0}$ with the NP volume fraction is rather quick: Already at the modest volume fraction of $\phi_N = 0.1$, the diffusion coefficient is reduced by $\approx 40\%$ for NPs of diameter $\sigma_N = 2$ and 3, and by $\approx 30\%$ for NPs of diameter $\sigma_N \ge 3$ (Fig. 13a). The most interesting characteristic of $D_N/D_{N0}$ is however the dependence on $\sigma_N$ at fixed $\phi_N$. To better visualize this, we have interpolated between the points in Fig. 13a in order to obtain approximately the reduced NP diffusion coefficient as a function of the NP diameter $\sigma_N$ at constant $\phi_N$ (Fig. 13b). The ratio $D_N(\sigma_N)/D_{N0}$ shows an initial increase with increasing $\sigma_N$, then an inflection point at $\sigma_N \approx 2.5$, and finally it reaches a plateau for $\sigma_N \gtrsim 4$. Such a peculiar behavior can be interpreted in the following manner: At $\sigma_N \approx 1$, increasing the NP diameter at fixed volume fraction has the effect of reducing the polymer-NP interface, and therefore decreasing the total interaction energy between polymers and NPs, resulting in an enhanced NP diffusion. When the NP size becomes larger than the mesh size $\xi \approx 1.3$ (the exact value depends on $\phi_m$), the motion of NPs starts to be geometrically hindered by the polymer segments [11, 37], and as a result the dependence of $D_N/D_{N0}$ on the NP size weakens. Then, when the NPs become large enough, since the surface-to-volume ratio becomes smaller, the importance of the energetic contribution to the diffusion coefficient starts to decrease, resulting in another increase in the diffusion coefficient. Finally, for large NPs energy becomes irrelevant and $D_N/D_{N0}$ is completely controlled by geometry, and therefore is constant at constant NP volume fraction.

## V. HIGHER NANOPARTICLE VOLUME FRACTIONS

If we keep increasing the number of NPs $N_N$ while keeping pressure and number of polymers constant, the volume of the simulation box will eventually start to increase proportionally to $N_N$. As a consequence, the NP volume fraction $\phi_N$ will reach a plateau, $\phi_N = \phi_N^{\max}(\sigma_N)$, which corresponds to the value of $\phi_N$ for a pure NP system at temperature $T = 1.0$ and pressure $P = 0.1$. This situation corresponds approximately to the one depicted in Fig. 2d. If a standard LJ potential was used for the NP-NP interaction, $\phi_N^{\max}$ would not depend on $\sigma_N$, since the interaction potential would only depend on the ratio $\sigma_N/r$ and all systems would be equivalent apart from a trivial distance rescaling. However, the expanded LJ potential, Eq. (3), does not simply depend on $\sigma_N/r$; therefore, pure NPs systems with the same $T$ and $P$ are not equivalent.

In Fig. 14, we show the reduced volume of the simulation box $V/V_0$, where $V_0 = V(\phi_N = 0)$, as a function of $N_N$, for different values of $\sigma_N$. For $\sigma_N \ge 3$, the volume increases monotonically with the NP number. For $\sigma_N = 1$ and 2, on the other hand, there is a range of $N_N$ values in which we observe a decrease in volume (see Sec. III B). In the inset of Fig. 14 we show the NP volume fraction $\phi_N$ as a function of $N_N$. Initially, the volume is almost constant and therefore as a good approximation $\phi_N = \pi\sigma_N^3 N_N / 6V_0$. Then, for larger values of $N_N$, the volume starts to increase proportionally to $N_N$ and $\phi_N$ reaches the plateau $\phi_N = \phi_N^{\max}$.

<figure><img src="image_15.png" /><figcaption>FIG. 14. Reduced volume of the simulation box versus NP number $N_N$. Inset: NP volume fraction $\phi_N$ versus $N_N$. Filled symbols: good NP dispersion. Open symbols: poor NP dispersion.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIG. 15. Normalized diffusion coefficient of NPs and polymers as a function of NP volume fraction for different values of the NP diameter $\sigma_N$. Filled symbols: good NP dispersion. Open symbols: poor NP dispersion. Dashed line: approximate value of $\phi_N^{\max}$ for $\sigma_N = 3$. The data for $D_N/D_{N0}$ have been shifted up by 0.6 for clarity.</figcaption></figure>

Since the number of polymer chains is constant, the increase of the volume at large $N_N$ results in a decrease of polymer volume fraction $\phi_m$, and therefore in a decrease of the polymer-NP interface per unit volume. This
in turn causes a weakening of the polymer-mediated attractive interaction between NPs and consequently an increase of the free volume $V_{\text{free}} = V[1 - (\phi_N + \phi_m)]$. Both of these mechanisms result in an increase in the polymer and NP diffusivities. We can observe this effect in Fig. 15, where we show the normalized diffusion coefficient of both polymers and NPs for all the simulated systems, including those where good NP dispersion is not realized (open symbols): $D/D_0$ reaches a minimum corresponding to the value of $\phi_N$ at which the volume starts to increase, and it continues to grow as $\phi_N^{\text{max}}$ is approached. We note however that for $\sigma_N = 7$, $D/D_0$ shows a monotonic decrease. While in all the other cases we found that the pure NP system at $T = 1.0$ and $P = 0.1$ is a liquid, for $\sigma_N = 7$ it is a crystal, which means that as we approach $\phi_N^{\text{max}}$ the ratio $D/D_0$ will decrease and eventually settle to a very small value.

## VI. SUMMARY AND CONCLUSIONS

We have carried out molecular dynamics simulations of a mixture of polymers and spherical nanoparticles (NPs) in a wide range of NP volume fractions ($\phi_N$) and for different NP diameters ($\sigma_N = 1, 2, 3, 4, 5$, and $7\sigma$, where $\sigma$ is the monomer size).

We have studied the structural properties of polymers and NPs, identifying the range of values of $\phi_N$ in which the NPs are well dispersed in the solution. In agreement with previous studies [51, 52, 68–70], we have found that the NPs of diameter $\sigma_N < 2R_g$ act like a good solvent, swelling the polymers. Surprisingly, however, we have also observed that NPs of the same size as that of the monomers, $\sigma_N = 1$, have the opposite effect, acting like a poor solvent. This is due to the high surface-to-volume ratio of small NPs, which causes the energetic contribution (which promotes chain contraction) to become stronger than the excluded volume contribution (which promotes chain expansion). Therefore, this effect is expected to depend on the strength of the monomer-NP interaction.

We have then analyzed the dynamical properties of the system, and in particular the diffusion coefficient of the NPs, $D_N$, and of the centers of mass of the polymer chains, $D_p$. We found that in the presence of good NP dispersion both $D_N$ and $D_p$ decrease monotonically with increasing NP volume fraction.

The reduction of $D_p$ in the presence of NPs can be well described, at fixed temperature, with a modification of the confinement parameter approach proposed by Composto and coworkers[19–22]. Having observed that the conventional “static” confinement parameter, defined as the ratio between the interparticle distance $h$ and $2R_g$, fails to capture the dynamics of the chains at different temperatures, we propose a *dynamic* confinement parameter $h/\lambda_d(T)$, where $\lambda_d(T)$ represents a cooperativity length scale which increases with decreasing temperature. When plotted against $h/\lambda_d(T)$, all the data for $\sigma_N \ge 2$

fall on the same master curve, $D_p = D_{p0}[1 - \exp(-h/\lambda_d)]$. Deviations from the master curve only appear at $\sigma_N = 1$.

The behavior of the reduced NP diffusion coefficient shows a complex dependence on $\sigma_N$, with an initial increase followed by an inflection point around $\sigma_N = 2.5$ and finally a plateau for $\sigma_N \gtrsim 4$. We speculate that this behavior results from an interplay between energetic and entropic contributions, with the latter depending on the mesh size of the solution, $\xi$.

We have also studied the single NP diffusion coefficient, $D_{N0} = D_N(\phi_N \to 0)$, performing additional simulations with $\sigma_N = 10, 12, 14$ at low NP volume fraction. For $\sigma_N < 2R_g$ the diffusion is faster than what is predicted by the Stokes-Einstein formula, according to which $D_{N0} \propto \sigma_N^{-1}$. Stokes-Einstein behavior is recovered when $\sigma_N \approx 2R_g$, in agreement with previous studies [12, 17, 18, 24, 27]. Theoretical studies predicted that for $\sigma_N$ sufficiently smaller than $2R_g$ the NP diffusion coefficient should decrease as $D_{N0} \propto \sigma_N^{-3}$ [11, 13, 37]; Cai et al. [11], in particular, predicted that this behavior should be observed in the range $\xi < \sigma_N < 2R_g$. However, it is not clear whether our data confirm these predictions or not: The main obstacle to find a conclusive answer is that $D_N$ shows a dependence on the mass density of the NPs for $\sigma_N \lesssim 7$, a dependence which has been reported in previous studies [27, 38, 39]. This mass dependence becomes stronger with decreasing $\sigma_N$, making the interpretation of the data difficult. In light of this observation, we believe that care has to be taken when comparing the results of simulations to those of experiment and to theory, especially at low monomer volume fractions.

In conclusion, the behavior of a polymer-NP solution is very rich, and it is dictated by the value of many different parameters, most importantly the polymer-NP interaction and the NP size. In particular, since changing the polymer-NP interaction can lead to a qualitatively very different dynamical behavior, we believe that systematic studies should be conducted in order to clarify the role of interactions in the dynamics of nanocomposites. Comparison of simulation data with the available theories is complicated by many factors, such as subtleties in the definition of important quantities, e.g. the effective hydrodynamic radius and the mesh size, and the presence of broad cross-overs and inertial effects (mass dependence of the NP diffusivity). Other fundamental issues that remain to be fully clarified are the relevance of hydrodynamic interactions (recently investigated in Ref. 60) and finite-size effects. We thus believe that the investigation of polymer-NP composites will remain an important topic of research also in the future.
# APPENDIX

## AI. INTERPARTICLE DISTANCE AND PORE-SIZE DISTRIBUTION

In a polymer nanocomposite, the interparticle distance $h$ is the average distance between the surfaces of neighboring NPs in the system. In the literature [19–22], $h$ has often been defined using Eq. (7), which we reproduce here:

$$h^{\text{th.}} = \sigma_N \left[ \left( \frac{\phi_N^M}{\phi_N} \right)^{1/3} - 1 \right],$$

where $\phi_N^M$ is the maximum achievable NP volume fraction [20]. There is however an evident problem with the above expression: It is a priori not clear at all what value should be used for $\phi_N^M$. Taking the NPs to behave approximately as hard spheres (an approximation that in our case is justified by the very steep NP-NP potential), there are several possibilities, like the close-packing value [91] $\phi_{cp} = \pi/\sqrt{18} \approx 0.740$, corresponding to an fcc or hcp lattice, or the volume fraction of some other crystal lattice, like the bcc ($\phi_{bcc} = \pi\sqrt{3}/8 \approx 0.680$) or the simple cubic ($\phi_{sc} = \pi/6 \approx 0.524$). The value $\phi_N^M = \pi/6$ was used in Ref. 92, one of the first to apply Eq. (7) to polymer nanocomposites. More recently, $\phi_{rcp} \approx 0.637 \approx 2/\pi$, that should correspond to a random close packing (RCP) of hard spheres, has often been invoked in the definition of $h$ [19–22]. However, it has been shown by Torquato et al. that the concept of RCP is ill-defined, and that different procedures can result in different values for $\phi_{rcp}$, ranging from 0.6 to 0.68 [64]. This issue could be solved, as suggested by Torquato et al., by replacing the ill-defined concept of RCP with that of maximally random jammed (MRJ) structure [64, 93]. For monodispere hard spheres, this redefinition should lead to a unique value, $\phi_{MRJ} = 0.642$ [93]. However, the problem of a priori assigning a certain value to $\phi_N^M$ remains.

We propose therefore a different way to define $h$, which relies on the concept of pore-size distribution [65]. The pore-size probability density function (PDF) $P(\delta)$ of a system consisting of two phases is defined such that $P(\delta)d\delta$ represents the probability that a randomly chosen point in the phase of interest lies at a distance between $\delta$ and $\delta + d\delta$ of the nearest point on the interface between the two phases [65]. It is clear from this definition that the typical interparticle distance $h$ for the NPs in a polymer nanocomposite should correspond approximately to the typical pore size.

There is another clue that suggests the identification of $h$ with some quantity derived from $P(\delta)$. For a system of randomly distributed overlapping spheres of radius $R$ with number density $\rho$, the pore-size PDF can be computed explicitly [94]:

<figure><img src="image_17.png" /><figcaption>FIG. 16. (color online) Pore-size probability density function of the NPs for $\sigma_N = 4$. Cartoon: schematic representation of how the pore-size distribution is calculated.</figcaption></figure>

$$P^{\text{os}}(\delta) = \frac{3\phi}{e^{-\phi}R} \left(\frac{\delta}{R} + 1\right)^2 \exp\left[-\phi \left(\frac{\delta}{R} + 1\right)^3\right]. \quad (16)$$

In this expression, $\phi = \pi(2R)^3\rho/6$ is the “volume fraction” of the spheres (although since the spheres can overlap, this does not correspond to their real volume fraction).

Let us now define $\delta_{max}$ as the value of $\delta$ for which $P(\delta)$ is maximized, i.e., $P(\delta_{max}) = \max_{\delta}\{P(\delta)\}$. For this model, we have $\delta_{max}^{os} = R[(2/3\phi)^{1/3} - 1]$, and therefore we can define a typical pore diameter as

$$2\delta_{\text{max}}^{\text{os}} = 2R \left[ \left( \frac{2}{3\phi} \right)^{1/3} - 1 \right]. \quad (17)$$

By comparing Eqs. (7) and (17) and making the identifications $2R = \sigma_N$ and $\phi = \phi_N$, we note that $2\delta_{max}^{os} = h^{\text{th.}}$, provided that we choose $\phi_N^M = 2/3 \approx 0.667$. Given this quite remarkable connection between $2\delta_{max}^{os}$ and $h^{\text{th.}}$, and given that the definition of $\delta_{max}$ does not present the same problems that affect Eq. (7), we are naturally lead to define the interparticle distance as $h = 2\delta_{max}$, where $2\delta_{max}$ is not computed using Eq. (17), which is strictly valid only for the overlapping spheres system, but rather evaluated directly from the pore-size PDF $P(\delta)$ obtained from the data.

The algorithm employed to obtain $P(\delta)$ from the data is described in Ref. 65: (1) A random point $(x, y, z)$ in the NP-free phase is chosen (the point must lay at distance $r > \sigma_N/2$ from every NP). (2) The radius $\delta$ of the largest sphere centered in $(x, y, z)$ that does not intersect any NP is recorded. (3) Step 1 and 2 are repeated many times, and a histogram of the radii is created. (4) The pore-size PDF $P(\delta)$ is obtained by normalizing the histogram. We note that in order for step 2 to be well-defined, we need to approximate the NPs as hard spheres of diameter
<figure><img src="image_18.png" /><figcaption>FIG. 17. (color online) Mean squared displacement of the NPs when varying the NP mass $m_N$ at fixed diameter $\sigma_N$, in the case $\sigma_N = 2, \phi_N = 0.024$ (a) and $\sigma_N = 7, \phi_N = 0.092$ (b).</figcaption></figure>

$\sigma_N$; this approximation is justified by the fact that the NP-NP potential is very steep.

In Fig. 16 we show the pore-size PDF for $\sigma_N = 4$. As already shown in Fig. 5, the interparticle distance obtained from the pore-size PDF and the one calculated using Eq. 7 are very similar. This means that Eq. 7 can be used to obtain an estimate of the “true” interparticle distance, despite the problems affecting its definition. We also mention that the pore-size PDF can also be extracted from experimental data: For example, Rintoul et al. used X-ray microtomography to obtain a three dimensional digitalized image of a porous magnetic gel and computed from it the pore-size PDF and other statistical correlation functions [94].

## AII. MASS DEPENDENCE OF NP AND POLYMER DIFFUSIVITIES

In previous studies of polymer nanocomposites [27, 39] and binary soft-sphere liquids [38] it has been shown that outside of the Stokes-Einstein regime not only the diameter, but also the mass density of a particle can affect its dynamics. In order to study this effect, we performed some simulations at low NP volume fraction and changed

<figure><img src="image_19.png" /><figcaption>FIG. 18. Diffusion coefficient of the NPs (a) and of the polymers (b) when varying the NP mass $m_N$ at fixed diameter $\sigma_N$. The values of the NP volume fraction $\phi_N$ are 0.024, 0.054 and 0.092 respectively for $\sigma_N = 2, 4$ and 7. Open symbols: simulation without rescaling of the viscous friction coefficient $\Gamma_N$. Dashed lines: slopes $-0.25$ (a) and $-0.03$ (b) for reference.</figcaption></figure>

the NP mass while leaving the diameter fixed.

The results for the mean squared displacement (MSD) $\langle r^2(t) \rangle$ of the particles are reported in Fig. 17 for the cases $\sigma_N = 2, \phi_N = 0.024$ and $\sigma_N = 7, \phi_N = 0.092$. At short times, the motion of the particles is ballistic and the mass dependence of the MSD is trivial:

$$ \langle r_N^2(t) \rangle = \langle v_N^2 \rangle t^2 = \frac{3k_B T}{m_N} t^2, \quad (18) $$

where $\langle v_N^2 \rangle$ is the mean squared speed of the NPs. At longer times, when the motion becomes diffusive, i.e. $\langle r^2(t) \rangle \propto t$, we can observe a much more interesting effect: While the motion of the larger NPs is unaffected by a change in the mass (Fig. 17b), the motion of smaller NPs presents a clear mass dependence (Fig. 17a).

This result can be better appreciated in Fig. 18a, where we show the NP diffusion coefficient $D_N$ as a function of the NP mass for different values of $\sigma_N$ and $\phi_N$ (here also the case $\sigma_N = 4, \phi_N = 0.054$ is shown). One recognizes that the mass dependence of the diffusion coefficient becomes weaker when the NP diameter is increased, and for
$\sigma_N = 7$ it has disappeared almost completely. The result that the mass dependence of the long time diffusion coefficient is stronger for smaller particles is in agreement with previous studies [27, 38].

Our interpretation of this result is that large particles are forced to wait for the polymers to relax in order to diffuse through the solution: The constraints on their motion are purely geometric, and mass plays no role. Smaller particles, on the other hand, can “slip” between the polymers; the probability that they find a passage to slip through in a certain time interval increases with their average velocity, which at equilibrium decreases as $m_N^{-1/2}$, as it follows from the equipartition of energy: $\langle v_N^2 \rangle = 3k_B T/m_N$. We expect this effect to be strongly suppressed when the polymer density is increased, and to be negligible for the motion of NPs of diameter $\sigma_N \gtrsim \sigma$ in a melt.

Although in Refs. 27, 38 it was claimed that the mass dependence disappears for high values of $m_N$, in our case for $\sigma_N = 2$ there is no hint that the mass dependence vanishes for larger mass values. In the case of Ref. 38, however, the disappearance of the mass dependence for high values of the mass is only apparent, as it becomes clear once the data are plotted in double logarithmic instead of linear scale (not shown). We believe that this could also be the case for Ref. 27, since also there the data are only reported in linear scale. From the analysis of our data and of those of Ref. 38, it seems that not only the mass dependence does not disappear when $m_N$ is increased, but on the contrary it becomes stronger.

The diffusion coefficient of the polymer chains $D_p$ is almost unaffected by changes in the NP mass density, as we can observe from Fig. 18b, although it is possible to see a very weak decrease of $D_p$ for $\sigma_N = 2$. It is possible that the slowing down of the NPs has an effect on the dynamics of the polymers, but this effect is very weak at low NP concentrations. Further study should be dedicated to clarifying this point.

Finally, we have made some tests to determine whether the observed mass dependence is an artifact resulting from the scaling of the friction coefficient of the Langevin thermostat (see Sec. SI in the S.I.). To this aim, we ran some simulations where the friction coefficient of the NPs, $\Gamma_N$, was kept constant and equal to that of the monomers: $\Gamma_N = \Gamma_m = 0.1$. The result is included in Fig. 18 (open blue circles). One sees that the effect of mass density is still present, but using $\Gamma_N = 0.1$ (fixed) has the effect of reducing the diffusion coefficient of the NPs, as one expects since this means that more massive NPs experience a higher solvent viscosity (we recall that $\eta_s \propto \Gamma_N m_N / \sigma_N$).

In conclusion, we have shown that the effect of mass density on the dynamical properties should be taken into careful consideration when performing molecular dynamics simulations of multi-component systems, such as polymer mixtures, binary fluids and solutions with explicit solvents. We think that the dependence of this effect on polymer and NP density and on polymer-NP interaction

should be more thoroughly investigated in order to gain a better understanding of polymer-NP mixtures from the point of view of molecular dynamics simulations.

## ACKNOWLEDGMENTS

We thank K. Schweizer, J. Oberdisse, L. Cipelletti, and M. Ozawa for fruitful discussions. 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). Part of the simulations were performed at the Center of High Performance Computing MESO@LR in Montpellier. The snapshots in Fig. 2 have been realized with VMD [95].

## SUPPLEMENTARY INFORMATION

### SI. SCALING OF THE FRICTION COEFFICIENT OF THE THERMOSTAT

<figure><img src="image_20.png" /><figcaption>FIG. S1. (color online) Mean-squared displacement (divided by time) of a particle of diameter $\sigma$ and mass $m$ diffusing in pure (implicit) solvent. The long-time diffusion coefficient depends only on $\sigma$.</figcaption></figure>

As discussed in the main text, we have used for our simulations a Langevin thermostat, so that the force experienced by particle $i$ is

$$m_{\alpha} \ddot{\mathbf{r}}_{\alpha,i} = -\nabla_i U(\{\mathbf{r}_{\alpha,k}\}) - m_{\alpha} \Gamma_{\alpha} \dot{\mathbf{r}}_{\alpha,i} + \sqrt{2m_{\alpha} \Gamma_{\alpha} k_B T} \zeta(t), \quad (S1)$$

where $\alpha$ denotes the particle type ($m$ = monomer, $N$ = NP). The mass and diameter of the monomers are respectively $m_m = m$ and $\sigma_m = \sigma$. $U(\{\mathbf{r}_{\alpha,k}\})$ is the total interaction potential acting on the particle, with $\{\mathbf{r}_{\alpha,k}\}$ representing the set of coordinates of all the particles in the system. The term $\sqrt{2m_{\alpha} \Gamma_{\alpha} k_B T} \zeta(t)$ is a random force which represents collisions with solvent molecules, while
$\Gamma_\alpha$ is a viscous friction coefficient, which is related to the viscosity of the implicit solvent $\eta_s$ by

$$ \Gamma_{\alpha} = \frac{C\eta_s\sigma_{\alpha}}{m_{\alpha}}, \qquad (S2) $$

where $\sigma_\alpha$ is the diameter of the particle and C a co-efficient which depends on the hydrodynamic boundary conditions.

In our simulations we made the assumption that every particle (monomer or NP) interacts with the implicit solvent as if it was a continuum with fixed $\eta_s$. Therefore, from Eq. (S2) we get

$$ \eta_s = \frac{\Gamma_m m}{C\sigma} = \frac{\Gamma_N m_N}{C\sigma_N} = \text{const.} \Rightarrow \Gamma_N = \Gamma_m \cdot \left( \frac{\sigma_N m}{m_N \sigma} \right). \qquad (S3) $$

With this choice of the friction coefficient, the long-time diffusion coefficient of a particle in the pure (implicit) solvent $D_{s,\alpha}$ follows the Stokes-Einstein law with viscosity $\eta_s$:

$$ D_{s,\alpha} = \frac{k_B T}{\Gamma_{\alpha} m_{\alpha}} = \frac{k_B T}{C \eta_s \sigma_{\alpha}}. \qquad (S4) $$

This way, $D_{s,\alpha}$ depends only on the diameter of the particle, and not on its mass, as we can see from Fig. S1, where the mean-squared displacement (divided by time) of particles of various masses and diameters in the pure solvent is shown. We note that since $\Gamma_m = 0.1$ and $m, \sigma = 1$, the numerical value of the solvent viscosity is $\eta_s = 0.1/3\pi = 1.06 \cdot 10^{-2}$.

## SII. MONOMER-MONOMER RADIAL DISTRIBUTION FUNCTION

In Figs. S2a-b, we report the monomer-monomer radial distribution function $g_{mm}(r)$ for $\sigma_N = 2$ and 4. One sees that at low NP volume fraction the most prominent features of $g_{mm}(r)$ are a sharp peak at $r \approx r_b = 0.96$ (first nearest neighbor distance in a chain) and a smaller one at $r \approx 2r_b = 1.92$ (second nearest neighbor). When the NP volume fraction is increased, the height of these two peaks increases. The reason is that, while the structure of the chain at the length scale $r \lesssim 2r_b$ remains almost unchanged when $\phi_N$ increases, the monomer density $\rho_m$ decreases, because the volume increases and the number of monomers is fixed. Since $g_{mm}(r)$ contains a factor $\rho_m^{-1}$ (see Eq. (5) in the main text), this results in an increase of this function for $r \lesssim 2r_b$. Incidentally, this is why the radial distribution function of a single polymer, used to derive the function $p(r)$ (Sec. III B in the main text) is a more sensible quantity when studying chain conformation.

Since $g_{mm}(r)$ is dominated by the two peaks at $r_b$ and $2r_b$ and shows only very small fluctuations around

$g_{mm}(r) = 1$ at larger $r$, we plot in Figs. S2c-d the function $|g_{mm}(r) - 1|$, hence allowing to detect more easily the structure at large $r$. For $\sigma_N = 4$, Fig. S2d, $|g_{mm}(r) - 1|$ shows very clearly at intermediate and high $\phi_N$ a long-range modulation with typical wavelength $\sigma_N$, due to the presence of the NPs. Notice that, since we have taken the absolute value, the wavelength must be calculated as the distance between the n-th peak and the $(n+2)$th. For $\sigma_N = 2$, the presence of this modulation is less clear, because the size of the NPs is close to the monomer size and as a consequence the signal coming from the NPs cannot be distinguished well from the one coming from the monomers themselves.

## SIII. SYSTEMATIC STUDY OF THE DIFFERENT DYNAMICAL REGIMES

In Sec. IV B of the main text we have studied the mean-squared displacement (MSD) $\langle r^2(t) \rangle$ of the NPs and of the centers of mass of the polymers, identifying three different regimes: a short-time ballistic regime, $\langle r^2(t) \rangle \propto t^2$, an intermediate-time subdiffusive regime, $\langle r^2(t) \rangle \propto t^{\beta}$, with $\beta < 1$, and a long-time diffusive regime, $\langle r^2(t) \rangle \propto t$. In order to study in a more systematic way these regimes and the transition between them, we introduce the function

$$ \beta_{\alpha}(t) = \frac{d \log \langle r_{\alpha}^{2}(t) \rangle}{d \log t} \qquad (S5) $$

where $\alpha = N, p$ respectively for the NPs and for the polymers. This quantity represents the slope of the MSD in a log-log scale, and it is a generalization of the above mentioned subdiffusive exponent $\beta$.

In Fig. S3a we show $\beta_p$ for the simulated systems. The presence of the transient subdiffusive regime with $\beta_p \approx 0.87$, already mentioned in the main text, is clearly visible. The dynamics of the polymers becomes again diffusive ($\beta_p = 1$) after a time $t \approx 10^4$, which we identify with the relaxation time of the chain [33]. We can observe that $\beta_p$ is not much affected by the NP volume fraction $\phi_N$ and size $\sigma_N$.

In Fig. S3b we report $\beta_N$. As already discussed in the main text, the dynamics of the NPs is more strongly influenced by the values of $\phi_N$ and $\sigma_N$ than that of the polymers. For $\sigma_N = 5, 7$, we can observe a subdiffusive transient ($\beta_N < 1$) appearing very clearly at intermediate times. For $\sigma_N < 5$, the effect is much smaller.

## SIV. VAN HOVE FUNCTION

The van Hove function $G(\mathbf{r}, t)$ [62] can be written as the sum of a self and a distinct part: $G(\mathbf{r}, t) = G_s(\mathbf{r}, t) + G_d(\mathbf{r}, t)$, where
<figure><img src="image_21.png" /><figcaption>FIG. S2. (color online) Monomer-monomer radial distribution function for $\sigma_N = 2$ (a) and $\sigma_N = 4$ (b) and different values of the NP volume fraction $\phi_N$.</figcaption></figure>

$$G_s(\mathbf{r}, t) = \frac{1}{M} \sum_{i=1}^{M} \langle \delta[\mathbf{r} - \mathbf{r}_i(t) + \mathbf{r}_i(0)] \rangle \quad (S6)$$

and

$$G_d(\mathbf{r}, t) = \frac{1}{M} \sum_{\substack{i=1 \\ j \neq i}}^{M} \langle \delta[\mathbf{r} - \mathbf{r}_j(t) + \mathbf{r}_i(0)] \rangle. \quad (S7)$$

The self part of the van Hove function represents the time-dependent spatial autocorrelation of a particle, while the distinct part represents the time-dependent spatial pair correlation. As usual, we will consider the spherical average of these two quantities: $G_s(r, t)$ and $G_d(r, t)$. The function $4\pi r^2 G_s(r, t)$ represents the probability to find a particle at time $t$ a distance $r$ from its original position. We note that $G_s(\mathbf{r}, 0) = \delta(\mathbf{r})$ and $G_d(\mathbf{r}, 0) = \rho g(\mathbf{r})$, where $g(\mathbf{r})$ is the pair correlation function [62].

For both small and large values of $t$, the self part of the van Hove function is a Gaussian [62]:

$$G_s(r, t) = \Gamma_s(r, t) = \left( \frac{3}{2\pi\langle r^2(t) \rangle} \right)^{3/2} \exp \left( -\frac{3r^2}{2\langle r^2(t) \rangle} \right). \quad (S8)$$

We can therefore define a rescaled self van Hove function which also preserves the probability with the following change of variables:

$$\begin{align} r &\rightarrow r' = r \cdot \left( \frac{3}{2\langle r^2(t) \rangle} \right)^{1/2} \tag{S9} \\ G_s &\rightarrow G'_s = G_s \cdot \left( \frac{3}{2\langle r^2(t) \rangle} \right)^{-3/2} \end{align}$$

If $G_s$ is Gaussian, the result of the transformation (S9) is

$$\Gamma'_s(r') = \pi^{-3/2} e^{-r'^2}, \quad (S10)$$

i.e., the distribution is independent of time.

In Fig. S4a, we show the rescaled van Hove function for the case $\sigma_N = 1$, $\phi_N = 0.05$ (but for other parameters we find the same qualitative behavior). In Fig. S4b, we report the ratio between the rescaled self van Hove function
<figure><img src="image_22.png" /><figcaption>FIG. S3. Slope of the MSD in log-log scale, $\beta_\alpha$ (Eq. (S5)), for the centers of mass of the polymers (a) and for the NPs (b). For the sake of clarity, for $\sigma_N < 7$ every set of curves has been shifted with respect to the previous one, by 0.5 in (a) and by 0.75 in (b). Different curves in the same set correspond to different NP volume fractions $\phi_N$. The dashed curves correspond to the lowest $\phi_N$, while the thick continuous curves correspond to the highest $\phi_N$.</figcaption></figure>

$G'_s(r, t)$ of the NPs and the same quantity in the Gaussian approximation $\Gamma'_s(r, t)$, for the same system. We observe that the shape of $G'_s(r, t)$ at short and long times is indeed very close to a Gaussian ($G'_s(r, t) = \Gamma'_s(r, t)$), and that the largest deviation from Gaussian behavior occurs when the dynamics of the NPs starts to be diffusive (in this case, at $t \approx 1$). These deviations are found to be most pronounced at large $r$, i.e., the NPs move a bit further than expected from a Gaussian approximation.

To better quantify how dissimilar $G_s(r, t)$ is from a Gaussian, it is customary to define a non-Gaussian parameter [80],

$$ \alpha_2(t) = \frac{3\langle r^4(t) \rangle}{5\langle r^2(t) \rangle^2} - 1, \quad (S11) $$

where

$$ \begin{aligned} \langle r^n(t) \rangle &= \frac{1}{M} \sum_{i}^{M} \langle |\mathbf{r}_i(t) - \mathbf{r}_i(0)|^n \rangle \\ &= 4\pi \int_{0}^{\infty} G_s(r, t) r^{2+n} dr. \end{aligned} \quad (S12) $$

If $G_s(r, t)$ is Gaussian, Eq. (S8), we have $\alpha_2 = 0$. Therefore, high values of $|\alpha_2|$ indicate a significant non-Gaussian behavior.

In Fig. S5, we show $\alpha_2(t)$ of the NPs for several values of the NP diameter $\sigma_N$ and of the NP volume fraction $\phi_N$. The largest departure from Gaussian behavior happens when the dynamics of the NPs starts to be diffusive, in agreement with what is observed from $G'_s(r, t)$. Both at short and long times $\alpha_2 \approx 0$, as expected. We notice that the maximum deviation from Gaussian behavior (the maximum of the curves in Fig. S5) becomes larger when $\phi_N$ is increased. This trend shows that the structure of the surrounding polymer mesh and the presence of nearby NPs both contribute to the non-Gaussian behavior. Moreover, increasing the NP size at fixed $\phi_N$ generally reduces the magnitude of $\alpha_2$. This is reasonable, since a large NP interacts with a large number of monomers and other NPs and thus “feels” an averaged interaction, which results in a reduction of the dynamical fluctuations and therefore of $\alpha_2$. One exception to this trend is $\sigma_N = 7$ at high NP volume fraction (Fig. S5d). The reason for this is likely that the system is approaching crystallization. Apart from the case $\sigma_N = 7$, $\phi_N = 0.467$, we always have $|\alpha_2| < 0.2$, and we can therefore state that the dynamics of the NPs is, to a good approximation, Gaussian. The non-Gaussian parameter of the polymer chains (not shown) always satisfies $|\alpha_2| < 0.1$, therefore also the dynamics of the polymers is approximately Gaussian. We also mention that in this case no clear dependence of $\alpha_2$ on $\phi_N$ and $\sigma_N$ is observed.

Finally, in Fig. S6 we show the distinct part of the van Hove function, $G_d(r, t)$, of the NPs for some selected systems. The relaxation happens in a way which is very similar to that observed in simple, non-supercooled liquids [80], in that we observe that in all cases the correlation hole at $r = 0$ is slowly filled as $t$ is increased. Since there is no evidence for the presence of a peak at $r = 0$, we conclude that hopping dynamics is absent in the studied systems [96].

## SV. POLYMER DIFFUSION COEFFICIENT AS A FUNCTION OF NP VOLUME FRACTION

In Sec. IV B of the main text we have seen how, in the regime of good NP dispersion, the decrease of the polymer diffusion coefficient $D_p$ relative to the diffusion coefficient in the pure polymer solution $D_{p0}$ can be empirically described by the function $D_p = D_{p0}[1 - (\phi_N/\phi_{p0})^\alpha]$
<figure><img src="image_23.png" /><figcaption>FIG. S4. (a) Rescaled self part of the van Hove function (Eqs.(S6)-(S9)) for $\sigma_N = 1$, $\phi_N = 0.05$ (b) Ratio between the rescaled self part of the van Hove function and the same quantity in the Gaussian approximation (Eq. (S10)), as a function of $r'^2$, for the same system.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>FIG. S5. (color online) Non-Gaussian parameter $\alpha_2$, Eq. (S11), for various NP diameters $\sigma_N$ and NP volume fractions.</figcaption></figure>

(see Fig. 10a in the main text). In Table S1, we report the fit parameters for this empirical relation.

## SVI. COMPARISON OF POLYMER AND NP DIFFUSIVITIES

In Fig. S7a, we compare the diffusion coefficient of the NPs, $D_N$, and of the CM of the polymers, $D_p$. One sees that, in almost all the systems we considered, $D_N > D_p$, i.e. the NPs move faster than the polymer chains. An
<figure><img src="image_25.png" /><figcaption>FIG. S6. (color online) Distinct part of the van Hove function, Eq. (S7), for some selected systems. The curves are equally spaced in logarithmic time ($t$). In (a) and (c), we have indicated some of the $t$ values associated to the curves. The minimum value of $t$ is for all the figures $t_{min} = \delta t = 0.003$ (integration time step), while the maximum value, $t_{max}$, is reported in the legend.</figcaption></figure>

<table><caption>TABLE S1. Fit parameters for $D_p = D_{p0}[1 - (\phi_N/\phi_{p0})^\alpha]$.</caption><thead><tr><th>σ<sub>N</sub></th><th>φ<sub>p0</sub></th><th>α</th></tr></thead><tbody><tr><td>1</td><td>0.104</td><td>0.763</td></tr><tr><td>2</td><td>0.231</td><td>0.768</td></tr><tr><td>3</td><td>0.364</td><td>0.812</td></tr><tr><td>4</td><td>0.471</td><td>0.851</td></tr><tr><td>5</td><td>0.519</td><td>1.006</td></tr><tr><td>7</td><td>0.594</td><td>1.150</td></tr></tbody></table>

## SVII. ADDITIONAL DETAILS ON THE SIMULATED SYSTEMS

In Tab. S2 we report some additional details on the simulated systems. $N_N$: number of NPs; $\sigma_N$: NP diameter; $L$: side of the cubic simulation box; $\phi_m$: monomer volume fraction; $\phi_N$: NP volume fraction; $R_g$: average radius of gyration of the polymers; $h$: interparticle distance.

exception is $\sigma_N = 7$ at high densities; however, we know that at higher NP volume fraction the NPs form in this case a crystal and therefore $D_N$ becomes very small.

Figure S7b shows the NP diffusion coefficient $D_N$ as a function of the polymer center of mass diffusion coefficient $D_p$. One can observe that there is a strong correlation between $D_N$ and $D_p$, which can be empirically described via a power law: $D_N \propto D_p^\alpha$, where $\alpha$ increases with increasing NP diameter (inset of Fig. S7b). This correlation suggests that there is a coupling between the long time diffusivities of NPs and the centers of mass of the polymers, as it was recently proposed by Chen *et al.* [79].
<figure><img src="image_26.png" /><figcaption>FIG. S7. (a) Ratio between the NP and polymer diffusion coefficients as a function of the NP volume fraction and for different values of the NP diameter $σ_N$. (b) NP diffusion coefficient $D_N$ versus the chain CM diffusion coefficient $D_p$. The correlation between $D_N$ and $D_p$ takes the form of a power law, $D_N \propto D_p^\alpha$, where $\alpha$ increases with the NP diameter. Inset: $\alpha$ as a function of the NP diameter.</figcaption></figure>
<table><caption>TABLE S2. Details of the simulated systems. Systems in which there is poor NP dispersion (p.d.) are denoted by an asterisk in the last column.</caption><thead><tr><th>N<sub>N</sub></th><th>σ<sub>N</sub></th><th>L</th><th>φ<sub>m</sub></th><th>φ<sub>N</sub></th><th>D<sub>N</sub></th><th>D<sub>p</sub></th><th>R<sub>g</sub></th><th>h</th><th>p.d.</th></tr></thead><tbody><tr><td>0</td><td>-</td><td>56.32</td><td>0.1466</td><td>0.0000</td><td>-</td><td>0.0114</td><td>6.282</td><td>-</td><td>-</td></tr><tr><td>30000</td><td>1</td><td>57.08</td><td>0.1408</td><td>0.0845</td><td>0.2793</td><td>0.0046</td><td>6.4681</td><td>0.8980</td><td>*</td></tr><tr><td>15000</td><td>1</td><td>54.03</td><td>0.1660</td><td>0.0498</td><td>0.2824</td><td>0.0050</td><td>6.2133</td><td>1.2760</td><td></td></tr><tr><td>7500</td><td>1</td><td>54.57</td><td>0.1611</td><td>0.0242</td><td>0.3477</td><td>0.0072</td><td>6.2358</td><td>1.9194</td><td></td></tr><tr><td>3000</td><td>1</td><td>55.64</td><td>0.1520</td><td>0.0091</td><td>0.3949</td><td>0.0096</td><td>6.2452</td><td>3.1421</td><td></td></tr><tr><td>1000</td><td>1</td><td>56.11</td><td>0.1482</td><td>0.0030</td><td>0.4221</td><td>0.0111</td><td>6.2688</td><td>5.0730</td><td></td></tr><tr><td>300</td><td>1</td><td>56.28</td><td>0.1468</td><td>0.0009</td><td>0.4199</td><td>0.0114</td><td>6.2535</td><td>8.1108</td><td></td></tr><tr><td>100</td><td>1</td><td>56.11</td><td>0.1482</td><td>0.0003</td><td>0.4380</td><td>0.0111</td><td>6.2621</td><td>12.0997</td><td></td></tr><tr><td>50</td><td>1</td><td>56.18</td><td>0.1476</td><td>0.0001</td><td>0.4400</td><td>0.0109</td><td>6.2639</td><td>15.5146</td><td></td></tr><tr><td>20000</td><td>2</td><td>73.92</td><td>0.0648</td><td>0.2074</td><td>0.1174</td><td>0.0063</td><td>7.1314</td><td>0.7926</td><td>*</td></tr><tr><td>15000</td><td>2</td><td>67.76</td><td>0.0842</td><td>0.2020</td><td>0.0964</td><td>0.0049</td><td>6.9420</td><td>0.8338</td><td>*</td></tr><tr><td>10000</td><td>2</td><td>61.64</td><td>0.1118</td><td>0.1789</td><td>0.0805</td><td>0.0040</td><td>6.6559</td><td>0.9663</td><td>*</td></tr><tr><td>7500</td><td>2</td><td>59.00</td><td>0.1275</td><td>0.1530</td><td>0.0807</td><td>0.0043</td><td>6.4846</td><td>1.1229</td><td>*</td></tr><tr><td>5000</td><td>2</td><td>57.13</td><td>0.1404</td><td>0.1123</td><td>0.0886</td><td>0.0050</td><td>6.3783</td><td>1.4177</td><td></td></tr><tr><td>4000</td><td>2</td><td>56.61</td><td>0.1443</td><td>0.0923</td><td>0.0971</td><td>0.0057</td><td>6.3390</td><td>1.6361</td><td></td></tr><tr><td>3000</td><td>2</td><td>56.31</td><td>0.1466</td><td>0.0704</td><td>0.1085</td><td>0.0067</td><td>6.3218</td><td>1.9785</td><td></td></tr><tr><td>1800</td><td>2</td><td>56.05</td><td>0.1487</td><td>0.0428</td><td>0.1222</td><td>0.0078</td><td>6.2599</td><td>2.7735</td><td></td></tr><tr><td>1000</td><td>2</td><td>56.18</td><td>0.1477</td><td>0.0236</td><td>0.1392</td><td>0.0095</td><td>6.2879</td><td>3.9275</td><td></td></tr><tr><td>500</td><td>2</td><td>56.23</td><td>0.1472</td><td>0.0118</td><td>0.1526</td><td>0.0104</td><td>6.2609</td><td>5.5744</td><td></td></tr><tr><td>100</td><td>2</td><td>56.30</td><td>0.1467</td><td>0.0023</td><td>0.1566</td><td>0.0114</td><td>6.2632</td><td>11.2290</td><td></td></tr><tr><td>50</td><td>2</td><td>56.31</td><td>0.1466</td><td>0.0012</td><td>0.1473</td><td>0.0114</td><td>6.2605</td><td>14.5222</td><td></td></tr><tr><td>20000</td><td>3</td><td>95.84</td><td>0.0297</td><td>0.3212</td><td>0.0589</td><td>0.0071</td><td>7.6393</td><td>0.6161</td><td>*</td></tr><tr><td>10000</td><td>3</td><td>76.81</td><td>0.0578</td><td>0.3120</td><td>0.0410</td><td>0.0048</td><td>7.2650</td><td>0.6705</td><td>*</td></tr><tr><td>5000</td><td>3</td><td>64.66</td><td>0.0968</td><td>0.2614</td><td>0.0329</td><td>0.0040</td><td>6.7890</td><td>0.9414</td><td>*</td></tr><tr><td>4000</td><td>3</td><td>62.21</td><td>0.1087</td><td>0.2348</td><td>0.0334</td><td>0.0040</td><td>6.5920</td><td>1.0718</td><td>*</td></tr><tr><td>3000</td><td>3</td><td>60.14</td><td>0.1204</td><td>0.1950</td><td>0.0382</td><td>0.0046</td><td>6.5102</td><td>1.3057</td><td></td></tr><tr><td>2300</td><td>3</td><td>58.75</td><td>0.1291</td><td>0.1604</td><td>0.0428</td><td>0.0054</td><td>6.4311</td><td>1.5429</td><td></td></tr><tr><td>1600</td><td>3</td><td>57.87</td><td>0.1351</td><td>0.1167</td><td>0.0524</td><td>0.0070</td><td>6.3742</td><td>2.0149</td><td></td></tr><tr><td>1000</td><td>3</td><td>57.07</td><td>0.1408</td><td>0.0761</td><td>0.0618</td><td>0.0080</td><td>6.3179</td><td>2.8246</td><td></td></tr><tr><td>500</td><td>3</td><td>56.64</td><td>0.1441</td><td>0.0389</td><td>0.0731</td><td>0.0095</td><td>6.2917</td><td>4.4527</td><td></td></tr><tr><td>100</td><td>3</td><td>56.37</td><td>0.1461</td><td>0.0079</td><td>0.0908</td><td>0.0110</td><td>6.2921</td><td>10.0294</td><td></td></tr><tr><td>50</td><td>3</td><td>56.35</td><td>0.1463</td><td>0.0040</td><td>0.0875</td><td>0.0112</td><td>6.2624</td><td>13.5295</td><td></td></tr><tr><td>10000</td><td>4</td><td>93.57</td><td>0.0320</td><td>0.4090</td><td>0.0200</td><td>0.0047</td><td>7.7135</td><td>0.5071</td><td>*</td></tr><tr><td>5000</td><td>4</td><td>75.96</td><td>0.0597</td><td>0.3823</td><td>0.0150</td><td>0.0036</td><td>7.2264</td><td>0.6233</td><td>*</td></tr><tr><td>3000</td><td>4</td><td>67.13</td><td>0.0865</td><td>0.3323</td><td>0.0145</td><td>0.0034</td><td>6.8100</td><td>0.8860</td><td>*</td></tr><tr><td>2000</td><td>4</td><td>62.90</td><td>0.1052</td><td>0.2694</td><td>0.0176</td><td>0.0045</td><td>6.6007</td><td>1.2091</td><td></td></tr><tr><td>1600</td><td>4</td><td>61.18</td><td>0.1143</td><td>0.2342</td><td>0.0205</td><td>0.0051</td><td>6.5100</td><td>1.3929</td><td></td></tr><tr><td>1000</td><td>4</td><td>59.15</td><td>0.1265</td><td>0.1619</td><td>0.0306</td><td>0.0067</td><td>6.3726</td><td>1.9771</td><td></td></tr><tr><td>700</td><td>4</td><td>58.12</td><td>0.1334</td><td>0.1195</td><td>0.0343</td><td>0.0076</td><td>6.3358</td><td>2.6185</td><td></td></tr><tr><td>500</td><td>4</td><td>57.45</td><td>0.1381</td><td>0.0884</td><td>0.0372</td><td>0.0085</td><td>6.3197</td><td>3.3626</td><td></td></tr><tr><td>300</td><td>4</td><td>57.21</td><td>0.1398</td><td>0.0537</td><td>0.0458</td><td>0.0099</td><td>6.3056</td><td>4.8379</td><td></td></tr><tr><td>100</td><td>4</td><td>56.56</td><td>0.1447</td><td>0.0185</td><td>0.0481</td><td>0.0110</td><td>6.2800</td><td>8.8842</td><td></td></tr><tr><td>2000</td><td>5</td><td>69.65</td><td>0.0775</td><td>0.3874</td><td>0.0071</td><td>0.0034</td><td>6.9309</td><td>0.8357</td><td>*</td></tr><tr><td>1600</td><td>5</td><td>66.72</td><td>0.0882</td><td>0.3526</td><td>0.0081</td><td>0.0040</td><td>6.7538</td><td>1.0106</td><td></td></tr><tr><td>1000</td><td>5</td><td>62.36</td><td>0.1080</td><td>0.2699</td><td>0.0130</td><td>0.0052</td><td>6.5122</td><td>1.4406</td><td></td></tr><tr><td>750</td><td>5</td><td>60.60</td><td>0.1176</td><td>0.2206</td><td>0.0165</td><td>0.0061</td><td>6.4122</td><td>1.7888</td><td></td></tr><tr><td>500</td><td>5</td><td>59.15</td><td>0.1265</td><td>0.1581</td><td>0.0222</td><td>0.0082</td><td>6.3587</td><td>2.5069</td><td></td></tr><tr><td>300</td><td>5</td><td>57.99</td><td>0.1342</td><td>0.1007</td><td>0.0273</td><td>0.0092</td><td>6.3333</td><td>3.8062</td><td></td></tr><tr><td>200</td><td>5</td><td>57.44</td><td>0.1381</td><td>0.0691</td><td>0.0301</td><td>0.0101</td><td>6.3126</td><td>5.1301</td><td></td></tr><tr><td>100</td><td>5</td><td>56.87</td><td>0.1423</td><td>0.0356</td><td>0.0329</td><td>0.0105</td><td>6.2973</td><td>7.7382</td><td></td></tr><tr><td>50</td><td>5</td><td>56.59</td><td>0.1444</td><td>0.0181</td><td>0.0354</td><td>0.0111</td><td>6.2820</td><td>11.2616</td><td></td></tr><tr><td>30</td><td>5</td><td>56.48</td><td>0.1453</td><td>0.0109</td><td>0.0383</td><td>0.0112</td><td>6.2879</td><td>14.4058</td><td></td></tr><tr><td>15</td><td>5</td><td>56.40</td><td>0.1459</td><td>0.0055</td><td>0.0389</td><td>0.0114</td><td>6.2622</td><td>19.9009</td><td></td></tr><tr><td>1000</td><td>7</td><td>72.71</td><td>0.0681</td><td>0.4672</td><td>0.0009</td><td>0.0029</td><td>6.9195</td><td>0.7881</td><td>*</td></tr><tr><td>750</td><td>7</td><td>68.71</td><td>0.0807</td><td>0.4153</td><td>0.0022</td><td>0.0039</td><td>6.7000</td><td>1.0633</td><td></td></tr><tr><td>500</td><td>7</td><td>64.51</td><td>0.0975</td><td>0.3345</td><td>0.0049</td><td>0.0054</td><td>6.5015</td><td>1.5040</td><td></td></tr><tr><td>300</td><td>7</td><td>61.26</td><td>0.1139</td><td>0.2344</td><td>0.0097</td><td>0.0073</td><td>6.3621</td><td>2.2690</td><td></td></tr><tr><td>200</td><td>7</td><td>59.65</td><td>0.1234</td><td>0.1692</td><td>0.0123</td><td>0.0088</td><td>6.3070</td><td>3.2795</td><td></td></tr><tr><td>100</td><td>7</td><td>58.00</td><td>0.1341</td><td>0.0920</td><td>0.0162</td><td>0.0101</td><td>6.2774</td><td>5.7758</td><td></td></tr><tr><td>50</td><td>7</td><td>57.17</td><td>0.1401</td><td>0.0481</td><td>0.0177</td><td>0.0107</td><td>6.2980</td><td>9.1238</td><td></td></tr><tr><td>30</td><td>7</td><td>56.83</td><td>0.1426</td><td>0.0294</td><td>0.0191</td><td>0.0109</td><td>6.2728</td><td>12.4074</td><td></td></tr><tr><td>15</td><td>7</td><td>56.57</td><td>0.1446</td><td>0.0149</td><td>0.0202</td><td>0.0115</td><td>6.2612</td><td>17.5854</td><td></td></tr></tbody></table>
[1] A. C. Balazs, T. Emrick, and T. P. Russell, Science **314**, 1107 (2006).

[2] K. I. Winey and R. A. Vaia, MRS Bull. **32**, 314 (2007).

[3] J. Jancar, J. Douglas, F. W. Starr, S. Kumar, P. Cas-
sagnau, A. Lesser, S. S. Sternstein, and M. Buehler,
Polymer **51**, 3321 (2010).

[4] H.-X. Zhou, G. Rivas, and A. P. Minton, Annu. Rev. Biophys. **37**, 375 (2008).

[5] K. A. Woodrow, Y. Cu, C. J. Booth, J. K. Saucier-Sawyer, M. J. Wood, and W. M. Saltzman, Nat. Mater. **8**, 526 (2009).

[6] H.-X. Zhou, FEBS Lett. **587**, 1053 (2013).

[7] K. S. Soppimath, T. M. Aminabhavi, A. R. Kulkarni,
and W. E. Rudzinski, J. Control. release **70**, 1 (2001).

[8] S. K. Lai, D. E. O'Hanlon, S. Harrold, S. T. Man, Y.-Y. Wang, R. Cone, and J. Hanes, Proc. Natl. Acad. Sci. **104**, 1482 (2007).

[9] I. Brigger, C. Dubernet, and P. Couvreur, Adv. Drug Deliv. Rev. **54**, 631 (2002).

[10] K. Cho, X. U. Wang, S. Nie, and D. M. Shin, Clin. Cancer Res. **14**, 1310 (2008).

[11] L.-H. Cai, S. Panyukov, and M. Rubinstein, Macromolecules **44**, 7853 (2011).

[12] S. A. Egorov, J. Chem. Phys. **134**, 84903 (2011).

[13] U. Yamamoto and K. S. Schweizer, J. Chem. Phys. **135**, 224902 (2011).

[14] U. Yamamoto and K. S. Schweizer, Macromolecules **48**, 152 (2014).

[15] Y. Dong, X. Feng, N. Zhao, and Z. Hou, J. Chem. Phys. **143**, 24903 (2015).

[16] A. Tuteja, M. E. Mackay, S. Narayanan, S. Asokan, and M. S. Wong, Nano Lett. **7**, 1276 (2007).

[17] C. A. Grabowski, B. Adhikary, and A. Mukhopadhyay, Appl. Phys. Lett. **94**, 21903 (2009).

[18] I. Kohli and A. Mukhopadhyay, Macromolecules **45**, 6143 (2012).

[19] S. Gam, J. S. Meth, S. G. Zane, C. Chi, B. A. Wood, M. E. Seitz, K. I. Winey, N. Clarke, and R. J. Composto, Macromolecules **44**, 3494 (2011).

[20] S. Gam, J. S. Meth, S. G. Zane, C. Chi, B. A. Wood, K. I. Winey, N. Clarke, and R. J. Composto, Soft Matter **8**, 6512 (2012).

[21] C. Lin, S. Gam, J. S. Meth, N. Clarke, K. I. Winey, and R. J. Composto, Macromolecules **46**, 4502 (2013).

[22] J. Choi, M. J. A. Hore, J. S. Meth, N. Clarke, K. I. Winey, and R. J. Composto, ACS Macro Lett. **2**, 485 (2013).

[23] S. Y. Kim, H. W. Meyer, K. Saalwächter, and C. F. Zukoski, Macromolecules **45**, 4225 (2012).

[24] R. Poling-Skutvik, R. Krishnamoorti, and J. C. Conrad, ACS Macro Lett. **4**, 1169 (2015).

[25] F. Babaye Khorasani, R. Poling-Skutvik, R. Krish-
namoorti, and J. C. Conrad, Macromolecules **47**, 5328
(2014).

[26] D. Bedrov, G. D. Smith, and J. S. Smith, J. Chem. Phys. **119**, 10438 (2003).

[27] J. Liu, D. Cao, and L. Zhang, J. Phys. Chem. C **112**, 6653 (2008).

[28] J. T. Kalathi, U. Yamamoto, K. S. Schweizer, G. S. Grest, and S. K. Kumar, Phys. Rev. Lett. **112**, 108301 (2014).

[29] A. Patti, J. Phys. Chem. B **118**, 3731 (2014).

[30] Y. Li, M. Kröger, and W. K. Liu, Soft Matter **10**, 1723 (2014).

[31] I. V. Volgin, S. V. Larin, E. Abad, and S. V. Lyulin, Macromolecules **50**, 2207 (2017).

[32] A. Karatrantos, R. J. Composto, K. I. Winey, and N. Clarke, J. Chem. Phys. **146**, 203331 (2017).

[33] M. Rubinstein and R. H. Colby, *Polymer physics* (Oxford University Press New York, 2003).

[34] I. Teraoka, *Polymer Solutions: An Introduction to Physical Properties* (John Wiley & Sons, Inc., 2002).

[35] D. S. Banks and C. Fradin, Biophys. J. **89**, 2960 (2005).

[36] R. Metzler, J.-H. Jeon, and A. G. Cherstvy, Biochim. Biophys. Acta (BBA)-Biomembranes **1858**, 2451 (2016).

[37] F. B. Wyart and P.-G. de Gennes, Eur. Phys. J. E **1**, 93 (2000).

[38] F. Ould-Kaddour and D. Levesque, Phys. Rev. E **63**, 11205 (2000).

[39] T. Desai, P. Keblinski, and S. K. Kumar, J. Chem. Phys. **122**, 134910 (2005).

[40] J. Liu, Y. Gao, D. Cao, L. Zhang, and Z. Guo, Langmuir **27**, 7926 (2011).

[41] J. B. Hooper and K. S. Schweizer, Macromolecules **39**, 5133 (2006).

[42] D. Meng, S. K. Kumar, S. Cheng, and G. S. Grest, Soft Matter **9**, 5417 (2013).

[43] P. G. Bolhuis, E. J. Meijer, and A. A. Louis, Phys. Rev. Lett. **90**, 68304 (2003).

[44] S. Asakura and F. Oosawa, J. Chem. Phys. **22**, 1255 (1954).

[45] W. C. Koon, J. Phys. Condens. Matter **14**, R859 (2002).

[46] E. J. Meijer and D. Frenkel, J. Chem. Phys. **100**, 6873 (1994).

[47] M. Dijkstra, J. M. Brader, and R. Evans, J. Phys. Condens. Matter **11**, 10079 (1999).

[48] K. Kremer and G. S. Grest, J. Chem. Phys. **92**, 5057 (1990).

[49] J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. **54**, 5237 (1971).

[50] S. Gong, Q. Chen, J. F. Moll, S. K. Kumar, and R. H. Colby, ACS Macro Lett. **3**, 773 (2014).

[51] A. L. Frischknecht, E. S. McGarrity, and M. E. Mackay, J. Chem. Phys. **132**, 204901 (2010).

[52] A. Karatrantos, N. Clarke, R. J. Composto, and K. I. Winey, Soft Matter **11**, 382 (2015).

[53] R. S. Hoy, K. Foteinopoulou, and M. Kröger, Phys. Rev. E **80**, 31803 (2009).

[54] L. Fetters, D. Lohse, D. Richter, T. Witten, and A. Zirkel, Macromolecules **27**, 4639 (1994).

[55] L. J. Fetters, D. J. Lohse, S. T. Milner, and W. W. Graessley, Macromolecules **32**, 6847 (1999).

[56] M. Kröger and S. Hess, Phys. Rev. Lett. **85**, 1128 (2000).

[57] In our simulations $N = 100$ and $\rho_m < 0.32$, therefore a crude estimate gives $N_e \gtrsim 85 \cdot (0.85/0.32)^2 = 600$.

[58] S. Plimpton, J. Comput. Phys. **117**, 1 (1995).

[59] R. Auhl, R. Everaers, G. S. Grest, K. Kremer, and S. J. Plimpton, J. Chem. Phys. **119**, 12718 (2003).

[60] A. Chen, N. Zhao, and Z. Hou, Soft Matter **13**, 8625 (2017).

[61] M. Tuckerman, *Statistical mechanics: theory and molecular simulation* (Oxford University Press, 2010).
[62] J.-P. Hansen and I. R. McDonald, *Theory of simple liquids* (Elsevier, 1990).

[63] L. Verlet, Phys. Rev. **165**, 201 (1968).

[64] S. Torquato, T. M. Truskett, and P. G. Debenedetti, Phys. Rev. Lett. **84**, 2064 (2000).

[65] S. Torquato, *Random heterogeneous materials: microstructure and macroscopic properties* (Springer Science & Business Media, 2013).

[66] V. Testard, L. Berthier, and W. Kob, J. Chem. Phys. **140**, 164502 (2014).

[67] The factor $1/\sqrt{3}$ is justified by considering that when $\rho_m$ is not too high the mesh size can be extracted from the low-$q$ limit of the monomer-monomer structure factor, which is given by the Ornstein-Zernike relation $S_{mm}(q) = S_{mm}(0)/[1 + (\xi q)^2]$, and that at low monomer densities and low wavevectors $S_{mm}(q) \approx S_1(q) = N/[1 + (R_g^*q/\sqrt{3})^2]$, where $S_1(q)$ is the single chain monomer-monomer structure factor [33]. From the previous relations, it follows that for $\rho_m < \rho_m^*$ we have $\xi = R_g^*/\sqrt{3}$.

[68] A. I. Nakatani, W. Chen, R. G. Schmidt, G. V. Gordon, and C. C. Han, Polymer (Guildf). **42**, 3713 (2001).

[69] M. E. Mackay, A. Tuteja, P. M. Duxbury, C. J. Hawker, B. Van Horn, Z. Guan, G. Chen, and R. S. Krishnan, Science (80-.). **311**, 1740 (2006).

[70] A. Tuteja, P. M. Duxbury, and M. E. Mackay, Phys. Rev. Lett. **100**, 77801 (2008).

[71] A. Karatrantos, N. Clarke, and M. Kröger, Polym. Rev. **56**, 385 (2016).

[72] As a first approximation, if $\phi_N$ is low, the excluded volume can be estimated as $V_{ex} = \pi N_N (\sigma_N + 1)^3 / 6V = \phi_N(1 + \sigma_N^{-3} + 3\sigma_N^{-2} + 3\sigma_N^{-1})$.

[73] G. D. Smith, W. Paul, M. Monkenbusch, and D. Richter, J. Chem. Phys. **114**, 4285 (2001).

[74] W. Paul, G. Smith, D. Y. Yoon, B. Farago, S. Rathgeber, A. Zirkel, L. Willner, and D. Richter, Phys. Rev. Lett. **80**, 2346 (1998).

[75] J. Padding and W. J. Briels, J. Chem. Phys **114**, 8685 (2001).

[76] W. Paul, K. Binder, D. W. Heermann, and K. Kremer, J. Chem. Phys. **95**, 7726 (1991).

[77] A. Kopf, B. Dünweg, and W. Paul, J. Chem. Phys. **107**, 6945 (1997).

[78] T. Ge, J. T. Kalathi, J. D. Halverson, G. S. Grest, and M. Rubinstein, Macromolecules **50**, 1749 (2017).

[79] R. Chen, R. Poling-Skutvik, A. Nikoubashman, M. P. Howard, J. C. Conrad, and J. C. Palmer, Macromolecules **51**, 1865 (2018).

[80] W. Kob and H. C. Andersen, Phys. Rev. E **51**, 4626 (1995).

[81] B. Dünweg and K. Kremer, J. Chem. Phys. **99**, 6983 (1993).

[82] I.-C. Yeh and G. Hummer, J. Phys. Chem. B **108**, 15873 (2004).

[83] A. Cavagna, Phys. Rep. **476**, 51 (2009).

[84] K. Binder and W. Kob, *Glassy materials and disordered solids: An introduction to their statistical mechanics* (World Scientific, 2011).

[85] L. Berthier and G. Biroli, Rev. Mod. Phys. **83**, 587 (2011).

[86] G. Adam and J. H. Gibbs, J. Chem. Phys **43**, 139 (1965).

[87] B. U. Felderhof, J. Phys. A. Math. Gen. **11**, 929 (1978).

[88] L. Bocquet, J.-P. Hansen, and J. Piasecki, J. Stat. Phys. **76**, 527 (1994).

[89] J. R. Schmidt and J. L. Skinner, J. Chem. Phys. **119**, 8062 (2003).

[90] J. T. Hynes, R. Kapral, and M. Weinberg, J. Chem. Phys. **70**, 1456 (1979).

[91] T. C. Hales, Ann. Math., 1065 (2005).

[92] S. Wu, Polymer **26**, 1855 (1985).

[93] Y. Jiao and S. Torquato, Phys. Rev. E **84**, 041309 (2011).

[94] M. Rintoul, S. Torquato, C. Yeong, D. Keane, S. Erramilli, Y. Jun, D. Dabbs, and I. Aksay, Phy. Rev. E **54**, 2663 (1996).

[95] W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graph. **14**, 33 (1996).

[96] J. Horbach and W. Kob, Phys. Rev. B **60**, 3169 (1999).