Rouse Mode Analysis of Chain Relaxation in Homopolymer Melts.

We use molecular dynamics simulations of the Kremer-Grest (KG) bead-spring model of polymer chains of length between 10 and 500, and a closely related analogue that allows for chain crossing, to clearly delineate the effects of entanglements on the length-scale-dependent chain relaxation in polymer melts. We analyze the resulting trajectories using the Rouse modes of the chains and find that entanglements strongly affect these modes. The relaxation rates of the chains show two limiting effective monomeric frictions, with the local modes experiencing much lower effective friction than the longer modes. The monomeric relaxation rates of longer modes vary approximately inversely with chain length due to kinetic confinement effects. The time-dependent relaxation of Rouse modes has a stretched exponential character with a minimum of stretching exponent in the vicinity of the entanglement chain length. None of these trends are found in models that allow for chain crossing. These facts, in combination, argue for the confined motion of chains for time scales between the entanglement time and their ultimate free diffusion.

Introduction

The structural relaxation of polymer melts has been the subject of considerable expn>erimental and theoretical investigation due to its relevance to the processing of these viscoelastic materials. The relaxation of chains strongly depends on the degree of polymerization (or chain length), N, in a manner that is qualitatively captured by the Rouse–Zimm and reptation models for short and long chains, respectively.[1] The question that logically follows is how the internal relaxations of the chains depend on the length of the subchain under consideration. For long enough chains, it is thought that short subchains relax according to Rouse dynamics and the longest by reptation with an intermediate crossover behavior. The crossover between these two limiting behaviors occurs around the critical molecular weight Nc ∼ Ne – 2Ne, which depends on the details of the polymer local architecture and packing[2,3] and is generally attributed to the topological constraints of entanglement that arise for large N.

Perez-Aparicio et al.[4] studied the dynamics of unentangled PEP melts using neutron spin echo (NSE) spectroscopy and molecular dynamics (MD) simulations. They found that the relaxation of the chains deviates from Rouse behavior at short length scales (i.e., below 2 nm), which are still larger than the Kuhn length of the chains (∼1.1 nm). Analogous deviations from Rouse behavior at time scales shorter than Rouse time were also observed for 1,4-polybutadiene, polyethylene, and poly(ethylene oxide).[5−7] These results are now thought to be due to chain stiffness, other local packing effects, and local frictional variations. For longer chain lengths, Richter et al.[8,9] found that there is a crossover where chain relaxation is significantly slowed down due to the presence of entanglements. Clearly, it is important to understand how topological constraints control the relaxation of different subsections of a chain in this crossover regime. This is the focus of the present work.

To provide context for our studies, we begin with the simplest model for the dynamics of short polymer chains in a homopolymer melt, the Rouse model.[10] A chain is represented by a sequence of beads connected by harmonic springs, each spring with mean-squared bond length b2. This model is known to describe the dynamics of short, unentangled melts reasonably well, though deviations appear at shorter length scales. While these deviations have been thought to be due to local excluded volume interactions and chain stiffness,[11] more recent work has suggested that these effects could also arise from hydrodynamic interactions operating over a range of length scales.[12] For long, entangled chains, the Rouse model describes the dynamics at intermediate time/length scales even though the longer scale dynamics are strongly affected by constraints formed by surrounding chains. The Rouse modes (p = 0, 1, 2, ..., N – 1) of a chain of length N are defined as[13]X⃗ = (2/N)1/2∑r⃗ cos[(pπ/N)(i – 1/2)]. The p = 0 mode describes the motion of the chain center-of-mass, while the modes with 1 ≤ p ≤ N – 1 describe internal relaxations with a mode number p corresponding to a subchain of (N – 1)/p segments. The autocorrelation of the Rouse modes is predicted to be ⟨X⃗(t)·X⃗(0)⟩ = ⟨X⃗2⟩e–; i.e., each mode should decay exponentially and the modes are independent. Here the square of the pth mode amplitude is ⟨X⃗2⟩ = b2/[4 sin2(pπ/2N)], and its reciprocal relaxation time is τ–1 = (12kBT/ζb2) sin2(pπ/2N), where ζ is the monomeric friction coefficient. In contrast to these predictions, the simulation determined autocorrelation function of Rouse modes for chains in a homopolymer melt has been found to be well fit by a stretched exponential relaxation:[14−16] ⟨X⃗(t)·X⃗(0)⟩ = ⟨X⃗2⟩e–(. Here, we note that intramolecular bond correlations, excluded volume effects, and chain stiffness cause ⟨X⃗2⟩ to deviate from the expected p–2 scaling, especially for large p.[17] Such effects along with additional hydrodynamic coupling between different chain sections recently discovered by Farago et al.[12,18] may be the source of these stretched exponential time correlations of the Rouse modes. The effective relaxation time of mode p can be obtained by integrating this relaxation function:[14−16] τeff = ∫0∞e–( dt = (τ/β)Γ(1/β), where Γ(x) is the gamma function. The effective monomeric relaxation rate (monomeric Rouse rate) can be calculated using Weff = 3kBT/ζb2 = 1/[4τeff sin2(pπ/2N)]. For the Rouse model, this quantity should be independent of mode number and only depend on the monomer friction ζ, temperature, and statistical segment length b.

Padding and Briels[15,19] have conducted MD simulations of polyethylene using a coarse-grained model in which one bead represented 20 methylenes. Using this model, they simulated melts as long as 50 beads (or 1000 methylenes) and found that the effective exponent β is dependent on mode number—for the shortest modes (largest p) they found β ∼ 0.8. With decreasing mode number β decreases, reaching a minimum ∼0.5 for modes that are in the vicinity of N/p ∼ Ne/3, independent of chain length. For small p, they find an asymptotic value of β ∼ 0.8. Li et al.[15] found similar results for the value of β but with its minimum value occurring close to N/p ∼ Ne for a bead–spring model with chain length N = 500. Shaffer suggested that this minimum (as well as the presence of stretched exponential relaxations) results from kinetic constraints experienced by long chains.[14] Several questions become apparent from this set of simulations: (i) Are these results general? (ii) Is the minimum of 0.5 in the vicinity of ∼Ne a general result, and how does it correlate with the concept of entanglements?

More recently, Likhtman[20] has examined several aspects relative to the applicability of the Rouse model and a Rouse mode analysis to melts of short, unentangled chains and longer entangled analogues. His results reiterate several important points that are critical to the current work. First, the normalized amplitude of the static Rouse modes, ⟨X2⟩ sin2(pπ/2N), are dependent on mode number, p—this result reflects the power law decay of intramolecular correlations in melts. Second, there is coupling between the dynamics of the different modes, for example the first and the third modes (of a chain of length N), and this coupling only relaxes on the time scale of the slowest mode. While these results thus clearly violate the fundamental principle of mode decoupling in the Rouse model, we shall continue to use this description for the following two reasons. First, as we shall discuss below, experimentalists still tend to model chain dynamics in the language of the Rouse model. Understanding experimental results therefore require us to analyze the simulations in the same manner. Second, our overall goal is to understand the role of nanoparticles in chain melts. Exploring these concepts based on the Rouse model thus seems to be the easiest means, at least in the context of current experimental practice.

Here, we study the internal relaxation of a series of monodisperse polymer melts with different lengths (N = 10–500) using molecular dynamics (MD) simulations. We analyze the results of the simulations using Rouse modes. Results for two representations of polymer chains, namely the Kremer–Grest (KG) bead–spring model and chain-crossing model (in which the chains are unentangled for all N), are compared. Comparison of these two models shows unequivocally that chain uncrossability for long-enough chains gives rise to the minimum in β and also a significant slowing down in the relaxation of the longest chain modes.

Model

Polymer chains are represented by the coarse-grained bead–spring Kremer–Grest (KG) model.[21] Nonbonded monomers interact through the Lennard-Jones (LJ) 12–6 potential: U(r) = 4ε[(σ/r)12 – (σ/r)6] for r ≤ rc, where ε is the LJ energy scale and σ is the monomer diameter. The LJ time scale is τ = (mσ2/ε)1/2, where m is the mass of a monomer. For most of the studies the LJ interaction was cut off at rc = 2.5σ, though some simulations were run with shorter cutoffs for comparison. Two successive segments in a chain are connected by a finitely extensible nonlinear elastic (FENE)[21] potential with k = 30ε/σ2 and R0 = 1.5σ.[21,22] These parameter values ensure noncrossing of the chains. In addition, a three-body bending potential of the form Ubend = kθ(1 + cos θ) is used to control the stiffness of the chains with kθ = 0, 0.75, and 1.5ε in three separate sets of simulations. The entanglement length Ne for these three values of kθ are ≈85, 45, and 28, respectively.[23−25] We have considered chain lengths N = 10–500 to study the relaxation of chains in both unentangled and entangled homopolymer melts.

To evaluate the role of entanglements, we also performed a series of simulations that allow for chain crossing (CC)—these latter class of models have no entanglement effects, even though their static properties (e.g., chain dimensions) closely match their entangled analogues. Following Likhtman[20] and Duering et al.,[26] this is done by “softening” the nonbonded polymer segment pair interactions to Useg–seg(r) = A[1 + cos(πr/rc)], with A = 6.5ε and rc = 1.6σ and switching the bonds from their FENE form to a softer harmonic form Vspring(r) = (k/2)(r – r0)2 with k = 20ε/σ2 and r0 = 1.222σ. These parameter values were chosen so that there is no significant effect on the static properties of the melts while still enabling chain crossing. (A detailed calculation for two bonds crossing at 90° yields a minimum energy of ∼2kBT per bead, corresponding to a bond crossing probability of ∼0.15.)

Most of the simulations are carried out on the NERSC Cray XE6 Hopper using the large-scale atomic/molecular massively parallel simulator (LAMMpan class="Chemical">PS).[27] The initial configurations of the systems are prepared at random at a constant number density while allowing for overlaps among beads. The overlaps are removed by initially using a soft potential between monomers and then by gradually increasing the strength of the potential. After all overlaps are removed, the LJ interactions between monomers are turned on and the volume of the simulation cell is allowed to adjust at a constant pressure pan class="Chemical">P* = 0. Systems of chain length N = 100–500 are equilibrated following the double-bridging algorithm.[22] The shorter N melts are equilibrated by running isobarically and then at constant volume until the chains have moved their own size.

After equilibration, the systems are run at constant volume with a Langevin thermostat with damping constant Γ = 0.1τ–1. All simulations are run at temperature T* = kBT/ε = 1.0 with a time step of 0.01τ. For longer chain lengths we find that the average pressure P* = 0 ± 0.05 ε/σ3, whereas for shorter chains P* = 0 ± 0.1 ε/σ3. We simulated homopolymer melts of M chains of length N for [M, N] as shown in Table 1 for kθ = 0.75. Additional systems with kθ = 0.0, 1.5ε for chain length N = 500 were also simulated to study the effect of bending stiffness on the relaxation.

Table 1

Details of Simulations for kθ = 0.75

			⟨R2g⟩1/2
chain length N	no. of chains M	length of simulation box L/σ	KG model	chain crossing model
10	2000	28.49	1.6	1.5
20	1000	28.37	2.4	2.3
40	500	28.28	3.4	3.4
60	500	32.38	4.4
80	500	35.62	5.1	5.0
100	500	38.37	5.7	5.6
150	500	43.91	6.9
200	500	48.33	8.2	8.0
400	500	60.80	11.4	11.4
500	500	65.56	13.1

Normalized autocorrelation function of different Rouse modes p for a melt of chain length N = 500 with varying chain stiffness: (a) kθ = 0 with rc = 21/6σ, (b) kθ = 0, (c) kθ = 0.75ε, and (d) kθ = 1.5ε with rc = 2.5σ. The dashed black lines are pan class="Disease">fits to stretched expn>onentials.

Results

Chain Stiffness Effects

We first present results for the effect of chain stiffness on the structure and relaxation of entangled chains of length N = 500. For rc = 2.5σ, the monomer number density ρ = 0.89σ–3 for pan class="Chemical">P* = 0, N = 500 and for all three kθ = 0, 0.75, and 1.5. These simulations were run for a total time of 5 × 107τ for kθ = 0 and 0.75ε and 8 × 107τ for kθ = 1.5ε. For comparison, we also ran the system of chains of length N = 500 with a purely repulsive LJ potential (rc = 21/6) for kθ = 0 at the same density.

All autocorrelation functions, i.e., for the different Rouse modes of index p, different stiffness and cutoffs (Figure 1) decay to zero, implying that the chains are fully relaxed. Attraction seems to have a relatively minor effect (Figures 1a vs 1b, also open squares vs closed squares in Figure 2), but the curves for a given mode p are shifted toward longer time with increasing stiffness (from Figures 1b to 1d). The scaled amplitudes of the Rouse modes, ⟨X⃗2⟩ sin2(pπ/2N) = b2/4 representing the mean-squared bond distances of chains, asymptotically approach the characteristic ratio, which increases with chain stiffness (Figure 2a). For small p, ⟨X⃗2⟩ sin2(pπ/2N) = AR[1 – c(N/p)1/2],[28,29] where the constant of proportionality, AR, is found to be linearly proportional to the characteristic ratio, C∞. These results simply reflect the importance of local intramolecular correlations in the chains and that the chains follow Gaussian statistics at long length scales.

Figure 1

Normalized autocorrelation function of different Rouse modes p for a melt of chain length N = 500 with varying chain stiffness: (a) kθ = 0 with rc = 21/6σ, (b) kθ = 0, (c) kθ = 0.75ε, and (d) kθ = 1.5ε with rc = 2.5σ. The dashed black lines are fits to stretched exponentials.

Figure 2

(a) Amplitudes of the autocorrelation function of the Rouse modes for chains of length N = 500. Lines correspond to the scaling relationship ⟨X⃗2⟩ sin2(pπ/2N) = AR(1 – c(N/p)1/2. (b) Effective monomeric relaxation rates of melts of chain length N = 500 with different stiffness. (c) Effective monomeric relaxation rates for chains of different lengths with stiffness kθ = 0.75ε for high N/p scaling. (d) Exponent β from fitting a stretched exponential to the autocorrelation function of the Rouse modes for N = 500.

Our results in Figure 2b for the monomeric relaxation rates are not consistent with the Rouse model prediction that Weff should be independent of N/p. Indeed, only on the shortest length scales (N/p < 20) is the pan class="Chemical">polymer relaxation Rouse-like. The effective monomeric relaxation rates reach a plateau value at large N/p, which we postulate as being due to the crossover to reptation-like scaling arising from constraints to chain motion from the neighboring chains. To lend credence to our postulate for smaller p modes, we note that the relaxation time for the longer chain modes should scale as τ ∼ [N/p]2[1 + (N/Ne)1.4]. From here we can show for small p that the Rouse rates follow Weff ∼ [1 + (N/Ne)1.4]−1. This scaling is well followed by long chain segments (Figure 2c). For intermediate length scales (i.e., for 20 < N/p < 90), we conjecture that there is a crossover between these two behaviors.

The stretching exponent β from fitting a stretched exponential form to the autocorrelation function is a nonmonotonic function of N/p with larger values in the Rouse and reptation regimes and a minimum at the crossover around the entanglement length Ne. For large N/p, β decreases from 0.9 to a value 0.5 around Ne (Figure 2d) with decreasing N/p, independent of chain stiffness. Li et al.[16] also found that the minimum in β occurs around N/p ∼ Ne for kθ = 0. Previous works by Padding and Briels[15] and by Shaffer[14] suggest that this minimum in β is due to kinetic constraints on the chains.[30,31] Our results agree with this general conclusion, and the location of the minimum of β does track qualitatively with the reduction in entanglement length with increasing stiffness. That is the minimum occurs for N/p ∼ 40 for kθ = 1.5 and N/p ∼ 90 for kθ = 0. (For comparison, the entanglement lengths are ∼45 and ∼85, respectively.)

pan class="Chemical">Previously, almost all simulations of this model have used a purely repn>ulsive (rc = 21/6σ) interaction between beads to reduce compn>utational resources necessary to reach long times. Our results in Figure 2 clearly show that the interaction cutoff has almost no effect on the Rouse modes of the chain. We also carried out a primitive path analysis following Everaers et al.[23] and found that the entanglement length Ne is indepn>endent of rc.

Chain Length Effects

The effective monomeric relaxation rate of chains with different N is presented in Figure 3a for kθ = 0.75ε. For N < 40 the rate Weff is essentially independent of N/p (or that the relaxation time follows Rouse scaling, Figure 3b), except for the largest p where chain stiffness presumably enters. Chains of length N = 100 show that Weff changes by a factor of ∼2 over the p values studied. As Ne ∼ 45, it appears that entanglement effects on Weff start around Ne. The behavior of longer chains is reminiscent of the behavior that we found for N = 500 in Figure 2a. In addition to a plateau for large p (i.e., small N/p) there is a plateau for small p (i.e., large N/p). In this longer chain length small p regime, the effective monomeric relaxation rate is expected to scale as Weff ∼ [1 + (N/Ne)1.4]−1, which is what is found when chain length is varied (Figures 3a and 2c). The longest relaxation time (p = 1) (Figure 3d) shows that its dependence on chain length changes from N2 (unentangled, Rouse-like chains) for short chains to N3.4 (entangled) for long chains consistent with experimentally observed behavior.[32]

Figure 3

(a) Effective monomeric relaxation rates for chains of different lengths with stiffness kθ = 0.75ε. (b, c) Effective relaxation times of melt chains of different length for kθ = 0.75ε. (d) Longest relaxation time (p = 1) of chains of different length is compared with CC model. Lines are provided for guidance. (e) Stretching parameter β for chains of different lengths.

Consistent with these ideas, the effective relaxation times of the different modes, especially for large N/p fall on a reptation-motivated master curve for different chain lengths as shown in Figure 3c, while the shorter length scale modes follow Rouse scaling (Figure 3b). Figure 3c also suggests that this constraint dominated regime is applicable for N3.4/p2 > 107, which implies that only chains longer than N = 100 (∼ 2Ne) are in this regime. The stretching exponents β (Figure 3e) increase monotonically with N/p for short chains, while for longer chains there is a well-defined minimum in the vicinity of Ne.

Chain Crossability

Results for the chain-crossing (CC) model are shown in Figure 4. In the CC model, chain relaxations are Rouse-like for any given N as topological constraints are eliminated by chain crossing. The effective monomeric relaxation rates are constant for N/p > 5 (Figure 4a). For smaller N/p, the results are similar to those for the KG model (Figure 4b). The effective relaxation times for different length scales fall on a master curve which scales as τ ∼ (N/p)2 (especially for small p) as shown in Figure 4c. The stretching exponent β monotonically increases with N/p (see Figure 4d) as for unentangled noncrossing chains (Figure 3e). Comparison of the CC chains to the chains that cannot cross (Figure 4b) clearly shows that entanglements affect the relaxation starting at the intermediate length scales around Ne. The extent of the effect depends on N as Weff decreases gradually with N/p in the KG model, whereas it is constant in the CC model. For large N/p, the plateau in Weff for the CC model shows that it better satisfies the assumptions inherent in the Rouse model than the KG model. The much lower plateau values of Weff for the entangled chains clearly show that the presence of topological constraints serve to significantly slow down chain motion.

Figure 4

Chain crossing (CC) model: (a) Effective monomeric relaxation rates of melts of different chain lengths with kθ = 0.75ε. (b) Comparison of effective monomeric relaxation rates of chains with different lengths between CC model (upper set of curves) and KG model (lower set of curves). (c) Effective relaxation times of chains in melts for different chain length for kθ = 0.75ε. (d) Stretching parameter βp for melts with CC chains of different lengths.

Discussion and Conclusion

The major point of our analysis is that the Rouse modes of the chains show characteristic trends in the regime of chain lengths where entanglement effects become prominent. The relaxation of chain segments smaller than this topologically defined size follow Rouse-like behavior, although local correlations also lead to departure from Rouse behavior at the shortest length scales. Our results, especially on the monomeric relaxation rates Weff of the chains, and their mode number dependence, are in good agreement with the experimental results of Richter et al.[8,9] These workers found that the effective relaxation rates decrease dramatically when one goes through the entanglement crossover (see Figure 3a). However, these workers did not see any minimum in the apparent stretched exponents as a function of mode number—rather, the mode relaxations of these very modestly entangled chains were described by a normal exponential relaxation. The reasons for this difference in behavior between our simulations and the experiments might arise from the fact that the experiments were conducted for chains which are only weakly entangled. (pan class="Chemical">Previous theories on short chains, below the entanglement length, also show deviations from Rouse behavior due to local packing effects and local friction. The expn>erimental evidence for such non-Rouse behavior is weaker, and we do not understand the source of these discrepn>ancies between simulations and expn>eriments at this time.)

One question is why we choose to use the Rouse modes of the chains to describe motion rather than true normal modes,[33] which are guaranteed to be orthogonal to each other. Our logic has a twofold rationale. First, as discussed above, experiments are normally analyzed in terms of the Rouse model. Using the same language as the experiments is a prudent means to make connections of our models to reality. Second, previous works by Hess,[34] Ronca,[35] and Edwards[36] have shown that the equations of motions for long, entangled chains in a pan class="Chemical">polymer melt can be reduced to a generalized Rouse form, but with spectra at short and high frequency showing distinctly different depn>endences. This is precisely what is seen in Figure 3a, in good agreement with previous theoretical works.

Our Rouse mode results for pan class="Chemical">polymer melts show very interesting trends which greatly expn>and our understanding of relaxation of different length scales of chain and the role of entanglements on the dynamics of pan class="Chemical">polymers. For all chain lengths, the effective monomeric relaxation rate (friction) for p ∼ N is not affected as Weff is same for all N. The relaxation of intermediate length scales is greatly affected by entanglements and depends on N in the studied crossover regime.