Review: Mean-Square Displacements of Simulated Polymers

1. Introduction

In this paper, we consider a few aspects of polymer motion in simulated bead-bond polymer melts. Our interest is the time-dependent mean-square displacement of individual beads, of polymer centers of mass, and of polymer beads with respect to the center of mass of their polymer chain. We advance via a review of representative papers from the literature, applying a quantitative approach based on standard numerical methods. A short letter demonstrating the approach has been published previously [1].

The behaviors of the mean-square displacements are predicted by tube-reptation-scaling models of polymer dynamics [2,3]. For long polymers in a melt, the tube model proposes that the polymer chains surrounding a chain of interest form a transient pseudolattice that restricts the chain of interest’s lateral motions. The transient pseudolattice is described as a tube, an irregular curving barrier having statistical diameter a and contour length L, within which the chain of interest is momentarily confined. According to the tube model, the polymer chain of interest is said to be constrained to move more-or-less parallel to its tube.

For long polymer molecules, the tube model introduces three characteristic times. The hypothesized times are claimed to separate four characteristic time scales for polymer motion. At times shorter than an entanglement time ${\tau }_e$, beads are said not to have encountered the walls of the tube, so bead motions are claimed to be described by the simple Rouse model. At times between ${\tau }_e$ and a later time ${\tau }_R$, Rouse relaxations continue, but are perturbed by the tube walls until they completely relax at a time ${\tau }_R$. ${\tau }_R$, the Rouse time, is the time scale on which the polymer’s longest-lived Rouse mode would relax, if the polymer’s dynamics were described by the Rouse model for a bead-spring chain [4]. ${\tau }_R$ is defined in terms of other chain properties as

$${\tau }_R={\displaystyle \frac{\zeta N\langle R^2\rangle }{3{\pi }^2{k}_{B}T}}.$$

(1)

Here $\zeta$ is a bead drag coefficient, N is the number of beads in the polymer, $\langle R^2\rangle$ is the polymer’s mean-square radius of gyration, $k_B$ is Boltzmann’s constant, and T is the absolute temperature.

Between ${\tau }_R$ and a disengagement time ${\tau }_d$, the polymer diffuses back and forth along the length of its tube. The tube’s ends fray and disappear whenever the polymer backs away from them. Over times greater than ${\tau }_d$, the polymer is no longer confined by the tube and instead performs simple diffusion. The tube model indicates that it does not describe the motions of short chains. Instead, for short chains there are only two time regimes, as described by the Rouse model, the time regimes being separated by the time ${\tau }_R$.

Tube models of polymer dynamics make a series of predictions that are accessible to simulational tests. Of particular interest are predictions of mean-square displacements, including displacements of individual beads and displacements of the chain center of mass. The atomic mean-square displacement for a chain is described by

$$g_1\left(t\right)={\displaystyle \frac{1}{N}}\sum _{i=1}^N\langle |{\mathbf{r}}_i(t+t_0)-{\mathbf{r}}_i\left(t_0\right){|}^2\rangle .$$

(2)

Here $t_0$ and $t+t_0$ are times. The sum may be over the positions ${\mathbf{r}}_i\left(t_0\right)$ of all N beads of a chain, or over the positions of a subgroup of n beads in a longer chain. Motions of a subgroup centered on the chain center are typically considered. The average $\langle \cdots \rangle$ includes averages over all chains and over all values of the initial time $t_0$.

The center-of-mass mean-square displacement is

$$g_3\left(t\right)=\langle |{\mathbf{r}}_{cm}(t+t_0)-{\mathbf{r}}_{cm}\left(t_0\right){|}^2\rangle ,$$

(3)

where ${\mathbf{r}}_{cm}\left(t\right)$ is the position of a chain center-of-mass at time t, the average being over all chains and all initial times $t_0$. Center-of-mass displacements have been reported using the center of mass of all beads in a chain or separately by using the center of mass of the central beads of a chain.

The mean-square displacement of monomer beads relative to their polymer’s center of mass ${\mathbf{r}}_{cm}\left(t\right)$ is described by

$$g_2\left(t\right)={\displaystyle \frac{1}{N}}\sum _{i=1}^N\langle |{\mathbf{r}}_i(t+t_0)-{\mathbf{r}}_{cm}(t+t_0)-({\mathbf{r}}_i\left(t_0\right)-{\mathbf{r}}_{cm}\left(t_0\right)){|}^2\rangle ,$$

(4)

where n is the number of beads in a chain.

When need be, we refer to the mean-square displacement functions collectively as the $g\left(t\right)$.

The Rouse and tube models predict the time dependences of these mean-square displacements, the predictions all having the form $g\left(t\right)\sim t^{\alpha }$. The exponent $\alpha$ depends on the polymer chain length and the time scale over which the time dependence is being described.

For melts of short chains the tube model predicts that polymers have the properties predicted by the Rouse model, namely

$$\begin{array}{ccc}\hfill g_1\left(t\right)& \sim & t^{1/2}\mathrm{for}t\le {\tau }_R,\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill g_1\left(t\right)& \sim & t\mathrm{for}\ge {\tau }_R,\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill g_3\left(t\right)& \sim & t\mathrm{for}t\le {\tau }_R\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill g_3\left(t\right)& \sim & t\mathrm{for}t\ge {\tau }_R.\hfill \end{array}$$

(8)

At times $t<{\tau }_R$, chain motions in a polymer melt are said to be described by the relaxation of pure Rouse modes, leading to $\alpha =1/2$ for $g_1\left(t\right)$ and $\alpha =1$ for $g_3\left(t\right)$. At longer times $t>{\tau }_R$, $\alpha =1$ is predicted both for $g_1\left(t\right)$ and for $g_3\left(t\right)$.

For long chains, the tube model predicts for single-bead mean-square displacements that

$$\begin{array}{ccc}\hfill g_1\left(t\right)& \sim & t^{1/2}\mathrm{for}t\le {\tau }_e,\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill g_1\left(t\right)& \sim & t^{1/4}\mathrm{for}{\tau }_e\le t\le {\tau }_R,\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill g_1\left(t\right)& \sim & t^{1/2}\mathrm{for}{\tau }_R\le t\le {\tau }_d,\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill g_1\left(t\right)& \sim & t\mathrm{for}{\tau }_e\le t.\hfill \end{array}$$

(12)

$g_2\left(t\right)$ and $g_1\left(t\right)$ have the same behavior until $t=t_d$, following which $g_2\left(t\right)$ reaches a plateau at $g_2\left(t\right)=R_g^2$.

For long chains, the tube model predicts for center-of-mass mean-square displacements that

$$\begin{array}{ccc}\hfill g_3\left(t\right)& \sim & t\mathrm{for}t\le {\tau }_e,\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill g_3\left(t\right)& \sim & t^{1/2}\mathrm{for}{\tau }_e\le t\le {\tau }_R\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill g_3\left(t\right)& \sim & t\mathrm{for}{\tau }_R\le t.\hfill \end{array}$$

(15)

The $t^{1/4}$ regime for long chains during ${\tau }_e\le t\le {\tau }_R$ is often viewed as a signature that the tube model is correct. However, as noted by Puetz, et al. [5], in the same time regime the Schweizer mode-coupling model [6,7] predicts a very similar $t^{0.28}$ time dependence for $g_1\left(t\right)$, even though Schweizer’s model gives an entirely different description of bead motion. Similarly, in the following regime, $g_1\left(t\right)$ is predicted by the tube model to be proportional to $t^{1/2}$, while the Schweizer model predicts $t^{0.56}$. From $g_1\left(t\right)$, it may therefore be difficult to distinguish between the tube and Schweizer models.

A few authors calculate the full distribution $P(r,t)$ of displacements r at each time t. When $P(r,t)$ is non-Gaussian, a non-Gaussianity parameter ${\alpha }_2\left(t\right)$ has been introduced. ${\alpha }_2\left(t\right)$ is defined as

$${\alpha }_2\left(t\right)={\displaystyle \frac{3\langle r^4\left(t\right)\rangle -5{\langle r^2\left(t\right)\rangle }^2}{5{\langle r^2\left(t\right)\rangle }^2}}.$$

(16)

If ${\alpha }_2\left(t\right)\ne 0$, interpreting $g_1\left(t\right)$ as supplying a (time-dependent) diffusion coefficient is invalid, because particle motion does not have the Gaussian displacement distribution required by Doob’s theorem [8] if the particle motion is Brownian diffusion. In addition, if ${\alpha }_2\left(t\right)\ne 0$, it would be incorrect to relate $g_1\left(t\right)$ to the dynamic structure factor $S(q,t)$ by using a Gaussian approximation $S(q,t)\sim exp(-q^2{g}_1\left(t\right)/2)$, because this form is only correct if particle displacements have a Gaussian distribution of values at all times [9]. If $P(r,t)$ is not a Gaussian, the Gaussian approximation for $S(q,t)$ may be incorrect even as an order-of-magnitude approximation.

Here we consider simulated determinations of polymer mean-square displacements and what they say about the validity of the deGennes-Doi-Edwards tube model. Our focus is careful numerical study of the time-dependent mean-square displacements, comparing with predictions, covering chains of all lengths, that there are regimes in which the mean-square displacements $g_1\left(t\right)$, $g_2\left(t\right)$, and $g_3\left(t\right)$ have power-law dependences $t^{\alpha }$ on time. Our interest here is less what exponents are found and more whether or not a power-law dependence is present at all.

How might power-law dependences be demonstrated? In the simplest case, if there were a power-law dependence $g\left(t\right)\sim t^{\alpha }$, a graph of $log\left(g\right(t\left)\right)$ against $log\left(t\right)$ would reveal a straight line having slope $\alpha$. Multiple straight-line regions might well be separated by smooth curves marking crossover regimes. Here we replace this crude graphical approach with a systematic quantitative method.

This paper represents an extension of our prior reviews on the phenomenology of the dynamics of polymeric liquids. In our prior volume Phenomenology of Polymer Solution Dynamics [10], we considered the concentration, molecular weight, and time dependences of polymer solution properties, including centrifugation, electrophoresis, polarized and depolarized light scattering spectroscopy, solvent diffusion, segmental relaxation, dielectric relaxation, mutual diffusion, probe diffusion, viscosity, viscoelasticity, non-linear viscoelastic behavior, and outcomes of the same techniques as applied to solutions of spherical colloids. A systematic comparison of experiments on polymer solutions with theoretical models found that hydrodynamic scaling models [11] worked well, but reptation-scaling forms [2], e.g., $D\sim c^{\nu }{M}^{\delta }$ for the self-diffusion coefficient, $\nu$ and $\delta$ being constants, were uniformly inconsistent with experiment. At that time, we were challenged to discuss whether our results applied to melts but were not prepared to do so. This paper is part of a series now responding to that challenge. In a recent paper [12], we extended our reviews to simulations of polymer melts, notably the behavior of the nominal Rouse modes, finding that Rouse mode amplitudes of polymers in melts do not have the properties predicted by Rouse. Here we continue this study, turning to the behavior of the mean-square displacements $g\left(t\right)$ as found in polymer simulations.

The remainder of this paper has three sections. The next Section presents our quantitative method for analyzing $g\left(t\right)$ to determine actual time dependences. We have previously presented tests of the method against representative data. [1]. The following Section applies these quantitative methods to results from the simulation literature. A final Section discusses what has been found. The verisimilitude of our analysis arises from the graphical presentations of our numerical analyses, so the Supplemental Figures present at full scale all of our results.

2. Methods

In this paper, we re-analyze literature reports of $g\left(t\right)$ of simulated polymeric systems. Most measurements in the literature are reported as smooth curves, not as the discrete points normally generated by molecular dynamics simulations. Digitization of reported $g\left(t\right)$ determinations was accomplished using UN-SCAN-IT [13] under manual control. The significant obstacle, at least in some papers, was that the authors conserved space by plotting multiple sets of results on a single graph. In some cases data points were sufficiently overlapped that determinations of $g\left(t\right)$ became impossible.

We advance by curve fitting. Identifying $y\left(t\right)=log\left(g\right(t\left)\right)$ and $x=log\left(t\right)$, we fit $log\left(g\right(t\left)\right)$ to polynomials

$$y\left(t\right)=\sum _{n=0}^N{k}_n{x}^n.$$

(17)

Here the $k_n$ are the fitting parameters, to be obtained from linear-least-mean-squares analysis, and N is the truncation limit, the largest value of n in a given fit. What value of N is appropriate? We repeatedly fit measurements to equation 17 with progressively increasing values of N. The mean-square difference between the fitted curve and the measurements was obtained. We consistently found that there was an upper limit on N above which further increases in N gave at most a very limited reduction in the mean-square difference between the simulational measurements and eq. 17. $N=8$ was generally used as the cutoff in the studies here. As seen below and in teh Supplemental Material, our fitted polynomials accurately follow simulational determinations of the $g\left(t\right)$, with no indication of a need to fit different time regimes with distinct polynomials employing different sets of $k_n$.

In our analysis,

$$K_1={\displaystyle \frac{dy\left(t\right)}{dx}}.$$

(18)

is the logarithmic derivative $dlog\left(g\right(t\left)\right)/dlog\left(t\right)$ of our polynomial. For each t, $K_1$ gives the instantaneous value of $\alpha$ in a local-in-time scaling law $g\left(t\right)\sim t^{\alpha }$. If power-law behavior were present in some time regime, within that regime $K_1$ would be constant up to noise. To test for constancy, we also calculated the second logarithmic derivative

$$K_2={\displaystyle \frac{d^{2}y\left(t\right)}{d{x}^2}}.$$

(19)

In general, $K_1$ and $K_2$ both depend on t. In a power-law regime, one expects $K_2\approx 0$ over an extended range of times. Noise in the simulation, the digitization process, and the numerical fitting steps will lead to noise in $K_2$, so in a power-law regime $K_2$ will in general not be exactly zero. $K_2\approx 0$ will also be found over a narrow range of times wherever $K_1$ has a maximum, a minimum, or a saddle point, but these cases are readily distinguished from power-law behavior.

Our analytical approach appears to be more sophisticated than has oft-times been applied to identify power-law regimes. This observation is not meant as a criticism of earlier work. The earlier processes were entirely adequate for the needs of their users. We only claim that our more powerful mathematical technique reveals hitherto obscure features that were actually present in the original studies. In particular, we do not attempt to force the interpretation that $g\left(t\right)\sim t^{\alpha }$, i.e., that $log\left(g\right(t\left)\right)$ actually has a linear dependence on $log\left(t\right)$ over some range of times. We instead determine the dependence of $log\left(g\right(t\left)\right)$ on $log\left(t\right)$, and inquire as to whether that dependence is linear in $log\left(t\right)$, as would be the case if $g\left(t\right)\sim t^{\alpha }$.

Our procedure creates smooth curves that, at least in the studies discussed here, always visibly pass within noise through the measured $g\left(t\right)$. The logarithmic derivative $K_1$ of the fitted function provides, for each time, the local value of the exponent, as could in principle be obtained by a fit of a power law to $log\left(g\right(t\left)\right)$ over a narrow range of $log\left(t\right)$. As a practical matter, we found that direct fits of $t^{\alpha }$ to $g\left(t\right)$ over short ranges in time were quite noisy, so we tried but did not pursue examining $K_1$ as obtained from a multitude of local fits. In contrast, logarithmic derivatives of equation 17 are nearly noiseless. This distinction between our approach and local fits is the well-known issue that, while numerical integration tends to suppress noise, numerical derivatives tend to amplify noise. We emphasize that polynomial fits are interpolants, not extrapolants. They fill in the space between measured points, but are unreliable as predictors outside of the range of times over which the fit was made.

3. Analysis

In this Section, we use our numerical analysis approach to consider $g\left(t\right)$ as reported by a series of authors. Strong conclusions, as will be presented here, need strong support. We therefore include eighty determinations of $g\left(t\right)$. While we could have selected from these an exemplary set of a few systems, this path would have led to suggestions that we had cherry-picked data to support our conclusions. On the other hand, there is a certain uniformity in the time dependences of the $g\left(t\right)$ and of $K_1$. In almost every case: $g\left(t\right)$ increases monotonically, on a log-log plot more rapidly at early and late times and less rapidly at intermediate times. With increasing time, the first logarithmic derivative $K_1$ decreases more or less smoothly to a minimum, corresponding to an inflection point of $g\left(t\right)$, and then increases more or less smoothly again. In a very few cases, we find that $K_1$ has a single extended region in which it is constant, corresponding to a power-law dependence of $g\left(t\right)$ on t. In no system do we observe more then one power-law regime for $g\left(t\right)$, other than at the short- or long-time limits.

Behbahani and Schmid [14] report simulations of a Grest-Kremer bead-spring polymer model [15]. Beads had mass m and a purely repulsive Weeks-Chandler-Anderson potential [16] with length scale $\sigma$ and energy $ϵ$. Bead-bead bonds were represented with a FENE (Finite Extensible Nonlinear Elastic) potential [17]. Studies were made for polymers having between 5 and 1000 beads with a minimum of 192 chains in a simulation box. Time was measured in natural units $\tau ={(m{\sigma }^2/ϵ)}^{1/2}$, simulations being extended out to ${10}^5\tau$ to ${10}^8\tau$ depending on chain lengths. Behbahani and Schmid note a series of estimates of the entanglement length, one group of estimates in a cluster near 50 beads and another in a cluster of 80 or a few more beads. The authors use the former number to estimate that their thousand-bead polymer has on average some 20 entanglements. The estimated distance between entanglements, which was identified as the tube diameter a, was estimated to be approximately $9\sigma$. From their simulations, Behbahani and Schmid calculate the mean-square bead displacements, the mean-square displacements of the central beads, the mean-square center-of-mass displacements, the time dependence of the polymer end-to-end vector, the single-chain dynamic structure factor, the stress relaxation modulus, and the zero-shear viscosity.

We first consider single-bead motions. Figure 1 and Figure 2 show Behbahani and Schmid’s determinations of mean-square single bead displacements $g_1\left(t\right)$ in their systems, and our fits of $log\left(g_1\left(t\right)\right)$ to a power series in $log\left(t\right)$. We used fits to eighth-order polynomials, except thatin one system (the 400-bead polymer, Figure 2c) the eighth-order fit to $g_1\left(t\right)$ was not entirely satisfactory, so a tenth-order fit was used instead. What features are revealed by $K_1$? For the 10-bead chains, $K_1$ has a minimum, and correspondingly $g_1\left(t\right)$ has an inflection point near $t=22$. For longer chains, there appear to be two features that emerge and separate with increasing chain length. The 1000-bead chains show these clearly. First, $K_1$ has a minimum, so $g_1\left(t\right)$ has an inflection point. Second, at shorter times before reaching its minimum, $K_1$ has an extended region in which it decreases linearly with $log\left(t\right)$ (and, correspondingly, $K_2\approx 0$). For shorter chains, the near-linear region acquires weak distortions from linearity. Finally, for the 100- and 50-bead polymers, the linear region and the minimum in $K_1$ blend into each other, so that one could reasonably describe the region near the apparent inflection point as showing power-law behavior. The value of $K_1$ at its minimum, the slope of $log\left(g_1\left(t\right)\right)$ at the inflection point, decreases with increasing chain length, from 0.65 for the 10-bead chains to 0.31 for the 1000-bead chains. Over the same series of simulations, the location of the minimum increases, from $t/\tau =22$ for the 10-bead chains to $t/\tau =8.1\times {10}^4$ for the 1000-bead chains. Table 1 shows the full dependence of $K_1$ on chain length.

Close to the shortest and longest times studied, we sometimes see an additional feature in the time dependence of $K_1$. As seen in Figure 1a,d and Figure 2a,c, the feature is a local maximum in $K_1$, in these figures located just before the largest times studied. We note the feature, but are not certain whether it is real or is a consequence of the operation of our polynomial fitting process near end points.

Behbahani and Schmid also report $g_1\left(t\right)$ as calculated using only the displacements of the central beads of the polymer. It is generally stated that polymer end beads are more mobile than beads near the center, for reasons not arising from core parts of the reptation-tube model, so a study of central beads may well be a better test of that model. Figure 3 shows their results for the 1000-bead polymer, and our fit to their results. $log\left(g_1\left(t\right)\right)$ for central beads indeed has an inflection point (a minimum in $K_1$) for $K_1\approx 0.25$, this being the slope predicted by the reptation-tube theory as a central outcome of the model. Those authors justifiably propose that they have observe the predicted $t^{1/4}$ region emerging from surrounding transition regions. One notes, however, that the minimum of $K_1$ progressively decreases with increasing polymer chain length, so a needed experiment would study somewhat longer polymers, say 1300 or 1500 beads, to determine whether for longer chains the minimum in $K_1$ remains near $K_1\approx 0.25$, or whether with increasing chain length the minimum in $K_1$ continues to decrease, in which case the observed 0.25 would be a coincidence arising from a fortunate choice of chain length. As with $g_1\left(t\right)$, at times shorter than that at which the inflection point is seen, there is an extended region in which $log\left(g_1\left(t\right)\right)$ decreases linearly with $log\left(t\right)$. The physical interpretation of this observation is unclear.

We finally consider Behbahani and Schmid’s determinations of center-of-mass motions. Figure 4 and Figure 5 show their measurements of mean-square center-of-mass displacements $g_3\left(t\right)$. Once again, the time dependence of $K_1$ shows that for each polymer length $g_3\left(t\right)$ has an inflection point. The slope $K_1$ at the inflection point decreases with increasing polymer length, from $K_1=0.87$ for the 10-bead polymers to $K_1=0.50$ for the 1000-bead polymers. With increasing chain length, the time at which the inflection point is found also increases, from $t\approx 12$ for the 10-bead polymer to $t\approx 8.8\times {10}^4$ for the 1000-bead polymer. At times before the inflection point, $K_1$ has distinct short- and long-chain behaviors. For the shorter chains, $\le 150$ beads per chain, $K_1$ decreases smoothly with decreasing slope until it reaches the inflection point. For the longer chains, 200 beads or more per chain, $K_1$ acquires a shoulder covering 2-3 decades in time before it reaches its minimum.

Chang and Yethiraj [18] explore a polymer model in which polymer molecules may or may not be allowed to pass through each other, depending on the value of a single parameter. They model polymers as chains of hard spheres having diameter $\sigma$ coupled by bonds whose lengths can fluctuate freely over a range of lengths $(1\pm {\delta }_B)\sigma$. Beads interact via elastic collisions. In different simulations, chains had lengths N between 8 and 512 beads at bead volume fractions of 0.3 or 0.4, with ${\delta }_B$ ranging from 0.05 to 0.9. Chain crossings are blocked for ${\delta }_B\le 0.4$, marginally possible for ${\delta }_B=0.45$, and become more frequent as ${\delta }_B$ is further increased. They calculated a diffusion coefficient D for the chains, finding that $DN$ was approximately constant for $N<64$ but decreased at larger N, and inferred that $N\approx 64$ marked the onset of entanglement. Static chain properties, the pair correlation function and the static structure factor, are nearly independent of ${\delta }_B$. Chain dynamic properties are sensitive to the chain crossing parameter, showing one behavior for ${\delta }_B\ge 0.5$ and a different behavior for ${\delta }_B\le 0.45$. In particular, for long chains the chain self-diffusion coefficient scales as $1/N$ for large ${\delta }_B$ and as $1/N^2$ for small ${\delta }_B$.

Of particular interest here, Chang and Yethiraj calculated the mean-square displacements of the chain centers of mass and the chain central monomers for ${\delta }_B=0.4$, $0.45$, $0.5$ and $0.6$. The 0.4 and 0.45 results are nearly indistinguishable. These are the chain-nearly-non-crossing simulations, so we will consider only the ${\delta }_B=0.45$ simulation.

As seen in Figure 6, Chang and Yethiraj’s determinations of $g\left(t\right)$ are all described well by polynomial curves. The first derivatives are equally well-behaved. Beginning with the center-of-mass displacement measurements as described by $g_3\left(t\right)$: As seen in Figure 6a, for ${\delta }_B=0.45$, $K_1$ increases smoothly from approximately 0.7 to 0.8 with increasing t, and at the largest t appears to increase to 0.9 or so. Figure 6b reveals for ${\delta }_B=0.5$ that earlier times $K_1\approx 0.74$ is nearly constant, increasing by 0.03 over its near-constant range; at $t/\tau >4000$, $K_1$ increases progressively toward $K_1=1$. From Figure 6c, for ${\delta }_B=0.6$, $K_1$ at first increases rapidly, and then rises gradually from 0.8 to a plateau near 0.95.

Chang and Yethiraj also calculated $g_1\left(t\right)$ for the central monomers of their chains. For ${\delta }_B=0.45$, the corresponding $K_1$ lies in the range 0.35-0.40, except at the longest and shortest times reported, as seen in Figure 6d. Figure 6e gives results for ${\delta }_B=0.5$, for which $K_1$ is consistently near 0.39. except at the two extrema, where it rises to 0.45 or 0.55. On increasing ${\delta }_B$ to 0.6, $K_1$ becomes larger, and increases smoothly from 0.45 to 0.5 over most of the observed time range, with a modest further increase to about 0.6 at the largest times studied.

Chang and Yethiraj propose that the power-law time regimes of the reptation-scaling model are seen in their central-bead results for ${\delta }_B=0.4$, these results being shown in Figure 6d. In this figure we find that $g_3\left(t\right)$ follows a smooth curve with a not-quite-constant slope in the range 0.35-0.4, the slope never decreasing to the 0.25 of the scaling model.

Tsalikis, et al. [19] report simulations of ring and linear polyethylene molecules in the melt. They examined monodisperse systems of chains having 120, 227, or 455 monomers, corresponding to molecular weights of 5330, 10,044, or 20,086 Da, which they refer to as the PEO-5K, PEO-10K, and PEO-20K melts. For their systems, they estimate $N_e\approx 45$. Simulations were carried out using Fischer’s united-atom force field [20,21]. For the PEO-20K melt, the simulation cell contained more than 50,000 united atoms. Tsalikis, et al., report detailed studies of static and dynamic properties of their polymers, including $g_1\left(t\right)$. Figure 7 shows $g_1\left(t\right)$ from their Figure 13, our fits to their simulation data, and the corresponding first derivatives. The ring data correspond to mean-square displacements of all monomers, while the linear chain data are for the ’most central’ (their words) monomers. They terminated their measurements at the Rouse times for each system. At each molecular weight, the linear chain’s central monomers move more rapidly than the monomers of the corresponding ring. Tsalikis, et al., interpret their results on the PEO-5K and PEO-10K systems as showing as showing two distinct breaks in $g_1\left(t\right)$, the lower being near 2nS and the higher being at 150 or 490 nS.

Figure 7 shows Tsalikis, et al.,’s measurements of $g_1\left(t\right)$ in their six systems. Points are their measurements as digitized here; the original measurements were reported as solid lines. In Figure 7a, the 120-bead linear molecules, the slope $K_1$ falls from 0.6 at 0.2 nS to 0.38 at 22 nS and then rises again to 0.8 at the longest time reported. $K_1=0.38$ is visibly a local minimum. Figure 7b shows the corresponding results for the 120-bead ring polymer. As a function of time $K_1$ has a minimum at 20 nS, very nearly the same time as for the linear polymer, but here $K_1$ is 0.44. At longer times, $K_1$ rises to 0.95, which does not appear to be a maximum. Figure 7c,d show the corresponding results for the 227-bead linear and ring polymers. In these two systems, $K_1$ has minima at 0.30 and 0.43, respectively. For the ring polymer, a case can be made that there is a region within which $K_1$ is nearly constant, namely $K_1$ is within 1% of its 0.43 minimum for $10\le t\le 80$ nS. Figure 7e,f, for the 455-bead linear and ring polymers, somewhat repeat the behavior seen for the 227-bead polymers. For the linear 455-bead polymer, $K_1$ falls from 0.57 at 0.33 nS to a minimum of 0.2 near 450 nS and then rises steeply at later time delays. The 455-bead ring polymer, like its 227-bead ring counterpart, has an extended minimum between 70 and 250 nS with $K_1\approx 0.39$.

For each of the larger rings, we have therefore found a single power-law region, though not with the same exponent in the two cases. With increasing chain length, the time dependence of $g_1\left(t\right)$ softens, but key characteristic features only change numerically. In describing their results, Tsalikis, et al., instead say that there are clear breaks in $g_1\left(t\right)$ separating extended regions in which $g_1\left(t\right)$ follows one power law or another. We do not see these multiple power-law regions in our analysis.

Takahashi, et al. [22], report extended simulations of a united-atom model for polyethylene. They used the TraPPE-UA model, including Lennard-Jones potentials between non-bonded atoms, rigid links using the LINCS algorithm, bond bending with a harmonic potential, and torsional potentials along the chain backbones. Chain molecular weights ranged from 423 to 2807 Da, with 150 to 1000 chains in a periodic box. They estimated $N_e\approx 80$ or equivalently $M_e\approx 1100$, implying that their chains were unentangled or slightly entangled, with at most 2-3 entanglements on a single chain. The mean-square displacement $g_1\left(t\right)$ of the chain central bead was calculated as a function of time. From the time correlation function of the chain end-to-end vector, a nominal relaxation time ${\tau }_R$ was inferred. A transition in the molecular weight dependence of ${\tau }_R$ was observed near a molecular weight of 1000, with ${\tau }_R\sim M^{2.1}$ at smaller molecular weights and ${\tau }_R\sim M^{2.7}$ at larger molecular weights.

Figure 8 shows our digitization of their measurements of $g_1\left(t\right)$, our polynomial fits to those measurements, and the computed logarithmic first derivatives $K_1$. Agreement between the measured $g_1\left(t\right)$ values and the polynomial fits is uniformly excellent. Of particular interest is Figure 8b. At longer times $>{10}^3$ pS, $log\left(g_1\left(t\right)\right)$ lies on a straight line having slope unity. Correspondingly, the plot of $K_1$ against time shows a straight horizontal line, confirming that $g_1\left(t\right)$ shows power-law behavior with $g_1\left(t\right)\sim t^1$. This result confirms that our method does reveal power-law behavior when it occurs. At shorter times, the slope of $g_1\left(t\right)$ in each of these systems has a local minimum, corresponding to $g_1\left(t\right)$ having a shorter-time inflection point. The slope at the inflection point decreases with increasing polymer molecular weight, the slopes at the six minima in order of increasing polymer length being 0.73, 0.64, 0.60, 0.53, 0.46, and 0.44, respectively.

For the 2807 Da polymer, Takahashi, et al., note that the sequence of straight-line slopes from the reptation-scaling model approximately resembles the curving $g\left(t\right)$, except that the “…threshold values …”, the moments in time at which one transitions from one power law to the next, are “…unclear …” in the $g_1\left(t\right)$ graph. We agree with their conclusion, namely that except at the longest time $g_1\left(t\right)$ follows a smooth curve with no ’transitions’.

Peng, et al. [23] report simulations of dynamically asymmetric A-B blends as functions of chain lengths and concentrations. Under some conditions the A and B chains become immiscible. The authors studied both phase-separated systems and systems which do not phase separate over the duration of the simulations at the simulated temperature. Their polymer beads interact via a Lennard-Jones potential and, along each chain, a FENE potential. A and B chains all use the same Lennard-Jones and FENE parameters. The A chains were made stiffer by endowing them with a harmonic bond-bending potential and a torsional potential. B chains had a fixed length $N_B=50$; A chains had lengths in the range 3-100. Nominal entanglement lengths were given as 28-45 for the A chains and 85 for the B chains, so the B chains were not entangled and the A chains would at most have been weakly entangled. The simulation cells all contained approximately 40,000 beads, all having the same mass. Simulations calculated an extended series of dynamic properties including bead and center-of-mass mean square displacements, Rouse mode temporal autocorrelation functions, and deviations from a Gaussian distribution of displacements, as represented by the non-Gaussianity parameter ${\alpha }_2\left(t\right)$ of equation 16.

Figure 9 shows a representative sample of Peng, et al.’s determinations of the mean-square displacements of B beads in a pure-B melt and in A-B mixtures, both in phase-separated and phase-non-separated systems. The A beads, when present, were at a mole fraction of 0.1. The dependence of $g_1\left(t\right)$ on the chain length of the A beads was very weak.

These are rather short polymers. In each system, at the shortest times $K_1$ approached 2.0. $K_1$ then declined rapidly to an inflection point, the minimum of $K_1$ being in the range 0.15-0.28 at times 2-2.5. $K_1$ then increased again to 0.55-0.59 in the four displayed systems. From the measurements, $g_1\left(t\right)$ in these systems show: shorter-time asymptotic behavior, an inflection point at intermediate times, and near-power-law behavior at the longest times studied. There are no indications of additional power-law-regimes at intermediate times. Our finding of a long-time power-law regime is in quantitative agreement with the analysis of Peng, et al., [23] who report power laws with $\alpha \in 0.56-0.58$ at times $t>{10}^3$.

Peng, et al., also report that the distribution function for displacements of the B beads is substantially non-Gaussian at intermediate times $1\le t\le 1000$. In the pure-B melt, ${\alpha }_2\left(t\right)$ at maximum was $>0.4$, while in the mixture ${\alpha }_2\left(t\right)$ at maximum was $>0.9$, the maxima in the two systems being seen for $10<t<100$. ${\alpha }_2\left(t\right)$ did not go to zero out to the largest observed times. It follows from Doob’s theorem that the forces on the individual beads do not have uncorrelated Gaussian random distributions [8].

Hsu and Kremer [24] examined static and dynamic properties of long-chain linear polymers, lengths 500, 1000, or 2000 beads, in simulated melts. Simulations each had 1000 chains in a periodic box, so simulated systems contained between 0.5 and 2 million beads, depending on the chain length. Polymer beads interacted with a truncated Lennard-Jones potential, a FENE potential representing chemical bonds, and a bond-bond bending potential. The authors determined the mean-square displacements $g_1\left(t\right)$ of individual beads in the center half of each chain, the mean-square displacements $g_2\left(t\right)$ of beads relative to the centers of mass of their respective chains, and the mean-square displacements $g_3\left(t\right)$ of the centers of mass of the chains. Hsu and Kremer infer two estimates of $N_e$, viz, 26 or 32 beads, from which one finds that the chains are nominally strongly entangled.

Figure 10 shows their measurements of the mean-square displacements for the 500 and 2000 bead chains, and our analysis of their measurements. Several features are noteworthy. First, the polynomial fits to the measurements visibly and accurate represent the reported mean-square displacements over a full range of reported times and displacements. Second, considering the first derivatives $K_1$ of the fitted curves, at early times the $g\left(t\right)$ all increase rapidly, so $K_1$ is close to 2.0. $K_1$ decreases rapidly until it reaches times in the range 1-10, following which it decreases less rapidly until it reaches a single minimum. Correspondingly, each $g\left(t\right)$ while increasing is concave downward until it reaches an inflection point.

An interesting exception is seen in Figure 10d, showing the mean-square displacements of central single beads of the 2000-bead chain. At long times, $K_1$ is in the near-constant range 0.26-0.24, while the second derivative of the fitted curve becomes $K_2=0.0\pm 0.01$. That is, there is a region at times $2.7·{10}^4-1·{10}^6$ in which $K_1\approx 0.25$ is very nearly constant. Correspondingly, $g\left(t\right)$ in this region follows a power-law with $g_1\left(t\right)\sim t^{0.25\pm 0.01}$.

The slope of $g\left(t\right)$ at an inflection point does not correspond to a power-law exponent, but it might nonetheless be of interest as reflecting a ghostly shadow of a power law almost but not quite buried under two adjoining non-power-law regimes. The slopes at the inflection points and the single power-law regime appear as Table 2.

The exception in Figure 10d, the revealed power-law behavior, has general significance for the analysis in this paper. Polynomials describe smooth curves, not (in general) straight lines. One might therefore be concerned that, because we use polynomial fits, our fitting function might somehow mask the presence of legitimate power-law behavior. As seen in this figure, a power-law dependence of $g\left(t\right)$ on t is captured by our fitting process, when the dependence is present.

Hsu and Kremer interpret their mean-square measurements by assuming that the reptation-scaling model predictions of power laws $g\left(t\right)\sim t^{\alpha }$ and the model-predicted exponents are correct. They were able to identify times near which power-law curves are tangents to the $g\left(t\right)$ results. Good agreement between the reptation-scaling power laws and the simulated curves was therefore found.

The results of Brodeck, et al., [25] are of general interest for this paper, because Brodeck, et al.,’s $g\left(t\right)$ do have power-law regimes that are revealed by our approach, thus confirming that our approach will find power-law behavior when it is present. Brodeck, et al., report simulations of polyethylene oxide/polymethylmethacrylate (PEO/PMMA) and simple bead-spring blends. PEO/PMMA is notable for its high degree of dynamical asymmetry, the glass temperatures of the two components being separated by 200 K. The PEO/PMMA simulations were based on the COMPASS force-field; the corresponding simulation cell contained equal numbers of PMMA and PEO chains, each chain containing 21 monomers.

Figure 11 and Figure 12 show mean-square displacements of single hydrogen atoms and of chain centers of mass, respectively, for the PEO chains in the PEO-PMMA mixtures. Simulations were made at each of four temperatures and times 0.01-${10}^5$ in natural units.

Figure 11 shows $g_1\left(t\right)$ for PEO in a PMMA melt at four temperatures. In each case, $g_1\left(t\right)$ clearly has a power-law region, as confirmed by our quantitative analysis. For the four Figure 11a-d, that being in order of decreasing temperature, $K_1$ takes on values 0.48, 0.42, 0.36, and 0.31, respectively. For each $g_1\left(t\right)$, there is an extended region from $t\approx {10}^2$ to t in the range 3000-8000 in which the second derivative satisfies $\mid K_2\mid \le 0.01$. Between the shortest and longest times at which the power-law behavior is seen, $K_1$ either is constant or perhaps increases by 0.02 over the full time range.

In each Figure, the ability of our method to capture minor features of $g_1\left(t\right)$, features that would not necessarily be apparent to the eye, is evident. As remarked by Brodeck, et al., in all four studies of $g_1\left(t\right)$, at the longest times $g_1\left(t\right)$ deviates from a power law. These deviations are captured by our fitting procedure. In Figure 11b, at long times $g\left(t\right)$ increases rapidly, while $K_1$ increases to unity. In Figure 11c.d, $g_1\left(t\right)$ appears to approach a plateau, so that $K_1\to 0$. Less obvious is the behavior of $g_1\left(t\right)$ at early times, times before the power-law regime is reached. In each case, at early times, $g_1\left(t\right)$ increases rapidly, with $K_1$ as large as $1.2$. With increasing t, $K_1$ then rolls over until it reaches the slope associated with the power-law regime. However, in each Figure, before $K_1$ reaches its power-law value, it has a hitherto-unmentioned local minimum, so that $g_1\left(t\right)$ has an inflection point and $K_1$ actually approaches its power-law behavior from values below its value in the power-law regime.

Brodeck, et al., report that the measured distribution functions for the displacements of the hydrogen atoms are not Gaussian in form. The non-Gaussianity parameter ${\alpha }_2\left(t\right)$ is small at early times but increases with increasing time.

Figure 12 shows $g_3\left(t\right)$ for PEO molecules in PMMA melts, these being the same systems as those described in Figure 11. Of these, as seen in Figure 12a, at 500K there is a somewhat narrow power-law region in which $K_1\approx 0.82$. Figure 12b–d show a fundamentally different behavior for $g_3\left(t\right)$, namely $g_3\left(t\right)$ has an initial rapid increase, then a near-plateau and corresponding local minimum in $K_1\left(t\right)$, and finally a continuously increasing slope until a long-time behavior is reached. The local minimum in $K_1$ becomes deeper with decreasing T, namely from $K_1=0.53$ at 500 K to $K_1=0.11$ at 300 K. The long-time limiting behaviors in Figure 12b–d are unique, namely in the three figures $g_3\left(t\right)$ at long times shows a progressive increase, a plateau, and an inflection point, respectively.

Brodeck, et al., [25] continued their study by simulating a two-component bead-spring model of a polymer melt. The polymer A and B components differ in their interaction diameters, the slower A component having larger diameter beads. Studies were made at each of six temperatures between 0.4 and 1.5 in natural units. Their results were qualitatively the same as those they had already found in their PEO:PMMA simulations, confirming that the behaviors seen above are not peculiarities unique to PEO-PMMA blend melts.

Stephanou, et al., [26] report an extensive effort to apply Everaers, et al.,’s Primitive Path Analysis [27] approach to a discussion of the dynamics of polymer bead-spring models. Their approach was to reanalyze previously published simulations [28,29] of linear polyethylenes and linear polybutadienes. The primary effort was a consideration of the average rate at which a polymer molecule diffusively escapes from the hypothesized tube of the deGennes-Doi-Edwards [2] model of polymer dynamics. However, as part of the study, Stephanou, et al., report the atomistic mean-square displacement $g_1\left(t\right)$ of the center section of polyethylenes in a melt. They also report for the same chains the mean displacement $N\left(t\right)$ of primitive path segments perpendicular to a nominal primitive path.

Stephanou, et al.,’s measurements of $g_1\left(t\right)$ and $N\left(t\right)$ over times 0.4-900 appear in Figure 13 together with our polynomial fits to their results. Stephanou, et al., assert for $g_1\left(t\right)$“Four distinct regimes are seen in the figure exactly as proposed by the reptation theory.". We would describe $g_1\left(t\right)$ and $N\left(t\right)$ differently, namely to our eyes $g_1\left(t\right)$ and $N\left(t\right)$ both follow smooth, monotonically increasing curves. Our analysis thus does not find the four power-law regimes that Stephanou, et al., see in their data. Namely, for $g_1\left(t\right)$, we find that $K_1$ has one pronounced local minimum, and correspondingly an inflection point in $g_1\left(t\right)$, and one time domain of significant extent within which $K_1$ is constant, namely $K_1\approx 0.5$ for times $t\approx 1$, corresponding to a single $t^{0.5}$ power-law regime at early times. Stephanou, et al.,’s determinations of $N\left(t\right)$ also find a smooth curve, whose slope $K_1$ increases monotonically with time. We do not find power-law behavior in $N\left(t\right)$.

Figure 13. (a) Mean-square displacement $g_1\left(t\right)$ of the central beads of a chain, and (b) mean displacement of individual primitive path segments perpendicular to the chain’s primitive path, based on Stephanou, et al. [26]. Heavy lines are digitizations of Stephanou, et al.’s data, circles are from the two eighth-order least-mean-square polynomial fits, and the thin solid line represents the first derivative $K_1$ of the polynomial fit.

Figure 13. (a) Mean-square displacement $g_1\left(t\right)$ of the central beads of a chain, and (b) mean displacement of individual primitive path segments perpendicular to the chain’s primitive path, based on Stephanou, et al. [26]. Heavy lines are digitizations of Stephanou, et al.’s data, circles are from the two eighth-order least-mean-square polynomial fits, and the thin solid line represents the first derivative $K_1$ of the polynomial fit.

Likhtman, et al., [30] report a study of bead-spring polymers as described by the Kremer-Grest model. Beads have Lennard-Jones potential energies; springs are described by the FENE potential; springs are tuned to prevent chain crossing. Their primary interest was calculating the linear viscoelastic spectra of linear polymer melts, but they also calculated mean-square polymer displacements. Their measurements cover more than six orders of magnitude in time. We consider here their results for the mean-square displacement of chain central monomers of chains containing 50, 100, 200, or 350 beads. Likhtman, et al., estimate that even their longest chains are only modestly entangled, with no more than seven entanglements per chain. They explored the temporal form of $g_1\left(t\right)$ by dividing their results for $g_1\left(t\right)$ by a time dependence $t^{1/2}$ identified with the Rouse model and separately with a time dependence $t^{1/4}$ associated with the reptation model, and examined graphs of $g_1\left(t\right)/t^{1/2}$ and $g_1\left(t\right)/t^{1/4}$. They found that plots of these functions were “…very rich and full of features.” Here we confine ourselves to $g_1\left(t\right)$ of the central monomer, of which they say “…we find it very difficult to identify any clear power-laws in this plot …”, a conclusion with which we agree.

Figure 14 shows their determinations of $g_1\left(t\right)$, our polynomial fits, and its logarithmic derivative $K_1$ for each polymer size. Qualitatively, the four figures are very similar. In each case, the polynomial fit (dark line) accurately describes the measurements. The first derivative curves are all qualitatively the same. $K_1$ slowly decreases to a local minimum, and then at longer times increases more rapidly. The minimum value of $K_1$ decreases with increasing polymer length, namely to 0.47, 0.40, 0.34, and 0.30, respectively. These are the only minima, the only regions where the slope of $K_1$ changes sign from negative to positive. Correspondingly, each $g_1\left(t\right)$ curve has a single inflection point. The time at which the inflection point is found increases from $\approx 470$ to $\approx 37,000$ with increasing chain length. There are other regions, notably at times $10\le t\le 100$, where $K_1$ is relatively slow-changing, but in all those regions the time derivative $K_2$ of $K_1$ keeps a single sign at all points.

Zhou and Larson [31] calculate $g_2\left(t\right)$, the mean-square monomer motion relative to the chain center of mass, for a bead-spring polymer model incorporating a Lennard-Jones potential, a FENE potential for each bead-bead bond, and a three-bead bending potential. Their polymer chains contained either 150 or 300 beads; chains were estimated to be subject to an average of 6.5 or 13 entanglements, respectively. Bead motions were calculated using molecular dynamics; the constant-temperature condition was maintained by applying to each bead a random force and a velocity-determined friction force. The calculations of $g_2\left(t\right)$ were part of a larger paper in which the relaxation modulus $G\left(t\right)$ and the potential energy of the confining tube, as determined from the excursion distance of each bead from its primitive path, were calculated.

Figure 15 shows $g_2\left(t\right)$ and its first logarithmic derivative for the two polymers studied by Zhou and Larson. $K_1$ has the same dominant features for the $N=150$ and $N=300$ polymers, namely at short times $K_1$ is relatively large (0.69 or 0.61, respectively), as time advances $K_1$ falls and reaches a minimum (0.08 or 0.23, respectively), and at longer times $K_1$ increases. For the $N=150$ polymer, near time $t=5000$, $K_1$ has a local minimum and is then nearly constant, with $K_1\in (0.31-0.34)$. There are several other narrow regions where $K_1$ could be said to be nearly constant, notably for the $N=300$ polymer near times 400-900, with $K_1\approx 0.41$.

Zhou and Larson compare their results with a theoretical model that predicts power-law behavior. They fit sections of $g_2\left(t\right)$ to $t^{1/4}$ and $t^{1/2}$ power laws, finding regions where the measured $g_2\left(t\right)$ is not transparently inconsistent with these power laws. However, on a log-log plot power laws are straight lines, while as seen above $g_2\left(t\right)$ for each polymer follows a smooth curve with a continuously changing slope. From our analysis, their interpretation of $g_2\left(t\right)$ corresponds to a valid best fit of their results to a particular theoretical model, even though the results might better be described elsewise.

Moreno and Colmenero [32] report simulations of a binary blend of two species of bead-spring polymers. Each polymer molecule contained ten beads. Beads interacted via an $r^{-12}$ and an $r^2$ potential with a cutoff distance. Springs were represented with a FENE potential. The A polymer beads were 60% larger than the B polymer beads and thus had appreciably slower dynamics. Mean-square displacements were reported for the pure B polymer, and for the B polymer in a blend with B mole fraction $x_B=0.3$, at temperatures between 0.4 and 1.5 in natural units. Periodic boxes contained between 250 and 600 chains.

Our fits to Moreno, et al.,’s data on the pure B polymer melts, as seen in Figure 16, show systematic changes with decreasing temperature. At the largest temperature studied, $T=1.5$, $g_1\left(t\right)$ as seen in Figure 16a has an initial rapid increase and then climbs approximately as a power-law in t, the exponent being $0.64\pm 0.02$. At the longest times reported, $g_1\left(t\right)$ increases more rapidly, the exponent of an inferred power-law description reaching 0.95 at long times. At lower temperatures, there are qualitative changes in the form of $g_1\left(t\right)$. At temperatures $T=1.0$ and below, the slope $K_1$ decreases to a local minimum and then increases again. At the local minimum, for $T=1.0$, $K_1=0.45$; with decreasing temperature $K_1$ at the minimum decreases monotonically. At times longer than the minimum in $K_1$, at larger temperatures $K_1\approx 2/3$ is nearly constant. At lower temperatures and times beyond its minimum, $K_1$ is not a constant; it instead increases monotonically with increasing time until a maximum is reached.

Moreno, et al.’s description of their results on the pure B polymer is qualitatively the same as ours. With decreasing temperature, at shorter times "...a plateau is developed...", that is, there is a region of small slope that is increasingly manifested. At longer times, we find an area with a slope close to 0.65, as already described by Moreno, et al., though at the lowest temperatures the increase in slope from the plateau area is extended, so that the region with $K_1\approx 2/3$ is reduced to a local maximum in $K_1$ corresponding to an inflection point in $g_1\left(t\right)$. Around that local maximum, a slope close to $2/3$ that could be approximated as a naked-eye near-straight line is seen.

The time dependence of $g_1\left(t\right)$ in the $x_B=0.3$ blend, as seen in Figure 17, has a somewhat simpler structure than does $g_1\left(t\right)$ in the pure melt. After an initial increase, for temperatures $T\ge 1.0$, $K_1$ becomes nearly constant, changing almost not at all. In this temperature range, for times $3\le t\le 100$$g_1\left(t\right)$ does follow a power law, with exponents 0.63 or 0.56, respectively. At temperatures $T\le 1.0$ and a broad range of times, we find that $K_1$ increases by perhaps 0.1 or 0.2 with increasing time, but is clearly not a constant. At lower temperatures, $g_1\left(t\right)$ therefore does not follow a power law in time.

Moreno, et al.’s [32] description of their results on the $x_{\mathrm{B}}=0.3$ blend is broadly consistent with what is seen in Figure 17. There is an early ballistic regime followed by a rollover to a much shallower slope. Moreno, et al., present straight-line descriptions of $g_1\left(t\right)$ in the shallower regime with slopes decreasing from 0.635 to 0.415 as temperature was reduced. Our results are quantitatively fairly close to theirs, except that we find that not-quite-straight lines having small curvatures give a more accurate description of $g_1\left(t\right)$ in this regime.

Padding and Briels [33,34,35,36] report simulations of linear polyethylene melts with chain lengths extending from $C_{80}$ to $C_{1000}$, including calculations of $g_1(q,t)$, the dynamic structure factor $S(q,t)$, and Rouse mode amplitudes and time correlation functions. They used molecular dynamics and a united atom model with harmonic bead-bead stretch and bend potential energies, a dihedral angle potential, and van der Waals forces between atoms that were not adjoining atoms along each polymer chain. For some purposes united atoms were aggregated into blobs, typically of 20 atoms. Blob dynamics were inferred from united-atom dynamics; the purpose of the blobs was to allow simulation of much longer polymers than would otherwise be possible with then-available computing power. Of relevance here, Padding and Briels studied the time-dependent mean-square displacements of individual atoms $g_1\left(t\right)$, chain centers of mass $g_3\left(t\right)$, and individual blobs $g_{\mathrm{bl}}\left(t\right)$. Other dynamic properties were also considered but are outside the scope of this paper.

Figure 17. A First part of Figure 17. Figure has been split to spread across two pages. (a) Mean-square displacement $g_1\left(t\right)$ of the B beads of a melt of an A-B polymer blend with B mole fraction $x_B=0.3$, based on Moreno, et al. [32] at temperatures (a) $1.5$, (b) $1.0$, (c) $0.75$, and (d) $0.60$. Heavy lines represent Moreno, et al.’s data, circles are eighth-order polynomial fits, and thin solid lines are first derivatives $K_1$ of the fits.

Figure 17. A First part of Figure 17. Figure has been split to spread across two pages. (a) Mean-square displacement $g_1\left(t\right)$ of the B beads of a melt of an A-B polymer blend with B mole fraction $x_B=0.3$, based on Moreno, et al. [32] at temperatures (a) $1.5$, (b) $1.0$, (c) $0.75$, and (d) $0.60$. Heavy lines represent Moreno, et al.’s data, circles are eighth-order polynomial fits, and thin solid lines are first derivatives $K_1$ of the fits.

Figure 17. B Second part of Figure 17. Figure has been split to spread across two pages. Mean-square displacement $g_1\left(t\right)$ of the B beads in an A-B polymer blend melt with $x_B=0.3$, based on Moreno, et al. [32] at temperatures (e) $T=0.50$, (f) $T=0.45$ , and (g) $T=0.40$. Heavy line show Moreno, et al.’s data, circles represent the eighth-order polynomial fits, and thin lines are the first derivatives $K_1$ of the polynomials.

Figure 17. B Second part of Figure 17. Figure has been split to spread across two pages. Mean-square displacement $g_1\left(t\right)$ of the B beads in an A-B polymer blend melt with $x_B=0.3$, based on Moreno, et al. [32] at temperatures (e) $T=0.50$, (f) $T=0.45$ , and (g) $T=0.40$. Heavy line show Moreno, et al.’s data, circles represent the eighth-order polynomial fits, and thin lines are the first derivatives $K_1$ of the polynomials.

Their first paper [33] considered simulations of an n-C₁₂₀H₂₄₂ melt. Their results for $g_1\left(t\right)$, $g_3\left(t\right)$, and $g_{\mathrm{bl}}\left(t\right)$ appear in Figure 18. Figure 18, showing $g_1\left(t\right)$, is noteworthy for identifying a temporal region in which the mean-square displacement clearly does increase as a power law in time, namely, for $6\le t\le 220$, $g_1\left(t\right)\sim t^{\alpha }$ for $\alpha \approx 0.66$. At larger t, $\alpha$ first decreases to $\approx 0.53$ and then increases to $\approx 0.68$. In contrast, $g_3\left(t\right)$ only approaches power-law behavior, namely for center-of-mass motion $K_1$ drops rapidly to 0.75, then over close to two decades in time increases to 0.84 and finally increases to 1.0 or a bit larger. The behavior of $g_{\mathrm{bl}}\left(t\right)$ is qualitatively similar to that of $g_1\left(t\right)$, namely after an initial rapid decrease in $K_1$, there is an extended region in which $\alpha \approx 0.78$ is encountered.

Figure 18. Mean-square displacements of (a) individual atoms, (b) chain centers of mass, and (c) individual blobs $g_{\mathrm{bl}}\left(t\right)$ for a simulated $T=450$K melt of n-C₁₂₀H₂₄₂. Heavy lines indicate Padding and Briel’s data [33], circles represent eighth-order polynomial fits, and thin solid lines represent first derivatives $K_1$ of the polynomial fits.

Figure 18. Mean-square displacements of (a) individual atoms, (b) chain centers of mass, and (c) individual blobs $g_{\mathrm{bl}}\left(t\right)$ for a simulated $T=450$K melt of n-C₁₂₀H₂₄₂. Heavy lines indicate Padding and Briel’s data [33], circles represent eighth-order polynomial fits, and thin solid lines represent first derivatives $K_1$ of the polynomial fits.

In a further paper, Padding and Briels [35] report mean-square displacements of blobs – coarse-grained groups of 20 adjoining monomers – in simulated melts of linear polyethylene. Coarse-graining introduces complications. Most of the mechanical variables in the monomers of each blob must become bath variables, in the sense of the Mori-Zwanzig formalism, so that the blob equations of motion gain frictional and random-thermal forces. The blob-blob potential energy, determined as the blob-blob potential of average force from a united-atom simulation of the same system, is soft, so an alternative procedure was introduced to create the chain non-crossing constraint.

Figure 19 shows their results for $g_{\mathrm{bl}}\left(t\right)$ for the 4-, 6-, and 50-blob chains. Their original figure (Figure 13 in their paper) also shows results for 10- and 20-bead polymer melts, but we encounter a limit of our procedure, namely if the data points are sufficiently overlapped, it becomes impossible to digitize the graph reporting the data. Having said that, $g_{\mathrm{bl}}\left(t\right)$ curves for the three chain lengths are actually rather similar. The logarithmic derivative $K_1$ for each of the three $g_{\mathrm{bl}}\left(t\right)$ curves begins at short times with a value close to 1.0, decreases to a minimum at times near $t=1000$, and then increases monotonically to a value close to 0.9. $K_1$ at the minimum is 0.60, 0.53, or 0.38, respectively, for the three chain lengths. Corresponding to the minimum in $K_1$, $g_{\mathrm{bl}}\left(t\right)$ has an inflection point. Except perhaps for the 4-blob chains at early times ($t<60$), there are clearly no power-law regions at any time for any chain length. At times beyond the inflection point, $K_1\left(t\right)$ of the two longer chains shows some additional structure, namely a region where it increases relatively slowly with increasing time.

In presenting their results, Padding and Briels [35] add to their log-log plot four straight lines to guide the eye, straight lines corresponding to power-law regimes, the labeled slopes being 1, 0.5, 0.4 (for the 50-blob chain) and 1. We find that the initial and terminal slopes are in reasonable agreement with Padding and Briels, namely $K_1$ approaches 1. Where Padding and Briels suggest a power-law region with a slope of 0.4, we identify an inflection point whose minimum slope is $K_1=0.38$. The suggested $K_1=0.5$ line is in a region where there are few data points, and in which we find a slope that changes continuously, but their suggested line would be a tangent at some ${10}^2\le t\le {10}^3$.

4. Discussion

This review considers simulations of polymer melts that determined polymer mean-square displacements $g\left(t\right)$. The time dependence of $g\left(t\right)$ is of specific interest because widely-used theoretical models have as a core prediction that $g\left(t\right)$ has power-law behavior $g\left(t\right)\sim t^{\alpha }$. Equations 5–12 show the predicted $\alpha$ in different time regimes. A particular focus of simulations has been a search for the $\alpha =1/4$ exponent of 10.

We applied a new mathematical approach to test these predictions. As described in Section II, for each simulated system we determined the (time-dependent) logarithmic derivative $K_1$ of $g\left(t\right)$. At each t, $K_1$ from our approach gives an effective local value for $\alpha$. In any of the hypothetical power-law regimes, $K_1$ would a constant over an appreciable range of times. With very few exceptions, such a phenomenon—a constancy of $K_1$—is not found. Instead, $K_1$ consistently shows a single more-or-less deep minimum, corresponding to a saddle point in the time dependence of $g\left(t\right)$.

We conclude that computer simulations serve to reject core predictions of the deGennes-Doi-Edwards description of a polymer melt. Efforts to extract power-law exponents $\alpha$ from simulations of $g\left(t\right)$ have been challenging. Our analysis indicates that the efforts failed because the hypothesized power-law behavior was not there to begin with.

An obvious objection to our conclusion is that, even with modern computational techniques, simulations can only examine relatively short polymers, polymers that are not heavily entangled. Indeed, the longest polymers examined here have as any as 2000 beads, perhaps corresponding to as many as forty entanglements. Of course, the deGennes-Doi-Edwards description does not actually indicate how many entanglements per chain are needed for the tube-reptation model to be applicable. There are also literature discussions as to what an entanglement is. A critic might therefore propose that the chains considered here were all unentangled. Nonetheless, as seen in equations 5–12, the deGennes-Doi-Edwards description also predicts the behavior of short chains, where there are few or no entanglements, and the behavior of long chains at short times, times too early for chain entanglements to affect polymer motion. These short-chain ( and perhaps short-time) predictions are indeed tested by our analysis. The predictions are incorrect; there is almost never power-law behavior; when there is, the exponent does not have the predicted value.

It should be entirely unsurprising that $g_1\left(t\right)$ and $g_3\left(t\right)$ of short polymers do not have the properties predicted in equations 5–8. These equations describe polymers that move in accordance with the Rouse model. However, as we have previously shown in an extensive review [12], for polymers in simulated melts the Rouse amplitudes $X_p\left(t\right)$ and their time correlation functions $\langle X_p\left(0\right)X_q\left(t\right)\rangle$ have none of the properties associated with a Rouse-model polymer. In particular:

(i) Mean-square Rouse mode amplitudes $\langle {\left(X_p\left(0\right)\right)}^2\rangle$ do not scale as predicted by the Rouse model;

(ii) Rouse-mode time correlation functions $\langle X_p\left(t\right)X_p\left(0\right)\rangle$ do not decay with time as the exponentials predicted by the Rouse model;

(iii) Polymer bead displacements are not described by the independent Gaussian random processes that are assumed by the Rouse model; and

(iv) Rouse mode amplitudes are not always independent, i.e., it is not always true that $\langle X_p\left(0\right)X_q\left(t\right)\rangle =0$ for $p\ne q$, as required by the Rouse model.

Proposals that polymer motions in the melt are described by the Rouse model, perhaps only at short times, are systematically inconsistent with extensive published simulations. This result should not be surprising, because the underlying Hamiltonian (properly, the underlying Rayleigh dissipation function) for short polymer chains in their melt is not the function assigned by Rouse. Rouse assumed that the forces on the beads of a polymer could be modelled as a series of independent Gaussian random processes, one for each bead. This assumption has been checked simulationally [12]. It is false. The displacement distribution function $P(\Delta r,t)$ for the individual beads, the probability of finding a displacement $\Delta r$ during a time t, is in general not a Gaussian; instead, the non-Gaussianity parameter ${\alpha }_2\left(t\right)$ for the polymer bead displacements is found to be non-zero and time-dependent.

In almost all cases, we find that power-law behavior is absent. Theoretical models that predict that power-law behavior is important are rejected by our analysis. Now, it might be proposed that the predicted power-law regimes are separated by transition regimes in which power-law behavior is absent. It could be proposed that the transition regimes between power-law regimes are broad. Indeed, the latter proposal could be said to be supported by our analysis, namely, the transition regimes could have been more-or-less completely swallowed up the regions where power-law activity is found, reducing the power-law behavior to the slope at an inflection point. In that case, however, attempts to calculate any hypothesized power-law exponents are pursuing a theoretical fata morgana. A satisfactory theory of polymer melt dynamics would in this case need to focus on the broad transition regimes, not on the hypothesized regimes where power laws are encountered.