Virtual experimentation of atomistic entangled polyethylene melts undergoing planar elongational flow revealed an amazingly detailed depiction of individual macromolecular dynamics and the resulting effect on bistable configurational states. A clear coil-stretch transition was evident, in much the same form as first envisioned by de Gennes for dilute solutions of high polymers, resulting in an associated hysteresis in the configurational flow profile over the range of strain rates predicted by theory. Simulations conducted at steady state revealed bimodal distribution functions, in which equilibrium configurational states were simultaneously populated by relatively coiled and stretched molecules which could transition from one conformational mode to the other over a relatively long time scale at critical values of strain rates. The implication of such behavior points to a double-well conformational free energy potential with an activation barrier between the two configurational minima.

Understanding the behavior and dynamics of polymeric liquids under flow has been a daunting challenge for well over a hundred years. Elongational flows, in particular, have proven to be especially difficult to model since the imposed velocity gradient highly stretches the flexible polymer macromolecules in the direction of flow while compressing them in the orthogonal directions. These stretched chain configurations are far from the equilibrium random coil conformations assumed by the individual macromolecules and often lead to non-Newtonian properties being exhibited by these liquids. The prevailing view over much of the history of polymer fluid dynamics was that the extent of individual chain stretching could be characterized at steady state by a unimodal distribution function whose peak narrowed and shifted to higher degrees of stretching as the elongation rate increased—see, for general reference, Fig. 3. This generic behavior was assumed to hold true for all types of flexible macromolecular liquids, ranging from unentangled dilute solutions to highly entangled concentrated solutions and melts.

The prevailing view began to change in 1974 with the publication of the seminal work of de Gennes^{1} on what soon thereafter became known as the “coil-stretch transition.” de Gennes used the elastic dumbbell model of kinetic theory to demonstrate that dilute solutions undergoing planar elongational flow could exhibit multiple steady-state configurations at strain rates in the vicinity of $\tau R\u22121$ which is the reciprocal of the Rouse time scale that governs the stretching of chain conformation. Two of these states are physically stable (i.e., a bistable equilibrium), corresponding to a relatively coiled (C) configuration and a relatively stretched (S) state, with a possible third state representing an unstable (P) condition.^{1} Whether or not multiple solution states appear is believed to depend on the chain length of the macromolecules comprising the dilute solution. This general behavior is illustrated in Fig. 1(a) on a plot of extension, *ℓ*, vs. elongation rate, $\epsilon \u0307$, for three different molecular weight linear polymers, denoted and ranked according to MW1 < MW2 < MW3. For shorter chains, MW1 and MW2, fractional extension of the macromolecule increases monotonically with increasing $\u03f5\u0307$, with the steepness of the increase becoming much more dramatic in the neighborhood of $\epsilon \u0307\u2248\tau R\u22121$. In the case of the longest macromolecule, three distinct values of chain extension are possible, as denoted by the points C, P, and S at $\epsilon \u0307ref\u2248\tau R\u22121$ in Fig. 1(a). An interesting consequence of the S-shaped profile of MW3 is the possibility that the configurational state of the system is determined by the imposed kinematics: i.e., increasing $\epsilon \u0307$ up to $\epsilon \u0307ref$ from below would most likely result in the coiled state, C, whereas decreasing $\epsilon \u0307$ down to $\epsilon \u0307ref$ from above would result in the stretched state S. This implies the possibility of the existence of a hysteresis loop in the macromolecular configuration of the elongational flow profile; such a possibility was noted by both de Gennes^{1} and Tanner.^{2}

To make sense of the physical behavior described above, it is instructive to consider the free energy profile of macromolecular chain configurations at $\epsilon \u0307ref$. Free energy is characterized vs. extension in Fig. 1(b) for three different chain lengths at the strain rate corresponding to the vertical line in Fig. 1(a). For shorter chains, there is only a single global minimum in the free energy which shifts to higher values of extension as either the chain length or the strain rate increases. For the longest chain, however, a second minimum appears at a much higher value of extension, which is significantly separated from the extrapolation of the original coiled state. Each minimum, as well as the maximum point separating them, is a valid state point for the system, although the maximum is unstable. Hence the configurational state of the macromolecular solution will be determined, in principle, at least at steady state, by the minimum point with the lower value of free energy (i.e., the global minimum). Which state, C or S, is the global minimum depends on the value of the strain rate, $\epsilon \u0307$, with the point P representing the activation barrier between the two possible configurational states. The critical value of $\u03f5\u0307$ where the two minima have the same depth, $\epsilon \u0307*$, is typically called the “coil-stretch transition” point, and it represents a value of the strain rate at which, even at infinite time, there is no preference for one state over the other. The critical extension rate, expressed in terms of the dimensionless Deborah number, $De=\epsilon \u0307\tau $, was later found by Larson and Magda^{3} to be 0.5, where *τ* is the longest relaxation time of the liquid (usually the Rouse time for dilute solutions). For $\epsilon \u0307\u2260\epsilon \u0307*$, either the C or the S state is preferred, and although the macromolecular configuration may be trapped within the local minimum for a period of time, eventually the chain configuration should find its way over the activation barrier to the global minimum at steady state. Unfortunately, this temporary trapping can potentially lead to ambiguous experimental observations since it can be very demanding to achieve a steady-state condition under elongational flow.

In spite of the initial theoretical evidence suggesting the coil-stretch transition, a long debate ensued regarding the existence of the S-shaped profile, which is evident in Fig. 1(a), and its implications. Fan *et al.*^{4} argued against the S-shaped profile on theoretical grounds, even going as far as to present a mathematical proof of its impossibility. Further computer simulations of bead-spring chains by Wiest *et al.*^{5} revealed a completely reversible monotonic flow profile with no signs of hysteresis. It should be noted, however, that these authors did not examine the impact of macromolecular free energy on the configurational flow profile; i.e., the authors’ choice of diffusion equation (itself merely an approximation of real fluid behavior) was only consistent with a single-well free energy potential function. Furthermore, even if their analysis were consistent with a double-well potential function, their solution would evolve to the infinite-time global minimum at the steady-state limit; consequently, their arguments overlook the possibility that the macromolecular conformation could become trapped in a distinct configuration other than the one corresponding to the global free energy minimum over the (insufficiently long) duration of the experiment.

With equally credible evidence suggesting both the existence and non-existence of the coil-stretch transition, the debate continued until experimental instrumentation evolved to the point where direct visualization of macromolecular configurations could be implemented under elongational flow. Perkins *et al.*^{6,7} directly visualized individual fluorescently labeled DNA molecules in dilute solution subjected to planar elongational flow and observed that different but identical molecules deformed under flow over greatly varying time scales, which were effectively governed by the specific configuration of the random conformation of the chain at the instant of flow inception. The implication of these observations was that any apparent dramatic increase in chain extension at a critical value of the strain rate was primarily a manifestation of steady-state behavior in the infinite-time limit (i.e., after all the chains had eventually extended) and was not an indication that each individual chain had suddenly experienced an abrupt unravelling. When only those molecules that had attained a steady-state extension during the duration of the experiment were considered, it was possible to identify a hypothetical coil-stretch transition near the predicted value of *De* = 0.5.

Schroeder *et al.*^{8} also performed direct visualization experiments on fluorescently dyed DNA molecules with contour lengths of 575 *μ*m and 1300 *μ*m in dilute solution under planar extensional flow over a wide range of strain rates. At each applied strain rate (quantified in terms of *De*), two distinct experiments were conducted to determine if hysteresis was evident in the macromolecular configurational response: one in which the strain rate was imposed from the quiescent state and another in which the strain rate was reduced to the equivalent value of the first experiment after pre-straining the material at a very high rate (*De* ≥ 5). For the case of shorter molecules, no hysteresis was evident as the initially extended molecules gradually recoiled to the same extension value as the molecules which had been slowly stretched from the quiescent condition. In the case of longer molecules, however, clear signs of a bistable steady state were evident in the range of *De* ∈ (0.3, 0.6), implying a hysteresis behavior very similar to that predicted by de Gennes three decades earlier. Furthermore, the authors performed Brownian dynamics simulations of a bead-spring chain in dilute solution and observed a similar hysteresis effect and also determined the conformational free energy landscape of the flow system, which exhibited a double-well potential similar to that shown in Fig. 1(b) within the same range of *De* in which bistable steady states had been observed in the experiments.

Dilute solutions are just one type of polymeric liquid and are not necessarily the most ubiquitous. Industrial processing of plastics and fibers is primarily concerned with concentrated solutions and melts, both of which are composed of densely packed and highly entangled macromolecules. To date, most theoretical and experimental studies of the coil-stretch transition have examined dilute solutions simply because they are the easiest to analyze theoretically and to visualize experimentally and primarily because the individual macromolecules do not interact with each other. Indeed, for entangled concentrated solutions or melts, there exists no conclusive evidence to date, either theoretically or experimentally, of a coil-stretch transition in polymer melts or concentrated solutions. From an experimental perspective, visualization of individual chains in a dense liquid is fraught with difficulties that must be overcome in order to collect meaningful statistical data, including the fact that the individual polymer chains mutually interact and thereby greatly influence each other’s dynamical behavior.

In this study, we examined the dynamics of individual molecules comprising entangled polyethylene melts undergoing planar elongational flow via “virtual experimentation,” i.e., using equilibrium and nonequilibrium molecular dynamics (NEMD) simulations of realistic atomistic chains at actual experimental state points. NEMD simulations of monodisperse linear C_{700}H_{1402} and C_{1000}H_{2002} melts were performed in the NVT statistical ensemble at 450 K and the experimental density of 0.766 g/cm^{3} using the Siepmann-Karaboni-Smit (SKS) united-atom potential model,^{9} as discussed by Nafar Sefiddashti *et al.*^{10} This model has been used in many studies of alkane and polyethylene fluid behavior with outstanding results; indeed, the original article^{9} has been cited over 475 times to date. The NEMD equations of motion^{11–13} with a Nosé-Hoover thermostat^{14,15} were integrated using the reversible reference system propagator algorithm^{16} subject to Kraynik and Reinelt periodic boundary conditions, which allowed for unrestricted simulation time of the deforming simulation box^{17} undergoing planar elongational flow over the full range of imposed strain rates. Further details of the simulation method can be found in the supplementary material.

Table I presents some of the characteristic quiescent properties of the polyethylene melts examined herein. Equilibrium properties of the C_{700}H_{1402} melt were comprehensively described in previous work;^{18} the same methods were used herein to calculate the corresponding properties of C_{1000}H_{2002} which are presented in Table I. The disengagement time, *τ*_{d}, is typically the longest characteristic time for polymer melts. It represents the time scale over which a macromolecule will diffuse from a tube formed by neighboring entangled chains. *τ*_{R} is the Rouse time, ⟨*R*^{2}⟩^{1/2} and $Rg21/2$ are the root mean squares of the molecular end-to-end distance and radius of gyration, respectively, ⟨*Z*⟩ is the ensemble average number of entanglements (per chain) as calculated from the Z1 code developed by Kröger,^{19} and *L* is the theoretical maximum chain length. The prevailing view prescribes the disengagement time as the characteristic time scale relevant to chain orientation, whereas the Rouse time scale is characteristic of macromolecular chain stretching.^{20} This study is primarily concerned with the stretching response of the molecules, and hence, the dimensionless extension rate (i.e., the Deborah number) is characterized in terms of the Rouse time scale as $De=\epsilon \u0307\tau R$, similar to the case of dilute solutions.

System . | τ_{d} (ns)
. | τ_{R} (ns)
. | ⟨R^{2}⟩^{1/2} (Å)
. | $Rg21/2$ (Å) . | ⟨Z⟩
. | L (Å)
. |
---|---|---|---|---|---|---|

C_{700}H_{1402} | 1216 | 66 | 123.5 | 49.8 | 8.6 | 902.8 |

C_{1000}H_{2002} | 5270 | 137 | 141.8 | 56.41 | 12.8 | 1290.2 |

System . | τ_{d} (ns)
. | τ_{R} (ns)
. | ⟨R^{2}⟩^{1/2} (Å)
. | $Rg21/2$ (Å) . | ⟨Z⟩
. | L (Å)
. |
---|---|---|---|---|---|---|

C_{700}H_{1402} | 1216 | 66 | 123.5 | 49.8 | 8.6 | 902.8 |

C_{1000}H_{2002} | 5270 | 137 | 141.8 | 56.41 | 12.8 | 1290.2 |

Figure 2 displays the fractional chain extension, $x=R21/2L$, as a function of the Hencky strain, $\epsilon H=\epsilon \u0307t$, at six values of *De* for each macromolecule in the C_{700}H_{1402} melt. The Hencky strain provides a quantitative measure of the relative deformation of the material; i.e., the relative separation of two fluid elements increases as exp(*ε*_{H}). Obtaining a Hencky strain of 10 in an actual experiment can be quite a daunting challenge; this is one advantage of virtual experimentation, where simulations can be run indefinitely within the limits of the available resources. Colored lines in Fig. 2 display the fractional extension of each individual molecule and the black line represents the ensemble average. At the lowest value of *De*, the random coils are only mildly extended under flow with a rather wide distribution—see the corresponding distribution plots in Fig. 3. As *De* increases to intermediate values, it becomes increasingly evident that a very broad distribution of chain extension is observed and that the individual chains experience very different dynamics; some of the chains stretch to their steady-state extensions very quickly, whereas others respond on a time scale that is a full order of magnitude slower. This behavior is entirely consistent with the molecular individuality observed in the visualization experiments of planar extensional flows of dilute polymer solutions.^{6–8}

Despite the different flow dynamics of individual chains undergoing transient elongational flow, their steady-state behavior looks fairly normal. Figure 3 presents the probability distribution functions of fractional chain extension at steady state for the same values of *De* displayed in Fig. 2. For this purpose, the steady state begins at the instant (or equivalently the Hencky strain) that the time derivative of the ensemble average chain extension becomes negligible, as marked by the vertical lines and arrows in Fig. 2; therefore, the distribution functions of Fig. 3 are calculated based on the data to the right of the vertical lines. These plots resemble Gaussian distributions with a single peak that shifts to higher values of chain stretch, *x*, as *De* increases. In general, the distributions are wider at low *De*, with a standard deviation close to that of the quiescent condition, and become narrower with increasing *De*. This behavior essentially matches common expectations based on intuition and mean-field theory. However, at *De* = 0.16 and *De* = 0.54, the distributions are wider on the left side than would normally be expected, with *De* = 0.54 even exhibiting a slight bimodal character, as shown in the inset. Such a distribution could possibly indicate a coil-stretch transition between a bistable steady state with an associated hysteresis behavior over a narrow range of *De*, although this conjecture cannot be supported conclusively based on the evidence presented thus far. Hence a plot of steady-state extension vs. strain rate would resemble the MW2 curve in Fig. 1(a), with no discernable signs of the S-shaped profile predicted by de Gennes.

A similar analysis was performed for the C_{1000}H_{2002} melt as shown in Figs. 4 and 5, which is composed of molecules significantly longer than those of the C_{700}H_{1402} liquid. Figure 4 displays the individual and ensemble average chain fractional extensions, and Fig. 5 presents the steady-state probability distribution functions of fractional extension for six values of *De*. It is immediately evident that the extensional behavior at low (*De* = 0.06) and high (*De* ≥ 3.0) flow strength is qualitatively similar to that of the C_{700}H_{1402} melt. Specifically, different molecules exhibit very different transient flow dynamics in response to the flow, with some molecules requiring up to 10 Hencky strain units to achieve a steady state. The probability distributions are Gaussian-like, unimodal, and rather wide at low *De* = 0.06 and have shifted to the right and narrowed at high extension rates (*De* ≥ 3.0), as all would be expected from the conventional wisdom. On the contrary, the dynamics are very peculiar within the intermediate range of flow strength, 0.3 ≤ *De* ≤ 1.5. Specifically, plots of individual chain extension versus strain in this region show a very diverse response for different chains. Whereas some chains significantly stretch only after 2 strain units, others remain tightly coiled virtually permanently. Furthermore, some of the chains actually cycle back and forth between highly extended and tightly coiled configuration states at time scales that seem only mildly correlated with the flow strength. This behavior is remarkably similar to what would be expected from a bistable steady-state configuration produced by a double-well potential, such as that displayed in Fig. 1(b), with a low activation barrier between the two states. In this type of scenario, minor conformational changes induced by the flow could possibly induce a transition over the energy barrier to different configurational state.

At the extremes of the strain rate regime examined, the distribution of fractional extension is very similar to that of the C_{700}H_{1402} liquid; i.e., a unimodal distribution that narrows and shifts to higher values of extension as *De* increases. Within the intermediate range of 0.3 ≤ *De* ≤ 1.5, however, an extremely wide bimodal distribution is observed, with peaks corresponding to relatively coiled and stretched states. The peak at lower values of extension is fairly close to the average extension of the chains under quiescent conditions throughout the range 0.3 ≤ *De* ≤ 1.5, whereas the position of the second peak is a strong function of flow strength and shifts to higher values of extension as *De* increases. In this intermediate *De* range, there are apparently two stable states (coil and stretched) that do not disappear even after about 40 strain units. Although it cannot be stated conclusively that further simulations would not result in a unimodal distribution indicating a single steady-state configuration profile, 40 strain units represent a deformation of exp(40), which is highly suggestive that a statistical steady state has been achieved.

The coexistence of coil and stretched configurations has been observed both theoretically and experimentally for dilute solutions, but these simulations provide the first evidence of the possibility of a similar phenomenon being present in dense polymer melts. Additional support is provided by examination of similar behavior exhibited by the distributions of the radius of gyration and the entanglement number, as presented in the supplementary material. By analogy, it is reasonable to postulate that a similar coil-stretch transition and conformational hysteresis can also exist in entangled polymeric liquids. It remains to be discussed whether such physical behavior can result in hysteresis of the configurational flow profile, as has been observed in the case of dilute solutions.

At *De* = 3.0 and 14.9, Figs. 4 and 5 reveal that all the macromolecules within the C_{1000}H_{2002} melt are highly stretched and fluctuating mildly about a definitive steady state. These steady states can be used as initial conditions for further simulations that reduce the strain rate to a lower value, thus allowing the macromolecules to relax into the apparent bistable regime. The left panels of Fig. 6 display simulated data for the case where a steady state flow system at *De* = 3.0 was suddenly reduced to *De* = 0.6. The macromolecules react very rapidly to lower their average fractional extension slightly and then tend to fluctuate mildly about the new steady state. Furthermore, the associated distribution is unimodal and contains only the stretch peak at high values of extension. This behavior is very different from that of the corresponding panels of Figs. 4 and 5, where clear signs of a bistable steady state were evident. At very high values of strain, the data suggest that a few molecules may fluctuate enough given sufficient time to surmount the activation barrier to achieve the second configurational state, but these events represent only a minor fraction of the total number of chains and apparently only occur at rather high strain values. Similar behavior is observed after a reduction from *De* = 14.9 to *De* = 1.5, as displayed in the panels on the right side of Fig. 6. Note that the extension values of the single peaks at both *De* in the distribution plots of Fig. 6 coincide with the stretch peaks of Fig. 5, indicating that they are associated with a stable steady-state highly stretched configuration that is energetically distinct from the coiled conformation. Also evident is the verification of the premise that the configurational state of the system is highly dependent on the flow deformation history within the intermediate flow strength regime.

Under the hypothesis that each peak of the distributions shown in Fig. 5 corresponds to a distinct configurational steady state of the C_{1000}H_{2002} melt undergoing planar elongational flow, a hysteresis flow profile can be constructed. Figure 7 shows the steady-state fractional extension as a function of *De* for the C_{1000}H_{2002} melt as extracted from the peaks of Fig. 5. Note that since the quiescent (*De* = 0) value cannot be shown on a semi-log scale, it is assigned at the value of *De* = 0.01. The S-shaped profile, as envisioned by de Gennes, was drawn over the data as a guide to the eye. This figure provides a clear indication that this entangled liquid exhibits a configurational hysteresis when subject to planar extensional flow, accompanied by a fairly dramatic coil-stretch transition at a critical value of the strain rate which occurs roughly at *De* ≈ 1.5. This value is higher than the theoretical value of *De* = 0.5 predicted for dilute solutions, and the width of the bistable steady-state region (0.2 < *De* < 1.5) is relatively wide compared to that observed for dilute solutions. The implication of this conclusion is that the abruptness of the coil-stretch transition in entangled polymer melts is less severe than that which occurs in dilute polymer solutions, but the overall system configurational state is more likely to be biphasic within the applicable strain rate regime.

The configurational hysteresis evident in Fig. 7 can also lead to a similar phenomenon with respect to the rheological properties of the fluid. For example, the primary extensional viscosity, *η*_{1}, will take on different values depending on whether the liquid is subjected to a step-strain from the quiescent condition [Figs. 4 and 5: *η*_{1}(*De*) = *η*_{1}(0.6) = 418 cP; *η*_{1}(1.5) = 393 cP] or from a pre-strained state [Fig. 6: *η*_{1}(0.6) = 834 cP; *η*_{1}(1.5) = 487 cP]. Therefore, which value is actually observed under experiment apparently depends on the deformation history of the sample as well as whether or not a true steady state condition has actually been attained.

See supplementary material for further details concerning the simulation parameters and additional results related to the radius of gyration and the entanglement number.

Computational resources for this project were provided by the National Science Foundation under Grant No. CBET-0742679 to the PolyHub Engineering Virtual Organization (EVO) as well as allocation of advanced computational resources by the National Institute for Computational Sciences (NICS) and the Oak Ridge National Laboratory Joint Institute for Computational Sciences. Financial support was provided by the National Science Foundation under Grant No. CBET-1602890.