J. Phys. Chem. C 2008, 112, 6653-6661
6653
Molecular Dynamics Study on Nanoparticle Diffusion in Polymer Melts: A Test of the Stokes-Einstein Law Jun Liu,†,‡ Dapeng Cao,*,† and Liqun Zhang*,‡ DiVision of Molecular and Materials Simulation, Key Lab for Nanomaterials, Ministry of Education, and Key Laboratory of Beijing City on Preparation and Processing of NoVel Polymer Materials, Beijing UniVersity of Chemical Technology, Beijing 100029, P. R. China ReceiVed: January 17, 2008; In Final Form: February 27, 2008
Molecular dynamics simulations are used to investigate the diffusion process of nanoparticles in polymer melts. The effects of size, concentration, and mass of the particle, chain length, and polymer-particle interaction on the diffusion of particles in polymer melts are also explored. Our simulated results indicate that the gyration radius of the polymer chain is the key factor determining the validity of the Stokes-Eintsein(SE) relation in describing the particle diffusion at infinite dilution. When the particle size is larger than the gyration radius of the polymer chain, the SE formula can predict the particle diffusion in the polymer melts correctly, and the particle diffusion does not show mass dependence. When the particle size is smaller than the gyration radius of the polymer chain, however, the SE equation becomes invalid in the prediction of particle diffusion, because particle diffusion is exactly related to the nanoviscosity rather than the macroviscosity used in the SE formula. Furthermore, it is also found that the diffusion coefficient of the particle is inversely proportional to the cube of its hydrodynamic radius. In this regime where the particle size is smaller than the gyration radius of the polymer, particle diffusion is independent of the chain length or molecular weight of the polymer, but dependent on the particle mass. In addition, we also observe the transition process of the particle experiencing macroviscosity to nanoviscosity of the polymer melts by gradually increasing the chain length. The concentration dependence of the particle diffusion is similar to the results from Heyes et al., and at high volume fractions, the negative deviations from the SE formula are found. By exploring the effect of particlepolymer interaction on the diffusion of the particle larger than the gyration radius of the polymer chain, it is found that the condition where the SE formula becomes valid in the prediction of the particle diffusion is that where different reasonably defined effective hydrodynamic radii must be used.
1. Introduction It is known that the dynamics and rheology of colloidal suspensions have been investigated well by using the classic continuum Stokes-Einstein(SE) formula.1 For a large and massive solute molecule of radius R in a solvent of much smaller and lighter molecules, the diffusion coefficient D of the solute is described as follows.2
D)
kBT fπηR
(1)
where η is the pure solvent viscosity, T is the absolute temperature, kB is the Boltzmann constant, and f is a constant determined by the choice of stick (f ) 6) or slip (f ) 4) hydrodynamic boundary conditions at the solute surface. However, in the case of the suspensions of particles in polymer melts and solutions, the validity of eq 1 becomes questionable, and it is affected at least by three length scales, including the size of the nanoparticle, the gyration radius (Rg) of the polymers, and the correlation length (ξ) for the polymer (for a polymer melt, ξ ) a, where a is the size of the monomer/polymer segment).3 In particular, as the particle becomes smaller than various characteristic lengths of the polymer chain, the extremely * Corresponding author.
[email protected] or zhanglq@mail. buct.edu.cn. † Division of Molecular and Materials Simulation. ‡ Key Laboratory of Beijing City on Preparation and Processing of Novel Polymer Materials.
great deviation from simple hydrodynamic motion appears. For example, Mackay et al.4 recently reported that the nanoparticles diffuse approximately 200 times faster in a polymer liquid than the prediction from the SE relation, possibly attributed to the nanoparticles being smaller than the entanglement mesh. Meanwhile, Narayanan et al.5 used X-ray photon correlation spectroscopy in conjunction with resonance-enhanced grazingincidence small-angle X-ray scattering to probe the particle dynamics in thin films, and also found that the particle dynamics differ from the SE Brownian motion, which is caused by the viscoelastic effects and interparticle interactions. Moreover, it is reported in the literature that the viscosity experienced by a large object, i.e., the macroviscosity, is different from the viscosity experienced by a small object, which is particular in the context of polymer melts and solutions.6-11 This point of view is firmly established in the rouse model, which assumes that the macroviscosity differs from the “segmental friction”, i.e., the nanoviscosity.12 In addition, Somoza et al.11 experimentally studied the anthracene rotation in poly(dimethylsiloxane) and poly(isobutylene) by gradually increasing the chain length, and found that the diffusivity of the particles exhibits a sharp transition with the increase of Rg, becoming just dependent on the nanoviscosity rather than the macroviscosity for small R/Rg ratios (R is the particle size). Accordingly, understanding nanoparticle diffusion in the polymer matrix is of fundamental importance. For instance, the preparation of the functional nanocomposites, such as the assembly of inorganic nanoparticles
10.1021/jp800474t CCC: $40.75 © 2008 American Chemical Society Published on Web 04/05/2008
6654 J. Phys. Chem. C, Vol. 112, No. 17, 2008 in polymer matrices and diblock copolymers, involves the dynamics of the particles.13-16 The “self-healing” nanoparticles are now being considered as a potential material based on the diffusion of the nanoparticles to the various defect points in the nanocomposites to retain their poperties.17-19 Furthermore, the mobility of the nanoparticles has a significant effect on the performance of the toughness of the nanocomposites.20-22 For example, by using molecular dynamics (MD) simulation, Gersappe20 found that the higher the particle mobility, the better the toughness of the nanocomposites. They also pointed out that particle mobility is a complex function of the size of the filler, the attraction between the polymer and the filler, and the thermodynamic state of the matrix. Besides, Gersappe et al.23 also observed that the dewetting dynamics of nanofiller-polymer films on a substrate is related to the particle mobility. More importantly, at present, many experimental techniques have moved toward regimes on the smaller length scales such as microrheology,24-30 which uses the generalized SE relation to assess the shear modulus and viscosity. Then, the accurate description of nanoscale dynamics becomes essential for the interpretation of experimental data. Besides the recent findings mentioned above,4,5 deviations from the SE law have also been observed elsewhere. Liu and his co-workers31 reported the invalidity of the SE relation for particle size comparable to the polymer segments in solution. However, Vergeles and co-workers32 found that the SE relation can serve as a qualitative approximation in the prediction of the particle diffusion when the particle size becomes comparable to the that of the solvent. Recently, extensive investigations into the breakdown of the SE law in some other cases have also been published in the literature.33-40 For instance, Kumar et al.35 used MD simulation to study hard-sphere fluids of high density, and found a clear breakdown of the Stokes-Einstein (SE) equation, attributed to the hopping particles in analogy with the unusual diffusive behavior of supercooled liquids. As for the simple fluids, the diffusion of the solute particle in the molecular liquid has been widely studied based on the SE formula by employing MD simulations.41-50 Moreover, some reviews have also been reported about the particle dynamics in viscoelastic media, such as the polymer solutions.30,51 However, the diffusion of the particle immersed in the polymer melts is less well investigated. Therefore, in light of this significant issue and the possible occurrence of the potentially unusual effects illustrated above, we use the MD method to investigate the diffusion of particles in the polymer melts. The emphasis is mainly placed on the extensive testing of the applicability of the classic continuum SE formula in describing the particle dynamics in polymer melts. Moreover, some comparisons are also made to recent experiments. To the best of our knowledge, the validity of the SE formula in describing particle diffusion in the polymer matrices is not systematically elucidated yet. The rest of this work is organized as follows. First, the models and simulation details are depicted. Then, the effects of the size, concentration, and mass of the particle, the chain length of the polymer and the polymer-particle interaction are explored. Finally, the important conclusions are drawn, and the discussion is addressed. 2. Models and Simulation Methods 2.1. Models. In our molecular dynamics (MD) simulation, we use the standard bead-spring model proposed by Kremer and Grest52 to represent the polymeric chain. The interactions between all pairs of monomers are described by a purely
Liu et al. repulsive potential derived by truncating and shifting LennardJones (TSLJ) potential at its minimum, given by
U(r) )
{
4
[(σr ) - (σr ) + 41] r < 2 12
6
1/6
σ
(2)
r g 21/6σ
0
The interaction between the adjacent bonded monomers is represented by a stiff finite extensible nonlinear elastic (FENE) potential
[ ( )]
VFENE ) -0.5kR02 ln 1 -
r R0
2
(3)
where k ) 30 /σ2 and R0 ) 1.5σ. The combined effects of the TSLJ and FENE potentials prevent chains from crossing each other. Since we do not intend to study a specific polymer, we use reduced units, in which , m, and σ are assumed to be unit ( is the LJ energy parameter, m and σ are the mass and diameter of the monomer respectively). The nanoparticles are modeled as Lennard-Jones spheres of radius Rn. The interactions between polymer segments and smooth nanoparticles, and nanoparticles and nanoparticles are modeled using modified Lennard-Jones functions, which account for the excluded volume of the beads and particles by offsetting the interaction range by REV53 as follows:
Unp(r) )
{
4np ∞
[(
σ r - REV
12
Unn(r) )
{
4nn ∞
[(
) ]
) (
σ r - REV
-
6 σ 1 + r > REV ) Rn - σ/2 r - REV 4 r e REV ) Rn - σ/2 (4a)
) ( 12
-
) ]
6 1 σ + r > REV ) 2Rn - σ r - REV 4 r e REV ) 2Rn - σ (4b)
where Unp and Unn denote the polymer-nanoparticle and nanoparticle-nanoparticle interactions, respectively. Rn is the radius of nanoparticle and r is the separation distance between interacting sites. Both the interaction parameters np and nn are equal to 1.0, unless otherwise noted. The polymer-particle interaction in the system is truncated and shifted at the separation r ) REV + 2.5σ, while the interaction between particles is truncated and shifted at r ) REV + 21/6σ. In this work, the reduced units of all quantities are adopted, and defined as below:
T* ) kBT/ t* ) t/τ; τ ) σxm/
(5)
F* ) Fσ3 η* ) η/kBTσ-3τ where T* is the reduced temperature, and t* is the reduced time, and η* is the reduced viscosity. All MD runs are carried out by using the large-scale atomic/molecular massively parallel simulator (LAMMPS) that was developed by Sandia National Laboratories.54 In particular, LAMMPS has potentials for soft materials (biomolecules, polymers), solid-state materials (metals, semiconductors), and coarse-grained systems. It can be used to model atoms or, more generically, as a parallel particle simulator at the mesoscale or continuum levels. 2.2. Methods. In our MD simulations, the NVT ensemble is used, where the temperature is fixed at T* ) 1.0 by using the Nose-Hoover thermostat, and the reduced number density of
Nanoparticle Diffusion in Polymer Melts
J. Phys. Chem. C, Vol. 112, No. 17, 2008 6655
TABLE 1: Diffusion Coefficient D of the Particle as a Function of the Size Ratio Rn/Rga Dn Rn/Rg 1.0 2.0 3.0 4.0 6.0 8.0 9.0
Mn
Nn
D (MD)
0.125 1.0 3 0.03522 ( 0.00162 0.25 8.0 3 0.008417 ( 0.0000398 0.375 27.0 3 0.003845 ( 0.000452 0.5 64.0 3 0.001893 ( 0.000105 0.75 216.0 1 0.0009133 ( 0.000063 1.0 512.0 1 0.0005525 ( 0.000057 1.125 729.0 1 0.000505 ( 0.000045
Dslip (SE) Dstick (SE) 0.001873 0.001249 0.0009366 0.00075 0.0005352 0.0004162 0.0003747
0.001249 0.0008326 0.0006244 0.0005 0.0003568 0.0002775 0.0002498
a Dn and Mn denote the diameter and mass of the nanoparticle, respectively, and Nn is the number of particles considered in each case.
monomer is fixed at F* ) 0.84, which corresponds to that of a dense melt. Periodic boundary conditions are employed in all three directions. For the filled systems, the volume of the simulation box is increased so that the volume fraction of the polymers is the same as that of the neat systems. The velocityVerlet55 algorithm is used to integrate the equations of motion, with a time step δt ) 0.001, where the time is reduced by τ. To generate the initial configurations, we place the polymer chains of stretching conformation in a large box, and obtain a system of low density. Then, the NPT simulation is used to compress the system of low density to the melt density. We equilibrate all structures over time such that each chain has moved at least 2Rg. Such equilibrated structures are then used as starting configurations for production runs of (10-20) × 106 MD steps, over which we collect dynamics data. As for the investigation of the particle diffusion in the dilute limit, in this simulation, we used one nanoparticle (i.e., Nn ) 1) for its size Rn larger than 2.0σ, and three nanoparticles (i.e., Nn ) 3) for Rn smaller than or equal to 2.0σ. The diffusion coefficient D of the nanoaprticles is obtained using the Einstein equation
〈(rcm(t) - rcm(0))2〉 tf∞ 6t
D ) lim
(6)
The particle velocity autocorrelation function (VACF) is given by
b(t)‚V b(0)〉 CV(t) ) 〈V
(7)
3. Results and Discussion 3.1. Size Effect. 3.1.1. Effect of Particle Size. First, we investigate the effect of nanoparticle size on its dynamics in the dilute limit. The diffusion coefficients of the nanoparticles obtained from the MD simulation are listed in Table 1. For clear illustration, the simulated results are shown in Figure 1. It is noted that the mass of the nanoparticle is proportional to its volume. So, the mass density of the nanoparticle remains unchanged. For simplicity, it is assumed that the mass density of the nanoparticle is equal to that of the monomers. In this way, the nanoparticle can be varied in a physically meaningful manner, i.e., gradually approaching to the Brownian limit (the nanoparticle mass Mn f ∞, and so the particle size Rn f ∞). In this simulation, we consider 100 chains of length N ) 60 with its radius of the gyration Rg ) 4.0σ, which is unperturbed by the addition of the nanoparticles. This finding agrees with many reports in the simulation literature,56 although it remains experimentally unverified.57,58 It is noted that the simulation of larger particles than those examined here would be challenging because of the increasing role of finite-size effects. The diffusion coefficients (D) of the nanoparticles are obtained by performing several parallel simulations with dif-
Figure 1. The diffusion coefficient D of the particle as a function of the size ratio R/Rg. The particle mass is proportional to its volume. Open squares represent MD data; full dots represent the SE prediction with slip boundary condition; some statistical errors from MD data are smaller than symbol size.
ferent initial configurations. The calculated error bars are also shown in Figure 1. The main objective of this paper is to investigate the applicability of the SE formula given by eq 1, which can be written in another form as
D)
T* / fπη*σ12
(8)
where the constants f ) 4 and 6 for the slip and stick / hydrodynamic boundary conditions, respectively, and σ12 is the / cross radius of the monomer-nanoparticle, σ12 ) r′ + r ) Rn + σ/2,43 where r′ and r are the radii of monomer and / is considered as the nanoparticle, respectively. The fact that σ12 hydrodynamic radius can be justified intuitively by noting that r′ + r is the contact distance between the particle and the monomers, and it would be invalid to apply the hydrodynamic boundary conditions at a smaller distance than the contact one in the case of the attractive interaction between the polymer and the particles. As for the polymer melt with monomer number density equal to 0.84, we can get the reliable viscosity value in the literature, η* ≈ 42.5,59 in which both the polymer model and the simulated number density of the monomer are the same as these used here. Although the total number of simulated monomers is a little different (in the ref 59, the number is 8400, while it is 6000 in this work for different chain lengths), the system size dependence of the zero-shear viscosity is relatively weak for systems in excess of 500 particles.60 On the basis of the viscosity value, we can use the SE formula in eq 8 to calculate the diffusion coefficient of the nanoparticles in the polymer melts. For comparison with the MD simulation, the diffusion coefficients from the SE formula are also added into Figure 1. It is found that, with the increase of the size ratio of R/Rg, the SE diffusion coefficient gradually approximates the MD data under the slip (dotted curve) boundary condition. As the size ratio of R/Rg increases to 1, the SE prediction and MD diffusion coefficient are the same basically. Obviously, the smaller the size ratio of R/Rg, the larger is the deviation between the SE and MD diffusion coefficients. Actually, in the small size ratio, for example, R/Rg ) 0.25, the SE formula becomes invalid and cannot describe the real diffusion of nanoparticles in the polymer melts, because the diffusion coefficient from the SE formula is 1 order of magnitude smaller than that from the MD simulation. This observation qualitatively agrees with the experimental results,4 where the size ratio of R/Rg is also 0.25 with 200 times faster diffusion compared with the SE equation. However, in the experiment, the polymer chain is
6656 J. Phys. Chem. C, Vol. 112, No. 17, 2008
Figure 2. The plot of ln D vs ln σ12, where D is the diffusion coefficient of the particle and σ12 is the hydrodynamics radius or the cross radius of the monomer particle. The slop of the fitted line is about -3, which suggests that the diffusion coefficient of the particle is inversely proportional to the cube of the hydrodynamic radius in the range of the size ratio R/Rg < 1, in contrast with the continuum SE prediction that the diffusion coefficient is inversely proportional to the hydrodynamic radius of the particle.
entangled in the reptation regime, while the polymer length used here is smaller than the entanglement length of the polymer chain.20 This simulated result is also analogous to the enhancement of the diffusion of small solutes in simple dense liquids, compared to the SE value.61 The reason may be that, in the SE formula, the diffusion coefficient is obtained on the basis of the macroscopic viscosity of the pure polymer melt. In fact, however, the diffusion coefficient of the nanoparticle is closely related to the microscopic viscosity of the polymer melt. That is to say, the local viscosity experienced by nanoparticles is much smaller than the macroscopic viscosity, as speculated by Wyart and de Gennes6 and other researchers,7-10 which therefore leads to underestimation of the diffusion coefficient of the nanoparticle in polymer melts by the SE equation related to the macroscopic viscosity in the case of the small size ratio of R/Rg. We believe that the mechanism for the small nanoparticle to sense the nanoviscosity is that the small nanoparticle can move through the polymer melts without necessarily waiting for chains to relax their conformations, which is coupled with the macroviscosity of the polymer melts. With the increase of the size ratio of R/Rg, the chainlike solvent becomes a continuum on the length scale of the chain size (i.e., the radius of the gyration Rg). Certainly, the big particles experience the macroviscosity of polymer melts and return to Brownian behavior. As a result, in the case of R/Rg ) 1, the diffusion coefficient from the MD simulation is basically the same as that from the macroviscosity-dependent SE formula under the slip boundary condition, as shown in Figure 1. In order to further investigate the regime in which the particle size is smaller than the gyration radius of the chains, we present in Figure 2 the diffusion coefficient of the particle as a function of the hydrodynamic radius of the particle or the cross radius of the particle-monomer, σ12. In the logarithm coordinates, the diffusion coefficient of the particle and the hydrodynamic radius are related by a linear function with the slop of the line equal to -3.0. It indicates that the diffusion coefficient of the particle is inversely proportional to the cube of its hydrodynamic radius in the case of the particle size smaller than the gyration radius of polymer. This is quite different from the SE prediction, in which the diffusion coefficient of the particle is inversely proportional to the hydrodynamic radius. This finding is in good agreement with the point of view from ref 4. For small nanoparticles, the friction between the polymer and the particle is caused by the monomer rubbing the nanoparticle surface. Therefore, the resulting friction is proportional to the particle
Liu et al.
Figure 3. The velocity autocorrelation function Cv(t) of the particle with respect to time t. The dashed line stands for R/Rg ) 0.125 and the solid line for R/Rg ) 0.25. τ denotes the LJ time (τ ) σxm/).
TABLE 2: Gyration Radius of the Polymer Chains of Different Lengths, Where the Number Density of the Monomers Is Equal to 0.84 chain length N
Rg/σ
10 20 30 40 50 60 80 120 200
1.487 ( 0.001 2.157 ( 0.003 2.753 ( 0.004 3.226 ( 0.004 3.542 ( 0.005 4.028 ( 0.003 4.682 ( 0.007 5.925 ( 0.021 6.465 ( 0.002
surface. The local viscosity (i.e., nanoviscosity) attributed to this friction is supposed to scale as R2 (R is the hydrodynamic radius). Applying the local viscosity rather than the macrovisosity into the SE formula, we can get the result illustrated in Figure 2. Next, we investigate the effect of particle size on the velocity autocorrelation function (VACF), CV(t). Two representative size ratios are shown in Figure 3. It can be observed that the corresponding VACF are strongly size dependent, which is consistent with the great size-dependent diffusion. As the particle size or equivalently its mass decreases, the VACF becomes more oscillatory, which indicates a stronger “backscattering” of the particle as it “rattles” in the chainlike solvent cage. This observation qualitatively agrees with the simulated results in the case of solute diffusion in the dense simple liquids.45 3.1.2. Effect of Chain Length. In this section, we consider the effect of the chain length of the polymer on the diffusion of particles. Therefore, we keep the particle size fixed, and vary the chain length of the polymer from N ) 10 to 200. In particular, the choice of the chain length N ) 200 is to address the effect of entanglement on the particle dynamics.62 In the following simulations, we use a total of 6000 monomers embedded in the periodic simulation box for all cases. First, we simulate the pure system only containing polymer chains. The simulated gyration radius of different chain lengths at the number density of monomer equal to 0.84 is listed in Table 2. For clear illustration, logarithm plots of the mean-squared gyration radius and the end-to-end distance of the polymers versus the chain length are shown in Figure 4. The simulated results are in excellent agreement with the data from Kumar et al.63 and qualitatively agree with the results from Kremer and Grest.52 It should be pointed out that the number density of their systems is F* ) 0.85, while the density of our system is F* ) 0.84. In the following simulation, we explore the case of the addition of particles into the polymer melts. In the case where
Nanoparticle Diffusion in Polymer Melts
Figure 4. Logarithm plot of the mean-square gyration radius 〈Rg2〉, and the end-to-end distance 〈Reed2〉 of the chains as a function of the chain length (N). The full square is from ref 63.
Figure 5. (a) The diffusion coefficient (D) of the particle with respect to the different chain length (N), where the number density is equal to 0.84 and the particle radius is set to 1.0σ, which is smaller than the gyration radius (Rg) of all the simulated chains (also see Table 2). The particle mass is proportional to its volume. The simulated results reveal that the particle diffusion (D) is independent of the chain length (N) within statistical errors. The dot horizontal line is a best fit for the simulated results. (b) The velocity autocorrelation function (VACF) of the particle for polymer chains of different lengths. The results shown in this figure indirectly suggest that the chain length of the polymer has no significant effect on the diffusion coefficient of the particle, which is related to the VACF by the Green-Kubo integral formula.
the radii of the particles are fixed at 1.0σ and 3.0σ, respectively, the corresponding numbers of particles in the simulation box are 3 and 1 in the dilute limit. As for the particles with radius equal to 1.0σ, which is smaller than the gyration radius of all the simulated chains (also see Table 2), the diffusion coefficients of the particle with respect to different chain lengths is shown in Figure 5a. It can be found from the figure that the diffusion coefficient of the particle shows no polymer weight or polymer chain length dependence within the statistical uncertainty at a fixed number density for different chain lengths. This reason
J. Phys. Chem. C, Vol. 112, No. 17, 2008 6657
Figure 6. The diffusion coefficient of the particle with respect to different chain length (N), where the number density is equal to 0.84 and the particle radius is set to 3.0σ. Open squares represent MD data; full dots represent the SE prediction with stick boundary condition for chains only at the rouse regime. Some statistical errors from MD data are smaller than symbol size. The great deviations from the SE formula are observed with increasing chain length.
may be the same as that obtained from Figure 1, where for particle size smaller than the gyration radius of the polymer, the particles experience local viscosity related to the monomeric friction coefficient and the surface area of the particle. However, for the chains in the rouse regime simulated here, it is wellknown that the monomeric friction coefficient is independent of the molecular weight.64 Therefore, as for the small particle, we observe that the molecular weight of the polymer has no significant influence on the particle diffusion. In fact, it is also well-established experimentally that the diffusivity of a small probe molecule in polymer solutions is independent of the molecular weight of the polymer matrix at the Rouse regime.65-67 This simulated result also qualitatively agrees well with the experimental result where the particle rotation becomes independent of the chain length at high molecular weight.10,11 Meanwhile, the VACF of the particle for several representative polymer chain lengths is also shown in Figure 5b. Interestingly, in this case, the simulated result reveals that the VACF of the particle is also independent of the chain length, even with the same initial values at t ) 0 in the figure. This also indirectly validates the simulated results in Figure 5a, because the diffusion coefficient is directly related to the VACF according to the Green-Kubo formula.68 When the radius of the particle increases to 3.0σ, the simulated results are shown in Figure 6. It should be pointed out that the simulated chain lengths from 10 to 80 are in the Rouse regime, and the macroviscosity values of the polymer melts, which are used to predict the diffusion coefficients of the particles according to the SE formula, are proportional to the molecular weight or chain length of the polymer.63 For the chains containing 10 and 20 monomers, respectively, we found that the simulated data agrees well with the SE prediction with the stick boundary condition within statistical error, because in this case, the particle radius is larger than the gyration radius of the chain. When the chain length increases from N ) 20, the SE equation starts to exhibit underestimation of the diffusion coefficient of the particle, compared to the MD simulation, because the gyration radius of the polymer becomes comparable to the particle size in this case (also see Table 2). For the chain lengths from N ) 20 to 50, the decrease amplitude of the diffusion coefficient from MD simulations is smaller than that from the SE prediction. Since it is well-established that the molecular weight dependence of the particle disappears at the high molecular weight of polymer melts, one reasonable explanation is that, in the regime, the motion of the entire chain
6658 J. Phys. Chem. C, Vol. 112, No. 17, 2008 is still correlated with the particle motion to some extent and the weak but finite molecular weight dependence still exists. With the further increase of the chain length from N ) 50, accompanying the fact that the particle size becomes much smaller than the gyration radius of the chain, the particle motion becomes decoupled with the motion of the entire polymer chain, only being dependent on the local viscosity of the polymer melts, i.e., the particle diffusion becomes independent of the chain length and the diffusion coefficient of the particle remains nearly unchanged in this case. This point is coincident with the result concluded from Figure 5a. The whole process from the gradual decrease of the particle diffusion to becoming polymer molecular weight independent is also rationalized by the experimental observation11 in which the nanoviscosity experienced by the anthracene increases slightly before reaching an asymptotic value with the increase of the chain length. Here, it is pointed out that the nanoviscosity is equal to the corresponding macroviscosity for short chains in the experiment. In addition, it should be noted that, in the case of the chain length N ) 200, which is in the entanglement regime,62 we can find from the simulated result in Figure 6 that the entanglement of the polymer chains has no significant effect on the particle diffusion. This analysis is manifested as well in Figure 5a, and this conclusion is also observed in the experimental results from Somoza et al.,11 in which it is claimed that the chain entanglement has no effect on the nanoviscosity. 3.2. Effect of Finite Particle Concentration. It is known that the application of the SE formula is initially linking the viscosity and diffusion at infinite dilution, while the concentration dependence of the particle diffusion is less well understood.69-71 Here, we explore the effect of particle concentration on its diffusion in the polymer melts. Under the same simulated conditions, in comparison with the dilute limit case with several particles before, we plot the diffusion coefficient of the particle versus the volume fraction of the particle in Figure 7. It is noted that, for convenience, we only simulate the case in which the radius of the particle is 1.0σ and the length of the polymer is 60. The number of particles is increased from 32 to 640, corresponding to the volume fraction of the particle from φ ) 0.0184 to 0.2728. We can see from Figure 7a that the diffusion coefficient of the particle decreases monotonically with the increase of the volume fraction of the particle. This is attributed to interparticle interactions related to the polymer chain-mediated hydrodynamic interactions (HI) because of the formation of chain bridges between neighboring particles with increasing volume fraction of the particle. The bridge networks develop gradually with the increase of the particle loadings, and it accordingly retards the particle diffusion. The dotted line in Figure 7a is the SE prediction value with the slip boundary condition at infinite dilution. As illustrated previously, in this case the SE formula becomes invalid, because the particle diffusion is only related to the local viscosity. The simulated results indicate that, with the network structure of chain bridges strengthening, the particle diffusion decreased monotonically. Furthermore, at high volume fraction, i.e., φ ) 0.2728, the diffusion coefficient of the particle even crosses over the prediction value from the SE formula at infinite dilution, and is only 1/4 of the SE prediction value. In this situation, it is claimed that the SE formula overestimates the particle diffusion at high volume fraction. Our results to some extent also confirm the experimental observation from Cole and co-workers,72 where it was demonstrated that the particle mobility was decreased 2 to 3 orders of magnitude compared with the predictions by the SE theory. It may be attributed to the formation of bridges by
Liu et al.
Figure 7. (a) The plot of the diffusion coefficient (D) of the particle with respect to the volume fraction of the particle (φ), where the particle radius is fixed at 1.0σ and the chain length of the polymer is N ) 60. The open squares represent the MD data, while the full line represents a polynomial fit including linear and quadratic terms. The dotted line stands for the SE prediction at infinite dilute. The figure reveals that the particle diffusion decreases monotonically with increasing the volume fraction of the particle. (b) The average number of particles adsorbed on each chain (An) as a function of the volume fraction of the particle φ. The solid line stands for a linear fit. This indicates that the number of bridge segments per chain is increased with the increase of the particle loadings.
polymer chains adsorbed on different particles due to the polymer-particle attraction. Additionally, we also found that the concentration dependence of the particle diffusion can be well-fitted by a polynomial function including linear and quadratic terms. By extrapolation, the polynomial function gives a value of 0.0081 for D at infinite dilution, which is in good agreement with our simulated data for the particle size with Rn ) 1.0, as listed in Table 1. This quantitative analysis agrees with the concentration dependence of the solute diffusion at high volume fractions in the case of dense simple liquids.70 In order to quantify the development of the network structure of chain bridges by increasing the particle loadings, we also plot the average amount of particles adsorbed on each chain represented by An versus the volume fraction φ of the particle in Figure 7b. It is noted that a polymer chain is considered to be adsorbed on the particle with at least a polymer segment located in its interface shells, which is defined as the shell of width 2σ surrounding the particle surface, i.e., within the cutoff distance of the polymer-particle interaction. It can be observed from the figure that the average number of particles adsorbed on each chain increases linearly with the volume fraction of the particle in the simulated case. Simultaneously, we also calculate the number of chains adsorbed on each particle, and find that the number of the adsorbed chains remains nearly unaltered and is equal to approximately 12 with the increase of the volume fraction of the particle. On the basis of the fixed number of simulated chains, this result just validates the linear
Nanoparticle Diffusion in Polymer Melts
Figure 8. The particle diffusion coefficient (D) as a function of the particle mass Mn. The chain length is set to 10, and the corresponding gyration radius is Rg ) 1.487σ.
relationship between the number of particles adsorbed on each chain and the volume fraction of the particle, as seen in Figure 7b. For the fixed chain length, the length of the bridge segments between neighboring particles shortens with the increase of the adsorption amount of particles on each chain. It accordingly results in stronger chain-mediated bridging interactions between particles. The interaction retards the particle diffusion, which supports the observation in Figure 7a. 3.3. Mass Dependence. After investigating the size and concentration dependence of the particle diffusion, in this section we are motivated to study the mass dependence of particle diffusion, since the ideas embodied in the SE equation imply that the mass of the particle should be an irrelevant variable. For simplicity, we choose the polymer chain with 10 monomers and two types of particles with radii of 0.5σ and 2.0σ. The number of particles is fixed at 3. For the two particle sizes, we investigate several particle masses with Mn equal to 1.0, 8.0, 64.0, 100.0, 150.0, and 216.0, respectively. Figure 8 shows the mass dependence of the particle diffusion in polymer melts for the two particle sizes. It can be observed that, for particle size Rn ) 2.0σ, which is greater than the gyration radius of the polymer, the particle diffusion exhibits mass-independent behavior, and the diffusion coefficient of the particle approximates the value 0.004496, which is predicted by the SE formula using the slip boundary condition. This suggests that, in this case, particle diffusion obeys SE behavior, as mentioned above, and also shows no mass dependence. However, for particle size Rn ) 0.5, which is smaller than the gyration radius of the polymer chain, the particle diffusion does not exhibit Brownian behavior and violates the SE law, exhibiting a decrease of the diffusion coefficient with increasing particle mass. However, for high particle mass, e.g., 100.0, the particle diffusion becomes independent of the particle mass, and the value is still larger than the SE prediction of 0.01124 with the slip boundary condition. The behavior in which the particle diffusion becomes mass-independent at heavier mass for the fixed particle size is similar to the diffusion of the small solutes in simple liquids.43 This can be explained as follows. For particle size Rn ) 0.5σ, which is smaller than the gyration radius of the polymer chain Rg ≈ 1.48σ, for small particle mass the time scale associated with particle diffusion is smaller than the relaxation time scale of polymer chains, contrary to the underlying assumption embodied in the SE formula. Therefore, particle diffusion presents mass dependence. When the particle mass becomes heavy enough, the relaxation time of the particle becomes larger than that of the polymer chains, and the particle diffusion no longer shows mass dependence, while the existing deviation from the SE formula is attributed to the fact that the
J. Phys. Chem. C, Vol. 112, No. 17, 2008 6659
Figure 9. Plot of the diffusion coefficient of the particle vs the polymer-particle interaction, where the chain length is 10 and the particle size is fixed at 3.0σ. The dash horizontal line represents the SE prediction with the effective hydrodynamic radius (EHR) equal to the bare particle radius, the dot horizontal line stands for the prediction with EHR equal to the cross radius of the monomer particle, and the dash-dot horizontal line means the prediction with EHR equal to Rn + 2Rg.
particle experiences the local viscosity of the polymer melts. This explanation suggests that the mass and size of the particle to some extent are coupled with each other to influence its diffusion, but the dominant factor for particle diffusion is the particle size. The observation in Figure 8 also agrees well with the results from Kumar et al.73 They found that the particle diffusion exhibits mass dependence, and is larger for the lighter particle for the fixed particle size in the simulated range. But, it should be noted that their simulation was performed in the regime where the particle radius is smaller than the gyration radius of the polymer chain, i.e., the radius of the nanoparticle is 2.5σ while the gyration radius of the chains is 4.7σ. Therefore, their simulated situation is similar to the case of our simulations, and these simulation results qualitatively agree well with each other. 3.4. Polymer-Particle Interaction Effect. For the SE formula in eq 1, the effect of the solute-solvent interaction on the particle diffusion is usually incorporated within the constant f and the effective hydrodynamic radius R. Here, we investigate the effect of the polymer-particle interaction on particle diffusion by adjusting the interaction parameters np from 0.1 to 8.0. For simplicity, we choose the polymer chain containing 10 monomers, and the particle of radius Rn ) 3.0σ with only one particle. In this case, the particle diffusion experiences macroviscosity and the SE equation is valid for predicting the particle diffusion, as mentioned earlier. All the other simulated conditions are the same as before except the interaction parameter of the polymer-particle. It should be noted that it is assumed that the change in the polymer-particle interaction has no significant effect on the configuration of the polymer chain in the dilute limit. The simulated diffusion coefficients are shown in Figure 9. Apparently, it can be found that the particle diffusion decreases monotonically with the polymerparticle interaction increasing. For large np, the particle diffusion becomes unchanged basically. In Figure 9, we also present three different horizontal lines predicted by the SE formula at differently effective hydrodynamic radius (EHR). Here, we choose the stick boundary condition, because the particle size is larger than the gyration radius of the polymer.46 As for the dashed horizontal line predicted by using the bare particle radius as the EHR, the SE predicted value approximates the simulated diffusion coefficient at weak interaction, i.e., np ) 0.1, while for the dotted horizontal line predicted by using the cross radius of the particle-monomer as the EHR, the SE predicted value
6660 J. Phys. Chem. C, Vol. 112, No. 17, 2008
Liu et al.
agrees well with the particle diffusion at interaction parameter np equal to 1.0. This also justifies the validity of our choice of cross radius of the particle-monomer as the hydrodynamic radius with np ) 1.0 in the simulation before. For the case of the strong particle-polymer interaction, the nearest-neighbor chainlike solvent molecules adhere to the particle, which leads to the simultaneous movement of the particle and polymer. Therefore, the resulting hydrodynamic radius Rh of this chainlike solvent-berg should be Rh ) R + 2Rg. By using Rh as the EHR, the prediction by the SE formula is also shown in Figure 9, denoted by the dash-dot horizontal line. This predicted value just agrees with the MD diffusion coefficient in the limit of strong interaction strength, i.e., np ) 8.0 within statistical error. In general, for interactions of different strengths, the condition that the SE relation validates is that the different reasonably defined hydrodynamic radii must be used.
significant theoretical information on the diffusion process of the particle in dense polymer melts or polymer solutions, and of course shed light on the regimes of applicability of experimental techniques such as microrheology24-30 and the preparation of advanced nanocomposites possessing high toughness.20-22
4. Summary and Conclusions
References and Notes
In this work, we used MD simulations to fully study the diffusion of the nanoparticle in the polymer melts, and to explore the effects of the size, concentration, and mass of the particle, chain length, and polymer-particle interaction on the particle diffusion. Furthermore, the validity of the Stokes-Einstein law in describing particle diffusion in the polymer melts was also tested systematically. The results from MD simulations indicated that, when the size of the particle is larger than the gyration radius of the polymer chain, the SE formula can predict well the diffusion of the particle in the polymer melts under the slip or stick boundary conditions. Moreover, the particle diffusion shows no mass dependence behavior. However, when the size of the particle is smaller than the gyration radius of the polymer chain, the SE formula becomes invalid in the prediction of the diffusion coefficient of the particle. It is because the local viscosity experienced by particles is smaller than the corresponding macroviscosity. Moreover, it was also found that the diffusion coefficient of the particle is inversely proportional to the cube of the particle size, and becomes independent of the molecular weight or chain length of the polymer. In addition, we also observed the transition process for the particle to experience from the macroviscosity related to the motion of the entire chain to the nanoviscosity relevant to the monomeric friction coefficient and the particle surface. In this regime of R < Rg, the particle diffusion also exhibits mass-dependent behavior and the particle diffusion coefficient approaches a constant at heavier mass. The velocity autocorrelation function (VACF) reflects different motion behavior at the microscopic scale for different particle sizes. By investigating the concentration dependence of the particle diffusion, we found that, at high particle loadings, the SE formula overestimates the diffusion coefficient of the particle, which is coincident with the corresponding experimental research. By exploring different particle-polymer interactions, it was found that, when the particle size is larger than the gyration radius of the polymer chain, the condition where the SE formula can be applied to describe the diffusion of the particle in the polymer melts is that in which the different reasonably defined effective hydrodynamic radii (EHR) must be used. However, for the finite size effect, we do not intend to simulate the broad range of particle sizes. More experiments and simulation investigations are expected to elucidate this complicated issue: for example, the determination of hydrodynamic boundary conditions in different situations in order to more accurately study the diffusion process based on the SE formula. It is expected that our MD simulation data provide
(1) Batchelor, G. K. An Introduction to Fluid Dynamics; Cambridge University Press: Cambridge, 1973. (2) Zwanzig, R.; Harrison, A. K. J. Chem. Phys. 1985, 83, 5861. (3) Ganesan, V.; Pryamitsyn, V.; Surve, M.; Narayanan, B. J. Chem. Phys. 2006, 124, 221102. (4) Tuteja, A.; Mackay, M. E.; Narayanan, S.; Asokan, S.; Wong, M. S. Nano Lett. 2007, 7, 1276. (5) Narayanan, S.; Lee, D. R.; Hagman, A.; Li, X. F.; Wang, J. Phys. ReV. Lett. 2007, 98, 185506. (6) Wyart, F. B.; de Gennes, P. G. Euro. Phys. J. E 2000, 1, 93. (7) Ye, X.; Tong, P.; Fetters, L. J. Macromolecules 1998, 31, 5785. (8) Cheng, Y.; Prud’homme, R. K.; Thomas, J. L. Macromolecules 2002, 35, 8111. (9) Won, J.; Onyenemezu, C.; Miller, W. G.; Lodge, T. P. Macromolecules 1994, 27, 7389. (10) Sluch, M. I.; Somoza, M. M.; Berg, M. A. J. Phys. Chem. B 2002, 106, 7385. (11) Somoza, M. M.; Sluch, M. I.; Berg, M. A. Macromolecules 2003, 36, 2721. (12) Ferry, J. D. Viscoelastic Properties of Polymers; John Wiley & Sons: New York, 1980. (13) Lopes, W. A.; Jaeger, H. M. Nature 2001, 414, 735. (14) Li, M.; Schnablegger, H.; Mann, S. Nature 1999, 402, 393. (15) Lin, Y.; Boker, A.; He, J. B.; Sill, K.; Xiang, H. Q.; Abetz, C.; Li, X. F.; Wang, J.; Emrick, T.; Long, S.; Wang, Q.; Balazs, A.; Russell, T. P. Nature 2005, 434, 55. (16) Bockstaller, M. R.; Lapetnikov, Y.; Margel, S.; Thomas, E. L. J. Am. Chem. Soc. 2003, 125, 5276. (17) Gupta, S.; Zhang, Q. L.; Emrick, T.; Balazs, A. C.; Russell, T. P. Nat. Mater. 2006, 5, 229. (18) Smith, K. A.; Tyagi, S.; Balazs, A. C. Macromolecules 2005, 38, 10138. (19) Lee, J. Y.; Buxton, G. A.; Balazs, A. C. J. Chem. Phys. 2004, 121, 5531. (20) Gersappe, D. Phys. ReV. Lett. 2002, 89, 058301. (21) Shah, D.; Maiti, P.; Jiang, D. D.; Batt, C. A.; Giannelis, E. P. AdV. Mater. 2005, 17, 525. (22) Zhou, T. H.; Ruan, W. H.; Rong, M. Z.; Zhang, M. Q.; Mai, Y. L. AdV. Mater. 2007, 19, 2667. (23) Luo, H. B.; Gersappe, D. Macromolecules 2004, 37, 5792. (24) Lu, Q.; Solomon, M. J. Phys. ReV. E 2002, 66, 061504. (25) Chen, D. T.; Weeks, E. R.; Crocker, J. C.; Islam, M. F.; Verma, R.; Gruber, J.; Levine, A. J.; Lubensky, T. C.; Yodh, A. G. Phys. ReV. Lett. 2003, 90, 108301. (26) Papagiannopoulos, A.; Waigh, T. A.; Fluerasu, A.; Fernyhough, C.; Madsen, A. J. Phys.: Condens. Matter 2005, 17, L279. (27) Oppong, F. K.; Rubatat, L.; Frisken, B. J.; Bailey, A. E.; de Bruyn, J. R. Phys. ReV. E 2006, 73, 041405. (28) Furukawa, A. J. Chem. Phys. 2004, 121, 9716. (29) Cheng, Z.; Mason, T. G. Phys. ReV. Lett. 2003, 90, 018304. (30) Solomon, M. J.; Lu, Q. Curr. Opin. Colloid Interface Sci. 2001, 6, 430. (31) Liu, J.; Gardel, M. L.; Kroy, K.; Frey, E.; Hoffman, B. D.; Crocker, J. C.; Bausch, A. R.; Weitz, D. A. Phys. ReV. Lett. 2006, 96, 118104. (32) Vergeles, M.; Keblinski, P.; Koplik, J.; Banavar, J. R. Phys. ReV. E 1996, 53, 4852. (33) Berthier, L.; Chandler, D.; Garrahan, J. P. Europhys. Lett. 2005, 69, 320. (34) Dunstan, D. E.; Stokes, J. Macromolecules 2000, 33, 193.
Acknowledgment. This work is supported by the National Natural Science Foundation (NSF) of China (20776005, 20736002), The Outstanding Young Scientists Foundation of NSF of China (50725310), Beijing Novel Program (2006B17), the Program for New Century Excellent Talents (NCET-060095) from Ministry of Education, and “Chemical Grid Project” and Excellent Talents Funding of Beijing University of Chemical Technology. J. Liu is also greatly thankful to Mr. Ran Ni for persistent help and significant discussion.
Nanoparticle Diffusion in Polymer Melts (35) Kumar, S. K.; Szamel, G.; Douglas, J. F. J. Chem. Phys. 2006, 124, 214501. (36) May, H. O.; Mausbach, P. Phys. ReV. E 2007, 76, 031201. (37) Mazza, M. G.; Giovambattista, N.; Stanley, H. E.; Starr, F. W. Phys. ReV. E 2007, 76, 031203. (38) Rajian, J. R.; Quitevis, E. L. J. Chem. Phys. 2007, 126, 224506. (39) Tarjus, G.; Kivelson, D. J. Chem. Phys. 1995, 103, 3071. (40) Zangi, R.; Kaufman, L. J. Phys. ReV. E 2007, 75, 051501. (41) Ouldkaddour, F.; Barrat, J. L. Phys. ReV. A 1992, 45, 2308. (42) Srinivas, G.; Bhattacharyya, S.; Bagchi, B. J. Chem. Phys. 1999, 110, 4477. (43) Ould-Kaddour, F.; Levesque, D. Phys. ReV. E 2000, 63, 011205. (44) Murarka, R. K.; Bhattacharyya, S.; Bagchi, B. J. Chem. Phys. 2002, 117, 10730. (45) Schmidt, J. R.; Skinner, J. L. J. Chem. Phys. 2003, 119, 8062. (46) Schmidt, J. R.; Skinner, J. L. J. Phys. Chem. B 2004, 108, 6767. (47) Sokolovskii, R. O.; Thachuk, M.; Patey, G. N. J. Chem. Phys. 2006, 125, 204502. (48) Sharma, M.; Yashonath, S. J. Phys. Chem. B 2006, 110, 17207. (49) Cappelezzo, M.; Capellari, C. A.; Pezzin, S. H.; Coelho, L. A. F. J. Chem. Phys. 2007, 126, 224516. (50) Ould-Kaddour, F.; Levesque, D. J. Chem. Phys. 2007, 127, 154514. (51) Masaro, L.; Zhu, X. X. Prog. Polym. Sci. 1999, 24, 731. (52) Kremer, K.; Grest, G. S. J. Chem. Phys. 1990, 92, 5057. (53) Smith, J. S.; Bedrov, D.; Smith, G. D. Compos. Sci. Technol. 2003, 63, 1599. (54) Plimpton, S. J. J. Comput. Phys. 1995, 117, 1. (55) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, 1987.
J. Phys. Chem. C, Vol. 112, No. 17, 2008 6661 (56) Vacatello, M. Macromolecules 2001, 34, 1946. (57) Nakatani, A. I.; Chen, W.; Schmidt, R. G.; Gordon, G. V.; Han, C. C. Polymer 2001, 42, 3713. (58) Botti, A.; Pyckhout-Hintzen, W.; Richter, D.; Straube, E.; Urban, V.; Kohlbrecher, J. Phys. B 2000, 276, 371. (59) Kroger, M.; Loose, W.; Hess, S. J. Rheol. 1993, 37, 1057. (60) Heyes, D. M. J. Phys.: Condens. Matter 2007, 19, 376106. (61) Bhattacharyya, S.; Bagchi, B. J. Chem. Phys. 1997, 106, 1757. (62) Picu, R. C.; Rakshit, A. J. Chem. Phys. 2007, 126, 144909. (63) Sen, S.; Kumar, S. K.; Keblinski, P. Macromolecules 2005, 38, 650. (64) de Gennes, P. G. Scaling Concepts in Polymer Physics; Cornell University Press: Ithaca, New York, 1979. (65) Gisser, D. J.; Johnson, B. S.; Ediger, M. D.; Vonmeerwall, E. D. Macromolecules 1993, 26, 512. (66) Vonmeerwall, E. D.; Amis, E. J.; Ferry, J. D. Macromolecules 1985, 18, 260. (67) Landry, M. R.; Gu, Q. J.; Yu, H. Macromolecules 1988, 21, 1158. (68) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids, 2nd ed; Academic: London, 1986. (69) McPhie, M. G.; Daivis, P. J.; Snook, I. K. Phys. ReV. E 2006, 74, 031201. (70) Nuevo, M. J.; Morales, J. J.; Heyes, D. M. Phys. ReV. E 1998, 58, 5845. (71) Snook, I.; O’Malley, B.; McPhie, M.; Daivis, P. J. Mol. Liq. 2003, 103, 405. (72) Cole, D. H.; Shull, K. R.; Rehn, L. E.; Baldo, P. Phys. ReV. Lett. 1997, 78, 5006. (73) Desai, T.; Keblinski, P.; Kumar, S. K. J. Chem. Phys. 2005, 122, 134910.
