# Effect of Chain Polydispersity on the Elasticity of Disordered Polymer Networks

Valerio Sorichetti, Andrea Ninarello, José M Ruiz-Franco, Virginie

Hugouvieux, Walter Kob, Emanuela Zaccarelli, Lorenzo Rovigatti

► To cite this version:

Valerio Sorichetti, Andrea Ninarello, José M Ruiz-Franco, Virginie Hugouvieux, Walter Kob, et al. Effect of Chain Polydispersity on the Elasticity of Disordered Polymer Networks. Macromolecules, 2021, 54 (8), pp.3769 - 3779. 10.1021/acs.macromol.1c00176. [https://doi.org/10.1021/acs.macromol.1c00176](https://doi.org/10.1021/acs.macromol.1c00176)

HAL Id: hal-03234741

https://hal.inrae.fr/hal-03234741

Submitted on 25 May 2021

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>
# Effect of Chain Polydispersity on the Elasticity of Disordered Polymer Networks

Valerio Sorichetti,*,¶ Andrea Ninarello,*,¶ José M. Ruiz-Franco, Virginie Hugouvieux, Walter Kob, Emanuela Zaccarelli, and Lorenzo Rovigatti*

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

Cite This: *Macromolecules* 2021, 54, 3769–3779

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

Read Online

## ACCESS

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

Metrics & More

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

Article Recommendations

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

Supporting Information

**ABSTRACT:** Due to their unique structural and mechanical properties, randomly cross-linked polymer networks play an important role in many different fields, ranging from cellular biology to industrial processes. In order to elucidate how these properties are controlled by the physical details of the network (e.g., chain-length and end-to-end distributions), we generate disordered phantom networks with different cross-linker concentrations C and initial densities $ρ_{init}$ and evaluate their elastic properties. We find that the shear modulus computed at the same strand concentration for networks with the same C, which determines the number of chains and the chain-length distribution, depends strongly on the preparation protocol of the network, here controlled by $ρ_{init}$. We rationalize this dependence by employing a generic stress–strain relation for polymer networks that does not rely on the specific form of the polymer end-to-end distance distribution. We find that the shear modulus of the networks is a nonmonotonic function of the density of elastically active strands, and that this behavior has a purely entropic origin. Our results show that if short chains are abundant, as it is always the case for randomly cross-linked polymer networks, the knowledge of the exact chain conformation distribution is essential for correctly predicting the elastic properties. Finally, we apply our theoretical approach to literature experimental data, qualitatively confirming our interpretations.

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

## 1. INTRODUCTION

For many applications, the elasticity of a cross-linked polymer network is one of its most important macroscopic properties.¹ It is thus not surprising that a lot of effort has been devoted to understanding how the features of a network, such as the fraction and functionality of cross-linkers or the details of the microscopic interactions between chain segments, contribute to generate its elastic response.²⁻⁶ The macroscopic behavior of a real polymer network (be it a rubber or a hydrogel) depends on many quantities, such as the properties of the polymer and of the solvent, the synthesis protocol, and the thermodynamic parameters. However, in experiments, it is difficult to disentangle how these different elements contribute to the elastic properties of the material. This task becomes easier in simulations because all the relevant parameters can be controlled in detail.⁷⁻¹⁶ In this regard, an important feature of real polymer networks that can be exploited is that their elasticity can be described approximately as the sum of two contributions: one due to the cross-linkers and one due to the entanglements.¹⁵⁻¹⁸ The former can be approximated well by the elastic contribution of the corresponding phantom network,¹⁹ that is, when the excluded volume between the strands is not taken into account.¹⁵,¹⁶ It is, therefore, very important to understand the role that the chain conformation

distribution plays in determining the dynamics and elasticity of phantom polymer models.

The distribution of the chemical lengths of the strands connecting any two cross-linkers in a network, that is, the chains (chain-length distribution for short), depends on the chemical details and on the synthesis protocol. In randomly cross-linked networks, this distribution is typically exponential,⁷,²⁰ whereas monodisperse or quasimonodisperse networks can be obtained by using specific end-linking protocols, for example, via the assembly of tetra-PEG macromers with a small polydispersity.²¹ Regardless of the synthesis route, the presence of short or stretched chains is common, although the exact form of the chain conformation fluctuations is highly nontrivial. From a theoretical viewpoint, however, the majority of the results on the elasticity of polymer networks have been obtained within the mean-field realm, in which scaling assumptions and chain Gaussianity are assumed.¹⁹,²² There-

Received: January 25, 2021

Revised: March 20, 2021

Published: April 14, 2021

<figure><img src="image_8.png" /></figure>
fore, simulations can be extremely helpful to clarify the exact role played by the chain-length distribution and better understand the experimental results. However, most simulation studies have focused on melt densities, where random or end-cross-linking can be carried out efficiently,⁷⁻¹⁰,¹²,¹⁴⁻¹⁶ or have employed idealized lattice networks.¹¹,¹³,²³⁻²⁵ This makes it challenging to compare the results from such simulations with common experimental systems such as hydrogels, which are both low-density and disordered.²⁶

In the present paper, we show that the knowledge of the exact chain end-to-end distribution is essential to correctly predict the linear elastic response of low-density polymer networks. We do so by simulating disordered phantom networks generated with different cross-linker concentrations C and initial monomer densities $\rho_{\text{init}}$. In our systems, the former parameter controls the number of chains and the chain-length distribution, while the latter determines the initial end-to-end distance distribution of the chains and, therefore, plays a similar role as the solvent quality in an experimental synthesis. To generate the gels, we exploit a recently introduced technique based on the self-assembly of patchy particles, which has been found to be correctly reproduce structural properties of experimental microgels.²⁷⁻²⁹ This method allows us to obtain systems at densities comparable with those of experimental hydrogels, that is, giving access to swelling regimes inaccessible through the previously employed techniques based on numerical vulcanization of high-density polymer melts.⁸⁻¹³,¹⁵,¹⁶ We first demonstrate that systems generated with the same C but at different values of $\rho_{\text{init}}$ can display very different elastic properties even when probed at the same strand concentration, despite having the same chain length distribution. Second, we compare the numerical results to the phantom network theory.¹⁹ In order to do so, we determine the theoretical relation between the shear modulus G and the single-chain entropy for generic non-Gaussian chains. We find a good agreement between theory and simulation only for the case in which the exact chain end-to-end distribution is given as an input to the theory, with some quantitative deviations appearing at low densities. On the other hand, assuming a Gaussian behavior of the chains leads to qualitatively wrong predictions for all the investigated systems except the highest density ones. Overall, our analysis shows that for low-density polymer networks and in the presence of short chains, the knowledge of the exact chain conformational fluctuations is crucial to predict the system elastic properties reliably. Notably, we validate our approach against recently published experimental data,²¹,³⁰ showing that the behavior of systems where short chains are present cannot be modeled without precise knowledge of the chain-size-dependent end-to-end distribution.

## 2. THEORETICAL BACKGROUND

In this section, we review some theoretical results on the elasticity of polymer networks, for the most part available in the literature,¹,¹⁹,²² by reorganizing them and introducing the terminology and notation that will be employed in the rest of the paper. We will consider a polydisperse polymer network made of cross-linkers of valence $\phi$ connected by $N_s$ strands. Here and in the following, we will assume the network to be composed of $N_s$ elastically active strands, defined as strands with the two ends connected to distinct cross-linkers, that is, that are neither dangling ends nor closed loops (i.e., loops of order one). Moreover, for those strands, which are part of

higher-order loops, we assume their elasticity to be independent of the loop order (see Zhong et al.³¹ and Lin et al.³²). We will focus on evaluating the shear modulus G of the gel, which relates a pure shear strain to the corresponding stress in the linear elastic regime.³³ One can theoretically compute G by considering uniaxial deformations of strain $\lambda$ along, for instance, the x axis. We assume the system to be isotropic; moreover, since we are interested in systems with no excluded volume interactions, we assume a volume-preserving transformation,ª that is, $\lambda_x = \lambda$ and $\lambda_y = \lambda_z = \lambda^{-1/2}$ as extents of deformation along the three axes.

The starting point to calculate the shear modulus is the single-chain entropy, which is a function of the chain's end-to-end distance.²² In general, we can write the instantaneous end-to-end vector of a single chain, which connects any two cross-linkers as $\mathbf{r}(t) = \mathbf{R} + \mathbf{u}(t)$, where $\mathbf{R} \equiv \mathbf{r}(t)$ represents the time-averaged end-to-end vector and $\mathbf{u}(t)$ the fluctuation term (see Figure 1 for a cartoon depicting these quantities). We also

<figure><img src="image_9.png" /><figcaption>**Figure 1.** Cartoon providing a visual explanation of some of the quantities used throughout the text. In the cartoon, the black dots are the average positions of two cross-linkers i and j, the green dots are their instantaneous positions at time t, and the gray dots are their positions at some other times. The time-averaged end-to-end vector **R** is the vector connecting the black dots, the instantaneous end-to-end vector **r**(t) is the vector connecting the green dots, and its instantaneous fluctuation is the difference between the instantaneous fluctuations of the positions of the two cross-linkers, **u**ᵢ(t) and **u**ⱼ(t).</figcaption></figure>

assume that there are no excluded volume interactions, so that the chains can freely cross each other. We thus have $\overline{r^2} = R^2 + \overline{u^2}$,ᵇ because $\overline{\mathbf{R} \cdot \mathbf{u}(t)} = 0$, the position and fluctuations of cross-linkers being uncorrelated.¹⁹

The entropy of a chain with end-to-end vector $\mathbf{r} = (r_{x0}, r_{y0}, r_{z0})$ is $S_n(\mathbf{r}) = k_B \log W_n(\mathbf{r}) + A_n$,³⁴ where $W_n(\mathbf{r})$ is the end-to-end probability density of $\mathbf{r}$ and $A_n$ is a temperature-dependent parameter that can be set to zero in this context. If the three spatial directions are independent (which is the case, e.g., if $W_n(\mathbf{r})$ is Gaussian) then $W_n(\mathbf{r})$ can be written as the product of three functions of $r_{x0}, r_{y0}$ and $r_{z0}$, so that $S_n(\mathbf{r}) = s_n(r_x) + s_n(r_y) + s_n(r_z)$, where $s_n$ is the entropy of a one-dimensional chain. Building upon this result, we can assume that each chain in the network can be replaced by three independent one-dimensional chains parallel to the axes using the so-called three-chain approximation.¹,³⁵ This assumption is exact for Gaussian chains, although for non-Gaussian chains the associated error is small if the strain is not too large.¹

We will also assume (i) that the length of each chain in the unstrained state ($\lambda = 1$) is $\tilde{r} \equiv (\overline{r^2})^{1/2} \equiv (R^2 + \overline{u^2})^{1/2}$, and (ii) that, upon deformation, the chains deform affinely with the network, so that the length of the chain oriented along the x axis becomes $\tilde{r}_{\lambda}$ and those of the chains oriented along the y and z axes become $r_{\lambda}^{-1/2}$. With those assumptions, the single-chain entropy $S_n(\lambda)$ becomes¹
$$S_n(\lambda) = \frac{s_n(\tilde{r}_\lambda) + 2s_n(\tilde{r}_\lambda^{-1/2})}{3} \quad (1)$$

where we need to divide by three because we are replacing each unstrained chain with end-to-end distance $\tilde{r}$ by three fictitious chains of the same size. Usually, the $\lambda$-dependence of $\tilde{r}_\lambda$ is controlled by the microscopic model and by the macroscopic conditions (density, temperature, etc.). Two well-known limiting cases are the affine network model,¹⁹ in which both the average positions and fluctuations of the cross-linkers deform affinely, $\tilde{r}_\lambda = \lambda\tilde{r}$, and the phantom network model,¹⁹ in which the fluctuations are independent of the extent of the deformation, so that

$$\tilde{r}_{\lambda} = [(\lambda^2 R^2 + \overline{u^2})]^{1/2}, \quad \text{and thus} \quad \tilde{r}_{\lambda^{-1/2}} = [(R^2/\lambda + \overline{u^2})]^{1/2} \quad (2)$$

The free-energy difference between the deformed and undeformed state of a generic chain is $\Delta F = -T[S_n(\lambda) - S_n(1)]$, and thus the $x$ component of the tensile force is given by

$$f_x(\lambda) = \frac{1}{L_{x0}} \frac{d\Delta F}{d\lambda} = - \frac{T}{L_{x0}} \frac{dS_n(\lambda)}{d\lambda} \quad (3)$$

The latter quantity divided by the section $L_{y0}L_{z0}$ yields the $xx$ component of the stress tensor, which thus reads

$$\sigma_{xx} = -\frac{TR^2}{3V} \left[ \frac{\lambda}{\tilde{r}_\lambda} \frac{ds_n(\tilde{r}_\lambda)}{d\tilde{r}_\lambda} - \frac{1}{\lambda^2 \tilde{r}_\lambda^{-1/2}} \frac{ds_n(\tilde{r}_\lambda^{-1/2})}{d\tilde{r}_\lambda^{-1/2}} \right] \quad (4)$$

where we have used eq 2.

Because the volume is kept constant, the Poisson ratio is $1/2$ and hence the single-chain shear modulus $g$ is connected to the Young's modulus $Y = \frac{d\sigma_{xx}}{d\lambda}|_{\lambda=1}$ by $g = Y/3$, which implies that

$$g = -\frac{TR^2}{3V} \left[ \frac{ds_n(\tilde{r})}{d\tilde{r}} \left( \frac{1}{\tilde{r}} - \frac{R^2}{2\tilde{r}^3} \right) + \frac{d^2 s_n(\tilde{r})}{d\tilde{r}^2} \frac{R^2}{2\tilde{r}^2} \right] \quad (5)$$

We note that, although similar equations can be found in Smith³⁵ and Treloar,¹ to the best of our knowledge eq 5 has not been reported in the literature in this form. In order to obtain the total shear modulus $G$ of the network, and under the assumption that the effect of higher-order loops can be neglected,³¹,³² one has to sum over the $N_s$ elastically active chains. Of course, the result will depend on the specific form chosen for the entropy $s_n$. We stress that a closed-form expression of the end-to-end probability density $W_n(\mathbf{r})$ is not needed because only its derivatives play a role in the calculation. Hence, it is sufficient to know the force-extension relation for the chain, because, as discussed above, the component of the force along the pulling direction satisfies eq 3 (see also Appendix A).

For a freely jointed chain (FJC)²² of $n$ bonds of length $b$, $W_n(\mathbf{r})$ has the following form¹,³⁶

$$W_n(\mathbf{r}) = \left[ \frac{n(n-1)}{8\pi r b^2} \right] \sum_{t=0}^{\tau} \frac{(-1)^t}{t!(n-t)!} \left[ \left( \frac{nb-r}{2b} \right) - t \right]^{n-2} \quad (6)$$

where $\tau = [(nb - r)/2b]$, that is, the largest integer smaller than $(nb - r)/2b$.

In the limit of large $n$, eq 6 reduces to a Gaussian³⁶

$$W_n^G(\mathbf{r}) = \left(\frac{3}{2\pi n b^2}\right)^{3/2} \exp\left(-\frac{3r^2}{2nb^2}\right) \quad (7)$$

Under this approximation, the shear modulus takes the well-known form

$$G^G = \frac{k_B T}{V} \sum_i^{N_s} \frac{R_i^2}{n_i b^2} = \left\langle \frac{R^2}{nb^2} \right\rangle k_B T \nu \equiv A k_B T \nu \quad (8)$$

where $\nu = N_s/V$ is the number density of elastically active strands and $A$ is often called the front factor.³⁵,³⁷⁻³⁹ We have also introduced the notation $\langle \cdot \rangle = N_s^{-1} \sum_i^{N_s} \cdot$ for the average over all the strands in the system. In the particular case that the $\tilde{r}_i^2$ values of the different chains are Gaussian-distributed (a distinct assumption from the one that $W_n(\mathbf{r})$ is Gaussian), which is the case, for example, for end-cross-linking starting from a melt of precursor chains,¹⁰ it can be shown that $A = 1 - \frac{2}{\phi}$ (we recall that $\phi$ is the cross-linker valence), so that one obtains the commonly reported expression (see also Appendix A)¹⁹,²²

$$G^G = \left(1 - \frac{2}{\phi}\right) k_B T \nu \quad (9)$$

Equation 8 was derived from eq 5, which assumes the validity of the phantom network model. If one assumes, on the other hand, that the affine network model is valid, a different expression for $G$ is obtained (see Supporting Information).

To obtain a more accurate description of the end-to-end probability distribution for strained polymer networks, one has to go beyond the Gaussian model and introduce more refined theoretical assumptions. Among other approaches, the Langevin-FJC¹ (L-FJC), the extensible-FJC⁴⁰ (ex-FJC), and the worm-like chain⁴¹ (WLC) have been extensively used in the literature. In the L-FJC model, the force-extension relation is approximated using an inverse Langevin function, whereas in the ex-FJC model bonds are modeled as harmonic springs. These models give a better description of the system's elasticity when large deformations are considered. The WLC model, in which chains are represented as continuously flexible rods, is useful when modeling polymers with a high persistence length (compared to the Kuhn length). More details about these models can be found in Appendix A.

### 3. MODELS AND METHODS

We build the polymer networks by employing the method reported in Gnan et al.,²⁷ which makes use of the self-assembly of a binary mixture of limited-valence particles. Particles of species A can form up to four bonds (valence $\phi = 4$) and bond only to B particles, thus acting as cross-linkers. Particles of species B can form up to two bonds ($\phi = 2$) and can bond to A and B particles. We carry out the assembly of $N_{init} = N_A + N_B = 5 \times 10^4$ particles at different number densities $\rho_{init} = N_{init}/V$, with $V$ the volume of the simulation box, and different cross-linker concentrations $C = N_A/(N_B + N_A)$. We consider two initial densities $\rho_{init} = 0.1, 0.85$, and $C = 1, 5$, and 10%. The results are averaged over two system realizations for each pair of $\rho_{init}$ C values.

The assembly proceeds until an almost fully bonded percolating network is attained, that is, the fraction of formed bonds is at least $N_{bond}/N_{bond}^{max} = 99.9\%$, where $N_{bond}^{max} = (4N_A + 2N_B)/2$ is the maximum number of bonds. The self-assembly process is greatly accelerated thanks to an efficient bond-swapping mechanism.⁴² When the desired fraction $N_{bond}/N_{bond}^{max}$ is reached, we stop the assembly, identify the percolating network, and remove all particles or clusters that do not belong to it. Because some particles are removed, at the end of the
procedure the values of $\rho_{init}$ and C change slightly. However, these changes are small (at most 10%) and in the following we will hence use the nominal (initial) values of $\rho_{init}$ and C to refer to the different networks.

The normalized distribution of the chemical lengths n of the chains, $P(n)/P(1)$, which constitute the network is shown in Figure 2. Here,

<figure><img src="image_10.png" /><figcaption>Figure 2. Rescaled distribution of chain lengths for all the simulated systems. We report the data for two samples (S1 and S2) generated with two values of initial density each. The orange dashed lines are the theoretical prediction of eq 10.</figcaption></figure>

the chemical chain length is defined as the number of particles in a chain, excluding the cross-linkers, so that a chain with n + 1 bonds has length n.⁶ In all cases, the distribution decays exponentially, as it is also the case for random-cross-linking from a melt of precursor chains.⁷ This exponential behavior can be explained by the Flory–Huggins polymerization theory⁴³ or, equivalently, by the Wertheim theory for associating fluids applied to patchy particles,⁴⁴ and it is ultimately due to the independence of the bonding events. We also note that $P(n)$ does not depend on the initial density²⁷,⁴⁵ and, as one expects given the equilibrium nature of the assembly protocol, it is fully reproducible. This distribution can be estimated from the nominal values of $\phi$ and C via the well-known formula of Flory⁴³

$$ \frac{P(n)}{P(1)} = \left(1 - \frac{1}{\langle n \rangle}\right)^{n-1} \qquad (10) $$

where $\langle n \rangle = 2(1 - C)/\phi C$ is the mean chain length,⁴⁵ which, using the nominal cross-linker valence ($\phi = 4$) and concentrations, takes the values 49.5, 9.5, and 4.5 for C = 1, 5, and 10%, respectively. The

parameter-free theoretical probability distribution is shown as orange-dashed lines in Figure 2 and reproduces almost perfectly the numerical data.

The network contains a few defects in the form of dangling ends (chains that are connected to the percolating network by one cross-linker only and first-order loops, that is, chains having both ends connected to the same cross-linker.³² Because there are no excluded volume interactions, these defects are elastically inactive and, therefore, do not influence the elastic properties of the network.³²,⁴⁶,⁴⁷ For the configurations assembled at C = 1%, the percentage of particles belonging to the dangling ends is ≈10% for $\rho_{init} = 0.1$ and ≈6% for $\rho_{init} = 0.85$. For higher values of C, the percentages are much smaller (e.g., ≈2% for C = 5%, $\rho_{init} = 0.1$ and ≈1% for C = 10%, $\rho_{init} = 0.1$). In order to obtain an ideal fully bonded network, the dangling ends are removed. We note that during this procedure, the cross-linkers connected to dangling ends have their valence reduced from $\phi = 4$ to $\phi = 3$ or 2 (in the latter case, they become type B particles). The percentage of the so-created three-valent cross-linkers remains small: for $\rho_{init} = 0.1$, it is ≈15, 4, and 2% for C = 1, 5, and 10%, respectively. The presence of these cross-linkers slightly changes the average cross-linker valence, but does not influence the main results of this work.

Once the network is formed, we change the interaction potential, making the bonds permanent and thus fixing the topology of the network. Since we are interested in understanding the roles that topology and chain size distribution of a polymer network play in determining its elasticity, we consider interactions only between bonded neighbors, similar to what has been done in Duering et al.¹⁰ Particles that do not share a bond do not feel any mutual interaction, and hence chains can freely cross each other (whence the name phantom network). Two bonded particles interact through the widely used Kremer–Grest potential,⁴⁸ which is given by the sum of a repulsive term accounting for steric interactions, modeled via the Weeks–Chandler–Andersen potential,⁴⁹ and of a finite extensible nonlinear elastic (FENE) potential modeling the intramolecular bonds. The former is given by

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

where $\sigma$ is the monomer diameter and $\epsilon$ sets the scale of the repulsion, while the latter is

<figure><img src="image_11.png" /><figcaption>Figure 3. Snapshots of a network with cross-linker concentration C = 5% and assembly density $\rho_{init} = 0.1$ simulated at (a, b, e, and f) $\rho = 0.1$ and (c, d, g, and h) $\rho = 1.5$. We show configurations in equilibrium (a, c, e, and g) and subject to a uniaxial deformation along the vertical direction with $\lambda = 1.2$ (b, d, f, and h). Top row (panels a–d): turquoise and red particles indicate cross-linkers and monomers, respectively. The orange scale bars in the top left corners are 10σ long. Bottom row (panels e–h): the same configurations are represented in a way such that the chains are colored according to the ratio between their end-to-end distance and contour length (see the legend on the right).</figcaption></figure>
$$ \mathcal{U}_{\text{FENE}}(r) = -\frac{kr_0^2}{2}\ln\left[1 - \left(\frac{r}{r_0}\right)^2\right] \quad (12) $$

where $k = 30\epsilon/\sigma^2$ and $r_0 = 1.5\sigma$. Here and 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 11 and $m$ is the mass of a particle, which is the same for A and B particles. The units of temperature, density, time, and elastic moduli are, respectively, $[T] = \epsilon/k_B$, $[\rho] = \sigma^{-3}$, $[t] = \sqrt{m\sigma^2/\epsilon}$, and $[G] = \epsilon\sigma^{-3}$. In these units, the Kuhn length of the model is $b = 0.97$.⁴⁸

We run molecular dynamics simulations in the NVT ensemble at constant temperature $T = 1.0$ by employing a Nosé-Hoover thermostat.⁵⁰ Simulations are carried out using the LAMMPS simulation package,⁵¹ with a simulation time step $\delta t = 0.003$.

In order to study the effects of the density on the elastic properties, the initial configurations are slowly and isotropically compressed or expanded to reach the target densities $\rho = 0.1, 0.2, 0.5, 0.85$, and $1.5$. Then, a short annealing of $10^6$ steps and subsequently a production run of $10^7$ steps are carried out. Even for the system with the longest chains, the mean-squared displacement of the single particles reaches a plateau, indicating that the chains have equilibrated (see the Supporting Information).

For each final density value, we run several simulations for which we perform a uniaxial deformation in the range $\lambda_\alpha \in [0.8, 1.2]$ along a direction $\alpha$, where $\lambda_\alpha = L_\alpha/L_{\alpha,0}$ is the extent of the deformation and $L_{\alpha,0}$ and $L_\alpha$ are the initial and final box lengths along $\alpha$, respectively. The deformation is carried out at constant volume with a deformation rate of $10^{-1}$. To confirm that the system is isotropic, we perform the deformation along different spatial directions $\alpha$.

Figure 3 shows representative snapshots of the $C = 5\%$, $\rho_{\text{init}} = 0.1$ system at low ($\rho = 0.1$) and high ($\rho = 1.5$) density, in equilibrium and subject to a uniaxial deformation along the vertical direction. In panels 3a–d, we show the particles (monomers and cross-linkers), highlighting the highly disordered nature of the systems and their structural heterogeneity, which is especially evident at low density. The same systems are also shown in panels 3e–h, where we use a gradient to color chains according to the ratio between their end-to-end distance and contour length. These panels highlight the effect that density has on the heterogeneous elastic response of these systems when they are subject to deformations.

Once the system acquires the target value of $\lambda_z$, we determine the diagonal elements of the stress tensor $\sigma_{aa}$ and compute the engineering stress $\sigma_{\text{eng}}$ as²²

$$ \sigma_{\text{eng}} = \frac{\sigma_{\text{tr}}}{\lambda_z} = \frac{1}{\lambda_z} \left[ \sigma_{zz} - \frac{1}{2}(\sigma_{xx} + \sigma_{yy}) \right] \quad (13) $$

where $\sigma_{\text{tr}}$ is the so-called true stress.¹,⁵² The shear modulus $G$ is then the quantity that connects the engineering stress and the strain through the following relation²²

$$ \sigma_{\text{eng}} = G \left[ (\lambda_z - \lambda_{\text{ref}}) - \frac{1}{(\lambda_z - \lambda_{\text{ref}})^2} \right] \quad (14) $$

In eq 14, $\lambda_{\text{ref}}$ is an extra fit parameter that we add to take into account the fact that in some cases $\sigma_{\text{eng}} \neq 0$ for $\lambda_z = 1$, which signals the presence of some prestrain in our configurations. The stress–strain curves we use to estimate $G$ are averaged over 10 independent configurations obtained by the randomization of the particle velocities, prior to deformation, with a Gaussian distribution of mean value $T = 1.0\epsilon/k_B$ in order to reduce the statistical noise.

Figure 4 shows the numerical data for the stress–strain curves for the $C = 1\%$, $\rho_{\text{init}} = 0.1$ system. We also report the associated theoretical curves, fitted to eq 14, through which we obtain an estimate of the shear modulus.

<figure><img src="image_12.png" /><figcaption>Figure 4. Example of stress–strain curves for the $C = 1\%$, $\rho_{\text{init}} = 0.1$ system. Symbols are simulation data, lines are fits with eq 14.</figcaption></figure>

## 4. RESULTS AND DISCUSSION

We use the simulation data to estimate $\tilde{r} \equiv (\overline{r^2})^{1/2}$ (RMS end-to-end distance) and $R \equiv \overline{r}$ for each chain to compute the elastic moduli of the networks through eq 5. In the following, we will refer to the elastic moduli computed in this way with the term “theoretical”.

Figure 5a shows the shear modulus as computed in simulations for all investigated systems as a function of $\nu$, the density of elastically active strands. We note that, to a first approximation, one would have $\nu = \phi C\rho/2$. However, the actual value of $\nu$ is slightly smaller than this figure because of the presence of elastically inactive strands, as discussed in Section 3. First of all, we observe that systems generated at the same $C$ but with different values of $\rho_{\text{init}}$ exhibit markedly different values of the shear modulus when probed under the same conditions (i.e., the same strand density). This result highlights the fundamental role of the cross-linking process, which greatly affects the initial distribution of the chains’ end-to-end distances even when the number of chains and their chemical length distribution, being dependent only on $C$ (see Figure 2), are left unaltered. Thus, the echo of the difference between the initial end-to-end distributions gives rise to distinct elastic properties of the phantom networks even at the same strand density.

In Figure 5a, we also plot the behavior predicted by eq 9 (blue line), which assumes Gaussian-distributed end-to-end distances. Even though the numerical data seem to approach this limit at very large values of the density, they do so with a slope that is clearly smaller than unity. For the $C = 1\%$ sample, this slope is almost exactly 1/3, and it is also very close to this value for the $C = 5\%$ and $C = 10\%$ samples assembled at $\rho_{\text{init}} = 0.1$. This behavior can be understood at the qualitative level from eq 8: $R$ is the average distance between cross-linkers and, therefore, it changes affinely upon compression or expansion, thereby scaling as $R \propto \nu^{-1/3}$.⁵³,⁵⁴ As a result, in the Gaussian limit, the shear modulus scales as $G^G \propto \nu^{1/3}$.²¹,⁵³⁻⁵⁵ As discussed above, our results show that the way this limiting regime is approached depends on the cross-linker concentration $C$ and on the preparation state, which is here controlled by $\rho_{\text{init}}$.

The quantitative differences in the elastic response of systems with different $C$ and $\rho_{\text{init}}$ can be partially rationalized by looking at the scaling properties of the end-to-end distances. We notice that the RMS equilibrium end-to-end distance $R(n)$ of the strands for different values of $\rho_{\text{init}}$ and $C$ nearly collapses on a master curve when divided by the initial cross-linker
<figure><img src="image_13.png" /><figcaption>**Figure 5.** (a) Shear modulus as a function of the elastically active strand number density $\nu$ for all the investigated systems. Solid blue line: eq 9 with $\phi = 4$. Solid orange line: slope 1/3. (b) Same as (a), with both G and $\nu$ rescaled by $\gamma\rho_{init}^A$ where $\gamma = 0.74$ for $\rho_{init} = 0.85$, $\gamma = 1$ for $\rho_{init} = 0.1$ is a fit parameter, and $\rho_{init}^A = C\rho_{init}$ is the initial cross-linker density.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 6.** Comparison between the shear moduli obtained through eq 5 and three different approximations and the numerical ones (G) for the simulated systems (see legends). Dashed-dotted line/stars in panel b: FJC approximation with no overstretched chains (chains with $r \ge 0.95 \cdot nb$).</figcaption></figure>

density, $\rho_{init}^A = C\rho_{init}$ (see Supporting Information). A slightly better agreement is found if the heuristic factor $\gamma\rho_{init}^A$ with $\gamma = 0.74$ for $\rho_{init} = 0.85$ and $\gamma = 1$ for $\rho_{init} = 0.1$, is used. Based on this observation, we rescale the data of Figure 5a multiplying both G and $\nu$ by $\gamma\rho_{init}^A$. The result is shown in Figure 5b: one can see that the shear modulus of systems with the same C but different values of $\rho_{init}$ nicely fall on the same curve. Moreover, in the large-$\nu$ limit, where all the curves tend to have the same slope, a good collapse of the data of systems with different C is also observed.

The differences arising between systems at different C can be explained by noting that the cross-linker concentration controls the relative abundance of chains with different $n$, whose elastic response cannot be rescaled on top of each other by using $n$ but depends on their specific end-to-end distribution (see e.g., Appendix A). As a result, the elasticity of networks generated at different C cannot be rescaled on top

of each other. In particular, systems with more cross-linkers, and hence more short chains, will deviate earlier and more strongly from the Gaussian behavior.

Interestingly, G exhibits a nonmonotonic behavior as a function of $\nu$; this feature appears for all but the lowest C and $\rho_{init}$ values. This behavior, which has also been observed in hydrogels,$^{21,30,53,54,56}$ cannot be explained assuming that the chains are Gaussian because in this case one has for all $\nu$ that G $\propto \nu^{1/3}$, as discussed above. Given that our model features stretchable bonds, at large strains it cannot be considered to be a FJC, being more akin to an ex-FJC.$^{40}$ Therefore, one might be tempted to ascribe the increase of G upon decreasing $\nu$ to the energetic contribution. For this reason, in addition to the Gaussian and FJC descriptions, we also plot in Figure 6b the shear modulus estimated by neglecting the contributions of those chains that have $r \ge 0.95 \cdot nb$. We note that chains with $r < 0.95 \cdot nb$ behave essentially as FJCs, and that in any case the
number of such chains is minuscule in all but the lowest-density $C = 10\%$ systems (see the Supporting Information). Because the sets of data with and without the overstretched chains overlap almost perfectly, we confirm that the energetic contribution due to the few overstretched chains is negligible: we can thus conclude that the nonmonotonicity we observe has a purely entropic origin. This holds true for all the systems investigated except for the $C = 10\%$, $\rho_{init} = 0.85$ system, which contains the largest number of short, overstretched chains (see Supporting Information).

In Figure 6, we compare the numerical shear modulus for all investigated systems with estimates as predicted by different theories, with the common assumption that the three-chain model remains valid (see Section 2). In particular, we show the results obtained with the FJC (eq 6), Gaussian (eq 7), and ex-FJC (see Appendix A) models. One can see that the agreement between the theoretical and numerical results is always better for larger values of $\nu$, that is, when chains are less stretched. Moreover, the agreement between data and theory is better for systems generated at smaller $\rho_{init}$. We note that the Gaussian approximation, which predicts a monotonically increasing dependence on $\nu$, fails to reproduce the qualitative behavior of G, whereas the ex-FJC systematically overestimates G. The FJC description is the one that consistently achieves the best results, although it fails (dramatically at large $C$) at small densities. We ascribe this qualitative behavior to the progressive failure of the three-chain assumption as the density decreases. Because the three-chain model is known to overestimate the stress at large strains compared to more complex and realistic approximations such as the tetrahedral model,¹ the resulting single-chain contribution to the elastic modulus for stretched chains is most likely overestimated as well. Regardless of the specific model used, our results suggest that when the samples are strongly swollen, something that can be achieved in experiments,²¹ any description that attempts to model the network as a set of independent chains gives rise to an unreliable estimate of the overall elasticity even when energetic contributions due to stretched bonds do not play a role.

In addition to providing the best comparison with the numerical data in the whole density range, the FJC description also captures the presence and (although only in a semi-quantitative fashion) the position of the minimum. This is the case for all the investigated systems, highlighting the role played by the short chains, whose strong non-Gaussian character heavily influences the overall elasticity of the network.

Although real short chains do not follow the exact end-to-end probability distribution we use here (see eq 6), they are surely far from the scaling regime and hence they should never be regarded as Gaussian chains, even in the melt or close to the theta point. This aspect has important consequences for the analysis of experimental randomly cross-linked polymer networks, for which one may attempt to extract some microscopic parameter (such as the contour length or the average end-to-end distance) by fitting the measured elastic properties to some theoretical relations such as the ones we discuss here. Unfortunately, such an approach will most likely yield unreliable estimates. We demonstrate that this happens even with a very idealized system such as the one we consider here. In order to follow the procedure usually applied to these polymeric systems,²¹ we make the assumption that the network can be considered as composed of $N_s$ strands of $\langle n \rangle$ segments.

We then compare in Figure 7a-c, the numerically estimated values of $G$ with those obtained with the L-FJC model (see

<figure><img src="image_15.png" /><figcaption>**Figure 7.** Fitting results to (a-c) simulation data (a) $C = 1\%$, (b) $C = 5\%$, and (c) $C = 10\%$ and (d,e) experimental data: (d) Young's modulus taken from Hoshino et al.²¹ and (e) shear modulus taken from Matsuda et al.³⁰ The quality of the fits in panel e does not depend on whether the point at $Q = 1$ is considered or not or if we restrict the fit to swelling ratios $Q \lesssim 15$.</figcaption></figure>

Appendix A). The expression we employ contains two quantities that can be either fixed or fitted to the data: the average end-to-end distance in a specific state (e.g., the preparation state), $R_0$, and the average strand length $\langle n \rangle$ (or, equivalently, the contour length $r_{max} = \langle n \rangle b$). Together with the numerical data, in Figure 7 we present three sets of theoretical curves: $G$ as estimated by using the simulation values of $R_0$ and $\langle n \rangle$ or fitted by using either $R_0$ or both quantities as free parameters. If $C$ is small (and hence $\langle n \rangle$ is large), the difference between the parameter-free expression and the numerical data is small (10-15%). However, as $C$ becomes comparable with the values that are often used in real randomly cross-linked hydrogels (≈5%), the difference between the theoretical and simulation data becomes very significant: for instance, for $C = 10\%$ the parameter-free expression fails to even capture the presence of the minimum. Fitting the numerical data makes it possible to achieve an excellent agreement, although the values of the parameters come out to be sensibly different (sometimes more than 50%) from the real values (see Appendix A). Our results thus show that even in the simplest randomly cross-linked system—a phantom network of freely jointed chains—neglecting the shortness of the majority of the chains, which dominate the elastic response, can lead to a dramatic loss of accuracy. Randomly cross-linked polymer networks contain short chains, which are inevitably quite far from the scaling regime, and hence even their qualitative behavior can become elusive if
looked through the lens of polymer theories that rest too heavily on the Gaussianity of the chains.

To conclude our analysis, we also apply our theoretical expressions to two sets of data, which have been recently published. Both experiments have been carried out in the group of Gong.²¹,³⁰ The first system is a tetra-PEG hydrogel composed of monodisperse long chains that can be greatly swollen by using a combined approach of adding a molecular stent and applying a PEG dehydration method.²¹ Because the system is monodisperse and the chains are quite long, we expect the theoretical expressions derived here to work well. Indeed, as shown in Figure 7d, the resulting Young's modulus is a nonmonotonic function of the swelling ratio Q, defined as the ratio between the volume at which the measurements are performed and the volume at which the sample was synthesized. The experimental data can be fitted with both the L-FJC and WLC expressions (see Appendix A) because both models reproduce the data with high accuracy when fitted with the two free parameters introduced above (R₀ and ⟨n⟩). However, better results are obtained with the WLC model, which fits well when ⟨n⟩ is fixed to its experimentally estimated value, yielding a value R₀ = 7.2 nm, which is very close to the independently estimated value of 8.1 nm,²¹ in agreement with what reported in Hoshino et al.²¹

The second system we compare to is a randomly cross-linked PNaAMPS network, for which the shear modulus as a function of Q has been reported.³⁰ As shown in Figure 7e, the theoretical expressions reported here cannot go beyond a qualitative agreement with the experimental data, even if two parameters are left free and we only fit to the experimental data in a narrow range of swelling ratios. In addition, the fitting procedure always yields unphysical values for the two parameters (e.g., R₀ comes out to be smaller than 1 nm, see Appendix A). Although part of the discrepancy might be due to the charged nature of the polymers involved,⁵⁷ we believe that the disagreement between the theoretical and experimental behaviors can be partially ascribed to the randomly cross-linked nature of the network, and hence to the abundance of short chains. Because the end-to-end distribution of such short chains is not known and depends on the chemical and physical details, there is no way of taking into account their contribution to the overall elasticity in a realistic way. These results thus highlight the difficulty of deriving a theoretical expression to assess the elastic behavior of randomly cross-linked real networks.

## 5. SUMMARY AND CONCLUSIONS

We have used numerical simulations of disordered phantom polymer networks to understand the role played by the chain size distribution in determining their elastic properties. In order to do so, we employed an in silico synthesis technique by means of which we can independently control the number and chemical size of the chains, set by the cross-linker concentration, as well as the distribution of their end-to-end distances, which can be controlled by varying the initial monomer concentration. We found that networks composed of chains of equal contour length can have shear moduli that depend strongly on the end-to-end distance even when probed at the same strand concentration. Hence, this shows that even in simple systems the synthesis protocol can have a large impact on the final material properties of the network even when it does not affect the chemical properties of its basic constituents, as recently highlighted in a microgel system.⁵⁸

We then compared the results from the simulations of the phantom network polymer theory, which was revisited to obtain explicit expressions for the shear modulus assuming three different chain conformation fluctuations, namely, the exact freely jointed chain, Gaussian, and extensible freely jointed chain models. We observed a nonmonotonic behavior of G as a function of the strand density that, thanks to a comparison with the theoretical results, can be completely ascribed to entropic effects that cannot be accounted for within a Gaussian description. We thus conclude that the role played by short-stretched chains in the mechanical description of polymer networks is fundamental and should not be overlooked. This insight is supported by an analysis of experimental data of the elastic moduli of hydrogels reported in the literature. We are confident that the numerical and analytical tools employed here can be used to address similar and other open questions concerning both the dynamics and the topology in systems in which excluded-volume effects are also taken into account, and hence entanglements effects may be relevant. Investigations in this direction are underway.

# APPENDIX A

## Shear Modulus of a System with Gaussian-Distributed End-to-End Distances

In this section, we show how eq 9 can be derived for a polydisperse network. Similar derivations for the case of monodisperse networks can be found in standard textbooks.¹⁹,²² We start by noting something that is sometimes overlooked: the Gaussian distribution $W_n^G(r)$, eq 7, applies to a single chain. However, at the ensemble level, the distribution of end-to-end vectors, which we may call $\Omega[r(n)]$, is no Gaussian in general. However, if we assume, for example, that the system has been obtained through end-cross-linking starting from a melt of precursor chains,¹⁰ then $\Omega[r(n)] = W_n^G(r)$,³⁴ so that the magnitudes r of the r vectors will be Gaussian-distributed. Under this assumption, one has

$$ \left\langle \frac{r^2}{nb^2} \right\rangle = \left\langle \frac{R^2}{nb^2} \right\rangle + \left\langle \frac{u^2}{nb^2} \right\rangle = 1 \qquad (15) $$

To evaluate the term in brackets in eq 8, we thus only need to evaluate the fluctuation term $\left\langle \frac{u^2}{nb^2} \right\rangle$. This term can be computed using the equipartition theorem (the same derivation can be found, e.g., in ref 59). The total energy of the fluctuations is $U_{\text{fluct}} = \frac{3}{2}k_B T N_x$, where $N_x$ is the number of active cross-linkers because there is one mode for each node and to each mode it is associated an energy of $3/2k_B T$.⁵⁹ Moreover, $N_x = 2N_s/\phi$, with $N_s$ the number of elastically active strands. Therefore, the mean energy per strand is

$$ \frac{U_{\text{fluct}}}{N_s} = \frac{3k_B T}{\phi} = \frac{3k_B T}{2N_s} \sum_i^{N_s} \frac{\overline{u_i^2}}{n_i b^2} = \frac{3k_B T}{2} \left\langle \frac{\overline{u^2}}{nb^2} \right\rangle \qquad (16) $$

where the sum extends over all the elastically active strands. Therefore, we get $\left\langle \frac{u^2}{nb^2} \right\rangle = \frac{2}{\phi}$ (a generalization to the polydisperse case of a well-known result for the phantom network²⁰,⁶⁰), from which we finally obtain

$$ \left\langle \frac{R^2}{nb^2} \right\rangle = 1 - \frac{2}{\phi} \qquad (17) $$
From eq 8, we obtain eq 9, that is, $G^G = (1 - \frac{2}{\phi})k_B T \nu$.

The validity of eq 17 depends not only on the cross-linking procedure but also on the macroscopic thermodynamic parameters (such as solvent quality, density, or pressure). For example, if the chain-size distribution is such that short chains are abundant, as it is the case for randomly cross-linked networks,⁷,²⁰ the front factor A (see eq 8) will depend on the chain size distribution because short chains are non-Gaussian. Another example is when the cross-linking procedure is performed in a state in which the chains are non-Gaussian, for example, under good solvent conditions, where the chains behave as self-avoiding random walks.²²

## Models of Chain Statistics

The freely jointed chain approximation, eq 6, describes an inextensible chain because each component $f_α$ where $α = x, y$, and $z$, of the force $f$ required to stretch the chain, for example

$$f_x = -T \frac{dS_n(\mathbf{r})}{dx} = - \frac{k_B T}{W_n(\mathbf{r})} \frac{dW_n(\mathbf{r})}{dx} \quad (18)$$

diverges in the limit $r \to nb$, that is, when the end-to-end distance approaches the contour length. In the limit of large strains and large degree of polymerization $n$, a better approximation is provided by the well-known Langevin dependence of the elongation on the exerted force $f$, which yields for the end-to-end probability distribution function¹

$$W_n^L(\mathbf{r}) = A \exp\left[-\frac{r}{b}\mathcal{L}^{-1}(r/nb)\right] \left[\frac{\mathcal{L}^{-1}(r/nb)}{\sinh \mathcal{L}^{-1}(r/nb)}\right]^{-n} \quad (19)$$

where $β = 1/k_B T$, $T$ is the temperature, $k_B$ is the Boltzmann constant, $A$ is a normalization constant, and $\mathcal{L}(x)^{-1}$ is the inverse Langevin function. The latter is defined as $\mathcal{L}(x) = \coth(x) - 1/x$ and it turns out to be equal to the ratio between the end-to-end distance and the contour length of a chain that is subject to a force $f$

$$\frac{r}{nb} = \coth\left(\frac{fb}{k_B T}\right) - \frac{k_B T}{fb} = \mathcal{L}(\beta fb) \quad (20)$$

We note on passing that $\mathcal{L}^{-1}(\cdot)$ cannot be written in a closed form⁶¹ and hence must be evaluated numerically.

Phantom Kremer-Grest chains behave exactly as freely jointed chains up to end-to-end distances that are very close to the contour length. However, beyond those values the Langevin description is no longer valid, and one has to resort to the ex-FJC model, for which the following analytical form of the force-extension curve has been recently derived⁴⁰

$$r = nb \mathcal{L}(\beta fb) + \frac{nf}{k} \left[ 1 + \frac{1 - \mathcal{L}(\beta fb) \coth(\beta fb)}{1 + \frac{f}{kb} \coth(\beta fb)} \right] \quad (21)$$

where $k$ is the monomer-monomer force constant in the harmonic approximation. In principle, it is possible to integrate the inverse of this relation to get the end-to-end probability density. However, as it is clear from eq 5, we only need the derivatives of the inverse function, and hence there is no need to obtain $W(r)$ explicitly. We set the value of $k$ to the value of the second derivative of the Kremer-Grest potential as computed in the minimum, $k \approx 867 \epsilon/\sigma^2$. We have also fitted the force-extension curves as obtained in simulations of single

chains, obtaining values that are compatible with this latter estimate.

We also report the following expression, which approximates the force–elongation relation for a WLC⁴² and is used in the main text to fit the experimental data shown in Figure 7c

$$\frac{fb}{k_B T} = \frac{1}{4} \left(1 - \frac{r}{nb}\right)^{-2} - \frac{1}{4} + \frac{r}{nb} - 0.8 \left(\frac{r}{nb}\right)^{2.15} \quad (22)$$

Figure 8 shows the force–extension curve for polymer chains described with different models. Because the force is

<figure><img src="image_16.png" /><figcaption>Figure 8. Force required to extend a chain by $r$ (scaled by its contour length $nb$) for the Gaussian (dashed line, eq 7), L-FJC and ex-JFC (solid lines, eq 21), and FJC (dashed-dotted lines, eq 6) models at $T = 1$. Note that being plotted in this way, $f$ depends on $n$ only for the latter case.</figcaption></figure>

plotted as a function of the distance scaled by the FJC contour length $nb$, the Gaussian, L-FJC and ex-FJC descriptions become independent of $n$. By contrast, the extensional force of the exact FJC model, whose end-to-end probability distribution is given by eq 6, retains an $n$-dependence that is very strong for small $n$ and decreases upon increasing the chain size. Indeed, for $n \gtrsim 40$, the resulting force is essentially $n$-independent and overlaps almost completely with the FJC curve.

## Fitting Procedure and Additional Results

As discussed in the main text in Section 4, we fit the experimental Young’s and shear moduli reported in Hoshino et al.²¹ (system A) and Matsuda et al.³⁰ (system B). As commonly done in the analysis of experimental systems, in both cases, we will consider the networks to be formed by strands of average size $\langle n \rangle$ and use relations based on eq 5: for system A, we use

$$G_{\text{exp}} = -\frac{TR^2\nu}{3} \left[ \frac{d s_n(\tilde{r})}{d\tilde{r}} \left( \frac{1}{\tilde{r}} - \frac{R^2}{2\tilde{r}^3} \right) + \frac{d^2 s_n(\tilde{r})}{d\tilde{r}^2} \frac{R^2}{2\tilde{r}^2} \right] \quad (23)$$

while for system B, we fit $Y_{\text{exp}} = 3G_{\text{exp}}$. We fit the experimental data with the L-FJC and WLC models by using eqs 19 and 22 to numerically evaluate the derivatives of $s_n(\tilde{r})$.

We evaluate $R$, $\tilde{r}$, and $\nu$ as follows. Let $V_0$, $\nu_0$, and $R_0$ be the volume, chain number density, and average end-to-end distance of the polymer chains in the preparation state, respectively, and $V$, $\nu$, and $R$ be the same quantities for the generic state point at which the experimental measurements are carried out. We define the swelling ratio $Q = V/V_0 = \nu_0/\nu$

so that $\nu = \nu_0/Q$ and $R = \left(\frac{V}{V_0}\right)^{1/3} R_0 = Q^{1/3}R_0$. Since
$\bar{r}^2 = \bar{r}^2 = R^2 + \bar{u}^2$ and, as shown in Appendix A, $\bar{u}^2 = \frac{2}{\phi}nb^2$,
for systems with $\phi = 4$ we also have that $\bar{r} = \sqrt{R^2 + r_{\max}b/2}$,
where $r_{\max} = \langle n \rangle b$ is the contour length of the strands. Using
the numbers reported in the original papers, we found that $\nu_0^A = 1.97 \times 10^{-3} \text{ nm}^{-3}$ and $\nu_0^B = 0.245 \text{ nm}^{-3}$. For system A, the
authors also report independent estimates for $R_0$ ($R_0^A = 8.1$ nm) and $r_{\max}$ ($r_{\max}^A = 82$ nm). For system A, we show the fitting
results obtained by either fixing $r_{\max}$ and using $R_0$ as a fitting
parameter or by leaving both quantities as fitting parameters.
For system B, we do the latter.

For system A, we obtain $R_0^{\text{FJC}} = 7.46$ nm and $R_0^{\text{WLC}} = 7.2$ nm for the one-parameter fits and $R_0^{\text{FJC}} = 5.08$ nm, $r_{\max}^{\text{FJC}} = 42.75$ nm and $R_0^{\text{WLC}} = 7.2$, $r_{\max}^{\text{WLC}} = 63.87$ nm for the two-parameter fits.

For system B, we obtain unphysical values ($R_0^{\text{FJC}} = 0.11$ nm, $r_{\max}^{\text{FJC}} = 0.65$ nm and $R_0^{\text{WLC}} = 0.17$, $r_{\max}^{\text{WLC}} = 0.96$ nm). The results do not improve if we restrict the fitting range to narrower Q-ranges.

## ■ ASSOCIATED CONTENT

### ● Supporting Information

The Supporting Information is available free of charge at
https://pubs.acs.org/doi/10.1021/acs.macromol.1c00176.

Shear modulus in the affine network model, monomer mean-squared displacement during equilibration, and density scaling of the RMS equilibrium end-to-end distance (PDF)

## ■ AUTHOR INFORMATION

### Corresponding Authors

**Valerio Sorichetti** – Laboratoire de Physique Théorique et Modèles Statistiques (LPTMS), CNRS, Université Paris-Saclay, F-91405 Orsay, France; Laboratoire Charles Coulomb (L2C), University of Montpellier, CNRS, F-34095 Montpellier, France; IATE, University of Montpellier, INRAE, Institut Agro, F-34060 Montpellier, France;
Email: valerio.sorichetti@universite-paris-saclay.fr

**Andrea Ninarello** – CNR-ISC Uos Sapienza, IT-00185 Roma, Italy; Department of Physics, Sapienza Università di Roma, IT-00185 Roma, Italy; ● orcid.org/0000-0001-5817-3084; Email: andrea.ninarello@roma1.infn.it

**Lorenzo Rovigatti** – Department of Physics, Sapienza Università di Roma, IT-00185 Roma, Italy; CNR-ISC Uos Sapienza, IT-00185 Roma, Italy; ● orcid.org/0000-0001-5017-2829; Email: lorenzo.rovigatti@uniroma1.it

### Authors

**José M. Ruiz-Franco** – CNR-ISC Uos Sapienza, IT-00185 Roma, Italy; Department of Physics, Sapienza Università di Roma, IT-00185 Roma, Italy

**Virginie Hugouvieux** – IATE, University of Montpellier, INRAE, Institut Agro, F-34060 Montpellier, France;
● orcid.org/0000-0001-6340-6487

**Walter Kob** – Laboratoire Charles Coulomb (L2C), University of Montpellier, CNRS, F-34095 Montpellier, France; Institut Universitaire de France, 75005 Paris, France

**Emanuela Zaccarelli** – CNR-ISC Uos Sapienza, IT-00185 Roma, Italy; Department of Physics, Sapienza Università di Roma, IT-00185 Roma, Italy

Complete contact information is available at:
https://pubs.acs.org/10.1021/acs.macromol.1c00176

## Author Contributions

†V.S. and A.N. contributed equally.

## Notes

The authors declare no competing financial interest.

## ■ ACKNOWLEDGMENTS

We thank F. Sciortino and D. Truzzolillo for helpful discussions and J. P. Gong, T. Nakajima, and T. Matsuda for sharing the data reported in Figure 7d,e. We acknowledge financial support from the European Research Council (ERC Consolidator Grant 681597, MIMIC). W.K. is a senior member of the Institut universitaire de France. This project was supported by the LabEx NUMEX (ANR-10-LABX-0020) within the I-Site MUSE.

## ■ ADDITIONAL NOTES

<sup>a</sup>In the absence of excluded volume interactions, the pressure of the system is negative and would, therefore, collapse if the volume was not kept constant.

<sup>b</sup>Here and in the following, italic is used to indicate the magnitude of the vector.

<sup>c</sup>We recall that because two cross-linkers cannot bind to each other, the minimum chain length is $n = 1$ (for two bonds)

## ■ REFERENCES

(1) Treloar, L. R. G. *The Physics of Rubber Elasticity*; Oxford University Press: USA, 1975.

(2) Treloar, L. R. G. The elasticity of a network of long-chain molecules I. *Trans. Faraday Soc.* **1943**, *39*, 36–41.

(3) Flory, P. J.; Rehner, J., Jr. Statistical mechanics of cross-linked polymer networks I. Rubberlike elasticity. *J. Chem. Phys.* **1943**, *11*, 512–520.

(4) Flory, P. J.; Rehner, J. Statistical mechanics of cross-linked polymer networks II. Swelling. *Chem. Phys.* **1943**, *11*, 521–526.

(5) Flory, P. J. Molecular theory of rubber elasticity. *Polym. J.* **1985**, *17*, 1.

(6) Broedersz, C. P.; MacKintosh, F. C. Modeling semiflexible polymer networks. *Rev. Mod. Phys.* **2014**, *86*, 995.

(7) Grest, G. S.; Kremer, K. Statistical properties of random cross-linked rubbers. *Macromolecules* **1990**, *23*, 4994–5000.

(8) Duering, E. R.; Kremer, K.; Grest, G. S. Relaxation of randomly cross-linked polymer melts. *Phys. Rev. Lett.* **1991**, *67*, 3531.

(9) Grest, G. S.; Kremer, K.; Duering, E. R. Kinetics of end crosslinking in dense polymer melts. *Europhys. Lett.* **1992**, *19*, 195.

(10) Duering, E. R.; Kremer, K.; Grest, G. S. Structure and relaxation of end-linked polymer networks. *J. Chem. Phys.* **1994**, *101*, 8169–8192.

(11) Everaers, R.; Kremer, K. Topological interactions in model polymer networks. *Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top.* **1996**, *53*, R37.

(12) Kenkare, N. R.; Smith, S. W.; Hall, C. K.; Khan, S. A. Discontinuous molecular dynamics studies of end-linked polymer networks. *Macromolecules* **1998**, *31*, 5861–5879.

(13) Everaers, R. Entanglement effects in defect-free model polymer networks. *New J. Phys.* **1999**, *1*, 12.

(14) Kenkare, N. R.; Hall, C. K.; Khan, S. A. Theory and simulation of the swelling of polymer gels. *J. Chem. Phys.* **2000**, *113*, 404–418.

(15) Svaneborg, C.; Everaers, R.; Grest, G. S.; Curro, J. G. Connectivity and entanglement stress contributions in strained polymer networks. *Macromolecules* **2008**, *41*, 4920–4928.

(16) Gula, I. A.; Karimi-Varzaneh, H. A.; Svaneborg, C. Computational Study of Cross-Link and Entanglement Contributions to the Elastic Properties of Model PDMS Networks. *Macromolecules* **2020**, *53*, 6907–6927.

(17) Rubinstein, M.; Panyukov, S. Nonaffine deformation and elasticity of polymer networks. *Macromolecules* **1997**, *30*, 8036–8044.
(18) Rubinstein, M.; Panyukov, S. Elasticity of polymer networks. *Macromolecules* **2002**, *35*, 6670–6686.

(19) Mark, J. E. *Physical Properties of Polymers Handbook*; Springer, 2007; Vol. 1076.

(20) Higgs, P. G.; Ball, R. C. Polydisperse polymer networks: elasticity, orientational properties, and small angle neutron scattering. *J. Phys. 1988*, *49*, 1785–1811.

(21) Hoshino, K.-I.; Nakajima, T.; Matsuda, T.; Sakai, T.; Gong, J. P. Network elasticity of a model hydrogel as a function of swelling ratio: from shrinking to extreme swelling states. *Soft Matter* **2018**, *14*, 9693–9701.

(22) Rubinstein, M.; Colby, R. H. *Polymer Physics*; Oxford University Press: New York, 2003.

(23) Everaers, R.; Kremer, K. Test of the foundations of classical rubber elasticity. *Macromolecules* **1995**, *28*, 7291–7294.

(24) Escobedo, F. A.; de Pablo, J. J. Simulation and theory of the swelling of athermal gels. *J. Chem. Phys.* **1997**, *106*, 793–810.

(25) Escobedo, F. A.; De Pablo, J. J. Molecular simulation of polymeric networks and gels: phase behavior and swelling. *Phys. Rep.* **1999**, *318*, 85–112.

(26) Richbourg, N. R.; Peppas, N. A. The swollen polymer network hypothesis: Quantitative models of hydrogel swelling, stiffness, and solute transport. *Prog. Polym. Sci.* **2020**, *105*, 101243.

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

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

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

(30) Matsuda, T.; Nakajima, T.; Gong, J. P. Fabrication of tough and stretchable hybrid double-network elastomers using ionic dissociation of polyelectrolyte in nonaqueous media. *Chem. Mater.* **2019**, *31*, 3766–3776.

(31) Zhong, M.; Wang, R.; Kawamoto, K.; Olsen, B. D.; Johnson, J. A. Quantifying the impact of molecular defects on polymer network elasticity. *Science* **2016**, *353*, 1264–1268.

(32) Lin, T.-S.; Wang, R.; Johnson, J. A.; Olsen, B. D. Revisiting the elasticity theory for real Gaussian phantom networks. *Macromolecules* **2019**, *52*, 1685–1694.

(33) Landau, L. D.; Lifshitz, E. M. *Course of Theoretical Physics: Theory of Elasticity*, 2nd ed.; Pergamon Press, 1970; Vol. 7.

(34) Flory, P. J. Statistical thermodynamics of random networks. *Proc. Roy. Soc. Lond. Math. Phys. Sci.* **1976**, *351*, 351–380.

(35) Smith, T. L. Modulus of tightly crosslinked polymers related to concentration and length of chains. *J. Polym. Sci., Polym. Symp.* **1974**, *46*, 97–114.

(36) Jernigan, R. L.; Flory, P. J. Distribution functions for chain molecules. *J. Chem. Phys.* **1969**, *50*, 4185–4200.

(37) James, H. M.; Guth, E. Theory of the elastic properties of rubber. *J. Chem. Phys.* **1943**, *11*, 455–481.

(38) Tobolsky, A. V.; Carlson, D. W.; Indictor, N. Rubber elasticity and chain configuration. *J. Polym. Sci.* **1961**, *54*, 175–192.

(39) Toda, M.; Morita, H. Rubber elasticity of realizable ideal networks. *AIP Adv.* **2018**, *8*, 125005.

(40) Fiasconaro, A.; Falo, F. Analytical results of the extensible freely jointed chain model. *Phys. Stat. Mech. Appl.* **2019**, *532*, 121929.

(41) Petrosyan, R. Improved approximations for some polymer extension models. *Rheol. Acta* **2017**, *56*, 21–26.

(42) Sciortino, F. Three-body potential for simulating bond swaps in molecular dynamics. *Eur. Phys. J. E: Soft Matter Biol. Phys.* **2017**, *40*, 3.

(43) Flory, P. J. *Principles of Polymer Chemistry*; Cornell University Press, 1953.

(44) Sciortino, F.; Bianchi, E.; Douglas, J. F.; Tartaglia, P. Self-assembly of patchy particles into polymer chains: A parameter-free comparison between Wertheim theory and Monte Carlo simulation. *J. Chem. Phys.* **2007**, *126*, 194903.

(45) Rovigatti, L.; Gnan, N.; Zaccarelli, E. Internal structure and swelling behaviour of in silico microgel particles. *J. Phys.: Condens. Matter* **2017**, *30*, 044001.

(46) Lang, M. Elasticity of phantom model networks with cyclic defects. *ACS Macro Lett.* **2018**, *7*, 536–539.

(47) Panyukov, S. Loops in polymer networks. *Macromolecules* **2019**, *52*, 4145–4153.

(48) Kremer, K.; Grest, G. S. Dynamics of entangled linear polymer melts: A molecular-dynamics simulation. *J. Chem. Phys.* **1990**, *92*, 5057–5086.

(49) Weeks, J. D.; Chandler, D.; Andersen, H. C. Role of repulsive forces in determining the equilibrium structure of simple liquids. *J. Chem. Phys.* **1971**, *54*, 5237–5247.

(50) Martyna, G. J.; Klein, M. L.; Tuckerman, M. Nosé–Hoover chains: The canonical ensemble via continuous dynamics. *J. Chem. Phys.* **1992**, *97*, 2635–2643.

(51) Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. *J. Comput. Phys.* **1995**, *117*, 1–19.

(52) Doi, M. *Introduction to Polymer Physics*; Oxford university Press, 1996.

(53) Gundogan, N.; Melekaslan, D.; Okay, O. Rubber elasticity of poly (N-isopropylacrylamide) gels at various charge densities. *Macromolecules* **2002**, *35*, 5616–5622.

(54) Horkay, F.; Tasaki, I.; Basser, P. J. Osmotic swelling of polyacrylate hydrogels in physiological salt solutions. *Biomacromolecules* **2000**, *1*, 84–90.

(55) Panyukov, S. Scaling theory of high elasticity. *Sov. Phys. JETP* **1990**, *71*, 372–379.

(56) Itagaki, H.; Kurokawa, T.; Furukawa, H.; Nakajima, T.; Katsumoto, Y.; Gong, J. P. Water-induced brittle-ductile transition of double network hydrogels. *Macromolecules* **2010**, *43*, 9495–9500.

(57) Fisher, L. W.; Sochor, A. R.; Tan, J. S. Chain characteristics of poly (2-acrylamido-2-methylpropanesulfonate) polymers. *I. Light-scattering and intrinsic-viscosity studies. Macromolecules* **1977**, *10*, 949–954.

(58) Freeman, K. G.; Adamczyk, J.; Streletzky, K. A. Effect of Synthesis Temperature on Size, Structure, and Volume Phase Transition of Polysaccharide Microgels. *Macromolecules* **2020**, *53*, 9244–9253.

(59) Everaers, R. Constrained fluctuation theories of rubber elasticity: General results and an exactly solvable model. *Eur. Phys. J. B* **1998**, *4*, 341–350.

(60) Graessley, W. W. Statistical mechanics of random coil networks. *Macromolecules* **1975**, *8*, 186–190.

(61) Jedynak, R. Approximation of the inverse Langevin function revisited. *Rheol. Acta* **2015**, *54*, 29–39.