1. Introduction {#sec1-polymers-11-00476}
===============

The study of flow properties of polymeric solutions and melts has a rich history of perplexing the physicists and engineers who have endeavored to understand and model the many and varied physical responses of these complex fluids to an imposed flow field. In particular, the description of fast flows of macromolecular fluids has proven to be a difficult challenge. Although many continuum level theories have proven capable of describing gross rheological data in the linear and weakly nonlinear viscoelastic flow regimes (i.e., at low to intermediate values of the strain rate relative to a characteristic relaxation time of the fluid), most of these have not been able to provide a quantitative description of the flow properties of solutions and melts at high flow strength. There are many possible reasons one could cite to explain this state of dysfunction, but the overall reason is abundantly clear: for polymeric fluids experiencing strong flow conditions, all of the physical and dynamical phenomena occurring within these materials have not been understood and accounted for in the prevailing mathematical models.

Developing reliable mathematical models necessarily depends upon complementary experimentation. A debilitating feature of rheological experimentation, however, is that these seemingly simple experiments typically only provide bulk-scale measurements that have effectively been averaged over macroscopic length and time scales. As a consequence, any dynamic behavior that is of much shorter length and time scale than those of the measuring instrument are effectively washed out of the system response, even though they contribute to the overall response. Therefore, for much of the 20th century, rheologists had little in the way of small length and time scale information to guide attempts at improved mathematical modeling.

The 21st century is proving to be a golden age of rheological discovery. New experimental methods have been developed which are beginning to tap into small time and length scale phenomena that have a dramatic impact on the bulk rheological response of a polymeric liquid, particularly under conditions of strong flow. Furthermore, the present century has seen the rise of a new form of scientific exploration; i.e., *virtual experimentation*. Advances in computational algorithms and efficiency have led to a new paradigm in experimentation that, under the right circumstances, can lead to a powerful new means to probe the small length and time scale phenomena that dominate the bulk rheological responses of polymeric fluids under strong flow conditions.

The primary advantage of virtual experimentation of an atomistically detailed polymer chain over experiment is that every chain within the sample can be examined individually, not simply the bulk rheological or microstructural response. This allows much more detailed information to be gleaned from the simulation with respect to the experiment, as statistically meaningful correlations can be established via ensemble averaging of the dynamical behavior of each individual chain. Additionally, simulations are readily amenable to topological analysis, extending equilibrium properties such as tube diameter, primitive path length, and number of entanglements to nonequilibrium flow situations \[[@B1-polymers-11-00476],[@B2-polymers-11-00476],[@B3-polymers-11-00476],[@B4-polymers-11-00476],[@B5-polymers-11-00476],[@B6-polymers-11-00476],[@B7-polymers-11-00476]\]. Certainly, bulk-averaged properties, such as the conformation and stress tensors, can still be calculated, but also with the ability to examine the effects of short timescale individual chain dynamics upon them. Ultimately, more and better information at the microscopic scale should lead to better rheological and microstructural models of polymeric liquids under flow.

Recent evidence collected via virtual experimentation of monodisperse atomistic melts has demonstrated that a flow-induced disentanglement of polymer macromolecules occurs at high strain rates in steady shearing flow. This reduction in interchain constraints leads to the onset of individual molecular retraction and rotation cycles, which occur within oriented tube-like structures composed of the highly-extended surrounding chain molecules. Eventually, the tube network disintegrates as the chains become effectively disentangled, allowing them to tumble with characteristic frequencies similarly to corresponding macromolecules in dilute solution. This new phenomenon has been observed via nonequilibrium molecular dynamics (NEMD) simulations of molten polyethylenes in the unentangled and moderately-entangled molecular-weight regimes (i.e., polyethylenes ranging up to C~700~H~1402~) \[[@B1-polymers-11-00476],[@B2-polymers-11-00476],[@B3-polymers-11-00476],[@B4-polymers-11-00476],[@B8-polymers-11-00476],[@B9-polymers-11-00476],[@B10-polymers-11-00476],[@B11-polymers-11-00476],[@B12-polymers-11-00476],[@B13-polymers-11-00476],[@B14-polymers-11-00476],[@B15-polymers-11-00476]\]. This unexpected observation from atomistic simulations has already been hypothesized to explain some of the difficulties that manifest in flow models for high strain-rate flows \[[@B16-polymers-11-00476],[@B17-polymers-11-00476],[@B18-polymers-11-00476]\].

In the present contribution, prior results of unentangled (liquids ranging in molecular weight roughly up to C~250~H~502~), mildly entangled (C~400~H~802~), and moderately entangled (C~700~H~1402~) polyethylene melts are extended to a highly entangled system, C~1000~H~2002~, thus completing the entire suite of virtual experiments of flexible, monodisperse linear macromolecular fluids ranging from unentangled alkane liquids to highly entangled polyethylene melts. Hence this publication presents the final piece of the puzzle to those that preceded it, providing a full description of the rich, complex dynamical behavior and the underlying physical mechanisms that give rise to it as macromolecular chain length increases from several carbon units up to 1000. Moreover, this work extends prior studies that were focused on steady-state dynamics to the transient response of these entangled liquids under startup of shear conditions. The data and analysis presented in the remainder of this article will help enable physicists and engineers to develop new and improved models for the bulk rheological behavior of these macromolecular fluids covering all the relevant length and time scales.

2. Simulation Methodology {#sec2-polymers-11-00476}
=========================

Equilibrium and nonequilibrium molecular dynamics simulations of a monodisperse, linear, C~1000~H~2002~ melt were performed in the NVT ensemble at a constant density of 0.766 g/cm^3^ (corresponding to a pressure of 1 atm) and constant temperature of 450 K. Four different rectangular simulation cells were chosen for different shear rate ranges in order to minimize the computational cost by optimizing the simulation box size and number of particles. [Table 1](#polymers-11-00476-t001){ref-type="table"} summarizes the cell sizes in various directions as well as the number of particles and applicable $Wi$ range. In the nonlinear viscoelastic regime ($Wi > 1$), the box dimension in the flow direction ($x$) was larger than the dimensions in the gradient ($y$) and neutral ($z$) directions to ensure minimal system size effects at high shear rates where chains orient and stretch in the direction of flow. These dimensions were chosen based on the same considerations in terms of the chain end-to-end distance at different $Wi~$which were employed for a shorter C~700~H~1402~ chain liquid in prior work \[[@B2-polymers-11-00476]\]. The smallest simulation cell, containing 20,000 particles, was equilibrated for more than 8 times the longest relaxation (disengagement) time before any data were gathered for analysis. The simulation cells containing 40,000 and 60,000 particles were created by replicating the equilibrated small simulation cell respectively once and twice in the x-direction, then equilibrated for one disengagement time. The longest cell was created by replicating the equilibrated cell containing 60,000 particles, twice in the x-direction, and then equilibrated for 0.8 disengagement time. It should be mentioned that the transient data were obtained using only a single independent initial equilibrium configuration to minimize the computational cost; however, this configuration varied from one simulation cell size to another---see [Table 1](#polymers-11-00476-t001){ref-type="table"}. Although ideally such data should be collected using more than one independent initial configuration at each $Wi$, based on prior experience, this has only a slight effect on the data presented in this work, considering the sufficiently large number of particles in the simulation cells.

The Siepmann--Karaboni--Smit (SKS) united-atom potential model \[[@B19-polymers-11-00476]\] was used to quantify the energetic interactions between the atomistic components of the polyethylene liquid. This is the same potential model employed in many other prior simulation studies \[[@B1-polymers-11-00476],[@B2-polymers-11-00476],[@B3-polymers-11-00476],[@B8-polymers-11-00476],[@B10-polymers-11-00476],[@B11-polymers-11-00476],[@B12-polymers-11-00476],[@B13-polymers-11-00476],[@B14-polymers-11-00476],[@B17-polymers-11-00476],[@B18-polymers-11-00476],[@B20-polymers-11-00476],[@B21-polymers-11-00476],[@B22-polymers-11-00476],[@B23-polymers-11-00476]\] to represent energetic interactions between either -CH~3~ for the end-groups of the chains or -CH~2~- groups for interior carbon atoms along the chain backbone. (Please refer to one of the references cited above for a detailed discussion of the SKS model equations and parameters.)

The NEMD equations of motion were used to perform the NEMD simulations, which were maintained at a constant temperature of 450 K using a Nosé--Hoover thermostat \[[@B24-polymers-11-00476],[@B25-polymers-11-00476],[@B26-polymers-11-00476],[@B27-polymers-11-00476],[@B28-polymers-11-00476],[@B29-polymers-11-00476],[@B30-polymers-11-00476],[@B31-polymers-11-00476]\]. The set of evolution equations for the particle positions and momenta were integrated within the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) environment, which is implemented using the p-SLLOD equations of motion \[[@B29-polymers-11-00476],[@B30-polymers-11-00476],[@B31-polymers-11-00476]\]. (Note that for steady-state and startup shear flow as considered herein the SLLOD and p-SLLOD algorithms are the same.) Boundary conditions were periodic at all box surfaces with a deforming simulation cell in the $x$ direction. The equations were integrated using the reversible-Reference System Propagator Algorithm (r-RESPA) \[[@B32-polymers-11-00476]\] with two different time steps. The long time step was 4.70 fs, which was used for the slowly varying nonbonded Lennard-Jones interactions, and the short time step was 1.176 fs (one-fourth of the long time step) for the rapidly varying forces including bond-bending, bond-stretching, and bond-torsional interactions. The relaxation time of the thermostat was set equal to 100 times the long time step. These time steps are longer than those used in many of the prior studies \[[@B1-polymers-11-00476],[@B2-polymers-11-00476],[@B3-polymers-11-00476],[@B4-polymers-11-00476],[@B8-polymers-11-00476],[@B11-polymers-11-00476],[@B12-polymers-11-00476],[@B13-polymers-11-00476],[@B14-polymers-11-00476],[@B17-polymers-11-00476],[@B18-polymers-11-00476],[@B21-polymers-11-00476],[@B22-polymers-11-00476]\]; however, a series of test simulations were performed at various *Wi* to ensure that the new (longer) time steps produced statically equivalent results as the prior (shorter) time steps. Furthermore, these time steps have been used successfully in recent NEMD studies of planar elongational flows of entangled polyethylene melts \[[@B20-polymers-11-00476],[@B23-polymers-11-00476]\]. Without this modification of the time steps, the simulations reported in this article would have been computationally intractable.

A wide range of Weissenberg numbers was examined over the interval \[0; 0.01, 11,700\], corresponding to the quiescent system and shear rates within the range $2.2 \times 10^{3}\mathbf{s}^{- 1} \leq \overset{˙}{\gamma} \leq 2.2 \times 10^{9}\mathbf{s}^{- 1}$. The topological analysis was performed using the Z1 code developed by Kröger \[[@B5-polymers-11-00476]\], which reduces atomistic configurations to a primitive path network in which the chains are not allowed to cross each other as the algorithm simultaneously minimizes the contour length of each polymer molecule \[[@B6-polymers-11-00476]\]. This method uses geometrical methods rather than dynamical algorithms to minimize the contour lengths of primitive paths in the most computationally efficient manner. The code further defines positions of kinks along the 3-dimensional primitive path of each chain, which are assumed to be roughly proportional to the number of entanglements per chain. Results of the code can be used to interpret other important reptative parameters, such as the effective tube diameter and entanglement strand length. The Z1 code has been compared with other topological analysis techniques by Shanbhag and Kröger \[[@B7-polymers-11-00476]\].

3. Results and Discussion {#sec3-polymers-11-00476}
=========================

3.1. Quiescent Properties {#sec3dot1-polymers-11-00476}
-------------------------

Equilibrium properties of the system can be calculated from the simulation results and compared with the predictions of reptation theory. The ensemble average squared end-to-end distance, $\left\langle R^{2} \right\rangle$, and radius of gyration, $\left\langle R_{g}{}^{2} \right\rangle$, were calculated as 20,107 $\mathring{\mathrm{A}}^{2}$ and 3353 $\mathring{\mathrm{A}}^{2}$, respectively, directly from the equilibrium simulation data. The theoretical fully extended chain end-to-end distance, $\left| R \right|_{\max}$, for a C~1000~H~2002~ molecule is 1290.2 $\mathring{\mathrm{A}}$. These values may be used to approximate the Kuhn length as $= \left\langle R^{2} \right\rangle/\left| R \right|_{\max} = 15.58\ \mathring{\mathrm{A}}$, and the number of Kuhn segments as $N = \left\langle R^{2} \right\rangle/b^{2} = 82.79 \approx 83$. Entanglement network properties were evaluated using the Z1-code \[[@B5-polymers-11-00476]\]. Specifically, the average primitive chain contour length, $\left\langle L \right\rangle = 508.6\ \mathring{\mathrm{A}}$, was obtained based on this analysis. These basic properties of the entangled liquid can be used in conjunction with reptation theory to estimate other (theoretical) system properties. The ensemble average entanglement density is thereby estimated as $$Z = \frac{\left\langle L \right\rangle}{a~} = \frac{Nb^{2}}{a^{2}} = \frac{\left\langle L \right\rangle^{2}}{\left\langle R^{2} \right\rangle~},$$ where $a$ is the tube diameter. From this expression, $Z = 12.9$ and $a = \left\langle L \right\rangle/Z = \left\langle R^{2} \right\rangle{}^{1/2}/Z^{1/2} = 39.5\ \mathring{\mathrm{A}}$. All of these values are in good agreement with the values estimated for C~400~H~802~ and C~700~H~1402~ molecules in prior work \[[@B2-polymers-11-00476],[@B3-polymers-11-00476],[@B18-polymers-11-00476]\]. The entanglement molecular weight for polyethylene at 443 K of $M_{e} = 1150$ g/mol was reported by Fetters, et al. \[[@B33-polymers-11-00476]\], which can be used to estimate an experimental entanglement density of $Z = M/M_{e} = 12.9$, and tube diameter $a = \left\langle R^{2} \right\rangle{}^{1/2}/Z^{1/2} = 40.6\ \mathring{\mathrm{A}}$. These values are in excellent agreement with the simulation results.

The diffusivity of the liquid can be evaluated from the slope of the chain center-of-mass mean-squared displacement (MSD) versus time. According to its definition, $D_{G}$ is $1/6$ of this slope at long times:$$D_{G} = \lim\limits_{t\rightarrow\infty}~\frac{1}{6t}\left\langle \left( {R_{G}\left( {t + \tau} \right) - R_{G}\left( \tau \right)} \right) \right.^{2}\rangle,$$where $R_{G}\left( t \right)$ is the position of the chain center of mass at time $t$. Using this method, the diffusivity of C~1000~H~2002~ is calculated as $1.28 \times 10^{- 12}{\ m}^{2}/s$. Note that this value makes it possible to calculate a key model parameter of the theory, the friction coefficient, $\xi$ \[[@B34-polymers-11-00476]\], $$D_{G} = \frac{k_{B}T~a^{2}}{3N^{2}\xi b^{2}},$$ where $k_{B}$ is the Boltzmann constant and $T$ is the absolute temperature.

According to reptation theory, the Rouse and disengagement timescales are governed by the expressions \[[@B34-polymers-11-00476]\] $$\tau_{R} = \frac{\xi N^{2}b^{2}}{3\pi^{2}k_{B}T~},$$ $$\tau_{d} = \frac{1}{\pi^{2}}\frac{\xi N^{3}b^{4}}{k_{B}T~a^{2}}.$$ Substituting $\xi$ from Equation (3) into Equation (4),$~\tau_{d}$ can be expressed as a function of $D_{G}$:$$\tau_{d} = \frac{1}{3\pi^{2}}\frac{R^{2}}{D_{G}}.$$ Hence, the theoretical value of the disengagement time (i.e., according to the equations of reptation theory after $D_{G}$ has been estimated from the simulations) is calculated as $5305$ ns. Note that from Equations (1), (4) and (5), the ratio $\tau_{d}/\tau_{R} = 3Z$, which leads to a theoretical Rouse time of $137$ ns.

The entanglement time is governed by the reptation-based Equation \[[@B3-polymers-11-00476]\] $$\tau_{e} = \frac{\pi}{36}\frac{R^{2}}{D_{G}Z^{3}},$$ which gives a value of 6.4 ns. This is very close to the values calculated for the C~400~H~802~ and C~700~H~1402~ liquids (5.1 ns and 6.4 ns, respectively). These values are consistent with theoretical arguments suggesting that the entanglement time is independent of the molecular weight of the polymeric liquid.

The characteristic relaxation times can also be estimated directly from the equilibrium simulation results using the characteristic breaks in the segmental mean-square displacement (MSD) plot versus time \[[@B34-polymers-11-00476]\]. The segmental MSD is defined as $\phi\left( t \right) = \langle\left( {r_{n}\left( {t + \tau} \right) - r_{n}\left( \tau \right)} \right)^{2}\rangle$, where $r_{n}$ is the position vector of the $n$-th monomer (i.e., the $n$-th -CH2- unit). In order to minimize chain-end effects, only the 500 monomers in the middle of each chain were included in these calculations. The details of the calculations are explained in prior publications \[[@B2-polymers-11-00476],[@B3-polymers-11-00476]\]. [Figure 1](#polymers-11-00476-f001){ref-type="fig"} displays these plots for very short times (a) and long times (b). As shown in the figure, the disengagement, Rouse, and entanglement times turn out to be, respectively, 5834 ns, 194 ns, and 2.7 ns. Both the Rouse and disengagement times are mildly overpredicted as compared to the theoretical values. Also, the ratio $\tau_{d}/3\tau_{R} = 10$, which is smaller than the expected theoretical value of *Z* = 12.9; this suggests that the Rouse time is significantly overpredicted by this method. On the other hand, $\tau_{e}$ is underpredicted compared to the theoretical value; however, it is in good agreement with the entanglement times calculated for the C~400~H~802~ and C~700~H~1402~ melts using the same method.

Another robust method for direct calculation of the disengagement time from the simulation data is to fit a sum of exponential functions to the autocorrelation function of the chain end-to-end unit vector, $\langle u_{i}\left( \tau \right) \cdot u_{i}{\left( {\tau + t} \right)\rangle} = \sum\limits_{i = 1}^{p}c_{i}\exp\left( {- t/\tau_{i}} \right)$, where the longest value of $\tau_{i}$ is considered as the disengagement time. $p$ is the minimum number of exponential terms (5 in this case) that results in the best fit (i.e., the closest coefficient of determination, R-squared, to unity, using a nonlinear least-squares method), and $c_{i}$ are fitting constants of order unity. The disengagement time, based on this method, is calculated to be 5270 ns, which agrees very well with the theoretical prediction (5305 ns).

[Figure 2](#polymers-11-00476-f002){ref-type="fig"} displays a log-log plot of the Rouse and disengagement times versus chain length (number of monomers per chain), as well as their relevant power-law fitting and exponents. The data for the C~400~H~802~ and C~700~H~1402~ liquids were obtained from prior work \[[@B2-polymers-11-00476],[@B3-polymers-11-00476]\]. These plots show that the power-law exponents for the disengagement time calculated from either the theoretical method or the fitting method are about $3.3 \pm 0.1$, in good agreement with experimental measurements for entangled polymers. This suggests that all physical phenomena, including contour length fluctuations (CLF), constraint release (CR), and of course reptation, are captured well by the simulations under quiescent conditions. One may expect a power-law exponent of 3.0 for the theoretical values of the disengagement time based on reptation model predictions alone; however, it should be noted that although Equations (3)--(5) were used for calculation of the theoretical characteristic times, the diffusivity (or equivalently the friction coefficient) was calculated from the simulation results and consequently includes all important physical phenomena, as explained earlier. In fact, the power-law exponent for the diffusivity itself is $- 2.3$, in excellent agreement with experimentally observed values \[[@B34-polymers-11-00476],[@B35-polymers-11-00476]\]. The same analysis is valid for the power-law exponent of the theoretical Rouse time, which scales as $N_{m}^{2.2}$ rather than the theoretical exponent of 2, which again should be attributed to CLF and CR effects \[[@B2-polymers-11-00476]\]. [Table 2](#polymers-11-00476-t002){ref-type="table"} summarizes the results of the calculations for the equilibrium characteristic relaxation times of the C~1000~H~2002~ melt obtained from the theoretical and MSD methods. In the rest of this chapter, we use the values $\tau_{d} = 5270$ ns (exponential method), $\tau_{R} = 137$ ns (theoretical method), and $\tau_{e} = 6.4$ ns (theoretical method) for the characteristic time scales of the C~1000~H~2002~ liquid (unless otherwise noted).

3.2. Steady Shear Flow Properties {#sec3dot2-polymers-11-00476}
---------------------------------

### 3.2.1. Steady-State Structural and Topological Properties {#sec3dot2dot1-polymers-11-00476}

The steady-state microstructural and topological properties of a C~1000~H~2002~ melt undergoing simple shear flow are qualitatively very similar to those of the C~400~H~802~ and C~700~H~1402~ liquids, which were discussed in detail in prior publications \[[@B2-polymers-11-00476],[@B3-polymers-11-00476],[@B4-polymers-11-00476],[@B18-polymers-11-00476]\]. These results are presented concisely herein; interested readers can refer to the cited references for more comprehensive discussions. Overall, steady-state shear properties of the C~1000~H~2002~ melt exhibit four distinct regions of behavior ($\overset{˙}{\gamma} < \tau_{d}^{- 1}$, $\tau_{d}^{- 1} < \overset{˙}{\gamma} < \tau_{R}^{- 1}$, $\tau_{R}^{- 1} < \overset{˙}{\gamma} < \tau_{e}^{- 1}$, $\tau_{e}^{- 1} < \overset{˙}{\gamma}$), as noted previously for the C~700~H~1402~ liquid \[[@B2-polymers-11-00476]\].

The probability distribution functions (PDFs) of the normalized end-to-end distance and the chain size (measured in terms of ensemble averages of chain end-to-end distance and six times the radius of gyration, respectively) are displayed for various values of $Wi$ in [Figure 3](#polymers-11-00476-f003){ref-type="fig"}a,b. In the linear viscoelastic regime ($Wi \leq 1$), the PDFs are Gaussian and remain essentially unchanged from the quiescent state. The ensemble averages of the squared end-to-end distance and (6 times the) radius of gyration also remain constant and almost equal to each other in this regime. This suggests that the flow is too weak to significantly perturb the global molecular sizes. Keep in mind that the timescale of the flow is larger than the disengagement time (i.e., $\overset{˙}{\gamma} < \tau_{d}^{- 1}$), implying that the constituent macromolecules have ample time for diffusive action to maintain their quiescent configurational properties even though the overall tube network begins to orient along a preferred direction in the shear plane relative to the direction of flow. Note that the ratio $\langle R^{2}\rangle/\langle R_{g}^{2}\rangle$ approaches the theoretical value of 6 for long flexible Gaussian chains.

[Figure 3](#polymers-11-00476-f003){ref-type="fig"}c displays the ensemble average orientation angle, $\langle\theta\rangle$, as a function of $Wi$. $\langle\theta\rangle$ is calculated as the angle between the principal eigenvector of the ensemble average of the unit end-to-end vector dyadic product, $\langle u_{i}u_{i}\rangle$, and the flow ($x$) direction. The orientation angle decreases from the zero-shear-rate limit of 45° (not shown in the figure) to about 30° at $Wi = 1$. Finally, the tube stretch is shown as a function of $Wi$ in [Figure 3](#polymers-11-00476-f003){ref-type="fig"}d. The tube stretch is defined as the ratio $\lambda = \langle L\rangle/L_{0}$, where $L_{0}$ is the quiescent primitive path length. Both $\langle L\rangle$ and $L_{0}$ are calculated using the Z1 code. No chain stretch is observed in the linear viscoelastic region, as expected.

As the flow enters the weakly nonlinear regime, $\tau_{d}^{- 1} < \overset{˙}{\gamma} < \tau_{R}^{- 1}$ (or equivalantly $1 < Wi \leq 38$), the orientation angle drops dramatically to values smaller than 5° and plateaus around 1--2° at higher $Wi$. The PDF of the end-to-end distance begins deviating from the equilibrium Gaussian distribution by developing a tail at higher values of $\left| R \right|/\left| R \right|_{\max}$, indicating that a portion of the macromolecules have become partially extended by the applied flow. Notably, the PDF peak is still approximately at the same location as the equilibrium distribution, which suggests that the overall conformation of a significant number of chains has not yet been perturbed. The growth of molecular size and the deviation from Gaussian behavior can also be inferred from [Figure 3](#polymers-11-00476-f003){ref-type="fig"}b, especially for $Wi > 10$ where $\langle R^{2}\rangle~$and $6\langle R_{g}^{2}\rangle$ begin to diverge. (Note that there is no theory which indicates these two quantities are equivalent under flow conditions.) Interestingly, the tube network also begins to extend moderately in in this shear-rate region ([Figure 3](#polymers-11-00476-f003){ref-type="fig"}d). This is an important observation because it contradicts the common notion of tube-based models that no stretching occurs for $\overset{˙}{\gamma} < \tau_{R}^{- 1}$. Quantitatively, [Figure 3](#polymers-11-00476-f003){ref-type="fig"}d indicates that tubes are stretched about 16% at $\left. \overset{˙}{\gamma} \right.\sim\tau_{R}^{- 1}$, which is not negligible although just a fraction of the maximum theoretical tube stretch, $\lambda_{\max} = 2.77$.

The third shear-rate regime of dynamical behavior is the range $\tau_{R}^{- 1} < \overset{˙}{\gamma} < \tau_{e}^{- 1}$ (approximately $40 \leq Wi < 800$). Within this region, vorticity excursions start playing an important role in the system properties. Brownian fluctuations caused by the vorticity of the shear field lead to random excursions of the chain ends outside of the confining tubes; some of these excursions, especially those with shear-plane projections that possess negative orientation angles relative to the flow direction, induce rotation and retraction quasi-periodic tumbling cycles of the individual molecules at moderate and high shear rates similar to those observed in previous work \[[@B2-polymers-11-00476],[@B3-polymers-11-00476],[@B8-polymers-11-00476],[@B14-polymers-11-00476],[@B18-polymers-11-00476],[@B22-polymers-11-00476]\]. A typical cycle begins as a chain molecule stretches and aligns in the flow direction (see [Figure 4](#polymers-11-00476-f004){ref-type="fig"}c). At this point, due to the flow vorticity, chain ends fold backward along the spine of the molecule and slide toward the middle of the chain until the molecule collapses into a compressed configuration. Then the orientation of the chain flips as the chain ends cross and the molecule unravels until it adopts a stretched conformation again that concludes a half cycle. At the lower end of this range ($Wi = 40$ to 100), the cycle is very irregular, almost chaotic. Here, the macromolecules will reside in the compressed state for a long period of time (see [Figure 4](#polymers-11-00476-f004){ref-type="fig"}a). Under this condition, the chain ends, which are typically very close to each other, exhibit a wagging behavior due to Brownian motion, passing each other back and forth multiple times before the molecule begins to reextend. As a consequence, the orientation angle of the chain end-to-end vector, $\theta_{ete}$, oscillates haphazardly between $- 90{^\circ}$ and $90{^\circ}$ as evident in [Figure 4](#polymers-11-00476-f004){ref-type="fig"}a. The orientation angle of the chain primary axis, $\theta_{pa}$, however, does not oscillate as much as $\theta_{ete}$ suggesting that the body of the molecule does not wag like a solid object. (The primary axis of the molecule is defined as the eigenvector corresponding to the largest eigenvalue of the molecule gyration tensor.) Yet, $\theta_{pa}$ changes rapidly between positive and negative values when the molecule is in a collapsed and highly compressed state, indicating that the coiled chains wag for some indefinite period of time before they begin to unravel.

In the uppper portion of the range $\tau_{R}^{- 1} < \overset{˙}{\gamma} < \tau_{e}^{- 1}$ (i.e., $100 < Wi < 800$), the dynamical behavior of the macromolecules is much more regular and resembles the tumbling behavior observed in prior work \[[@B2-polymers-11-00476],[@B3-polymers-11-00476],[@B8-polymers-11-00476],[@B14-polymers-11-00476],[@B18-polymers-11-00476],[@B22-polymers-11-00476]\]. During a typical cycle, the chain end-to-end distance varies dramatically from high values associated with the stretched configurations to values that are even smaller than the average equilibrium end-to-end distance. This is manifested in the wide non-Gaussian bimodal probability distribution function at this flow regime, as displayed in [Figure 3](#polymers-11-00476-f003){ref-type="fig"}a. Specifically, the peak at low values of $\left| R \right|/\left| R \right|_{\max}$ shifts to the left as $Wi$ increases and occurs at extensions smaller than the equilibrium peak, indicating the increasing population of the collapsed configurations during the course of the tumbling cycle. At the same time, the ensemble average molecule size ([Figure 3](#polymers-11-00476-f003){ref-type="fig"}b) and tube stretch ([Figure 3](#polymers-11-00476-f003){ref-type="fig"}d) increase with $Wi$ in this flow region. Based on theoretical arguments, this is the region wherein tube stretch becomes significant. As mentioned earlier, [Figure 3](#polymers-11-00476-f003){ref-type="fig"}d shows that tube stretching begins at lower flow strength than theoretically expected; however, tube stretch in the third flow region is apparently of a different nature than within the second flow regime. In the third region, $\lambda$ scales as $\lambda^{0.04}$ while this power-law exponent is 0.03 in the second range, $\tau_{d}^{- 1} < \overset{˙}{\gamma} < \tau_{R}^{- 1}$. This suggests that the tube stretch is influenced by the tumbling dynamics of the individual macromolecules, and that it has an influential contribution to the shear stress and constitutes a major relaxation mechanism in this intermediate flow strength regime. Note that the time-average orientation angle of the molecules is very close to its plateau value in the third region and does not change significantly, indicating that the chain end-to-end vectors are almost completely aligned in the flow direction on the molecular length scale, although not necessarily on the tube segment length scale.

The fourth and final flow regime is the strong flow region where $\overset{˙}{\gamma} > \tau_{e}^{- 1}$, approximately $Wi > 800$. Although the molecules continue to stretch in this region ([Figure 3](#polymers-11-00476-f003){ref-type="fig"}b,d), the molecular size and the tube stretch ultimately attain plateau values, which are significantly smaller than their corresponding maximum theoretical values. The tube stretch profile has an inflection point around $\left. \overset{˙}{\gamma} \right.\sim\tau_{e}^{- 1}$ where the curvature changes from positive to negative. This signals a new regime where the tube stretch becomes saturated as chain rotation becomes the more dominant dynamic mechanism. The shape of the end-to-end distance distribution curve is also very different in this high $Wi$ region compared to that at lower $Wi$ regimes. Specifically, the distributions become relatively flat with a characteristic rotational peak at low $\left| R \right|/\left| R_{\max} \right|$ and a stretch peak that emerges at very high $Wi$. (See [Figure 3](#polymers-11-00476-f003){ref-type="fig"}a. The stretch peaks can also be easily recognized in C~400~H~802~ and C~700~H~1402~ systems \[[@B18-polymers-11-00476]\].) These flat distributions, which become wider as $Wi$ increases, are attributed to the more regular molecular rotation cycles at very high shear rates, as discussed by Nafar Sefiddashti, et al. \[[@B3-polymers-11-00476]\]. The skewed distributions within the intermediate $Wi$ regime suggest that during a rotation cycle individual molecules spend on average a longer time at collapsed (or less stretched) configurations than they do at relatively stretched configurations \[[@B18-polymers-11-00476]\], or that some of the chains have not yet stretched enough to begin their rotation cycles (see [Figure 4](#polymers-11-00476-f004){ref-type="fig"}a). Both cases lead to unbalanced lifetimes for various configurations, and consequently irregular rotation cycles. Within the high $Wi$ regime, on the other hand, molecules undergo more regular periodic cycles. Hence various configurations between a highly stretched chain and a tightly packed coil have fairly similar lifetimes or probabilities (see [Figure 4](#polymers-11-00476-f004){ref-type="fig"}b), which manifest in the flat probability distribution of the end-to-end distance \[[@B18-polymers-11-00476]\].

[Figure 5](#polymers-11-00476-f005){ref-type="fig"} displays the entanglement network properties of the C~1000~H~2002~ melt at various $Wi$. The ensemble average entanglement density and the probability distribution function for the entanglement density are displayed in [Figure 5](#polymers-11-00476-f005){ref-type="fig"}a,c. [Figure 5](#polymers-11-00476-f005){ref-type="fig"}b shows the tube diameter, determined as the step length of the primitive path $a = \langle L\rangle/\left( {\langle Z_{k}\rangle/2} \right)$ \[[@B2-polymers-11-00476],[@B3-polymers-11-00476]\]. The probability distribution function of the primitive path contour length is also shown in [Figure 5](#polymers-11-00476-f005){ref-type="fig"}d for various values of $Wi$. Note that the primitive path contour length, $\langle L\rangle$, is essentially commensurate with the tube stretch, $\lambda$ (see [Figure 3](#polymers-11-00476-f003){ref-type="fig"}d), which is the normalized primitive path contour length. These plots show that, within the linear viscoelastic regime, the entanglement network is practically unperturbed as compared to quiescent conditions. Specifically, the entanglement density and tube diameter do not change as the flow strength increases. The probability distribution function for the entanglement density, $P\left( Z_{k} \right)$, follows a Poisson distribution and is independent of $Wi$ in this regime. $P\left( L_{pp} \right)$ exhibits a similar behavior, except that it follows a Gaussian distribution. As $Wi$ increases and shear rate enters into the nonlinear viscoelastic regime, the tube network begins to lose entanglements. Notably, there is no sharp boundary between the second ($1 < Wi \leq 38$) and the third $58 \leq Wi < 800$) flow regimes, as discussed for the structural properties of the system. Rather, there is an initial stage of convective constraint release wherein the chains disentangle at a moderate rate in the region $1 < Wi < 500$ such that $\left. Z_{k} \right.\sim Wi^{- 0.07}$. Accordingly, the tube diameter increases moderately in this region. The probability distribution function for the entanglement density, $P\left( Z_{k} \right)$, shifts to the left with increasing $Wi$ as the chains disentangle. The shape of the distribution, however, remains approximately similar to that of the linear viscoelastic regime and still follows a Poisson distribution. On the other hand, $P\left( L_{pp} \right)$ shifts to the right and becomes wider (i.e., with a higher standard deviation); nevertheless, the distribution continues to follow a Gaussian distribution. These results suggest that, although by the end of this flow regime the system loses about 30% of its entanglements, the nature of the entanglement network does not change radically. Note that even at the highest shear rate within this regime, none of the chains has lost all of its entanglements. For instance, the curve for $Wi = 117$ in [Figure 5](#polymers-11-00476-f005){ref-type="fig"}c shows that all molecules possess 5 or more kinks.

At higher flow strength, (i.e., $Wi \geq 585$), the entanglement density begins to drop dramatically as $Wi$ increases with a power-law exponent of $\langle Z_{k}\rangle{}^{- 0.35}$. The tube diameter also increases substantially in this region, such that at $Wi = 2340$ the tube diameter grows almost as large as the molecular radius of gyration. This means that a molecule could effectively diffuse as far as its size without feeling the confining tube. This interpretation essentially questions the existence of the tube concept and of an entangled system. These subtleties can be understood by examining the probability distribution of the entanglement density. [Figure 5](#polymers-11-00476-f005){ref-type="fig"}c shows that the distributions begin to deviate somewhat from Poisson distributions. More importantly, these distributions suggest that, unlike before, in this flow region some of the molecules have lost all their entanglements and have become virtually unentangled. The distribution of the primitive path contour length also deviates considerably from the Gaussian distribution in this region. All of these observations suggest that the entanglement network is effectively destroyed due to the strong flow. This also might explain why the system behavior at such high shear rates resembles that of a dilute solution, as has been argued for shorter chain C~400~H~802~ and C~700~H~1402~ liquids \[[@B3-polymers-11-00476],[@B18-polymers-11-00476]\]. In this regime, the tumbling cycles are comparatively more regular, similar to those of dilute solutions. Tube stretch approaches its plateau value as macromolecular tumbling becomes the dominant dynamic mechanism.

An important characteristic of the entanglement network is the tube orientation tensor, $S$, that is one of the principal variables of many tube-based constitutive models. [Figure 6](#polymers-11-00476-f006){ref-type="fig"} displays the nonzero components of $S$ as functions of $Wi$ at steady state obtained from the NEMD data for the C~1000~H~2002~ melt. The average orientation tensor of the tube segments in this figure is defined as $S = \langle u^{t}u^{t}\rangle$, where $u^{t}$ is the unit end-to-end vector of an entanglement strand: knowing the positions of the entanglements (kinks) along the chain from the Z1 code, the end-to-end vectors of the entanglement strands can be easily identified and the appropriate ensemble averages of the components of the orientation tensor can then be readily calculated from the NEMD data. The component, $S_{xy}$, begins to increase at very low $Wi$ within the linear viscoelastic regime. This segmental orientation leads to an increase in the shear stress in this region, in agreement with theory. $S_{xy}$ passes through a maximum in the range $3 < Wi < 12$, which is somewhat higher than the theoretical prediction of $\left. Wi \right.\sim 1$. At higher shear rates, $S_{xy}$ decreases almost monotonically. Such behavior can lead to excessive shear thinning, as observed in vesrsions of the tube model that do not incorporate CCR, especially within the shear rate range $\tau_{d}^{- 1} < \overset{˙}{\gamma} < \tau_{R}^{- 1}$ or the plateau region where the tube stretch is insignificant. In fact, models that incorporate CCR predict a nearly constant $S_{xy}$, and consequently constant shear stress, in the plateau region in agreement with typical experiments. Hence, the decrease in $S_{xy}$ observed in the NEMD data calls into question the theoretical mechanism of CCR in some tube-based models like MLD. The diagonal components of $\mathbf{S}$ remain nearly constant in the linear viscoelastic regime and then diverge from their equilibrium value (\~0.33) as $Wi$ increases. At very high shear rates, i.e., $\overset{˙}{\gamma} > \tau_{e}^{- 1}$, the rate of change in these components increases significantly. This is the shear rate range wherein the entanglement network begins to disintegrate. Generally, the features of these plots are qualitatively similar to those of $\lambda$ and $\langle Z_{k}\rangle$ (see [Figure 3](#polymers-11-00476-f003){ref-type="fig"}d and [Figure 5](#polymers-11-00476-f005){ref-type="fig"}a), which could be indicative of an inherent connection between $\lambda$ and $\langle Z_{k}\rangle$ with the normal components of the tube segmental orientation tensor. This is specifically important from a modeling perspective as it suggests that the evolution equations for the tube stretch and entanglement density should be expressed in terms of the diagonal components of $S$ rather than the shear component.

[Figure 7](#polymers-11-00476-f007){ref-type="fig"} displays the important characteristic timescales of the C~1000~H~2002~ liquid as functions of $Wi$. These timescales are calculated based on fitting the autocorrelation function data of the end-to-end vector with the functional form $\langle u_{i}\left( \tau \right) \cdot u_{i}{\left( {\tau + t} \right)\rangle} = \exp\left( {- t/\tau_{d}} \right)\cos\left( {2\pi t/\tau_{rot}} \right)$. Hence, $\tau_{d}$ is the decorrelation time of the end-to-end vector, which is equal to the longest relaxation time (i.e., the disengagement time under quiescent conditions and within the linear viscoelastic regime). $\tau_{rot}$ quantifies the period of the rotation and retraction cycle of the macromolecules, assuming that the cycles are quasi-periodic. A characteristic time for the tumbling period can be defined conceptually as $\tau_{r} = \tau_{rot}/2\pi$ \[[@B2-polymers-11-00476]\], displayed as diamonds in [Figure 7](#polymers-11-00476-f007){ref-type="fig"}.

[Figure 7](#polymers-11-00476-f007){ref-type="fig"} indicates that $\tau_{d}$ does not change significantly within not only the linear viscoelastic region ($Wi \leq 1$) but also in the nonlinear regime for $1 < Wi \leq 12$. At higher shear rates, the longest relaxation time decreases with a power-law exponent of $- 0.71 \pm 0.06$; this is consistent with the scaling exponents of the C~400~H~802~ and C~700~H~1402~ liquids at high shear rates \[[@B2-polymers-11-00476],[@B3-polymers-11-00476]\]. Unlike for the C~1000~H~2002~ liquid, the relaxation times of C~400~H~802~ and C~700~H~1402~ decreased with shear rate at all $Wi > 1$, and hence a separate power-law exponent for the $\tau_{d}^{- 1} < \overset{˙}{\gamma} < \tau_{e}^{- 1}$ regime was reported in prior work \[[@B2-polymers-11-00476],[@B3-polymers-11-00476],[@B18-polymers-11-00476]\]. Nevertheless, a separate scaling factor for low $Wi$ looks to be irrelevant here. This is perhaps caused by the higher entanglement density of the C~1000~H~2002~ melt, which possibly delays any meaningful change in the relaxation time until approximately $Wi = 10$ when $\langle Z_{k}\rangle$ begins to decrease---see [Figure 5](#polymers-11-00476-f005){ref-type="fig"}a.$~\tau_{rot}$ also exhibits a power-law behavior that scales as ${\overset{˙}{\gamma}}^{- 0.7 \pm 0.07\ }$with flow strength. Although this value is slightly smaller than those of the C~400~H~802~ and C~700~H~1402~ melts (−0.78 and −0.75 respectively), they are all in reasonable agreement within statistical bounds. The ratio $\tau_{rot}/\tau_{d}$ averages about 7.3 over all $Wi \geq 50$, which is reasonably close to $2\pi$, similarly to the prior cases \[[@B2-polymers-11-00476],[@B3-polymers-11-00476]\]. This suggests that a single timescale, one associated with the period of the molecular tumbling cycles, is the sole configurational relaxation mechanism of the C~1000~H~2002~ chains for $\overset{˙}{\gamma} \geq \tau_{R}^{- 1}$.

### 3.2.2. Rheological Response {#sec3dot2dot2-polymers-11-00476}

[Figure 8](#polymers-11-00476-f008){ref-type="fig"} displays the steady-state rheological properties of the C~1000~H~2002~ liquid as functions of $Wi$. As expected, the shear stress scales as $\overset{˙}{\gamma}$ in the linear viscoelastic regime; however, at higher shear rates, the system's response is quite different from typical experimental observations, as evident from [Figure 8](#polymers-11-00476-f008){ref-type="fig"}a. Specifically, the shear stress passes through a maximum in the shear rate range $3 < Wi < 12$ and a subsequent minimum in the range $58 < Wi < 117$, in contradiction with the experimentally observed plateau region where the shear stress remains approximately constant or increases slightly as shear rate increases, usually within the shear rate ranges $\tau_{d}^{- 1} < \overset{˙}{\gamma} < \tau_{R}^{- 1}$ or $\tau_{d}^{- 1} < \overset{˙}{\gamma} < \tau_{e}^{- 1}$. Considering the uncertainties of the calculations, it appears that the local maximum and minimum in the shear stress profile occur roughly at about $\left. \overset{˙}{\gamma} \right.\sim\tau_{d}^{- 1}$ and $\left. \overset{˙}{\gamma} \right.\sim\tau_{R}^{- 1}$, respectively, and the shear stress surpasses the local maximum value at a shear rate of approximately $\tau_{e}^{- 1}$. This possibly implies that the flow is unstable over a fairly wide range of shear rates. Such behavior is enticingly consistent with the discussion of Doi and Edwards (see Figure 7.22 of Reference \[[@B34-polymers-11-00476]\]) concerning the DE model predictions at high shear rates, who argued that the power-law exponent of the shear stress is very sensitive to the relaxation spectrum of the linear relaxation modulus. They argued that the absolute value of the exponent becomes smaller (closer to zero) as the relaxation spectrum becomes broader. Therefore, the shear stress should be approximately independent of the shear rate for polydisperse samples that are commonly used in experiments (hence the plateau), whereas a maximum in the shear stress profile could result from a completely monodisperse sample. Nevertheless, even for monodisperse samples, multiple relaxation processes tend to broaden the relaxation spectrum and weaken the shear rate dependence of the stress. However, as evident from [Figure 7](#polymers-11-00476-f007){ref-type="fig"}, the number of timescales becomes effectively unity for $\overset{˙}{\gamma} \geq \tau_{R}^{- 1}$. Based on the DE model, $\left. \sigma_{xy} \right.\sim{\overset{˙}{\gamma}}^{- 0.5}$ for $\tau_{d}^{- 1} \ll \overset{˙}{\gamma} \ll \tau_{R}^{- 1}$ and as $\overset{˙}{\gamma}\tau_{R}^{- 1}$ becomes close to unity, the shear stress increases due to tube stretching \[[@B34-polymers-11-00476]\]. This implies that if the number of entanglements is not large enough (i.e., $\tau_{d}/\tau_{R}$ is not high enough), the shear rate dependence weakens. In [Figure 8](#polymers-11-00476-f008){ref-type="fig"}a, $\left. \sigma_{xy} \right.\sim{\overset{˙}{\gamma}}^{- 0.2}$ for $12 < Wi < 58$, which is consistent with this argument.

The plateau region in the shear stress profile has also been postulated to result from the onset of the molecular tumbling cycles that begin to manifest in this shear rate regime \[[@B2-polymers-11-00476],[@B3-polymers-11-00476]\]. This hypothesis led to further investigations which indicated the possibility that shear banding, caused by the molecular periodicity, was a possible cause of the experimentally observed plateau in the shear stress profile \[[@B36-polymers-11-00476],[@B37-polymers-11-00476],[@B38-polymers-11-00476]\]; however, it is unlikely that shear banding occurs in the present simulations since the p-SLLOD equations of motion impose a homogeneous linear velocity profile throughout the simulation cell in the NEMD simulations. That being said, however, recent DPD simulations have demonstrated shear banding in monodisperse polymers in the same range of molecular weight where the flow curve is non-monotonic \[[@B36-polymers-11-00476],[@B37-polymers-11-00476],[@B38-polymers-11-00476],[@B39-polymers-11-00476]\].

For $\overset{˙}{\gamma} > \tau_{e}^{- 1}$, the shear stress scales as ${\overset{˙}{\gamma}}^{0.3}$. The power-law exponents for the C~400~H~802~ and C~700~H~1402~ melts over the same range of shear rates are approximately $- 0.5$ and $- 0.4$, respectively \[[@B2-polymers-11-00476],[@B3-polymers-11-00476]\], which suggest a molecular weight dependence for the shear stress at these high shear rates. [Figure 8](#polymers-11-00476-f008){ref-type="fig"}b,c show the first normal stress coefficient $\mathsf{\Psi}_{1} = N_{1}/{\overset{˙}{\gamma}}^{2}$, and the second normal stress coefficient $\mathsf{\Psi}_{2} = N_{2}/{\overset{˙}{\gamma}}^{2}$, where $N_{1} = \sigma_{xx} - \sigma_{yy}$ and $N_{2} = \sigma_{yy} - \sigma_{zz}$. Both coefficients exhibit strong shear thinning behavior in the nonlinear regime with power-law exponents of --1.7 and --1.8, in agreement with those of the C~700~H~1402~ melt \[[@B2-polymers-11-00476]\]. The ratio $- \mathsf{\Psi}_{2}/\mathsf{\Psi}_{1}$ ranges over $0.04 < - \mathsf{\Psi}_{2}/\mathsf{\Psi}_{1} < 0.27$ in the nonlinear regime, again in reasonable agreement with C~700~H~1402~ melt \[[@B2-polymers-11-00476]\] and typical experimental values \[[@B40-polymers-11-00476]\].

3.3. Startup of Shear Flow Properties {#sec3dot3-polymers-11-00476}
-------------------------------------

### 3.3.1. Transient Behavior {#sec3dot3dot1-polymers-11-00476}

The time-dependent microstructural and rheological properties of the C~1000~H~2002~ melt were investigated under startup of simple shear flow similarly to those of the C~700~H~1402~ liquid presented in a prior publication \[[@B17-polymers-11-00476]\]. [Figure 9](#polymers-11-00476-f009){ref-type="fig"} displays the transient shear viscosity, $S_{xy}$, $\lambda$, and $\langle Z_{k}\rangle$ as functions of time for various $Wi$ obtained from the NEMD simulations. The data for the transient viscosity and $S_{xy}$ have been smoothed using a running time average over a number of successive sample times spanning 0.05--0.1 relaxation time at the corresponding $Wi$, as represented by the circles in [Figure 7](#polymers-11-00476-f007){ref-type="fig"}. It should be noted that $\lambda$ is very sensitive to the box shape when calculated using the Z1 code; since the box shape continuously changes during the simulation due to the Lagrangian rhomboid periodic boundary conditions, it is difficult to calculate transient tube stretch using the Z1 code. A solution to this problem is to calculate the tube stretch only at time steps when the box is rectangular or slightly (e.g., less than 5%) tilted. The tube stretch profiles displayed in [Figure 9](#polymers-11-00476-f009){ref-type="fig"} were obtained using this method. A major disadvantage of this method is that it significantly reduces the resolution of data, which could lead to the loss of important dynamical features, such as an overshoot or undershoot. However, unlike tube stretch, the entanglement density is not very sensitive to the simulation box shape; since the entanglement density has essentially similar dynamics as the primitive path contour length (and equivalently, the tube stretch---see [Figure 9](#polymers-11-00476-f009){ref-type="fig"}), it can be used to estimate the overshoot and undershoot times of the tube stretch. Note that an overshoot in tube stretch corresponds an undershoot in entanglement density, and vice-versa.

The transient first and second normal stress differences are shown in [Figure 10](#polymers-11-00476-f010){ref-type="fig"} for various $Wi$ as functions of time. This figure also displays the tube orientation tensor differences $S_{xx} - S_{yy}$ corresponding to $N_{1}$ and $S_{zz} - S_{yy}$ corresponding to $- N_{2}$ for comparison. Note that the data have been smoothed using the same method as discussed above. The data are displayed at three $Wi$ numbers, which were chosen to represent the three distinct nonlinear viscoelastic flow regimes: $\tau_{d}^{- 1} \leq \overset{˙}{\gamma} < \tau_{R}^{- 1}$, $\tau_{R}^{- 1} \leq \overset{˙}{\gamma} < \tau_{e}^{- 1}$, and $\overset{˙}{\gamma} \geq \tau_{e}^{- 1}$.

It is evident from [Figure 9](#polymers-11-00476-f009){ref-type="fig"} and [Figure 10](#polymers-11-00476-f010){ref-type="fig"} that the transient viscosity and normal stresses are in qualitative agreement with typical experimental data. Specifically, except for $N_{2}$ at $Wi = 12$, they all exhibit an overshoot for $Wi \geq 12$ before they attain steady state. Additionally, the overshoot in shear viscosity is followed by an undershoot at least for $Wi \geq \tau_{R}^{- 1}$, again in agreement with typical experiments. These overshoots and undershoots (if any) also occur in the entanglement network variables, as shown in [Figure 9](#polymers-11-00476-f009){ref-type="fig"} and [Figure 10](#polymers-11-00476-f010){ref-type="fig"}. These figures make it possible to compare the dynamics of the stress tensor with those of tube variables (i.e., the tube segmental orientation tensor $\mathbf{S}$, and the tube stretch $\lambda$) to investigate the origins of these phenomena, as discussed in the next section. It is worth mentioning that steady or transient shear banding might occur in the range $\tau_{d}^{- 1} \leq \overset{˙}{\gamma} < \tau_{e}^{- 1}$. This phenomenon cannot be investigated here due to the use of p-SLLOD equations of motion, as discussed in [Section 3.2.2](#sec3dot2dot2-polymers-11-00476){ref-type="sec"}. As a consequence, the quantities presented in this section could be affected, assuming shear banding occurs. However, we do not expect a significant change, especially in ensemble-averaged quantities such as the stress tensor and tube variables. Cao and Likhtman \[[@B41-polymers-11-00476]\] compared the startup shear behavior of entangled melts obtained from NEMD simulations using the SLLOD equations and a Langevin thermostat with those of boundary-driven DPD simulations. These comparisons suggested that the ensemble average shear stresses obtained from these two methods were consistent (although not identical) despite the presence of shear banding at the examined shear rates.

### 3.3.2. Stress Overshoot and Undershoot {#sec3dot3dot2-polymers-11-00476}

From [Figure 9](#polymers-11-00476-f009){ref-type="fig"}, it appears that the dynamic response of shear viscosity and $S_{xy}$ are roughly synchronized over a wide range of $Wi$. Specifically, the overshoot and undershoot of transient viscosity (if any) occur approximately at the same time as those of the $S_{xy}$ component of the tube orientation tensor. On the contrary, tube stretch and entanglement density respond to the applied flow field with a notable lag as compared to $\eta^{+}$ and $S_{xy}$. It is worth noting that the displayed $Wi$ numbers represent various flow regimes: $Wi = 12$ lies in $Wi_{R} < 1$ regime where tube stretch is negligible; $Wi = 58$ is within the regime where tube stretch is significant, and $Wi = 1170$ is within the regime where molecular tumbling is dominant. It should also be noted that this classification is based on the steady-state responses and that it might not necessarily remain valid in transient situations. For instance, whereas the tube stretch is minor at $Wi_{R} < 1$, it could exhibit an overshoot in transient situations. Although the magnitude of the shear viscosity (and stress) is a function of both tube orientation, $S_{xy}$, and stretch, $\lambda$, these plots suggest that the dynamics of shear viscosity are mainly influenced by the tube segment orientation, $S_{xy}$, which indicates that the principal origin of stress overshoot and undershoot is possibly tube segmental orientation. These plots also show that there is no significant undershoot in $\lambda$ (or equivalently, an overshoot in $\langle Z_{k}\rangle$) at any $Wi$. This observation, that also applies to other shear rates (not shown in [Figure 9](#polymers-11-00476-f009){ref-type="fig"}), practically rules out tube stretch as the origin of the stress undershoot at high shear rates. It is also evident from [Figure 10](#polymers-11-00476-f010){ref-type="fig"} that the dynamics of $N_{1}$ and $N_{2}$ are also in good agreement with their corresponding components of the tube orientation tensor, i.e., $S_{xx} - S_{yy}$ and $S_{yy} - S_{zz}$, respectively, suggesting that the overshoot in normal stresses arises from the tube segment orientation.

[Figure 11](#polymers-11-00476-f011){ref-type="fig"} shows the overshoot (panel (a)) and undershoot (panel (b)) times for the transient viscosity and $S_{xy}$ component of the tube orientation tensor as functions of $Wi$. It also displays the undershoot time for the entanglement density in both panels for comparison. Note that an undershoot in $\langle Z_{k}\rangle$ corresponds to an overshoot in $\lambda$, as discussed before. It is evident that the transient viscosity overshoot and undershoot times effectively overlap with those of $S_{xy}$ at all $Wi < 585$. At higher $Wi$, although these two curves look to be diverging, the difference between the two times is not significant, considering the error associated with extracting these small values from the noisy data, as shown in [Figure 9](#polymers-11-00476-f009){ref-type="fig"}. On the other hand, it is obvious that there is a significant difference between the undershoot time in $\langle Z_{k}\rangle$ and either the overshoot and undershoot times of shear viscosity. These results again imply that both the stress overshoot and undershoot are originated from similar phenomena in the tube segmental orientation. This conclusion is in agreement with observations for a C~700~H~1402~ melt at high shear rates \[[@B17-polymers-11-00476]\]. It also agrees with the results of Cao and Likhtman \[[@B42-polymers-11-00476]\] for unentangled and mildly entangled systems, indicating that the origin of the stress overshoot at low shear rates is the orientation of the tube network rather than chain stretching. Jeong, et al. \[[@B9-polymers-11-00476]\] also attributed the stress overshoot to the segmental orientation in a wide range of flow strength for a mildly entangled C~400~H~802~ polyethylene melt. However, unlike the current results, they did not observe a clear overshoot in the primitive path of the contour length (and hence in the tube stretch) even at very strong flow fields. This may be due to the relatively low entanglement density of the C~400~H~802~ molecules used in their NEMD simulations. Masubuchi et al. \[[@B43-polymers-11-00476]\] also investigated the origin of the stress undershoot at high shear rates using primitive chain network simulations. They examined segmental orientation, tube stretch, and the ensemble average squared sine of the chain end-to-end orientation angle (representing the tumbling motion) and showed that all these variables exhibit undershoots, although not synchronized with shear stress. Masubuchi et al. \[[@B43-polymers-11-00476]\] concluded that their results supported the mechanism proposed by Costanzo et al. \[[@B16-polymers-11-00476]\]

[Figure 12](#polymers-11-00476-f012){ref-type="fig"} displays the shear strain, $\lambda$, at the times of overshoot (panel a) and undershoot (panel b) in the transient viscosity and $S_{xy}$ component of the tube orientation tensor, as functions of $Wi$. Shear strain expresses the theoretical deformation due to the applied flow field, which is calculated as $\gamma = \overset{˙}{\gamma}t$. This figure also shows the strain at the undershoot time for the entanglement density in both panels for comparison. The agreement between $\eta^{+}$ and $S_{xy}$ curves in the region $Wi < 585$ is not surprising, considering the results of [Figure 11](#polymers-11-00476-f011){ref-type="fig"} and how the strain is calculated. The important point to notice is that up to very high$~Wi$ the overshoot in $S_{xy}$ occurs at about $\gamma = 2$ consistently. This suggests that regardless of flow strength, the material deforms affinely during the initial 2 strain units until $S_{xy}$ attains a maximum. However, it does not look to be the case at later times that $S_{xy}$ passes beyond its minimum, especially at intermediate and high $Wi$ as evident from [Figure 12](#polymers-11-00476-f012){ref-type="fig"}b. It is also worth noting that the strain at the maximum of viscosity is roughly 2 at low and intermediate $Wi$, in agreement with the prediction of the Doi-Edwards model \[[@B44-polymers-11-00476]\], and shifts to higher values for $Wi \geq 585$. The experimental value of strain at stress overshoot is also about 2 at low shear rates; however, it shifts to higher strains, as the shear rate exceeds $\tau_{R}^{- 1}$ \[[@B45-polymers-11-00476]\].

The discussion concerning the overshoot and undershoot dynamics in the last few paragraphs should not lead to misinterpretation about the role of tube stretch in the stress overshoot and undershoot. [Figure 13](#polymers-11-00476-f013){ref-type="fig"}a shows the magnitude of the overshoots in the normalized shear stress and tube orientation $S_{xy}^{os}$, versus $Wi$. The shear stress is normalized with the plateau modulus, $G_{N}^{0}$. [Figure 13](#polymers-11-00476-f013){ref-type="fig"} also displays the magnitude of the tube stretch at the time of stress overshoot. Note that this quantity is different from the magnitude of the tube stretch overshoot. It is evident from this figure that for $Wi < 585$ the shear stress closely mimics $S_{xy}^{os}$, while the tube stretch is fairly close to the equilibrium value of unity, or only mildly greater. This suggests that in this region tube stretch has a minor or negligible contribution to the stress overshoot, $\sigma_{xy}^{os}$. At higher shear rates, whereas $S_{xy}^{os}$ looks to become saturated and remain roughly constant, $\sigma_{xy}^{os}$ increases quickly as $Wi$ increases. The tube stretch in this region also begins to increase and diverge from its equilibrium value. This indicates that although the dynamics of the stress overshoot are essentially controlled by the tube segmental orientation (as discussed in the preceding paragraph), its magnitude is significantly influenced by the tube stretch at high flow strength. A similar argument can be made for the stress undershoot---see [Figure 13](#polymers-11-00476-f013){ref-type="fig"}b. This conclusion can be rationalized by hypothesizing that the tube stretch is itself originated from the tube orientation or another dynamic variable. It is however immediately evident from [Figure 9](#polymers-11-00476-f009){ref-type="fig"} that the $S_{xy}$ component could not be that variable, considering the significant differences between the features of $S_{xy}$ and $\lambda$ plots in this figure.

[Figure 14](#polymers-11-00476-f014){ref-type="fig"} shows the undershoot time for the entanglement density as well as the overshoot time for the ensemble average squared end-to-end distance, $\langle R^{2}\rangle$, and the normal (diagonal) components, $S_{xx}$, $S_{yy}$, and $S_{zz}$, of the tube orientation tensor. Overall, this figure shows that the overshoot times of these variables roughly overlap, within a wide range of $Wi$ including mildly to strongly nonlinear viscoelastic flow regimes. Specifically, there is a good agreement between the undershoot time of $\langle Z_{k}\rangle$ and that of $\langle R^{2}\rangle$. It should be emphasized that $\langle R^{2}\rangle$ is essentially the trace of the ensemble average chain conformation tensor and represents the overall average extensional state of the molecules. [Figure 14](#polymers-11-00476-f014){ref-type="fig"} implies that the entanglement network, and hence tube stretch dynamics, are mainly influenced by the diagonal components of the orientation tensor or by the overall extensional properties of the molecules rather than the shear component.

4. Conclusions {#sec4-polymers-11-00476}
==============

Transient and steady-state dynamic responses of an entangled C~1000~H~2002~ polyethylene melt were examined via virtual experimentation using NEMD simulations. Under quiescent conditions, reptation theory could explain equilibrium properties fairly well. Under steady shear flow conditions, four flow regimes were recognized in agreement with prior results for moderately and mildly entangled C~700~H~1402~ and C~400~H~802~ liquids \[[@B2-polymers-11-00476],[@B3-polymers-11-00476]\]. The first regime was the linear viscoelastic regime ($\overset{˙}{\gamma} < \tau_{d}^{- 1}$) where most of the structural and topological properties of the system remain unperturbed compared to the quiescent conditions. Orientation effects dominated the rheological response in this flow regime, although they are quite weak. In the second regime ($\tau_{d}^{- 1} < \overset{˙}{\gamma} < \tau_{R}^{- 1}$), the molecules began to align with the flow direction and a significant degree of chain orientation was observed as $Wi$ increased. Additionally, the tube segments began to stretch mildly and chain molecules partially unraveled and disentagled as flow strength increased. However, the dominant relaxation mechanism in this region was the orientation of the tube segments. In the third regime ($\tau_{R}^{- 1} < \overset{˙}{\gamma} < \tau_{e}^{- 1}$), while on average the chains were fully aligned with the flow direction, the molecular disentangling continued and tube stretching dominated the rheological response. Additionally, the rotation of molecules became a significant source of the overall system dynamics. In the fourth regime ($\overset{˙}{\gamma} > \tau_{e}^{- 1}$), the chain stretching decelerated, and tube stretch approached a plateau value. At the same time, flow-induced disentanglement continued and the entanglement network began to deteriorate such that some molecules became completely devoid of entanglements. The molecular tumbling, on the other hand, gradually became the dominant relaxation mechanism, and molecular configurations followed more regular cycles when compared to similar behavior at lower flow strength.

The comparison of transient shear viscosity, $\eta^{+}$, with the dynamic responses of key variables of the tube model, including the tube segmental orientation, $S$, and tube stretch, $\lambda$, revealed that the stress overshoot and undershoot in steady shear flow of entangled liquids were essentially originated and dynamically controlled by the $S_{xy}$ component of the tube orientation tensor, rather than the tube stretch $\lambda$, over a wide range of flow strengths (including shear rates faster than $\tau_{R}^{- 1}$). Nevertheless, the magnitude of the stress is significantly affected by $\lambda$ at high shear rates.

Computational resources for this project were provided by the NSF under CBET-0742679 to the PolyHub Engineering Virtual Organization as well as allocation of advanced computational resources by the National Institute for Computational Sciences and the ORNL Joint Institute for Computational Sciences. Financial support was provided by NSF under CBET-1602890. This work also used the Extreme Science and Engineering Discovery Environment, which is supported by NSF (ACI-1548562) using Bridges (TG-CTS150054).

Conceptualization, B.E.; Formal analysis, M.H.N.S.; Funding acquisition, B.E.; Methodology, M.H.N.S.; Project administration, B.E. and B.K.; Resources, B.K.; Supervision, B.E. and B.K.; Writing -- original draft, M.H.N.S. and B.E.

The authors declare no conflicts of interest.

![Short (**a**) and long (**b**) timescale segmental MSD versus time of the (500) centermost chain atomic units.](polymers-11-00476-g001){#polymers-11-00476-f001}

![Characteristic timescales as functions of molecular length. Power-law exponents agree well with experimental observations.](polymers-11-00476-g002){#polymers-11-00476-f002}

![Probability distribution of the normalized magnitude of the end-to-end vector $\left| R \right|/\left| R \right|_{\max}$ at various values of $Wi$ (**a**), the mean-square of the end-to-end vector and 6 times the radius of gyration vs. $Wi$ (**b**), ensemble average chain orientation angle with respect to the flow direction vs $Wi$ (**c**), and tube stretch vs $Wi$ (**d**).](polymers-11-00476-g003){#polymers-11-00476-f003}

![Normalized end-to-end vector and orientation angle vs. time for a random chain at $Wi = 58$ (**a**) and $Wi = 1170$ (**b**). A time progression of configurations of a random chain at $Wi = 117$ is depicted in panel (**c**). $\theta_{ete}$ is the angle between the chain end-to-end vector and the flow direction ($x$) and $\theta_{pa}$ is the angle between the primary axis of the chain and the flow direction.](polymers-11-00476-g004){#polymers-11-00476-f004}

![Average number of entanglements per chain versus $Wi$ (**a**), the effective tube diameter versus $Wi$ (**b**), and probability distributions of the number of entaglements per chain (**c**) and the primitive path contour length (**d**) at various values of $Wi$.](polymers-11-00476-g005){#polymers-11-00476-f005}

![The $S_{xy}$ (**a**) and diagonal (**b**) components of the tube segmental orientation tensor $\mathbf{S}$, as functions of $Wi$ at steady state, as obtained from the NEMD data for the C~1000~H~2002~ melt.](polymers-11-00476-g006){#polymers-11-00476-f006}

![The three characteristic times ($\tau_{d}$, $\tau_{rot}$, and $\tau_{r}$) of the C~1000~H~2002~ melt versus $Wi$.](polymers-11-00476-g007){#polymers-11-00476-f007}

![Shear stress (**a**), first and second normal stress coefficients (**b**), (**c**), and the ratio of the normal stress coefficients (**d**), all as functions of $Wi$.](polymers-11-00476-g008){#polymers-11-00476-f008}

![Shear viscosity (blue), the $S_{xy}$ component of the tube segmental orientation tensor (red), tube stretch *λ* (purple), and normalized entanglement density (green) versus time upon startup of shear flow at $Wi = 12$ (**a**), $Wi = 58$ (**b**), and $Wi = 1170$ (**c**). The dynamics of $\langle Z_{k}\rangle$ are similar to those of $\lambda$ except that the minimum in $\langle Z_{k}\rangle$ corresponds to a maximum in $\lambda$. Note that in panel (**c**) there appear small gaps in some of the data profiles at long times where simulation data was accidentally deleted. Since these data points had no bearing on the present discussion, the simulations were not repeated.](polymers-11-00476-g009){#polymers-11-00476-f009}

![Transient first (left panels) and second (right panels) normal stress differences as well as their corresponding tube orientation tensor differences as functions of time upon startup of shear flow. Normal stress differences are normalized with respect to the plateau modulus. Weissenberg numbers are 12, 58, and 1170 from top to bottom rows, respectively.](polymers-11-00476-g010){#polymers-11-00476-f010}

![Overshoot (**a**) and undershoot (**b**) times for the transient viscosity, $\eta^{+}$, and $S_{xy}$ component of the tube orientation tensor, as well as the undershoot time for the entanglement density as a function of $Wi$ in both panels. Note that the undershoot time in $\langle Z_{k}\rangle$ (which corresponds to the overshoot time in tube stretch) does not coincide with either the overshoot or undershoot time in $\eta^{+}$ or $S_{xy}$.](polymers-11-00476-g011){#polymers-11-00476-f011}

![Shear strain, $\gamma$, at the times of overshoot (**a**) and undershoot (**b**) in the transient viscosity, $\eta^{+}$, and $S_{xy}$ component of the tube orientation tensor, as well as strain at the undershoot time for the entanglement density as a function of $Wi$ in both panels. Note that the undershoot time in $\langle Z_{k}\rangle$ does not coincide with either the overshoot or undershoot strain in $\eta^{+}$ or $S_{xy}$. The strain overshoots for $\eta^{+}$ and $\langle Z_{k}\rangle$ scale as $Wi^{0.17}$ and $Wi^{0.33}$ respectively, for $Wi > 58$.](polymers-11-00476-g012){#polymers-11-00476-f012}

![Magnitudes of the normalized shear stress, $\sigma_{xy}/G_{N}^{0}$, the $S_{xy}$ component of the tube orientation tensor, and the tube stretch $\lambda$ at the stress overshoot (**a**) and undershoot (**b**) times versus $Wi$. The stress overshoot scales as $Wi^{0.3}$ in the range $3.5 \leq Wi \leq 2340$, and the tube stretch at the stress overshoot time scales as $Wi^{0.1}$ for $Wi > 58$.](polymers-11-00476-g013){#polymers-11-00476-f013}

![Comparison of the undershoot time for the entanglement density with the overshoot times for the ensemble average squared end-to-end distance, $\langle R^{2}\rangle$, and the normal (diagonal) components, $S_{xx}$, $S_{yy}$, and $S_{zz}$, of the tube orientation tensor.](polymers-11-00476-g014){#polymers-11-00476-f014}

polymers-11-00476-t001_Table 1

###### 

Details of the simulation cells: $L_{x}$, $L_{y}$, and $L_{z}$ are box lengths in the $x$, $y$, and $z\ $dimensions, respectively.

  $\mathbf{Wi}$   $\mathbf{L}_{\mathbf{x}}\left( \mathring{\mathrm{A}} \right)$   $\mathbf{L}_{\mathbf{y}}~\mathbf{and}~\mathbf{L}_{\mathbf{z}}\left( \mathring{\mathrm{A}} \right)$   Number of Particles
  --------------- --------------------------------------------------------------- ---------------------------------------------------------------------------------------------------- ---------------------
  0--1.2          84.7                                                            84.7                                                                                                 20,000
  3.5--12         169.3                                                           84.7                                                                                                 40,000
  41--58          254.0                                                           84.7                                                                                                 60,000
  117--11,700     508.0                                                           84.7                                                                                                 120,000

polymers-11-00476-t002_Table 2

###### 

The three relaxation times for the C~1000~H~2002~ melt calculated according to reptation theory and segmental mean-square displacement (MSD) data.

  Relaxation Time   Theory   MSD
  ----------------- -------- ------
  $\tau_{e}$ (ns)   6.4      2.7
  $\tau_{R}$ (ns)   137      194
  $\tau_{d}$ (ns)   5305     5834
