Disentangling Self-Atomic Motions in Polyisobutylene by Molecular Dynamics Simulations

We present fully atomistic molecular dynamics simulations on polyisobutylene (PIB) in a wide temperature range above the glass transition. The cell is validated by direct comparison of magnitudes computed from the simulation and measured by neutron scattering on protonated samples reported in previous works. Once the reliability of the simulation is assured, we exploit the information in the atomic trajectories to characterize the dynamics of the different kinds of atoms in PIB. All of them, including main-chain carbons, show a crossover from Gaussian to non-Gaussian behavior in the intermediate scattering function that can be described in terms of the anomalous jump diffusion model. The full characterization of the methyl-group hydrogen motions requires accounting for rotational motions. We show that the usually assumed statistically independence of rotational and segmental motions fails in this case. We apply the rotational rate distribution model to correlation functions calculated for the relative positions of methyl-group hydrogens with respect to the carbon atom at which they are linked. The contributions to the vibrational density of states are also discussed. We conclude that methyl-group rotations are coupled with the main-chain dynamics. Finally, we revise in the light of the simulations the hypothesis and conclusions made in previously reported neutron scattering investigations on protonated samples trying to address the origin of the dielectric β-process.


Introduction
During the last decades, neutron scattering (NS) has proved to be an extremely valuable microscopic technique to decipher the structural and dynamical properties of glass-forming systems, including polymeric materials [1][2][3]. In particular, its sensitivity to hydrogen motions allows characterizing the molecular dynamics of protonated systems at length scales in the Å to nm range. At such scales, well above the glass-transition temperature T g the structural or α-relaxation is the most relevant dynamical process. In the early 90's, by using backscattering techniques on protonated samples [poly(vinyl methyl ether) (PVME), phenoxy (PH) and poly(vinyl chloride) (PVC)], the correlation function corresponding to self-atomic (hydrogen) motions in this regime was characterized by stretched exponentials with similar stretching exponents (β-parameter) as those obtained by dielectric relaxation for the α-process [4]. The momentum transfer (Q) dependence of the characteristic times was found to be correlated with the shape parameter β such that the atomic probability distributions could be approximately described by Gaussian functions. The underlying motions were qualified as anomalous diffusion-like processes where the atomic mean squared displacement r 2 increases sublinearly with time ( r 2 ∝ t β ). These results were later confirmed by similar experiments on 'canonical' polymers [like e.g., 1,4-polybutadiene (PB), poly(vinyl ethylene) (PVE), polyisoprene (PI) and polyisobutylene (PIB)] [5], and also on other systems like ionic liquids [6,7], alcohols [8], nanocomposites [9,10], hydration water [11] or other complex macromolecules like lignin [12]. Backscattering measurements are however restricted to a rather limited dynamic window and the results are affected by the instrumental resolution function, which hampers an accurate determination of the spectral shape. Deconvolution from the instrumental resolution is possible by the neutron spin echo (NSE) technique, which in addition accesses the widest dynamic range for NS instruments. The first NSE results on a fully protonated polymer were reported in Ref. [13]. The sample was PIB. Though still subjected to uncertainties in the determination of the exact value of the shape parameter β, the NSE data revealed unambiguously stretched functional forms and Q-dependences of the characteristic times compatible with the Gaussian approximation, confirming the previous general conclusions. In the high Q-range investigated, however, indications for deviations from Gaussian behavior were found, suggesting a crossover from Gaussian to non-Gaussian behavior in a Q-range of about 1 Å −1 .
Such deviations were investigated in more detail on PI with deuterated methyl groups by combining NS (backscattering and NSE) experiments with molecular dynamics (MD) simulations [14,15]. This combined approach results to be an extremely powerful tool, as it has been shown in a series of systems during the last years (see, e.g., [3,16]). To explain the behavior found for the hydrogen motions in PI, the so-called anomalous jump diffusion model was proposed. This model is a generalization of the well-known jump-diffusion model [17] to the case of sublinear diffusion. Similar crossover from Gaussian to non-Gaussian behavior was also later reported for other glass-forming polymers like PB [18], PVE [18,19], PVME [20], poly(methyl methacrylate) (PMMA) [21], poly(ethylene propylene) (PEP) [22], head-to-head polypropylene (hh-PP) [22], poly(vinyl acetate) (PVAc) [23], poly(vinyl pirrolidone) (PVP) [24], poly(ethyl methacrylate) (PEMA) [25] and poly(tetrahydrofurane) (PTHF) [26], as well as for other systems as e.g., n-alkanes [27] or H-bonding liquids [28][29][30]. The anomalous jump diffusion model has also been applied to some of these polymers and seems to be successful in capturing the essence of the caging effect occurring in these systems in a very simple way. The limited Q-range of the NSE experiments on PIB [13] prevented however a reliable application of such model.
In addition to the α-relaxation, other processes like the β-relaxation might be found in polymers. In two cases, PB and PIB, there were attempts to characterize by NS the molecular motions behind the β-process. This is an extremely difficult task, since in the relatively high frequency window accessed by NS, contributions not only from the β-process but also from the α-relaxation are expected to be present. For PB, coherent scattering data obtained on a deuterated sample were analyzed assuming a scenario of statistically independent α and β-processes [31,32]. The localized motions behind the β-relaxation could be characterized as rotational jumps of the chain building blocks around their center of mass. This was corroborated by MD-simulations some years later [33]. Such an approach was also applied to describe the NS results on PIB. From the description of the self-correlation function accessed on a protonated sample, the length of the jumps involved in the localized motions giving rise to the dielectric β-process was determined to be of about 2.7 Å [34]. However, NSE experiments on a fully deuterated sample addressing collective relaxation in PIB and evaluated in the same framework delivered a much smaller value (0.5-0.9 Å) for such jump distance [35]. Conversely, the dielectrically detected βrelaxation showed very similar features as the so-called δ-process reported from NMR experiments [36][37][38], that was interpreted as rotations of methyl groups (which are not dielectrically active). Interestingly, the NS spectra of PIB in the vicinity of T g do not show any clear hint of classical hopping of its methyl groups [39,40], as it is usually the case in glass-forming polymers [41]. The NS experimental signature of methyl group dynamics in PIB is the presence of a pronounced, broad and structured peak in the vibrational density of states (VDOS) measured on fully protonated samples [39,42]. Such data in fact reveal the total-hydrogen averaged VDOS, containing contributions from both, methyl-group and main-chain hydrogens. The reported peak is located at an energy of about 40 meV, i.e., a rather high value as compared to the typical energies corresponding to the methyl-group librational peaks in polymers [41]. This observation was explained in Ref. [39] in terms of a strong coupling between the methyl groups, which also may be related to the strong steric hindrance of the polymer chain in PIB. Adams et al. [42] proposed that the peak at ≈40 meV was due to a broadened methyl group torsional frequency, and the structure observed at 25 K (a kind of shoulder or additional peak at higher frequencies) was due mainly to the C-C-C skeletal vibrations being in the same region of the spectrum. They remarked the importance of using selective deuteration to experimentally resolve possible different components of the structured broad peak observed.
Polyisobutylene is thus one of the glass-forming polymers that has been most extensively investigated by NS [13,34,35,39,40,[42][43][44][45][46][47][48][49][50][51][52][53][54]. At the same time, however, there are still a number of open questions regarding the interpretation of these results. In addition, we cannot forget that PIB is a very interesting polymer also for industrial applications due to its well-known and unusual physical properties; remarkably low T g ≈ 200 K resulting in low cost processability, high packing efficiency and low gas permeability. Any contribution to the microscopic understanding of the underlying dynamical processes in this system is always important.
With these ideas in mind, we have performed fully atomistic MD-simulations on PIB in a wide temperature range above T g . We have followed the strategy that has proved to be highly successful in previous studies that combined NS and MD-simulations [3,16]. First, the cell is validated by direct comparison of magnitudes computed from the simulation and measured by NS. We addressed problems related with collective features in a previous work [55]. There, the realism of the simulated cell was checked by comparison of data corresponding to the static and dynamic structure factor. Here we focus on self atomic motions. Consistently, the validation will be performed against NS results on protonated samples mainly determined by the self-correlation function of the hydrogens. Once the reliability of the simulated cell is assured, with all the information contained in the simulations at hand, we can ask questions that an experimentalist cannot answer due to technical limitations. In particular, we can address how good are the approximations involved in the experimental data analysis (always a crucial problem).
This work is focused on the characterization of the dynamics of the different kinds of atoms in PIB. In theory, partial deuteration of the sample allows selectively investigating by NS the motions of the remaining hydrogens. However, partially deuterated PIB samples are not easy to synthesize and consequently there are no available NS results that allow isolating main-chain from methyl-group hydrogen motions. Moreover, since the incoherent scattering cross-section of carbon atoms is 0, their self-motions are inherently inaccessible by NS. One particularly interesting question is how different are the motions of the total subensamble of hydrogens from those of the main-chain atoms. Do the latter also show deviations from Gaussian behavior? Along this paper, we first discuss the dynamic heterogeneity associated to the different atomic species and the applicability of the anomalous jump diffusion model to the different kinds of atoms in PIB. Thereafter, the characterization of the methyl-group hydrogen motions is performed. We note that these motions are very commonly found to contribute the the quasielastic spectra together with other processes like the segmental relaxation (see, e.g., [56][57][58][59][60][61][62]) and these contributions are not easy to disentangle. We show that the usually assumed statistically independence of rotational and segmental motions is inappropriate for PIB. Exploiting the possibility offered by the simulations of calculating correlation functions corresponding to relative motions of atoms (including those defined in real space), we apply the rotational rate distribution model and discuss the particularities of PIB methyl-group dynamics. The contributions to the VDOS are also disentangled. Finally, we revise in the light of the simulations the hypotheses and conclusions made in previously reported neutron scattering investigations on protonated samples trying to address the origin of the dielectric β-process of this polymer.

Molecular Dynamics Simulations
Fully atomistic molecular dynamics simulations were carried out using the COMPASS forcefield implemented within the commercial software package MATERIALS STUDIO 5.1. The cubic cell contained 20 PIB chains of 70 monomers (M w = 3923 g/mol, i.e., smaller than the entanglement mass M e ≈ 7000 g/mol [63], and a total number of atoms N = 16,840). The system was simulated at the temperatures of 470, 390, 365, 335 and 320 K, following the procedure explained in Ref. [55], to carry out 100 ns dynamics.

Computed Magnitudes
The trajectory in space of a single atom is described by the self part of the van Hove correlation function G s ( r, t). G s ( r, t) is the probability to find an atom at time t at a position r if it was at r=0 for t=0: Here r i is the position of the atom i and N the number of nuclei in the sample. The brackets denote the ensemble average. The Fourier transform of G s ( r, t) in Q-space is the so-called intermediate scattering function F s ( Q, t), For some simple cases-free nuclei in a gas, harmonic crystals, simple diffusion at long times; G s ( r, t) is a Gaussian function [64,65]; in an isotropic system this implies In the Gaussian approximation the even moments of G s (r, t) can be straightforwardly calculated. For instance, the mean squared displacement of the atom r 2 (t) is given by r 2 (t) = 3/[2α(t)]. Moreover, in such Gaussian case F s (Q, t) is entirely determined by r 2 (t) : In more general cases, deviations of G s (r, t) from the Gaussian form [Equation ((3)] may be expected. F s (Q, t) can then be expressed in terms of its expansion in Q (see, e.g., Ref. [66]) where α 2 (t) giving the leading correction is the so called second order non-Gaussian parameter. It is defined as [67] α 2 (t) = 3 5 Evidently, in the Gaussian case α 2 (t) = 0.

Validation
As we will compare MD-simulation with experimental neutron scattering (NS) results, we will introduce some magnitudes and functions related with NS that will be invoked. In general the neutron intensity scattered by a sample contains an incoherent contribution reflecting the self motion of atoms, as above defined, and a coherent contribution related to the atomic pair correlations. These contributions are weighted by the factors I inc and I coh (Q) respectively. These factors depend on the incoherent and coherent scattering lengths b inc i and b coh i of the nuclei in the sample (i refers to the nucleus considered) as: and Here, r i and r j are taken at the same time and the brackets mean the thermal average. Hydrogen presents a very large value for the incoherent scattering length (b inc H = 25.274 fm). In a fully protonated polymer sample consisting of H and C, the scattered intensity is completely dominated by the incoherent contribution from the protons (C only scatters coherently with b coh C = 6.6511 fm, and b coh H = −3.7406 fm). Therefore, in such case I inc >> I coh (Q) for all Q-values and NS results are predominantly determined by F s (Q, t) related to the self motion of the hydrogens.
In this work, to validate the simulated cell we have used previously published neutron spin echo (NSE) results on fully hydrogenated PIB [13]. NSE is distinguished from all other dynamic neutron scattering techniques in that it measures time-dependent correlation functions. The NSE signal is given by [68,69]: whereF(Q, t) and F s (Q, t) are the intermediate pair correlation function (normalized to its value at t = 0) and the self correlation function above defined respectively. Figure 1 shows the comparison between the experimental results measured by NSE on a fully hydrogenated sample and computed from the atomic trajectories in the simulation following Equation (10). The agreement is excellent. The simulations fairly reproduce the Qand T-dependence of the characteristic times and spectral shape experimentally observed. It is noteworthy that no shift in timescale has been imposed for a direct matching.
Though the main focus of this work is on the dynamical behavior at times longer than the picosecond, we will also discuss some aspects related with the vibrational processes. Therefore, we have also considered the vibrational density of states (VDOS, Z(E)) to validate the simulated cell. The VDOS of PIB was investigated applying NS to fully protonated samples (revealing the VDOS corresponding to all the hydrogens in the system) by Frick et al. [39] at 240 K and later by Adams et al. [42] at 25 K, using better resolution. These data are reproduced in Figure 2 in the region below 150 meV.
From the simulations, the Z(E)-function can be calculated from the spectral density of the velocity autocorrelation function as where the velocity autocorrelation function is calculated in terms of the velocity autocorrelation function of each (Hydrogen) atom as: Figure 1. Comparison of the neutron spin echo results obtained from measurements on a fully hydrogenated PIB sample [13] (symbols) and calculated from the simulations (solid lines) at 335 K (a), 365 K (b) and 390 K (c) at the Q-values indicated in (a). Bandpass corrections have been applied to the experimental data (see, e.g., [19]). The value considered for the band pass time was 0.2 ps. For comparison, the F s (Q, t) calculated for all hydrogens at 365 K is shown by the dotted lines in (b).
The Z(E) was calculated running a NVT dynamics of 800 ps at 25 K. However, we note that equilibrating a simulated cell at 25 K is just impossible. Thus the dynamics was run at 25 K but in the cell previously equilibrated at 320 K, i.e., with the density and the structure corresponding to 320 K. Data were recorded every 0.01 ps. The comparison with the results of Ref. [42] is presented in Figure 2. The overall agreement found is rather good, supporting the realism of the simulations also regarding vibrational aspects. We note that there are no quantum mechanic corrections in MD-simulations, which should mainly affect the high-energy region of the spectrum. Therefore, the comparison has to be taken with care in such a region.

Results and Discussion
Once the cell has been properly validated, from the MD-trajectories it is possible to calculate also magnitudes that are not easily (or not at all) experimentally accessible, like e.g., correlation functions in real space, or analyzing separately the different kinds of nuclei (including the 'experimentally invisible' carbons) in the sample. This information can be of utmost help to interpret experimental observations, and/or to check different theoretical frameworks. Such capabilities of the simulations are exploited in this section.

Distinguishing the Different Atomic Species: Chemical Dynamic Heterogeneity in PIB
For a Q-value of 1.5 Å −1 and T = 335 K, Figure 3a shows the intermediate incoherent scattering function of the different kinds of atoms in PIB: main-chain carbons (cC), mainchain hydrogens (cH), methyl-group carbons (mC) and methyl-group hydrogens (mH) (see inset in Figure 3a for the definitions). With tH we denote the total hydrogens (main-chain and methyl-group hydrogens). Figure 3a reveals that: (i) for any type of atom the scattering function decays in two main steps, before and after ≈1 ps; (ii) the curves completely decay to 0 in the dynamic window investigated in this temperature range. The achievement of such a full decorrelation can be attributed to the dynamics of the α-relaxation; (iii) the atomic motions in PIB are very heterogeneous. Defining a global characteristic time τ 0.2 as that at which the function reaches the value of 0.2 (see dotted line in the figure), at 335 K we find a difference of more than one order of magnitude between τ mH 0.2 and τ cC 0.2 . (iv) the decay at short times is more pronounced for methyl-group atoms than for main-chain atoms, and (v) the stretching of the second decay is much stronger for mH than for the other atoms in the system. Finally, as expected, the behavior of tH is very close to that of mH -6 out of the 8 hydrogens in PIB are located in the methyl groups. We have characterized the second decay (times longer than 2 ps) of the incoherent scattering functions of the different atomic species by using the commonly invoked stretched exponential or Kohlrausch-Williams-Watts functional form [1]: Here A is the amplitude accounting for the decay of the correlations in the microscopic regime, 0 < β ≤ 1 is the stretching exponent characterizing the deviations with respect to a single exponential decay and τ w (Q, T) is the characteristic time. In a first step, we fitted Equation (13) to the simulated curves allowing the three parameters to float. For a given species and temperature, we did not find a systematic variation with Q of the values obtained for the β-parameter. Therefore, we calculated the Q-averaged values of this parameter β Q and considered them as representative for each temperature and kind of atom. They are presented in Figure 3b. The highest values of this parameter are found for cC and the lowest for mH (and tH). A tendency to increase with increasing temperature is observed for all the species, in particular for mH. In a second step, we fitted again Equation (13) by fixing the values of the β-parameters to the average values β Q . The descriptions obtained were reasonably good. Given the change of the shape parameter with temperature and considered species, in the following we will discuss the results corresponding to the characteristic times in terms of the average characteristic times τ , that in the case of KWW functions are related with τ w through the expression Figure 4 shows the Q-dependence of τ obtained for the different species. In the low-Q limit (large length scales), all the times roughly coincide. With increasing Q they decrease, showing a certain curvature in the log-log representation.
The values corresponding to cC, mC and cH tend to coincide in the high-Q range (local length scales) where those of mH remain faster.
(see dotted lines). Here, the β-value used at each temperature is the Q-averaged value shown in Figure 3b. This observation implies the Gaussian form of the scattering function (Equation (5) in such a Q-range [4], as can be deduced when Equation (14) is considered together with Equation (13). Note that since the value of the stretching parameter used in the fits is common for all Q-values, considering τ w or τ is equivalent for the discussion of the Q-dependence of the characteristic times. The only difference is a prefactor. The asymptotic Gaussian laws τ ∝ Q −2/β cease to describe the Q-dependence of the characteristic times in a Q-range that varies from one species to the other but for a given kind of atom hardly depends on temperature (see shadowed areas in Figure 5). At higher Q-values, the Q-dependence is weaker than that given by Equation (14), indicating deviations from the Gaussian approximation. Thus, the crossover from Gaussian to non-Gaussian behavior is found not only for the total hydrogens (as it was reported from the NSE experiments [13]), but for all kinds of atoms in PIB.  Table 1).
The anomalous jump diffusion model [14,15] Table 1. For the main-chain atoms, they increase from cC to cH (0.3 to 0.4 Å, respectively); hydrogens undergo larger jumps than carbons in the decaging process. The value obtained for cH is comparable to that reported for the main-chain hydrogens in e.g., polyisoprene (0.42 Å) [14,15]. For methyl-group atoms this length is larger than for the main-chain counterparts and rather similar for both mH and mC (0.50 and 0.54 Å). For the total hydrogens, the value obtained for o is largest. This result can be understood considering that in this case the Q-dependence of the characteristic times reflects deviations from Gaussian behavior which are a consequence not only of the discrete nature of the underlying mechanism for sublinear diffusion (basic ingredient of the anomalous jump diffusion model) but also of the heterogeneous dynamics developed by cH and mH. The more pronounced deviations from Gaussian behavior are captured in this simple model by an effectively larger value of the jump distance. A similar situation was discussed for PVE [19], PVME [20], styrene-butadiene-rubber (SBR) and polystyrene (PS) [70]. Additional ingredients were also included to describe the behavior found in complex systems like e.g., confined liquid crystals [71] or proton diffusion in polymer eletrolyte fuel cell membranes [72]. Furthermore, this intrinsic heterogeneity -together with the occurrence of mH-rotations as we will show in the following section-is also the cause of large stretching for tH. On the other hand, the values of τ o corresponding to cC, cH and mC are rather similar and show nearly identical values of the activation energy E AJD a (see Table 1). For mH the obtained values are smaller and the activation energy weaker. Thus, methyl-group rotations lead to a quicker delocalization of mH atoms. It is noteworthy that in essence the anomalous jump diffusion model incorporates the ingredient of caging by considering a distribution of jumps underlying the diffusive-like motion of atoms in the α process. The cage effect is the main concept involved also in the mode coupling theory (MCT) [73,74]. The applicability of MCT to the case of PIB was explored in a previous work [75]. The localization length r sc involved in the MCT was found to be r sc = 0.30 Å for cC and r sc = 0.50 Å for mH. These values are very similar to those obtained in the anomalous jump diffusion model for the preferred jump distance. We note that this was also the case for cH in PI [15].
From the simulations we thus confirm the experimental observation of stretching of F s (Q, t), its Gaussian behavior at low Qs and deviations at higher Q-values reported from experiments. The larger window covered, the better statistics of the data and the capability of isolating the incoherent contribution has made it possible to determine with more accuracy the value of the shape parameter and the limits of the Gaussian approximation. It is meaningful to compare the outcome of the analysis of the NSE results [13] and those of the simulations. In the interpretation of the NSE data, it is assumed that S NSE (Q, t) = F tH s (Q, t). Figure 1b shows that this is a very good approximation, though the coherent contribution can be noticeable in the Q-range of the first structure factor peak (≈1 Å −1 ). Figure 4b compares the average relaxation times for tH from the simulations with those reported from the experiments. The latter are systematically faster. This is mainly due to the choice of the shape parameter β: the NSE results were fitted by assuming a fixed value of 0.55 deduced from the analysis of collective data. This value is sensitively larger than that obtained from the simulations (see Figure 3b). Such relatively small values of β are due to the heterogeneous motions of cH and mH, i.e., ultimately to the occurrence of methyl-group rotations. These motions are discussed in more detail in the following.

Methyl Group Dynamics
In the high temperature range here investigated, the mH are expected to undergo classical hopping in a rotational potential. We have found indications for the rotations in the low values of the β-parameter used to describe the second decay of F s (Q, t) for mH and in the faster corresponding times observed in the high Q-range with respect to those of the other species. However, since the segmental dynamics (α-relaxation) is superimposed, a two-steps decay of F s (Q, t) after the microscopic regime is not observed for mH in the whole temperature range investigated (see e.g., Figure 3a), preventing a direct characterization of these rotational motions.
Exploring the correlation functions directly in real space (the self-part of the van Hove correlation functions G s (r, t)) can be of utmost help to characterize the atomic motions. For the two extreme temperatures here investigated, Figure 6 represents these functions for mH at different times. At 470 K, the features of the probability distribution function correspond to a diffusive-like process: increasing broadening and shift toward larger distances with increasing time. At 320 K some hints of a jump-like process with the appropriate characteristic length (r HH = 1.78 Å is the distance between the hydrogens in the methyl group) could be envisaged in the range 100...1000 ps. However, the overall behavior is also rather close to a diffusional process.
Thus, in order to characterize the methyl-group rotations the data have to be 'cleaned out' from the diffusive-like motions of the α-relaxation. In a first approximation, it could be assumed that both kinds of motions occur simultaneously and are statistically independent. We have checked this hypothesis for PIB. Under such conditions, the incoherent scattering functions for methyl group hydrogens F mH s (Q, t) would be given by the product of the incoherent scattering functions corresponding to the purely rotational motion and to the α-relaxation process. The latter can be approximated with the incoherent scattering function calculated for main-chain hydrogens F cH s (Q, t). This approach was successfully used in the case of PI [76] and is in fact inspired in the way an experimentalist would face this kind of problem, provided the availability of NS measurements on fully protonated and partially deuterated samples. This allows to obtain in a straightforward way the deconvoluted rotational function by simple division of the two known functions, F Deconv s (Q, t) = F mH s (Q, t)/F cH s (Q, t). If the invoked assumptions are correct, consistent results should be obtained. First of all, its functional form should be that of a localized motion. This implies that its long-time limit should be a constant value, the so-called Elastic Incoherent Incoherent Structure Factor (EISF). The EISF gives account for the geometry of the motion. Figure 7 shows the results of applying the deconvolution approach to PIB. A pronounced first decay in the microscopic regime below 1 ps is seen, that is due to the larger motional amplitude at shorter times of the mH with respect to the cH. A second decay towards a plateau is clearly visible at high Q-values for the low temperature (Figure 7a), that is reached at about 2 ns. At longer times, mainly in the low-Q range, the obtained curves do not really show a constant value, as it should happen if motions are independent. At the highest temperature of 470 K (Figure 7b), no plateau can even be identified in the curves. We can thus infer that the initial assumption of statistical independence fails. If the segmental contribution is approximated by considering the scattering function of main-chain carbons instead of that of main-chain hydrogens, the deconvolution does not work better. The conclusion is that the deconvolution assumption implying the statistical independence of α-process and methyl-group rotations fails for PIB, not only for long times and high temperatures, but also at the lowest temperatures investigated at intermediate times and intermolecular length scales, as we show in detail in the SM.  In order to characterize the rotational dynamics of the methyl groups in PIB, we have calculated from the simulations the correlation functions corresponding to the relative motions of mH with respect to the mC to which they are linked. Figure 8a shows as an example the intermediate scattering function F rel s (Q, t) calculated at 335 K and different Q-values. We applied the rotational rate distribution model (RRDM) [41,77,78] Figure 8b shows the values obtained for the EISF parameter at 335 K (empty diamonds) as compared with the theoretical prediction Equation (16).  As can be seen in Figure 8b, the EISF-values show a more pronounced decay than that corresponding to a 3-fold potential for Q-values > 1.5 Å −1 approx. Obviously, the geometry of the motions carried out by mH is not fully describable by a rigid rotor in a three-fold potential. The mH motions lead to almost a complete decay of EISF at high Q-values. Now let us focus on the parameters describing the distribution function of rotational times. As expected for a localized motion, the characteristic time τ MG o is Q-independent (see Figure 9b). The width of the distribution function σ also displays an approximately constant value for Q < 2 Å −1 but a slight tendency to increase is observed at higher Q-values. We have thus considered the average value of σ in the low-Q region as that representing the actual width of the distribution of rotational times at a given temperature (lines in Figure 9a). The results corresponding to σ are inversely proportional to temperature, implying that the data are compatible with an underlying Gaussian distribution of activation energies E MG a (see SM). The value obtained for σ E is 38.5 meV (447 K). The average activation energy E MG a is obtained from the Arrhenius fit of τ MG leading to E MG a = 0.21 eV (2452 K) and τ ∞ = 0.08 ps. Compared with other glass-forming systems [41], PIB shows a very large value of E MG a (typical values are around 1000 K) but comparable to that displayed by methyl groups located in bisphenol-A moeties (like in polycarbonate, polysulfone and phenoxy [41]).
The direct insight of the atomic motions in real space facilitated from the simulations can help to unveil more details of mH dynamics. Figure 10 shows the radial probability distribution function calculated from the mH trajectories relative to those of mC. At times below ≈ 1 ps the probability shows only one clear maximum broadened by vibrational motions. With increasing time a second maximum is clearly developed which position is separated from that of the first one by the distance r HH = 1.78 Å. The relative populations of the two peaks continuously change with increasing time. At the longest times represented, the radial probability distribution function continuously increases with distance in the physically accessible spatial range. This behavior is qualitatively different from that expected for the 3-fold model of a rigid rotor. This observation, together with the previously commented increase of σ and decrease of EISF detected at high Q-values could be attributed to a superimposed translational component of mH atoms. In order to proof this hypothesis, we generated some random distributions of the accessible positions of such a proton (accordingly to the geometric configuration of a methyl group) allowing for some minor fluctuations for the methyl carbon-proton bond length and also for some stochastic vibrations around the three minima of the methyl rotor. The result is depicted by the filled circles in the insert of Figure 10 that is the one to be expected for a (semi-)rigid rotor of the corresponding geometry with a symmetry axis which direction is fixed in space. Now, if one releases this constraint and allows the axis of the rotor to vary its direction in space, the curve with empty symbols is obtained, which is very similar to those obtained in Figure 10b for the longest time.
To determine the time at which the probability distribution function 4πr 2 G rel s (r, t) deviates from the expected behavior for a rotor with fixed axis, we have considered the value of this function at r = r HH /2 = 0.89 Å. This parameter represents the way how the minimum in the probability distribution function is getting populated. In principle, its value should remain constant in time once the stationary regime for rotations has been achieved. Figure 11b shows the results obtained at 335 K. The time at which deviations start is about 200 ps. For comparison, Figure 11a shows with circles the mean squared displacement r 2 (t) and the non-Gaussian parameter α 2 (t) calculated for main-chain carbons at the same temperature. A characteristic time for the processes developed in the cageing regime t is usually defined as that where the α 2 -parameter reaches its maximum. For main-chain carbons, t is very close to 200 ps at this temperature. Thus, the non-Gaussian events undergone by backbone atoms seem to be at the origin of the peculiarities observed for methyl-group rotations.
The inverted triangles in Figure 11a show the results corresponding to the relative motion of mH with respect to mC at 335 K. At this temperature, for t ≥ 1 ns mH explore the full available space around mC. There, the function F rel s (Q, t) reaches the limit t → ∞ (see Figure 8b). From the comparison between the mean squared displacements in Figure 11a we can see that the development of the localized motions by the methyl groups' hydrogens takes place just during the occurrence of the decaging events in the backbone atoms. At longer times, in the region where r 2 rel (t) is constant, r 2 cC (t) follows the sublinear increase r 2 (t) ∝ t β characteristic for the α-relaxation.
We can now consider the information provided by the distribution functions of rotational times for mH, H(log τ). From these functions, the values of the average residence time for methyl-group hydrogens τ o R can be obtained for the different temperatures (since τ R is related with τ through τ R = 3τ/2, see SM, τ o R = 3τ o /2). This time is represented in the insert of Figure 11b together with the characteristic time for decaging processes of backbone carbons t cC . Both times almost perfectly agree for all temperatures. Thus, methyl-group rotations and main-chain decaging are clearly coupled for PIB in the temperature range investigated. A similar conclusion was deduced from the MD-simulations reported in Ref. [38]. Finally, in the same figure we have also included (dashed line) the characteristic times reported for the so-called δ-relaxation from NMR [36,38] and ESR [37] investigations. The δ-process had been assigned by Slichter [36] to originate undoubtedly from a methyl group rotation. The excellent agreement found between our results and this process definitely supports our analysis.
Furthermore, interestingly, we note that at short times (around 0.08 ps) the r 2 (t) and α 2 (t) functions calculated for the relative motions of mH with respect to mC display marked peaks suggesting a kind of 'recoil' effect (see Figure 11a). This feature occurs just when the second maximum in the 4πr 2 G rel s (r, t) function starts to be populated (see Figure 10a), coinciding with the value of the prefactor τ ∞ in Equation (17). We now address another manifestation of methyl-group dynamics in NS experiments, namely the presence of librational peaks in the inelastic spectra. A peak in the VDOS reveals the energy of the transition between the ground and first excited level E 01 (see SM). We have calculated the VDOS corresponding to mH, cH and tH [see Figure 12a]. In this way, the 'selective deuteration' invoked by Adams et al. [42] was realized. In the region where the NS experiments showed the broad peak for tH we can clearly distinguish two peaks. The first (and mean) peak indeed corresponds to a mH peak (presumably that corresponding to E 01 ). The second peak of the tH-VDOS at about 50 meV, which in the experiments was not really resolved from the first one, is also mainly due to a pronounced peak in the VDOS of mH, but contains additional cH-contributions in the range 40-55 meV. We note that a double-peak structure in the range ≈34-44 meV was also observed in the VDOS of polysulfone, polycarbonate and phenoxy [41]. Thus, there are striking similarities in the behavior of PIB and those systems, regarding both, the presence of a double peak in this energy region and the above commented unusually large value of E MG a (and accordingly, of E 01 ).
The VDOS of the carbons in the system has also been calculated, distinguishing even the four different carbons. The peak at about 40 meV in the VDOS of mH appears as main signature of the VDOS in mC. A peak at about 50 meV is present in the VDOS of all the carbon species, and is specially pronounced in the methylene cC. Thus, the appearance of a second peak in the VDOS of mH at about 50 meV seems to be a consequence of the strong coupling of mH and cC dynamics, as proposed by Adams et al. [42] and is consistent with our conclusions from the analysis of the rotational dynamics of mH at high temperature.
As pointed out in Ref. [42], the peak attributed to methyl-group librations is rather broad. In principle, such broadening could be understood as consequence of the disorder in the glassy material. In fact, this ingredient has already been taken into account analysis of the F rel s (Q, t) functions by invoking the RRDM. From the obtained distribution of activation energies f (E MG a ) and using the proper relationships (see SM) the corresponding distribution function of librational energies F(E 01 ) can be obtained under the assumption of an underlying 3-fold potential. The result is included in Figure 12a. Though the deduced average value E 01 = 36 meV is slightly shifted toward lower energies with respect to the simulated one, the broadening is well reproduced. The discrepancy observed for E 01 could be due to the fact that the analysis of F rel s (Q, t) was made for temperatures above T g .
There, the potential 'seen' by the methyl groups could be slightly different from that in the glassy state, due to the softening of the non-bonded interactions at high temperature.

Revisiting Neutron Scattering Studies
As it was mentioned in the Introduction, to shed light on the origin of the dielectric βprocess, in Ref. [34] backscattering NS-measurements were performed on fully protonated PIB in the temperature range 260-280 K. There, the α-relaxation was not expected to sensitively contribute in the experimental window accessed. Even so, its contribution was considered in the data analysis assuming the merging of α and β-relaxations as statistically independent processes (i.e., the corresponding intermediate scattering functions are combined as a simple product). The β-process was represented by a localized motion involving atomic jumps between two equivalent positions separated by a distance d. The corresponding EISF is then given by The resulting value obtained for d was d = 2.7 Å. If protonic jumps over such large distances actually occurred in PIB, we could have a chance to identify them in the radial probability distribution function. The best conditions would be at the lowest temperature investigated (320 K), where we could expect to have the least contribution from the diffusional motions involved in the α-process. According to the dielectric results, at 320 K the occurrence of jumps should manifest at most at about 1 ns (where the distribution function of the corresponding characteristic times reaches its maximum). As can be seen in the insert of Figure 6a, at such time the function 4πr 2 G tH s (r, t) seems indeed to show a kind of double-peak structure, but the location of the peak (nearly a hump, see arrow in the figure) centered at larger distances would be around, at most, 2 Å. In fact, the inspection of this function for mH and cH separately reveals that the cH atoms do not show any clear hint of jumps. Conversely, we have successfully characterized the localized motions of mH as methyl-group rotations. These are the events producing the second feature at ≈2 Å in the distribution function of tH.
It is worth or remark that the result of an apparently large motional amplitude for localized motions arises as a consequence of the deconvolution procedure assuming statistically independence of localized and diffusive-like processes. To illustrate this, we may define an 'apparent EISF' from the values of the above defined deconvoluted function at t ≈2 ns. At this time, F Deconv s (Q, t) seems to approach a plateau at low tempera-tures, as shown in Figure 7a. To take into account the different vibrational dynamics of mH and cH, we have corrected the results by the value of this function at t = 2 ps, EISF app ≡ F Deconv s (Q, t = 2ns)/F Deconv s (Q, t = 2ps). The resulting values are shown as circles in Figure 8a. This kind of data treatment for mH motions results in apparent EISF-values that are lower than those corresponding to a methyl-group rotation. The fit of Equation (18) to the such obtained data in the experimentally accessed Q-window deliver a jump distance of 2.6 Å (see solid line in Figure 8a). Thus, the coupling between methyl group rotations and main-chain dynamics in PIB invalidates for this polymer the hypothesis of diffusive and localized motions being statistically independent processes.

Conclusions
In this work we have shown the power of properly validated fully atomistic molecular dynamics simulations to unravel the details of the atomic motions at inter-and intramolecular length scales in polymers. First of all, the realism of the simulated cell has been proved by direct comparison of magnitudes calculated from the simulated atomic trajectories and determined on the real polymer by neutron scattering experiments. Taking advantage of the information available from the atomic trajectories, the simulations have been exploited to disentangle the self-motions of the different kinds of nuclei in PIB. The different atomic species display distinct motions, being rotations of methyl groups coupled with main-chain dynamics; in all cases, a crossover from Gaussian to non-Gaussian behavior is observed. In more detail, we have shown that:

•
The motions of the different atomic species are distinct. All atoms manifest Gaussianlike behavior at low Q-values and deviations from it at high Qs.

•
The anomalous jump diffusion model captures this behavior by invoking an underlying distribution of elementary jumps at the origin of the sublinear diffusion of atoms in the α-relaxation regime. The typical lengths of the jumps deduced are ≈0.4 Å for main-chain atoms, in accordance with results reported for other polymers. • Due to rotations, methyl-group hydrogens display an enhanced mobility in the region where the other atoms show decaging processes. These motions are manifested in low values of the shape parameter and shorter characteristic times at small length scales. • The usually invoked statistically independence of rotation and translational motions does not work in polyisobutylene. The methyl-group rotations could be fully characterized by analyzing the relative motion of these atoms with respect to the carbon at which the are linked. The rotational rate distribution model gives account for these motions, though deviations are found at high Q-values. They can be rationalized as due to non-Gaussian processes taking place during the cageing regime of backbone atoms, which allow a superimposed motion of the rotation axis. • The methyl-group rotations are coupled with the main-chain dynamics, as deduced from the coincidence of the characteristic times for rotations and non-Gaussian events within the cage of main-chain carbons. This coupling was also concluded from a study combining MD-simulations and NMR experiments [38] and has also been deduced from the analysis of the coherent structure factor of the present simulations [55]. This shall be the origin of the high value observed for the average potential barrier for methyl-group dynamics, leading consequently to a high value of the librational frequency. In fact, a double peak appears in the VDOS of methyl-group hydrogens, where the peak at higher energies occurs where main-chain carbons also present a characteristic peak. • The information provided by the simulations in real space rules out the occurrence of large amplitude jumps (≈3 Å) associated to the β-process, as proposed from the analysis of neutron scattering data on protonated samples [34]. The deduction of such a motion is an artifact arising when the hypothesis of statistically independence of α-relaxation and localized processes is assumed. The reason is the coupling of methylgroup rotations and main-chain dynamics, possibility that was already suggested by the authors of Ref. [34] as an explanation of the apparently incompatible observations from incoherent and coherent scattering studies.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1. Institutional Review Board Statement: Not applicable.

Author
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are openly available in: https:// ehubox.ehu.eus/s/6zpwz274ntyw5tt.

Acknowledgments:
We thank all our collaborators in the neutron scattering facilities all over the world who, along the past decades, have participated in the experiments. We also would like to thank the rest of the members of our research group in San Sebastián as well as Prof. Dieter Richter and his group in Jülich for the long-standing and fruitful collaboration.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: