Detailed analysis of Rouse mode and dynamic scattering function of highly entangled polymer melts in equilibrium

We present large-scale molecular dynamics simulations for a coarse-grained model of polymer melts in equilibrium. From detailed Rouse mode analysis we show that the time-dependent relaxation of the autocorrelation function (ACF) of modes p can be well described by the effective stretched exponential function due to the crossover from Rouse to reptation regime. The ACF is independent of chain sizes N for N/p < Ne (Ne is the entanglement length), and there exists a minimum of the stretching exponent as N/p → Ne. As N/p increases, we verify the crossover scaling behavior of the effective relaxation time τeff,p from the Rouse regime to the reptation regime. We have also provided evidence that the incoherent dynamic scattering function follows the same crossover scaling behavior of the mean square displacement of monomers at the corresponding characteristic time scales. The decay of the coherent dynamic scattering function is slowed down and a plateau develops as chain sizes increase at the intermediate time and wave length scales. The tube diameter extracted from the coherent dynamic scattering function is equivalent to the previous estimate from the mean square displacement of monomers.


Introduction
The dynamics of polymer chains in a melt is a complicated many-body problem where the motion of chains depends not only on different length scales but also time scales. It is well known that for short unentangled chains in a melt, excluded volume and hydrodynamic interactions are screened, and the viscoelastic properties of chains can be approximately described by the Rouse model [1][2][3]. If the polymer chains become long enough, the topological constraints dominate the dynamics of the chains. At intermediate time and length scales, the chains are assumed to move back and forth (reptation) within a tube-like region created by surrounding entangled chains and depending on the entanglement length N e . The dynamic behavior within this time frame is well described by the reptation theory of de Gennes, Doi and Edwards [2][3][4]. a e-mail: hsu@mpip-mainz.mpg.de b e-mail: kremer@mpip-mainz.mpg. de Rouse mode analysis provides a straightforward way of understanding the dynamics of single chains in a melt by mapping the trajectories of the chains into orthogonal Rouse modes. The chains are assumed to be Gaussian in this analysis. This has been applied widely to the analysis of experimental data and simulation data. Besides the original predictions of this model, in the literature [5][6][7][8][9] there also exist modified theoretical predictions including the excluded volume interaction, the intrinsic stiffness of chains, the intramolecular correlations in chains, and the topological constraints for analyzing the relaxation of Rouse modes.
In our recent work [10], we have investigated the chain conformations of fully equilibrated and highly entangled polymer melts and compared our simulation results to the related theoretical predictions in detail. The fully equilibrated and highly entangled polymer melts were generated by a novel and very efficient methodology through a sequential backmapping of soft-sphere coarse-grained configurations from low resolution to high resolution, and finally the application of molecular dynamics (MD) simulations of the underlying bead-spring model [11][12][13][14]. We have also studied the dynamics of fully equilibrated polymer melts, characterized by the mean square displacement of monomers, and determined the characteristic time scales: the characteristic time τ 0 , the entanglement time τ e ≈ τ 0 N 2 e , the Rouse time τ R ≈ τ 0 N 2 , and the disentanglement time τ d ≈ τ 0 N 2 (N/N e ) 1.4 , according to the predictions given by the Rouse model and the reptation theory, where N e is the entanglement length and N is the chain size. For N < N e in the Rouse regime, there exist exact solutions of almost all physical observables. Therefore, based on this work, we are interested in understanding to what extent the dynamics of single chains can be analyzed through the Rouse mode analysis and check the scaling predictions of the relaxation of the Rouse modes in the literature [1][2][3][4][5][6][7][8][9][14][15][16] whenever it is possible. On the other hand, Rouse mode decay should display a similar slowing down of the dynamic structure factor. The tube diameter introduced in the reptation theory can be extracted from the single chain dynamic structure factor measured via neutron spin echo (NSE) experiments [17][18][19]. Therefore we also study the dynamic scattering function from single chains in a melt, and check the consistency of the tube diameter estimated from different physical quantities [3,14,15,20,21]. We use the ESPResSo++ package [22] to perform the standard MD simulations with Langevin thermostat at the temperature T = 1ε/k B where k B is the Boltzmann factor to study fully equilibrated polymer melts consisting of 1000 chains of sizes N = 500, 2000 for k θ = 1.5 (τ 0 ≈ 2.89τ , N e ≈ 28 [10,12]) and of sizes N = 1000, 2000 for k θ = 0.0 (τ 0 ≈ 1.5τ , N e ≈ 87 [12]) in the framework of the standard bead-spring model [14] with a bond bending interaction parameter k θ at a volume fraction φ = 0.85.
The outline of the paper is as follows: Section II describes the theoretical background of the Rouse model and the Rouse mode analysis of highly entangled polymer melts. Section III describes the scaling behavior of dynamic structure factors and the comparison between theory and simulation. Finally, our conclusions are summarized in Section IV.

Rouse mode analysis
In the Rouse model, neglecting inertia effects Rouse chains undergo Brownian motion and therefore the Langevin equation of motion for the ith monomer is thus given by with the Rouse potential where b is the effective bond length, r i is the position vector of the ith monomer, and f R i (t) is a random force. The friction ζ and the random force f R i (t) are related by the fluctuation dissipation theorem, i.e., The Rouse modes X p (t) are defined as the cosine transforms of position vectors r i at time t for i = 1, . . . , N as given in reference [23], Here the p > 0 modes describe the internal relaxation of a chain of N/p monomers while the 0th (p = 0) mode corresponds to the motion of the center-of-mass of chain.
Since each Rouse mode X p (t) for p > 0 performs a Brownian motion in a harmonic potential, independent of each other, the cross-correlations should vanish, the normalized autocorrelation function is expected to follow an exponential decay, with and For p = N , τ N = τ 0 is the shortest relaxation time of the Rouse model, where τ 0 = ζb 2 /(3k B T π 2 ) is the characteristic relaxation time, while τ 1 for p = 1 is the longest relaxation time equal to the Rouse time, i.e. τ 1 = τ R . However, as the excluded volume interaction and topological constraint are taken into account, it has been pointed out in the literature [5,6,24] that simulation results of the normalized autocorrelation function are well described by the stretched exponential Kohlrausch-Williams-Watts (KWW) function, i.e.
where the KWW characteristic relaxation time τ * p and the stretching exponent β p depend on the mode p, and both are measures of importance of excluded volume interactions and topological constraints. The effective Rouse time of mode p is thus given by where Γ (x) is the gamma function. Figure 1 shows the typical relaxation of the time-dependent normalized autocorrelation function of modes p, X p (t)X p (0) / X 2 p , according to the definition of Rouse

696
The European Physical Journal Special Topics  modes X p (t) given in Equation (4). Data are for polymer melts of chain sizes N = 500 and 2000 with k θ = 1.5. The data sets shown in Figure 1a from left to right correspond to N/p = 5, 10, 25, 50, 125, 250, and 500. We see that for N/p < N e (N e ≈ 28 for k θ = 1.5) the relaxation of X p (t)X p (0) / X 2 p is independent of chain size N , namely, data for N = 500 and N = 2000 are on top of each other. Similar results are also observed for polymer melts of chain sizes N = 1000 and 2000 with k θ = 0.0 where the entanglement length N e ≈ 87 [12] (not shown). In the regime where N/p > N e the deviation between two data sets corresponding to the same value of N/p becomes more prominent as N/p increases since the entanglement effect between chain segments becomes more important. In Figure 1b, the plot of ln[ X p (t)X p (0) / X 2 p ] versus t/τ p with τ p = (N/p) 2 (see Eq. (5)) shows that the exponential decay is only valid for small values of N/p at initial relaxation time t. As p becomes small and t increases, one sees systematic deviations from the exponential decay due to the crossover from Rouse to reptation behavior. In Figure 1c, the curves indicate the best fit of our data for N = 500 to the theoretical prediction (Eq. (8)) with two fitting parameters β p and τ * p . Polymer chains containing N/p monomers are relaxed completely as X p (t)X p (0) / X 2 p → 0 for t 1. Therefore, one can estimate roughly the required relaxation time to relax very long chains in a melt through this curve fitting procedure. Fitted values of the stretching exponent β p and the estimates of the effective relaxation time τ eff,p from β p and τ * p (Eq. (9)) are shown in Figure 2.   We see that β p reaches a minimum around N/p ≈ N e for both cases k θ = 1.5 and 0.0 in Figure 2a. For chains having the same intrinsic stiffness, β p is independent of size N for N/p < N e while β p decreases with increasing size N for N > N e . Our results agree with the argument that the development of a minimum in β p is due to kinetic constraints on the chains [5,6], and the recent work [16] using the same simulation model but different cut-off for the LJ potential. At short length scales, the local constraints on the motion of monomers is related to the chain connectivity and the excluded volume interactions between monomers while at long length scales, the entanglement effect sets in that the motion of entangled chains is strongly hindered by topological constraints. Results of the effective relaxation time τ eff,p show the crossover behavior from the Rouse regime (τ eff,p ∼ (N/p) 2 ) to the reptation regime (τ eff,p ∼ (N/p) 3.4 ) as N/p increases.
In previous Monte Carlo simulations of polymer melts based on the bond fluctuation model [7], the authors showed that the Rouse model overestimates the correlation as p/N > O(10 −1 ) due to the lack of considering the intrinsic stiffness of the chains. Replacing random walk chains by freely rotating chains with a specific bond angle θ, the analytical expression of the amplitude X p (0)X p (0) is as follows, Figure 3a presents the rescaled amplitude of the autocorrelation function of Rouse mode p, X p (0)X p (0) /b 2 , plotted versus p/N . Here b 2 is determined by the best fit of our data to Equation (6) for small values of p/N . b 2 = 2.74σ 2 for k θ = 1.5 and b 2 = 1.74σ 2 for k θ = 0.0. Our fitted values of b 2 for both cases satisfy the relation b 2 = 2 b C ∞ predicted for freely rotating chains, where C ∞ is Flory's characteristic ratio (C ∞ = 2.88, 1.83 for k θ = 1.5, 0.0, respectively), and b = 0.964σ is the mean bond length [10,12]. Theoretical predictions given in Equations (6) and (10) are also shown for comparison. The deviation from the Rouse prediction for p/N > O(10 −1 ) is indeed seen as shown in reference [7]. Taking the estimates of cos θ from our simulations,  (6) and (10) are shown for comparison. b 2 = 2.74σ 2 for k θ = 1.5 and b 2 = 1.74σ 2 for k θ = 0.0. In (b), Equation (11) with C = 0.23 and 0.32 for k θ = 1.5 and k θ = 0, respectively, are also shown for comparison.
our data are described quite well by Equation (10) for p/N > O(10 −1 ). For small p/N (large N/p), one should expect that 4 sin 2 (pπ/(2N )) X 2 p reaches a plateau value b 2 . However, since at short length scale the local intramolecular correlations in the chains (the correlation hole effect) are important, a correction term [8,9,16] O((N/p) −1/2 ) is needed to be considered as follows, where C is a fitting parameter. The prediction is also verified as shown in Figure 3b.
Since the entanglement effect already sets in at N/p ≈ 28 for k θ = 1.5, we see that for chains of size N = 2000, the data for k θ = 1.5 fluctuate more than that for k θ = 0.0.

Dynamic structure factors
Dynamic behavior of polymer chains in a melt can also be described by the dynamic scattering from single chains. The coherent and incoherent dynamic structure factors S coh (q, t) and S inc (q, t) for single chains are defined by and The average · · · denotes an average over all chains, many starting states (t = 0), as well as over orientations of the wave vector q having the same wave length. Note that S coh (q) is the q-space representation of the Rouse modes. In the Rouse model, the displacement between monomer positions is Gaussian distributed since the force has a Gaussian probability distribution. Therefore, Equations (12) and (13) can be written as and where [r i (t) − r i (0)] 2 ∼ g 1 (t) is simply the mean square displacement of monomers. For short chains (N < N e ), the scaling predictions of S coh (q, t) and S inc (q, t) for t < τ e (N < N e ) and q > 2π/R g (N ) (R g (N ) is the radius of gyration of chains containing N monomers) predicted by the Rouse model are as follows, and where the factor W = 12kB T b 4 πζ and R g is the radius of gyration of the chain of size N . For t τ e one expects the standard diffusion behavior, i.e.
For long chains (N > N e ), the entanglement effect due to the topological constraints between chains in a melt becomes important. According to the reptation theory, local reptation processes for short time and escape processes from the tube (creep motion) for longer times and small values of q should be considered. Thus a pronounced plateau in S coh (q, t) is predicted and can loosely be interpreted as a Debye-Waller factor for τ e t τ d , S coh (q, t) S coh (q, 0) = 1 − q 2 d 2 /36 . (19) Note that here the tube diameter is defined by d = R e (N e ) where R e (N e ) is the end-toend distance of chains of size N e [3,15]. In our previous work [10], the definition of tube diameter d T is different by a factor of √ 3, i.e. our simulation estimate of tube diameter d T = (2 R 2 g (N e ) ) 1/2 = ( R 2 e (N e ) /3) 1/2 = d/ √ 3. In the deep-reptation regime, an analytic expression of the coherent dynamic structure used often in the neutronspin-echo (NSE) measurements for the determination of the tube diameter is given by [3,14,15,20,21] S coh (q, t) S coh (q, 0) with f (u) = exp(u 2 /36)erfc(u/6). Figure 4 shows the results of the incoherent dynamic structure factor S inc (q, t) according to Equations (13) and (17) for polymer melts of sizes N =500 and 2000 with k θ =1.5. The characteristic time scales, τ 0 , τ e , and τ R,N =500 taken from reference [10] are indicated by arrows. Since ln S inc (q, t) ∼ g 1 (t), the scaling laws of g 1 (t) showing  the crossover behavior from the Rouse regime to the reptation regime as t increases are also shown for comparison. In the Rouse regime where q > 2π/d T ≈ 1.25σ −1 , we see that ln S inc (q, t) ∼ t for t < τ 0 while ln S inc (q, t) ∼ t 1/2 for τ 0 < t < τ e . In the reptation regime one should expect that ln S inc (q, t) ∼ t 1/4 for τ e < t < τ R and 2π/R g (N ) < q < 2π/d T . We see that indeed S inc ∼ t 1/4 for 0.4σ −1 < q < 1.25σ −1 in Figure 4a and for 0.2σ −1 < q < 1.25σ −1 in Figure 4b. Here for k θ = 1.5, R g ≈ √ 0.4839N σ [10]. In Figure 4a, we have also observed that ln S inc (q) ∼ t 1/2 for t > τ R and q < 0.4σ −1 . Therefore, our results are in perfect agreement with the theoretical predictions.
For checking the scaling behavior of the normalized coherent dynamic structure factor predicted by the Rouse model, we plot ln[S coh (q, t)/S coh (q, 0)] versus q 2 t 1/2 /6 for polymer melts of size N = 2000 with k θ = 0.0 and 1.5 in Figure 5. We see that in the intermediate time regime τ 0 < t < τ e , our data are also in perfect agreement with the scaling law t 1/2 for q > 2π/d T ∝ N −1/2 e σ −1 . For τ d t τ e , one would expect that S coh (q, t)/S coh (q, 0) reaches a plateau predicted by the reptation theory    (Eq. (19)). Therefore, in Figure 6, we show the similar data as shown in Figure 5 but plot S coh (q, t)/S coh (q, 0) versus q 2 t 1/2 /6. We see that first data tend to collapse onto a single master curve for τ 0 < t < τ e as q increases and S coh (q, t)/S coh (q, 0) is independent of chain size N . As t > τ e , S coh (q, t)/S coh (q, 0) for two different chain sizes N start to deviate from each other. For polymer chains of size N = 2000 in both cases (k θ = 0.0 and k θ = 1.5), S coh (q, t)/S coh (q, 0) slows down as t τ e . It gives the first evidence from simulations that a pronounced plateau in S coh (q, t) shall occur as chain size N increases. Finally, we determine the tube diameter d T = d/ √ 3 for k θ = 1.5 by fitting our simulation data of S coh (q, t)/S coh (q, 0) to Equation (20) (see Fig. 7). Our results show that d T ≈ 7.10σ for N = 500, and d T ≈ 5.95σ for N = 2000 which are compatible to our previous estimate of d T ≈ 5.02σ [10]. incoherent structure factors for chains of two different sizes and stiffnesses. The relaxation of time-dependent autocorrelation functions of Rouse modes p is independent of chain size N for N/p < N e (Fig. 1a) where the entanglement length N e is obtained through the primitive path analysis (PPA) [10,12,25]. Estimates of the stretching exponent β p also show that the minimum of β p occurs in the vicinity of N/p ≈ N e (Fig. 2(a)). The crossover behavior of effective relaxation times τ eff,p from the Rouse regime to the reptation chain as N/p increases is verified as well. Since all these estimated quantities for N/p < N e behave differently from that for N/p > N e , we see that N e can also be determined roughly via the Rouse mode analysis of large polymer melt systems of two different chain sizes, and the value is consistent with that obtained through PPA. Our results are also in perfect agreement with the extended theoretical predictions considering the excluded interaction, topological constraint, the intramolecular interactions, and chain stiffness (Fig. 3).
The scaling behavior of coherent and incoherent dynamic structure factors strongly depends on the time t and wave length q. Therefore, it is a delicate matter to analyze the dynamic structure factors. However, we have provided evidence that the scaling behavior of S inc (q, t) ( Fig. 4) is compatible with the mean square displacement of monomers g 1 (t), and the crossover points characterized by the characteristic time scales τ 0 , τ e , and τ R and the corresponding wave length scales (the inverse of length scales) are consistent with each other. The slowing down of S coh (q, t) (Fig. 6) gives the first evidence that S coh (q, t) exhibits a plateau for τ e t τ d as the chain size N increases. The tube diameter d T = d/ √ 3 extracted from S coh (Fig. 7) is also in perfect agreement with d T obtained from g 1 (t) in reference [10].
We hope that the present work showing the detailed analysis of the dynamic properties of highly entangled chains covering the scaling regimes from Rouse to reptation will help for the further understanding of the dynamic behavior of the deformed polymer melts, polydisperse polymer melts, and the related experiments.