Entanglements in polymer nanocomposites containing spherical nanoparticles †

We investigate the polymer packing around nanoparticles and polymer/nanoparticle topological constraints (entanglements) in nanocomposites containing spherical nanoparticles in comparison to pure polymer melts using molecular dynamics (MD) simulations. The polymer–nanoparticle attraction leads to good dispersion of nanoparticles. We observe an increase in the number of topological constraints (decrease of total entanglement length N e with nanoparticle loading in the polymer matrix) in nanocomposites due to nanoparticles, as evidenced by larger contour lengths of the primitive paths. An increase of the nanoparticle radius reduces the polymer–particle entanglements. These studies demonstrate that the interaction between polymers and nanoparticles does not affect the total entanglement length because in nanocomposites with small nanoparticles, the polymer–nanoparticles topological constraints dominate.


Introduction
The dynamics of long polymers is controlled by entanglements, which are topological constraints imposed by the other chains. These can dramatically change the polymer viscosity, dynamics, mechanical and tribological properties. In this paper we explore how spherical nanoparticles affect rheology by studying the entanglements in polymer-nanoparticle composites in the cases when the polymer radius of gyration (R g ) is larger than [1][2][3][4] or of the order of the nanoparticle diameter (D). [5][6][7][8][9] The quality of nanoparticle dispersion 10,11 can play an important role on the polymer structure thus on polymer entanglements. In a well dispersed polystyrene (PS) chains/(PS) nanoparticles nanocomposite, 12 neutron scattering showed polymer chain expansion for polymer chains with radius of gyration larger than the nanoparticle radius (R), similar to the study of poly-(dimethylsiloxane)/polysilicate nanocomposites. 13 This has also been observed recently by simulations 14,15 and a thermodynamic model 16 but it is contrary to other recent studies of PS/silica 10,11,17 nanocomposites where polymers are unperturbed, and in poly-(ethylene-propylene)(PEP)/silica nanocomposites 5 where the polymer chains are contracted at very high nanoparticle loading; however, in some of the previous studies 10, 17 good nanoparticle dispersion has not been achieved and in some others 5,13 transmission electron microscopy (TEM) data were not reported, so the extent of dispersion is unknown.
Spherical nanoparticles can affect the primitive path and entanglement network of long polymers. 18 An increase of the entanglement polymer density is the origin of mechanical reinforcement in nanocomposites. [19][20][21][22][23] In particular it was shown using a slip-link model, 24 that in nanocomposites with ''bare'' fillers, a relatively small level of reinforcement was evidenced, 23 which is not verified in PS/silica and poly(methyl-methacrylate)/silica nanocomposites. 25 In addition, in ref. 23 was shown that the viscosity of the nanocomposite, Z, seems to be independent of the state of dispersion and can be predicted by the classical Einstein law: Z = Z 0 (1 + 2.5)f (where Z 0 is the viscosity of the pure polymer melt, and f is the nanoparticle volume fraction). However, the reinforcement is considerably higher when nanocomposites contain nanoparticles with grafted chains. 23 Such an observation was also reported recently by molecular dynamics simulations. 26 In addition, other parameters may also play a role on mechanical reinforcement such as size and shape of fillers, polymer matrix, interaction between fillers and matrix, and computer simulations 20,21,[27][28][29][30][31] have been used extensively to answer that fundamental problem.
The geometry of the nanoparticle (such as buckyball, graphene, nanodiamond 32,33 or nanorod [34][35][36] can also affect the polymer/ nanoparticle entanglement network. However in most of the entanglement network studies, a dilute nanoparticle regime [32][33][34][36][37][38] has been investigated, in which the polymer entanglement network, excluding the nanoparticles, remains unaffected. Only the MD study of nanoparticles R = 5 (where R g E D) 39 and Monte Carlo (MC) study by Termonia,40,41 have been performed at nanoparticle loadings above percolation. It is worthy to note that the MC study 40,41 is based on a body-centered-cubic (bcc) lattice in which the nanoparticles are fixed, and the polymer free volume remains constant irrespective of the nanoparticle size, which is not the real case in an experiment. 5,42,43 In this article, we investigate the polymer packing (free volume) in nanocomposites which contain spherical nanoparticles as fillers. We calculate the number of monomers between entanglements, entanglements per chain 44 and primitive path (the shortest path connecting the two ends of the polymer chain subject to the topological constraints) 45 in both polymer melts and nanocomposites of oligomers and weakly entangled polymers by using topological algorithms 44,[46][47][48] and applying different entanglements estimators. 47 The rest of this paper is organized as follows. In Section 2, the theoretical background is given for the entanglement analysis that is implemented in polymer melts and polymer nanocomposites. In Section 3, we discuss first the primitive path and entanglements of the polymer model used in this study. In polymer nanocomposites, we investigate the free volume and entanglements as a function of polymer molecular weight, volume fraction of fillers, interaction of polymers with fillers and nanoparticle radius in comparison to theoretical relations. Finally, in Section 4, conclusions are presented.

Estimators for entanglement length N e
In polymer melts of sufficiently long flexible chain molecules, neighboring chains strongly interpenetrate and entangle with each other. 49 Thus, the motion of polymers whose degree of polymerization is greater than the ''entanglement length'' N e is confined to a tube-like region.
The N e is determined by the estimator of Everaers et al. 45 (which we denote as classical S-coil), evaluated using the geometrical Z1 algorithm. 44,[46][47][48] This N e estimator is determined by statistical properties of the primitive path as a whole coil and evaluated for a given number of monomers N in the polymer chain as follows: where R ee is the end-to-end vector distance of a polymer chain and L pp is the contour length of its primitive path, the averages are taken over the ensemble of chains. Another estimator for the entanglement length can be used by measuring the number of interior ''kinks'' 44,46 which is considered to be proportional to the number of entanglements. The estimator on the number of ''kinks'', hZi, is denoted here as classical S-kink is given by: 44 In addition, there are modified estimators that provide an upper bound for N e , such as the modified S-coil, 47 but they tend to overestimate N e for weakly entangled chains: In order to eliminate the systematic errors that appear in the previous estimators 47 and to obtain an accurate N-independent value, we use an ideal N e estimator (M-coil), 47 which requires simulation of multiple systems of different chain lengths, using coil properties: where C(x) hR ee 2 i/R RW 2 (x) is the characteristic ratio 50 for a chain with x monomers, and R RW 2 (x) = (x À 1)r 0 2 is the reference mean squared end-to-end distance of a random walk, where r 0 = 0.967. Because of the dependence of C(x) on the number of monomers x, this estimator can be applicable to non-Gaussian chains. This non-Gaussian statistics of chains and primitive paths produces systematic errors in the old estimators for N e such as eqn (1). 51 The M-coil estimator converges faster than eqn (1)

Polymer melt
The chain and primitive path dimensions as calculated from the Z1 algorithm 44,46-48 for the polymer melts, of the semiflexible model used in this study, are presented in Table 1.
We depict the behavior of the M-coil estimator (eqn (4)) for the semiflexible Kremer-Grest (KG) polymer model studied, in Fig. 1, in comparison with results 47 of the fully flexible KG model. 53 From Fig. 1, it can be extracted that for the polymer model used in this study, N e E 59.7, while the value obtained from the S-coil estimator (eqn (1)) is 54.9. Other methodologies than the primitive path analysis, such as mean square displacement measurements (MSD), 51,54 can predict a much different N e value. 51,54 The MSD methodology assumes validity of the reptation model and assumptions made to come up with the i/L c = 2.02, where L c is the contour length of the polymer chain) and the packing length p 45 (the characteristic length at which polymers start to interpenetrate) decreases, thus the N e value unavoidly decreases for a semiflexible model 55 in comparison to the fully flexible Kremer-Grest model. 53 The glass transition of a polymer model which contains a bending potential (but not a torsional potential) is T g = 0.4. 52

Nanocomposites
For nanocomposites, we consider systems of spherical nanoparticles in a dense polymer melt. In the nanocomposite systems studied, a total number of N t = 23 600 monomers were used in a cubic cell with nanoparticles of radius R = 1 or 2 (in nanocomposites with nanoparticles R = 4 and polymer matrices N r 160, N t = 9440 monomers were used, whereas for R = 4 and polymer matrices N = 200 the total number of monomers was N t = 23 600). We define the nanoparticle (filler) volume fraction f in our simulations where D is the nanoparticle diameter and hVi is the total average volume of the nanocomposite simulation box during the NPT simulation. The mass of nanoparticle is m = 0.85pD 3 /6. Details of the nanocomposite systems studied (equilibration, length of the simulations) are given in ref. 14 and in Table 2.
In the next sections we investigate the effect of nanoparticle volume fraction, polymer-nanoparticle interaction and nanoparticle radius on the entanglement length N e and primitive path. In the nanocomposite systems studied we consider the case of the primitive path analysis for both the frozen particle limit, where nanoparticles are held fixed in space during the primitive path analysis, and the phantom particle limit. In the phantom particle limit nanoparticles are unable to restrict polymer motion on the time scales relevant to reptation dynamics. We approximate this limit, by removing the nanoparticles from the simulated system prior to the primitive path analysis. In the frozen particle limit, the N e contains two types of entanglements: polymerpolymer and polymer-nanoparticle entanglements, whereas in the phantom limit, it contains only polymer-polymer entanglements.
3.2.1 Effect of nanoparticle loading on polymer density. Nanocomposites with nanoparticles of R = 1 are similar to experimental systems of POSS nanoparticles (R = 1 nm) dispersed in polymer matrix (such as poly(ethylene-alt-propylene)). 56 Since we perform NPT simulations where the volume of the simulation box hVi fluctuates, the monomer density fluctuates and we can consider free volume effects. The quantification of free volume can be measured by calculating the net packing fraction (NPF): where hVi is the simulation average volume, N n the number of nanoparticles, N t is the total number of monomers, and s m is the monomer diameter. As can be seen in Fig. 2, the net packing fraction can increase (free volume decreases) as nanoparticle loading increases. This effect becomes stronger by dispersing large nanoparticles as also reported in ref. 57. While in ref. 57 only the dilute nanoparticle loading has been explored, we show the net packing fraction for a wider nanoparticle concentration range. The polymer nanoparticle interaction alters the free volume dramatically in the case of small nanoparticles R = 1. Since the smaller nanoparticles have a high interfacial area which increases with the nanoparticle loading, it is shown that the interfacial area is a factor that controls the packing fraction of such nanocomposite. However in nanocomposites containing larger attractive nanoparticles such as of size R = 4, the interfacial area diminishes and does not contribute to the net packing fraction. The polymer matrix has a slight effect on NPF calculations (shown in Fig. 2), an unentangled polymer matrix can reduce slightly the net packing fraction of the nanocomposites for all nanoparticle loadings. By dispersing attractive nanoparticles in the polymer matrix the polymer density around the nanoparticles increases. 14 By tuning the polymer-nanoparticle interaction e mp , not only the interface but also the dispersion and aggregation behavior  (which have been explored by molecular dynamics simulations 14,21,58,59 ) of the nanoparticles in nanocomposites is altered. It was shown that attractive nanoparticles of radii R = 1-4 can be dispersed in an unentangled or entangled polymer matrix (R g /R E 0.4-8) 14 in agreement to previous simulation studies. 21 However, for repulsive polymer nanoparticle interaction, poor dispersion is observed for R g /R E 0.4-4 (in our case radii R = 2 to R = 4), 14 which has also been observed experimentally for systems with weak interactions such as polystyrene-silica nanocomposite 10,17 and possibly for a repulsive nanoparticle nanocomposite such as PEP-silica. 5 3.2.2 Effect of nanoparticle loading on entanglement length N e . The modified S-coil, classical S-coil, and S-kink estimators (eqn (1)-(3)) for determining N e are applied on the nanocomposite systems, and shown in Fig. 3 for two nanoparticle loadings and nanoparticle size of R = 1. In general, both S-coil, S-kink (eqn (1) and (2) respectively) and modified S-coil estimators (eqn (3)) cannot correctly predict the N e value for all N exceeding N e . Eqn (1) and (2) converge to N e value slowly with increasing N, whereas eqn (3) tends to approach N e from above rather than from below, but overestimating its value for weakly entangled chains. However, they clearly show the effect of the volume fraction on the behaviour of the N e . In particular, at a nanoparticle loading of 5.5% (inset of Fig. 3), all the estimators show no changes in the N e value, between the phantom and frozen particle limits, for N = 200. In this case, the nanoparticles do not affect the primitive path of the polymer chains. By increasing the nanoparticle loading above than 5.5%, to 14.5%, it can be seen that the values of all three estimators reduce in the frozen particle limit, in comparison to phantom particle limit. This shows that at such loading, nanoparticles are additional topological constraints, that increase the length of the primitive path and unavoidably decrease the entanglement length. Moreover, (in particular in Fig. 3 from 5.5% to 14.5% loading), it can be seen that the N e in the phantom limit (red symbols) increases for all polymers.
This indicates that small attractive nanoparticles can alter the polymer primitive network. A similar trend, but to a smaller extent, has been also observed at a nanoparticle loading of 10.3% (results not shown). By increasing the nanoparticle loading, there is a higher deviation between the phantom and frozen particle limits, and there is a further decrease of the total entanglement length. In addition to previous estimators, the M-coil estimator (eqn (4)) is used for nanocomposites with small nanoparticles (R = 1) in order to have an N-independent estimation of N e . Clearly as can be seen in Fig. 4, by increasing the volume fraction of nanoparticles dispersed in the polymer matrix, the N e is reduced due to the contour length of primitive path, L pp , increase. It is noted that these N e values contain two types of entanglements: polymer-polymer and polymer-nanoparticle entanglements. Results from molecular dynamics 36 and dissipative dynamics 34 simulations of nanorods embedded in a polymer matrix indicate that the inclusion of nanorods into the polymer matrix does not significantly alter the polymerpolymer entanglements network at low nanoparticle loading, instead it creates additional topological constraints of polymernanorod origin.
Moreover, we depict in Table 3, the N e values (and number of ''kinks'' hZi) in the frozen particle limit as calculated by the S-coil, S-kink and S-modified estimators for long polymers N = 200, in nanocomposites at different nanoparticle loading, for both repulsive (Re) and attractive (A) nanoparticles. It can be seen that outside the error margin there is no difference in these N e (and hZi) values. The N e values estimated from M-coil agree with those estimated from S-coil for long polymers (N = 200), within the error margin. Similarly, for nanocomposites consisted of a matrix of polymers N = 100, the type of nanoparticles do not alter the N e in the frozen particle limit (results not shown). Thus in the case where small nanoparticles are dispersed in a polymer matrix, the polymer-nanoparticle interaction does not  play any role on the primitive path. The topological constraints created by nanoparticles seem to dominate the entanglement network even if the polymer dimensions can be altered 14 by the polymer nanoparticle interaction. We also report in Table 3 the number of kinks hZi in the phantom limit. We can see that in the phantom limit, hZi decreases with nanoparticle loading (whereas it increases in the frozen limit due to polymer nanoparticle entanglements). This shows that nanoparticle loading reduces the polymer-polymer entanglements for nanocomposites containing nanoparticles of radius R = 1.
The concept of entanglement length is useful because it relates changes in structure to rheological properties. 45,55,60 In polymer melts and semidilute solutions, a temperature and concentration dependent material constant, the plateau shear modulus G 0 N , which is of the order of 10 6 Pa, or five orders of magnitude smaller than the shear modulus of ordinary solids, is related to rheological entanglement length, N rheol e , by eqn (6): 45,61 where, r is the monomer density, and k B T is the thermal energy. The rheological entanglement length N rheo e should be equal to N e calculated from eqn (1) (the classical S-coil estimator) for loosely entangled polymer chains. 61 In polymer nanocomposites the validity of eqn (6) is unclear, and especially at high nanoparticle loading; however, a dependence of the plateau modulus G 0 56,62 on the filler degree f has been observed for repulsive nanoparticle nanocomposites when R g 4 R filler (such as PEP-POSS, PI-POSS 56 ), [63][64][65] where f (f) is given by: 56,64 f (f) = 1 + [Z]bf + a 2 (bf) 2 + a 3 (bf) 3 + Á Á Á where, Z = 2.5, 66,67 a 2 = 14.1, 68 and b is an effectiveness factor. 56 For b = 1 and a 3 = 0, eqn (7) leads to the Guth-Gold relation, 68 while if additionally a 2 = 0 the Einstein-Smallwood relation 64,66,67 is obtained. Another model for estimating the plateau modulus has been proposed by Eilers 69 where, The addition of small nanoparticles in the polymer matrix decreases the N e value, as shown in Fig. 4, thus the plateau shear modulus G 0 N increases, according to eqn (6), since it is approximately inverse proportional to N e . In Fig. 5 we depict a comparison between the plateau modulus experimental measurements, 56 theoretical predictions 56,69 and simulation data for the N e (f = 0)/N e (f) ratio, at different nanoparticle loadings. The ratio of N e (f)/N e (f = 0), as calculated by our simulations, decreases with the nanoparticle volume fraction. At a volume fraction, f = 24.2%, there is approximately 60% decrease in N e mainly due to the polymer-nanoparticle entanglements. Instead in a polymer nanorod composite such a decrease in N e appears at a much smaller nanorod volume fraction, f nanorod E 11%. 34 It seems from Fig. 5 that the N e decrease in nanocomposites with small nanoparticles follows quantitatively the theoretical trends of Guth-Gold relation (eqn (7)), however eqn (7) is not necessarily proportional to the ratio of N e values. The Einstein relation 66 is invalid for a such small nanoparticle composite, in contrast to nanocomposites with bare spherical nanoparticles studied through the slip-link model. 23 Small nanoparticles, such as R = 1, can reinforce polymers effectively. All three estimators in Fig. 5 show that the mechanical reinforcement effect in nanocomposites can be induced by the change of primitive path network due to the additional topological constraints created by small nanoparticles.
3.2.3 Effect of nanoparticle radius on entanglement length N e . Increasing the radius of the nanoparticles at a constant volume fraction decreases the surface area to volume ratio of the nanoparticles. The effect of the nanoparticle radius on the N e (f) from our simulations is depicted in Fig. 6. We observe that by increasing the nanoparticle radius to R = 2 we can see a decrease to the discrepancy between the phantom and frozen particle limit. In particular, at f = 25.8%, the S-coil estimator predicts for polymers (N = 200) a value of N e = 66.1 AE 2.4 in the frozen limit (from M-coil estimator N e = 64) and N e = 72.5 AE 1.9 in the phantom limit. Specifically, for nanocomposites with nanoparticles R = 4, the N e estimators in the phantom and frozen limit are indistinguishable. This implies that nanoparticles of R = 4 do not alter the underlying polymer network for polymer lengths used in our study. Tuteja also found that N e Fig. 4 Dependence of N e using eqn (4) in the frozen particle limit in nanocomposites with attractive small nanoparticles of R = 1 for different nanoparticle loading: (i) 10.3% (triangles), (ii) 14.5% (squares), (iii) 18.2% (circles), (iv) 24.2% (diamonds), (ii) 36% (star symbols). Inset: Dependence of L pp for different polymers in the frozen particle limit. The vertical lines show the N e extracted values at each nanoparticle loading.   is unaffected by nanoparticles of R = 5 nm at low nanoparticle loading (f = 8%). 70,71 The interparticle distance of nanoparticles is: ID = D((2/pf) 0.333 À 1). For all the nanoparticle volume fractions studied ID o R g (R g E 8 for N = 200). In that regime, there is no change in N e , in the frozen limit, when R g E D, whereas it changes only if R g c R (figure in ESI †). Also we can observe that in nanocomposites with nanoparticles of R = 2 in a polymer matrix N = 200 (see Fig. 6) the N e in the phantom limit is enhanced with respect to its polymer melt value. In order to investigate further the polymer path network, we calculated the polymer tube diameter 37,72 and depict it in Fig. 7: in which L pp (ph) is extracted in the phantom limit. As can be seen in Fig. 7, for attractive small nanoparticles, the tube diameter increases with the nanoparticle loading, whereas for repulsive nanoparticles it remains constant in agreement to the data by Li et al. 39 This increase means that such small nanoparticles (R = 1) do alter the polymer network, and the polymers disentangle with nanoparticle loading. This implies that in the case of attractive nanoparticles, the N e in the phantom limit is increased, as observed in Fig. 3. While the total N e (in the frozen limit) of the nanocomposite decreases with the nanoparticle loading (see Fig. 4), the polymer-polymer entanglements convert to polymer-nanoparticle entanglements approximately for nanoparticle loading f Z 15%. This disentanglement effect does also appear in thin polymer films, 73 under cylindrical confinement, 74,75 on a bare flat surface 76 and on the vicinity of large repulsive spherical nanoparticles at a high nanoparticle loading. 39 Furthermore, the polymer chain dynamics can also be affected by the nanoparticle volume fraction. Since by increasing the nanoparticles loading the tube diameter increases, that can enhance the polymer chain's diffusivity. However, direct studies of diffusion remain quite difficult to study the slow reptational dynamics of nanocomposites using molecular dynamics simulations.

Conclusions
The polymer density, polymer/polymer and polymer/nanoparticle topological constraints (entanglements) of polymers in melts and nanocomposites containing spherical nanoparticles were   investigated by means of molecular dynamics simulations. We applied different N e (N) estimators for the calculation of the number of entanglements in our systems, and extracted the N-independent entanglement length N e . We observe that the total N e decreases even with low volume fraction of small nanoparticles, and significantly for f Z 25%. This decrease of N e , in the nanocomposite, originates from the polymer/ nanoparticle entanglements, because the contour length of the primitive path, L pp , increases with the addition of nanoparticles. In order for polymer nanoparticles entanglements to be formed, the polymers need to be substantially larger than the nanoparticles in order to wrap around them, and in that case the nanoparticles act as topological constraints. Interaction between polymers and nanoparticles does not affect the total entanglement length when there is good nanoparticle dispersion. For the case of attractive small nanoparticles (such as R = 1) the polymer-polymer entanglements decrease (increase of tube diameter) due to the expansion of the polymer chains for f Z 15%. This effect on the polymer network is enhanced by the nanoparticle loading. Instead for the case of repulsive nanoparticles the tube diameter remains the same up to E24% nanoparticle loading in agreement with previous studies. 39