Characterizing the mesh size of polymer solutions via the pore size distribution

TL;DR

This paper introduces a novel method to characterize the mesh size of polymer solutions using pore size distribution, providing a more accurate and accessible measure than traditional correlation length methods.

Contribution

The authors propose using pore size distribution as a new way to quantify mesh size, validating it through simulations and analytical mapping to sphere systems.

Findings

Mean pore size correlates with mesh size predicted by scaling theory.

Pore size distribution can be accurately estimated from polymer radius of gyration and density.

Mapping polymer solutions to sphere systems reproduces pore size distribution results.

Abstract

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. (Ref.1) and one by Gubbins et al. (Ref.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.

Tables

Table 1. Table 1: Details of the simulated systems: chain length N 𝑁 N , number of chains N c subscript 𝑁 𝑐 N_{c} , number of monomers M 𝑀 M , monomer density ρ 𝜌 \rho , equilibration time t e subscript 𝑡 𝑒 t_{e} .

$N$	$N_c$	$M=N_{c}N$	$\rho$ (range)	$t_e$
$50$	$200$	$1\cdot {10}^4$	$[0.001,1.00]$	$1.5\cdot {10}^5$
$200$	$50$	$1\cdot {10}^4$	$[0.001,0.85]$	$4.5\cdot {10}^5$
$1000$	$10$	$1\cdot {10}^4$	$[0.001,0.26]$	$3.0\cdot {10}^6$
$1000$	$50$	$5\cdot {10}^4$	$[0.001,0.26]$	$3.0\cdot {10}^6$

Table 2. Table 2: N 𝑁 N -dependent properties of the simulated systems: Radius of gyration in the dilute limit, R g ​ 0 subscript 𝑅 𝑔 0 R_{g0} and approximate overlap concentration (Eq. ( 20 )).

$N$	$R_{g0}$	$3N/4\pi R_{g0}^3$
$50$	$4.83\pm 0.02$	$1.1\cdot {10}^{-1}$
$200$	$11.44\pm 0.07$	$3.2\cdot {10}^{-2}$
$1000$	$29.9\pm 0.2$	$8.9\cdot {10}^{-3}$

Equations

$$\rho \approx \frac{N_b}{{\xi }_b^3},$$

$${\xi }_b\approx b{N}_b^{\nu },$$

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

$${\xi }_b\approx \{\begin{array}{ll}{R}_{g0}& \rho <{\rho }^{\ast }\\ b(b^3\rho {)}^{-\nu /(3\nu -1)}& {\rho }^{\ast }<\rho <{\rho }^{\ast \ast }\\ b& \rho >{\rho }^{\ast \ast },\end{array}$$

$${\rho }^{\ast }\approx N{R}_{g0}^{-3}\approx N^{1-3\nu }{b}^{-3}$$

$${\rho }^{\ast \ast }\approx b^{-3}\approx {\rho }_{\text{melt}}.$$

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

$$\frac{{\xi }_b}{R_{g0}}\approx \{\begin{array}{ll}1& \chi <{\chi }^{\ast }\\ {\chi }^{-\nu }& {\chi }^{\ast }<\chi <{\chi }^{\ast \ast }\\ N^{-\nu }& \chi >{\chi }^{\ast \ast }.\end{array}$$

$$g(r)=1+\frac{A}{r}{e}^{-r/{\xi }_c},$$

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

$$S(q)=1+4\pi \rho {\int }_0^{\infty }[g(r)-1]\frac{\sin (qr)}{qr}{r}^{2}dr.$$

$$S(q)=\frac{S(0)}{1+(q{\xi }_c{)}^2},$$

$$S_1(q)\simeq N\left[1-{\left(\frac{q{R}_g}{\sqrt{3}}\right)}^2\right]\simeq \frac{N}{1+(q{R}_g/\sqrt{3}{)}^2},$$

$${\xi }_c\approx b(b^3\rho {)}^{-\nu /(3\nu -1)}({\rho }^{\ast }<\rho <{\rho }^{\ast \ast }),$$

$${\xi }_c\approx \{\begin{array}{ll}{R}_{g0}& \rho <{\rho }^{\ast }\\ b(b^3\rho {)}^{-\nu /(3\nu -1)}& {\rho }^{\ast }<\rho <{\rho }^{\ast \ast }\\ h(\rho )& \rho >{\rho }^{\ast \ast },\end{array}$$

$$\xi \approx {\xi }_b\approx {\xi }_c({\rho }^{\ast }<\rho <{\rho }^{\ast \ast }).$$

$$\xi \approx \{\begin{array}{ll}(\rho /N{)}^{-1/3}& \rho <{\rho }^{\ast }\\ b(b^3\rho {)}^{-\nu /(3\nu -1)}& {\rho }^{\ast }<\rho <{\rho }^{\ast \ast }\\ f(\rho )& \rho >{\rho }^{\ast \ast },\end{array}$$

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

$$\xi \simeq R_{g0}{\left(\frac{\rho }{{\rho }^{\ast }}\right)}^{-\nu /(3\nu -1)},$$

$${\rho }^{\ast }\simeq \frac{3N}{4\pi R_{g0}^3}.$$

$$\mathcal{U}(r)=\{\begin{array}{ll}4ϵ\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{array}$$

$${\mathcal{U}}_{\text{bond}}(r)=-\frac{k{r}_0^2}{2}\ln [1-(r/r_0{)}^2],$$

$$m\ddot{\mathbf{r}}=-\nabla {\mathcal{U}}_{\text{tot}}(\{\mathbf{r}\})-m\Gamma \dot{\mathbf{r}}+\sqrt{2m\Gamma k_{B}T}\mathbf{\zeta }(t),$$

$$S(\mathbf{q})\equiv \frac{1}{M}\sum _{k,j=1}^M⟨\exp [-i\mathbf{q}\cdot ({\mathbf{r}}_k-{\mathbf{r}}_j)]⟩,$$

$$S(q)\simeq S_1(q)=\{\begin{array}{ll}N/(1+q^2{R}_{g0}^2/3)& q<2\pi /R_{g0}\\ A{q}^{-1/\nu }& 2\pi /R_{g0}<q<2\pi /b\\ \mathcal{O}(1)& q>2\pi /b,\end{array}$$

$$S(q)\simeq \{\begin{array}{ll}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{array}$$

$$g(r)\equiv \frac{1}{4\pi M\rho r^2}\sum _{\begin{array}{c}k=1\\ j\ne k\end{array}}^M⟨\delta (|\mathbf{r}+{\mathbf{r}}_k-{\mathbf{r}}_j|)⟩.$$

$$F_T(0)=1;\underset{r\to \infty }{\lim }{F}_T(r)=0.$$

$$P_T(r)\equiv -\frac{d{F}_T(r)}{dr}.$$

$$F_T(r)={\int }_r^{\infty }{P}_T(x)dx,$$

Full text

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 De Gennes (1979); Rubinstein and Colby (2003). 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$ De Gennes (1979).

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 Zhou et al. (2008); Zhou (2013) and medicine, e.g. for drug delivery Soppimath et al. (2001); Lai et al. (2007); Schuster et al. (2013). Many theoretical approaches have been proposed to describe the diffusion of particles in polymer solutions and gels, such as geometric obstruction models Ogston et al. (1973); Altenberger and Tirrell (1984); Johansson et al. (1991), hydrodynamic models Phillies et al. (1985); Phillies (1986); Cukier (1984); Odijk (2000); Fan et al. (2007); Tuinier and Fan (2008) and mode-coupling theory Egorov (2011); Yamamoto and Schweizer (2011); Dong et al. (2015) (for a review of some of these models, see Ref. 22). Although the predictions of different 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 Cai et al. (2011). This model, which seems to be the one in best agreement with recent experiments Kohli and Mukhopadhyay (2012); Babaye Khorasani et al. (2014); Poling-Skutvik et al. (2015) and simulations Chen et al. (2018, 2019), 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 Lustig and Peppas (1988); Masaro and Zhu (1999); Dell and Schweizer (2013); Rose et al. (2013); Cai et al. (2015); Parrish et al. (2017, 2018).

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 De Gennes (1979), 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 De Gennes (1979), or measurements of the monomer density fluctuation correlation length, $\xi_{c}$ De Gennes (1979); Rubinstein and Colby (2003); Koshy et al. (2003). 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 Koshy et al. (2003).

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 Scheidegger (1957); Prager (1963); Torquato et al. (1990); Torquato and Avellaneda (1991); Torquato (1995, 2013); Gelb and Gubbins (1999); Thomson and Gubbins (2000); Pikunic et al. (2003); Bhattacharya and Gubbins (2006). 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 Parrish et al. (2017), 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 concepts of blob De Gennes (1979); Rubinstein and Colby (2003), monomer density fluctuation correlation length De Gennes (1979); Rubinstein and Colby (2003); Koshy et al. (2003) and geometrical mesh size. Most of what is discussed in this section can be found in standard texts De Gennes (1979); Doi and Edwards (1988); Teraoka (2002); Rubinstein and Colby (2003). 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 Rubinstein and Colby (2003). This choice is also motivated by the fact that the model chosen for the simulations is athermal, in the sense that all the nonbonded 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\simeq 0.588$ Rubinstein and Colby (2003); Clisby (2010).

II.1 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 De Gennes (1979). 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 De Gennes (1979); Rubinstein and Colby (2003). 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 Rubinstein and Colby (2003)

[TABLE]

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 De Gennes (1979); Rubinstein and Colby (2003):

[TABLE]

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

[TABLE]

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 De Gennes (1979); Rubinstein and Colby (2003). If the chains are ideal at all length scales $r>b$, the solution is called concentrated and $\xi_{b}\approx b$ Rubinstein and Colby (2003); Teraoka (2002) 111We 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$.. 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:

[TABLE]

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 De Gennes (1979); Rubinstein and Colby (2003), we obtain

[TABLE]

and

[TABLE]

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 Paul et al. (1991); Gerroff et al. (1993) 222We choose to use $\chi$ and not $\rho/\rho^{*}$ as a scaling variable since $\rho^{*}$ is not a sharply defined quantity.

[TABLE]

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

[TABLE]

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.

II.2 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)$ Hansen and McDonald (1990), 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 Fisher (1964); De Gennes (1979); Doi and Edwards (1988)

[TABLE]

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 Verlet (1968); Evans et al. (1994, 1993); Koshy et al. (2003),

[TABLE]

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)$ Hansen and McDonald (1990); Rubinstein and Colby (2003), which is related to $g(r)$ through a Fourier transform Hansen and McDonald (1990):

[TABLE]

From Eqs. (9) and (11), we obtain the low $q$ behavior of $S(q)$, which is given by the Ornstein-Zernike function Fisher (1964); De Gennes (1979); Doi and Edwards (1988) 333To be more precise, Eq. (12) follows formally from Eq. (11) if $g(\mathbf{r})=1+Ae^{-r/\xi_{c}}/r-\delta(\mathbf{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.,

[TABLE]

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)$ Teraoka (2002); Rubinstein and Colby (2003). Since we know that for $q<2\pi/R_{g}$ Teraoka (2002); Rubinstein and Colby (2003)

[TABLE]

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 Rubinstein and Colby (2003); Teraoka (2002)

[TABLE]

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 densityKoshy et al. (2003). Therefore one can summarize the behavior of $\xi_{c}$ as follows:

[TABLE]

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

II.3 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):

[TABLE]

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$:

[TABLE]

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

[TABLE]

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))

[TABLE]

where Doi and Edwards (1988); Teraoka (2002)

[TABLE]

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 Kremer and Grest (1990). The excluded volume interactions between all monomers are modeled as a purely repulsive Weeks-Chandler-Andersen (WCA) potential Weeks et al. (1971),

[TABLE]

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

[TABLE]

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$ Kremer and Grest (1990). Since non-bonded interactions are purely repulsive, this model mimics the behavior of polymers in an athermal solvent Koshy et al. (2003); Rubinstein and Colby (2003).

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_{c}N=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. 1. 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 Plimpton (1995). The temperature was kept constant at $T=1.0$ by means of a Langevin thermostat Schneider and Stoll (1978), so that the force experienced by a monomer is

[TABLE]

where $\mathbf{r}$ is the position vector, and $\mathcal{U}_{\text{tot}}(\{\mathbf{r}\})=\mathcal{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 Auhl et al. (2003)). 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 Rubinstein and Colby (2003), 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. 1.

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.

IV.1 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 Hansen and McDonald (1990)

[TABLE]

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 Rubinstein and Colby (2003); Teraoka (2002)

[TABLE]

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. SII 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 Rubinstein and Colby (2003)

[TABLE]

where $B>0$ is a constant. From our data at low $\rho$, we can observe at intermediate $q$ the slope $-1/\nu\simeq-1.70$; 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\geq 0.85$, the structure becomes virtually indistinguishable from that of a dense liquid Hansen and McDonald (1990); Verlet (1968). Comparing $S(q)$ for different 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$ Rubinstein and Colby (2003).

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 Hansen and McDonald (1990)

[TABLE]

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 Koshy et al. (2003). 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 two estimates give almost identical results for $\rho<0.2$. For $\rho\geq 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 Koshy et al. (2003).

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$.

IV.2 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) Scheidegger (1957); Prager (1963); Torquato (2013); Torquato et al. (1990); Torquato and Avellaneda (1991); Torquato (1995); Gelb and Gubbins (1999); Thomson and Gubbins (2000); Pikunic et al. (2003); Bhattacharya and Gubbins (2006).

We will consider two different definitions of the PSD: The first one, although introduced by others Scheidegger (1957); Prager (1963), was formalized and extensively studied by Torquato and coworkers Torquato (2013); Torquato et al. (1990); Torquato and Avellaneda (1991); Torquato (1995), while the second one was introduced by Gubbins and coworkers Gelb and Gubbins (1999); Thomson and Gubbins (2000); Pikunic et al. (2003); Bhattacharya and Gubbins (2006). 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 Torquato (2013); Torquato et al. (1990); Torquato and Avellaneda (1991); Torquato (1995). 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)\leq 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

[TABLE]

Torquato’s PSD is then defined as the probability density function Torquato (2013); Torquato et al. (1990); Torquato and Avellaneda (1991); Torquato (1995)

[TABLE]

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

[TABLE]

and that $P_{T}(r)$ is normalized to unity:

[TABLE]

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 Gelb and Gubbins (1999) (Fig. 6). The surface $S_{G}$ in which $V_{G}$ is enclosed is sometimes called Connolly surface or reentrant surface Connolly (1983); Thomson and Gubbins (2000). 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$ Gelb and Gubbins (1999) 444Other 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 $\mathcal{D}$ such that if the center of the probe particle is in $\mathcal{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).. 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 Lee and Richards (1971); Richards (1977); Connolly (1983); Thomson and Gubbins (2000).

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 Prager (1963):

[TABLE]

with $r\leq R$, while

[TABLE]

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.,

[TABLE]

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$ 555We note that different choices of $r_{m}$ will lead to different PSDs. However, given two different values $r_{m}$ and $r^{\prime}_{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^{\prime}_{m}|$.. 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 $\operatorname{sgn}(\tilde{r})\cdot{\tilde{r}}^{2}$ (with $\operatorname{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 better appreciated by considering the relative standard deviation,

[TABLE]

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.4. 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.

IV.3 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 with the theoretical expectations for $\xi$. The most natural choice is to consider the mean pore radius,

[TABLE]

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\left(\rho/\rho^{*}\right)^{-{\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). 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$. 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_{\text{max}}$ at which $P(r_{\text{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_{\text{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.

IV.4 $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 Torquato (2013); Torquato et al. (1990); Torquato and Avellaneda (1991); Torquato (1995), 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 Torquato (2013); Torquato et al. (1990); Torquato and Avellaneda (1991); Torquato (1995)

[TABLE]

with

[TABLE]

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 Torquato and Avellaneda (1991); Torquato (2013)

[TABLE]

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}$ Torquato and Avellaneda (1991); Torquato (2013). For small values of $\eta$, Eq. (37) reduces to Eq. (39), i.e., the HS model is equivalent to the OS model Torquato and Avellaneda (1991); Torquato (2013).

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{\tiny{HS}}}(\rho,N)$ (resp. $R_{\text{\tiny{OS}}}(\rho,N)$) and density $\rho_{\text{\tiny{HS}}}(\rho,N)$ (resp. $\rho_{\text{\tiny{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{\tiny{HS}}}(r;R_{\text{\tiny{HS}}},\rho_{\text{\tiny{HS}}})$, and with the OS model $P_{T}^{\text{\tiny{OS}}}(r;R_{\text{\tiny{OS}}},\rho_{\text{\tiny{OS}}})$. It is virtually impossible to distinguish $P_{T}^{\text{\tiny{HS}}}(r;R_{\text{\tiny{HS}}},\rho_{\text{\tiny{HS}}})$ from $P_{T}^{\text{\tiny{OS}}}(r;R_{\text{\tiny{OS}}},\rho_{\text{\tiny{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{\tiny{HS}}},R_{\text{\tiny{OS}}}$ and $\rho_{\text{\tiny{HS}}}/\rho,\rho_{\text{\tiny{OS}}}/\rho$ as function of the monomer density. Through an extrapolation of the low density behavior, we deduce that in the limit $\rho\to 0$, $R_{\text{\tiny{HS}}},R_{\text{\tiny{OS}}}\to R_{g0}(N)$, and $\rho_{\text{\tiny{HS}}},\rho_{\text{\tiny{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:

[TABLE]

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

[TABLE]

as a function of $\rho$. We note that both $\eta_{\text{\tiny{HS}}}$ and $\eta_{\text{\tiny{OS}}}$ depend basically only on $(\rho,N)$, with $\eta_{\text{\tiny{OS}}}$ being slightly smaller than $\eta_{\text{\tiny{HS}}}$. At low $\rho$, both $\eta_{\text{\tiny{OS}}}$ and $\eta_{\text{\tiny{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

[TABLE]

i.e., the dimensionless density the chains would have were they spheres of radius $R_{g}$. Importantly, we observe that when $N$ increases, $\eta_{\text{\tiny{OS}}}$ and $\eta_{\text{\tiny{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

[TABLE]

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

[TABLE]

where $E_{n}(\eta)$ is an exponential integral Abramowitz and Stegun (1965),

[TABLE]

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

[TABLE]

In Fig. 17, we compare the measured average pore size $\langle r\rangle_{T}$ with the OC model result $\langle r\rangle_{T}^{\text{\tiny{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 De Gennes (1979); Rubinstein and Colby (2003); Teraoka (2002)

[TABLE]

This implies, given the definition of $\eta_{c}$, that

[TABLE]

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 Davenport et al. (1993); the leading terms are

[TABLE]

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

[TABLE]

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

[TABLE]

We thus find that for very long chains, ${\langle r\rangle_{T}^{\text{\tiny{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}^{\text{\tiny{OC}}}$ will always break down for $\rho\geq\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}^{\text{\tiny{OC}}}$ are so similar even if $P_{T}(r)$ and $P_{T}^{\text{\tiny{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.3), 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 Rintoul et al. (1996).

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 Torquato (2013) and by Gubbins Gelb and Gubbins (1999). 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 Parrish et al. (2017).

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$ Tsuji et al. (2018). 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 Duering et al. (1994). 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 Doi and Edwards (1988), since on timescales shorter than the entanglement relaxation time the entanglements effectively act as crosslinks Dell and Schweizer (2013); Cai et al. (2015).

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 Chen et al. (2019). 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

[TABLE]

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 De Gennes (1979); Rubinstein and Colby (2003); Teraoka (2002)

[TABLE]

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. 2, 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 Wittmer et al. (2004)

[TABLE]

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 Schäfer and Elsner (2004); Wittmer et al. (2004); Hsu et al. (2010):

[TABLE]

where $l_{p}$ is the persistence length Rubinstein and Colby (2003), $l_{b}=\sqrt{\langle\mathbf{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}$ Hsu and Kremer (2016). 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 Hsu et al. (2010), 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

AII.1 Torquato’s PSD

The algorithm to compute Torquato’s PSD is described in Ref. 40. First of all, as mentioned in Sec. IV.2, 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:

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

The smallest distance $r_{\text{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_{\text{min}}-\sigma/2$.

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

AII.2 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

[TABLE]

subject to the constraint

[TABLE]

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 Press et al. (2007). We have used the open-source Sbplx algorithm of the NLopt library Johnson , which is a re-implementation of Subplex Rowan (1990). Other choices are possible, although one has to make sure that the chosen algorithm can handle discontinuous objective functions Bhattacharya and Gubbins (2006).

Acknowledgements.

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 Hsu and Kremer (2016)

[TABLE]

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

[TABLE]

where $\mathbf{r}_{\text{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 Hsu and Kremer (2016).

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 Doi and Edwards (1988); Rubinstein and Colby (2003) gives a good approximation of chain dynamics at all densities, whereas in the presence of entanglements, the dynamics is described by the reptation model Doi and Edwards (1988); Rubinstein and Colby (2003). In both cases, when $t>t_{\text{rel}}$, with $t_{\text{rel}}$ the longest relaxation time of the system, we expect $g_{1}(t)=g_{3}(t)\propto t$ Doi and Edwards (1988); Hsu and Kremer (2016): 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 $\simeq 85$ at $\rho=0.85$ ($T=1.0$) Hoy et al. (2009), and it decreases with decreasing $\rho$ Fetters et al. (1994, 1999); Kröger and Hess (2000). 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 Doi and Edwards (1988); Rubinstein and Colby (2003). 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)\simeq g_{3}(t)$, and thus we conclude that the equilibration time $t_{e}$ (see Tab. 1 in the main text) is larger than the longest relaxation time of the system, $t_{\text{rel}}$, and therefore all the systems are well equilibrated.

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 $\rho$.

We recall that the theoretical prediction for $S_{1}(q)/N$ is Rubinstein and Colby (2003); Teraoka (2002); Paul et al. (1991)

[TABLE]

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}$ De Gennes (1979); Rubinstein and Colby (2003). 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 $\rho$ 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. 2 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 Rubinstein and Colby (2003).

SIV Bond angle distribution

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

[TABLE]

where $\mathbf{r}_{i,j}\equiv\mathbf{r}_{i}-\mathbf{r}_{j}$ and $\mathbf{r}_{i},\ i=2,\dots,N-1$ is the monomer’s position vector. We have therefore $\theta=\pi-\theta_{1}$, where $\theta_{s}$ is defined in Eq. (54) of the main text.

At low density, $P(\theta)$ shows a maximum at $\theta\simeq 112^{\circ}$ and falls to zero rather sharply at $\theta\simeq 60^{\circ}$ because of the excluded volume interaction. When $\rho$ is increased past $\rho\simeq 0.3$, the shape of $P(\theta)$ starts to change significantly in that it develops a peak at $\theta\simeq 70^{\circ}$, signaling that the chains are compressed. Overall, the average bond angle decreases by $\simeq 8\%$ when going from the dilute regime to density $\rho\gtrsim 0.85$. This effect is also observable in the average bond length (not shown), although the decrease is in this case only $\simeq 1\%$. Although it is possible to conceive chain conformations with a small $\langle\theta\rangle$ but a high persistence length (e.g., an “accordion-like” rigid polymer), the fact that up to very high $\rho$ the form factor is compatible with the scaling predictions for chains in a concentrated solution (Sec. SII) suggests that the reduction of $\langle\theta\rangle$ corresponds to a reduction of the persistence length, as also observed from the bond-bond correlation function (Sec. AI in the main text).