Self-motion and thea relaxation in a simulated glass-forming polymer: Crossover from Gaussian to non-Gaussian dynamic behavior

J. Colmenero, 1,2,3,* F. Alvarez, 1,2 and A. Arbe Departamento de Fı ́sica de Materiales, Universidad del Paı ́s Vasco UPV/EHU, Apartado 1072, 20080 San Sebastia ́n, Spain Unidad de Fı ́sica de Materiales (CSIC-UPV/EHU), Apartado 1072, 20080 San Sebastia ́n, Spain Donostia International Physics Center, Apartado 1072, 20080 San Sebastia ́n, Sp in ~Received 6 July 2001; revised manuscript received 10 December 2001; published 4 April 2002 !


I. INTRODUCTION
The understanding of how atoms or molecules move within a supercooled liquid and the way this liquid becomes a glass-the glass transition-is still one of the main challenges in the field of condensed matter. Nowadays there is an increasing interest in the so-called ''dynamical heterogeneity'' of the main dynamical process in supercooled liquids: the ␣ relaxation. Computer simulations ͑see ͓1-5͔ as representative references͒ as well as several experimental techniques ͓6͔ have provided different evidences of heterogeneous behavior. Dynamical heterogeneities have also been directly observed in colloidal models of glass-forming systems ͓7͔. Moreover, the heterogeneous dynamics in polymer films of poly͑methylacrylate͒ near their glass-transition temperature has recently been proved by means of singlemolecule spectroscopy ͓8͔. Many different theoretical concepts of dynamical heterogeneity are usually invoked ͓9͔. However, from a general point of view, a system is considered as dynamically heterogeneous if a dynamically distinguishable subensemble ͑e.g., ''fast'' or ''slow'' particles͒ can be isolated by computer simulation or experiment. Computer simulations of model systems ͓3,5͔, as well as the experiments in colloidal systems mentioned above, show that the most mobile particles are spatially correlated and move cooperatively along ''stringlike paths.'' Recent neutron scattering ͑NS͒ results in a fragile glass, Ca 0.4 K 0.6 ͑NO 3 ͒ 1.4 also seem to suggest this interpretation ͓10͔. Dynamical heterogeneity is often discussed in terms of the self-part of the Van Hove correlation function G s (r,t) ͓2-4,7͔. This function gives the probability distribution of the position of an atom at time t relative to its position at tϭ0. In the case of a simple atomic diffusion, a Gaussian distribution is obtained. However, deviations of G s (r,t) from the Gaussian form can be expected in the case of more complex dynamic processes as, for instance, a heterogeneous dynamics. G s (r,t) can easily be evaluated in a computer simulation ''experiment'' where the atomic trajectories are directly recorded. Strong deviations from the Gaussian behavior, which depend on temperature, have been reported for simple model liquids ͑see, for example, ͓2͔͒. Similar results have also been recently reported for colloidal models of supercooled liquids where G s (r,t) can be experimentally calculated ͓7͔. These deviations have been interpreted in terms of a heterogeneous dynamic behavior of the slow atomic motions that give rise to the ␣ relaxation. In the case of real atomic or molecular glass-forming systems, G s (r,t) is not experimentally accessible. The only experimental way to obtain information about the Gaussian behavior is through the momentum transfer Q dependence of the so-called intermediate incoherent scattering function F s (Q,t), which is defined as the Fourier trans-*Email address: wapcolej@sc.ehu.es form of G s (r,t). Because the Q range accessed by NS studies is too small to perform an inverse Fourier transform to space coordinates to derive G s (r,t) unambiguously, the validity of the Gaussian approximation is checked by studying the Q dependence of F s (Q,t). This has been extensively investigated in the ␣-relaxation regime of glass-forming polymer systems ͓11-16͔. NS results of many different polymers ͓11-13͔ show that in the time and Q range where the ␣ relaxation is observed by these techniques ͑typically: 5 ps ՇtՇ2 ns; 0.2ՇQՇ1.5 Å Ϫ1 ͒ F s (Q,t) shows an approximate Gaussian behavior. These results are, in principle, in disagreement with those obtained by other techniques e.g., 4D-NMR or selective photobleaching of probe molecules ͓6͔. However, we have to point out that the time/frequency range usually explored by these techniques-and thereby the temperature range-is rather different from that accessible by neutron scattering. Why does dynamical heterogeneity not seem to be visible for NS in glass-forming polymers? How can results from different techniques be compared? Is there any characteristic time for the dynamical heterogeneity in these systems? And, if so, how does it compare with the ␣-relaxation times? Computer simulations are unique techniques to shed new light onto these questions because, not only G s (r,t) and F s (Q,t), but also other different correlators can be computed from the atomic trajectories. Moreover, the Q range available can be extended beyond the actual possibilities of experimental neutron scattering techniques. However, in spite of these possibilities, the Q dependence of F s (Q,t) has only been explored in a reduced number of molecular dynamics ͑MD͒ simulations studies on Lennard-Jones liquids ͓17͔, and water ͓18͔, and most of the reported works only deal with the temperature dependence of F s (Q,t) at the Q value of the first maximum of the static structure factor S(Q) ͑see, for example, ͓2͔͒. In the case of polymer systems, the Q dependence of F s (Q,t) has been previously investigated for a simple ''bead-spring'' model of a polymer melt ͓19,20͔ and for a ''united atom'' model as well ͓21͔. Although the main focus of these papers was to check the mode-coupling theory, deviations from the Gaussian behavior of F s (Q,t) were also reported in both cases.
With these ideas in mind, we have carried out a fully atomistic molecular dynamics simulation in a realistic model of polyisoprene ͑PI͒. Our results show that there is a crossover from Gaussian to non-Gaussian behavior of the ␣ relaxation, which takes place in the Q range of the first maximum of the static structure factor S(Q). The microscopic origin of such a crossover and its possible relationship with the abovementioned apparent contradiction between NS and relaxation techniques is discussed.

A. Model and simulation
The simulations were carried out by using the INSIGHT ͑INSIGHT II 4.0.0 P version͒ and the Discover-3 module from Molecular Simulations Inc. with the Polymer Consortium Force Field ͓22͔. Most parameters of this field were derived based on ab initio data using a least-squared-fit technique developed by Hagler and co-workers ͓22͔. The functional form includes terms that can be divided into two categories-valence terms including diagonal and offdiagonal cross-coupling terms and nonbonded interactions terms. The valence terms represent internal coordinates of the bond, angle, torsion angle, and out of plane angle, and the cross-coupling terms include combinations of two or three internal coordinates. The cross-coupling terms are important for predicting vibration frequencies and structural variations associated with conformational changes. The analytical expression employs quartic polynomials for bond stretching and angle bending and a three-term Fourier expansion for torsions. The nonbonded interaction terms include a Coulombic function for the electrostatic interaction and a Lennard-Jones 9-6 potential function rather than the more customary Lennard-Jones 12-6 potential for the van der Waals term. More information about this kind of force fields, including the complete analytical expression for the functional form, can be found in ͓22,23͔. The model system was built by means of the well-known amorphous cell protocol, which was proposed for the first time by Theodorou and Suter ͓24͔. In this work, a cubic cell containing one polymer chain of 100 monomer units ͓uCH 2 uCHvC(CH 3 )uCH 2 u͔ 100 was constructed at 363 K and a density (ϭ0.869 g/cm 3 ), which was extrapolated to 363 K from the available experimental data ͓25͔. Such a density leads to a cell dimension of 23.53 Å of side. Periodic boundary conditions were assumed in order to model the bulk system. Standard minimization procedures ͑Polak-Ribiere conjugate gradients method͒ were followed in order to minimize the so-obtained energy structure, and a subsequent dynamics was run for 1 ns in order to equilibrate the sample. The chosen temperature is high enough to allow local structural equilibration of the sample in this time ͓26͔. The system obtained in this way was used as a starting point for collecting data every 0.01 ps during a MD run of 1 ns. As integration method we have used the velocity-Verlet algorithm with a time step of 1 fs. The simulations were carried out in the constant number of atoms, volume, and temperature ͑NVT͒ ensemble. However, instead of a real temperature-bath coupling ͑Nosé-Hoover or Berendsen thermostats, for instance͒ in order to control the temperature we have followed a rather crude velocity scaling procedure but with a wide temperature window of 10 K. In these conditions, greater temperature fluctuations are allowed but the trajectory is disturbed less. In fact, we have checked that by following this simple procedure we obtain results similar to those obtained with a constant number of atoms, volume, and energy ͑NVE͒ ensemble, which has the proper Newtonian dynamics ͑see Sec. III͒. Moreover, it is worth noting that in a previous work ͓27͔ we checked that the temperature control method used in this work also gives a vibrational density of states similar to that obtained by using the Nosé-Hoover thermostat in this kind of polyisoprene models. After the first 1-ns MD run, two more successive runs of 2 and 20 ns were carried out, collecting data every 0.05 and 2 ps, respectively. As it will be shown in Sec. III, nearly indistinguishable results were obtained from the different simulation runs. Thus, no signature of any aging process was observed during the successive runs confirming local equilibration of the sample. In addition, a different cell was also constructed by the same protocol and equilibration procedure but starting from a different conformation of the parent chain. The results obtained for the dynamic magnitudes, as e.g., F s (Q,t), were similar to those obtained with the first cell within the estimated uncertainties ͑see the Results͒. Finally, in order to compare the structure and vibrational properties of our simulation cell to actual neutron scattering data ͑see Sec. II B͒ the system was suddenly quenched to the glassy state at a temperature of 100 K, similar to the temperature at which the reported neutron scattering measurements were carried out. The density of the system was then adjusted to the estimated value at 100 K (ϭ0.976 g/cm 3 ) by changing the cubic cell edge to 22.63 Å and correspondingly scaling all atomic coordinates. After a minimization procedure similar to that used at 363 K, a dynamic run of 10 ns was carried out in order to accommodate the change in density.

B. Validation
For an amorphous system the structural information is contained in the so-called radial distribution function or in its Fourier transform counterpart: the static structure factor S(Q). In the case of amorphous polymers S(Q) can be measured by neutron diffraction using a fully deuterated sample, i.e., a sample where all hydrogens are replaced by deuterons. The coherent differential cross section d coh /d⍀ measured by neutron scattering can be defined for an isotropic sample and taking the orientational averaging as ͓28͔ where the ͗b i ͘ stand for the nuclear scattering lengths for neutrons and N is the number of atoms ͓the cross sections like coh are given in units of barns/atom, and I coh (Q) is in units of barns/͑sr atom͒, 1 barn ϭ10 Ϫ28 m 2 ]. The scattering lengths of carbon and deuteron atoms are very similar ͑͗b c ͘ϭ0.6648ϫ10 Ϫ14 and ͗b D ͘ϭ0.6674ϫ10 Ϫ14 m͒ thereby these two atoms become almost indistinguishable for neutrons. In this context, the coherent intensity measured results to be just proportional to the static structure factor S(Q). However, what it is usually measured in a neutron diffraction experiment is the total differential scattering cross section given by where d inc /d⍀ is a Q-independent incoherent term defined as For a fully deuterated sample of PI, the value of I inc is 0.10 ͑barns/sr atom͒.
Starting from the atomic coordinates of our simulation runs at 363 K, we have calculated I(Q) by means of expressions ͑1͒-͑3͒ and averaging it for a large number of frames throughout the atomic trajectories. The result obtained is shown in Fig. 1. It is essentially the same result that was previously reported in Ref. ͓29͔ ͑see Fig. 2 of that reference͒. In that paper, the I(Q) obtained from MD simulations at 363 K was compared with the experimental I(Q) measured by neutron diffraction and reported in Ref. ͓30͔, which is also included here in Fig. 1. The conclusion of such comparison was that although the shape of the experimental I(Q) is qualitatively reproduced by the simulation data, there is a clear disagreement concerning the relative intensities of the first and the second peak ͑see Fig. 1͒. Various possibilities for this discrepancy were outlined in Ref. ͓29͔, in particular, the different microstructure of the simulated and actual PI. However, it is worth remarking that neutron scattering measurements of I(Q) in polymers are usually carried out at low temperatures ͑10-100 K͒ well below the glass transition in order to minimize the Debye-Waller factor influence and inelasticity effects. Therefore, in order to properly compare simulation and neutron data of I(Q) and to see to what extent the above-mentioned discrepancy is due to the temperature, we have ''quenched'' our simulation cell to TϷ100 K following the procedure described above. The results obtained for I(Q) at 100 K are also shown in Fig. 1. Now the agreement between the experimental and simulation data is quite good even taking into account the possible nonequilibrium effects of the quenched cell at 100 K and other experimental uncertainties as multiple scattering and inelasticity effects. It is worth remarking that the Q-range included in our comparison corresponds to both intermacromolecular and intramacromolecular correlations in polymer systems. While the first two peaks of I(Q) are mainly dominated by intermacromolecular correlations, the third one at Q ϭ3 Å Ϫ1 is a common feature of all polymers and its Q-position does not change with temperature, thereby suggesting a pure intramacromolecular origin ͓31͔. Therefore, from these results we can conclude that the intermacromolecular and intramacromolecular structure obtained in the simulation cell constitutes a quite reasonable mimic of the actual structure of PI within the limitations of the simulation FIG. 1. Total differential scattering cross section obtained for a fully deuterated sample by neutron diffraction at 100 K ͑᭺͒ ͓30͔ and calculated from our simulations at the same temperature ͑solid line͒ and 363 K ͑dashed-dotted line͒. method used. Although a complete structural characterization of the simulation cell is beyond the aim of this paper and it will be published elsewhere, we can anticipate that there is also a good agreement between the calculated partial static structure factors-corresponding to different partially deuterated PI samples-and those measured in real samples by means of neutron scattering with polarization analysis ͑D7 instrument, ILL, Grenoble, France͒. Some preliminary results can be found in Ref. ͓32͔.
In addition to the structural features, and taking advantage of the ''quenched'' cell, we can also check whether the force field used in the simulations reproduces the main vibrational properties of PI at a low temperature. The vibrational density of states ͓VDOS, Z(E)͔ can be evaluated from inelastic neutron scattering measurements by applying the incoherent Gaussian approximation ͓33͔. Due to the high value of the incoherent scattering cross section of the hydrogen ͑80.27 barns͒ as compared to its coherent cross section ͑1.76 barns͒ and to the cross sections of other typical nuclei in polymers ͑carbon, oxygen, deuterium, etc.͒, the VDOS obtained from inelastic neutron scattering measurements in polymers corresponds to the subset of hydrogens in the sample ͑hydrogen weighted VDOS͒. In Ref. ͓33͔ the VDOS corresponding to two partially deuterated PI samples ͑PId3 and PId5͒ were evaluated from neutron scattering results. They are reproduced in Fig. 2. These two VDOS correspond to the vibrational properties of methyl group hydrogens in the case of PId5 sample ͑where the hydrogens of the main chain have been replaced by deuterons͒ and to the main chain hydrogen vibrations in the case of PId3 sample ͑where the hydrogens of the methyl groups have been replaced by deuterons͒. From the MD simulation results the vibrational density of states can be calculated in general as the spectral density of the velocity autocorrelation function ͓34͔ where the velocity autocorrelation function is calculated in terms of the velocity autocorrelation function of each atom as and N is the number of atoms considered in the calculations. This approach has already been used to characterize the vibrational properties of different molecular crystals ͑see, e.g., Ref. ͓35͔͒. By means of this procedure we have computed Z(E) from the atomic trajectories of the hydrogens corresponding to the ''quenched'' cell at TϷ100 K. In order to compare our results with the experimental ones, we have used in our calculations either the methyl group hydrogens ͑corresponding to the PId5 sample͒ or the main chain hydrogens ͑corresponding to the PId3 sample͒. The results obtained in both cases are shown in Fig. 2. As can be seen, there is a good agreement between the simulation and experimental data at least in the energy range (Eр40 meV) where experimental data are available. This again validates our MD-simulation method and the force field used as well. On the other hand, the calculated Z(E) from MD simulations also shows other maxima at higher energies, which cannot be directly compared with inelastic neutron scattering results. However, although quantum effects certainly affect the highenergy range of Z(E), the energies of the above-mentioned maxima roughly correspond to those of the different infrared bands reported for PI ͑see, for example, Ref. ͓25͔͒.
In addition to the structural and high frequency dynamical aspects discussed above, the simulated slow segmental dynamics, which will be described in the Results section, also reproduces the main experimental features measured by neutron scattering. The reasonable agreement found gives us confidence that those results obtained from simulation, which are not easily accessible by the current experimental techniques, are also reasonably realistic. This is in fact the ultimate goal of any simulation exercise: once a realistic system is simulated, take advantage of it to go beyond the experimental possibilities.

III. RESULTS
From the atomic trajectories obtained in the simulations we have calculated the self-part of the Van Hove correlation function. Assuming isotropic behavior, this function is given by where r is the radial distance from a given particle, N is the number of particles, and r ជ i is the position vector of the ith particle. The angular brackets denote canonical averaging. G s (r,t) is directly related to the incoherent scattering function measured by NS ͓28͔. In the case of polymers, and due to the scattering cross-section values of the different atoms, the NS is dominated by the self-motions of the protons ͓28͔. Therefore, since we want to connect with NS results, we have calculated, first of all, G s (r,t) from the proton trajectories. Moreover, in order to avoid the ''contamination'' from the fast rotational motion of methyl side groups ͓36͔, only main chain protons were considered in the calculations.
In the simplest case, the self-part of the Van Hove correlation function can be approximated by a Gaussian function This form holds rigorously for an ideal gas, for a harmonic crystal and for a system where the motion of the atoms is governed by Langevin's equation. It also holds in any case for t→0, because, under such conditions, the atoms behave as if they were free. However, in more complicated systems at longer times deviations from Eq. ͑7͒ can be found. These deviations can be quantified in a first approximation by the so-called second-order non-Gaussian parameter ␣ 2 (t) ͓37͔ with the even moments of G s (r,t) given by The results obtained for the mean-squared displacement ͗r 2 (t)͘ and for ␣ 2 (t) are shown in Fig. 3͑a͒. ͗r 2 (t)͘ displays three typical dynamic ranges: ͑i͒ a microscopic regime until about 0.8 -1 ps; ͑ii͒ a crossover regime until about 10 ps; and ͑iii͒ a sublinear time dependence that extends until the limit of our simulations ͑20 ns͒. ␣ 2 (t) has a double peak structure where the short-time maximum corresponds to the microscopic regime of ͗r 2 (t)͘ and the other is centered in the crossover regime of ͗r 2 (t)͘. The short-time regime of both ͗r 2 (t)͘ and ␣ 2 (t) strongly depends on the kind of atom considered. This implies that the short-time behavior observed has to be very different in a realistic polymer and in a simple model liquid. For instance, the short-time peak of ␣ 2 (t) almost vanishes when ␣ 2 (t) is calculated from the trajectories of the main-chain carbons instead of main-chain protons ͓see Fig. 3͑a͔͒. This indicates that this peak is related to the fast librational motions of C-H bonds and explains why it is not observed in simple model systems as ''bead-spring'' or ''united-atom'' polymer models. However, the second peak of ␣ 2 (t) shows a similar behavior to that observed in computer simulations of simple Lennard-Jones systems ͓2͔ or in the experiments with colloidal glass-forming systems ͓7͔. It shows a maximum at a time t*Ϸ4 ps centered in the crossover regime of ͗r 2 (t)͘. Once the sublinear behavior of ͗r 2 (t)͘ is well established, ␣ 2 (t) decreases to its long-time limit, zero. Moreover, this second peak is also present in the ␣ 2 (t) corresponding to main-chain carbons, i.e., it is not related to any additional C-H bond motion. The apparent shift of the main peak of ␣ 2 (t) for protons with respect to carbons may be due to the superposition of C-H bond librations and carbonlike motions. When atoms ͑protons in this case͒ are participating simultaneously in different dynamical processes, the resulting non-Gaussian parameter is not the simple addition of those ␣ 2 (t) corresponding to the processes involved, but is given by a more complicated expres- sion ͓14,36͔. Since this ␣ 2 peak is rather broad, we have characterized it not only by t*, but also by the full width at half maximum ͑FWHM͒. This time regime is represented by the shaded area in Fig. 3. The hatched area in the upper part of this figure shows for comparison the ␣ 2 values estimated by means of the procedure described in ͓14͔ from NS data of PI at 340 K. As can be seen, an accurate determination of ␣ 2 (t) from experimental data is rather difficult ͑see ͓14͔ and ͓15͔ for a detailed discussion about this question͒. From the calculated G s (r,t) we have computed the selfcorrelation function F s (Q,t) as This corresponds to the incoherent intermediate scattering function measured by quasielastic NS techniques in protonated polymer samples. Let us remember that if the Gaussian approximation is fulfilled the self-correlation function is expressed just in terms of the mean-squared displacement The results obtained for our polymer are plotted in Fig. 3͑b͒ for different Q values. F s (Q,t) exhibits the two-step decay, which is characteristic of glass-forming supercooled liquids in general. This two-step feature is not present when F s (Q,t) is calculated at a very high temperature ͑Ϸ513 K͒ where the system behaves as a simple ''polymer liquid'' ͓see Fig. 3͑b͔͒. In this case, the second peak of ␣ 2 (t) also vanishes ͓see Fig.  3͑a͔͒ suggesting that this peak is a main signature of supercooled liquid dynamics. Here we will focus on the slower decay of F s (Q,t) usually known as the ␣ relaxation. We can immediately see that in the low Q range, the time range where ␣ 2 has significant values ͑FWHM͒ only covers the initial part of the slow decay of F s (Q,t). However, as Q increases, this time range starts to cover almost completely the slow decay. Thereby, even without any analysis, we should expect a different influence of the non-Gaussian events on the decay of the correlations in the time region close to its characteristic time scale depending on the Q value. Before analyzing the slow decay of F s (Q,t), we will comment on the main uncertainties affecting the calculated F s (Q,t). The accuracy of the calculated values of F s (Q,t) at longer times, close to the limit of the MD-simulation run t run , depends on the number of time origins that are available for calculating the ensemble average of the Van Hove correlation function. Thus, the uncertainties increase with the value of the time t at which the function is calculated ͑see, e.g., ͓38͔ as a general reference͒. Figure 4͑a͒ shows for two representative Q values the intermediate scattering function corresponding to three successive dynamic runs of different duration: t run ϭ1, 2, and 20 ns. We can realize that, for example, the points calculated from the 2 ns run in the time range tՇt run slightly deviate from the corresponding points calculated from the 20 ns run. However, apart from this tendency, which can be understood in terms of the inherent uncertainties of the calculation, the points corresponding to the three runs nicely superimpose, indicating that there is not any signature of aging processes. This confirms the local equilibrium state of our sample in the meaning that the density-density correlation function at the intermacromolecular level ͑structural ␣ relaxation͒ decays to zero within our simulation time window ͑see Ref. ͓26͔͒. However, it is likely that this time does not allow the long-length scale conformational properties to equilibrate. This is evident in the behav- dϳ5 Å is the average intermacromolecular chain distance defined as dϭ2/Q max ͓Q max is the Q value at which S(Q) shows the first maximum͔. According to the results described in Ref. ͓39͔, this implies that the sublinear behavior of ͗r 2 (t)͘ observed from about 10 ps to the limit of our simulations, actually corresponds to the small-scale motions involved in the ␣ relaxation. A full equilibration of the large scales would give rise to a long-time linear behavior of ͗r 2 (t)͘ corresponding to the diffusion of the chain center of mass. This could only be observed in a fully atomistic simulation at a very high-unrealistic-temperature or by simulating very short chains. Thereby, we can consider that in our sample the large-scale conformational fluctuations are frozen out during the simulation runs. Although this can, in principle, mean a potentially serious problem of ergodicity, it is likely that the ␣ relaxation is mainly controlled by the packing density and the intermolecular structure rather than by large-scale conformational properties. For instance, in Refs. ͓29,40͔ it is reported that chain dimension seems to be not a critical parameter for local segmental dynamics ͑␣ relax-ation͒ in a polyisoprene model similar to that used in this work. With respect to this and as an example, we show in Fig. 4͑b͒ the F s (Q,t) obtained at three different Q values from similar dynamic runs (t run ϭ2 ns) but in two different cells, which were constructed starting from a different conformation of the parent chain. As can be seen, there is a very good agreement between the two sets of results. Finally, as mentioned in the preceding section, we have also checked whether or not our results depend on the statistical ensemble ͑NVT͒ and the temperature control used. To do this we have carried out an additional 2 ns run starting from the final system configuration and temperature obtained at the NVT conditions, but in the microcanonical NVE ensemble ͑which has the proper Newtonian dynamics͒. We have also taken care that the drift in temperature in this NVE run was kept less than 1%. The results obtained for F s (Q,t) at different Q values are shown in Fig. 4͑c͒ in comparison to the results obtained in the NVT ensemble. As can be seen, there is a very good agreement between the two sets of data. Moreover, as it was previously shown in Fig. 3͑a͒, the ␣ 2 (t) obtained in the NVE condition also agrees with that obtained in the NVT ensemble. Now, following the procedure used for analyzing NS data, we have fitted the slow decay of F s (Q,t) to a Kohlrausch-Williams-Watts ͑KWW͒ or stretched exponential function where ␤ is a parameter measuring the deviation from a single exponential form (0Ͻ␤Ͻ1) and a relaxation time, which, in principle, depends on Q. A is a generalized Lamb-Mössbauer factor ͑LMF͒ giving account for the first fast decay of F s (Q,t). The fitting was carried out in the typical time range covered by neutron backscattering and spin echo techniques ͑Ϸ5 to Ϸ20 ns͒. In a first fitting, A, , and ␤ were taken as free parameters. As can be seen in Fig. 5, ␤ resulted to be hardly dependent on Q around the value of 0.4, which is the value of ␤ usually obtained by NS and spectroscopic methods in PI ͓13,30͔. This gives additional support to our MD simulations. Therefore, in a second fitting step ␤ was fixed to 0.4. The resulting fitting curves are plotted in Fig.  3͑b͒. They describe quite well the F s (Q,t) values in the fitting time range. The Q dependence found for A follows a LMF-like behavior which gives a value of ͗u o 2 ͘Ϸ0.41 Å 2 , in the range of the values usually found by NS in polymers in such a temperature range. The obtained Q dependence of is shown in Fig.  6. In the Q-range Qр1 Å Ϫ1 , nicely follows the power law where a is a temperature-dependent prefactor. This kind of Q dependence has been experimentally observed in several polymers including PI in a similar Q range ͓11-13͔. This once again gives support to our simulation and polymer model. Introducing the results obtained for the Q dependencies of and A ͓Eqs. ͑13͒ and ͑14͔͒ in Eq. ͑12͒ one finds

͑15͒
This coincides just with the expression corresponding to the Gaussian approximation for F s (Q,t) ͓Eq. ͑11͔͒, where we can identify ͗r 2 (t)͘ϭ2͗u o 2 ͘ϩ6(t/a) ␤ . This implies that in the Q range QՇ1 Å Ϫ1 Gaussian behavior is found for the dynamics of the ␣ relaxation, without any previous assumption on the validity of such an approximation. However, as Q increases, (Q) deviates from the power law given by Eq. ͑14͒ and a clear crossover towards an approximate Q Ϫ2 law is evident in Fig. 6. This crossover takes place in a Q range that we will characterize by a crossover Q, Q*Ϸ1.3 Å Ϫ1 . We note that Q* lies very much in the range where S(Q) shows the first maximum, Q max . A similar breakdown of the Gaussian approximation was found in the neighborhood of Q max in Refs. ͓19,20͔.
The approximate law ϰQ Ϫ2 found for QϾQ*, together with the above-mentioned Q dependence of A ͓Eq. ͑13͔͒ means a strong deviation of F s (Q,t) from the Gaussian form. The crossover from Gaussian to non-Gaussian behavior at Q* can be understood as a consequence of the influence of the non-Gaussian events at the time scale tϷ(Q). This is evident in Fig. 6 where we have assigned to each Q value the value of ␣ 2 at tϭ(Q), ␣ 2 ͓(Q)͔. By means of this representation, we immediately realize that when Q approaches Q*, ␣ 2 ͓(Q)͔ takes significant values. In this representation, ␣ 2 ͓(Q)͔ reaches its maximum value at the Q value ͑Ϸ3 Å Ϫ1 ͒ at which (Q)Ϸt*. However, due to the large width of ␣ 2 (t), the crossover from Gaussian to non-Gaussian behavior already takes place when (Q)Ϸ15t*.

IV. DISCUSSION
First of all, it is worth remarking that the crossover described above takes place in a Q range that is close to the limit usually explored by incoherent neutron scattering measurements (QϷ1.8 Å Ϫ1 ). There, the uncertainties are usually higher, since the intensity of the signal is diminished by the Lamb-Mössbauer factor. This could be the reason why this crossover was never reported, although deviations from the Q Ϫ2/␤ law in the high Q range have already been ob-served in different glass-forming polymers ͓16,41͔. Now there is a chance to extend the Q range until QϷ5 Å Ϫ1 by means of, e.g., the IN13 spectrometer, ILL, Grenoble, France. This should allow, in principle, an experimental observation of this crossover. However, we have to point out that this is not an easy task. From an experimental point of view there are many difficulties. For instance, the full Q-range from about 0.3 to 5 Å Ϫ1 cannot be covered by one single spectrometer. Moreover, it is also very difficult to observe the second step of F s (Q,t) centered in the experimental dynamical range at the different Q values but at the same temperature. Finally, the amplitude of this second step is rather low in the high Q regime ͓see, e.g., Fig. 3͑b͔͒. In spite of these difficulties, experiments in this direction are now being planned in PI.
How do the results reported here compare with those previously published about MD simulations of polymer systems? As it has already been mentioned in the Introduction, the Q dependence of simulated F s (Q,t) corresponding to polymer models has been investigated in recent papers ͓19-21͔, though such papers were mainly focussed on checking the mode-coupling theory ͑MCT͒ ͓42͔. In Ref. ͓21͔, Van Zon and De Leeuw, using a ''united atom'' model for a polymer melt, reported a Gaussian behavior of F s (Q,t) in the low Q regime (QϽ1 Å Ϫ1 ). However, although they also reported deviations from this behavior for higher Q values, a direct comparison with our results is not possible because they found a ␤ value that strongly depends on Q. On the other hand, the Q dependence of the incoherent ␣-relaxation time in a simulated simple polymer model is discussed by Bennemann, Baschnagel, and Paul ͓19͔. They find a law compatible with ϳQ Ϫ2/␤ ͑Gaussian approximation͒ for the two lower Q values investigated as well as a crossover towards a law close to ϳQ Ϫ4/3 for larger Q ͑see also Ref. ͓20͔͒. The reported crossover is rather smooth and takes place at about Qϳ0.7Q max when S(Q) starts to have significant values ͑see Fig. 1 of Ref. ͓19͔͒. It is worth remarking that this crossover seems to have a behavior qualitatively similar to that reported here even though the system used by Bennemann, Baschnagel, and Paul is rather simple ͑''bead-spring'' model͒ and does not consider valence or rotational potentials. The authors of Ref. ͓19͔ try to understand this trend in terms of the MCT, which was the main subject of that paper. In the high Q limit, where the MCT equations can be analytically solved, the incoherent scattering function F s (Q,t) shows a KWW behavior with the stretching exponent ␤ equal to the von Schweidler exponent b ͑see Ref. ͓43͔͒. In this regime, the expected Q dependence of the KWW relaxation time is given by the power law ϳQ Ϫ1/b . The b value (bϳ0.75) reported in Refs. ͓19,20͔ seems to be in good agreement with the law found ϳQ Ϫ4/3 . However, it is worth noting that in the framework of the MCT the Q dependence ϳQ Ϫ1/b should be observed in the high Q limit. The authors of Ref.
͓19͔ themselves realized this problem. If we try to consider this MCT interpretation for the crossover reported in this paper, the corresponding b value should be of the order of 0.5, which can be reasonable for realistic polymer systems. It is worth remarking that a KWW functional form for F s (Q,t) ͓Eq. ͑12͔͒ with ␤ϭb and ϳQ Ϫ1/b reads as F s (Q,t) ϭA exp͓ϪQt b a(T)͔, where a(T) accounts for the temperature dependence of . The Q dependence of this equation strongly deviates from that expected in the Gaussian case ͓see Eq. ͑11͔͒. In fact, the time dependence of the non-Gaussian parameter ␣ 2 (t) calculated in the framework of the MCT for a hard-sphere system ͑HSS͒ ͓44͔ shows a qualitatively similar behavior to that reported here. In the framework of the MCT the non-Gaussian behavior is related to the decaging processes and ␣ 2 (t)Ͼ0 would indicate a more fuzzy cage boundary. In this framework, we can try to obtain a first estimation of the so-called mean characteristic localization length r sc from the value of the mean squared displacement at t*, ͗r 2 (t*)͘ while assuming that ͗r 2 (t*)͘ϳ6r sc 2 . This gives r sc ϳ0.45 Å. Taking into account that the average intermacromolecular chain distance d is given by d Ϸ2/Q max (Q max Ϸ1.3 Å Ϫ1 ), we obtain that r sc ϳ0.094d. It is worth remarking that this is a similar value as that obtained for the HSS ͓44͔ and for the above-mentioned ''beadspring'' polymer model as well ͓20͔. As in these cases, the value found here for r sc also approximately fits to the Lindemann criterion for melting.
On the other hand, as it has already been mentioned in the Introduction, dynamical heterogeneity of supercooled systems is often discussed in terms of the self-part of the Van Hove correlation function or its Fourier counterpart F s (Q,t). Deviations from the Gaussian behavior, similar to those reported here, have usually been considered as a signature of the dynamically heterogeneous behavior ͑see, e.g., ͓2,7͔͒. Positive ␣ 2 (t) suggests that the probability for the particle to move very far is enhanced relative to the one expected for a random-walk process. In the case of Lennard-Jones systems, it was shown by Kob et al. ͓2͔ that the number of these ''faster particles'' giving rise to ␣ 2 (t)Ͼ0 was only in the range of 5% of all particles in the system. In our case, taking into account the obtained values of ␣ 2 (t), we can estimate a similar percentage of faster protons in the system. In this framework, the crossover from Gaussian to non-Gaussian behavior at Q* could be interpreted as a crossover from homogeneous to heterogeneous dynamical behavior. For instance, a possible heterogeneous scenario giving account for the law ϳQ Ϫ2 found at QϾQ*, was already outlined in Ref. ͓13͔ in the simplified case of Aϭ1. It is also worth remarking that the behavior found is also compatible, at least qualitatively, with a simple sublinear jump diffusion model resulting from jumps with the jumping length distributed ͓16͔. In this framework, the asymptotic low Q dependence of is given by ϰQ Ϫ2/␤ implying a Gaussian behavior of F s (Q,t). However, in the high Q regime, deviations from this Q dependence take place ͑non-Gaussian behavior͒ due to the intrinsic jump-length distribution for the diffusive jumps. Finite jump lengths tend to cause a bending of (Q), which is qualitatively similar to that found in simulations at QϾQ* ͑see, e.g., Fig. 1 Since the characteristic time scale of ␣ 2 (t), t*, corresponds to the decaging region of ͗r 2 (t)͘, it is quite reasonable to think that the non-Gaussian or heterogeneous behavior is in fact due to slightly different atomic environments for each particle in this region. This could be of particular relevance for a complicated system, such as a polymer. There-fore one can expect that the heterogeneous or non-Gaussian effects would start to vanish once the particles move to larger distances and thereby leave the decaging region. In fact, Fig.  3͑a͒ shows that ␣ 2 (t)→0 at about tϷ2 ns, a time for which ͗r 2 (t)͘ approaches the square value of the average intermacromolecular distance dϭ2/Q max . This might explain in a natural way why the crossover takes place in the range of Q max . It is worth noting that this heterogeneous interpretation is, in principle, compatible with the MCT approach discussed above ͑see also ͓44͔͒. Although a microscopic characterization of the effect of the different possible atomic environments in our sample will certainly require much more simulation effort ͑it will be the subject of future work͒, we have calculated here first, as an example, the ␣ 2 (t) behavior for three different kind of protons in the monomer of PI. The results obtained are shown in Fig. 7 where it is clear that different protons exhibit different ␣ 2 (t) behavior. Now, in the framework of the homogeneous/ heterogeneous interpretation, we can come back to the question: why does the dynamics of the ␣ relaxation seem to be heterogeneous when it is observed by spectroscopic methods in general and in particular by NMR techniques? To answer this question we have to know what is the time scale of the ␣ process as observed by these methods. From the atomic trajectories of our MD-simulation runs we have also computed, as an example, the second-order orientational autocorrelation function M 2 (t)ϭ͗P 2 ͓cos (t)͔͘, where P 2 is the second Legendre polynomial and (t) describes the orientation of a given vector at time t relative to its orientation at tϭ0. This is the correlation function determining the NMR experimental observables as, for example, the spin-lattice-relaxation time T 1 measured by 13 C NMR. M 2 (t) was calculated for the different main chain C-H bonds of PI monomer. The resulting curves are depicted in Fig. 8. The slow decay of M 2 (t) is produced by the ␣ relaxation, and the corresponding characteristic time scale NMR is the characteristic time deduced from NMR results for the ␣ relaxation. In Fig. 8 we have also included F s (Q,t) at different Q values. The direct comparison of the functions shows that NMR is in the range of the time scale of F s (Q,t) at QϷ1 -1.2 Å Ϫ1 . These Q values lay in the region where the crossover of (Q) takes place. This can also be deduced from Fig. 6, where the range of time scales observed by NMR, NMR has been represented by a shadowed area giving account for the different possible choices of the C-H bond orientation followed. Thus, we can expect that the non-Gaussian effects influence to a certain extent the slow molecular dynamics in such time scales. This suggests why the ␣ relaxation observed by NMR techniques reveals heterogeneous contributions. It is worth emphasizing that the time scale of the ␣ process, which is probed by other relaxation techniques as, for instance, dielectric relaxation, is also in this range. Data reported for different polymers show that the dielectric relaxation time at high frequency coincides with the relaxation time corresponding to incoherent NS, (Q), at QϷ1 Å Ϫ1 ͓45͔. On the other hand, since usually universal behavior for the ␣ process is observed-the time scales deduced from different correlators show similar temperature dependencies-we can expect that the results obtained for the temperature investigated are extensible to other temperature ranges.
How do the different time scales discussed ͓t*,(Q), NMR ,...͔ compare with the actual relaxation time ␣ ? This characteristic time should be defined as the relaxation time of the dynamic structure factor F(Q,t) at the Q value Q max of the first intermolecular maximum of S(Q) ͓S(Q)ϵF(Q,tϭ0)͔. The so-defined time scale has the meaning of the relaxation time of the slow decay of the density-density correlation function at Q max . This is the way the ␣-relaxation time is defined, for instance, in the framework of the mode-coupling theories of the glass transition ͓42͔. From the atomic trajectories of the simulation runs, F(Q,t) can be calculated as the Fourier transform of the total Van Hove correlation function but now considering all the atoms in the simulation cell. By means of this procedure we have computed the densitydensity correlation function F(Q,t)/F(Q,0) for different Q values in the range 0.5рQр5 Å Ϫ1 . This is the correlation function measured by neutron-spin-echo ͑NSE͒ when a fully deuterated polymer sample is used. Some of the results obtained are shown in Fig. 9͑a͒ as an example. As F s (Q,t), F(Q,t)/F(Q,0) also shows the well-known two-step decay. The second step was also fitted by a KWW function. We followed a procedure similar to the one described in the Results section for F s (Q,t). In a first fitting, the three parameters A, ␤, and , were taken as free. ␤ resulted to be modulated mirroring S(Q) around similar values to those obtained in the case of F s (Q,t) ͑see Fig. 5͒. Therefore, as in the case of F s (Q,t) and following the usual way of analyzing NSE experimental data, the ␤ value was fixed to 0.4 in a second fitting series. As can be seen in Fig. 9͑a͒, the KWW fitting curves can describe the slow decay of F(Q,t)/F(Q,0), though the deviations from the imposed functional form become more evident than for the self-correlation function. In spite of this, it is clear that the time scales can be unambiguously determined from these fittings. The obtained Q dependence of the KWW time ͑we will call it coh -coherent-in the following͒ is shown in Fig. 9͑b͒ together with the (Q) corresponding to F s (Q,t). As can be seen, coh (Q) is modulated by S(Q) mainly in the Q range of the first intermolecular maximum of S(Q). A similar behavior has been found in other systems where the Q dependence of the collective time has been studied, e.g., diatomic molecules ͓46͔, water ͓47͔, and simple polymer models ͓20͔. This modulation can be understood as a consequence of some kind of ''de Gennes narrowing'' ͓48͔. In the high Q regime, coh (Q) tends to follow the Q dependence of (Q) corresponding to F s (Q,t). Due to the modulation of coh (Q), the value of ␣ ϭ coh (Q max ) is higher than the corresponding value of (Q) ͓ ␣ ϳ6(Q max )͔ and in fact is close to the time scale of F s (Q,t) at QϷ0.85 Å Ϫ1 . In order to compare ␣ with the other time scales discussed above, we have plotted it in Fig.  6. As can be seen, ␣ Ͼ NMR Ͼt*. Moreover, ␣ lies in the time regime at which the slow dynamics behaves as Gaussian. Therefore, in the framework of this interpretation of the crossover found, we can understand why, although being intrinsically heterogeneous, the slow dynamics in glassforming polymers looks homogeneouslike at the time scale of the structural ␣ relaxation.

V. CONCLUSIONS
We have shown that the Q dependence of the KWWincoherent relaxation time ϳQ Ϫ2/␤ , which has been ob-served in different glass-forming polymers in the range Q р1 Å Ϫ1 , is nicely reproduced by MD simulations in a fully atomistic model system of polyisoprene. In the Q range of the first maximum of S(Q) we have obtained a crossover to a Q dependence close to ϳQ Ϫ2 . We have shown that this crossover is a consequence of the effect of the non-Gaussian events at the time scale (Q). As Q increases, (Q) approaches the characteristic time t* of the non-Gaussian processes and (Q) deviates from the Gaussian behavior given by ϳQ Ϫ2/␤ . Various possible interpretations of the crossover from Gaussian to non-Gaussian behavior have been discussed. The interpretation in terms of an underlying crossover from homogeneous to heterogeneous dynamics opens a possible way for rationalizing the apparent contradiction between neutron scattering and relaxation techniques concerning dynamical heterogeneity of the ␣ relaxation. The closer the time scale of the ␣ process probed by a particular technique is to the t* range, the more sensible this technique is to the heterogeneous dynamics.