The Role of Intramolecular Barriers on the Glass Transition of Polymers: Computer Simulations vs. Mode Coupling Theory

We present computer simulations of a simple bead-spring model for polymer melts with intramolecular barriers. By systematically tuning the strength of the barriers, we investigate their role on the glass transition. Dynamic observables are analyzed within the framework of the Mode Coupling Theory (MCT). Critical nonergodicity parameters, critical temperatures and dynamic exponents are obtained from consistent fits of simulation data to MCT asymptotic laws. The so-obtained MCT $\lambda$-exponent increases from standard values for fully-flexible chains to values close to the upper limit for stiff chains. In analogy with systems exhibiting higher-order MCT transitions, we suggest that the observed large $\lambda$-values arise form the interplay between two distinct mechanisms for dynamic arrest: general packing effects and polymer-specific intramolecular barriers. We compare simulation results with numerical solutions of the MCT equations for polymer systems, within the polymer reference interaction site model (PRISM) for static correlations. We verify that the approximations introduced by the PRISM are fulfilled by simulations, with the same quality for all the range of investigated barrier strength. The numerical solutions reproduce the qualitative trends of simulations for the dependence of the nonergodicity parameters and critical temperatures on the barrier strength. In particular, the increase of the barrier strength at fixed density increases the localization length and the critical temperature. However the qualitative agreement between theory and simulation breaks in the limit of stiff chains. We discuss the possible origin of this feature.


I. INTRODUCTION
Since they do not easily crystallize, polymers are probably the most extensively studied systems in relation with the glass transition phenomenon. Having said this, their macromolecular character, and in particular chain connectivity, must not be forgotten. The most evident effect of chain connectivity is the sublinear increase of the mean squared displacement (Rouse-like) [1] arising after the decaging process, in contrast to the linear regime found in non-polymeric glass-formers. Moreover, in the case of strongly entangled polymer chains, the reptation model predicts other two sublinear regimes between the Rouse and linear regimes [1,2,3].
Another particular ingredient of polymers is that, apart from fast librations or methyl group rotations [4], every motion involves jumps over carbon-carbon rotational barriers and/or chain conformational changes. Intramolecular barriers play a decisive role in the physical properties of polymer systems. Thus, they are responsible of partial or total crystallization [5,6]. They also enhance dynamic features which are usually associated to reptation [7,8], which controls rheological properties [3]. Models for semiflexible polymers are of great interest, since they can be applied to many important biopolymers such as proteins, DNA, rodlike viruses, or actin filaments [9,10,11]. Moreover, chain stiffness seems to play an * Corresponding author: sckbernm@ehu.es important role in the absorption behavior of polymers at interfaces [12,13]. Thus, an understandig of the role of intramolecular barriers on structural, dynamic and rheological properties of polymers is of practical as well as of fundamental interest.
In this work we investigate, by means of molecular dynamics simulations, the role of intramolecular barriers on the glass transition of polymer melts, by systematically tuning barrier strength in a simple bead-spring model. We discuss the obtained results within the framework of the Mode Coupling Theory (MCT) of the glass transition [14,15,16,17,18,19]. We extend preliminary results reported by us in Ref. [20] by testing a large set of predictions, including the factorization theorem and time-temperature superposition principle. A consistent set of dynamic exponents associated to asymptotic scaling laws is obtained. By increasing the barrier strength a crossover is observed for the values of the so-called λexponent. In the limit of fully-flexible chains λ takes values ∼ 0.7, characteristic of simple fluids dominated by packing effects. On the contrary, for strong intramolecular barriers the λ-values approach the upper limit λ = 1 characteristic of higher-order MCT transitions. The latter arise in systems with different coexisting mechanisms for dynamic arrest [21,22,23]. In the system investigated here, the obtained results suggest an interplay between general packing effects and polymer-specific intramolecular barriers.
Chong and co-workers [24,25] have recently presented an extension of the MCT to simple fully-flexible bead-Typeset by REVT E X spring models of polymer systems, in the framework of the polymer reference interaction site model (PRISM) [26,27,28]. In this formalism each molecule is divided into interaction sites corresponding to monomers. A key assumption of the PRISM is the replacement of the sitespecific intermolecular surroundings of a monomer by an averaged one (equivalent site approximation), while keeping the fully intramolecular dependence. We have tested the PRISM approximations used by MCT in the polymer model here investigated, which incorporates intramolecular barriers. Likewise, we have solved the MCT equations for the location of the MCT 'glass transition' temperatures (MCT critical temperatures) and for the nonergodicity parameters, which quantify the stability of density fluctuations in the reciprocal space. We compare solutions of the MCT equations with the results obtained from the phenomenological analysis of the simulation data. We observe that the theory reproduces qualitative trends in the nonergodicity parameters and critical temperatures. However, the agreement breaks as the limit of stiff chains is approached. We discuss the possible origins of this feature.
The article is organized as follows. In Section II we describe the model and give simulation details. Static correlators are shown in Section III. Moreover the PRISM approximations are tested for representative values of the barrier strength. Section IV presents qualitative dynamic trends as a function of the barrier strength. In Section V we summarize the universal predictions of the MCT and the equations of motion of the version for polymer melts introduced by Chong and co-workers. In Section VI we perform a phenomenological analysis of simulation data within the MCT, by testing universal scaling laws and deriving their associated dynamic exponents. In Section VII we compare the results of the former analysis with numerical solutions of the MCT equations. We discuss the observed differences for stiff chains in Section VIII. Conclusions are given in Section IX.

II. MODEL AND SIMULATION DETAILS
We have performed molecular dynamics (MD) simulations of a bead-spring model for which we have implemented bending and torsional intramolecular barriers. The monomer-monomer interaction is given by a corrected soft-sphere potential V (r) = 4ǫ[(σ/r) 12 where ǫ = 1 and σ = 1. The potential V (r) is set to zero beyond the cutoff distance r ≥ cσ, with c = 1. 15. The values C 0 = 7c −12 and C 2 = 6c −14 guarantee continuity of potential and forces at 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. Along the chain backbone, of N monomers, an additional finitelyextensible nonlinear elastic (FENE) potential [29,30] is used to introduce bonds between consecutive monomers: where K F = 15 and R 0 = 1.5. The superposition of potentials (1) and (2) provides 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 a combined bending V B , and torsional potential V T . We have used the potentials proposed by Bulacu and van der Giessen in Refs. [8,31]. The bending potential acts on three consecutive monomers along the chain. The angle between adjacent pairs of bonds is mantained close to the equilibrium value θ 0 = 109.5 o by the cosine harmonic bending potential where θ i is the bending angle between consecutive monomers i − 1, i and i + 1 (with 2 ≤ i ≤ N − 1). The torsional potential constrains the dihedral angle φ i,i+1 , which is defined for the consecutive monomers i − 1, i, i + 1 and i + 2 (with 2 ≤ i ≤ N − 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 this 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 [8,31]. 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. [8,31], numerical instabilities arising when two consecutive bonds align are naturally eliminated by choosing the torsional potential (4), without the need of imposing rigid constraints on the bending angles. In the following, temperature T , time t, distance, wave vector q, and monomer density ρ are given respectively in units of ǫ/k B (with k B the Boltzmann constant), σ(m/ǫ) 1/2 (with m the monomer mass), σ, σ −1 , and σ −3 . We investigate, at fixed monomer density ρ = 1.0, the temperature dependence of the dynamics for different values of the bending and torsion strength, (K B ,K T ) = (0,0), (4,0.1), (8,0.2), (15,0.5), (25,1), (25,4), and (35,4). In the following, all the data presented in the figures and discussed in the main text will correspond to ρ = 1.0. This value will not be, in general, explicitly mentioned there. We have also studied the case (K B ,K T ) = (35,4) at density ρ = 0.93. The specific information of this case is given in Table I (see below). We investigate typically 8-10 different temperatures for each set of values (K B ,K T ).
We simulate 300 chains, each chain consisting of N = 10 monomers of mass m = 1, placed in a cubic simulation box of lenght L box = 14.4225 for ρ = 1.0, or L box = 14.7756 for ρ = 0.93, with periodic boundary conditions. Equations of motion are integrated by using the velocity Verlet scheme [32]. Computational expense is reduced by implementing a linked-cell method [32]. We use a time step ranging from 10 −4 to 5 × 10 −3 . We take shorter and longer steps for respectively higher and lower values of temperatures, bending and torsional constants. The system is prepared by placing and growing the chains randomly in the simulation box, with a constraint avoiding monomer core overlap. The initial monomer density is ρ = 0.375. Equilibration consists of a first run where the box is rescaled periodically by a factor 0.99 < f < 1 until the target density ρ is reached, and a second isochoric run at that ρ. Thermalization at the target T is achieved by periodic velocity rescaling. After reaching equilibrium, energy, pressure, chain radii of gyration, and end-to-end distances show no drift. Likewise, dynamic correlators show no aging effects. Once the system is equilibrated, a microcanonical run is performed for production of configurations, from which static and dynamic correlators are computed. Static correlators presented here are averaged over typically 300 equispaced configurations. Dynamic correlators are averaged over typically 40 equispaced time origins. The typical duration of a production run is of 40-200 million time steps for respectively high and low temperatures.

III. STATIC PROPERTIES a)Orientational correlations
Simulation results presented in this work 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 fullyflexible chains, (K B ,K T ) = (0,0), and for representative stiff 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. 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. Note that for small r cm differences in the represented data for different barrier strength may even be statistical artifacts, arising from

b)Static structure factors and chain form factors
Now we present results for static structure factors and chain form factors, both for fully-flexible chains and for a representative case of stiff chains. Let us consider an isotropic homogeneous system of volume V containing n identical chains of N monomers. The densities of chains and monomers are respectively denoted by ρ c = n/V and ρ = nN/V . Let us denote the location of a monomer along its chain by the index 1 ≤ a ≤ N . The site-site static structure factor for monomers of indices a and b is defined as: Brackets denote ensemble average. The monomer density distribution for wave vector q is given by In this expresion r a j is the position vector of the ath monomer in the jth chain (1 ≤ j ≤ n). The quantity S ab (q) can be splitted into intrachain and interchain a-b correlations: S ab (q) = ω ab (q) + ρ c h ab (q), (7) or in matrix form, S(q) = w(q) + ρ c h(q). In Eq. (7) ω ab (q) and h ab (q) respectively denote the intrachain and interchain correlations between monomers of type a and b. By averaging over all the possible pairs (a, b) we obtain the static correlators S(q), ω(q) and h(q), which are related through: In this expression S(q) is the total static structure factor, which equivalently can be obtained as S(q) = (nN ) −1 ρ(−q, 0)ρ(q, 0) , where ρ(q) = N a=1 ρ a (q) is the total monomer density distribution. In Eq. (8) the chain form factor, ω(q), accounts for all the static intrachain correlations, while h(q) accounts for all the static interchain correlations. In both cases, no signature of crystallization is present. Indeed no sharp Bragg peaks are observed. In both cases S(q) shows a maximum at q max ≈ 7.0. Since S(q max ) comes from the packing in the first shell around a monomer, the latter corresponds to a typical distance 2π/7.0 ≈ 0.90 in the real space between neighboring monomers. On cooling, the peak at q max ≈ 7.0 increases in intensity, which is a signature of increasing short-range order.
In Fig. 3 we show, for the former values of (K B , K T ), the corresponding results for the form factors ω(q). We note that in the case of fully-flexible chains the form factor is nearly independent on temperature. The form factor for stiff chains exhibits a certain T -dependence, which is however rather weak in comparison with that of S(q).
The T -dependence of ω(q) becomes more clear at low qvalues. The way the form factor behaves on lowering the temperature is directly connected with the values of the mean chain end-to-end radius R ee . Thus, by decreasing temperature from T = 2.0 to T = 0.96, the computed R ee increases from 4.8 to 5.5 for the selected stiff chains. This leads, for lower T , to a stronger decay in ω(q) at low-q. On the other hand, the value R ee = 3.6 for the fully-flexible chains is almost T -independent, leading to a negligible T -dependence of ω(q).

c)Test of the PRISM approximations
The MCT for polymer melts developed by Chong and co-workers [24,25] invokes several approximations of the PRISM theory [27]. In this subsection we summarize such approximations and test their validity for all the investigated range of barrier strength. The site-site direct  correlation function, c ab (q), is introduced via the generalized Ornstein-Zernike relation for polyatomic molecules, or 'reference interaction site model' (RISM) [33], in which intramolecular contributions are accounted by the form factor terms ω ab (q). By inserting (7) in Eq. (9), c ab (q) is related to S ab (q) and ω ab (q) as: Here ω −1 ab (q) and S −1 ab (q) are the elements of, respectively, the matrices w −1 (q) and S −1 (q), which are defined as the inverses of w(q) and S(q).
In the equivalent-site approximation (which is exact for polymer rings) of the PRISM, chain end effects are neglected and all sites are treated equivalently for interchain correlations. Thus, c ab is replaced by the average over all (a, b)-pairs: By introducing this approximation in Eq. (9) and averaging over all (a, b)-pairs we find h(q) = ω(q)c(q)[ω(q) + ρh(q)]. By introducing Eq. (8) in the latter expression we arrive to the scalar equation also known as PRISM equation. In Figs. 4 and 5 we test the validity of the equivalentsite approximation c ab (q) ≈ c(q). We calculate c ab (q) and c(q) respectively through Eqs. (10) and (12), by using the quantities ω −1 ab (q), S −1 ab (q), ω(q), and S(q) as computed from the simulations. Fig. 4 shows results for the fully-flexible case. Data for the case (K B , K T ) = (25, 1) are displayed in Fig. 5. Both data sets correspond to the respective lowest investigated temperatures. We use a representation analogous to that of Ref. [34]. Thus, top and bottom panels in both figures show the comparison of the averaged c(q) with respectively the matrix elements c aa (q) and c a5 (q). The data of Fig. 4 are consistent with results of Ref. [34] for a similar fully-flexible bead-spring model. Data in Fig. 5 constitute new results for the case of implemented intramolecular barriers. By looking at both figures we conclude that the quality of the equivalent-site approximation is not altered by the introduction of strong intramolecular barriers. Data in Fig. 5 display the same trends as in the fully-flexible case. Thus, c ab (q) ≈ c(q) is an excellent approximation except for correlations involving chain end monomers a = 1 (and a = N by symmetry). The latter show deviations from c(q) which are moderate around q max , this q-range being the dominating one in the MCT kernel.
An additional approximation of the PRISM is the ring approximation (which is again exact for polymer rings). First we define the quantitiesS a (q) = N b=1 S ab (q) and S −1 a (q) = N b=1 S −1 ab (q). By exploiting the fact that for a ring polymerS a (q) is a-independent, i.e.,S a (q) ≈ N −1 N a=1S a (q), we find From the definition of fully-flexible and stiff chains. Only for the end monomers a = 1 (and a = N by symmetry) significant differences between S(q) and 1/S −1 a (q) are observed around the wave vector q max .
With all these results we conclude that the approximations assumed by the PRISM theory and introduced in the MCT equations for polymer melts (see below) are fulfilled, with the same quality for all the investigated range of barrier strength.

IV. DYNAMIC PROPERTIES
In this section we show some phenomenological dynamic features induced by the introduction of intramolecular barriers in our model. Panels in Fig. 7 show the Tdependence of the monomer mean squared displacement (MSD) for fully-flexible and representative stiff chains with (K B , K T ) = (25, 1). We observe similar features in both cases, but also some differences. After the initial ballistic regime, a plateau extends over longer times with decreasing temperature. This plateau corresponds to the caging regime -i.e., the temporary trapping of each monomer in the shell of neighboring monomers around it -which is usually observed when approaching a liquidglass transition. At longer times, leaving the plateau, a crossover to a Rouse-like sublinear regime (∆r) 2 ∝ t 0.65 [30,35] is observed for the fully-flexible case. The final crossover to linear diffusion (∆r) 2 ∝ t is reached at long times only for the highest investigated temperatures. However, for the case of stiff chains it is difficult to discriminate power-law behavior over significant time windows. Apparently, the linear diffusive regime is not reached within the simulation time window. Fig. 8 shows the monomer MSD, for fixed values of density ρ = 1.0 and temperature T = 1.5, as a function of the barrier strength. Consistently with results in Ref. [31], we observe that increasing the strenght of the inter-  nal barriers at fixed ρ and T leads to slower dynamics. Fig. 9 shows simulation results at several temperatures, both for fully-flexible and stiff chains, for the normalized density-density correlator f (q, t). The latter is defined as f (q, t) = ρ(−q, 0)ρ(q, t) / ρ(−q, 0)ρ(q, 0) . In both cases the correlator is evaluated at the maximum, q max ≈ 7, of the static structure factor S(q). As in the case of the MSD, both the fully-flexible and stiff cases exhibit the standard behavior in the proximity of a glass transition [30,35]. After the initial transient regime, f (q, t) shows a first decay to a plateau connected with the caging regime. On lowering the temperature this plateau extends over longer time intervals. At long times, a second decay is observed from the plateau to zero. This second decay corresponds to the structural α-relaxation.
Let us define the relaxation time as a time scale probing the α-structural relaxation. This can be done by introducing the time τ x for which the correlator for q max takes the value f (q max , τ x ) = x, provided x is small in comparison with the plateau height. Here we use x = 0.2. Fig. 10 shows τ 0.2 as a funcion of T , for different values of the bending and torsional constants. As observed in the analysis of the mean squared displacements, increasing the chain stiffness slows down the dynamics. At fixed temperature, the relaxation time for the stiffest investigated chains increases up to three decades with respect to the fully-flexible case.
In this section we have demonstrated a main dynamic feature: the slowing down of the dynamics, at fixed density and temperature, by progressively increasing the strength of the intramolecular barriers. This feature strongly suggests that intramolecular barriers constitute and additional mechanism for dynamic arrest, coexisting with the general packing effects induced by density and temperature. In the following we summarize the main predictions of the Mode Coupling Theory and discuss simulation dynamic features within this theoretical framework.

V. MODE COUPLING THEORY: SUMMARY
In this section we briefly summarize universal dynamic scaling laws concerning the MCT liquid-glass dynamics, and test them in the simulated polymer melt for all the investigated range of barrier strength. Extensive reviews on MCT can be found, e.g., in Refs. [14,15,16,17,18,19,36,37]. Though initially derived for simple hard-sphere systems, these predictions follow as consequences of the mathematical structure of the MCT equations. More specifically, they are associated to the bilinear dependence of the memory kernel on the density correlators (see below). Thus, MCT predicts the same dynamic scaling laws of the monoatomic case if such a mathematical structure is retained in systems of polyatomic molecules. This is indeed the case of the MCT for polymer melts developed by Chong and coworkers [24,25] (see below). Therefore, the phenomenological analysis of our simulation results in terms of MCT dynamic scaling laws is justified within the theory.
By starting from the fundamental Liouville equation of motion and using the Mori-Zwanzig projection operator formalism one arrives to an integro-differential equation for the normalized density-density correlator: This equation is obtained by using projectors over the subspace spanned by the densities and the longitudinal currents. The memory kernel , where the quantities R q are, within the Mori-Zwanzig formalism, the associated fluctuating forces. Since the kernel cannot be exactly expressed in terms of f (q, t) and/or its time derivatives, Eq. (15) is not solvable. MCT introduces several approximations for the memory kernel, in order to provide a closed solvable form of Eq. (15). These approximations are: i) It is assumed that the long-time, slow dynamic regime of any observable coupled to density fluctuations can be expressed as a linear combination of 'mode pairs', ρ k ρ q−k . Since the exact expression of the correlator of the fluctuating forces contains a slow contribution which is a linear combination of mode pairs (see e.g., Ref. [17] for details), the former assumption is equivalent to neglecting the fast contribution of the fluctuating forces. In other words, it is equivalent to assuming a large separation between the time scales of the former contributions.
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, where the superscript Q denotes evolution with projected dynamics (see e.g., Ref. [17] for details), and F (q, t) = ρ(−q, 0)ρ(q, t) are just the unnormalized density-density correlators.
By making use of these three approximations, the memory kernel m(q, t) becomes a bilinear form in f (q, t), where the vertex V(q − k, k) is given by: In a monoatomic fluid the direct correlation function c(q) is related to the static structure factor via the exact Ornstein-Zernike relation [38] ρc(q) = 1 − S −1 (q). With all this, Eq. (15) has been reduced to 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). For the case of systems with molecular architecture, Chong and Hirata have obtained [39], by using projectors over site-densities and site-currents, generalized MCT equations of motion for site-site correlators F ab (q, t). The latter are defined as F ab (q, t) = n −1 ρ a (−q, 0)ρ b (q, t) . Note that F ab (q, 0) = S ab (q). The corresponding MCT equations of motion arë where the memory kernel is now given by with p = q − k. By comparing Eqs. (20,21) with Eqs. (15,18,19) we note that the general mathematical structure of the kernel (bilinear in site-site correlators), and of the MCT equations of motion is retained.
Except for very small values of N , numerical solution of the MCT equations for site-site correlators is extremely expensive, and further simplifications are needed in order to obtain a tractable set of equations. For the case of simple bead-spring chains, Chong and co-workers have reduced [24,25] Eqs. (20,21) to a scalar form for f (q, t). This is achieved by introducing in Eqs. (20,21) the equivalent site, Eq. (11), and ring, Eqs. (13,14), approximations of the PRISM theory. The so-obtained scalar MCT equations of motion, memory kernel and vertex for polymer chains are formally identical to Eqs. (15,18,19). The polymer character of the system only enters implicitly through the PRISM relation ρc(q) = 1/ω(q) − 1/S(q), which differs from the Ornstein-Zernike equation, ρc(q) = 1 − S −1 (q), for monoatomic systems. With this, general MCT predictions which originate from the mathematical structure of Eqs. (15,18,19) will be, due to the mentioned formal equivalence, analogous both for monoatomic systems and for polymer chains. Now we summarize such general predictions.
In MCT, nonergodic arrested states (glasses) are defined as those for which density correlators do not exhibit full relaxation. More specifically, if we introduce the nonergodicity parameters, defined as f q = lim t→∞ f (q, t), MCT discriminates between fluid states (f q = 0) and glassy states (f q > 0). At the MCT critical temperature T c , the nonergodicity parameters jump from zero to nonzero values [40]. In the following we use the notation f c q for referring to the critical nonergodicity parameters, i.e., the values of f q at T = T c .
By Laplace transform (t → z) of Eqs. (15,18) and taking the limit z → 0, one finds a coupled set of equations for the nonergodicity parameters: where F q ({f }) denotes a functional, whose explict expression is given in the right-hand side of the equation. Note that Eq. (22) always has the trivial solution f q = 0. Thus, glassy states take place when solutions f q > 0 also exist. Given a tagged chain (labeled s), the density distribution for the ath monomer of the tagged chain is defined as ρ s a (q) = exp[iq · r a s ]. The site-site intrachain correlator is defined as F s ab (q, t) = ρ s a (−q, 0)ρ s b (q, t) . Note that F s ab (q, 0) = ω ab (q). For the derivation of the MCT equations for F s ab (q, t) we refer to [25]. In this case the reduction to a scalar form is not possible. The corresponding nonergodicity parameters f s ab (q) = lim t→∞ F s ab (q, t) are obtained by solving the N × N -matrix equation [41] with I the identity matrix. The corresponding functional F s ab (q) is given by with the vertex The normalized self-correlator, usually introduced as f s (q, t) = (nN ) −1 n j=1 N a=1 exp[iq · (r a j (t) − r a j (0))] , can be equivalently obtained as f s (q, t) = N −1 N a=1 F s aa (q, t). Likewise, the corresponding nonergodicity parameters, defined as the long-time limit of f s (q, t), can be obtained as f s q = N −1 N a=1 f s aa (q). Thus, the solution of Eq. (23) also provides trivially the nonergodicity parameters for the self-correlator.
The separation parameter, ǫ = (T − T c )/T c , is introduced to quantify the relative distance to the critical temperature T c . We are interested in the behavior of f (q, t) in the ergodic fluid by approaching T c from above. Thus we express the long-time behavior of the density-density correlators as: where g q (t) quantifies (small) deviations around f c q for |ǫ| → 0. By introducing Eq. (26) in Eqs. (15,18), expanding the functional F q of Eq. (22) in a power series of |ǫ|, comparing the so-obtained resulting expressions and retaining the lower-order terms (see, e.g., Ref. [36] for a detailed exposition), one finds that g q (t) = h q G(t), where h q only depends on q, and G(t) is a q-independent term which contains the full time dependence of the deviations of f (q, t) around f c q . Thus, we rewrite Eq. (26) as: This expression is known as the first universality of the MCT or factorization theorem. It predicts a scaling function G(t) (known as the β-correlator) that is common for all the density correlators (since it is q-independent). Following the procedure mentioned in the previous paragraph [36], the function G(t) is found to obey the equation: whereG(z) and L[G 2 (t)] are the Laplace transform of respectively G(t) and G 2 (t). In this equation σ = c|ǫ|, with c a constant (see [36] for its explicit expression), and λ is another constant given by The quantities e q and e T q are respectively the eigenvectors of the so-called stability matrix C c (see below) and its traspose, with the normalization conditions q e T q e q = 1 and q e T q (1 − f c q )e 2 q = 1. The elements of the stability matrix are given by The terms C c (q, k, |q − k|) in Eq. (29) are given by: Eq. (28) for the β-correlator does not have an analytical solution. Still, asymptotic expressions can be obtained for different time windows. With this idea in mind the β-time scale is first defined as with t 0 a microscopic time scale and a an exponent. The β-correlator is then rewritten as G(t) = |σ| 1/2 g σ (t/τ β ). By introducing this expression in Eq. (28) and taking the limits t ≪ τ β and t ≫ τ β one finds the asymptotic solutions [16]: where B is a constant [36]. The exponents a and b follow the constraint where Γ denotes the Euler's Gamma function. According to Eqs. (33,34), one finds the asymptotic expressions for Eq. (27): The analysis of the long-time decay usually includes higher-order corrections [36] to Eq. (37): The latter is also known as the von Schweidler expansion. In these equations τ α is the α-time scale, defined as: The exponent γ follows the constraint: Another important prediction of the MCT for states approaching T c from above, is the second universality or time-temperature superposition principle (TTSP). This prediction arises as a long-time scaling property of the MCT equations of motion [16]. According to the TTSP, the long-time decay of any correlator f (q, t) (i.e. the final part of the α-relaxation) is invariant under scaling by the α-relaxation time τ α . In other words, for two temperatures T 1 and T 2 above T c one finds (41) where f (q,t) is a T -independent master function of the normalized timet. While G(t) is common to all correlators, the master function f (q,t) associated to the TTSP is different for each correlator f (q, t). The superposition principle implies that the estimated α-relaxation time, defined in this work as the time τ x where f (q max , t) takes a value x well below the plateau, is proportional to τ α . Thus, it also follows the asymptotic power law The α-decay from the plateau to zero is often well described by an empirical Kohlrausch-Williams-Watts (KWW) function, with A q , β q < 1. Note that the latter does not come out as an analytical solution of the MCT equations. However in the limit q → ∞ of the KWW time τ K q , MCT predicts that [42] where b is the von Schweidler exponent introduced above. The set of equations exposed in this section constitute a series of universal results which originate from the structure of the MCT equations of motion, Eqs. (15,18,19). As mentioned above, the latter were initially derived for simple hard-sphere systems, but the corresponding ones for polymer melts become formally identical following the derivation by Chong and co-workers [24,25]. With this, the scaling laws exposed in this section will also hold in the MCT for polymer melts. Thus, the phenomenological analysis of our simulation data in terms of such scaling laws is justified within the framework of MCT. This analysis is presented in the next section.

VI. MCT ANALYSIS OF SIMULATIONS
In order to test the factorization theorem, Eq. (27), we compute the ratio: where t ′ and t ′′ are arbitrary times in the β-regime.
The ratio for the self-correlators, R s q (t), is defined analogously. If the factorization theorem, and then also the right-hand side of Eq. (45), is fulfilled, the ratios R q (t) and R s q (t) do not depend on the specific correlator. Fig.  11 shows R q (t) and R s q (t) over a broad range of wave vectors 2.3 ≤ q ≤ 16.5. The data correspond to barrier strength (K B , K T ) = (15, 0.5) at T = 0.80. The fixed times t ′′ = 0.8 and t ′ = 100 roughly correspond to the beginning and the end of the plateau regime. There is an intermediate time window of about two decades where the data for density-density and self-correlators collapse onto a q-independent master curve, while they split at both short and late times. Fig. 12 demonstrates that the master curve is, moreover, the same for both density-density and self-correlators. Thus, Figs. 11 and 12 demonstrate the validity of the MCT first universality. Fig. 13 shows a test of the TTSP, Eq. (41), for the density-density correlator evaluated at q max (maximun of the static structure factor S(q)). The data correspond to the case (K B , K T ) = (15, 0.5) and cover a broad temperature range 0.80 ≤ T ≤ 1.5. Data collapse onto a master curve after rescaling the absolute time by the relaxation time τ 0.2 . Thus, the MCT second universality also holds for chains with strong intramolecular barriers.
Solving numerically the MCT equations and determining the dynamic exponents (a, b, γ, λ) is in general a difficult task. When numerical solutions are not available, nonergodicity parameters, prefactors and exponents in Eqs. (36,38,42,44) can be obtained as fit parameters from simulation or experimental data (see, e.g., Refs. [15,35,43,44,45]). Consistency of the analysis requires that dynamic correlators and relaxation times are described by a common set of exponents, all of them related to a single λ-parameter through Eqs. (35,40).
We have performed this consistency test for all the investigated range of barrier strength. The following figures in this section illustrate, for some representative cases, the analysis of simulation data in terms of MCT asymptotic laws. Fig. 14 shows for a broad q-range (2.0 ≤ q ≤ 14.4), fits to the von Schweidler expansion, Eq. (38) (up to second-order terms). Data correspond to density-density correlators f (q, t) for the state point (K B , K T ) = (15, 0.5), T = 0.80 (labelled S1), and to  self-correlators f s (q, t) for (K B , K T ) = (35, 4), T = 1.33 (labelled S2). A good description of the simulation data is achieved, for all the range of q-values and over almost four time decades, with a fixed b-exponent (b = 0.50 and 0.37 for respectively S1 and S2).   rier strength, the q-dependence of the so-obtained critical nonergodicity parameters (f c q for f (q, t) and f sc q for f s (q, t)). For comparison, we also include the fullyflexible case (K B , K T ) = (0, 0). As deduced from the stronger decay of f c q and f sc q for stronger barriers, the introduction of chain stiffness yields a weaker stability of density fluctuations. It also induces a weaker localization for self-motions at fixed density. Thus, by making an approximate fit of f sc q to Gaussian behavior, f sc q ≈ exp(−q 2 l 2 c /6), we estimate, at fixed ρ = 1.0, a localization length l c = 0.19, 0.21, and 0.23 for respectively (K B , K T ) = (0,0), (15,0.5), and (35,4).
Data of self-correlators from the plateau to the limit of the simulation window have been fitted to KWW functions, Eq. (43) (not shown). Fig. 16 shows the qdependence of the so-obtained KWW relaxation times τ K q for the former values of the barrier strength, at their respective lowest investigated temperatures. The lines represent tests of the MCT prediction τ K q ∝ q −1/b for large q. A good description of the data is obtained with the same b-exponents used for the independently obtained von Schweidler fits of Fig. 14. Fig. 17 shows, for the same values of (K B , K T ) in Fig.  16, a test of the power law τ 0.2 ∝ (T − T c ) −γ for the temperature dependence of the estimated α-relaxation times. The fit covers about three time decades. By representing the data in terms of the separation parameter T /T c − 1, clearly different γ-exponents are evidenced for different   Table I displays the results for the so-obtained λ-exponents and critical temperatures T c as a function of (K B , K T ). We also include the corresponding value of the mean endto-end radius (computed at T c ), which provides a qualitative characterization of chain stiffness. From the numerical values in Table I a [36] where dynamic arrest is driven by packing effects. The largest ones, λ 0.9, are similar to those observed in realistic models of polymer melts which incorporate the chemical structure of the chains. Some examples include poly(vinyl methylether) [46], polybutadiene [47] or poly(vinyl ethylene) [48], with respective values of λ = 0.87, 0.93, and 0.93. Thus, the analysis presented here rationalizes the difference in the MCT exponents between fully-flexible bead-spring models and real polymers. The systematic study performed by tuning the barrier strength suggests that large λ-exponents in real polymers arise from the interplay between two distinct mechanisms for dynamic arrest. These are general packing effects and polymerspecific intramolecular barriers. Large λ-values arising from the interplay between distinct arrest mechanisms have been observed in systems of very different nature, as short-ranged attractive colloids [22,49,50] (competition between hard-sphere repulsion and short-ranged reversible bonding), polymer blends [51,52] and colloidal mixtures with strong dynamic asymmetry [53,54] (bulklike caging and matrix-induced confinement), or densified silica [55] (presumably bonding and packing). Numerical solutions of the MCT equations in short-ranged attractive colloids [22,49] and quenched-annealed mixtures [23] have revealed the existence of higher-order MCT transitions, which are characterized by the upper limit λ = 1. Whether higher-order MCT transitions are present at some region of the control parameter space of the investigated model is an open question.
In this section we have performed a phenomenological analysis of the simulation data within the framework of MCT. In the next section the observed trends are compared with numerical solutions of the MCT equations.

VII. SOLUTION OF THE MCT EQUATIONS
We have solved Eqs. (22) and (23) for the nonergodicity parameters, for all the investigated range of barrier strength. In analogy with the procedure exposed in, e.g., Refs. [36,37], the integrals over the reciprocal space in the corresponding MCT functionals of Eqs. (22) and (24) are discretised to a grid of M = 600 equispaced points, with q-spacing ∆q = 0.1, leading to the expressions: and In these expressions the wavevectors are defined as q = x q ∆q, k = x k ∆q, and p = x p ∆q, with x q , x k , x p = 1/2, 3/2, ...1199/2. The prime at the sums over x p means that the latter are restricted to x p -values following the condition |x q − x k | + 1/2 ≤ x p ≤ x q + x k − 1/2. The solutions of Eq. (46) are found by a standard it- , with j the iteration step, and with the initial condition f 0 q = 1. It can be demonstrated that the stability matrix in Eq. (30) has always a maximum non-degenerate eigenvalue E ≤ 1, which takes the upper value E c = 1 at the critical point [36]. Thus, by following the drift of E with changing temperature it is possible to bracket the values of the critical nonergodicity parameters f c q , and the critical temperature T c , with very high precision. Once the values of f c q are obtained, they are fixed in the functional of Eq. (47), and a small number of iterations is needed to find the corresponding critical values f sc ab (q). Finally, the critical nonergodicity parameters for self-correlations are obtained as f sc q = N −1 N a=1 f sc aa (q). Following the procedure exposed above, we solved Eq. (46) by inserting as external inputs the structural quantities, S(q) and c(q), as directly computed from the simulations. However, as previously reported in Ref. [25] for fully-flexible chains, a MCT transition was not observed for any of the investigated barrier strength. This means that the theoretical critical temperature T c is below the lowest simulation temperature for which equilibration was possible. This result is different from the usual observation in non-polymeric systems, for which the theoretical critical point is accessible in simulation time scales. The reason of this difference is, in some way, related with the unability to crystallize of beadspring models, which avoids a fast growing of peaks under cooling in the static structure factor S(q), leading to MCT kernels which are not sufficiently strong to provide nonzero solutions of f q .
Since static correlations computed from our equilibrium simulations do not induce a MCT transition, we are forced to use a structural theory for estimating S(q) and c(q) at lower temperatures, which will allows us to insert them in the MCT equations and to search for the critical temperature. Thus, we solve numerically the PRISM equation with the Percus-Yevick (PY) closure relation [38] for the non-bonded potential V (r) of Eq. (1). The PY relation is given by: where c(r) and h(r) are the Fourier transforms in the real space of c(q) and h(q). The coupled set of nonlinear equations (48,49) is solved by a standard Picard iteration method [56] for the quantity Γ(r) = h(r) − c(r), which is a smooth function over all the range of r. The form factor ω(q) is an external input in this procedure. We observed (see above) that ω(q) exhibits a very weak temperature dependence in comparison to the total static structure factor S(q). Thus we just use for each barrier strength the ω(q), as computed from the simulations, at the lowest temperature for which equilibration was possible.
In Fig. 18 we show a comparison of the critical nonergodicity parameters f c q (top panel) and f sc q (bottom panel) as obtained from numerical solution of the MCT equations, with the results of the fitting procedure of simulation data (see above). The theoretical results qualitatively reproduce the simulation trends, and in particular the observation that at fixed density the intramolecular barriers induce a weaker localization length. Quantitatively, the MCT solutions overstimate the amplitude of the nonergodicity parameters, except in the low-q region of f c q , for which MCT clearly understimates the results. In figure 19 we show a representation of the critical temperature T c as a function of the end-to-end radius R c ee , which quantifies chain stiffness. Values of T c obtained from the phenomenological analysis of the simulations (T MD c ) and from the numerical solutions of the MCT equations (T MCT c ) are compared. We note that T MD c seems to grow monotonously with chain stiffness. This trend is well reproduced by the theory for low and moderate values of the internal barriers. Thus, for values of bending and torsional constants K B < 15 and K T < 0.5, the dependence of T MCT c on R c ee roughly displays the same slope as for T MD c , with a shift factor T MD c /T MCT c ≈ 1.25. Similar shifts between simulation and theory, which have their origin in the mean-field character of the MCT, are observed in other systems [43,57,58]. The range of barrier strength for which T MCT c and T MD c are roughly parallel is significant. Note that for (K B , K T ) = (8, 0.2) the end-to-end radius R c ee is a factor 1.3 longer than for fully-flexible chains.
By further increasing chain stiffness the differences between T MD c and T MCT c progressively increase. We observe a saturation of the theoretical T MCT c around ≈ 0.55, while the simulation T MD c grows up to a value of 1.23 for the stiffest investigated chains. Thus, the agreement between theory and simulation clearly breaks for stiff chains.
Finally, we have computed the corresponding theoretical λ-exponents according to the definitions of Eqs. (29,30,31). We find an almost constant value of λ ≈ 0.72 for all the investigated range of barrier strength. This result is clearly different from the observations in the phenomenological MCT analysis of simulation data (see data in Table I), which provides a strong dependence of λ on the barrier strength.

VIII. DISCUSSION
The results reported in the previous section show that, though reproducing some qualitative simulation trends for low and moderate barriers, numerical solutions of the MCT equations exhibit important differences with simulation values as the limit of stiff chains is approached. Another result to be understood is the clear disagreement between the almost constant value of the theoretical λ-exponent, and the observed strong dependence of simulation values on the barrier strength.
The observed disagreement between theory and simulation for strong barriers does not seem to be related with the failure of the PRISM approximations for stiff chains, which have been introduced in the derivation of the MCT equations. Indeed we have shown that the quality of the used PRISM approximations is the same for fully-flexible and stiff chains (Figs. 4, 5 and 6). Having said this, it might be argued that the theory is simply wrong: the phenomenological MCT analysis is apparently successful, but one finds that it has little to do with the theory, for stiff chains, when solving the MCT equations. However, we remind that the phenomenological analysis has shown, for all the investigated range of barrier strength: i) the validity of the two MCT universalities, i.e., the factorization theorem (Figs. 11 and 12) and the TTSP (Fig. 13) ii) the possibility of a good description of different dynamic observables (Figs. 14, 16 and 17) with a set of dynamic exponents which are consistently transformed, through Eqs. (35,40), to a single λ-exponent.
We believe that all these observations, for all the investigated cases, are not fortuitous. At this point it must be noted that the predictions referred to in points i) and ii) arise, within MCT, as a consequence of the mathematical structure of the equations of motion, more precisely they originate from the bilinear form of the memory kernel. The specific values of the numerical solutions clearly depend on the coefficients of the bilinear products (which enter through the vertices of the kernel), but the factorization theorem, the TTSP, and the asymptotic scaling laws are universal properties provided the kernel is bilinear. Thus, the results of the phenomenological analysis suggest that the underlying physics may be connected to a bilinear memory kernel, though for high barriers the actual coefficients strongly differ from those introduced by MCT through the vertices, thus leading to theoretical results which strongly differ from simulations.
In other words, the present results suggest that there may be relevant static contributions for the case of stiff chains which are missing in the MCT vertices. Thus, the inclusion of such contributions will increase the strength of the kernel and will induce the theoretical transition at higher values of T c , which might improve the comparison between T MCT c and T MD c of Fig. 19. Recalling the three main approximations of MCT, we suggest that the convolution approximation, Eq. (16), might break for stiff chains. Though possibly it is not the case for intermolecular contributions, its breakdown for intramolecular contributions in stiff chains is plausible. It is known that the convolution approximation fails when static correlations show a strong directionality at near-neighbor distances, as for e.g., network-forming liquids as silica [59]. This directionality is clearly enhanced for intrachain correlations by increasing the barrier strength, as evidenced by the progressively larger values of the end-to-end radius (see Table I). For the case of silica, it has been shown that the explicit inclusion of three-point static correlations in the MCT vertex improves significantly the quality of the comparison between theory and simulations [59]. A similar improvement might be achieved in the present case by similarly incorporating the intrachain three-point static contributions. Work in this direction is in progress.

IX. CONCLUSIONS
We have performed simulations on a simple beadspring model for polymer melts with intramolecular barriers. The role of such barriers on the glass transition has been investigated by systematically tuning the barrier strength. Dynamic correlators probing the structural relaxation have been analyzed in the framework of the Mode Coupling Theory. We have obtained critical nonergodicity parameters, critical temperatures and dynamic exponents of the theory from consistent fits of simula-tion data to MCT asymptotic laws. From the analysis of the critical nonergodicity parameters we deduce that the presence of the barriers induces a weaker localization length in the system at fixed density. The increase of the barrier strength at fixed density also induces a higher critical temperature T c . The values of the dynamic exponents, as obtained from the phenomenological analysis of the simulation data, exhibit significant differences between the limit of fully-flexible and stiff chains. In particular the so-called λ-exponent takes standard values λ ∼ 0.7 for the fully-flexible case and values approaching the upper limit λ = 1 for strong intramolecular barriers. While the former λ-values are characteristic of simple systems dominated by packing effects, transitions with λ ≈ 1 arise in systems with different competing mechanisms for dynamic arrest. In our systems these large λvalues suggest a competition between two distinct mechanisms: general packing effects and polymer-specific intramolecular barriers.
For a comparison between simulation and theory, we have numerically solved the MCT equations, following a recent extension of the MCT by Chong and co-workers for polymer melts. The approximations assumed by the structural PRISM theory, which are introduced in the MCT equations, are fulfilled for all the investigated values of the barrier strength. We have compared the critical nonergodicity parameters and critical temperatures T c , as obtained by solving the MCT equations, with the corresponding values from the phenomenological analysis of the simulation data. The theoretical calculations qualitatively reproduce the trends observed in the simulations for low and moderate barriers. However strong discrepancies are observed as the limit of high barriers is approached. The reason for such a disagreement possibly lies in the nature of the approximations made in the derivation of the MCT equations. In particular, the convolution approximation for three-point static correlations might be unadequate for stiff chains. We suggest that a reformulation of MCT equations for polymer melts, explicitly including intrachain three-point static correlations, might lead to a better agreement between simulations and theory. Work in this direction is in progress.