Methyl group dynamics above the glass transition temperature: a molecular dynamics simulation in polyisoprene

Chemical Physics 261 (2000) 47±59

www.elsevier.nl/locate/chemphys

Methyl group dynamics above the glass transition temperature: a molecular dynamics simulation in polyisoprene F. Alvarez, A. Arbe, J. Colmenero * Departamento de Fõsica de Materiales y Centro Mixto CSIC-UPV/EHU, Universidad del Paõs Vasco, Apartado 1072, 20080 San Sebasti an, Spain Received 11 February 2000; in ®nal form 11 May 2000

Abstract We have carried out molecular dynamics (MD) simulations of polyisoprene at 363 K, about 150 K above the experimental glass transition temperature, using the I N S I G H T and D I S C O V E R programs from MSI Inc. with the Polymer Consortium Force Field. The model system was built using the MSI amorphous cell construction protocol with periodic boundary conditions. Starting from the self part of the van Hove correlation function, the incoherent intermediate scattering function was calculated for the protons in the main chain and in the methyl groups (MGs). The dynamics of the latter ones can be well described assuming decoupled segmental dynamics and rotations in a threefold potential. The limits of such an approximation are also discussed. We ®nd the existence of a distribution of potential barriers for MG rotation which is very similar to that deduced from low temperature (150 K) MD-simulation results and inelastic neutron scattering measurements. The glass transition would thus hardly modify this distribution. Ó 2000 Elsevier Science B.V. All rights reserved. Keywords: Methyl group dynamics; Glass forming polymers

1. Introduction Recent neutron scattering experiments of quantum and classical methyl group (MG) dynamics in glassy polymers have shown the signi®cant role played by the distribution of potential barriers for MG rotation [1±3]. Due to the structural disorder inherent to the glassy state, the MGs in a glassy polymer have di€erent local environments arising, in principle, from both the lack of regularity of the main chain (MC) conformation and

*

Corresponding author. Tel.: +34-943-015962; fax: +34-943212236. E-mail address: [email protected] (J. Colmenero).

the di€erent local packing conditions. This manifests itself in a di€erent mean ®eld or single particle potential barrier for each rotating MG. The results for both quantum and classical temperature regimes of di€erent polymers can be described consistently in the framework of the so-called rotation rate distribution model (RRDM) which was ®rst proposed in 1994 [4]. This model assumes a pure single particle threefold potential fV …u† ˆ V3 ‰1 ÿ cos …3u†Š=2g for MG rotation, together with a Gaussian distribution of potential barriers G(V3 ). The distribution of V3 yields distributions of the parameters driving the MG dynamics in both the classical and quantum regimes, namely, the activation energy for hopping through the barriers, the torsional libration levels and the

0301-0104/00/$ - see front matter Ó 2000 Elsevier Science B.V. All rights reserved. PII: S 0 3 0 1 - 0 1 0 4 ( 0 0 ) 0 0 2 2 5 - 1

48

F. Alvarez et al. / Chemical Physics 261 (2000) 47±59

tunnelling frequencies. These are the parameters that can experimentally be obtained by means of neutron scattering spectroscopy [1±4] in particular. Due to the fast time scale of the MG dynamics, molecular dynamics (MD) simulations are also a suitable tool to gain insight into this problem. For example, the calculation of the density of librational states for the MGs is quite straightforward and does not require lengthy MD runs as it was shown in a previous paper [5]. From the atomic trajectories which are obtained in a MD simulation, di€erent correlators can be constructed, in particular, the incoherent and coherent scattering functions which are measured by neutron scattering. However, the advantage of MD simulations is that we can calculate directly the correlations in the real space (the van Hove correlation functions), while in a neutron scattering experiment only their Fourier transformations are obtained. Moreover, MD simulations allow to follow selectively the dynamics of particular molecular groups like the MGs. In principle, this can also be carried out in real neutron scattering experiments by chemical manipulation (partial deuteration) of the samples and taking advantage of the di€erent scattering cross-sections of protons and deuterons. However, in practice, this procedure is not always possible. In addition, MD simulations have the advantage that parameters de®ning the force ®eld used in the simulations can easily be manipulated in order to investigate their in¯uence on the MG dynamics. Therefore, MD simulations can be considered as a complementary tool to neutron scattering techniques. In previous works [5,6], we carried out a MD simulation of MG dynamics in polyisoprene (PI) in the glassy state at a temperature of 150 K. PI was chosen because it is a rather simple polymer which can be simulated in a realistic way by commercial force ®elds and well-established methods of MD simulations (see for example Refs. [7,8]). Moreover, there are neutron scattering data for MG dynamics already reported in the literature [9,10] which allow to validate MD simulation results. In particular, the calculated static structure factor S(Q) (see Fig. 1 of Ref. [6]) was found to agree quite well with the experimental data, validating the structure obtained in the simulation

Fig. 1. Comparison between the static structure factor S(Q) of protonated PI as measured by neutron scattering at 10 K (dots) and calculated at 363 K from the cell used in this simulation (line). Inset: S(Q) measured for 1,4-PB at 4 K (


) and 270 K (Ð) (from Ref. [11]).

cells, even if the simulation cells were in that work  dimension side). rather small (cubic cell of 15.5 A The main results obtained in such a ®rst MD investigation of MG dynamics in PI can be summarized as follows. First of all our results supported the generally assumed threefold approximation for the single particle MG potential. Moreover, the density of states for MG torsional librations agrees quite well with inelastic neutron scattering results on PI [10] and shows a broad feature re¯ecting a distribution of potential barriers. By changing the parameters of the force ®eld used, we concluded that the width of the distribution of potential barriers, which can be quanti®ed in the framework of the threefold approximation, is mainly controlled by the nonbonded interactions. As mentioned above, the MD simulations previously reported were performed at 150 K in the glassy state of PI where the segmental motions associated to the a-relaxation (glass transition) are frozen. Now, the question arising is: what is the behaviour of the MG dynamics in the high temperature regime well above the glass transition?

F. Alvarez et al. / Chemical Physics 261 (2000) 47±59

There are two main points which should be addressed. First of all, the relationship between segmental dynamics ± which is active at high temperature ± and the MG dynamics. Are these statistical independent processes in a good approximation? Second, does the distribution function of potential barriers (or activation energies) for MG rotation change with temperature? Trying to address these questions, we have carried out a MD simulation of both, MG and MC dynamics of PI, at a high temperature of 363 K. In this work, a  bigger cubic cell was constructed of about 23 A side by following similar methods to those used in our previous work. The force ®eld used was also the same. The results obtained are compared with the neutron scattering data which are available in the high T > Tg range (Tg : glass transition temperature). 2. Simulation method and validation The simulations were carried out by using the (Insight II 4.0.0 P version) and the D I S C O V E R -3 module from MSI (Molecular Simulations Inc. San Diego) with the Polymer Consortium Force Field (PCFF). The model system was built by means of the Amorphous Cell protocol. A cubic cell containing one polymer chain of 100 monomer units [±CH2 ±CH@C(CH3 )±CH2 ±]100 was constructed simulating a temperature of 363 K and subjected to the corresponding experimental density at that temperature. Such a density leads to  at the side. Periodic a cell dimension of 23.53 A boundary conditions were assumed in order to model the bulk system. The microstructure of the sample was 82% trans and 18% cis units. Standard (conjugate gradients) minimization procedures were followed in order to minimize the thus obtained structure, and subsequent dynamics were run for 1 ns in order to equilibrate the sample. The chosen temperature is high enough to allow structural equilibration of the sample in this time. The system obtained in this way was used as a starting point for collecting data. Initially, the atoms of the system are assigned velocities from a Boltzmann distribution at the corresponding temperature. Then, Newton's laws INSIGHT

49

of motion are solved numerically in 1 fs time steps, using the velocity Verlet algorithm to follow the evolution of the conformations. An NVT ensemble was used with the velocity scaling algorithm to keep the temperature constant. A ®rst MD at 363 K was run for 1 ns using the D I S C O V E R -3 program collecting data every 0.01 ps, and a subsequent one (taking the previous output sample as an input for the following dynamics) was run for 2 ns collecting data every 0.05 ps. The results of the second run agreed to those of the ®rst run, indicating that the sample was well equilibrated at this high temperature. Therefore, only one cell was used in this work. In order to validate the structure obtained in the simulation cell, the static structure factor S…Q† was calculated as in our previous work [6]. In Fig. 1, the S…Q† obtained now is compared with the experimental data measured with spin polarization analysis by means of the D7 spectrometer (ILL, Grenoble, France) [9]. These new simulation results show an even better agreement with the experimental data than that obtained from the smaller cell used in Ref. [6]. However, some differences between the D7 and the simulation results can still be observed. It is worth noting that the error bars shown for the D7 points, though quite large, only re¯ect the statistical errors involved in the measurement. There are other e€ects as multiple scattering that have not been corrected and could lead to slight changes in the shape of the experimental peak (the corrected data would show a sharper peak) which would improve the agreement. However, even taking into account these considerations, the ®rst peak of S…Q† in the case of the simulation shows a shift towards lower Q values and a broadening in the low Q-¯ank with respect to the D7 data. These di€erences can nevertheless be rationalized taking into account that the simulation was carried out at 363 K and the experiments at 10 K (this kind of experiments are usually performed at very low temperatures in order to minimize the Debye±Waller factor in¯uence on the data). The evolution with temperature of the ®rst S…Q† peak in polymers generally shows the features mentioned above, as can be seen in the inset of Fig. 1, where, a typical example of experimental data is presented. The curves, extracted

50

F. Alvarez et al. / Chemical Physics 261 (2000) 47±59

from Ref. [11], correspond to the archetypal polymer 1,4-PB at 4 and 270 K (around 90 K above Tg ). Thus, the apparent discrepancies found between our simulation results and the D7 points can be qualitatively understood and it can be concluded that the agreement is quite reasonable.

3. Results From the atomic trajectories obtained in the simulation the van Hove correlation function [12] 1 can be calculated. This consists of two parts, a self part Gs …~ r; t†, that relates to the single atom r; t†, that involves motion, and a distinct part Gd …~ correlations between atom pairs. These two parts are directly related to the incoherent and coherent scattering functions, respectively, which are accessible by neutron scattering. The incoherent scattering function dominates the scattering measured when the sample mainly consists of protons, as it is the case of polymers in general. This is due to the high value of the incoherent cross-section of the proton rH inc (80.27 barns) as compared to its coherent cross-section rH coh (1.76 barns) and to the cross-sections of other nuclei (carbon, oxygen, deuterium, ...) [16] that polymers usually consist of. For this reason, we will focus on the self part of the van Hove correlation function. This is given by * + N h i 1 X r; t† ˆ d~ r ‡~ ri …0† ÿ~ ri …t† ; …1† Gs …~ N iˆ1 where N is the number of atoms, and ~ ri …t† denotes the position of the ith atom at time t. The bracket r; t† is indicates standard ensemble averaging. Gs …~ thus the ensemble average probability of ®nding a given particle for a time t at the position ~ r if the same particle was at the origin for time t ˆ 0. In the simplest case, the self part of the van Hove correlation function can be approximated by a Gaussian function

1

The same basic ideas were published by Glauber [13±15].



Ggauss …~ r; t† s

a…t† ˆ p

3=2

  exp ÿ a…t†r2 :

…2†

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. (2) can be found. r; t† are given by the The even moments of Gs …~ following expression: Z

2n  r; t†d~ r r …t† ˆ r2n Gs …~ Z 1 r2n 4pr2 Gs …r; t† dr; …3† ˆ 0

where isotropic behaviour has been assumed. For n ˆ 1, the well-known mean squared displacement hr2 …t†i is obtained. It is straightforward to show that in the Gaussian case (Eq. (2)), hr2 …t†i is given by

2  3 : …4† r …t† ˆ 2a…t† Neutron scattering experiments provide information on the Fourier transform of the van Hove correlation function through the momentum transfer (Q) and the energy transfer dependence of the number of neutrons scattered by the sample. The incoherent intermediate scattering function is thus given by Z   ~ t† ˆ Gs …~ ~r d~ r; t† exp iQ~ r Fs …Q; Z 1 sin …Qr† dr; …5† 4pr2 Gs …r; t† ˆ Qr 0 where again an isotropic system has been considered. In the framework of the Gaussian approximation, this function can easily be calculated:   Q2 gauss Fs …Q; t† ˆ exp ÿ 4a…t†   hr2 …t†i 2 Q : …6† ˆ exp ÿ 6

F. Alvarez et al. / Chemical Physics 261 (2000) 47±59

51

In the general non-Gaussian case, a further approximation for Fs …Q; t† can be obtained by using the following expression: " # 2 hr2 …t†i 2 a2 …t†hr2 …t†i 4 Q Q ‡ Fs …Q; t† ˆ exp ÿ 72 6 …7† that results from the cumulant expansion [17] up to the order of Q4 . In this equation, we have introduced the so-called non-Gaussian parameter of order 2, a2 …t†, which is de®ned as a2 …t† ˆ

3 hr4 …t†i ÿ 1: 5 hr2 …t†i2

…8†

The deviations from the Gaussian approximation usually found in complex systems can be quanti®ed in a ®rst approximation by this parameter. In this work, we have calculated the self part of the van Hove correlation function given by Eq. (1) for two kinds of atoms in the simulated cell: the protons in the MC, and the protons in the MGs. In the following, we will distinguish the results corresponding to these subsystems by superscripts. In order to achieve adequate averaging, GMC s …r; t† and …r; t† were calculated by considering a maxiGMG s mum of 400 di€erent time origins. The other functions and parameters introduced above were thereafter obtained for both kinds of atoms. As commented before, the incoherent scattering from protons dominates the experimental spectra due to the high value of rH inc . For this reason, the incoherent scattering function FsMC …Q; t† calculated from the simulations would be the right quantity to be compared to results from neutron scattering experiments on a sample with deuterated MGs and the rest protonated. In analogy, FsMG …Q; t† would correspond to a real sample with protonated MGs and deuterated MC. Fig. 2 shows the results obtained for the func2 MG …r; t† tions 4pr2 GMC s …r; t† (Fig. 2(a)) and 4pr Gs (Fig. 2(b)) for di€erent times. In both cases, the maximum of the probability distribution function continuously shifts towards larger values of the distance from the origin when time increases. This is a signature of a di€usive process. As commented

Fig. 2. Probability distribution functions obtained from the self part of the van Hove correlation function calculated for the protons located (a) in the MC of PI and (b) in the MG at the di€erent times indicated. Solid lines show the functions obtained assuming the Gaussian approximation for the correlation function at 4 ps with the same mean squared displacement.

in Section 1, at the temperature here investigated the di€usive a-relaxation is expected to contribute to the dynamics of our system. Therefore, the behaviour found for the self part of the van Hove correlation function at longer times can be attributed to the main relaxation. On the other hand, Fig. 2 shows that, for a given time, the maximum of the probability distribution for the MG atoms is located at higher r values than for the MC atoms, i.e., the hydrogens located in the MGs show higher mobility than the MC hydrogens. Additional mobility associated to the MG dynamics with

52

F. Alvarez et al. / Chemical Physics 261 (2000) 47±59

respect to the dynamical processes in the rest of the protons can thus be observed. The mean squared displacement obtained for both kinds of atoms (Eq. (3) with n ˆ 1) is displayed in Fig. 3(a) as a function of time. In this ®gure, the qualitative features of the motion deduced by inspection of the probability distributions can be analysed in more detail. The hr2 …t†i of the MC protons increases almost linearly with time for t < 1 ps approximately, in the so-called regime of the fast dynamics. At this point, some cage effect prevents the atom to continue moving in that way and hr2 …t†i increases much slower until a time

Fig. 3. Time dependence of (a) the mean squared displacement and (b) the non-Gaussian parameter for the MC protons ( ) and for the MG protons (r). Dotted line in (a) shows the mean squared displacement for the carbons of the MGs.

of around 10 ps. After this time, a continuous sublinear dependence of hr2 …t†i on t can be found, which can be attributed to the a-relaxation [18]. Concerning the protons in the MG, it is clear that their hr2 …t†i is higher than that of the MC protons in the whole time range studied. However, an approximately parallel behaviour for hr2 …t†i of both subsystems is found for longer times, where the a-relaxation dominates the dynamics. Let us now concentrate on the shape of the van Hove correlation functions. As an example, we have included as solid lines in Fig. 2 the resulting curves for t ˆ 4 ps when the Gaussian approximation is assumed. These curves have been calculated from Eqs. (2) and (4) with the corresponding values of hr2 …t ˆ 4 ps†i found in each subsystem. Whereas for the MG protons, the Gaussian form for the self correlation function would provide a rather good approximation for the results at 4 ps, it is clear that strong deviations from this behaviour can be found in the case of the MC protons for the same value of the time. In order to quantify the possible deviations from the Gaussian approximation, the values of the nonGaussian parameter a2 …t† (Eq. (8)) were calculated for each subsystem. They are shown in Fig. 3(b). Two clear maxima can be found in a2 …t† for both cases, and, as it was expected, the limit of this quantity for t ! 0 is 0. For the MC protons, the ®rst maximum is located in the region of the fast dynamics and the second one in the time range where the so-called decaging processes (transient from ``fast'' to ``slow'' dynamics) take place. It is worthy of remark that, as soon as the sublinear dependence of hr2 …t†i is established, the values of a2 …t† decrease. This implies that the behaviour of the van Hove correlation function becomes closer and closer to Gaussian with increasing time in the region of the a-relaxation. A similar maximum in a2 …t† to that found for the MC atoms at shorter times is observed in the case of the protons in the MGs. However, the second maximum in the latter case is much less pronounced and shifted towards longer times. It is worth noting that the values of a2 …t† for both kinds of atoms are very close in the range where their hr2 …t†i run parallel, i.e., where the a-relaxation plays the most important role. The main di€erences in the non-Gaussian para-

F. Alvarez et al. / Chemical Physics 261 (2000) 47±59

meter are obviously found in the intermediate time range. In fact, it is at times around 4 ps ± the example depicted in Fig. 2 ± when the di€erence between the values of a2 …t† for the two subsystems reaches its maximum. The intermediate scattering functions were calculated through Fourier transformation (Eq. (5)). In Fig. 4, the results of these calculations for Q ˆ 1 ÿ1 are presented. It is worthy of remark that the A accuracy of the calculated values of Fs (Q,t) at longer times, close to the limit of the MD simulation run trun (trun ˆ 2 ns in our case) depends on the number of time origins which are available for calculating the ensemble average of the van Hove correlation function. To have an estimation of the uncertainties involved, we have followed the procedure described for the example in Ref. [19]. The standard error at a time t is given by r  1=2 ‰2sc =…N …trun ÿ t††Š , where sc is a ``¯uctuation'' correlation time and, N, the number of particles considered. In the error estimation, we have taken

53

sc of the order of the time at which Fs …Q; t†  1=e. The estimated errors are shown in the inset of Fig. 4 for the case of the MC protons where N ˆ 500 and sc  100 ps. Errors of the same order are obtained for the other curves shown in Fig. 4. As can be seen in Fig. 4, at shorter times, a more pronounced decay is observed for FsMG …Q; t† than for FsMC …Q; t†. We will ®rst compare our simulation results to previous neutron scattering results on PI [20]. In the experiments reported in that work, a fully protonated sample was investigated at 340 K by means of backscattering (BS) techniques, covering a dynamical range of approximately 5 ps < t < 2 ns and a Q-interval of 0.2 ÿ1 . At the temperature investiÿ1 < Q < 1:9 A A ÿ1 , the agated, and for Q-values below Q  1 A relaxation dominates the experimental spectra, that were described according to the expression (  b ) t Fs …Q; t† ˆ A…Q; T † exp ÿ : …9† sKWW …Q; T † This equation describes the decay due to the arelaxation by means of the Kohlrausch±Williams± Watts (KWW) function and gives account for the faster processes ± vibrations, MG dynamics, ... ± through the prefactor A…Q; T †. The parameters characterizing the a-process are thus the characteristic time sKWW …Q; T † and the stretching parameter b. For PI, a value of b ˆ 0:40 was ÿ1 , reported and the value of sKWW …Q ˆ 1 A T ˆ 340 K† can be estimated to range between 15 and 50 ps from Fig. 3(b) in that work. For a qualitative comparison with these neutron scattering results, we have calculated the intermediate scattering function which would be observed in the case of a fully protonated sample, Fstot …Q; t†. This curve can be easily constructed starting from FsMG …Q; t† and FsMC …Q; t†, taking into account that in each monomer, three out of eight protons are located in the MG and the rest in the MC:

Fig. 4. Incoherent intermediate scattering function obtained from the MD simulations for the MC protons ( ), for the MG protons (r), and for all protons (e). Solid lines correspond to Eq. (9) with b ˆ 0:40 and sKWW ˆ 120 ps and A ˆ 0:96, 0.84, and 0.65 (top to bottom) in the time range usually accessed by BS spectrometers. Dashed lines show the extrapolations to shorter times. The inset shows a magni®cation of the long time range of FsMC …Q; t† including the estimated error bars.

Fstot …Q; t† ˆ 58FsMC …Q; t† ‡ 38FsMG …Q; t†:

…10†

The resulting function has been included in Fig. 4. As one could expect, the decay of this function takes place in between the other two curves. We followed the same procedure for analysing Fstot …Q; t† as that used for the neutron scattering

54

F. Alvarez et al. / Chemical Physics 261 (2000) 47±59

data. A ®t of Fstot …Q; t† to Eq. (9) with a ®xed value of b ˆ 0:40 gives a reasonable description of the function for t longer than 5 ps, i.e., in the time range accessible by BS. This can be seen in Fig. 4. The time scale obtained from the ®t is sKWW …Q ˆ 1 ÿ1 , T ˆ 363 K† ˆ 120 ps, slightly slower than A that experimentally obtained. Taking into account that at such high temperatures above Tg , the apparent activation energy of sKWW is quite low, we can ®nally estimate a ratio of about 6 between the ÿ1 , T ˆ 363 K† as obtained value of sKWW …Q ˆ 1 A from our simulations and that extrapolated from the BS experiments. From this comparison, one can thus conclude that the simulations seem to give results for the a-relaxation which are quite close to the experimental ones. Concerning the prefactor A…Q; T †, a value of 0.84 is found in the ÿ1 . This value is case of the simulation for Q ˆ 1 A close to that estimated for the same conditions from the neutron scattering results [0:76  0:06] [21]. Starting from the description obtained for the a-relaxation in the ``fully protonated'' sample, we can now check whether the same behaviour is observed in the case of the MC protons and of the MG protons in this region. In Fig. 4, we have included the ®ts of FsMG …Q; t† and FsMC …Q; t† with Eq. (9), in which only the prefactor A…Q; T † has been allowed to vary. The values of sKWW and b have been ®xed to those previously used for describing Fstot …Q; t†. As can be seen, the agreement between the results of the simulations and the theoretical curves is quite acceptable in the BS region, although there are some deviations at longer times, in particular in the case of the MG protons. In fact, a free-®tting procedure gives, in the latter case, a b-value close to 0.5. In spite of these differences, these results suggest that both kinds of protons in the sample participate in the segmental relaxation in a similar way. The main di€erence in the dynamical behaviour seems to show up at shorter times, in the picosecond region. In this, as commented above, the MC protons undergo the fast dynamics crossing over to the a-relaxation after a decaging process; the MG protons, probably in addition to that dynamical processes, rotate with a characteristic time that lies in this time regime.

4. Discussion The ®nding of a similar behaviour in the range of the a-relaxation for both kinds of protons in the sample suggests that the rotational motion of the protons in the MG can be considered to a good approximation as decoupled from the slower structural a-relaxation process. This can be taken as a starting point for understanding the dynamics of the MG protons also at shorter times. Assuming that also in the short time regime the segmental motion and the MG rotation are independent dynamical processes, FsMG …Q; t† can be written as the product of the incoherent scattering function corresponding to a purely rotational motion Fsrot …Q; t† and the incoherent scattering function giving account for the segmental motions Fsseg …Q; t†, FsMG …Q; t† ˆ Fsrot …Q; t†Fsseg …Q; t†:

…11†

In the ®rst approach, we can identify Fsseg …Q; t† with the incoherent scattering function calculated for the protons in the MC, i.e., Fsseg …Q; t†  FsMC …Q; t†:

…12†

In this way, we are facing the problem as an experimentalist would do, since FsMC …Q; t† can be accessed by means of neutron scattering, performing measurements on samples with deuterated MGs. The expressions (11) and (12) allow to obtain in a straightforward way the function Fsrot …Q; t† by simple division of the two known functions. If our assumption of statistically independent rotation and segmental dynamics is correct, consistent results should be obtained for Fsrot …Q; t†. The resulting curves of such a ``deconvolution'' procedure are shown for di€erent Q-values in Fig. 5 up to 20 ps. At ®rst sight, it can be deduced that the behaviour found qualitatively ®ts to that expected for a rotational motion: (i) there is no apparent Q-dependence of the characteristic time of the decay and (ii) the correlation functions present clear elastic contributions, i.e., Fsrot …Q; t† get constant values, which depend on Q, at longer times. We checked whether these elastic levels corresponded to the elastic incoherent structure factor

F. Alvarez et al. / Chemical Physics 261 (2000) 47±59

Fig. 5. Incoherent intermediate scattering function for the rotation of the MG protons for di€erent Q-values, obtained under the assumption of being independent of the segmental motions. Solid lines show the ®t to the RRDM with the distribution function of potential barriers shown in Fig. 7 as a solid line. The ÿ1 data to a single dotted line corresponds to the ®t of the 1.6 A barrier model.

(EISF) expected for a rotation of the proton in a threefold potential. This is given by [22]   1 sin …QrHH † : …13† EISF ˆ 1 ‡ 2 3 QrHH  for the H±H distance rHH in with a value of 1.78 A a MG. Fig. 6 shows the comparison between this theoretical expression and the EISF obtained from the values of Fsrot …Q; 10 ps < t < 20 ps†, where the function seems to remain constant. As can be seen in this ®gure, the agreement between these data is excellent. This ®nding gives strong support to the assumption of statistic independent rotation and segmental motions in this time range. At times longer than 20 ps, the deconvolution procedure described above produces inconsistent results (the curves eventually bend upwards or downwards depending on the Q-value). We can envisage two di€erent origins for such a behaviour: (i) the statistical independence of the rotational and segmental motions fails as a consequence of some coupling between these dynamics; (ii) di€erences between the segmental dynamics at longer times

55

Fig. 6. Momentum transfer dependence of the EISF obtained from the long time limit of Fsrot …Q; t† compared to the theoretical prediction for a rotation in a threefold potential (Ð).

(a-relaxation) of the protons in the MC and protons in the MG, i.e., Eq. (12) would not be a good approximation. We will come back to this point when a quantitative characterization of the rotational contribution is done. In order to perform a quantitative analysis of Fsrot …Q; t† including the inelastic contribution, let us remember the expression for the incoherent intermediate scattering function of a rotational motion with an activation energy barrier Ea (see, e.g., Ref. [2]): FsEa …Q; t† ˆ EISF ‡ …1 ÿ EISF† 2 3 t   5;  exp 4 ÿ s0 exp kEBaT

…14†

where kB is the Boltzmann constant and s0 ˆ …2pm0 †ÿ1 (m0 : attempt frequency of the transition). In a threefold potential, the EISF is given by Eq. (13). Expression (14) implies a single exponential decay for the inelastic contribution to the scattering function. It is clear that in the obtained Fsrot …Q; t†, as already found in real polymers [1±4], this is not the case. This is evidenced in Fig. 5 where we have included as an example the ®t to Eq. (14) (dotted line) of the data corresponding to

56

F. Alvarez et al. / Chemical Physics 261 (2000) 47±59

ÿ1 . As commented in Section 1, the RRDM 1.6 A allowed to describe the Fs …Q; t† corresponding to MG rotations measured by neutron scattering in a series of polymers. In this framework, the existence of a distribution of potential barriers (or activation energy barriers) G…Ea † in the sample is assumed. The total scattering function is thus given by the superposition of the di€erent contributions weighted by the distribution function Z 1 RRDM …Q; t† ˆ EISF ‡ …1 ÿ EISF† G…Ea † Fs 0 3 2 t   5 d…Ea †:  exp 4 ÿ s0 exp kEBaT …15† The Fsrot …Q; t† found were analysed in terms of the RRDM. According to the assumptions of the RRDM, in the ®rst approximation, a Gaussian function could be chosen for describing G…Ea †, "  2 # 1 1 Ea ÿ E0 G…Ea † ˆ p exp ÿ ; …16† r 2 r 2p where E0 is the most probable value of the activation energy and r gives the width of the distribution. The number of free parameters in this description is three: E0 , r and s0 . Note that the parameters E0 and s0 are coupled if only one temperature is analysed, as this is the case. We have therefore taken advantage of the information about the distribution function G…Ea † provided by the low temperature results of the MD simulations on the torsional libration density of states [6]. The distribution G…Ea † at 150 K can easily be obtained from the distribution of MG torsional librations g…E† given in the previous papers [5,6], as G…Ea † ˆ g…E†dE=dEa . In the framework of the threefold approximation and assuming that the simulated (or measured) density of torsional states only corresponds to the ®rst librational state, a numerical relationship between the librational energy E and the activation energy barrier Ea can be established from the solutions of the corresponding Schr odinger equation. Tabulated solutions of this equation ± which for a pure sinusoidal threefold potential takes the form of the well-known Mat-

hieu equation ± can be found for example in Ref. [23] for di€erent values of the potential barrier V3 . These numerical solutions allow one to correlate the values of Ea and E. The correlation can be well described by a general power law like Ea …meV†  a‰ E …meV†Ša ÿ Ea0 …meV†, where the values of the parameters a, a and Ea0 depend on the range considered. In the range of Ea which is relevant for this work (E < 42 meV approximately), the values found for these parameters are a ˆ 0:182, a ˆ 1:981, and Ea0 ˆ 0. The numerical correlation between E and Ea as well as its description by the power law with such parameters are shown in Fig. 7 in the mentioned Ea interval. From this analytical function, dE=dEa and thereby G…Ea † can be easily calculated. The distribution function G…Ea † obtained by means of this procedure is shown in Fig. 7. This ®gure also includes the same distribution G…Ea † obtained by the same procedure but from the density of torsional states measured by inelastic neutron scattering in PI [5,10]. Both of them show the maximum value for Ea  100 meV. Therefore, in order to ®t the high temperature behaviour of Fsrot …Q; t† to Eq. (15), we have ®xed the value of E0 in Eq. (16) to 100 meV. This ®rst approximation is based on the idea that the average activation energy for MG rotation in polymers is mainly determined by the chemical structure of the monomer [2] and the very local non-bond interactions [6]. In our previous works [5,6], and in performing MD simulations under ``phantom-conditions'' (i.e., when the non-bond interactions are switched o€), we estimated that the contribution of the non-bond interactions to the mean value of the potential barrier was of the order of 30% for PI. We have checked (see later) what would be the in¯uence of such an uncertainty for E0 on the ®nal value found by ®tting for r. The curves Fsrot …Q; t† shown in Fig. 5 were then ®tted to Eq. (15) with Eqs. (13) and (16) and E0 ˆ 100 meV. Therefore, only the values of r and s0 were allowed to change. The resulting ®tting curves are shown in Fig. 5 as solid lines. They reproduce the Fsrot …Q; t† very nicely. Concerning the values of the parameters, we found s0 ˆ 1:9  10ÿ2 ps, corresponding to an attempt frequency of m0 ˆ 8  1012 Hz, and r ˆ 37 meV. The thus deduced distribution function G…Ea † is displayed in

F. Alvarez et al. / Chemical Physics 261 (2000) 47±59

Fig. 7. Distribution function of potential barriers obtained from the torsional libration density of states calculated by MD simulations [6] (


) and measured by inelastic neutron scattering [10] ( ). The solid line shows the distribution function deduced in this work from the ®ts of the high temperature simulation data to the RRDM (see text). Inset: correlation between E and Ea obtained from the numerical solutions of the Mathieu equation in the range of E relevant for this work. The solid line shows the description by the power law given in the text.

Fig. 7 together with the results from the MD simulations at low temperature [6] and from inelastic neutron scattering measurements [10]. As it can be seen, the agreement among these functions is quite good. The main di€erence between the distribution function G…Ea † found at 363 K and those corresponding to 150 K is the high energy tail shown by these last. In order to rationalize this di€erence, it is worthy of remark that the low temperature G…Ea † are calculated from the distribution of torsional librations g…E† obtained either from MD simulations or inelastic neutron scattering measurements. As it was described above, the transformation from g…E† to G…Ea † is based on the assumption that g…E† only corresponds to the ®rst librational state. However, second-order contributions cannot be discarded at 150 K which would be visible in the high frequency range of g…E†. These contributions would produce an arti®cial high energy tail in the calculated G…Ea †. On the other hand, it is worth noting that the value

57

found for r seems to be almost insensitive to possible changes of the value of E0 caused by a softening of the non-bond interactions at high temperature. Assuming the extreme case that these interactions become negligible, G…E† would be centered at about 70 meV (the estimated 30% contribution from the non-bond interactions to the activation energy is subtracted from the total value). Under these conditions, the value obtained by ®tting for r is 39 meV, i.e., it only di€ers by 5% from that reported above corresponding to E0 ˆ 100 meV. These ®ndings seem to indicate that the distribution function of potential barriers for MG rotation in polymers is not very sensitive to the glass transition. This result is in qualitative agreement with the idea suggested by our previous results [6] that the width of this distribution is mainly controlled by the very local non-bonded interactions. At this point, it is worthy of remark that the reduction of the value of the non-Gaussian parameter in the MG protons with respect to that of the MC protons in the long time range mentioned above ®nds a natural explanation in the framework proposed here. In the case of the combination of two independent motions (let us denote with superscripts a and b), the parameter a2 can be calculated as a function of the non-Gaussian parameters (aa2 and ab2 ) and the mean squared displacements …hr2 …t†ia and hr2 …t†ib † corresponding to the processes involved [24]. It can be shown that if aa2 takes vanishing or negative values, as in the case of a con®ned motion in a certain volume at long times, then the value of a2 is smaller than that of ab2 . This is just the situation found for the MG protons, where at the temperature investigated the non-Gaussian parameter corresponding to the rotational contribution vanishes at 2 ps approximately and reaches an asymptotic limit of ÿ0.1 for times longer than 10 ps. Thus, the behaviour observed for a2 can be understood, at least qualitatively, supporting the scenario proposed. However, it has also to be mentioned that for t > 20 ps the values calculated from the simulation for hr2 i of the MG protons are larger than those obtained by the addition of the corresponding ones to the pure rotational motion and the MC-protons. This ®nding is related to the inconsistencies above

58

F. Alvarez et al. / Chemical Physics 261 (2000) 47±59

commented that are observed in this time range when applying the deconvolution procedure given by Eqs. (11) and (12). In the following, the sources of these discrepancies will be discussed. The ®rst origin we have mentioned in advance is related to a possible coupling between rotation and segmental dynamics. In order to quantify this coupling, we have calculated the so-called connected intermediate scattering function, Fcon …Q; t†, in a similar way as it was introduced in Ref. [25] for water. In our case, this function, which describes the strength of the coupling, can be de®ned as Fcon …Q; t† ˆ FsMG …Q; t† ÿ Fsrot …Q; t†FsMC …Q; t†:

…17†

The functions involved in this calculation are ÿ1 (close to depicted in Fig. 8(a) for Q ˆ 1:4 A the ®rst maximum of S(Q)). For this Q-value, the function constructed assuming no coupling, Fsrot …Q; t†FsMC …Q; t†, is very close to FsMG …Q; t†, leading thus to very small values for Fcon …Q; t†. This is the case for most of the Q-values investigated, as can be deduced from Fig. 8(b). The maximum value of Fcon …Q; t†…0.1† is found at longer times (approximately for t > 20 ps) in the intermediate Q range. Therefore, the limit of validity of this approximation could be set around 20±30 ps, when Fcon …Q; t† ceases to be negligible for all Q-values. Note that negative values of Fcon …Q; t† correspond to faster decays of FsMG …Q; t† than that of Fsrot …Q; t†FsMC …Q; t†. The other source of uncertainties at longer times could be the deviations from the assumption expressed by Eq. (12). The purely segmental component of the MG protons dynamics is approached by that of the protons in the MC. As mentioned before, this would be the procedure followed by a neutron scattering experimentalist when analysing and interpreting the data. However, in principle, a more precise approximation could be to consider the intermediate scattering function of the carbon in the MG, which is not accessible by neutron scattering. Simulation allows nevertheless to follow the dynamical behaviour of these C-atoms and to compare it with that of the MC protons used in the deconvolution procedure. A qualitatively similar dynamical behaviour for

ÿ1 : Fig. 8. (a) Intermediate scattering function at Q ˆ 1:4 A MG Simulation results for MG protons, Fs …Q; t† (r) and MC protons, FsMC …Q; t† (O); calculated rotational contribution (±


±), the product, Fsrot …Q; t†FsMC …Q; t† (Ð) and the connected function, Fcon …Q; t† (‡). (b) Time dependence of Fcon …Q; t† for several Q-values investigated. Symbols are as in Fig. 5.

both kinds of atoms is observed, as can be appreciated in Fig. 2(a). The C-atoms however move slightly slower than the MC protons below 100 ps and faster above 100 ps. Therefore, in this long time regime the description of the MG dynamics, though not yet perfect, would improve by approaching the segmental part of the motion by the dynamics observed for the C-atom in the MG. For example, at 1 ns, the value of the mean squared displacement deduced for the MG proton in this

F. Alvarez et al. / Chemical Physics 261 (2000) 47±59

2 when the MC protons case would be 19 and 16 A are considered (the value obtained in the simula2 ). tion is 22 A Thus, the limits of the decomposition proposed here for the motion of the MG protons could be set around 20±30 ps, being the origin of this limitation either due to a small coupling between the processes or to di€erences between the segmental dynamics of MC protons and MG protons. Despite this ®nding, this scenario could be considered as a good approximation when dealing with experimental data corresponding to neutron scattering, taking into account the uncertainties usually involved in such experiments.

5. Conclusions From the results of MD simulations of PI at 363 K, we have compared the dynamical behaviour of the protons located in the MG and in the MC. It can be concluded that (i) the rotational and the segmental dynamics are decoupled to a good approximation and (ii) the energy landscape for the MG motion seems to be una€ected by the glass transition.

Acknowledgements We acknowledge support from the following projects: DGICYT, PB97-0638; GV, EX 1998-23; UPV/EHU, 206.215-G20/98. Support from ``Donostia International Physics Center'' is also acknowledged.

59

References [1] J. Colmenero, R. Mukhopadyay, A. Alegrõa, B. Frick, Phys. Rev. Lett. 80 (1998) 2350. [2] R. Mukhopadyay, A. Alegrõa, J. Colmenero, B. Frick, Macromolecules 31 (1998) 3985. [3] A. Moreno, A. Alegrõa, J. Colmenero, B. Frick, Phys. Rev. B 59 (1999) 5983. [4] A. Chaid, A. Alegrõa, J. Colmenero, Macromolecules 27 (1994) 3282. [5] F. Alvarez, A. Alegrõa, J. Colmenero, T.M. Nicholson, G.R. Davies, in: M.R. Johnson, G.J. Kearley, H.G. B uttner (Eds.), CP479 Neutrons and Numerical Methods ± N2 M, American Institute of Physics, New York, 1999, p. 201. [6] F. Alvarez, A. Alegrõa, J. Colmenero, T.M. Nicholson, G.R. Davies, Macromolecules, accepted for publication. [7] N.E. Moe, M.D. Ediger, Polymer 37 (1996) 1787. [8] N.E. Moe, M.D. Ediger, Phys. Rev. E 59 (1999) 623. [9] F. Alvarez, et al., submitted for publication. [10] B. Frick, L.J. Fetters, Macromolecules 27 (1994) 974. [11] B. Frick, D. Richter, Cl. Ritter, Europhys. Lett. 9 (1989) 557. [12] L. van Hove, Phys. Rev. 95 (1954) 249. [13] R.J. Glauber, Phys. Rev. 87 (1952) 189. [14] R.J. Glauber, Phys. Rev. 94 (1954) 751. [15] R.J. Glauber, Phys. Rev. 98 (1955) 1692. [16] V.F. Sears, Neutron News 3 (1992) 26. [17] A. Rahman, K.S. Singwi, A. Sj olander, Phys. Rev. 126 (1962) 986. [18] J. Colmenero, A. Alegrõa, A. Arbe, B. Frick, Phys. Rev. Lett. 69 (1992) 478. [19] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1997. [20] A. Arbe, J. Colmenero, M. Monkenbusch, D. Richter, Phys. Rev. Lett. 81 (1998) 590. [21] S. Ho€mann et al., private communication. [22] M. Bee, Quasielastic Neutron Scattering, Adam Hilger, Bristol, 1988. [23] M. Prager, A. Heidemann, Chem. Rev. 97 (1997) 2933. [24] R. Zorn, Phys. Rev. B 55 (1997) 6249. [25] S.-H. Chen, P. Gallo, F. Sciortino, P. Tartaglia, Phys. Rev. E 56 (1997) 4231.