From caging to Rouse dynamics in polymer melts with intramolecular barriers: a critical test of the Mode Coupling Theory

By means of computer simulations and solution of the equations of the Mode Coupling Theory (MCT), we investigate the role of the intramolecular barriers on several dynamic aspects of non-entangled polymers. The investigated dynamic range extends from the caging regime characteristic of glass-formers to the relaxation of the chain Rouse modes. We review our recent work on this question, provide new results and critically discuss the limitations of the theory. Solutions of the MCT for the structural relaxation reproduce qualitative trends of simulations for weak and moderate barriers. However a progressive discrepancy is revealed as the limit of stiff chains is approached. This disagreement does not seem related with dynamic heterogeneities, which indeed are not enhanced by increasing barrier strength. It is not connected either with the breakdown of the convolution approximation for three-point static correlations, which retains its validity for stiff chains. These findings suggest the need of an improvement of the MCT equations for polymer melts. Concerning the relaxation of the chain degrees of freedom, MCT provides a microscopic basis for time scales from chain reorientation down to the caging regime. It rationalizes, from first principles, the observed devations from the Rouse model on increasing the barrier strength. These include anomalous scaling of relaxation times, long-time plateaux, and non-monotonous wavelength dependence of the mode correlators.


I. INTRODUCTION
The different dynamic processes present in amorphous polymers cover a extremely broad range of characteristic time scales, spanning from about 100 femtoseconds up to years. There are two main reasons for this. First, polymers are usually good glass-formers, which inherently exhibit a dramatic increase of the viscosity and structural (α-) relaxation times on approaching the glass transition temperature T g . As in non-polymeric glassformers, localized dynamic processes are also present below T g [1]. Second, their macromolecular character introduces relaxation processes related to the dynamics of the internal chain degrees of freedom. In the case of lowmolecular weight, nonentangled, polymer chains a sublinear increase (Rouse-like) arises in the mean squared displacement prior to the linear diffusive regime. In the case of high-molecular weight, strongly entangled, chains further sublinear regimes are found between the Rouse and linear regimes, which are usually interpreted in terms of reptation dynamics [2][3][4]. Such processes are inherent to chain connectivity, and extend over more time decades on increasing chain length. This broad time window for chain dynamics is observed even for temperatures far * Corresponding author: wabmosea@ehu.es above T g , when the structural relaxation extends over just a few picoseconds.
Another particular ingredient of polymer systems is that, apart from fast librations or methyl group rotations [5], every motion involves jumps over carbon-carbon rotational barriers and/or chain conformational changes. The corresponding map of relaxation processes is largely influenced by the barrier strength. Intramolecular barriers play a decisive role in, e.g., crystallization [6,7], adsorption onto surfaces [8,9], viscoelastic properties [10], or phase behavior of block copolymers [11]. Models for semiflexible and stiff polymers are of great interest in biophysics, since they can be applied to many important biopolymers as DNA, rodlike viruses, or actin filaments [12][13][14]. Thus, an understandig of the role of the intramolecular barriers on structural and dynamic properties of polymer systems is of practical as well as of fundamental interest in many fields of research.
A possible theroretical approach to this problem is provided by the Mode Coupling Theory (MCT) [15]. MCT introduces a closed set of coupled Mori-Zwanzig equations for the time dependence of density correlators. Static correlations enter the memory kernel as external input. Since the former can be related to the interaction potential through liquid state theories, MCT constitutes a first-principle theory for slow dynamics in complex systems. MCT has been developed over the last years to include systems with intramolecular structure (see e.g., Typeset by REVT E X [16][17][18]). This includes the approach of Chong and coworkers for simple polymer melts [19,20], based on the polymer reference interaction site model (PRISM) [21] for the static correlations. This approach was applied to the specific case of fully-flexible chains [19,20], i.e., without intramolecular barriers. A major success was the derivation, from first-principles, of the scaling laws predicted by the phenomenological Rouse model [2] for chain dynamics in nonentangled polymer melts. Likewise, it provided a unified microscopic description of both chain dynamics and the structural relaxation associated to the glass transition [19,20]. Some of us have recently performed a systematic computational investigation of the role of intramolecular barriers on the glass transition in polymer systems [22,23]. Starting from fully-flexible bead-spring chains, we introduced stiffness by implementing intramolecular barriers with tunable bending and torsion terms. In Ref. [23] we discussed the glass transition within the framework of the MCT for polymer melts, comparing simulations with numerical solutions of the MCT equations, in the long-time limit, for a broad range of barrier strength. This was possible since the quality of the PRISM approximations observed for fully-flexible chains [24] was not affected at all by the introduction of internal barriers in all the investigated range [23]. Numerical solutions reproduced trends in the nonergodicity parameters and MCT critical temperatures for weak and moderate barriers. However, strong discrepancies were observed on approaching the limit of stiff chains [23].
In this article we briefly summarize the main points of Refs. [22,23] and present extensive new results. Thus, we solve the time-dependent MCT equations for density correlators, and compare simulation and theoretical trends in α-relaxation times. We critically discuss the limitations of the theory by analyzing the accuracy of the assumed approximations. We find that dynamic heterogenities, static three-point correlations, and chain packing effects not accounted by MCT, do not play a major role on increasing the barrier strength. Indeed their effects seem to be weaker that in the case of fully-flexible chains. The reason for the observed discrepancies between simulation and theory for very stiff chains remains to be understood.
We also present here a systematic investigation on the effect of intramolecular barriers on the internal chain dynamics of nonentangled polymers. We analyze correlators for the chain normal modes (Rouse modes) and for bond reorientation. The simulations reveal strong deviations from the Rouse model on increasing chain stiffness. These include anomalous scaling of relaxation times [25], long-time plateaux, and nonmonotonous wavelength dependence of the mode correlators. We show that these anomalous dynamic features are reproduced by the corresponding MCT equations for the Rouse modes. This generalizes the analysis of Ref. [20], which was limited to fully-flexible chains, to polymers with intramolecular barriers of arbitrary strength. Thus, beyond usual phe-nomenological models for chain dynamics, MCT provides a unified microscopic picture down to time scales around and before the α-process [26].
The article is organized as follows. In Section II we describe the model and give simulation details. In Section III we compare simulation results with MCT solutions for several dynamic correlators probing structural relaxation and chain dynamics. In Section IV we discuss the possible origin of the observed deviations from MCT predictions. Conclusions are given in Section V.

II. MODEL AND SIMULATION DETAILS
We have performed molecular dynamics (MD) simulations of a bead-spring model with tunable intramolecular barriers. All chains consist of N m = 10 identical monomers of mass m = 1. Non-bonded interactions between monomers are given by a corrected soft-sphere potential where ǫ = 1 and σ = 1. The potential V (r) is set to zero for r ≥ cσ, with c = 1. 15. The values C 0 = 7c −12 and C 2 = 6c −14 guarantee continuity of potential and forces at the cutoff distance r = cσ. The potential V (r) is purely repulsive. It does not show local minima within the interaction range r < cσ. Thus, it drives dynamic arrest only through packing effects. Chain connectivity is introduced by means of a finitely-extensible nonlinear elastic (FENE) potential [27,28] between consecutive monomers: where K F = 15 and R 0 = 1.5. The superposition of potentials (1) and (2) yields an effective bond potential for consecutive monomers with a sharp minimum at r ≈ 0.985, which makes bond crossing impossible. Intramolecular barriers are implemented by means of the combined bending and torsional potentials proposed by Bulacu and van der Giessen in Refs. [29,30]. The bending potential V B acts on three consecutive monomers along the chain and is defined as where θ i is the bending angle between consecutive monomers i − 1, i and i + 1 (with 2 ≤ i ≤ N m − 1). We use θ 0 = 109.5 o for the equilibrium bending angle. The torsional potential V T constrains the dihedral angle φ i,i+1 . The latter is defined for the consecutive monomers i − 1, i, i + 1 and i + 2 (with 2 ≤ i ≤ N m − 2), as the angle between the two planes defined by the sets (i − 1, i, i + 1) and (i, i + 1, i + 2). The form of the torsional potential is The values of the coefficients a n are a 0 = 3.00, a 1 = −5.90, a 2 = 2.06, and a 3 = 10.95 [29,30]. The torsional potential depends both on the dihedral angle φ i,i+1 and on the bending angles θ i and θ i+1 . As noted in Refs. [29,30], the functional form (4) avoids numerical instabilities arising when two consecutive bonds align, without the need of imposing rigid constraints on the bending angles.
The investigated range of barrier strength corresponds to a strong variation of the chain stiffness. This can be quantified by the average end-to-end radii, R ee , of the chains. Thus, for the representative values (K B , K T ) = (0, 0), (8,0.2), (25,1) and (35,4), which cover the range from fully-flexible chains to the stiffest investigated chains, we find R ee = 3.6, 4.7, 5.5, and 6.5 at the respective lowest investigated temperature.

III. RESULTS: SIMULATIONS vs. THEORY
Before addressing the dynamic aspects of the present system, we want to stress that the investigated state points correspond to isotropic phases. We do not observe signatures of global orientational order induced by chain stiffness for the investigated state points. Thus, by measuring the quantity P 2 (Θ) = (3 cos 2 Θ − 1)/2, where Θ is the angle between the end-to-end vectors of two chains, and averaging it over all pairs of distinct chains, we obtain in all cases values |P 2 (Θ)| < 3 × 10 −3 . This is illustrated in Fig. 1, which shows the time evolution of P 2 (Θ) along a typical simulation window, both for fully-flexible chains, (K B ,K T ) = (0,0), and for the stiffest investigated chains, (K B ,K T ) = (35,4).
Local orientational order is also negligible. This is evidenced by computing a similar correlator P 2 (Θ; r cm ). In this case the average is performed only over pairs of distinct chains for which the distance between their respective centers-of-mass is less than r cm . Fig. 1 displays, for the former cases of fully-flexible and stiff chains, data of P 2 (Θ; r cm ) for several values of r cm . Negligible values of P 2 (Θ; r cm ) are obtained for r cm ≥ 2.0. Thus, the time average over the simulation time window, t sim , provides values | P 2 (Θ; r cm ≥ 2.0) tsim | < 0.02. By comparing both panels we conclude that chain stiffness does not induce a significant increase, if any, of local orientational order in the investigated systems. Weak local orientational order | P 2 (Θ; r cm ) tsim | 0.1 is observed only for very small interchain distances (see data for P 2 (Θ; r cm = 1.4)). Again, the introduction of chain stiffness does not induce clear changes in the orientational order at this length scale.

A. Structural relaxation
Now we characterize dynamic features associated to the caging regime and the structural α-relaxation.
. The sum is done over the coordinates r j of all the N monomers in the system. In all the cases the correlator is evaluated at the maximum, q max ≈ 7 [23], of the static structure factor S(q) = N −1 N j,k=1 exp[iq · (r j (0) − r k (0))] . We observe that increasing the strenght of the internal barriers at fixed ρ and T leads to slower dynamics. In the fully flexible case f s (q, t) decays to zero in a single step. On increasing the strength of the internal barriers f s (q, t) exhibits the standard behavior in the proximity of a glass transition. After the initial transient regime, f s (q, t) shows a first decay to a plateau, which is associated to the caging regime, i.e., the temporary trapping of each particle by its neighbors. At long times, a second decay is observed from the plateau to zero. This corresponds to the structural α-relaxation. Similar trends are displayed by the density-density correlator (not shown), defined as f (q, t) = ρ(q, t)ρ(−q, 0) / ρ(q, 0)ρ(−q, 0) , with ρ(q, t) = N j=1 exp[iq · r j (t)]. Let us define τ 0.2 as the time for which f s (q max , τ 0.2 ) = 0.2, and τ K q as that obtained from fitting the αdecay to a Kohlrausch-Williams-Watts (KWW) function, . Both times correspond to a significant decay from the plateau, and therefore can be used as operational definitions of the αrelaxation time τ α . Fig. 2b shows τ 0.2 as a funcion of T , for different values of the bending and torsional constants (results for τ K q are analogous). As observed in the analysis of the self-correlators, increasing the chain stiffness slows down the dynamics. At fixed temperature, the relaxation time for the stiffest investigated chains increases by several decades with respect to the fully-flexible case.
The dynamic trends displayed in Fig. 2 demonstrate that intramolecular barriers constitute an additional mechanism for dynamic arrest, coexisting with the general packing effects induced by density and temperature. Now we discuss this scenario within the framework of the (ideal) MCT. We briefly summarize the basic concepts and predictions of the theory. Extensive reviews can be found, e.g., in Refs. [15,[31][32][33][34][35]. On approaching a glass transition from the ergodic phase, density fluctuations decay in a slower fashion, remaining frozen in amorphous configurations when the glass transition occurs. MCT describes this phenomenon as a feedback mechanism driven by the slow density fluctuations. By starting from the fundamental Liouville equation of motion and using the Mori-Zwanzig projection operator formalism, an integrodifferential equation is obtained for the density-density correlator:f The memory kernel m(q, , is expressed in terms of the associated fluctuating forces R f q [34]. In order to provide a closed solvable form of Eq. (5), MCT introduces several approximations for the memory kernel. These approximations are: i) The fluctuating force can be splitted in two terms: the regular (fast) contribution, linear in density fluctuations, and a second term which can be expressed as a a linear combination of 'mode pairs', ρ k ρ q−k . The latter provides the slow contribution relevant for the structural relaxation. Thus, the first approximation consists of neglecting the fast contribution.
ii) Convolution approximation: three-point static correlations are approximated as products of static structure factors, iii) Kawasaki approximation: dynamic four-point correlations are factorized in terms of products of dynamic two-point correlations (see e.g., Ref. [34] for details). Nowadays there is plenty of evidence that this approximation worsens on decreasing temperature, specially around the time scale of the α-relaxation. The breakdown of the former approximation is usually assigned to the emergence of strong dynamic heterogeneities in the proximity of the glass transition [36][37][38][39][40].
After applying the former approximations, the memory kernel m(q, t) becomes bilinear in f (q, t): The vertex V(q, q − k) is given by: where c(q) is the direct correlation function [41]. Eq. (5) constitutes a closed set of coupled equations which can be solved self-consistently, provided S(q) and c(q) are known. The latter are external inputs in the MCT equations. Since static correlators contained in the vertex vary with the control parameters (e.g., density, temperature, or barrier strength), the MCT equations (7) establish a direct connection between statics and dynamics. Moreover, the former static correlators can be related to the interaction potential through closure relations from liquid state theories [41]. With this, MCT provides a first-principle approach for the slow relaxation of density correlators.
Recently, Chong and co-workers have derived MCT equations for simple models of polymer melts [19,20]. By exploiting the polymer reference interaction site model (PRISM) [21], the MCT equations are considerable simplified. This is achieved by replacing site-specific intermolecular surroundings of a monomer by an averaged one (equivalent site approximation), whereas the full intramolecular dependence is retained in the MCT equations [19,20]. The so-obtained scalar MCT equations of motion, memory kernel and vertex for polymer chains are formally identical to Eqs. (5,7,8). The polymer character of the system only enters implicitly through the PRISM relation [21] ρc(q) = 1/ω(q) − 1/S(q), which differs from the Ornstein-Zernike equation [41], ρc(q) = 1 − S −1 (q), for monoatomic systems. The quantity ω(q) is the chain form factor, defined as where r I a are the coordinates of the ath monomer in the Ith chain. N c is the total number of chains. The use of the former MCT equations is a priori justified for polymers of variable stiffness. Indeed, it has been shown that the PRISM approximations retain their validity not only in the fully-flexible limit [24], but also when strong intramolecular barriers are present [23]. For the case of the self-density correlators f s (q, t), the MCT equations are different from the monoatomic case. The former are obtained by summation of the diagonal terms of the self site-site density correlators. The latter are given by with indices defined as in Eq. (9). The correlators F s ab (q, t) are determined by solving the corresponding MCT matrix equations (see Ref. [20]).
Ideal MCT predicts a sharp transition from an ergodic liquid to an arrested state (glass), at a given value of the relevant control parameter (temperature in the present case). At the transition (or 'critical') temperature T = T c , the non-ergodicity parameter, defined as f q = lim t→∞ f (q, t), jumps from zero to a nonzero value f c q . The latter is called the critical non-ergodicity parameter. By taking the limit t → ∞ in the MCT equations, one finds the relation Eq. (11) always has the trivial solution {f q } = 0. Glassy states take place when solutions f q > 0 also exist. The temperature at which the jump from zero to nonzero solutions occurs defines T c . The corresponding solutions define the critical non-ergodicity parameters.
The separation parameter, ǫ T = (T − T c )/T c measures the distance to the critical temperature. We are interested in the behavior of f (q, t) in the ergodic fluid, i.e., for ǫ T > 0. For small values of ǫ T , MCT predicts several asymptotic laws for dynamic observables [42], which are characterized by several dynamic exponents. Thus, the decay from the plateau follows a von Schweidler expan- , and diffusivities and α-relaxation times obey The exponents of these asymptotic laws are related to the so-called exponent parameter λ, which is the only independent one. The MCT expression for λ is determined by the static correlators S(q) and c(q) evaluated at T = T c , and by the critical non-ergodicity parameters f c q (see, e.g., [23,43]). We solved Eqs. (5) and (11) by combining simulation results of ω(q) with the PRISM equation ρc(q) = 1/ω(q) − 1/S(q) and the Percus-Yevick closure relation [41]. Details of the numerical procedure for solving (11) can be found in Ref. [23]. Numerical integration of the density correlators was performed following the method of Ref. [44]. It often happens in the analysis of experiments or simulations that numerical solutions of the MCT equations are not available. In such cases, a phenomenological analysis can be performed, and the values of T c and the former dynamic exponents can be obtained as fit parameters from the experimental or simulation data. Consistency of the analysis requires that the exponents, which are obtained from independent fits to different scaling laws, are related to the same λ-parameter, as predicted by the theory. This consistency test was done in the analysis of our simulation data (see Refs. [22,23] for a detailed explanation), providing different values of T c and λ for each barrier strength. These values obtained from simulations can be compared with the values provided by solution of the MCT equations.
This comparison is shown in Fig. 3. Superscripts 'MD' and 'MCT' are used respectively for simulation and theoretical values. The data are represented as a function of the end-to-end radius, which quantifies chain stiffness. A clear correlation between the barrier strength and the values of T MD c and λ MD is unambiguously demonstrated. The interplay between monomer packing effects and intramolecular barriers [23] induces a progressive increase of T MD c at fixed density. We note that from the fully-flexible limit (K B , K T ) = (0, 0) to barriers with (K B , K T ) = (15, 0.5), the data sets for T MD c and T MCT c roughly display the same slope. As usual, there is a shift factor between simulation and theoretical temperatures (here T MD c /T MCT c ≈ 1.25), which may originate from the mean-field character of MCT [15]. The range of barrier strength for which T MCT c and T MD c are roughly parallel is significant. Indeed, for (K B , K T ) = (8, 0.2) the endto-end radius R c ee is a 30% larger than for fully-flexible chains. However, a strong discrepancy between simulation and theory becomes evident on increasing the barrier strength from (K B , K T ) = (15, 0.5). While beyond this point T MCT c seems to approach an asymptotic limit, T MD c increases up to 1.23 for the stiffest investigated chains. Similar trends are observed for the λ-exponent. Simulation values increase from λ MD = 0.76 for fully-flexible chains to λ MD = 0.86 for the stiffest investigated chains. The first ones are typical of simple glass-formers as the archetype hard-sphere fluid (λ = 0.74 [43]), where dynamic arrest is driven by packing effects. The largest ones are similar to those observed in realistic models of polymer melts which incorporate the full chemical structure of the chains [45][46][47]. On the contrary, the theoretical exponent exhibits a very weak variation, 0.71 ≤ λ MCT ≤ 0.72, over the investigated range of barrier strength.
These discrepancies in the case of strong intramolecular barriers are also reflected in the q-dependence of density correlators computed from simulations and from solution of the MCT equations. In both cases we fitted the corresponding α-decay to a KWW function (see above). Fig. 4 compares the q-dependence, at fixed T , of the KWW time τ K q for the self-correlators f s (q, t), as obtained from simulations and from theory. Results are presented for the fully-flexible case and for the stiffest investigated chains.
Before discussing such results, some points must be clarified. As mentioned above, the mean-field character of MCT usually yields a temperature shift between simulation and theory (see Fig. 3a). Moreover, MCT times are affected by an undetermined constant factor [15]. Thus, a proper comparison between theory and simulation for time-dependent correlators can be done by rescaling t by some characteristic relaxation time, and using a common separation parameter ǫ T [20]. This is the case for the data of Fig. 4. Thus, each data set is rescaled by the respective KWW time τ qmax corresponding to f s (q max , t).
Temperatures of both panels correspond to ǫ T ≈ 0.04 and 0.08 for respectively fully-flexible and stiff chains. In the rest of the article we will present several comparisons between simulation data and MCT solutions. It will be understood that the times and temperatures of the compared data obey the former criteria.
A clear disagreement between simulation and theoretical trends becomes evident in Fig. 4. The two sets of KWW times obtained from simulation show a rather different q-dependence, which is more pronounced for the stiff chains. On the contrary, after rescaling by τ qmax , the theoretical sets become esentially identical. Fig 5 shows a similar comparison between simulation and theory for the rescaled KWW times of the density-density correlators, f (q, t). For q ≥ q max the theory reproduces qualitatively the shape of the relaxation times, which are modulated by the respective static structure factor S(q) (not shown). However, MCT fails at reproducing the broad peak at intermediate q C ≈ 4 which is present in the simulation data. This failure was already noted for fully-flexible chains in Ref. [20], and is confirmed here for the general case with intramolecular barriers. Apparently (note the error bars), the peak does not shift significantly and decreases its intensity as chains become stiffer, leading to a shoulder. In previous works [48] on similar fully-flexible bead-spring chains of N m = 10, the value of q C has been identified with 2π/R g , where R g is the chain radius of gyration. Data of Fig. 5a do not seem compatible with this assignment. Appart from results for fully-flexible and stiffest investigated chains of N m = 10, we include data for (K B , K T ) = (15, 0.5) of additional simulations with N m = 21. With this, the data sets of Fig. 5a cover a significant variation in R g . Namely, for (K B , K T ) = (0, 0), (35,4), and (15,0.5) we respectively find 2π/R g = 4.2, 2.8, and 2.0. Thus, the observation q C ≈ 2π/R g for fully-flexible chains is apparently fortuitous. The associated length scale 2π/q C ≈ 1.6σ rather seems to be a characteristic feature which does not depend significantly on the barrier strength. We will come back to these points in Section IVc. In this subsection we compare simulation and theoretical results for the dynamics of the Rouse modes. center-of-mass coincides with X 0 (t)/ √ N m . The mode correlators are defined as C pq (t) = [ X p (0) · X q (t) − δ 0,p×q X p (0) · X q (0) ]/3N m . For p, q > 0 we define the matrixĈ pq (0) ≡ C pq (0). In the Rouse model the stochastic forces are fully spatial and time uncorrelated, i.e., f j (t) · f k (t ′ ) = 6ζk B T δ jk δ(t − t ′ ). These properties of the random forces lead to orthogonality and exponentiality of the Rouse modes [2]. Thus, the mode correlators obey C pq (t) =Ĉ pq (0) exp[−t/τ p ], withĈ pq (0) = δ pq (b 2 /24N 2 m ) sin −2 [pπ/2N m ] and τ p = (ζb 2 /12k B T ) sin −2 [pπ/2N m ]. Accordingly, for p ≪ N m the quantitiesĈ pp (0) and τ p scale as ∼ p −2 .
In the following we show how the former scaling properties are strongly altered by the introduction of intramolecular barriers. This is demonstrated in Fig. 6 for the case of intrachain static correlations.
We show the off-diagonal terms of Ψ pq (0) = X p (0) · X q (0)/(X p (0)X q (0)) (the diagonal terms are trivially Ψ pp (0) ≡ 1). Data for fully-flexible chains exhibit small deviations from orthogonality, indeed |Ψ pq (0)| < 0.05 for all p = q, independently of T . Instead, orthogonality is clearly violated for strong intramolecular barriers. Offdiagonal terms can take values of even a 60% percent of the diagonal ones. Moreover, deviations are enhanced by decreasing temperature. Fig. 6c shows results for the unnormalized diagonal termsĈ pp (0) (see above). In the low p-range the data can be described by an effective power lawĈ pp (0) ∼ p −x . For fully-flexible chains we find approximate gaussian behavior,Ĉ pp (0) ∼ p −2.2 [28]. However, the introduction of internal barriers leads to strong non-gaussian behavior. On increasing the barrier strength, the effective exponent x increases up to a value of 3.8 for the stiffest case, (K B , K T ) = (35,4), at the lowest T . The most local effects of the intramolecular barriers are manifested by flattening ofĈ pp (0) at large p. The trends observed for intrachain static correlations have their dynamic counterparts. Fig. 7 shows the relaxation times τ p , of the normalized mode correlators Φ pp (t) = C pp (t)/Ĉ pp (0), as a function of the mode index p. We display data for several values of the bending and torsion constants (K B , K T ) and temperatures T . The relaxation times have been operationally defined as Φ pp (τ p ) = 0.3. Data can be again described at low-p by an effective power-law τ p ∼ p −x . The observed trends are analogous to those found for the static correlations (Fig. 6c). Rouse behavior (x = 2) is observed only in the fully-flexible limit. Again, as for the static amplitudeŝ C pp (0), x is weakly dependent on T [50] but strongly dependent on the barrier strength, taking higher values for stiffer chains. The x-values forĈ pp (0) and τ p at the same (K B , K T ) and T are similar. This suggests that the structural origin of the observed dynamic anomalies is mainly controlled by intrachain static correlations. Fig. 8a shows simulation results for the normalized mode correlators Φ pp (t), for (K B , K T ) = (35,4), at T = 1.48. Times are rescaled by the relaxation time of the first mode, τ 1 . Several salient features are revealed. First, the unambiguous presence of a long-time plateau for the modes p = 3 and p = 5, followed by an ultimate slow decay. It must be stressed that this feature is not related to the structural α-relaxation. Indeed, the plateau arises at times far beyond the α-time scale (τ α ∼ 5 × 10 −3 τ 1 for the considered T ). This feature is instead intimately connected to the relaxation of the internal torsional degrees of freedom of the chain. Indeed we observe (not shown) that for fixed bending constant K B , the long-time plateau tends to vanish as the value of the torsional constant K T is decreased.
The observed long-time plateau constitutes a clear breakdown of the Rouse model, which predicts single, purely exponential decays of the mode correlators (see above). Its origin can be temptatively understood as follows. The relaxation of the pth-mode is equivalent to the relaxation of a harmonic oscillation of wavelength N/p . In the case of strong torsional barriers, the wavelengths of some particular modes probe characteristic lengths over which chain deformation involves a strong energetic penalty (due to the presence of the barriers). Thus, at the time scales for which the barrier amplitudes are probed, the relaxation of such modes becomes strongly hindered, leading to the observed long-time plateau regime and ultimate slow relaxation. Another intriguing feature of Fig. 8a, also inconsistent with the Rouse model, is the non-monotonous p-dependence of the mode correlators at intermediate times prior to the long-time plateau (see data for p > 4). Now we demonstrate that all the former dynamic features can be rationalized in terms of the PRISMbased MCT approach of Chong et al.. As exposed in Refs. [19,20], the MCT equations for the unnormalized Rouse correlators C pq (t) are derived as the q → 0 limit of the equations for the self site-site density correlators F s ij (q, t). The equations for C pq (t) read [19,20]: [51][52][53]. Thus, prior to solve Eq. (12), we obtained the density-density correlators f (k, t) and self site-site density correlators F s ij (k, t) from their respective MCT equations (see Ref. [20]). The static quantitiesĈ pq (0) andĈ −1 pq (0), which also enter Eq. (12) as external inputs, were directly computed from the simulations at the respective lowest investigated temperature. Fig. 8 shows a comparison, at ǫ T ≈ 0.2, of the MCT solutions for the normalized mode correlators Φ pp (t) [panel (b)] of the stiffest investigated chains, with the respective simulation results previously discussed [panel (a)]. A full correspondence between MCT solutions and simulation trends is obtained. These include the long-time plateaux for p = 3 and p = 5, as well as the sequence in the complex, non-monotonous p-dependence for p > 4 at intermediate times. As previously done for the simulation data, we can obtain the theoretical relaxation times τ p from the condition Φ pp (τ p ) = 0.3 in the theoretical correlators. The p-dependence of the simulation and theoretical times are compared in Fig. 7, at common ǫ T ≈ 0.2, for several values of (K B , K T ). Again, MCT solutions are in semiquantitative agreement with the anomalous trends of simulations, with similar exponents for the effective power-laws.
As we observed in Fig. 6 for Ψ pq (0), there are offdiagonal terms of the intrachain static correlations which are non-orthogonal. This non-orthogonality persists over long time scales, as can be seen in Fig. 9. The latter shows simulation and theoretical results for normalized Rouse cross-correlators Φ pq (t), with p = 3 and q = 1, 3, 5, 7, 9. Data correspond to the same temperatures and barrier strength (the stiffest investigated case) of the diagonal correlators of Fig. 8. Again, MCT qualitatively reproduces simulations trends for the case of the off-diagonal terms. Finally, it is worth noting that the good agreement between simulations and MCT for the Rouse correlators is similar for other observables probing chain dynamics. The reason is that, through the transformation X p (t) = Nm j=1 P jp r j (t) (see above), such observables can be expressed in terms of the Rouse diagonal and cross-correlators [2]. An example is given by the orientational bond correlator P b (t) = b(0) · b(t) / b 2 (0) , where b(t) is the bond vector joining two consecutive monomers. Following the former transformation we find b(0) · b(t) ≡

Nm−1 j=1
Nm−1 p,q=0 [P −1 p,j+1 − P −1 p,j ][P −1 q,j+1 − P −1 q,j ] X p (t) · X q (0) , where P −1 is the inverse of the matrix of coefficients P jp . Note that this expression is exact (the Rouse model makes the approximation X p (t) · X q (0) = 0 for p = q). Since MCT solutions pro- vide the Rouse correlators for all (p, q), insertion of these in the former exact expression directly provides P b (t). Fig. 10 shows simulation and MCT results of P b (t) for several values of (K B , K T ) from the fully-flexible limit to the stiffest investigated chains. As in previous figures, times are rescaled by the respective τ 1 , and data in both panels correspond to a common separation parameter ǫ T ≈ 0.2. MCT reproduces semiquantitatively the observed simulation trends. These include, on increasing barrier strength, a relative speed up and slowing down (in terms of the scaled time t/τ 1 ) of respectively the short-time and long-time dynamics. MCT also accounts for the emergence, for strong barriers, of a plateau at t/τ 1 ∼ 10 −2 and a change in the concavity of the decay.

IV. DISCUSSION
On Section III we have shown that, concerning the critical temperature T c , MCT reproduces qualitative simulation trends for low and moderate barriers. However a strong disagreement is found on approaching the limit of stiff chains. We have also found a clear discrepancy in the trends of the λ-exponent, with a nearly constant value from theory and strongly barrier-dependent values from simulations. In this Section we discuss possible origins of these discrepancies.

A. Three-point static correlations
In Ref. [23] we showed that the failure of the MCT predictions for strong intramolecular barriers was not apparently related with the breakdown of the PRISM approximations, which are invoked in the derivation of the MCT equations for polymers. Indeed the quality of such approximations appeared to be the same for all the range of barrier strength here investigated by simulation and MCT. Despite the mentioned discrepancies between theory and simulation, the phenomenological analysis of simulation results in terms of a huge set of general asymptotic laws of MCT was consistent [22,23]. This means that the dynamic exponents involved in the different tested laws could be, in each case, related to a same λ MD . As we discussed in [23], such scaling laws are a mathematical consequence of the bilinear dependence of the memory kernel on the density correlators (see Eq. (7)). The specific numerical values of λ (and by transformation, of the other dynamic exponents) are determined by the static quantities entering the vertex (8) [23,43]. Given the consistency of the phenomenological analysis we speculated that, by retaining the bilinear form of the MCT memory kernel, there may be missing static contributions in the vertex which are not significant for low barriers, but become increasingly important as the limit of stiff chains is approached. Including them and solving the MCT equations accordingly, might raise the theoretical values of T c and λ, leading to a better agreement with the simulation trends.
Thus, we suggested that intrachain three-point static correlations should be explicitly included in the MCT vertex. Chain stiffness induces a strong directionality in the intrachain static correlations, at least at nearneighbor distances. It has been shown that directionality in static correlations can break the static convolution approximation of MCT, Eq. (6). A well-known example is given by silica, a network-forming system. For the latter the inclusion of three-point static correlations in the MCT vertex significantly improves the comparison between theory and simulation, with respect to the solutions obtained under the convolution approximation [54].
The calculation of the three-point static correlations involved in Eq. (6) is very demanding. This is because most of the computational time is consumed by the interchain three-point correlations. For intrachain three-point correlations the computation is not demanding. Fortunately, in the present case only the latter is necessary, since the directionality of correlations is only relevant along the chain. Thus, the convolution approximation is retained for interchain correlations, and it is modified only to include the intrachain three-point correlations.
Figs. 11 and 12 show representative tests of the convolution approximation for respectively fully-flexible and stiffest investigated chains. Following the scheme proposed in Ref. [20], the vectors q, k and p = q − k define the sides of a triangle, the first two enclosing an angle φ given by cos φ = (q 2 + k 2 − p 2 )/2qk. Panels (a) and (b) in Fig. 11 show a test of the corresponding expression for an equilateral triangle, ω 3 (q, q, q) = ω 3 (q). Panels (a) and (b) in Fig. 12 show a similar test for equal moduli k = q and all the relative orientations (given by cos φ) of q and k. In other words, we test the approximation ω 3 (q, q, p = q 2(1 − cos φ)) = ω 2 (q)ω(p). Data in Fig. 12 are represented as a function of cos φ for two characteristic wave vectors, corresponding to the first minimum and second maximum of the respective ω 3 (q, q, q) (see Fig. 11).
As already noted in Ref. [24], the convolution approximation for intrachain static correlations provides a good description of ω 3 in the fully-flexible limit. As expected, the quality of the approximation decreases by introducing intramolecular barriers. Still it constitutes a good approximation for all the investigated barrier strength. In the case of wave vectors around the first peak of S(q), q max ≈ 7, the quality is almost unaffected by the barrier strength, i.e., the terms c 3 (q max , q max − k) will be small even for the stiffest investigated chains. It must be noted that the MCT kernel is usually dominated by the contributions around q max . Thus, the former observations suggest that the inclusion of the three-point static correlations will modify weakly the MCT solutions obtained under the convolution approximation. We confirm this by obtaining numerical solutions with the vertex (13), for which we compute the input quantities involved in Eqs. (14,15) directly from the simulations. The so-obtained values of the critical temperature T c and λexponents raise by ≈ 1% as much, even for the stiffest chains, with respect to the previous values ( Fig. 3) found under the assumption c 3 = 0. With all this, we conclude that the observed discrepancies between simulation and theoretical trends of T c and λ are not related to the breakdown of the convolution approximation for static threepoint correlations. The latter indeed retains its validity for all the investigated range of barrier strength.

B. Dynamic heterogeneities
It is well-known that the quality of the Kawasaki approximation for dynamic correlations (see above) breaks on decreasing temperature. This feature is specially critical around the α-time scale [36][37][38][39][40], leading to the complete failure of the MCT predictions associated to it, as the power law behavior D −1 , τ α ∼ (T − T c ) −γ , or the time-temperature superposition of density correlators. This breakdown is usually assigned to the emergence of strong dynamic heterogeneities on approaching the glass transition [36][37][38][39][40]. Having noted this we may speculate that, for some reason to be understood, increasing the barrier strength strongly enhances dynamic heterogeneities. This might result in a lower quality of the MCT and might be the reason for the observed discrepancies between theory and simulation trends for T c . Now we show that this is not actually the case, and that there is no correlation between barrier strength and enhanced dynamic heterogeneity. Non-gaussian parameters provide a simple way of quantifying the strength of the dynamic heterogeneity. They display large positive values at the time scales for which the respective van Hove function strongly deviates from the gaussian limit. This occurs when a significant fraction of particles has performed displacements very different from the average. Here we discuss the two most popular non-gaussian parameters. The standard or 'fast' parameter is given by α 2 (t) = 3 5 (∆r(t)) 4 / (∆r(t)) 2 2 − 1. The 'slow' parameter introduced by Flenner and Szamel [56] is defined as γ 2 (t) = 1 3 (∆r(t)) 2 1/(∆r(t)) 2 − 1. By construction α 2 (t) and γ 2 (t) are identically zero for a gaussian form of the van Hove self-correlation function. However, as noted in [56], large positive values of these parameters have a very different microscopic origin, reflecting distinct aspects of dynamic heterogeneity. In the case of the fast parameter α 2 (t), large values originate from a significant population of particles which have performed much larger displacements than the average. This effect is generally maximum at the time scale t * around the end of the caging regime. Thus, α 2 (t) increases from zero at t = 0 up to a maximum at t * , and decays to zero at longer times. The increase of the maximum α 2 (t * ) on decreasing temperature reflects a progressive enhancement of dynamic heterogeneity, at the decaging process, on approaching the glass transition.
The slow parameter γ 2 (t) exhibits analogous trends for the temperature and time-dependence. However the maximum of γ 2 (t) takes place at much longer scales than t * , namely around the α-relaxation time τ α . This effect originates from a significant population of particles which at t ∼ τ α have performed much smaller displacements than the average [56]. Fig. 13 compares simulation results of α 2 (t) and γ 2 (t), for the fully-flexible case and for very stiff chains. For a fair comparison we have selected temperatures at which the respective α-relaxation times are similar. These are T = 0.50 and T = 1.05, for respectively fully-flexible and stiff chains, and correspond to a separation parameter ǫ T ∼ 0.04 (see Fig. 3). Let us remind that the αtime scale can be estimated, e.g., as f (q max , t) = 0.2. In Fig. 14 we display f (q max , t) for both systems at the former temperatures, showing that the respective α-time scales are roughly the same, τ α ∼ 10 4 . The decaging times, which can be estimated from the start of the decay from the plateau in f (q max , t), are also roughly the same, t * ∼ 500. This equivalence is indeed reflected in the trends of the non-gaussian parameters in Fig. 13. Thus, in both systems α 2 (t) is peaked at t * ∼ 500 and γ 2 (t) is peaked at τ α ∼ 10 4 . Having noted this equivalence of time scales, data in Fig. 13 do not reflect any enhancement of dynamic heterogeneity on increasing the barrier strength. Actually, the opposite effect is suggested by the lower values of α 2 (t) and γ 2 (t) for stiff chains with respect to the fullyflexible case. With this, we discard a major role of dynamic heterogeneities as the reason for the observed discrepancies between simulations and MCT solutions for very stiff chains.

C. Chain packing
As mentioned in Section III and shown in Fig. 5, MCT fails at reproducing the peak around q C ∼ 4 for the qdependence of the KWW times of density-density correlators. As noted in Ref. [20] for the fully-flexible case, the origin of this peak may be related to dynamic correlations between centers-of-mass of the chains. The latter might arise from the effective packing between the polymer coils, interacting as fully penetrable spheres of size R g . This interpretation is not clear in view of the results of Fig. 5, since the value of q C does not seem to be related with 2π/R g . Having noted this, Chong et al. found that the incorporation of the static correlations between the centers-of-mass in the MCT equations did not improve the description of the simulation results. As shown in [20], this is not unexpected due to the almost featureless form of the static structure factor of the centers-of-mass S CM (q). The inset of Fig. 5a shows simulation results of S CM (q), for the same barrier strength and temperatures of the KWW times in the main panel. The introduction of chain stiffness does not induce significant features in S CM (q), appart from a stronger signal at low q. The latter indeed suggests that packing effects between the polymer coils are even weaker than for the fully-flexible case. Within the former interpretation, this would be consistent with the lower intensity of the mentioned peak of τ K q /τ qmax at q C ∼ 4. All these results suggest that discrepancies between simulations and MCT on increasing chain stiffness are not related to a dynamic coupling, not accounted for within the theory, to the slow modes at q C ∼ 4. Indeed, this coupling seems to be weaker for stiff chains.

D. Outlook
In summary, in this section we have discussed possible origins for the discrepancies, concerning the structural relaxation, between simulations and MCT on increasing barrier strength. We discard a major role, in comparison with the fully-flexible case, of three-point static correlations, dynamic heterogeneities and chain packing. These effects become even weaker on increasing chain stiffness. We remind that such effects are indeed neglected in the derivation of the MCT equations used here (see Section III). Results in this section suggest that this is not less justified for very stiff chains than for fully-flexible ones.
How to improve the theory to account for dynamic trends in stiff chains is an open question. A way might be the reformulation or extension of the MCT equations, retaining the bilinear form of the kernel, in terms of new dynamic observables coupled to density fluctuations. Such observables can be adequate for describing particular dynamic features which are not captured by the usual observables, i.e., the number density fluctuations ρ(q, t). Some examples are roto-translational site fluctuations adapted to the molecular symmetry, as has been shown, e.g., for dumbbell-like molecules [17,57] or for a simple model of orthoterphenyl [18]. The inclusion of density fluctuations of centers-of-mass improve results for rigid molecules [18] concerning a peak in τ K q /τ qmax at intermediate q [58], similar to that observed here at q C ∼ 4. As discussed above, this is not the case for polymer chains. Though there is no characteristic symmetry in polymer chains, roto-translational density fluctuations can also be defined over sites a, b at some characteristic distance |a− b|, perhaps probing the relevant length scale 2π/q C , which according to the data of Fig. 5 seems to depend weakly on the barrier strength. Whether this procedure may improve the agreement between MCT and simulations remains to be solved.

V.CONCLUSIONS
By means of simulations and solution of the equations of the Mode Coupling Theory, we have studied the role of intramolecular barriers, of arbitrary strength, on several aspects of polymer dynamics. The investigated dynamic range extends from the caging regime characteristic of glass-formers to the relaxation of the chain Rouse modes. Solutions of the MCT for the structural relaxation reproduce qualitative trends of simulations for weak and moderate barriers. However, a progressive discrepancy between MCT and simulations is revealed as the limit of stiff chains is approached. We have tested the validity of several assumptions inherent to the theory. Deviations from the theoretical predictions do not seem related with dynamic heterogenities, which indeed are not enhanced by increasing the barrier strength. Moreover, the convolution approximation for three-point static correlations retains its validity for stiff chains. Even the role of slow modes at intermediate length scales, not accounted for by MCT, becomes less significant on increasing chain stiffness. At this point it is not clear how to improve the MCT equations in order to remove the mentioned discrepancies for the case of stiff chains. We have suggested the possibility of formulating the MCT equations in terms of roto-translational density fluctuations over specific length scales.
Concerning the relaxation of the chain degrees of freedom, MCT provides a microscopic basis for the observed deviations from the Rouse model on increasing the barrier strength. These include anomalous scaling of relaxation times, long-time plateaux, and non-monotonous wavelength dependence of the mode correlators. Beyond usual phenomenological models for chain dynamics (the Rouse model being the corresponding one for fullyflexible chains), MCT provides a unified microscopic picture down to time scales around and before the α-process, which is not accounted for within the mentioned models.