Polymer and spherical nanoparticle diffusion in nanocomposites

Long polymer chains diffuse in dense systems following the reptation model of de Gennes.¹ The dynamics of long polymers is controlled by entanglements, which are topological constraints imposed by the other chains. However, in polymer nanocomposites,² there are more topological constraints due to the dispersed nanoparticles (NPs) in the polymer matrix that can influence the polymer primitive path,^3–7 thus the polymer dynamics. Polymer dynamics can be altered by the characteristics of nanoparticles such as size, shape, aspect ratio, surface area, volume fraction of nanoparticles, and the nature of the interactions between the nanoparticle and the polymer matrix.^8,9 There are several theories that have been used to predict how the polymer diffusivity is impeded by impenetrable particles.^10–19 The first theory that was developed for a two phase system (particle and polymer) was introduced by Maxwell.¹⁰ It predicted a monotonic decrease in monomer diffusivity, as the particles act as obstacles. More diffusion models^11,12,14 were developed following the same obstruction approach which also included the effect of particle geometry on polymer diffusivity. While most models for polymer nanocomposites predict a decrease in polymer diffusivity due to tortuosity, only recent models^19–22 are based on the interaction of polymer chains with the nanoparticle surface. Polymer tracer diffusivity decreases in nanocomposites of large silica nanoparticles dispersed in polystyrene (PS)²³ or (PMMA).²⁴ The tracer diffusivity decreases across a wide range of polymer molecular weight and nanoparticle size^23,24 and unexpectedly decreases more strongly at higher temperatures.²⁵ In addition, Richter and co-workers explored the dynamics of entangled^26,27 and unentangled poly(ethylene-alt-propylene) (PEP) polymers²⁸ in nanocomposites by neutron spin echo (NSE) experiments. For entangled polymers, it was found that the initial Rouse relaxation rate was unaffected by the hydrophobic silica nanoparticles. A crossover from classical chain entanglement to nanoparticle entanglement, for which chains motion is hindered by nanoparticles rather than surrounding chains, has been observed experimentally.²⁶ The confined polymer was found to be surprisingly mobile and did not exhibit glassy behavior.^29,30 For unentangled poly(ethylene-alt-propylene) (PEP), there is a strong deceleration of the polymer center of mass motion, but even at the highest silica volume fraction, the segmental chain dynamics is not affected.²⁸ 

From a simulation point of view, there are several studies that use the Kremer-Grest model³¹ to investigate the polymer dynamics. However searching for the reptation power laws¹ as predicted by reptation theory, through molecular dynamics simulations, is difficult.^32–38 For example, by using the Kremer-Grest model,³¹ the reptation regime cannot be reached for very long polymers.³³ By using Monte Carlo (MC) simulations and the bond fluctuation model,^39,40 it was found that there is a really long crossover region from the nonentangled dynamics to reptation dynamics.^41,42 In a study of polymer chains which are confined between two walls,⁴³ it was observed that confinement influences both the static⁴⁴ and dynamic properties of the polymers.^43,45–50 For a self-avoiding walk chain that diffuses between randomly distributed impenetrable fixed obstacles, the chain dynamics is slower than even reptation.⁵¹ For such of random walk polymers of length N which diffuse in a network of obstacles, the center of mass diffusion, D, varies as N^−2.4.⁵² In a recent study on lattice MC simulations of long polymers in nanocomposites,^44,53 a decrease in polymer chain diffusion was observed when the nanoparticles are fixed in the polymer matrix, with radius (R) smaller than polymer radius of gyration, by both geometric (confinement by fillers)^44,53 and energetic (monomer-nanoparticle interaction) effects.^53–60 Nanoparticle diffusion in a medium is predicted theoretically by the Stokes–Einstein formula⁶¹ in a dilute nanoparticle regime. Nanoparticles smaller than the tube diameter^62–65 do not follow the Stokes-Einstein relation.^62,66–74 To the best of our knowledge, there is no simulation study that investigates the polymer (except Ref. 75 for unentangled polymers) and nanoparticle diffusion in nanocomposites, in a non-dilute nanoparticle regime. In this article, we explore how spherical attractive nanoparticles and polymers diffuse in dense melts in the cases when the polymer radius of gyration (R_g) is larger than^76–79 the nanoparticle diameter (2R), in comparison to experiments. The rest of this paper is organized as follows. In Section II, the methodology and simulation details of the present study are described. Subsequently in Section III A 1, the diffusion coefficient of unentangled polymers, normalized to its bulk value, is presented as a function of nanoparticle loading. In Section III A 2, the polymer diffusion of weakly entangled chains is compared to that of experiments.²⁴ In Section III B, the diffusion of small nanoparticles of different radii, in the various molecular weight polymer matrices, is calculated, for a range of nanoparticle loading, and compared to the Stokes-Einstein formula. Finally, in Section IV, conclusions are presented.

In this study, NPT molecular dynamics simulations were performed using the GROMACS package. A Kremer-Grest polymer model³¹ was used to incorporating an intrinsic stiffness through a cosine harmonic bending potential,^4,80 which acts on three consecutive beads along the chain

where $θijk$ is the bending angle between three consecutive bonds. The equilibrium value $θ0=109.5°$⁠, and the bending constant $kθ=25$⁠.⁸⁰ For the polymer model used in this study, N_e obtained from the S-coil estimator is 54.9,⁴ and the tube diameter is $dT≈10.3σm$⁴ using the primitive path analysis method.⁸¹ More details regarding the polymer model used can be found in Ref. 78. The nanocomposite systems contain nanoparticles represented by single spheres dispersed in the dense polymer melt. The Lennard-Jones (LJ) potential, $VijLJ$⁠, acting along the line between the centres of mass of two particles force⁸² is given by

where $ϵij$ is the interaction energy between particle i and particle j (for monomers: $ϵm=σm=mm=1$⁠, for nanoparticles: $ϵp=1$⁠, $σp=2R$⁠, $R/σm=1−2$⁠, $mp=0.8543πR3$⁠), and r_ij represents the distance between particles i and j. For monomer (m)- nanoparticle (p) interactions, the combination rules $ϵmp=(ϵmϵp)1/2$ and $σmp=(σm+σp)/2$⁸² were used. The monomers (or nanoparticles) are modelled with the repulsive only part of Eq. (2), shifted and truncated with a cutoff radius at $rc=21/6σm$ (⁠$rc=21/6σp$ for nanoparticles). For the monomer–nanoparticle interaction also the attractive part of the LJ potential (Eq. (2)) is considered and the cutoff radius is $rc=2.4σmp$⁠. The interaction energy between monomers and nanoparticles is $ϵ=kBT$⁠. A total number of N_t = 23 600 monomers were used in a cubic cell with nanoparticles of radius R = 1, 1.5, or 2. We define the nanoparticle (filler) volume fraction ϕ in our simulations as $φ=π(2R)3Nn6⟨V⟩$⁠, where 2R is the nanoparticle diameter and $⟨V⟩$ is the total average volume of the nanocomposite simulation box during the NPT simulation. The mass of nanoparticle is $m=0.85πD3/6$⁠. For equilibration of polymer melts, the fast push off method was applied.⁸³ The pressure calculated for the N = 200 polymer melt was $P∗=Pσm3/ϵm$ = 4.864. That pressure was used to perform all the nanocomposite system simulations in the NPT ensemble. In nanocomposites, the starting structures were created by an ensemble of polymers with N = 200 (or N = 400) and nanoparticles inserted at random positions within a large simulation box. Subsequently, the NPT ensemble was used at melt pressure $P∗=4.864$ to gradually squeeze the simulation box to a dense polymer nanocomposite. Specifically, for nanocomposites with polymer matrices N = 200 and 100, the equilibrium time was $6⋅105τ$ and $3⋅105τ$⁠, respectively, which is sufficient time to evolve the entanglement density. For nanocomposites with unentangled matrix (N = 50, 20, 10), the starting structure was the equilibrated nanocomposite with matrix N = 100, and then a further equilibration of $105τ$ and $1.5⋅103τ$ was performed for systems with N = 50 and N = 20, 10, respectively. The polymer radius of gyration reached a constant value with time. We also checked that the mean square internal distances for all nanocomposites reach a Gaussian limit.⁸⁴ The duration of the simulation production runs was between $6⋅103−104τ$ for unentangled polymers (N = 10, 20),$3⋅105τ$ for polymers with N = 50, and $4.5⋅105−8⋅105τ$ for polymers with N = 100. Details of the nanocomposite systems studied are given in Table I. A snapshot of the simulated polymer melt containing spherical nanoparticles is shown in Figure 1.

The length of the simulation cell, which is $L≈31−34.6σm$ (depending on the nanoparticle volume fraction ϕ: see Table I), was always larger than the end-to-end distance, R_e, of the polymer chains plus the cutoff distance (⁠$rc=21/6σm$⁠), thus the polymers do not interact with its periodic image. In a previous simulation study of Desai,⁸⁵ a smaller simulation cell length was used (⁠$L≈22σm$⁠) studying diffusion of polymers of similar length (N = 80).

To set the temperature at $T=kBT/ϵ=1$⁠, the Langevin thermostat⁸⁶ was used with a friction constant $Γ=0.5τ−1$⁠. To set pressure at $P∗=4.864$⁠, the Berendsen barostat was used with time constant $τp=2τ$⁠. The equations of motion were integrated using the leapfrog-algorithm⁸⁷ with a time step equal to 0.004τ, where τ is the Lennard-Jones time. The glass transition of a polymer model which contains a bending potential (but not a torsional potential) is T_g = 0.4.⁸⁰ The monomer–nanoparticle radial distribution function g_mp(r) exhibits a three layer structure.⁷⁸ The high monomer density of the layers establish a well defined interface between nanoparticles and polymers whose structure differs from that of the amorphous polymer melt. Also, the nanoparticle loading increases the monomer density of polymers in contact with the nanoparticle surface.⁷⁸ 

The dynamic properties of polymer chains are estimated through the self-diffusion coefficient D₀ which is calculated through the mean square displacement⁸² in the regime where the mean square displacement increases linearly with time

where $⟨|𝐫i(t)−𝐫i(0)|2⟩$ is the time dependent displacement of the monomer averaged over time and monomers of the ensemble (r_i is the position of a monomer i).

The diffusion of unentangled tracer polymers in an entangled matrix⁸⁸ follows the Rouse model. In a polymer nanocomposite containing non-attractive large nanoparticles⁷⁵ (where diameter 2R is larger than R_g), unentangled polymer diffusion decreases with nanoparticle loading. However in our study, we investigate spherical nanoparticles with attractive interactions with the polymers that are well dispersed⁷⁸ in the polymer matrix. We show in Figure 2 the polymer diffusivity, normalized to its bulk value, as a function of the nanoparticle loading for $Rg>2R$ (2R = 2). We can see that at nanoparticle loading higher than the dilute regime (⁠$ϕ>5.5%$⁠), the polymer diffusivity is lower than that predicted from the Maxwell model. We observe that the polymer diffusion decreases with the addition of nanoparticles in the matrix. This is due to, by adding nanoparticles in the matrix, the increased interfacial area around the nanoparticles and more monomers exist into the interfacial area.⁷⁸ The same trend appears for unentangled matrices of N = 10, 20, 50. Similarly, a reduction of unentangled polymer diffusion with the nanoparticle loading has also been observed in the case of very strong attractive polymer-filler interaction (⁠$ϵmp=4−10ϵm$^59,85,89). The polymer chain dynamics of the nanocomposites containing attractive polymer/nanoparticle interaction is largely affected by the nanoparticles,⁹⁰ in nanoparticle loading larger than the dilute regime (⁠$ϕ=5.5%$⁠) as in the experimental systems of poly(ethylene glycol) (PEG) and poly(butylene oxide) (PBO) chains⁹¹ or (poly(2-vinylpyridine) (P2VP)) chains⁹ with silica nanoparticles. Since the PEG polymers interact attractively with neat silica nanoparticles, a fraction of PEG chains was attached to the nanoparticle surface.³⁰ 

We show in Figure 3 the diffusivity (as calculated from Eq. (3)) of weakly entangled polymers (N = 100) as a function of the nanoparticle loading. We observe that the polymer diffusion decreases with the addition of nanoparticles in the matrix. In particular, in the nanocomposites studied at 24.2% nanoparticle loading, the ratio, of polymer diffusivity, D/D₀ = 0.3, whereas in a Monte Carlo (MC) study of athermal nanocomposites,⁴⁴ the ratio D/D₀ decreases to $≈0.28$⁠, for weakly entangled polymers, at nanoparticle loading 30%. However, it is worth noting that in the MC study, the nanoparticles are fixed into the simulation cell, thus they are immobile, contributing to the decrease of polymer diffusivity, whereas in our work, they are mobile and not fixed in the simulation cell. In addition, in the case of repulsive nanoparticles,^54,59,85 the polymer diffusion coefficient initially increases with increasing nanoparticle loading (up to 5%),^54,85 it reaches a maximum, and decreases further to a D/D₀ = 0.8 ratio, at nanoparticle loading (10%). While in the initial dilute nanoparticle loading, the increase of polymer diffusion (in comparison to its bulk value) reflects the repulsive polymer-nanoparticle interaction, at a higher loading (10%), a reduction of the polymer diffusion is observed, due to simply geometrical reasons, and is affected by the presence of tortuous path.⁸⁵ We can see that the simulation data agree with the experimental system of PMMA tracer diffusion,²⁴ at high nanoparticle loadings. However, in the dilute nanoparticle regime, there is discrepancy to the data. In the experimental nanocomposite, tracer diffusion²⁴ collapses onto a master curve when the neat homopolymer diffusion (D/D₀) is plotted versus a confinement parameter defined as the interparticle distance relative to the polymer radius of gyration (ID/2R_g)^23,24,92 where the interparticle distance ID is given as $ID=2R[(2πϕ)1/3−1]$⁠, where ϕ is the volume fraction of nanoparticles. However in our data, the mechanism of reduction of diffusion in our simulation data is not due to entropic barriers but due to nanoparticle interfacial area.^93,94 Although the, experimental, nanocomposite system²⁴ considers immobile nanoparticles (of radius 6.5-25 nm) that are different than our simulated small nanoparticle composites, this is the only studied nanocomposite system which contains an attractive nanoparticle polymer interaction, and thus the nanoparticles are well dispersed in the polymer matrix. In order to estimate the local segmental dynamics, the bond autocorrelation, C_b(t), was calculated according to the equation $Cb(t)=⟨P2[b(t)⋅b(0)]⟩$⁠, where P₂ is the second Legendre polynomial, b(t) is a unit vector aligned along the bond of a polymer, and the angular brackets indicate an average over all polymer bonds in the system. C_b(t) measures the local segmental dynamics of polymer chains, equally weighting reorientation motions in all spatial directions. It was found that the segmental dynamics decreases (thus the relaxation time increases) with nanoparticle loading and with increasing nanoparticle radius (see Figures S1 and S2 in the supplementary material).

In a simplified model where polymer chains were modeled as spheres diffusing down cylindrical pores,⁹⁵ it was shown that excluded volume effects can account for the observed reduction in polymer diffusion for values of $ID/2Rg>5$⁠. In that region, the excluded volume model predicted a pseudomaster curve for the reduced diffusion coefficients, in good agreement with experimental results^23,92 at higher nanoparticle loading. In addition, in nanocomposites of poly(ethylene glycol) (PEG)-tethered silica (SiO₂) nanoparticles and poly(methyl methacrylate) (PMMA), polymer chain dynamics undergoes a continuous transition from bulk-like behavior at low nanoparticle loading to confinement behavior at intermediate nanoparticle loading and ultimately to glassy behavior at high nanoparticle loading.⁹⁶ 

Nanocomposites containing nanoparticles have been investigated experimentally such as systems of polyhedral oligomeric silsesquioxane (POSS) nanoparticles (R = 1 nm) dispersed in poly(ethylene-alt-propylene) matrix²⁷ or cadmium selenide nanoparticles dispersed in PS matrix.⁶⁶ While in other theoretical studies,^{70,85,97–99} only the dilute nanoparticle loading has been explored, we investigate, for the first time, the nanoparticle diffusion for a wide nanoparticle concentration range. It is known that the dynamics of spherical particles (colloids) in a fluid can be predicted by the classic continuum Stokes-Einstein relation. For a spherical particle of radius R in a fluid of viscosity η, the diffusion coefficient D_SE of the particle is given by^59,61

where η is the viscosity of the dense polymer melt (Kremer Grest model), given by Ref. 100 and for long polymers (N = 200, 400) by Ref. 101. However, in the case of small nanoparticles dispersed in polymer melts^{27,79,102,103} and solutions,^104,105 the validity of Equation (4) is questionable as observed by experiments⁶⁶ and theory.^69,106 The nanoparticle diffusivities of the simulated nanocomposite systems are calculated from the mean square displacement measurements given by Equation 3, where $D0=16limt→∞ddt⟨|𝐫i(t)−𝐫i(0)|2⟩$⁠, where $⟨|𝐫i(t)−𝐫i(0)|2⟩$ is the time dependent displacement of the nanoparticle averaged over time and nanoparticles of the ensemble.

It can be seen from Figure 4 the nanoparticle diffusivity for different polymer matrices and nanoparticle loading. Clearly it can be seen that in the dilute nanoparticle regime (⁠$ϕ=5.5%$⁠), small nanoparticles (R = 1) diffuse much faster than the SE prediction⁶⁴ for matrices $N≥50$⁠. The Stokes-Einstein relation can be valid, for nanoparticles diffusivity, in nanocomposites containing very short (unentangled) polymers at the dilute regime in agreement with experiments.⁷² The diffusion of such small sized nanoparticles reaches a plateau in entangled polymer matrices for $N≥100$ in agreement with the Generalized Langevin Equation (GLE) theory^62,64 and theoretical predictions by de Gennes¹ who argued that bulk viscosity does not capture the behavior of surrounding flows near NPs, and thus the NPs diffusion is decoupled from the SE relation (Eq. (4)). In addition, by increasing the nanoparticle loading, nanoparticles diffusivity is decreased dramatically in all polymer matrices. The same plateau region is obtained for all nanoparticles loading in entangled polymer matrices. Since the nanoparticle diameter is smaller than the entanglement strand, such dispersed nanoparticles can diffuse inside the mesh of polymers, faster in comparison to the SE predictions, due to local viscosity.⁶⁴ Such small fillers provide the best reinforcement effects^{3,4,107–109} in all the loading range considered and diffuse faster than the monomers. In the inset of Figure 4, the nanoparticle diffusivity is scaled with the SE prediction (Eq. (4)) for different polymer matrices. It can be seen, from all the data, that such nanoparticles, smaller than the tube diameter d_T (⁠$dT≈10.3σm$ for our polymer model), deviate more from the SE predictions by increasing the molecular weight of the polymer matrix and nanoparticle loading. The discrepancy is large for weakly entangled (and even more for entangled) polymer matrices due to the viscosity scaling with molecular weight $η∝N3.4$ in the entangled regime. In an experimental study of silica nanoparticles in entangled poly(2-vinylpyridine) (P2VP) melts, in which 2R/d_T = 3, the diffusion of silica nanoparticles in P2VP⁹ is well-described by the hydrodynamic Stokes-Einstein relation as can be seen in the inset of Fig. 4.

The nanoparticle loading has a high impact on the nanoparticle diffusivity decrease, especially for unentangled polymers, as can be seen in Fig. 5. By increasing the nanoparticle radius, the nanoparticle diffusivity decreases as can be seen in the inset of Fig. 5 for different polymer matrix molecular weights. It is worth noting that a similar plateau region in the entangled regime appears in experiments, and nanoparticles diffusivity will approach the Stokes-Einstein predictions when the radius is approximately 4-5 times of the matrix tube diameter, as was evidenced in poly(n-butyl methacrylate)/gold nanocomposites.⁷³ 

The nanoparticle and unentangled polymer dynamics in nanocomposites containing spherical nanoparticles were investigated by means of molecular dynamics simulations. In conclusion, we have found that attractive interactions between the nanoparticles and the polymers appear to hinder chain motion, especially at higher nanoparticle volume fractions. Contrary to previously studied models containing unentangled polymers, we also calculate the diffusivity of weakly semiflexible entangled polymers and predict a better comparison to experimental data of nanocomposites, in which there is a good dispersion of silica nanoparticles in a polystyrene matrix. We reported that the polymer diffusivity decreases with nanoparticle loading due to the increase of the nanoparticle surface area. Also other parameters may affect polymer diffusion, such as polymer radius of gyration with nanoparticle loading,⁷⁸ polymer entanglements,⁴ nanoparticle diffusivity, strength of monomer-nanoparticle interaction, and confinement.

We also show that small sized nanoparticles can diffuse much faster than that predicted from the Stokes-Einstein relation in the dilute regime. However, for large volume fractions, the nanoparticles diffusivity decreases due to both nanoparticle–polymer surface area and nanoparticle loading.

See supplementary material for the bond autocorrelation decay function and mean square displacement of polymers in different nanocomposites.

This research was funded by the EPSRC/NSF Materials Network Program Nos. EP/5065373/1 (EPSRC: N.C. and A.K.) and DMR-1210379 (NSF: K.I.W. and R.J.C.). Support was also provided by ACS/PRF 54028-ND7 (RJC) and NSF/DMR 1507713 (RJC).