Structure and elasticity of model disordered, polydisperse, and defect-free polymer networks

Polymer networks are solids that can be obtained by cross-linking polymeric chains.¹ From rubbers² to hydrogels³ to biological networks,^4,5 these systems have countless industrial⁶ and biomedical⁷ applications. Understanding how the macroscopic properties of polymer networks, and in particular their elasticity, depend on their structure, chemistry, and topology is still a largely open problem.^8–13 To unravel these questions, numerical simulations^12–34 are an invaluable tool, as they allow a control on the structure and topology of the network that is impossible to achieve in experiments, making it possible to disentangle the effects of different microscopic contributions.

Since cross-linked polymer networks are systems with quenched (frozen-in) disorder, their properties depend on the way the network is formed. Therefore, when simulating a model polymer network, the first step is to choose an assembly protocol. One possibility is starting with a system of precursor polymers, which can be mono- or polydisperse, linear or branched, which are then cross-linked via some procedure. The cross-linking can be allowed to occur only between chain ends (end-linking)^13–20,25 or between any pair of monomers (random cross-linking).^12,21–24 Random cross-linking, however, is only efficient at melt densities, whereas end-linking suffers from kinetic limitations, as reaching a perfect (fully bonded) configuration requires a time that grows quickly as the length of chains grows.^14–16, Ad hoc methods can be used to increase the number of bonded sites.^15,16 Using these methods is, however, not completely satisfactory since the final structure in principle depends on the exact method that was used to force the formation of the bonds. Another option is to impose some lattice connectivity, like the diamond lattice.^26–31 Several of these “lattice networks” can then be randomly superimposed in order to obtain a disordered structure.^27,29,31 These systems, however, present an underlying ordered topology and monodisperse strand length, contrary to most experimental systems.

Here, we study a model of disordered, polydisperse, and defect-free (i.e., fully bonded) networks, which has originally been developed for the study of microgels^35–40 and later applied to the study of phantom,³² double,³³ and hyper-auxetic networks.³⁴ This model has been previously validated against experiments of microgels^35,40 and we therefore consider it as a viable computational counterpart of experimentally realizable hydrogel networks. The macroscopic properties of these networks, which are self-assembled via equilibrium simulations, depend only on a very small number of parameters. Moreover, the self-assembly procedure naturally gives rise to a disordered topology with a well-defined exponential chain length distribution, similar to that of randomly cross-linked networks^21,41 or resulting from step-growth polymerization.⁴² Here, we use molecular dynamics simulations to characterize the structure and elasticity of these systems, and show how these properties are related to each other.

In Sec. II, we give a detailed description of the simulation model and of the self-assembly procedure. In Sec. III, we study the structural properties of these networks, such as strand length distribution, radius of gyration of the chains, bond angle distribution, and structure factor. In Sec. IV, we analyze their elastic properties, connecting them to static observables. We conclude in Sec. V with a discussion of our results and of future perspectives.

We generate a polydisperse network via the method described in Ref. 35. This method was originally developed for MD simulations of microgels,^35–40 but it can be generalized to the case of bulk systems.^32–34 In contrast to network-generation approaches using direct simulations of cross-linking dynamics^12–25 or based on compenetrating lattices,^27,29,31 this method allows for an efficient generation of fully bonded networks by using a bottom-up self-assembly approach⁴³ based on a bond-swapping potential as detailed below.

The starting point is a mixture of two different species of patchy particles, i.e., spheres of identical size and mass decorated by a certain number of interaction sites (the “patches”) arranged in a regular configuration. The number of patches per particle is called the valence. We consider systems of volume V_init containing M_tot = M₂ + M_f particles, or monomers, with M₂ bivalent particles and M_f f-valent particles. Two patchy particles are reversibly attached to each other when at least two of their respective patches are overlapping. Only pairs formed by two bivalent particles or by a bivalent particle and an f-valent particle can be attached to each other, so that cross-links can be connected only through chains made of bivalent particles to form branched structures. In the following, we consider the cases f = 3 and f = 4.

The assembly is performed via molecular dynamics simulations (GPU implementation of the oxDNA software⁴⁶) run in the NVT ensemble at T = T_assembly = 0.05. At this low temperature, the system approaches the fully bonded ground state. Since the bonds can break and reform with an efficiency that is significantly improved by the bond-swapping potential, Eq. (4), the system can quickly reach equilibrium. Thus, for a given T_assembly, the properties of the final state are uniquely determined by only three parameters, as it will be shown below: the fraction of cross-links c = M_f/M_tot, the cross-link valence f, and the number density ρ_init = M_tot/V_init at which the assembly has been carried out. We note that the variation of the latter parameter can be regarded as a crude way for tuning solvent quality. In fact, from a mean-field-like perspective an increase in the assembly density corresponds to a stronger effective attractive monomer-monomer interaction.

Once the majority of the bonds $(>99.8%)$ are formed, the simulation is stopped and the particles that do not belong to the percolating cluster (at most 4% of the total in all the simulated systems) are removed. Although we chose for practical reasons to stop the reaction before reaching the fully bonded ground state of the system, reaching this state is in principle possible with a greater computational effort.

The so-obtained networks still contain a few dangling ends. As only a very small fraction of particles is part of these loose strands, we do not expect them to have a significant effect on the elastic properties of the network. However, due to their very long relaxation times, which are known to grow exponentially with the length of the strands,^47,48 they could lead to a slow relaxation of the monomer mean-squared displacement, which is analyzed below. To avoid these long relaxation times, we remove all the (simple and branched) dangling ends in the system, obtaining a perfect fully bonded network. We note that for systems with tetravalent cross-links, the removal of the dangling ends has the side effect of introducing some trivalent cross-links in the system, since when a dangling end is cut, a cross-link loses a bond. In this case, the cross-link fraction should be calculated as c = (M₃ + M₄)/(M₂ + M₃ + M₄). However, for all the systems considered, with the exception of c = 1%, where the number of dangling ends can become comparable to that of cross-linkers, we have M₃ ≪ M₄ and therefore the presence of the trivalent cross-links is negligible. By contrast, in the case of trivalent cross-links, the removal of dangling ends transforms cross-links into regular bifunctional particles. Regardless, in all cases at the end of the procedure, we obtain a network that is naturally disordered and almost fully bonded.^32,34

Molecular dynamics simulations of the network with frozen topology and connectivity are run using the LAMMPS package.⁵⁸ The system is initially allowed to relax to pressure P = 0 at constant temperature T = 1.0; then, NPT simulations are run at these T, P values. An equilibration run of 10⁷ time steps is followed by a production run of 2 · 10⁸ time steps. Temperature and pressure are kept constant by Nosé–Hoover chains (three thermostats and barostats),^59–61 insuring a correct sampling of the NPT ensemble. The integration time step is δt = 0.003 for all the simulations, and the relaxation time for the thermostat is chosen to be δt · T_damp = 0.3, whereas the relaxation time for the barostat is δt · P_damp = 3.⁵⁹ The three dimensions of the box, L_x, L_y, and L_z, are allowed to fluctuate independently, so that we have access to both volume and shape fluctuations, which are used to estimate the shear and bulk moduli G and K, respectively.³⁹ The total density ρ is thus defined as ρ ≡ M_tot/⟨V⟩, where V = L_xL_yL_z and ⟨·⟩ denotes an average over all configurations. We note that, here, M_tot refers to the number of particles after all the dangling ends have been removed (see Tables I and II).

In this section, we study the topological and structural properties of the networks.

After the relaxation at P = 0, the final mean density of the system, ρ ≡ M_tot/⟨V⟩, will in general be different from the assembly density ρ_init. In Fig. 2(a), we show ρ as a function of ρ_init for different values of c and f. One can see that for most of the systems with ρ_init < 0.85, we have ρ > ρ_init, i.e., the network contracts after switching the interaction from patchy to Kremer–Grest and letting the system relax to P = 0. This is due to the fact that the chains are stiffer with the patchy potential due to the directionality of the bonds, which is absent with the Kremer–Grest potential. The only exceptions are when f = 3, c = 5%, ρ_init = 0.5 and f = 3, c = 1%, ρ_init = 0.2, for which the network expands slightly (ρ < ρ_init). For ρ_init = 0.85, on the other hand, all the systems expand, as discussed below. Finally, we note that the curves for f = 4, c = 5% and f = 3, c = 10% superimpose perfectly, implying that the parameter controlling ρ is the mean valence F [Eq. (7)]. Indeed, it can be verified that both these systems have F = 2.1 (see also Tables I and II). This is expected for systems of patchy particles, where the mean valence F is known to control the equilibrium thermodynamic properties of the system.^53–56 

In order to study the spatial conformation of the strands, we consider the mean radius of gyration of single strands of length n, R_g(n). To calculate R_g, we consider as part of the strand only the bivalent particles and the bonds that connect them, excluding the cross-links. Thus, R_g(1) = 0, since for n = 1 we only have a single bead, and for n = 2 (two beads) we have R_g(2) = l_b/2, with l_b the bond length. In Fig. 3(a), we show R_g(n)/R_g(2) for all the simulated systems. As a consequence of Eqs. (8) and (9), networks with low c contain longer strands on average. In particular, for c = 1% we observe chain lengths up to n à 500. Fitting R_g(n) for all c values to the power law (n − 1)^ν in the range n − 1 ≥ 10 yields ν = 0.59 ± 0.02 for f = 3 and ν = 0.60 ± 0.02 for f = 4 [see Fig. 3(a)]. Within the accuracy of the data, these values are compatible with the Flory exponent for self-avoiding walks, i.e., ν = 0.59,¹ suggesting that in the considered range of n, the strands adopt conformations akin to those of linear chains in a good solvent. The inverse of the above metric exponent yields the fractal dimension of the individual strands, 1/ν = 1.7. Note that this fractal dimension is, in general, different from that of the whole network, and that for randomly branched polymers, this is predicted to be d_f = 2.⁶⁵ We also note that at fixed n, R_g(n) is basically independent of the cross-link concentration c and of the cross-link valence f, and it depends only weakly on the density of the system. This is shown in Fig. 3(b), where we plot the quantity $⟨Rg2(n)⟩n/⟨(n−1)2ν⟩n$⁠, which is proportional to the effective Kuhn length¹ of the chains. For all systems considered here, this is a monotonically increasing function of ρ and, discounting the numerical noise, a monotonically decreasing function of c. From Fig. 3(b), we also see that in the systems f = 4, c = 5% and f = 3, c = 10%, which have the same mean valence F, the strands have the same conformation. The same qualitative behavior is observed when considering the end-to-end distance (not shown).

The fact that R_g increases with density might be surprising, since in polymeric systems with free chains, R_g usually decreases with increasing ρ.¹ This behavior, however, can be understood qualitatively by analyzing the bond angle distribution of the strands, P(θ), where θ is the angle between bonded triplets of monomers. In Fig. 4, we report with solid lines P(θ) for c = 5% and f = 3 (a) and f = 4 (b) (the other systems show the same qualitative behavior). These are compared with the distribution obtained before changing the interaction potential from patchy to Kremer–Grest and allowing the system to relax to P = 0 (dotted lines). Before the system is brought to P = 0, P(θ) displays a peak at θ ≃ 142°, which originates from the patchy potential. Moreover, P(θ) drops to zero for θ ≲ 60°, which is the minimum allowed angle considering the excluded volume (WCA) interaction.

We note that before relaxation, P(θ) only depends on f and c and not on ρ_init. As the bonding potential is switched to the directional patchy one to the nondirectional FENE one and the system relaxes to P = 0, the curve becomes broader and the peak shifts to lower values of θ. At low ρ, the system can completely relax and P(θ) assumes a form identical to that of chains in a dilute solution (dashed lines; data from Ref. 66). By contrast, for high ρ, P(θ) retains a shape that is close to the shape it had before relaxation, signaling the presence of strong topological constraints, i.e., entanglements,¹ resulting from the non-crossability of the strands. These constraints prevent the bond angle from fully relaxing, as signaled by the survival of the θ ≃ 142° peak. We can see that for both systems, the average bond angle ⟨θ⟩ increases with density: This results in an increase in the effective persistence length¹ of the strands and of R_g(n), as shown in Fig. 3(b). This phenomenon is reminiscent of topological rigidification, when the system’s rigidity increases with increasing topological complexity as observed in knotted rings and star polymers,⁶⁷ although the limited range of strand lengths available here does not allow establishing the exact nature of this rigidification. We note that the persistence length of chains of patchy particles in solution decreases with ρ_init, but it does so only weakly (see  Appendix A); thus, this ρ_init-dependence cannot, by itself, explain the effect observed in P(θ).

The increase in the number of entanglements with increasing density can be more directly quantified using the method of primitive path analysis.^68,69 The procedure is as follows: (1) The cross-links are fixed in space, then (2) the intra-chain excluded volume interactions are turned off, while the interchain interactions are kept, and finally (3) the bonds of the segments are gradually shrunk, as when the system is cooled to T = 0. Since intra-chain excluded volume interactions are turned off, the chains are straightened; however, the interchain interaction still prevents the chains from passing through each other, and thus the topology is conserved. The network is thus reduced to a mesh of primitive paths. We then count the number of contacts between any two primitive paths and refer to it as the number of kinks, M_k. Under the assumption that two strands do not entangle more than once, M_k is proportional to the number of entanglements M_e, i.e., M_k ∝ M_e ≃ M_tot/4N_e, where N_e is the entanglement length.^1,68 One can see from Fig. 5 that the fraction of kinks M_k/M_tot increases with ρ and decreases with c and f. The ρ-dependence at high-density seems to be compatible with a power law of exponent 1.5, although the limited density range does not allow to draw any definite conclusion (see the dashed line reported in Fig. 5 as reference). It is important to stress that in our systems, an increase in ρ is not equivalent to the one obtainable by simply compressing the system (i.e., by collapsing the network), since the degree of entanglement of the network strongly depends on the assembly density ρ_init.

The results of the primitive path analysis are important for several reasons: First of all, from the methodology point of view, it is clear that in the protocol we employ to build disordered networks the number of topological kinks (and therefore of entanglements) can be controlled, in a statistical sense, by varying ρ_init. Leveraging this feature will make it possible to investigate fully entangled, polydisperse, disordered networks. Second, for the particular choice of parameters we made here, these results show that the investigated systems are at most lightly entangled as the order of magnitude of the number of entanglements (which is smaller than M_k) never exceeds that of the cross-links (also shown in Fig. 5). Finally, the increase in the number of entanglements with density may explain the apparent topological rigidification discussed above.

To gain insight into the global structure of the network, it is useful to study the structure factor.⁷⁰ We report the total structure factor S(q) in Fig. 6(a) and the structure factor of the cross-links S_f(q) in Fig. 6(b) for f = 3, c = 5% and different values of the monomer density ρ. The structure factors of the other systems display the same qualitative features. All the systems have rather large isothermal compressibility κ_T = S(0)/ρk_BT,⁷⁰ which grows with decreasing density. In a finite range of small wave-vectors q, S(q) is compatible with a power law, S(q) ∝ q^−α. This behavior results from the fractal structure and inherent (strand-)polydispersity of the network and also from the presence of holes with a wide size range, clearly visible in Fig. 1.^71–73 The value of α decreases with increasing monomer density; for the systems studied here, it was found that 0.5 ≲ α ≲ 1.3. For larger q, S(q) behaves similarly to that of a liquid, with a contact peak at q ≃ 8 ≃ 2π/σ, followed by periodic oscillations.

In Fig. 6(b), we show the structure factor of the cross-links, S_f(q), for the same systems. The qualitative behavior of S_f(q) is very similar to that of S(q), with a power-law decrease at small q, followed at larger q by a liquid-like oscillatory behavior (inset). However, the height of S_f(q) in the small-q limit is much lower than that of S(q) and also the oscillations have much smaller amplitude, showing that the location of the cross-links is more random than that of the monomers. Note that since two f-valent particles cannot bind, the distance of closest approach between two cross-links is ≃2σ and therefore the contact peak of S_f(q) is found at q ≃ π/σ.

In Fig. 6(c), we compare the total structure factors S(q) of the systems f = 4, c = 5% and f = 3, c = 10% for different values of the total density ρ. As noted above, these two systems have approximately the same mean valence F and therefore the same ratio ρ/ρ_init as well [Fig. 2(a)]. One can see that curves corresponding to similar densities superimpose almost perfectly, suggesting that S(q) is controlled either by ρ alone or by ρ and F. In Fig. 6(d), we compare S(q) for two systems with very similar densities but rather different F (2.05 and 2.20): The two curves are quite different, suggesting that the behavior of S(q) is determined by both F and ρ. These two structure factors are also compared with the S(q) of a solution of 50 chains of length 1000 at similar densities (dashed lines, data from Ref. 66), evidencing a rather different low-q behavior. The power-law behavior of the network’s S(q) at small q suggests that the structure of the network is significantly more heterogeneous than that of the solution at sufficiently large length scales or, equivalently, the inhomogeneities caused by the cross-links make the networks much more compressible than their chain-only counterparts.

In this section, we investigate the long-time limit of the mean-squared displacement of the particles and connect it to the shear modulus of the networks.

In Fig. 7(b), we show the localization length as a function of the monomer density ρ. The observed dependence of λ on ρ is roughly compatible with a power law with exponent −1, i.e., λ ∝ ρ⁻¹ although the curves display a steeper decrease at the highest densities, where local packing starts to play a more significant role. The same behavior is observed when considering the localization length of the cross-links or that of the bivalent particles, see Fig. 12 in the  Appendix. We note that also in this case, the curves for f = 4, c = 5% and f = 3, c = 10%, which have the same mean valence F, superimpose almost perfectly, confirming once more that the static properties of the system are controlled by F alone, as also shown in Sec. III. As the strand length distribution is independent of density, we ascribe the decrease of λ with ρ to topological constraints stemming from the non-crossability of the chains and possibly to a difference in the overall network topology.

Up to the present day, most simulation studies of polymer networks have considered structures that are ordered, monodisperse (with a unique strand length), or both.^26–31 Even in the case of networks produced through random cross-linking of precursor chains, it is very difficult, due to the slow chain dynamics, to reach a fully bonded state in which dangling ends are absent.^13–18,20 By contrast, here we used a method to generate disordered, polydisperse networks that are almost defect-free from the assembly of bivalent and f-valent (cross-links) patchy particles.^32,34,35 In the present study, all the networks satisfied more than 99.8% of all possible bonds and no more than 4% of the monomers were part of dangling ends, but it is in principle possible, with a larger computational effort, to obtain disordered fully bonded networks (100% of bonds formed) with no dangling ends. Moreover, the assembly process is an equilibrium one: As a result, the final structure of the network depends only on the percentage c of cross-links, on their valence f, and on the assembly density ρ_init. Another interesting property of these networks is that the distribution m_n of strand lengths is independent of density and only depends on f and c. However, we showed that the number of topological kinks depends rather strongly on density and increases approximately as a power law: This makes it possible to tune the entanglement length by changing the assembly density while keeping the average strand length N_s constant. We further note that the same model considered here can be used to study the influence that f and c have on the local network structure, e.g., ring statistics. Likewise, the model can be used to investigate how short-range attractions that mimic the solvent quality during the assembly will impact the final network.

By analyzing the dynamical properties of the networks on the level of the strands, we found that the dynamics of long strands is qualitatively consistent with the tube model. Moreover, we rationalized the behavior of the ratio between the localization length of the cross-links, $limt→∞⟨rf2(t)⟩$⁠, and the localization length of the middle monomers of the strands, $limt→∞⟨rmid2(t)⟩$⁠, finding a crossover between a phantom network model-like behavior at low density and a mean-field-like regime at high density where each particle can be described as being connected via an effective spring. Finally, we discussed the elasticity of the networks, showing that the values of the Poisson ratio we observe are in line with experimental values and that the peculiar ρ-dependence of the network shear modulus we observe, G ∝ ρ³, can be interpreted in light of the behavior of the localization length of the particles discussed above. We mention that the study of how the localization length and its relation to the shear modulus change when the network is subject to deformations is a promising direction for future work. Overall, the results presented here show that the assembly method we have used yields polymer networks that display realistic properties and thus can be used to model interesting phenomena where the polydisperse and disordered nature of the networks become important, such as polymer–nanoparticle composites⁸⁸ or double networks.⁸⁹ However, the peculiar scaling (or apparent scaling) behavior we found here cannot be straightforwardly understood by using classic polymer theories in view of the presence of short chains that are abundant in the systems investigated.³² Therefore, the results we have shown here should be taken as motivation for future theoretical efforts attempting to describe the behavior of realistic disordered polymer networks.

We thank Michael Lang for helpful discussions. We acknowledge financial support from the European Research Council (ERC Consolidator Grant No. 681597, MIMIC) and from LabEx NUMEV (Grant No. ANR-10-LABX-20) funded by the “Investissements d’Avenir” French Government program, managed by the French National Research Agency (ANR). W.K. is a senior member of the Institut Universitaire de France.

The authors have no conflicts to disclose.

Valerio Sorichetti: Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Andrea Ninarello: Data curation (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – review & editing (equal). José Manuel Ruiz Franco: Formal analysis (equal); Investigation (equal); Software (equal); Validation (equal); Writing – review & editing (equal). Virginie Hugouvieux: Supervision (equal); Writing – review & editing (equal). Emanuela Zaccarelli: Conceptualization (equal); Supervision (equal); Writing – review & editing (equal). Cristian Micheletti: Data curation (equal); Investigation (equal); Software (equal); Writing – review & editing (equal). Walter Kob: Conceptualization (equal); Supervision (equal); Writing – review & editing (equal). Lorenzo Rovigatti: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Methodology (equal); Project administration (equal); Software (equal); Supervision (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

In Fig. 12, we report the localization length of the bivalent particles [λ₂, panel (a)] of the cross-links [λ_f, panel (b)], and of the middle monomers of the odd-length strands [$⟨λmid⟩n$⁠, panel (c)]. We note that the total localization length λ is given by $λ2=cλf2+(1−c)λ22$ and that since c is small in the simulated systems, we have λ ≃ λ₂, i.e., the total localization length is approximately equal to the localization length of the bifunctional particles.