# Determining the mesh size of polymer solutions via the pore size distribution

Valerio Sorichetti, Virginie Hugouvieux, Walter Kob

► To cite this version:

Valerio Sorichetti, Virginie Hugouvieux, Walter Kob. Determining the mesh size of poly-mer solutions via the pore size distribution. Macromolecules, 2020, 53 (7), pp.2568-2581. 10.1021/acs.macromol.9b02166 . hal-02540128

HAL Id: hal-02540128

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

Submitted on 10 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Characterizing the mesh size of polymer solutions via the pore size distribution

Valerio Sorichetti

Laboratoire Charles Coulomb (L2C), Univ. Montpellier, CNRS,
F-34095, Montpellier, France and IATE, INRA, CIRAD,
Montpellier SupAgro, Univ. Montpellier, F-34060, Montpellier, France

Virginie Hugouvieux

IATE, INRA, CIRAD, Montpellier SupAgro, Univ. Montpellier, F-34060, Montpellier, France

Walter Kob

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

(Dated: August 6, 2019)

In order to characterize the geometrical mesh size $\xi$, we simulate a solution of coarse-grained polymers with densities ranging from the dilute to the concentrated regime and for different chain lengths. Conventional ways to estimate $\xi$ rely either on scaling assumptions which give $\xi$ only up to an unknown multiplicative factor, or on measurements of the monomer density fluctuation correlation length $\xi_c$. We determine $\xi_c$ from the monomer structure factor and from the radial distribution function, and find that the identification $\xi = \xi_c$ is not justified outside of the semidilute regime. In order to better characterize $\xi$, we compute the pore size distribution (PSD) following two different definitions, one by Torquato *et al.* [1] and one by Gubbins *et al.* [2]. We show that the mean values of the two distributions, $\langle r \rangle_T$ and $\langle r \rangle_G$, both display the behavior predicted for $\xi$ by scaling theory, and argue that $\xi$ can be identified with either one of these quantities. This identification allows to interpret the PSD as the distribution of mesh sizes, a quantity which conventional methods cannot access. Finally, we show that it is possible to map a polymer solution on a system of hard or overlapping spheres, for which Torquato's PSD can be computed analytically and reproduces accurately the PSD of the solution. We give an expression that allows $\langle r \rangle_T$ to be estimated with great accuracy in the semidilute regime by knowing only the radius of gyration and the density of the polymers.

## I. INTRODUCTION

Polymer solutions are outstanding liquids showing a rich behavior due to the exceptional variety of their structural and dynamical properties. While simple atomic and molecular liquids generally have only one relevant length scale –the size of the atoms/molecules– a polymer solution has many: the monomer size, the Kuhn length, the entanglement length, all the way up to the chain size [3, 4]. In semidilute polymer solutions, i.e., solutions which are dense enough that different chains start to overlap, the only physically relevant length scale is (in the absence of entanglement) the *mesh size*, $\xi$. The general idea behind the concept of mesh size is that a dense polymer solution observed at a certain instant in time looks very similar to an intricate network with a certain mesh size $\xi$ [3].

The relevance of the mesh size is most apparent when considering the problem of the diffusion of proteins or nanoparticles in polymer liquids, which has received a vast amount of attention in recent years because of its applications to biology [5, 6] and medicine, e.g. for drug delivery [7–9]. Many theoretical approaches have been proposed to describe the diffusion of particles in polymer solutions and gels, such as geometric obstruction models [10–12], hydrodynamic models [13–18] and mode-coupling theory [19–21] (for a review of some of these models, see Ref. 22). Although the predictions of differ-

ent models can vary widely, most of them agree upon the fact that $\xi$ is central in controlling the diffusion of the particles. In the last few years, a model based on scaling theory has been proposed to describe the diffusion of nonsticky spherical nanoparticles in polymer liquids [23]. This model, which seems to be the one in best agreement with recent experiments [24–26] and simulations [27, 28], predicts that particles with radius $R < \xi$ can slip through the mesh, experiencing the viscosity $\eta_s$ of the pure solvent, while larger particles will experience a higher effective viscosity. The mesh size represents therefore the length scale of the transition from free diffusion to obstructed diffusion in a polymer solution. Naturally, the mesh size is also a fundamental quantity when considering the diffusion of particles in polymer networks and gels [22, 29–34].

When considering the diffusion of a probe particle in a polymer solution, it is clear that the relevant mesh size is, to a first approximation, just the average geometrical size of the pores, or holes, in the solution, i.e., those regions which are filled only with solvent. In view of applications to diffusion problems, we will in the present work adopt this intuitive concept of $\xi$ as the average pore size of the polymer solution. This also allows to extend the intuitive definition of $\xi$ proposed by de Gennes, as the average size of the transient mesh in a semidilute solution [3], to the dilute and concentrated regimes: In the dilute regime the chains do not overlap, but a mesh size still exists as an average geometric distance between neighboring chains.
In the concentrated regime, the chains overlap so strongly that there is no clear mesh structure, but we can still measure the size of the pores.

We are, however, immediately confronted with the question: How can one determine $\xi$? Normally, two ways are used to estimate the mesh size: scaling theory [3], or measurements of the monomer density fluctuation correlation length, $\xi_c$ [3, 4, 35]. Scaling calculations are useful and allow to obtain easily an estimate of $\xi$. However, by construction they can only give the quantity of interest up to an unknown multiplicative factor, which is very inconvenient in cases where a precise knowledge of $\xi$ is required. Estimates based on the correlation length –which can be obtained from scattering experiments– are also useful; however, as we will discuss below, the correlation length is not the same quantity as the mesh size, and taking them to be the same quantity can lead to apparently nonsensical results, such as $\xi$ increasing when the polymer concentration is increased [35].

In the present paper, we use molecular dynamics simulations of coarse-grained polymers to perform a systematic comparison of the different conventional ways to measure $\xi$, i.e., those based on scaling theory and those based on static correlation functions in Fourier and real space. We show that these techniques have some limitations, can lead to apparently contradictory behaviors and provide at best an approximate value for the average mesh size in the semidilute regime. To overcome these problems, we propose a different method to estimate $\xi$, based on the concept of pore size distribution [1, 2, 36–43]. If the coordinates of the monomers are known, this method allows to measure $\xi$ reliably at any density, and to obtain not only an average value, but also the distribution of mesh sizes. This last feature can be particularly important, since it is known that particles diffusing in systems with the same average mesh size can display completely different dynamical behaviors [33], i.e., the heterogeneity of the polymer medium plays a relevant role in particle diffusion. Although we consider here only the case of a solution of linear chains, the proposed method is very general and can be applied to systems of polymers with different topologies (e.g. rings, randomly branched, dendrimers) and also to networks and gels.

The paper is organized as follows: In Sec. II we review the theoretical background, discussing the concept of blob, of density fluctuation correlation length and of geometrical mesh size. In Sec. III, the model and the details of the simulations are presented. Sec. IV deals with the main results: Measurements of the correlation length and scaling estimates of $\xi$ are compared with results from the pore size distribution. We conclude with a summary in Sec. V.

## II. THEORETICAL BACKGROUND

Before discussing the approach that we propose to characterize the geometrical mesh size, we review the

<figure><img src="image_2.png" /><figcaption>FIG. 1. Schematic representation of a polymer solution in the dilute (a) and semidilute (b) regime. Circles denote the blobs, of size $\xi_b$. The other relevant length scales are the geometrical mesh size, $\xi$, and the polymer correlation length, $\xi_c$. In the semidilute regime, these three length scales are proportional through a $O(1)$ dimensionless factor. In the dilute regime, $\xi$ is very different from the other two length scales.</figcaption></figure>

concepts of blob [3, 4], monomer density fluctuation correlation length [3, 4, 35] and geometrical mesh size. Most of what is discussed in this section can be found in standard texts [3, 4, 44, 45]. However, we feel it is useful to summarize these results, since they will be extensively used in what follows. Moreover, this section also serves to introduce the terminology and notations which will be used in the rest of this work.

To simplify the discussion, we limit ourselves to the case of polymers in athermal solvent [4]. This choice is also motivated by the fact that the model chosen for the simulations is athermal, in the sense that all the non-bonded interactions are purely repulsive, as discussed in Sec. III below. Since our discussion is limited to the athermal solvent case, we use the symbol $\nu$ exclusively to denote the Flory exponent in good/athermal solvent, i.e., $\nu \approx 0.588$ [4, 46].

### A. The concept of blob

The structure of a polymer solution can be rather complicated, but its main features can easily be obtained from the blob model, originally introduced by de Gennes [3]. A blob is defined as a region of space of size $\xi_b$ which contains almost exclusively monomers from a single chain, or equivalently inside of which inter-chain interactions are negligible [3, 4]. A graphical representation of a blob is given in Fig. 1. The above definition implies that in a semidilute solution, where the polymer chains overlap, the blobs are space-filling, i.e., if the number of monomers in a blob is $N_b$ we have [4]

$$ \rho \approx \frac{N_b}{\xi_b^3}, \qquad (1) $$
where $\rho$ is the bulk monomer density, see Fig. 1b. (In Eq. (1), and also in the following, we use the symbol $\approx$ to indicate that two quantities are proportional through a dimensionless factor of order 1, while we will use the symbol $\simeq$ when two quantities can be considered approximately equal.)

On length scales $r \approx \xi_b$ or smaller, the chain is by definition mostly surrounded by solvent, and therefore it assumes a swollen configuration [3, 4]:

$$ \xi_b \approx bN_b^\nu, \quad (2) $$

where $b$ is the effective bond length. From Eqs. (1) and (2) we obtain

$$ \xi_b \approx b(b^3\rho)^{-\nu/(3\nu-1)}. \quad (3) $$

We note that, since $\nu \simeq 0.588$, $-\nu/(3\nu - 1) \simeq -0.770$.

On length scales $r > \xi_b$, the excluded volume interactions are screened by the surrounding chains and therefore the structure of the chain is ideal [3, 4]. If the chains are ideal at all length scales $r > b$, the solution is called concentrated and $\xi_b \approx b$ [4, 45] [47]. In the opposite limit we have a dilute solution, where chains do not overlap and therefore $\xi_b \approx R_{g0}$, where $R_{g0} \equiv \lim_{\rho \to 0} R_g$ (Fig. 1a). Between these two regimes, we find the semidilute behavior given by Eq. (3). We can thus summarize the scaling behavior of $\xi_b$ as follows:

$$ \xi_b \approx \begin{cases} R_{g0} & \rho < \rho^* \\ b(b^3\rho)^{-\nu/(3\nu-1)} & \rho^* < \rho < \rho^{**} \\ b & \rho > \rho^{**}, \end{cases} \quad (4) $$

where $\rho^*$ (the overlap concentration) marks the onset of the semidilute regime and $\rho^{**}$ the one of the concentrated regime.

Imposing continuity between the three relations in Eq. (4), and recalling that $R_{g0} \approx bN^\nu$, where $N$ is the degree of polymerization [3, 4], we obtain

$$ \rho^* \approx N R_{g0}^{-3} \approx N^{1-3\nu} b^{-3} \quad (5) $$

and

$$ \rho^{**} \approx b^{-3} \approx \rho_{\text{melt}}. \quad (6) $$

We note that $\lim_{N \to \infty} \rho^* = 0$, i.e., infinitely long chains are never in the dilute regime.

Exploiting the relation $R_{g0} \approx bN^\nu$ and introducing the dimensionless scaling variable $\chi$, defined as [48, 49] [50]

$$ \chi \equiv N(b^3\rho)^{1/(3\nu-1)} \approx \left(\frac{\rho}{\rho^*}\right)^{1/(3\nu-1)}, \quad (7) $$

Eq. (4) can be rewritten in the following dimensionless form:

$$ \frac{\xi_b}{R_{g0}} \approx \begin{cases} 1 & \chi < \chi^* \\ \chi^{-\nu} & \chi^* < \chi < \chi^{**} \\ N^{-\nu} & \chi > \chi^{**}. \end{cases} \quad (8) $$

It is easy to verify using Eqs. (5)-(7), that $\chi^* \equiv \chi(\rho^*) = 1$ and $\chi^{**} \equiv \chi(\rho^{**}) = N$. We note that $\chi^*$ is independent of $N$, whereas $\chi^{**}$ depends on $N$. Conversely, $\rho^*$ depends on $N$, whereas $\rho^{**}$ is independent of $N$.

The concept of blob is extremely useful for scaling calculations; however, there is no quantity which can be measured in experiments or simulations which unambiguously corresponds to the blob size $\xi_b$. We will see below that it is possible to introduce two different quantities, the density fluctuation correlation length and the geometrical mesh size, which have the same scaling behavior as $\xi_b$ in the semidilute regime and are accessible from simulations or experimental data.

## B. Monomer density fluctuation correlation length

Due to chain connectivity, the local monomer density in a polymer has significant spatial fluctuations. One possibility to quantify the spatial extent of these fluctuations is to study the radial distribution function $g(r)$ [51], defined in such a way that $4\pi\rho r^2 g(r) dr$ represents the probability to find a monomer at distance $r < x < r + dr$ from a randomly chosen monomer. At low density, the large $r$ behavior of $g(r)$ is expected to have the form [3, 44, 52]

$$ g(r) = 1 + \frac{A}{r} e^{-r/\xi_c}, \quad (9) $$

where $A > 0$ and $\xi_c > 0$ are constants. In Eq. (9), $\xi_c$ has clearly the role of a correlation length, since it represents the typical decay length scale of the density fluctuations. Therefore $\xi_c$ is known as the monomer density fluctuation correlation length (or simply correlation length). At higher density, Eq. (9) must be modified by the introduction of a periodic modulation [35, 53–55],

$$ g(r) = 1 + \frac{B}{r} e^{-r/\xi_c} \sin\left(\frac{2\pi r}{\lambda} + \theta\right), \quad (10) $$

where $B > 0$, $\lambda > 0$, and $\theta$ are constants.

The radial distribution function is easily accessible in simulations; in experiments, however, it is easier to study the static structure factor $S(q)$ [4, 51], which is related to $g(r)$ through a Fourier transform [51]:

$$ S(q) = 1 + 4\pi\rho \int_{0}^{\infty} [g(r) - 1] \frac{\sin(qr)}{qr} r^2 dr. \quad (11) $$
From Eqs. (9) and (11), we obtain the low $q$ behavior of $S(q)$, which is given by the Ornstein-Zernike function [3, 44, 52] [56],

$$S(q) = \frac{S(0)}{1 + (q\xi_c)^2}, \quad (12)$$

where $S(0) \equiv \lim_{q \to 0} S(q)$.

In the dilute regime $\rho < \rho^*$, the monomer structure factor coincides approximately with the structure factor of a single chain: $S(q) \simeq S_1(q)$ [4, 45]. Since we know that for $q < 2\pi/R_g$ [4, 45]

$$S_1(q) \approx N \left[ 1 - \left( \frac{qR_g}{\sqrt{3}} \right)^2 \right] \approx \frac{N}{1 + (qR_g/\sqrt{3})^2}, \quad (13)$$

it follows from Eq. (12) that for $\rho < \rho^*$, $\xi_c \simeq R_{g0}/\sqrt{3}$. Therefore, in a dilute solution $\xi_c \approx \xi_b \approx R_{g0}$.

Since the structure of a semidilute solution of short chains is identical to the one of long chains, we expect that in the semidilute regime $\xi_c$ is independent of $N$. Using scaling arguments, it is then straightforward to prove that [4, 45]

$$\xi_c \approx b(b^3 \rho)^{-\nu/(3\nu-1)} \quad (\rho^* < \rho < \rho^{**}), \quad (14)$$

which is equivalent to (3). We can therefore conclude that $\xi_c \approx \xi_b$ in the whole range $\rho < \rho^{**}$.

The relation (14) loses its validity for $\rho > \rho^{**}$. The reason is that in this density regime chain connectivity becomes irrelevant, and the structure of the system is dominated by local packing constraints. The correlation length behaves therefore like that of a dense molecular liquid, i.e., it increases with increasing density[35]. Therefore one can summarize the behavior of $\xi_c$ as follows:

$$\xi_c \approx \begin{cases} R_{g0} & \rho < \rho^* \\ b(b^3 \rho)^{-\nu/(3\nu-1)} & \rho^* < \rho < \rho^{**} \\ h(\rho) & \rho > \rho^{**}, \end{cases} \quad (15)$$

where $h(\rho)$ is an increasing function of $\rho$ for $\rho > \rho^{**}$.

### C. Geometrical mesh size

In the present work we are interested in measuring the geometrical mesh size (or simply “mesh size”) $\xi$, which we define as the average size of the pores in the system. Intuitively, in the semidilute regime $\xi$ will be very close to $\xi_b$, and therefore to $\xi_c$ (Eq. (14) and Fig. 1b):

$$\xi \approx \xi_b \approx \xi_c \quad (\rho^* < \rho < \rho^{**}). \quad (16)$$

<figure><img src="image_3.png" /><figcaption>FIG. 2. Expected $\rho$-dependence of the blob size $\xi_b$, the correlation length $\xi_c$ and the geometrical mesh size $\xi$, Eqs. (4), (15), and (17).</figcaption></figure>

However, if we want to maintain our definition of $\xi$ as the average size of the pores, then it is clear that Eq. (16) cannot be extended outside of the density range $\rho^* < \rho < \rho^{**}$. Indeed, in the dilute regime $\rho < \rho^*$ the mesh size is nothing else than the average distance between neighboring chains, i.e., $\xi \approx (\rho/N)^{-1/3}$, and becomes infinite in the limit $\rho \to 0$, whereas $\xi_b \approx \xi_c \approx R_{g0}$. Moreover, in the concentrated regime $\rho > \rho^{**}$ the size of the pores becomes vanishingly small, $\xi \to 0$, while $\xi_b \approx b$, and $\xi_c$ increases. To summarize, we expect the following behavior for the geometrical mesh size $\xi$ as a function of $\rho$:

$$\xi \approx \begin{cases} (\rho/N)^{-1/3} & \rho < \rho^* \\ b(b^3 \rho)^{-\nu/(3\nu-1)} & \rho^* < \rho < \rho^{**} \\ f(\rho) & \rho > \rho^{**}, \end{cases} \quad (17)$$

where $f(\rho)$ is a decreasing function of $\rho$ for $\rho > \rho^{**}$. In terms of the scaling variable $\chi$,

$$\frac{\xi}{R_{g0}} \approx \begin{cases} \chi^{-\nu+1/3} & \chi < \chi^* \\ \chi^{-\nu} & \chi^* < \chi < \chi^{**} \\ f[\rho(\chi)] & \chi > \chi^{**}. \end{cases} \quad (18)$$

The scaling behaviors of $\xi_b$, $\xi_c$ and $\xi$ as a function of $\rho$ are schematically represented in Fig. 2. We mention that $\xi$ is often estimated using the relation (see Eqs. (3) and (5))

$$\xi \approx R_{g0} \left( \frac{\rho}{\rho^*} \right)^{-\nu/(3\nu-1)}, \quad (19)$$
where [44, 45]

$$ \rho^* \approx \frac{3N}{4\pi R_{g0}^3}. \qquad (20) $$

Although Eq. (20) is a reasonable estimate, in reality the transition between the dilute and the semidilute regimes is far from being sharp (see Sec. IV), meaning that the uncertainty on $\xi$ estimated from Eqs. (19)-(20) can be quite large.

### III. MODEL AND SIMULATION METHOD

We have performed NVT molecular dynamics simulations (MD) of a system of $N_c$ coarse-grained polymer chains of length (degree of polymerization) $N$ at different monomer densities $\rho = NN_c/V = M/V$. The model chosen to simulate the chains is the bead-spring model of Kremer and Grest [57]. The excluded volume interactions between all monomers are modeled as a purely repulsive Weeks-Chandler-Andersen (WCA) potential [58],

$$ U(r) = \begin{cases} 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 + \frac{1}{4} \right] & r \le 2^{1/6}\sigma \\ 0 & \text{otherwise.} \end{cases} \qquad (21) $$

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

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

where $k = 30\epsilon/\sigma^2$ and $r_0 = 1.5\sigma$. The combined effect of the FENE and the WCA potentials prevents the chains from crossing each other even at high density at temperature $T = 1.0$ [57]. Since non-bonded interactions are purely repulsive, this model mimics the behavior of polymers in an athermal solvent [4, 35].

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

We considered chain lengths $N = 50, 200$, and $1000$. For $N = 50$ and $N = 200$, we consider $N_c = 200$ and $N_c = 50$ chains, respectively, so that the total number of monomers is $M = N_cN = 10^4$. For $N = 1000$, we simulated two systems, one of 10 chains ($M = 10^4$) and one of 50 chains ($M = 5 \cdot 10^4$), in order to check for the presence of finite-size effects. For the system of 10 chains, five independent simulations were performed in order to improve the statistics. For none of the studied quantities any relevant difference was found between the systems with $M = 10^4$ and those with $M = 5 \cdot 10^4$.

The range of monomer densities $\rho$ for the different systems are reported in Tab. I. While many of these densities are meant to represent a polymer solution, we do not take into account hydrodynamic interactions between the monomers, since we focus on static quantities, which are unaffected by hydrodynamic interactions.

All the simulations were carried out using the LAMMPS software [59]. The temperature was kept constant at $T = 1.0$ by means of a Langevin thermostat [60], so that the force experienced by a monomer is

$$ m\ddot{r} = -\nabla U_{\text{tot}}(\{\mathbf{r}\}) - m\Gamma\dot{r} + \sqrt{2m\Gamma k_B T} \zeta(t), \qquad (23) $$

where $\mathbf{r}$ is the position vector, and $U_{\text{tot}}(\{\mathbf{r}\}) = U_{\text{tot}}(\mathbf{r}_1, \dots, \mathbf{r}_M)$ is the total interaction potential acting on the monomer. The second term on the right side of Eq. (23) represents a viscous damping, with $\Gamma = 0.1$ the friction coefficient. The last term is a random, uncorrelated force mimicking the collisions with “solvent” particles.

The simulation box is cubic and periodic boundary conditions are applied in all directions. In all the simulations, the MD integration time step is $\delta t = 3 \cdot 10^{-3}$. The initial configurations are prepared by randomly placing the polymers in the box; for every set of values ($N, N_c, \rho$) a different initial configuration is created. Initially, overlaps between the monomers 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 [61]). After the overlaps have been removed, we switch to the WCA potential, Eq. (21), and perform an equilibration run of duration $t_e$ (corresponding to $t_e/\delta t$ MD integration steps) before starting the production run. Since the relaxation time of a polymer chain scales as $N^2$ for unentangled chains and as $N^{3.4}$ for entangled chains [4], longer equilibration times were chosen for the systems with longer chains. For each system, we checked that the equilibration time $t_e$ was larger than the longest relaxation time of the system (see Sec. SI in the S.I.). The values of $t_e$ for the different systems are reported in Tab. I.

<table><caption>TABLE I. Details of the simulated systems: chain length $N$, number of chains $N_c$, number of monomers $M$, monomer density $\rho$, equilibration time $t_e$.</caption><thead><tr><th>N</th><th>N<sub>c</sub></th><th>M = N<sub>c</sub>N</th><th>ρ (range)</th><th>t<sub>e</sub></th></tr></thead><tbody><tr><td>50</td><td>200</td><td>1 · 10<sup>4</sup></td><td>[0.001, 1.00]</td><td>1.5 · 10<sup>5</sup></td></tr><tr><td>200</td><td>50</td><td>1 · 10<sup>4</sup></td><td>[0.001, 0.85]</td><td>4.5 · 10<sup>5</sup></td></tr><tr><td>1000</td><td>10</td><td>1 · 10<sup>4</sup></td><td>[0.001, 0.26]</td><td>3.0 · 10<sup>6</sup></td></tr><tr><td>1000</td><td>50</td><td>5 · 10<sup>4</sup></td><td>[0.001, 0.26]</td><td>3.0 · 10<sup>6</sup></td></tr></tbody></table>

### IV. RESULTS AND DISCUSSION

As a preliminary characterization of the structural properties of the simulated systems, we have analyzed
the distribution of bond angles and the radius of gyration $R_g$ as a function of density. A complete discussion can be found in Sec. AI of the Appendix. The main result is that the transition from the semidilute to the concentrated regime occurs at density $\rho \approx 0.3-0.4$. In order to simplify the discussion, we will in the following operatively set the onset of the concentrated regime at $\rho^{**} = 0.3$, which implies that $\chi^{**} \equiv \chi(\rho^{**}) = 0.3N$ (Eq. (7)).

In the following subsection, we perform a systematic study of the correlation length, obtained alternatively from the structure factor and from the radial distribution function, as a function of density.

### A. Density fluctuation correlation length

As discussed in Sec. II, the correlation length $\xi_c$ can be measured either from the monomer structure factor $S(q)$ or from the monomer radial distribution function $g(r)$. The monomer structure factor is defined as [51]

$$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 (24)$$

where $M$ is the number of monomers and $\langle \cdot \rangle$ denotes the thermodynamic average. Since the system is isotropic, we will consider the spherically averaged structure factor, $S(q)$, which depends only on $q \equiv |\mathbf{q}|$.

In Fig. 3 we show the monomer structure factor for $N = 50, 200$, and $1000$ and different densities. In the dilute regime, we expect [4, 45]

$$S(q) \approx S_1(q) = \begin{cases} N/(1+q^2 R_{g0}^2/3) & q < 2\pi/R_{g0} \\ Aq^{-1/\nu} & 2\pi/R_{g0} < q < 2\pi/b \\ \mathcal{O}(1) & q > 2\pi/b, \end{cases} \quad (25)$$

where $S_1(q)$ is the single chain structure factor, which can be obtained by computing the expression (24) for a single chain, and $A > 0$ is a constant (the $\rho$-dependence of $S_1(q)$ is reported in Sec. II in the S.I.). In Fig. 3a-c, we show as thin continuous lines $S_1(q)$ for $\rho = 10^{-3}$. For $N = 50$ and $200$, we see that at this density $S(q) = S_1(q)$, indicating that the system is already in the infinite dilution limit, whereas for $N = 1000$ this limit is still not yet reached.

In the semidilute regime, we expect [4]

$$S(q) \approx \begin{cases} S(0)/[1 + (q\xi_c)^2] & q < 2\pi/\xi_c \\ B/[1 + (q\xi_c)^{1/\nu}] & 2\pi/\xi_c < q < 2\pi/b \\ \mathcal{O}(1) & q > 2\pi/b, \end{cases} \quad (26)$$

where $B > 0$ is a constant. From our data at low $\rho$, we can observe at intermediate $q$ the slope $-1/\nu \approx -1.70$;

<figure><img src="image_4.png" /><figcaption>FIG. 3. Monomer structure factor $S(q)$, Eq. (24), for $N = 50$ (a), $N = 200$ (b) and $N = 1000$ (c), at different monomer densities. Thin continuous lines: $S_1(q)$ for $\rho = 10^{-3}$. Dashed lines: slope $-1/\nu \approx -1.70$ ($\nu \approx 0.588$). Dash-dotted curve in (a): fit with the Ornstein-Zernike function, Eq. (12), for $\rho = 0.20$.</figcaption></figure>

as expected, this regime is more clearly observable for the longest chains, $N = 1000$ (Fig. 3c).

As $\rho$ is increased, the isothermal compressibility $\kappa_T = S(0)/\rho k_B T$ decreases and eventually, for $\rho \ge 0.85$, the structure becomes virtually indistinguishable from that of a dense liquid [51, 53]. Comparing $S(q)$ for differ-
<figure><img src="image_5.png" /><figcaption>FIG. 4. The function $r|g(r) - 1|$, with $g(r) = \text{radial distribution function, for } N = 200$ and different monomer densities. To facilitate the visualization, every curve starting from $\rho = 0.08$ was shifted by a factor 100 with respect to the previous one, except for $\rho = 0.85$, which is shifted by a factor 200. The dashed lines are exponential fits.</figcaption></figure>

ent chain lengths, we note that for densities $\rho > 0.11 \simeq \rho^*(N = 50)$ (see Eq. (20)), $S(q)$ becomes independent of $N$ (see Fig. S3 in the S.I.). This is in agreement with the prediction that in the semidilute and concentrated regimes, $\rho > \rho^*$, the global structure of the system does not depend on $N$ [4].

The correlation length $\xi_c$ can be measured by fitting $S(q)$ in the low-$q$ range with the Ornstein-Zernike expression, Eq. (12). It is clear, however, that this expression only gives a good description of the data for densities $\rho \lesssim 0.3 = \rho^{**}$. To obtain $\xi_c$ at higher densities we consider the radial distribution function of the monomers, $g(r)$, defined as [51]

$$g(r) = \frac{1}{4\pi M \rho r^2} \sum_{\substack{k=1 \\ j \neq k}}^{M} \langle \delta(|\mathbf{r} + \mathbf{r}_k - \mathbf{r}_j|) \rangle. \quad (27)$$

The expected asymptotic expression of $g(r)$ at low and high $\rho$ is given by Eqs. (9) and (10). At low $\rho$, we can therefore obtain $\xi_c$ from $g(r)$ by fitting it with an exponential, while at higher $\rho$ we fit the exponential envelope of the function [35]. This method is illustrated in Fig. 4, where we show $r|g(r) - 1|$ for $N = 200$ and different values of $\rho$, with the slopes resulting from the fit (dashed lines). One can see that the low-$\rho$ form works well up to density $\rho \simeq 0.11$, while the high-$\rho$ one gives a good description of the data starting from $\rho \simeq 0.64$. For densities in-between, neither of these two expressions gives a really satisfactory description of the data. We decide nevertheless to fit the data with an exponential also in this intermediate range, as shown in Fig. 4, in order to have at least an estimate of $\xi_c$.

The values of $\xi_c$ obtained by fitting $S(q)$ and $g(r)$ are reported in Fig. 5a as a function of $\rho$. We note that the

<figure><img src="image_6.png" /><figcaption>FIG. 5. The length $\xi_c$ measured from $g(r)$ (dark symbols, continuous lines) and from $S(q)$ (light symbols, dashed lines), as a function of monomer density (a) and the scaling variable $\chi$, Eq. (7) (b). Continuous lines represent the slopes predicted by scaling theory for the semidilute regime. Open symbols denote densities $\rho > \rho^{**} = 0.3$</figcaption></figure>

two estimates give almost identical results for $\rho < 0.2$. For $\rho \ge 0.26$, i.e., when we start to approach $\rho^{**} = 0.3$ we cannot fit reliably the structure factor with the Ornstein-Zernike function, Eq. (12), and we have to rely on the radial distribution function. We note that $\xi_c$ reaches a minimum at $\rho \simeq 0.3 = \rho^{**}$ and increases for larger $\rho$. The increase of $\xi_c$ with $\rho$ is due to local packing constraints, and it is a typical behavior for any dense liquid [35].

In Fig. 5b, we plot $\sqrt{3} \xi_c / R_{g0}$ as a function of the scaling variable $\chi$, Eq. (7). At low densities, $\sqrt{3} \xi_c / R_{g0} \simeq 1$, in agreement with the theoretical predictions (Eqs. (12) and (13)). We observe that in the semidilute regime, $\xi_c$ follows approximately the predictions of scaling theory Eq. (15), i.e. $\xi_c \propto \rho^{-\nu/(3\nu-1)}$ (Fig. 5a), or equivalently $\xi_c \propto \chi^{-\nu}$ (Fig. 5b). The agreement with the theory becomes better with increasing $N$, as expected from the fact that $\lim_{N \to \infty} \rho^* = 0$, or analogously $\lim_{N \to \infty} \chi^{**} = \infty$.
<figure><img src="image_7.png" /><figcaption>FIG. 6. Example of a porous medium with a single star-shaped pore. $V_T(r)$ is the volume which is accessible to the center of a spherical probe particle of radius $r$. This volume is enclosed in the surface $S_T$. $V_G(r)$ is the volume which is coverable by spheres of radius $r$ or smaller, and it is enclosed in the Connolly (or reentrant) surface $S_G$.</figcaption></figure>

## B. Pore size distribution

We have seen in the previous subsection how to estimate the monomer correlation length $\xi_c$. As discussed in Sec. II, only in the semidilute regime we have $\xi_c \approx \xi$. In the dilute regime $\xi_c \approx R_{g0}$ while $\xi \propto (\rho/N)^{-1/3}$. In the concentrated regime, $\xi_c$ increases, while $\xi$ is expected to decrease. In order to get the actual geometrical mesh size, we therefore need a way to measure directly the size of the pores in the system. With this objective in mind, we turn to a concept that was developed to measure the distribution of pore sizes in solid porous media: the pore size distribution (PSD) [1, 2, 36–43].

We will consider two different definitions of the PSD: The first one, although introduced by others [36, 37], was formalized and extensively studied by Torquato and coworkers [1, 38–40], while the second one was introduced by Gubbins and coworkers [2, 41–43]. We will in the following consider a generic porous medium in three dimensions, but all the definitions can be extended to any two-phase system and any number of spatial dimensions. Our presentation closely follows the one given by Torquato [1, 38–40]. A porous medium is a three-dimensional domain of volume $V$ which is composed of two sub-domains: a void (or pore) region, of volume $V_p$ and volume fraction $\phi_p \equiv V_p/V$ (also called the porosity), and a solid region of volume $V_s = V - V_p$ and volume fraction $\phi_s \equiv V_s/V = 1 - \phi_p$. In Fig. 6 we show a simple example of a porous medium with a single star-shaped pore.

Following Torquato, we start by defining the quantity $V_T(r)$ as the pore volume accessible to the center of a

spherical particle with radius $r$. From Fig. 6 one recognizes that this is the volume inside the surface $S_T$ followed by the center of the particle when the particle slides over the void-solid interface. Naturally, $V_T(r) \le V_p$. We can then define the fraction of pore volume which is accessible to the particle as $F_T(r) \equiv V_T(r)/V_p$. It is clear from the definition that

$$F_T(0) = 1; \quad \lim_{r \to \infty} F_T(r) = 0. \qquad (28)$$

Torquato’s PSD is then defined as the probability density function [1, 38–40]

$$P_T(r) \equiv -\frac{dF_T(r)}{dr}. \qquad (29)$$

Thus, $P_T(r)dr$ represents the probability that a randomly chosen point in the pore region lies at a distance between $r$ and $r+dr$ from the nearest point on the pore-solid interface. From (28) and (29) it follows that

$$F_T(r) = \int_r^\infty P_T(x)dx, \qquad (30)$$

and that $P_T(r)$ is normalized to unity:

$$\int_0^\infty P_T(x)dx = 1. \qquad (31)$$

Moreover, $P_T(r)$ will vanish for $r \to \infty$.

The definition of Gubbins differs from that of Torquato by a simple, yet significant detail. Instead of considering the part of pore volume which is accessible to the center of a spherical particle of radius $r$, we consider the volume which is accessible to any point of the probe particle. In other words, one defines $V_G(r)$ as the pore volume coverable by spheres of radius $r$ or smaller [2] (Fig. 6). The surface $S_G$ in which $V_G$ is enclosed is sometimes called Connolly surface or reentrant surface [41, 62]. Therefore $P_G(r)dr$ represents the probability that a randomly chosen point in the pore space is coverable by spheres of radius $r$ but not by spheres of radius $r+dr$ [2] [63]. Once $V_G(r)$ is defined, all the other quantities can be defined exactly as done above for Torquato’s PSD. We note that analogous definitions of $V_T$ and $V_G$ have been given by other authors [41, 62, 64, 65].

The PSDs $P_T(r)$ and $P_G(r)$, obtained respectively from $V_T(r)$ and $V_G(r)$, can differ significantly from each other. To see this, let’s consider the simple case of a solid material containing a spherical pore of radius $R$. It is easy to see that in this case [37]:

$$F_T(r) = \left(\frac{R-r}{R}\right)^3; \quad P_T(r) = \frac{3(R-r)^2}{R^3}, \qquad (32)$$

with $r \le R$, while
$$F_G(r) = \begin{cases} 1 & r \le R; \\ 0 & r > R; \end{cases} \qquad P_G(r) = \delta(r - R), \qquad (33)$$

where $\delta$ is Dirac's delta distribution. It is clear that in this case the size $R$ of the pore is most readily identified by considering $P_G$, which is nothing else than a delta distribution centered at $R$, whereas $P_T$ goes to zero at $R$. We will see below that also when dealing with more realistic systems, $P_G$ often conveys the information regarding the typical size of the pore in a much more direct manner than $P_T$.

In order to compute $P_T$ and $P_G$ from the simulation data, we first need to divide the sample in pore and solid regions. For simplicity, we assume that the interaction potential between the probe particle of radius $r$ and the monomers is hard-sphere like, i.e.,

$$U(d) = \begin{cases} \infty & d < r + r_m \\ 0 & d \ge r + r_m, \end{cases} \qquad (34)$$

where $d$ is the distance between the probe particle and the monomer and $r_m$ is the radius of the monomer when it is approximated as a hard sphere. In the present work, we choose $r_m = \sigma/2 = 0.5$ [66]. Once this assumption is made, we can proceed to the calculation of the PSD, following the two definitions. The algorithms used to calculate $P_T$ and $P_G$ are described in Sec. AII.

In Fig. 7, we show the complementary cumulative distribution functions $F_G$ and $F_T$ (Fig. 7a) and the PSDs $P_G$ and $P_T$ (Fig. 7b) for $N = 200$ and different densities. In both graphs, $r$ has been normalized by the length of the simulation box $L = (M/\rho)^{1/3}$ in order to make the plot more readable. We find that $F_T$ drops to zero much faster than $F_G$, and that the shapes of the two functions are quite different. As a consequence, also the distributions $P_T$ and $P_G$ are very different.

In Fig. 8 we show $P_G(\tilde{r})$, where $\tilde{r} \equiv r/\langle r \rangle_G - 1$, for $N = 200$ and different densities. One recognizes that at intermediate and high densities the distribution is very similar to a Gaussian. This becomes clearer when plotting the probability distribution of the variable sgn($\tilde{r}$)·$\tilde{r}^2$ (with sgn the sign function) in semi-logarithmic scale, as shown in the inset of Fig. 8. Only at density $\rho < \rho^* \simeq 0.032$ we observe significant deviations from Gaussianity, with the left-side tail displaying a markedly exponential decay (dashed line in the inset of Fig. 8). The small shoulder observable at low $\rho$ for small values of $\tilde{r}$ (indicated by the arrow in the inset of Fig. 8) comes from distances inside the pervaded volume of the chain.

We note that the width of $P_G(\tilde{r})$ increases with increasing density, indicating that the pore space becomes more heterogeneous on the length scale of the mean pore size $\langle r \rangle_G$. Moreover, although the shape of $P_G(\tilde{r})$ changes significantly when $\rho$ is increased in the dilute regime, it remains almost unchanged for $\rho > \rho^*$. This effect can be

<figure><img src="image_8.png" /><figcaption>FIG. 7. Complementary cumulative distribution $F_\alpha(r)$ ($\alpha = T, G$) (a) and PSD $P_\alpha(r)$ (b), as obtained using the definitions of Gubbins (filled symbols, continuous lines) and of Torquato (open symbols, dashed lines).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 8. Gubbins's PSD, $P_G(\tilde{r})$, for $N = 200$ and at different densities $\rho$, as a function of $\tilde{r} \equiv r/\langle r \rangle_G - 1$. Inset: Probability distribution of the variable sgn($\tilde{r}$)·$\tilde{r}^2$ in semi-logarithmic scale. Dashed line: exponential fit. See main text for the explanation of the arrow.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 9. Relative standard deviation $s_\alpha$ ($\alpha = G, T$), Eq. (35), as a function of the scaling variable $\chi$. $s_G$ is shifted up by 0.2 to facilitate the comparison with $s_T$. Open symbols denote densities $\rho > \rho^{**} = 0.3$</figcaption></figure>

better appreciated by considering the relative standard deviation,

$$s_{\alpha} = \sqrt{\frac{\langle r^2 \rangle_{\alpha}}{\langle r \rangle_{\alpha}^2} - 1} \quad (\alpha = T, G), \quad (35)$$

which is a measure of the width of the PSD relative to its mean value.

In Fig. 9, we show $s_G$ and $s_T$ as a function of the scaling variable $\chi$. One sees that $s_G$ increases with increasing $\chi$ up to the beginning of the semidilute regime at $\chi = \chi^* = 1$, whereas it is a constant for $\chi > \chi^*$. This transition from dilute to semidilute behavior is quite sharp if compared with that of $\xi_c$ (Fig. 5b) or $R_g$ (Fig. 18 of Sec. AI), suggesting that $P_G$ can be used to determine $\rho^*$ in a more precise way. In addition, $s_G$ is almost independent of $N$ for a given value of $\chi$.

The relative standard deviation $s_T$ of Torquato's PSD behaves quite differently, showing a marked dependence on $N$ at fixed $\chi$ for $\chi > 1$. Only for $N = 1000$ a behavior similar to that of $s_G$, with an increase up to $\chi = 1$ followed by a plateau, is recovered. These observations suggest that $P_G$ may be more suitable than $P_T$ to characterize the pore structure of the system. This will also be argued below based on different considerations. Additional details on $P_T$ will be given in Sec. IV D. Finally we note that the value of $s_G$ is significantly smaller than the one of $s_T$ (in Fig. 9 $s_G$ is shifted upwards by 0.2) which shows that $P_G$ is much narrower than $P_T$ when both are normalized by the respective mean values.

### C. Estimating the geometrical mesh size from the pore size distribution

Once the PSD has been determined, we need to extract from it a quantity that we can compare with $\xi_c$ and

<figure><img src="image_11.png" /><figcaption>FIG. 10. Comparison between the average pore sizes, the density fluctuation correlation length $\xi_c$ and the scaling estimate of Eqs. (19)-(20) (with $\rho$ expressed as a function of $\chi$ using Eq. (7)), as a function of the scaling variable $\chi$. Dark filled symbols with continuous lines: $\xi_c/R_{g0}$. Thick-bordered open symbols with dashed lines: $\langle r \rangle_T/R_{g0}$. Light filled symbols with dash-dotted lines: $\langle r \rangle_G/R_{g0}$. Thin-bordered open symbols denote densities $\rho > \rho^{**} = 0.3$. Thick black lines: scaling estimate for $N = 50$ (dash-dotted) $N = 200$ (dashed) and $N = 1000$ (continuous).</figcaption></figure>

with the theoretical expectations for $\xi$. The most natural choice is to consider the mean pore radius,

$$\langle r \rangle_{\alpha} = \int_{0}^{\infty} x P_{\alpha}(x) dx = \int_{0}^{\infty} F_{\alpha}(x) dx \quad (\alpha = T, G). \quad (36)$$

Although other choices are possible, such as the position of the peak of $P(r)$, we found that the main results do not depend significantly on this choice.

In Fig. 10 we show $\langle r \rangle_G/R_{g0}$ and $\langle r \rangle_T/R_{g0}$ as a function of the scaling variable $\chi$ for different values of $N$. These two quantities are compared with $\xi_c/R_{g0}$, where $\xi_c$ is measured from the radial distribution function (points for $\rho > 0.3$ not shown), and with the scaling estimate $\xi/R_{g0} \simeq (\rho/\rho^*)^{-\nu/(3\nu-1)}$, Eq. (19), where $\rho^*$ was estimated using Eq. (20). We observe that, despite the fact that $\langle r \rangle_G > \langle r \rangle_T$, both quantities show a remarkably good agreement with the scaling prediction, Eq. (18). Moreover, above the dilute regime, i.e., $\chi \gtrsim 1$, $\langle r \rangle_T \simeq \xi_c$, likely because the calculation of both quantities involves averaging over distances which are of similar magnitude.

To better see that $\langle r \rangle_\alpha$ ($\alpha = T, G$) follow the scaling prediction, we show in Fig. 11a $\langle r \rangle_\alpha$ multiplied by $\chi^\nu$ (points corresponding to $\rho > 0.3$ not shown): For $\chi \gtrsim 1$, they reach a plateau, the length of which increases with increasing $N$ implying that $\langle r \rangle_\alpha$ is indeed proportional to $\chi^{-\nu}$. Moreover, we also see that for low $\chi$ the function tends towards the predicted behavior $\chi^{-\nu+1/3}$, corresponding to the average distance between neighboring chains, Eq. (18). In Fig. 11b we show the ratio $\langle r \rangle_G/\langle r \rangle_T$ as a function of $\chi$ (points corresponding to $\rho > 0.3$ not shown).
<figure><img src="image_12.png" /><figcaption>FIG. 11. (a) $\chi^\nu \cdot \langle r \rangle_\alpha / R_{g0}$ ($\alpha = G, T$) as a function of the scaling variable $\chi$. Dash-dotted lines: Predicted slopes in the dilute limit, Eq. (18). (b) Ratio $\langle r \rangle_G / \langle r \rangle_T$ as a function of $\chi$. Thick dashed line: $2.18 + 0.0458 \ln(\chi)$. Points corresponding to $\rho > \rho^{**} = 0.3$ are not shown.</figcaption></figure>

shown). For $\chi < 1$, the ratio increases logarithmically with increasing $\chi$ and is independent of $N$. This demonstrates that in the dilute regime $\langle r \rangle_G$ and $\langle r \rangle_T$ are directly related to each other via an $N$-independent function. For $\chi > 1$ the $N = 50$ data follows this logarithmic dependence whereas the larger systems show a plateau the width of which increases with $N$, and only for larger $\chi$ the ratio increases again. In summary one can conclude that this ratio displays a surprisingly weak dependence on $\chi$, changing only by $\simeq 15\%$ over several decades in $\chi$.

To summarize, we have two possible estimators for the mesh size, which are both in excellent agreement with the scaling prediction and differ from each other by a multiplicative factor $\simeq 2$. This factor is of course irrelevant if we are only interested in order-of-magnitude estimates of the mesh size, but it is relevant when more precise information about the size of the pores are needed. In other words, it is important to understand which PSD, $P_G$ or $P_T$, gives us a more precise information about the “real” size of the pores. In the following we will therefore discuss how to interpret the PSD, and which estimator to choose between $\langle r \rangle_G$ and $\langle r \rangle_T$ when more quantitative information about the mesh size is required.

When considering the case of a solid material containing disconnected spherical pores, Eqs. (32) and (33), one can see that it is much more straightforward to infer the value of $R$ from $P_G(r)$ than from $P_T(r)$. Indeed, using Eq. (32) we find $\langle r \rangle_T = R/4$, whereas from Eq. (33) one obtains $\langle r \rangle_G = R$. Thus already this simple example hints that $\langle r \rangle_G$ is a better indicator for the “real” pore size than $\langle r \rangle_T$. In order to consider a somewhat less artificial example, we simulated a simple model of a polymer gel: Polymer strands of length $N = 18$ are placed on the edges of a cubic lattice and connected to each other at the vertices (see snapshot in Fig. 12). Initially, the distance between each pair of bonded monomers is $r = 1$.

<figure><img src="image_13.png" /><figcaption>FIG. 12. PSD for a model polymer gel (cubic lattice with strand length $N = 18$; a snapshot of the system is reported in the main figure). Dashed blue line: Torquato’s PSD. Continuous red line: Gubbins’s PSD. Inset: double logarithmic plot of the same PSDs.</figcaption></figure>

The network is then allowed to relax at constant volume. In order to grant additional flexibility to the chains, we set the $k$ parameter in the FENE potential, Eq. (22), to $k = 20$. By construction, the mesh size of this system is $\xi \simeq N$, and therefore the “ideal” PSD should have a strong signal at $r \simeq N/2 = 9$. More precisely, taking into account the tri-dimensionality of the system and the finite diameter of the monomers, $2r_m = 1$, we expect a strong signal between $(N-1)/2 = 8.5$ (edge of the cubic cell) and $(\sqrt{2}N-1)/2 \simeq 12.2$ (face diagonal of the cubic cell).

In Fig. 12, we show $P_T(r)$ and $P_G(r)$ for this model gel. We observe that, while $P_T(r)$ displays a broad profile which peaks at $r \simeq 4$, $P_G(r)$ has a negligible value over the whole $r$ range, except for a very sharp peak at $r \simeq 11 \simeq (\sqrt{2}N-1)/2$. The information about the mesh size is not, however, completely absent from $P_T(r)$, but it is “hidden” in the value of $r_{max}$ at which $P(r_{max}) = 0$ (the largest pore size), as it becomes clear when plotting $P_T(r)$ in semi-logarithmic scale (inset of Fig. 12). From this semi-logarithmic plot, one can also see that $P_G(r)$ has a secondary peak (much smaller than the main one) at $r \simeq 8 \simeq (N-1)/2$, which originates from the portion of pore space in the vicinity of the surface of the lattice cells. Furthermore, we note that also for this idealized system the ratio $\langle r \rangle_G / \langle r \rangle_T$ is very close to $2 (\langle r \rangle_G / \langle r \rangle_T \simeq 2.6)$, and this remains true also when other densities are considered (not shown). This is a consequence of the fact that for all the systems studied $P_T(r)$ is fairly symmetrical, and $P_G(r)$ peaks close to $r_{max}$. The precise value of $\langle r \rangle_G / \langle r \rangle_T$ will however depend on the relative forms of $P_T$ and $P_G$.

From these two examples, we conclude that $P_G(r)$ is the quantity that gives a more immediate and easy to interpret information on the pores in the system. This does not mean that $P_T(r)$ should not be used, but only that care should be taken in its interpretation.
## D. $P_T(r)$: Comparison with analytical models

In the previous section we have shown how the PSD can be connected to the geometrical mesh size and have compared the PSD of Torquato, $P_T(r)$, to that of Gubbins, $P_G(r)$. In this section, we will illustrate some analytical results which can be used to gain insight into $P_T(r)$. We focus on this distribution because quasi-exact results are available for it [1, 38–40], whereas this is not the case for $P_G(r)$.

Torquato and coworkers have studied extensively the properties of porous media, and in particular the PSD of systems of identical particles interacting via an arbitrary potential. For a system of hard spheres (HS) of radius $R$ at density $\rho$, they demonstrated that [1, 38–40]

$$P_T^{\text{HS}}(x) = \frac{3\eta}{R} F_T^{\text{HS}}(x)(a_0x^2 + 2a_1x + 4a_3), \quad (37)$$

with

$$F_T^{\text{HS}}(x) = \exp[-\eta(a_0x^3 + 3a_1x^2 + 12a_2x + a_3)]. \quad (38)$$

Here $x \equiv (r + R)/R$, $\eta \equiv 4\pi R^3 \rho/3$ is a dimensionless density which for hard spheres is equivalent to the solid volume fraction $\phi_s$, and $a_i$, $i = 0, 1, 2, 3$ are functions of $\eta$ only, whose explicit expression can be found in Ref. 39.

For overlapping spheres (OS), i.e., spheres which can overlap with no energy penalty, they found [38, 40]

$$P_T^{\text{OS}}(x) = \frac{3\eta}{R} x^2 \exp[-\eta(x^3 - 1)]. \quad (39)$$

Note that in this case, $\eta$ is not equivalent to the particle volume fraction $\phi_s$, and indeed it can be shown that $\phi_s = 1 - e^{-\eta}$ [38, 40]. For small values of $\eta$, Eq. (37) reduces to Eq. (39), i.e., the HS model is equivalent to the OS model [38, 40].

We make the assumption that the measured PSD $P_T(r)$ can be fitted by one of the two functional forms (37)-(39), with $R$ and $\rho$ as fit parameters. In other words, we assume that it is possible to map our system on a system of hard (resp. overlapping) spheres with radius $R_{\text{HS}}(\rho, N)$ (resp. $R_{\text{OS}}(\rho, N)$) and density $\rho_{\text{HS}}(\rho, N)$ (resp. $\rho_{\text{OS}}(\rho, N)$).

In Fig. 13 we compare for $N = 200$ and different monomer densities the calculated PSD, with the result of the fits using the HS model, $P_T^{\text{HS}}(r; R_{\text{HS}}, \rho_{\text{HS}})$, and with the OS model $P_T^{\text{OS}}(r; R_{\text{OS}}, \rho_{\text{OS}})$. It is virtually impossible to distinguish $P_T^{\text{HS}}(r; R_{\text{HS}}, \rho_{\text{HS}})$ from $P_T^{\text{OS}}(r; R_{\text{OS}}, \rho_{\text{OS}})$ via a simple visual inspection, since the two curves superimpose almost perfectly. An analysis of the squared difference between the calculated and the fitted function (not shown) reveals however that the HS model fits the data slightly better for all densities, except $\rho = 0.85, 1.00$.

In Figs. 14a and b we show, respectively, $R_{\text{HS}}, R_{\text{OS}}$ and $\rho_{\text{HS}}/\rho, \rho_{\text{OS}}/\rho$ as function of the monomer density.

<figure><img src="image_14.png" /><figcaption>FIG. 13. Torquato’s PDS, $P_T(r)$, for $N = 200$ and at different densities (open symbols), compared with the fit results $P_T^{\alpha}(r)$ (with $\alpha = \text{OS}$ or HS) employing the HS (Eq. (37), points) and OS (Eq. (39), lines) models. The HS and OS fits superimpose almost perfectly.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 14. Fit parameters as a function of monomer density $\rho$ for the fit of the PSD $P_T(r)$ with the HS (Eq. (37), filled symbols) and OS (Eq. (39), open symbols) models. (a): Effective radius. (b): Effective density divided by $\rho$.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIG. 15. Dimensionless densities $\eta_{OS}$ (open symbols) and $\eta_{HS}$ (filled symbols), Eq. (41), calculated using the fit parameters from the OS and HS model fits, as a function of monomer density. Continuous lines: $\eta_c(0, N) = 4\pi R_{g0}^3 \rho/3N$ (see Eq. (42)). The approximate values of the overlap densities $\rho^*(N)$ have been estimated using Eq. (20).</figcaption></figure>

Through an extrapolation of the low density behavior, we deduce that in the limit $\rho \to 0$, $R_{HS}$, $R_{OS} \to R_{g0}(N)$, and $\rho_{HS}$, $\rho_{OS} \to \rho/N$ (dashed horizontal lines in Figs 14a and 14b). In other words, $P_T(r)$ is described well by the PSD of a system of spheres with radius $R_{g0}$, and density equal to the chain density $\rho/N$. This is valid independently of the model used (OS or HS), since, as mentioned above, both models give the same result in the limit $\eta \to 0$, or equivalently $\rho \to 0$, $R = \text{const}$.

In the semidilute regime $\rho^*(N) < \rho < \rho^{**}$, both $R_\alpha$ and $\rho_\alpha$ ($\alpha = \text{OS,HS}$) tend towards a power-law behavior with increasing $N$. Since in the semidilute regime the only relevant length scale is the mesh size $\xi$, we make the hypothesis that for large $N$ these quantities will show the following scaling behavior:

$$ 
\begin{aligned}
R_{\alpha} &\propto \xi \approx \rho^{-\nu/(3\nu-1)} \\
\rho_{\alpha} &\propto \xi^{-3} \approx \rho^{3\nu/(3\nu-1)}.
\end{aligned} 
\quad (\alpha = \text{OS,HS}) \quad (40) 
$$

This *ansatz* is motivated by the fact that the PSD is a purely geometrical object and hence has to scale with the intrinsic length scales of the system. In Fig. 15 we show the dimensionless density

$$ \eta_{\alpha} \equiv 4\pi R_{\alpha}^3 \rho_{\alpha} / 3 \quad (\alpha = \text{OS,HS}) \quad (41) $$

as a function of $\rho$. We note that both $\eta_{HS}$ and $\eta_{OS}$ depend basically only on $(\rho, N)$, with $\eta_{OS}$ being slightly smaller than $\eta_{HS}$. At low $\rho$, both $\eta_{OS}$ and $\eta_{HS}$ approach the asymptotic expression $\eta_c(0, N) = 4\pi R_{g0}^3 \rho/3N$ (continuous lines), where $\eta_c$ is in general defined as

$$ \eta_c(\rho, N) \equiv 4\pi R_g^3 \rho/3N, \quad (42) $$

<figure><img src="image_17.png" /><figcaption>FIG. 16. Torquato's PSD for $N = 200$ and at different densities (open symbols), compared with the OC model prediction, Eqs. (39)-(43).</figcaption></figure>

i.e., the dimensionless density the chains would have were they spheres of radius $R_g$. Importantly, we observe that when $N$ increases, $\eta_{OS}$ and $\eta_{HS}$ develop a plateau in the semidilute regime: This corroborates the hypothesis that for large $N$ Eq. (40) is valid, since this equation implies $\eta_\alpha \approx \text{const}$.

We have seen in the preceding paragraph that at low $\rho$ the polymer solution can be mapped on a system of overlapping spheres of radius $R_g$ and density equal to the chain density, $\rho/N$. Inspired by this result, we introduce the overlapping chains (OC) model, by assuming that the PSD of the system is approximated for all $\rho$ by Eq. (39), with

$$ 
\left\{ 
\begin{aligned}
R_{\text{OS}} &= R_g(\rho, N) \\
\rho_{\text{OS}} &= \rho/N \\
\eta_{\text{OS}} &= \eta_c.
\end{aligned} 
\right. 
\quad (43) 
$$

This approximation fails at intermediate and high $\rho$, as it is clear from Fig. 14a-b. This is also illustrated in Fig. 16, where the calculated PSD for $N = 200$ and different densities is compared with the one predicted by the OC model and one sees that already for $\rho$ as low as 0.02 the two PSDs are qualitatively very different. However, given the simplicity of the OC model, we decide nevertheless to compare the computed average pore size, $\langle r \rangle_T$, with the one that the model predicts.

From the PSD (39) one can compute the average pore size with the result

$$ \langle r \rangle_{T}^{\text{OS}} = \frac{R}{3} \exp(\eta) E_{2/3}(\eta), \quad (44) $$

where $E_n(\eta)$ is an exponential integral [67],

$$ E_n(\eta) = \int_1^\infty \frac{e^{-\eta x}}{x^n} dx. \quad (45) $$
<figure><img src="image_18.png" /><figcaption>FIG. 17. Reduced average pore size $\langle r \rangle_T/R_g$ (symbols) compared with the OC model prediction $\langle r \rangle_T^{OC}/R_g$ (thick line), Eq. (43), as a function of $\eta_c$. Open symbols correspond to densities $\rho > 0.3$ (concentrated regime). Dashed (dash-dotted) line: Low (high) density behavior predicted by the OC model, Eq. (49).</figcaption></figure>

The OC model result is obtained simply by applying Eqs. (43) to Eq. (44), i.e.,

$$ \langle r \rangle_T^{OC} = \frac{R_g(\rho, N)}{3} \exp(\eta_c) E_{2/3}(\eta_c). \quad (46) $$

In Fig. 17, we compare the measured average pore size $\langle r \rangle_T$ with the OC model result $\langle r \rangle_T^{OC}$. One observes that, except for the points corresponding to the concentrated regime ($\rho > 0.3$), the prediction of the OC model is in surprisingly good agreement with the data, without using any fit parameter. This is a striking result, since we have seen that the calculated PSD starts to diverge from the theoretical prediction already at relatively low density (Fig. 16).

Another related surprising result is the following. The scaling behavior of $R_g$ is [3, 4, 45]

$$ R_g \approx \begin{cases} R_{g0} & \chi < \chi^* \\ R_{g0} \chi^{-\nu+1/2} & \chi^* < \chi < \chi^{**}. \end{cases} \quad (47) $$

This implies, given the definition of $\eta_c$, that

$$ \eta_c \approx \begin{cases} \chi^{3\nu-1} & \chi < \chi^* \\ \chi^{1/2} & \chi^* < \chi < \chi^{**}. \end{cases} \quad (48) $$

The function appearing on the right hand side of Eq. (46) can be expanded for $\eta_c \to 0$ and $\eta_c \to \infty$ via a Puiseux series [68]; the leading terms are

$$ \exp(\eta_c) E_{2/3}(\eta_c) = \begin{cases} \Gamma(\frac{1}{3}) \eta_c^{-1/3} + O(1) & \eta_c \to 0 \\ \eta_c^{-1} + O(\eta_c^{-2}) & \eta_c \to \infty, \end{cases} \quad (49) $$

where $\Gamma$ is the Gamma function, and therefore

$$ \frac{\langle r \rangle_T^{OC}}{R_g} \approx \begin{cases} \eta_c^{-1/3} & \eta_c \to 0 \\ \eta_c^{-1} & \eta_c \to \infty. \end{cases} \quad (50) $$

From Eqs. (47), (48), (50), along with the fact that $\lim_{N \to \infty} \chi^{**}(N) = \infty$, we obtain for infinitely long chains

$$ \frac{\langle r \rangle_T^{OC}}{R_{g0}} = \frac{\langle r \rangle_T^{OC}}{R_g} \frac{R_g}{R_{g0}} \approx \begin{cases} \chi^{-\nu+1/3} & \chi < \chi^* \\ \chi^{-\nu} & \chi^* < \chi. \end{cases} \quad (51) $$

We thus find that for very long chains, $\langle r \rangle_T^{OC}/R_{g0}$ behaves exactly as $\xi/R_{g0}$, Eq. (18), for all values of $\chi$. We stress, however, that while $\lim_{N \to \infty} \chi^{**}(N) = \infty$, $\rho^{**}$ is independent of $N$: Therefore, when considering functions of density instead of functions of $\chi$, we must be aware that the approximation $\langle r \rangle_T \simeq \langle r \rangle_T^{OC}$ will always break down for $\rho \ge \rho^{**} = 0.3$.

It is not evident whether there is a deep reason behind the fact that $\langle r \rangle_T$ and $\langle r \rangle_T^{OC}$ are so similar even if $P_T(r)$ and $P_T^{OC}(r)$ are qualitatively very different. Nonetheless, this remarkable fact allows to estimate $\langle r \rangle_T$ for $\rho < \rho^{**}$ with good accuracy by knowing only $R_g$ and $\rho$. Since we have seen that $\langle r \rangle_G \simeq 2\langle r \rangle_T$ (Sec. IV C), this means that we are also able to estimate with reasonable accuracy the average value of the PSD of Gubbins. This is a useful result for estimating $\xi$ in real systems, since measuring the PSD experimentally is not an easy task [69].

## V. SUMMARY AND CONCLUSIONS

The geometrical mesh size $\xi$ is a key quantity in polymer solutions and networks, since it governs the diffusion of particles or molecules in these systems. Conventional methods to measure the mesh size rely on scaling estimates or on measurements of the monomer density fluctuation correlation length $\xi_c$. Scaling estimates, however, only give $\xi$ up to an unknown multiplicative factor, and identifying $\xi_c$ with $\xi$ works only in the semidilute regime, whereas it gives nonsensical results in dense systems.

We propose a method to directly probe the geometrical mesh size $\xi$, which we identify with the average size of the pores of the polymer solution, i.e., the space filled by the solvent. To obtain the distribution of the pore sizes, we make use of the definitions of pore size distribution (PSD) formalized by Torquato [40] and by Gubbins [2]. We find that in both cases, the average of the PSD, $\langle r \rangle$, follows the expected scaling behavior for $\xi$ with high accuracy, validating the hypothesis that $\langle r \rangle$ can be identified with the geometrical mesh size. An important feature of our approach is that it naturally gives access to the distribution of mesh (pore) sizes in a polymer medium: This aspect can be quite relevant, since we know from experiments on polymer networks that the degree of heterogeneity of the medium can strongly influence the diffusion of particles [33].
The range of applicability of the PSD method to evaluate the mesh size is very broad, and by no means limited to polymer solutions. Moreover, it can also be extended to polymer networks and gels. We point out, however, that in gels there exist another relevant “mesh size”, which is related to the average distance between neighboring crosslinks and is in general quite different from $\xi$ [70]. If, for example, a polymer melt is crosslinked, one expects the structural properties, and therefore $\xi$, to remain basically unchanged with respect to the un-crosslinked system [71]. However, the diffusivity of a spherical probe particle will change significantly if the particle diameter is comparable to the average distance between the crosslinks. The same is of course true for entangled polymer melts, where $\xi \simeq 0$ and the relevant parameter controlling the diffusion of particles and the mechanical properties is nothing else than the diameter of the reptation tube [44], since on timescales shorter than the entanglement relaxation time the entanglements effectively act as crosslinks [30, 32].

Another factor which should be taken into account is that the PSD does not depend on the flexibility of the polymers, although in reality this can have a very strong influence on the diffusion of nanoparticles in polymeric systems. For example, it has recently been shown in simulations that increasing the chain rigidity will decrease the diffusivity of nanoparticles in polymer solutions [28]. This rigidity will, however, also affect estimates based on scaling predictions and on measurements of the correlation length, which in addition only give access to an average mesh size, and not to the whole distribution.

To summarize, the method which we have proposed to measure the mesh size should represent a substantial improvement over conventional methods such as scaling estimates and methods which identify the mesh size with the correlation length. It will be interesting, in further studies, to address the problem of characterizing the relevant mesh size in networks and gels and understanding the role played by polymer flexibility.

# APPENDIX

## AI. RADIUS OF GYRATION AND BOND-BOND CORRELATION FUNCTION

As a preliminary analysis of the structural properties of the simulated systems, we have studied the behavior of the radius of gyration $R_g$ as a function of monomer density $\rho$. The radius of gyration is defined as

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

where $\mathbf{r}_{\text{cm}}$ is the position of the center of mass of the chain and $\langle \cdot \rangle$ represents the thermodynamic average. Scaling theory predicts the dependence of $R_g$ on $\rho$ to be [3, 4, 45]

<figure><img src="image_19.png" /><figcaption>FIG. 18. Normalized polymer radius of gyration, Eq. (52), as a function of monomer density $\rho$. The dashed lines have slope $-\frac{\nu-1/2}{3\nu-1} \simeq -0.115$ ($\nu \simeq 0.588$). The overlap densities have been estimated using Eq. (20).</figcaption></figure>

$$\frac{R_g}{R_{g0}} \approx \begin{cases} 1 & \rho < \rho^* \\ (\rho/\rho^*)^{-\frac{\nu-1/2}{3\nu-1}} & \rho^* < \rho < \rho^{**} \\ N^{1/2-\nu} & \rho > \rho^{**} \end{cases} \quad (53)$$

where $R_{g0} \equiv \lim_{\rho \to 0} R_g(\rho)$ and $\rho^*$ is the overlap density.

In Fig. 18, we report $R_g/R_{g0}$ as a function of $\rho$, with the dashed lines representing power laws with exponent $-\frac{\nu-1/2}{3\nu-1} \simeq -0.115$. The radius of gyration in the dilute limit, $R_{g0}$, was computed by simulating a single polymer chain at very low density ($\rho = 10^{-3}$ for $N = 50, 200$ and $\rho = 10^{-4}$ for $N = 1000$): The results are reported in Tab. II, alongside with the resulting estimate for the overlap density, given by Eq. (20). Figure 18 shows that $R_g/R_{g0}$ follows the scaling predictions for all values of $N$ up to the density $\rho \simeq 0.3$. For $\rho \gtrsim 0.3$, the decrease is steeper than what is predicted by scaling theory. This is due to the fact that in this density range the persistence length is density-dependent, an effect which scaling theory does not take into account.

The $\rho$-dependence of the persistence length can be understood by studying the bond-bond correlation function $\langle \cos(\theta_s) \rangle$, where $\cos(\theta_s)$ is defined as [72]

$$\cos(\theta_s) = \frac{\mathbf{b}_n \cdot \mathbf{b}_{n+s}}{|\mathbf{b}_n| |\mathbf{b}_{n+s}|} \quad (s = 1, \dots, N-n), \quad (54)$$

where $\mathbf{b}_n \equiv \mathbf{r}_{n+1} - \mathbf{r}_n$ is the $n$-th bond vector, with $\mathbf{r}_n$ the monomer's position vector. For very long chains, $\langle \cos(\theta_s) \rangle$ is expected to have the following behavior [72–74]:

$$\langle \cos(\theta_s) \rangle \propto \begin{cases} s^{-2(1-\nu)} & l_p/l_b < s \ll N_b(\rho) \\ s^{-3/2} & N_b(\rho) \ll s \ll N, \end{cases} \quad (55)$$
<figure><img src="image_20.png" /><figcaption>FIG. 19. Bond-bond correlation function $\langle\cos(\theta_s)\rangle$, Eq. (54), for $N = 200$ and 1000 and for different values of $\rho$. The dashed line indicates the predicted theoretical behavior for small $s$, whereas the continuous lines indicate the predicted high-$s$ behavior (see Eq. (55)). We note that $2(1-\nu) \simeq 0.824$.</figcaption></figure>

<table><caption>TABLE II. N-dependent properties of the simulated systems: Radius of gyration in the dilute limit, $R_{g0}$ and approximate overlap concentration (Eq. (20)).</caption><thead><tr><th>N</th><th>R<sub>g0</sub></th><th>3N/4πR<sub>g0</sub><sup>3</sup></th></tr></thead><tbody><tr><td>50</td><td>4.83 ± 0.02</td><td>1.1 · 10<sup>-1</sup></td></tr><tr><td>200</td><td>11.44 ± 0.07</td><td>3.2 · 10<sup>-2</sup></td></tr><tr><td>1000</td><td>29.9 ± 0.2</td><td>8.9 · 10<sup>-3</sup></td></tr></tbody></table>

where $l_p$ is the persistence length [4], $l_b = \sqrt{\langle b_n \rangle}$ is the RMS bond length and $N_b$ is the number of monomers in a blob (see Eq. (1)). For the fully flexible Kremer-Grest model considered here, $l_p \approx l_b$ [75]. Since in the dilute regime $N_b \approx N$, we expect in this regime $\langle\cos(\theta_s)\rangle \propto s^{-2(1-\nu)} \forall s > l_p/l_b \approx 1$. Analogously, since in the concentrated regime $N_b \approx 1$, we expect $\langle\cos(\theta_s)\rangle \propto s^{-3/2} \forall s \ll N$.

In Fig. 19 we report $\langle\cos(\theta_s)\rangle$ for $N = 200$ and 1000 and different densities (for $N = 50$ we find similar results). The chains considered here are too short to clearly observe the small-$s$ regime $s^{-2(1-\nu)}$ (dashed line). On the other hand, the large-$s$ behavior is compatible with the theoretical prediction $\langle\cos(\theta_s)\rangle \propto s^{-3/2}$, as illustrated by the continuous lines. Although it is difficult to precisely quantify the persistence length outside of the melt regime for highly flexible and rather short chains as those considered here [74], the marked reduction of $\langle\cos(\theta_s)\rangle$ as a function of $\rho$ observed for small value of $s$ is compatible with a reduction of the persistence length $l_p$. A reduction of $l_p$ with increasing density can also be inferred by considering the distribution of the bond angle $\theta = \pi - \theta_1$ (Sec. S4 in the S.I.) These observations explain qualitatively the high-$\rho$ behavior of $R_g$ observed in Fig. 18.

## AII. CALCULATION OF THE PORE SIZE DISTRIBUTION

### A. Torquato's PSD

The algorithm to compute Torquato's PSD is described in Ref. 40. First of all, as mentioned in Sec. IV B, we have to divide the sample in a “pore” region and a “solid” region. In order to do so, we approximate the monomers as hard spheres of diameter $\sigma$. The procedure to calculate $P_T(r)$ is then as follows:

1. A random point $\mathbf{r}$ is chosen in the pore phase.

2. The smallest distance $r_{\min}$ between $\mathbf{r}$ and the center of a monomer is calculated. The distance between $\mathbf{r}$ and the nearest pore-solid interface is calculated as $r = r_{\min} - \sigma/2$.

This procedure is repeated many times, until a large number $N$ of $r$ values is recorded. In the limit of large $N$, the normalized histogram of these values is equivalent to $P_T(r)$.

### B. Gubbins's PSD

In order to calculate $P_G(r)$, we have used the algorithm proposed by Bhattacharya and Gubbins in Ref. 43. This algorithm is based on the observation that the problem of finding the largest sphere containing $\mathbf{r}$ and which does not overlap with any monomer can be reformulated as the problem of maximizing the function

$$r(\mathbf{r}_c) \equiv \min_{i=1,...,M} \{s_i\} - \sigma/2, \quad (56)$$

subject to the constraint

$$|\mathbf{r}_c - \mathbf{r}| - r(\mathbf{r}_c) \le 0, \quad (57)$$

where $\mathbf{r}_c$ is the position of the sphere's center and $s_i$ is the distance between $\mathbf{r}_c$ and the centers $\mathbf{r}_i$ of the monomers: $s_i = |\mathbf{r}_c - \mathbf{r}_i|$. If the maximization of the function (56) is carried out for a large enough number of points $\mathbf{r}$, the resulting (normalized) histogram of $r$ values will converge to $P_G(r)$.

The problem of calculating $P_G(r)$ reduces therefore to a nonlinear optimization problem, which can be solved with a standard algorithm [76]. We have used the open-source Sbplx algorithm of the NLopt library [77], which is a re-implementation of Subplex [78]. Other choices are possible, although one has to make sure that the chosen algorithm can handle discontinuous objective functions [43].
ACKNOWLEDGMENTS

We thank J. Baschnagel and K. Schweizer for useful discussions, and L. Rovigatti for implementing the code which was used to calculate $P_G$. This work has been supported by LabEx NUMEV (ANR-10-LABX-20) funded by the Investissements d'Avenir French Government program, managed by the French National Research Agency (ANR).

SUPPLEMENTARY INFORMATION

SI. MEAN-SQUARED DISPLACEMENT

In order to check that the time we used to equilibrate the system is sufficiently long, we consider two quantities: the mean-squared displacement (MSD) of the monomers, which we compute as [75]

$$g_1(t) = \frac{1}{(N/2 + 1)} \sum_{i=N/4}^{3N/4} \langle [\mathbf{r}_i(t) - \mathbf{r}_i(0)]^2 \rangle, \quad (S1)$$

and the MSD of the centers of mass of the chains,

$$g_3(t) = \langle [\mathbf{r}_{cm}(t) - \mathbf{r}_{cm}(0)]^2 \rangle, \quad (S2)$$

where $\mathbf{r}_{cm}(t)$ is the position vector of the center of mass of the chain and $\langle \cdot \rangle$ denotes, as usual, the thermodynamic average. Note that in Eq. (S1), only the central half of the monomers belonging to each chain are considered, in order to suppress the fluctuations caused by chain ends [75].

In the present work, the motion of every monomer is governed by the Langevin equation, and hydrodynamic interactions between monomers are neglected. Therefore, in the absence of entanglements the Rouse model [4, 44] gives a good approximation of chain dynamics at all densities, whereas in the presence of entanglements, the dynamics is described by the reptation model [4, 44]. In both cases, when $t > t_{rel}$, with $t_{rel}$ the longest relaxation time of the system, we expect $g_1(t) = g_3(t) \propto t$ [44, 75]: Therefore, we can check that the system has equilibrated by verifying that this condition is met at long times.

In Fig. S1, we show $g_1(t)$ and $g_3(t)$ for different values of $N$ and $\rho$. For the fully-flexible Kremer-Grest model employed here, the entanglement length $N_e$ is $\approx 85$ at $\rho = 0.85$ ($T = 1.0$) [79], and it decreases with decreasing $\rho$ [80–82]. Therefore, for $N = 50$ (Fig. S1a) the system is unentangled, for $N = 200$ (Fig. S1b) it is lightly entangled, and for $N = 1000$ (Fig. S1c) it is entangled. This is confirmed by the fact that for $N = 50$ the behavior of $g_1$ and $g_3$ is in qualitative agreement with the predictions of the Rouse model, whereas for $N = 200$ and 1000 it follows approximately the reptation model [4, 44].

<figure><img src="image_21.png" /><figcaption>FIG. S1. MSD of the monomers, $g_1$ (Eq. (S1)), and of the centers of mass of the chains, $g_3$ (Eq. (S2)), for chain lengths $N = 50$ (a), $N = 200$ (b) and $N = 1000$ (c) and different monomer densities $\rho$. Slopes are reported for comparison with the Rouse model and the reptation model [4, 44].</figcaption></figure>

This is also evidenced by the comparison with the slopes reported in Fig. S1.

For all values of $N$ and $\rho$ considered, we observe that for large times, $g_1(t) \approx g_3(t)$, and thus we conclude that the equilibration time $t_e$ (see Tab. I in the main text) is larger than the longest relaxation time of the system, $t_{rel}$, and therefore all the systems are well equilibrated.
<figure><img src="image_22.png" /><figcaption>FIG. S2. Single chain form factor for chain lengths N = 50 (a), N = 200 (b) and N = 1000 (c) and different monomer densities ρ. Continuous lines: slope −2 (ideal chain). Dashed lines: slope −1/ν ≃ −1.70 (swollen chain).</figcaption></figure>

## SII. CHAIN FORM FACTOR

In Fig. S2, we show the chain form factor $S_1(q)/N$ for N = 50, 200, and 1000 and for different values of the monomer density ρ.

We recall that the theoretical prediction for $S_1(q)/N$ is [4, 45, 48]

$$ \frac{S_1(q)}{N} = \begin{cases} 1/(1+q^2 R_g^2/3) & q < 2\pi/R_g \\ Aq^{-2} & 2\pi/R_g < q < 2\pi/\xi_c \\ Bq^{-1/\nu} & 2\pi/\xi_c < q < 2\pi/b \\ O(1/N) & q > 2\pi/b, \end{cases} \quad (S3) $$

where A, B > 0 are constants. In the dilute regime, $\xi_c \approx R_{g0}$ and therefore the regime $S_1(q)/N \propto q^{-2}$ disappears; in the concentrated/melt regime, $\xi_c \approx b$ and therefore the regime $S_1(q)/N \propto q^{-1/\nu}$ is not present (here $\nu = 0.588$ is the Flory exponent in good solvent).

Swollen chain behavior $S_1(q)/N \propto q^{-1/\nu}$ and the ideal chain behavior $S_1(q)/N \propto q^{-2}$ are clearly observable in Fig. S2 respectively at low and high densities for N = 50 and N = 200. At intermediate densities, we expect the form factor to transition from a $q^{-2}$ dependence to a $q^{-1/\nu}$ dependence at $q \approx 1/\xi_c$ [3, 4]. Since $S_1(q)/N$ assumes values from $\approx 1/N$ to 1, this transition is more easily observable for the longer chains: N = 200 and N = 1000.

## SIII. COMPARISON OF $S(q)$ FOR DIFFERENT VALUES OF N

In Fig. S3, we compare the monomer structure factor $S(q)$ for systems with the same density ρ but with different chain lengths N. One can see that for densities $\rho \gtrsim 0.11$, which is approximately the overlap density of the N = 50 system, $\rho^*(N = 50) \simeq 0.106$ (see Tab. II in the Appendix of the main text), $S(q)$ is independent of N. This is expected, since for $\rho > \rho^*(N = 50)$ all the systems here considered are in the semidilute regime, where the global structure is independent of chain length [4].

<figure><img src="image_23.png" /><figcaption>FIG. S3. Comparison between the monomer structure factors $S(q)$, Eq. (24), of systems with different chain lengths N at the same monomer density ρ. For densities $\rho > 0.11 \simeq \rho^*(N = 50)$, $S(q)$ becomes independent of N.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>FIG. S4. Bond angle distribution for N = 50 and different monomer densities ρ . The bond angle θ is defined by Eq. (S4). Cartoon: schematic representation of how θ is defined.</figcaption></figure>

SIV. BOND ANGLE DISTRIBUTION

In Fig. S4 we report the bond angle distribution P ( θ ) for N = 50 and different densities (for larger N , the results are basically the same). The angle θ is defined as

$$
\theta \equiv \arccos \left( - \frac{\mathbf{r}_{i-1,i} \cdot \mathbf{r}_{i,i+1}}{|\mathbf{r}_{i-1,i}| |\mathbf{r}_{i,i+1}|} \right), \quad (S4)
$$

where r$_{i,j}$ ≡ r$_{i}$ − r$_{j}$ and r$_{i}$ , i = 2 , . . . , N − 1 is the monomer’s position vector. We have therefore θ = π − θ$_{1}$ , where θ$_{s}$ is defined in Eq. (54) of the main text.

At low density, P ( θ ) shows a maximum at θ ≃ 112 ◦ and falls to zero rather sharply at θ ≃ 60 ◦ because of the excluded volume interaction. When ρ is increased past ρ ≃ 0 . 3, the shape of P ( θ ) starts to change significantly in that it develops a peak at θ ≃ 70 ◦ , signaling that the chains are compressed. Overall, the average bond angle decreases by ≃ 8% when going from the dilute regime to density ρ ≳ 0 . 85. This effect is also observable in the average bond length (not shown), although the decrease is in this case only ≃ 1%. Although it is possible to con- ceive chain conformations with a small 〈 θ 〉 but a high persistence length (e.g., an “accordion-like” rigid poly- mer), the fact that up to very high ρ the form factor is compatible with the scaling predictions for chains in a concentrated solution (Sec. SII) suggests that the reduc- tion of 〈 θ 〉 corresponds to a reduction of the persistence length, as also observed from the bond-bond correlation function (Sec. AI in the main text).

[1] S. Torquato, B. Lu, and J. Rubinstein, Phys. Rev. A 41 , 2059 (1990).

[2] L. D. Gelb and K. Gubbins, Langmuir 15 , 305 (1999).

[3] P.-G. De Gennes, Scaling concepts in polymer physics (Cornell University Press, 1979).

[4] M. Rubinstein and R. H. Colby, Polymer physics (Oxford University Press, 2003).

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

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

[7] K. S. Soppimath, T. M. Aminabhavi, A. R. Kulkarni, and W. E. Rudzinski, J. Controlled 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] B. S. Schuster, J. S. Suk, G. F. Woodworth, and J. Hanes, Biomaterials 34 , 3439 (2013).

[10] A. G. Ogston, B. Preston, and J. Wells, Proc. R. Soc. Lond. A 333 , 297 (1973).

[11] A. R. Altenberger and M. Tirrell, J. Chem. Phys. 80 , 2208 (1984).

[12] L. Johansson, C. Elvingson, and J. E. Loefroth, Macro- molecules 24 , 6024 (1991).

[13] G. D. Phillies, G. S. Ullmann, K. Ullmann, and T.-H. Lin, J. Chem. Phys. 82 , 5242 (1985).

[14] G. D. Phillies, Macromolecules 19 , 2367 (1986).

[15] R. Cukier, Macromolecules 17 , 252 (1984).

[16] T. Odijk, Biophys. J. 79 , 2314 (2000).

[17] T.-H. Fan, J. K. Dhont, and R. Tuinier, Phys. Rev. E 75 , 011803 (2007).

[18] R. Tuinier and T.-H. Fan, Soft Matter 4 , 254 (2008).

[19] S. Egorov, J. Chem. Phys. 134 , 084903 (2011).

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

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

[22] L. Masaro and X. Zhu, Prog. Polym. Sci. 24 , 731 (1999).

[23] L.-H. Cai, S. Panyukov, and M. Rubinstein, Macro- molecules 44 , 7853 (2011).

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

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

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

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

[28] R. Chen, R. Poling-Skutvik, M. P. Howard, A. Nikoubashman, S. A. Egorov, J. C. Conrad, and J. C. Palmer, Soft Matter (2019).

[29] S. R. Lustig and N. A. Peppas, J. Appl. Polym. Sci. 36 , 735 (1988).

[30] Z. E. Dell and K. S. Schweizer, Macromolecules 47 , 405 (2013).
[31] S. Rose, A. Marcellan, D. Hourdet, C. Creton, and T. Narita, *Macromolecules* **46**, 4567 (2013).

[32] L.-H. Cai, S. Panyukov, and M. Rubinstein, *Macromolecules* **48**, 847 (2015).

[33] E. Parrish, M. A. Caporizzo, and R. J. Composto, *J. Chem. Phys.* **146**, 203318 (2017).

[34] E. Parrish, S. C. Seeger, and R. J. Composto, *Macromolecules* **51**, 3597 (2018).

[35] R. Koshy, T. Desai, P. Keblinski, J. Hooper, and K. Schweizer, *J. Chem. Phys.* **119**, 7599 (2003).

[36] A. Scheidegger, *The physics of flow through porous media* (Macmillan, 1957).

[37] S. Prager, *Chem. Eng. Sci.* **18**, 227 (1963).

[38] S. Torquato and M. Avellaneda, *J. Chem. Phys.* **95**, 6477 (1991).

[39] S. Torquato, *Phys. Rev. E* **51**, 3170 (1995).

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

[41] K. T. Thomson and K. E. Gubbins, *Langmuir* **16**, 5761 (2000).

[42] J. Pikunic, C. Clinard, N. Cohaut, K. E. Gubbins, J.-M. Guet, R. J.-M. Pellenq, I. Rannou, and J.-N. Rouzaud, *Langmuir* **19**, 8565 (2003).

[43] S. Bhattacharya and K. E. Gubbins, *Langmuir* **22**, 7726 (2006).

[44] M. Doi and S. F. Edwards, *The theory of polymer dynamics* (Oxford University Press, 1988).

[45] I. Teraoka, *Polymer solutions: an introduction to physical properties* (John Wiley & Sons, 2002).

[46] N. Clisby, *Phys. Rev. Lett.* **104**, 055702 (2010).

[47] We are assuming here that $b \approx \sigma$, with $\sigma$ the monomer size, and that the solvent is athermal, so that the size of the thermal blob is approximately equal to $b$.

[48] W. Paul, K. Binder, D. W. Heermann, and K. Kremer, *J. Phys. II* **1**, 37 (1991).

[49] I. Gerroff, A. Milchev, K. Binder, and W. Paul, *J. Chem. Phys.* **98**, 6526 (1993).

[50] We choose to use $\chi$ and not $\rho/\rho^*$ as a scaling variable since $\rho^*$ is not a sharply defined quantity.

[51] J.-P. Hansen and I. R. McDonald, *Theory of simple liquids* (Elsevier, 1990).

[52] M. E. Fisher, *J. Math. Phys.* **5**, 944 (1964).

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

[54] R. Evans, R. Leote de Carvalho, J. Henderson, and D. Hoyle, *J. Chem. Phys.* **100**, 591 (1994).

[55] R. Evans, J. Henderson, D. Hoyle, A. Parry, and Z. Sabeur, *Mol. Phys.* **80**, 755 (1993).

[56] To be more precise, Eq. (12) follows formally from Eq. (11) if $g(r) = 1 + Ae^{-r/\xi_c}/r - \delta(r)/\rho$. The presence of the delta distribution is however only an artifact stemming from the fact that we are connecting the large $r$ behavior of $g(r)$ with the small $q$ behavior of $S(q)$. In any case, since we are interested in the large $r$ behavior, this factor can be disregarded.

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

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

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

[60] T. Schneider and E. Stoll, *Phys. Rev. B* **17**, 1302 (1978).

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

[62] M. L. Connolly, *J. Appl. Crystallogr.* **16**, 548 (1983).

[63] Other definitions of the PSD are possible. In the definition of Do et al. (J. Colloid Interface Sci. **328**, 110 (2008)), for example, the void-solid interface is the boundary of the region $D$ such that if the center of the probe particle is in $D$, the solid-probe interaction energy is zero. If the solid-probe interaction potential is taken to be a hard-sphere interaction, this definition is equivalent to the one of Torquato. Despite the arbitrariness deriving from the choice of the solid-probe potential, this method can indeed be more suitable than the one of Torquato in situations where the solid-probe interaction cannot be approximated as a hard-sphere interaction (e.g., in the presence of very soft potentials).

[64] B. Lee and F. M. Richards, *J. Mol. Bio.* **55**, 379 (1971).

[65] F. M. Richards, *Annu. Rev. Biophys. Bio.* **6**, 151 (1977).

[66] We note that different choices of $r_m$ will lead to different PSDs. However, given two different values $r_m$ and $r'_m$, we expect the difference between the PSDs computed using these two values to be relevant only when the average pore size is comparable to, or smaller than, $|r_m - r'_m|$.

[67] M. Abramowitz and I. A. Stegun, *Handbook of mathematical functions: with formulas, graphs, and mathematical tables* (Courier Corporation, 1965).

[68] J. H. Davenport, Y. Siret, and E. Tournier, *Computer Algebra: Systems and Algorithms for Algebraic Computation (2nd ed.)* (Academic press, 1993).

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

[70] Y. Tsuji, X. Li, and M. Shibayama, *Gels* **4**, 50 (2018).

[71] E. R. Duering, K. Kremer, and G. S. Grest, *J. Chem. Phys.* **101**, 8169 (1994).

[72] J. Wittmer, H. Meyer, J. Baschnagel, A. Johner, S. Obukhov, L. Mattioni, M. Müller, and A. N. Semenov, *Physical review letters* **93**, 147801 (2004).

[73] L. Schäfer and K. Elsner, *The European Physical Journal E* **13**, 225 (2004).

[74] H.-P. Hsu, W. Paul, and K. Binder, *Macromolecules* **43**, 3094 (2010).

[75] H.-P. Hsu and K. Kremer, *J. Chem. Phys.* **144**, 154907 (2016).

[76] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, *Numerical recipes 3rd edition: The art of scientific computing* (Cambridge university press, 2007).

[77] S. G. Johnson, “The nlopt nonlinear-optimization package,” http://github.com/stevengj/nlopt.

[78] T. H. Rowan, *Functional Stability Analysis of Numerical Algorithms*, Ph.D. thesis, Department of Computer Sciences, University of Texas at Austin (1990).

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

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

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

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