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

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

Marco Bernabei sckbernm@ehu.es Donostia International Physics Center, Paseo Manuel de Lardizabal 4, 20018 San Sebastián, Spain.    Angel J. Moreno Centro de Física de Materiales (CSIC, UPV/EHU)-Materials Physics Center, Apartado 1072, 20080 San Sebastián, Spain.    Juan Colmenero Donostia International Physics Center, Paseo Manuel de Lardizabal 4, 20018 San Sebastián, Spain. Centro de Física de Materiales (CSIC, UPV/EHU)-Materials Physics Center, Apartado 1072, 20080 San Sebastián, Spain. Departamento de Física de Materiales, Universidad del País Vasco (UPV/EHU), Apdo. 1072, 20080 San Sebastián, Spain.
July 5, 2019
Abstract

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 -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 -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.

pacs:
64.70.pj, 64.70.qj, 61.20.Ja

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) doi () 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 doi (); degennes (); mcleish ().

Another particular ingredient of polymers is that, apart from fast librations or methyl group rotations methylrev (), 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 meyer (); vettorel (). They also enhance dynamic features which are usually associated to reptation faller (); bulacu1 (), which controls rheological properties mcleish (). 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 bustamante (); kas (); ober (). Moreover, chain stiffness seems to play an important role in the absorption behavior of polymers at interfaces semenov (); ubbink (). 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 mctrev1 (); mctrev2 (); gotzebook (); reichman (); das (); gotzenewbook (). We extend preliminary results reported by us in Ref. bernabei () 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 , characteristic of simple fluids dominated by packing effects. On the contrary, for strong intramolecular barriers the -values approach the upper limit characteristic of higher-order MCT transitions. The latter arise in systems with different coexisting mechanisms for dynamic arrest schematic (); sperl (); krakoviack (). In the system investigated here, the obtained results suggest an interplay between general packing effects and polymer-specific intramolecular barriers.

Chong and co-workers chongprl (); chongpre () have recently presented an extension of the MCT to simple fully-flexible bead-spring models of polymer systems, in the framework of the polymer reference interaction site model (PRISM) schweizer (); curro (); putz (). In this formalism each molecule is divided into interaction sites corresponding to monomers. A key assumption of the PRISM is the replacement of the site-specific 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

(1)

where and . The potential is set to zero beyond the cutoff distance , with . The values and guarantee continuity of potential and forces at . The potential is purely repulsive. It does not show local minima within the interaction range . Thus, it drives dynamic arrest only through packing effects. Along the chain backbone, of monomers, an additional finitely-extensible nonlinear elastic (FENE) potential grest (); bennemann () is used to introduce bonds between consecutive monomers:

(2)

where and . The superposition of potentials (1) and (2) provides an effective bond potential for consecutive monomers with a sharp minimum at , which makes bond crossing impossible.

Intramolecular barriers are implemented by means of a combined bending , and torsional potential . We have used the potentials proposed by Bulacu and van der Giessen in Refs. bulacu1 (); bulacu2 (). 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 by the cosine harmonic bending potential

(3)

where is the bending angle between consecutive monomers , and (with ).

The torsional potential constrains the dihedral angle , which is defined for the consecutive monomers , , and (with ), as the angle between the two planes defined by the sets (, , ) and (, , ). The form of this potential is

(4)

The values of the coefficients are , , , and bulacu1 (); bulacu2 (). The torsional potential depends both on the dihedral angle and on the bending angles and . As noted in Refs. bulacu1 (); bulacu2 (), 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 , time , distance, wave vector , and monomer density are given respectively in units of (with the Boltzmann constant), (with the monomer mass), , , and . We investigate, at fixed monomer density , the temperature dependence of the dynamics for different values of the bending and torsion strength, , (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 . This value will not be, in general, explicitly mentioned there. We have also studied the case , at density . The specific information of this case is given in Table 1 (see below). We investigate typically 8-10 different temperatures for each set of values ,.

We simulate 300 chains, each chain consisting of monomers of mass , placed in a cubic simulation box of lenght for , or for , with periodic boundary conditions. Equations of motion are integrated by using the velocity Verlet scheme frenkel (). Computational expense is reduced by implementing a linked-cell method frenkel (). We use a time step ranging from to . 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 . Equilibration consists of a first run where the box is rescaled periodically by a factor until the target density is reached, and a second isochoric run at that . Thermalization at the target 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 , 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 . This is illustrated in Fig. 1, which shows the time evolution of along a typical simulation window, both for fully-flexible chains, , (0,0), and for representative stiff chains, , (35,4).

Local orientational order is also negligible. This is evidenced by computing a similar correlator . 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 . Fig. 1 displays, for the former cases of fully-flexible and stiff chains, data of for several values of . Negligible values of are obtained for . Thus, the time average over the simulation time window, , provides values . 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 is observed only for very small interchain distances (see data for ). Again, the introduction of chain stiffness does not induce clear changes in the orientational order at this length scale. Note that for small differences in the represented data for different barrier strength may even be statistical artifacts, arising from the small number of neighboring chains within such distances and the limitted time of observation (see the amplitude of the fluctuations in data of Fig. 1).

Figure 1: Time evolution of the global and local orientational parameter (see text) for fully-flexible (top) and stiff chains with , (bottom), for two selected low temperatures.

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 containing identical chains of monomers. The densities of chains and monomers are respectively denoted by and . Let us denote the location of a monomer along its chain by the index . The site-site static structure factor for monomers of indices and is defined as:

(5)

Brackets denote ensemble average. The monomer density distribution for wave vector is given by

(6)

In this expresion is the position vector of the th monomer in the th chain (). The quantity can be splitted into intrachain and interchain - correlations:

(7)

or in matrix form, . In Eq. (7) and respectively denote the intrachain and interchain correlations between monomers of type and . By averaging over all the possible pairs we obtain the static correlators , and , which are related through:

(8)

In this expression is the total static structure factor, which equivalently can be obtained as , where is the total monomer density distribution. In Eq. (8) the chain form factor, , accounts for all the static intrachain correlations, while accounts for all the static interchain correlations.

Fig. 2 (top panel) shows simulation results for as a function of temperature for fully-flexible chains, . Data for representative stiff chains, , are shown in the bottom panel. In both cases, no signature of crystallization is present. Indeed no sharp Bragg peaks are observed. In both cases shows a maximum at . Since comes from the packing in the first shell around a monomer, the latter corresponds to a typical distance in the real space between neighboring monomers. On cooling, the peak at increases in intensity, which is a signature of increasing short-range order.

Figure 2: Temperature dependence of the static structure factor for fully-flexible chains (top panel) and for chains with barrier strength (bottom panel).

In Fig. 3 we show, for the former values of , the corresponding results for the form factors . 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 -dependence, which is however rather weak in comparison with that of . The -dependence of becomes more clear at low -values. The way the form factor behaves on lowering the temperature is directly connected with the values of the mean chain end-to-end radius . Thus, by decreasing temperature from to , the computed increases from 4.8 to 5.5 for the selected stiff chains. This leads, for lower , to a stronger decay in at low-. On the other hand, the value for the fully-flexible chains is almost -independent, leading to a negligible -dependence of .

Figure 3: Temperature dependence of the form factor for fully-flexible chains (lines) and for chains with barrier strength (symbols). For clarity, the inset shows results in the range of low-. Different colors correspond to different temperatures, following the legends of Fig. 2.

c)Test of the PRISM approximations

The MCT for polymer melts developed by Chong and co-workers chongprl (); chongpre () invokes several approximations of the PRISM theory curro (). 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, , is introduced via the generalized Ornstein-Zernike relation for polyatomic molecules, or ‘reference interaction site model’ (RISM) rism (),

(9)

in which intramolecular contributions are accounted by the form factor terms . By inserting (7) in Eq. (9), is related to and as:

(10)

Here and are the elements of, respectively, the matrices and , which are defined as the inverses of and .

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, is replaced by the average over all -pairs:

(11)

By introducing this approximation in Eq. (9) and averaging over all -pairs we find . By introducing Eq. (8) in the latter expression we arrive to the scalar equation

(12)

also known as PRISM equation.

In Figs. 4 and 5 we test the validity of the equivalent-site approximation . We calculate and respectively through Eqs. (10) and (12), by using the quantities , , , and as computed from the simulations. Fig. 4 shows results for the fully-flexible case. Data for the case are displayed in Fig. 5. Both data sets correspond to the respective lowest investigated temperatures. We use a representation analogous to that of Ref. aichelepre (). Thus, top and bottom panels in both figures show the comparison of the averaged with respectively the matrix elements and . The data of Fig. 4 are consistent with results of Ref. aichelepre () 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, is an excellent approximation except for correlations involving chain end monomers (and by symmetry). The latter show deviations from which are moderate around , this -range being the dominating one in the MCT kernel.

Figure 4: Test of the equivalent site approximation, Eq. (11), for fully-flexible chains at . Top and bottom panels compare with respectively matrix elements and . The insets enhance the region around the wave vector for the maximum of the static structure factor .
Figure 5: Test of the equivalent site approximation, Eq. (11), for stiff chains with , at . Top and bottom panels compare with respectively matrix elements and . The insets enhance the region around the wave vector for the maximum of the static structure factor .

An additional approximation of the PRISM is the ring approximation (which is again exact for polymer rings). First we define the quantities and . By exploiting the fact that for a ring polymer is -independent, i.e., , we find

(13)

From the definition of and the relation is exact. By introducing the ring approximation the former relation is transformed into:

(14)

Fig. 6 shows a test of the ring approximation of Eqs. (13) (main panels) and (14) (insets). This is done both for fully-flexible chains (top panel) and for stiff chains with (bottom panel). The comparison between , and as computed from simulations is in general excellent, with the same quality for fully-flexible and stiff chains. Only for the end monomers (and by symmetry) significant differences between and are observed around the wave vector .

Figure 6: Test of the ring approximation, Eqs. (13,14). Top panel: fully-flexible case at . Bottom panel: barrier strength , at . Main panels and insets compare (symbols) with respectively and 1/, for the sites 1,2 and 5 (lines).

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 -dependence of the monomer mean squared displacement (MSD) for fully-flexible and representative stiff chains with . 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 liquid-glass transition. At longer times, leaving the plateau, a crossover to a Rouse-like sublinear regime bennemann (); aichele () is observed for the fully-flexible case. The final crossover to linear diffusion 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.

Figure 7: Temperature dependence of the monomer mean squared displacement for fully-flexible (top) and stiff chains with (bottom). The solid and dashed lines indicate respectively sublinear () and linear behavior.

Fig. 8 shows the monomer MSD, for fixed values of density and temperature , as a function of the barrier strength. Consistently with results in Ref. bulacu2 (), we observe that increasing the strenght of the internal barriers at fixed and leads to slower dynamics.

Figure 8: Monomer mean squared displacement, for several values of the barrier strength, at fixed density and temperature .

Fig. 9 shows simulation results at several temperatures, both for fully-flexible and stiff chains, for the normalized density-density correlator . The latter is defined as . In both cases the correlator is evaluated at the maximum, , of the static structure factor . 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 bennemann (); aichele (). After the initial transient regime, 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.

Figure 9: Temperature dependence of for fully-flexible (top) and stiff chains with (bottom). The wave vector is , in both cases corresponding to the maximum of the static structure factor .

Let us define the relaxation time as a time scale probing the -structural relaxation. This can be done by introducing the time for which the correlator for takes the value , provided is small in comparison with the plateau height. Here we use . Fig. 10 shows as a funcion of , 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.

Figure 10: Temperature dependence of the relaxation times for several values of the barrier strength.

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. mctrev1 (); mctrev2 (); gotzebook (); reichman (); das (); gotzenewbook (); franosch (); fuchspre (). 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 co-workers chongprl (); chongpre () (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:

(15)

This equation is obtained by using projectors over the subspace spanned by the densities and the longitudinal currents. The memory kernel , where the quantities are, within the Mori-Zwanzig formalism, the associated fluctuating forces. Since the kernel cannot be exactly expressed in terms of 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’, . 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. reichman () 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,

(16)

iii) Kawasaki approximation: dynamic four-point correlations are factorized in terms of products of dynamic two-point correlations,

(17)

where the superscript denotes evolution with projected dynamics (see e.g., Ref. reichman () for details), and are just the unnormalized density-density correlators.

By making use of these three approximations, the memory kernel becomes a bilinear form in ,

(18)

where the vertex is given by:

(19)

In a monoatomic fluid the direct correlation function is related to the static structure factor via the exact Ornstein-Zernike relation hansen () . With all this, Eq. (15) has been reduced to a closed set of coupled equations which can be solved self-consistently, provided and are known (the latter are external inputs in the MCT equations).

For the case of systems with molecular architecture, Chong and Hirata have obtained chonghirata (), by using projectors over site-densities and site-currents, generalized MCT equations of motion for site-site correlators . The latter are defined as . Note that . The corresponding MCT equations of motion are

(20)

where the memory kernel is now given by

(21)

with . 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 , 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 chongprl (); chongpre () Eqs. (20,21) to a scalar form for . 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 , which differs from the Ornstein-Zernike equation, , 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 , MCT discriminates between fluid states () and glassy states (). At the MCT critical temperature , the nonergodicity parameters jump from zero to nonzero values hopping (). In the following we use the notation for referring to the critical nonergodicity parameters, i.e., the values of at .

By Laplace transform () of Eqs. (15,18) and taking the limit , one finds a coupled set of equations for the nonergodicity parameters:

(22)

where 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 . Thus, glassy states take place when solutions also exist.

Given a tagged chain (labeled s), the density distribution for the th monomer of the tagged chain is defined as . The site-site intrachain correlator is defined as . Note that . For the derivation of the MCT equations for we refer to chongpre (). In this case the reduction to a scalar form is not possible. The corresponding nonergodicity parameters are obtained by solving the -matrix equation chongdum ()

(23)

with the identity matrix. The corresponding functional is given by

(24)

with the vertex

(25)

The normalized self-correlator, usually introduced as , can be equivalently obtained as . Likewise, the corresponding nonergodicity parameters, defined as the long-time limit of , can be obtained as . Thus, the solution of Eq. (23) also provides trivially the nonergodicity parameters for the self-correlator.

The separation parameter, , is introduced to quantify the relative distance to the critical temperature . We are interested in the behavior of in the ergodic fluid by approaching from above. Thus we express the long-time behavior of the density-density correlators as:

(26)

where quantifies (small) deviations around for . By introducing Eq. (26) in Eqs. (15,18), expanding the functional of Eq. (22) in a power series of , comparing the so-obtained resulting expressions and retaining the lower-order terms (see, e.g., Ref. franosch () for a detailed exposition), one finds that , where only depends on , and is a -independent term which contains the full time dependence of the deviations of around . Thus, we rewrite Eq. (26) as:

(27)

This expression is known as the first universality of the MCT or factorization theorem. It predicts a scaling function (known as the -correlator) that is common for all the density correlators (since it is -independent). Following the procedure mentioned in the previous paragraph franosch (), the function is found to obey the equation:

(28)

where and are the Laplace transform of respectively and . In this equation , with a constant (see franosch () for its explicit expression), and is another constant given by

(29)

The quantities and are respectively the eigenvectors of the so-called stability matrix (see below) and its traspose, with the normalization conditions and . The elements of the stability matrix are given by

(30)

The terms in Eq. (29) are given by:

(31)

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

(32)

with a microscopic time scale and an exponent. The -correlator is then rewritten as . By introducing this expression in Eq. (28) and taking the limits and one finds the asymptotic solutions gotzebook ():

(33)
(34)

where is a constant franosch (). The exponents and follow the constraint

(35)

where denotes the Euler’s Gamma function. According to Eqs. (33,34), one finds the asymptotic expressions for Eq. (27):

(36)
(37)

The analysis of the long-time decay usually includes higher-order corrections franosch () to Eq. (37):

(38)

The latter is also known as the von Schweidler expansion. In these equations is the -time scale, defined as:

(39)

The exponent follows the constraint:

(40)

Another important prediction of the MCT for states approaching 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 gotzebook (). According to the TTSP, the long-time decay of any correlator (i.e. the final part of the -relaxation) is invariant under scaling by the -relaxation time . In other words, for two temperatures and above one finds

(41)

where is a -independent master function of the normalized time . While is common to all correlators, the master function associated to the TTSP is different for each correlator . The superposition principle implies that the estimated -relaxation time, defined in this work as the time where takes a value well below the plateau, is proportional to . Thus, it also follows the asymptotic power law

(42)

The -decay from the plateau to zero is often well described by an empirical Kohlrausch-Williams-Watts (KWW) function,

(43)

with . Note that the latter does not come out as an analytical solution of the MCT equations. However in the limit of the KWW time , MCT predicts that fuchsjcns ()

(44)

where 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 chongprl (); chongpre (). 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:

(45)

where and are arbitrary times in the -regime. The ratio for the self-correlators, , is defined analogously. If the factorization theorem, and then also the right-hand side of Eq. (45), is fulfilled, the ratios and do not depend on the specific correlator. Fig. 11 shows and over a broad range of wave vectors . The data correspond to barrier strength at . The fixed times and 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 -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.

Figure 11: Test of the factorization theorem, Eq. (45), for density-density (top) and self-correlators (bottom) at and . The different curves correspond to equispaced wave vectors in the range . The fixed times are and .
Figure 12: Common representation of (symbols) and (lines), for selected wavevectors (common colors correspond to common -values). As in Fig. 11, data correspond to and , and the selected fixed times are and .
Figure 13: Test of the time-temperature superposition principle, Eq. (41), for the density-density correlators at , for the case .

Fig. 13 shows a test of the TTSP, Eq. (41), for the density-density correlator evaluated at (maximun of the static structure factor ). The data correspond to the case and cover a broad temperature range . Data collapse onto a master curve after rescaling the absolute time by the relaxation time . Thus, the MCT second universality also holds for chains with strong intramolecular barriers.

Figure 14: Symbols: simulations results for density correlators. Top panel: for , at . Bottom panel: for , at . Identical symbols in both panels correspond to identical wave vectors [values are given in panel (a)]. Lines are fits to the von Schweidler expansion, Eq. (38) (up to second-order terms), with (top) and 0.37 (bottom).

Solving numerically the MCT equations and determining the dynamic exponents 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. mctrev2 (); aichele (); vanmegen (); koband2 (); baschrev ()). 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 -range (), fits to the von Schweidler expansion, Eq. (38) (up to second-order terms). Data correspond to density-density correlators for the state point , (labelled S1), and to self-correlators for , (labelled S2). A good description of the simulation data is achieved, for all the range of -values and over almost four time decades, with a fixed -exponent ( and 0.37 for respectively S1 and S2).

Fig. 15 displays, for the former values of the barrier strength, the -dependence of the so-obtained critical nonergodicity parameters ( for and for ). For comparison, we also include the fully-flexible case . As deduced from the stronger decay of and 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 to Gaussian behavior, , we estimate, at fixed , a localization length 0.19, 0.21, and 0.23 for respectively (0,0), (15,0.5), and (35,4).

Figure 15: Critical nonergodicity parameters, as determined from fits to Eq. (38), for different barrier strength. Top and bottom panels show data for respectively and .

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 -dependence of the so-obtained KWW relaxation times for the former values of the barrier strength, at their respective lowest investigated temperatures. The lines represent tests of the MCT prediction for large . A good description of the data is obtained with the same -exponents used for the independently obtained von Schweidler fits of Fig. 14.

Figure 16: Symbols: -dependence of KWW relaxation times for different barrier strength. Lines are fits to (see text). From top to bottom , 0.50 and 0.37.
Figure 17: Symbols: -dependence of relaxation times for different barrier strength. Lines are fits to (see text). From top to bottom , 2.74 and 3.43.

Fig. 17 shows, for the same values of in Fig. 16, a test of the power law 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 , clearly different -exponents are evidenced for different barrier strength. A good description of the data is obtained with the -values derived, through Eqs. (35, 40), from the -values used in Figs. 14 and 16. This result demonstrates the consistency of the MCT analysis for the representative examples showed here, which cover all the range of investigated barrier strength between fully-flexible and stiff chains.

1 0 0 3.6 0.48 0.761
1 4 0.1 4.4 0.54 0.767
1 8 0.2 4.7 0.67 0.773
1 15 0.5 5.2 0.75 0.785
1 25 1 5.5 0.82 0.827
1 25 4 6.4 1.02 0.845
1 35 4 6.5 1.23 0.862
0.93 35 4 6.9 1.02 0.885
Table 1: Values of the MCT -exponents and critical temperatures for different and barrier strength. Also included are the mean chain end-to-end radii at .

Similar consistent tests (not shown) have been performed for the rest of investigated systems. Table 1 displays the results for the so-obtained -exponents and critical temperatures as a function of . We also include the corresponding value of the mean end-to-end radius (computed at ), which provides a qualitative characterization of chain stiffness. From the numerical values in Table 1 a clear correlation between the strength of the internal barriers and the values of and is unambiguously demonstrated. The interplay between packing effects and intramolecular barriers induces a progressive increase of at fixed density. A similar effect is observed for the -exponent, which increases from for fully-flexible chains to for the stiffest investigated chains. The smallest -values in Table 1 are typical of simple glass-formers as the archetype hard-sphere fluid () franosch () 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 chemical structure of the chains. Some examples include poly(vinyl methylether) capponi (), polybutadiene narros1 () or poly(vinyl ethylene) narros2 (), 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 polymer-specific 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 sperl (); dawson (); zaccarelli () (competition between hard-sphere repulsion and short-ranged reversible bonding), polymer blends blends1 (); blends2 () and colloidal mixtures with strong dynamic asymmetry mixtures1 (); mixtures2 () (bulk-like caging and matrix-induced confinement), or densified silica voigtsildens () (presumably bonding and packing). Numerical solutions of the MCT equations in short-ranged attractive colloids sperl (); dawson () and quenched-annealed mixtures krakoviack () have revealed the existence of higher-order MCT transitions, which are characterized by the upper limit . 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. franosch (); fuchspre (), the integrals over the reciprocal space in the corresponding MCT functionals of Eqs. (22) and (24) are discretised to a grid of equispaced points, with -spacing , leading to the expressions: