Predictive relation for the α-relaxation time of a coarse-grained polymer melt under steady shear

See allHide authors and affiliations

Science Advances  24 Apr 2020:
Vol. 6, no. 17, eaaz0777
DOI: 10.1126/sciadv.aaz0777


We examine the influence of steady shear on structural relaxation in a simulated coarse-grained unentangled polymer melt over a wide range of temperature and shear rates. Shear is found to progressively suppress the α-relaxation process observed in the intermediate scattering function, leading ultimately to a purely inertially dominated β-relaxation at high shear rates, a trend similar to increasing temperature. On the basis of a scaling argument emphasizing dynamic heterogeneity in cooled liquids and its alteration under material deformation, we deduce and validate a parameter-free scaling relation for both the structural relaxation time τα from the intermediate scattering function and the “stretching exponent” β quantifying the extent of dynamic heterogeneity over the entire range of temperatures and shear rates that we can simulate.


Glass-forming liquids are often subjected to shear deformation and flow during their processing, leading to substantial changes in their relaxation dynamics in comparison to unsheared materials at equilibrium. Many measurements (14) and simulations (57) have shown that the structural relaxation time and effective viscosity of deformed complex fluids can decrease following material deformation, sometimes by orders of magnitude. Measurements have also shown that the relaxation dynamics of deformed glassy materials tend to recover to their equilibrium properties after long “aging” times (8), suggesting some sort of structural and dynamical recovery of the material. Despite the abundant experimental literature, there is limited understanding of these changes in dynamics from a theoretical perspective. This situation is perhaps not unexpected, since there is not even a consensus about the origin of the dynamics of glass-forming liquids under equilibrium conditions.

In the present work, we address the essential phenomenology of a model glass-forming polymer melt subjected to steady shear in the liquid regime at temperatures below the onset temperature TA, where non-Arrhenius α-relaxation becomes prevalent (9). We restrict our attention to relatively short chains to avoid complications arising from chain entanglement. First, we characterize the slowing down of the α-relaxation through the study of the α-relaxation time τα0 as a function of temperature in the absence of shear. τα0 was determined in a standard fashion from the intermediate scattering function Fself(q0, t) at a scale q0 corresponding to the intermolecular distance, as determined by the primary peak of the static structure factor S(q). These equilibrium data are compared to the relaxation time τα(γ·,T) of the fluid after the material has settled into a steady state (i.e., shear transients are not considered in the present work). We find that the α-relaxation component of Fself(q0, t) becomes progressively suppressed with increasing γ·, ultimately merging with the inertial fast β-relaxation process. The β-relaxation component of Fself(q0, t) at fixed T is insensitive to γ·. Specifically, the α-relaxation time τα decreases with increasing γ·, ultimately merging with the β-relaxation time, and τα follows a universal scaling relation as a function of a reduced stress, with apparently no adjustable parameters for both polymeric and small-molecule glass-forming liquids. The stretching exponent β of the α-relaxation in Fself(q0, t) also shows a general transition from a diffusion-dominated regime to an inertially dominated ballistic regime and can be described by its own reduced variable description. In the context of the critical phenomena of fluids near their critical point for phase separation, it is well known that shear causes a progressive breakdown of transient particle clusters arising in connection with the incipient phase separation process, so that the lifetime of these transient clusters becomes independent of the average cluster size, as measured by the correlation length ξ for density and composition fluctuations in one and two component phase separation processes (10, 11). This mechanism of shear thinning is probably very common in complex fluids exhibiting dynamic self-association, and we may likewise expect large changes in the stress relaxation time τα of glass-forming liquids under steady shear deformation to be accompanied by the breakdown of dynamical immobile particle clusters under steady shear that should occur generally in this broad class of fluids (12, 13). Correspondingly, our analysis of shear thinning in our model polymeric glass-forming liquid is based on the working hypothesis that shear thinning in glass-forming liquids arises from the breakdown of immobile particle clusters under flow rather than from an alteration of the activation energy as in the conventional Eyring model for relaxation in deformed fluids (3). Previous work has shown that the lifetime of the immobile clusters in both polymeric and nonpolymeric liquids can be identified with the structural relaxation time extracted from the intermediate scattering function under equilibrium conditions (12, 13), and the present work shows that this interpretation extends to our fluid under steady shear.

Our modeling of the reduction in τα under steady shear in terms of the breakdown of dynamic clusters of immobile particles accords with the mechanism of the reduction in the order parameter relaxation time under steady shear in critical fluids, which is one of the few complex fluids for which there is a rigorous, and extensively experimentally validated, theory of the nonlinear rheology of a real complex fluid (10, 11). Nonetheless, our modeling of the shear thinning in glass-forming liquids is just a working model of a complex phenomenon for which there is no rigorous theory, even in the absence of shear deformation.


The effective shear viscosity η, defined during the steady state as η=σ/γ·, is reported in the main panel of Fig. 1. The effective η measured during shear approaches the equilibrium η in the limit of low γ·. With increasing γ·, shear thinning occurs as the viscosity drops by orders of magnitude and the T dependence decreases up to the point where it is not appreciable for high shear rates.

Fig. 1 Shear rate dependence of viscosity and structural relaxation time.

Main panel: Effective viscosity of all systems derived as η=σ/γ·, where σ is the steady-state value of the stress reported in the inset of Fig. 5. The drop in viscosity η with increasing γ· is a signal of shear thinning for this model. A strong temperature dependence is observed for low shear rate, while it is not apparent for shear rates γ· values higher then γ·=101. Lines are a guide for the eye. Inset: α-Relaxation times of the systems derived from Fself(q0, t), as indicated in Fig. 2 and Eq. 2 and associated discussion.

Shear thinning is a common phenomenon in complex fluids (14) in which the structural organization of the fluid is broken down under shear. The shear thinning of solutions of self-assembled actin and tubulin molecules provides a protypical example of this phenomenon (15). As noted before, previous work has identified growing dynamical clusters of immobile particles forming in our model glass-forming liquid (12), so this interpretation of shear thinning in terms of “depolymerization” under flow seems plausible also for glass-forming liquids (we do not refer here to the breaking of chemical bonds, but rather to the transient associative bonds between molecules, as normally encountered in systems undergoing supramolecular assembly). We adopt this dynamic heterogeneity perspective below in our quantification of the shear-thinning effect that we observe.

To describe relaxation in our system under shear, we consider the self-part of the intermediate scattering functionFself(q0,t)=1Nj=1Nexp{iq·[rj(t)rj(0)]}(1)for all the systems, where q0 ≃ 7.14 is the value at which the maximum of the static structure factor is observed. Following previous protocols (7), only q vectors in the yz plane orthogonal to the direction of the shear are considered.

Fig. 2 shows the intermediate scattering functions at the lowest temperature investigated for all shear rates. It is apparent that with increasing shear rate, the α-relaxation is suppressed, the α-relaxation time decreases, and, eventually, the α- and fast β-relaxations merge for high shear rates, similarly to the effect achieved by increasing T. Similar results have been observed before in two-dimensional (2D) (6) and 3D (7) model atomic liquids. To quantify these results, we fit our data to a two-step relaxation function, successfully used to describe Fself in thin films and bulk materials at equilibrium (16)Fself(q0,t)=(1Aα)e(t/τf)βf+Aαe(t/τα)β(2)

Fig. 2 Intermediate scattering function over a range of shear rates showing the suppression of the α-relaxation process.

The self-part of the intermediate scattering function is reported for all γ· for the temperature T = 0.43. The brown curve labeled as γ·=0 is the reference equilibrium system before shearing. Circles are simulation data, and lines are fits to Eq. 2, with the fitted constants τf = 0.3, βf = 1.3, and A = 0.73. The dashed black line indicates the onset of a plateau following the initial fast decay, corresponding to the limit value of the intermediate scattering function in the limit of infinite α-relaxation time. The height of the plateau is 0.73, coinciding with the value of the estimated A parameter for this T.

In qualitative accord with Eq. 2, the mode coupling theory likewise indicates (17) an additive contribution of a fast β-relaxation, with a T-independent β-relaxation time, and α-relaxation process, whose relaxation time is strongly dependent on T and whose amplitude is independent of T above the “glass transition” of this model, Tg. Correspondingly, we tentatively take τf and βf to be fixed and determine α-relaxation time τα, Aα, and stretching exponent β as function of T and γ·. The equilibrium, zero shear reference case is discussed in the Supplementary Materials.

We find that the values of the parameters τf and βf of the initial “fast” relaxation process (see Fig. 2) can be fixed to τf = 0.3 and βf = 1.3, respectively, defining the fast β-relaxation time as the exponent βf quantifying the nonexponential nature of the fast relaxation process. The “mode amplitude” parameter Aα shows a relatively weak T and γ· dependence for the parameter range studied in the present work, and, correspondingly, we fix Aα to equal Aα = 0.7, except at the lowest temperature T = 0.43, where this parameter increases to a somewhat larger value, Aα = 0.73. The accuracy of the fit is not appreciably reduced in the T range investigated by fixing the three parameters (τf, βf, and Aα). The lines in Fig. 2 are the fit of Eq. 2 to the data (circles) for all shear rates. From the fits of Eq. 2 to the simulation data, values of τα and β can be derived for all T and γ·. It is apparent from Fig. 2 that the α contribution to Fself(q0, t) becomes increasingly “suppressed” with increasing γ·, an effect superficially similar to increasing the temperature of the fluid [see the Supplementary Materials for an extended discussion of the effect of varying T on Fself(q0, t)]. In particular, we observe that τα becomes progressively reduced by γ·, reaching the value τα0(TA) (the value of τα0 at equilibrium at T = TA, where TA ≈ 0.7 for this model) at a characteristic shear rate γ·c101. Previous work defined τβ as τα0(T=TA) at equilibrium (18). For the present model, we find τα0(TA)=2.12, or ≈2 ps in laboratory units. γ·c is slightly T dependent and can be roughly estimated from the relaxation time data shown in the inset of Fig. 1. Below, we use these values of γ·c to estimate a corresponding critical steady-state stress of the system σc=σ(γ·c) (see the inset of Fig. 5 in the Materials and Methods section), which plays an important role in our reduced variable description of τα(γ·). For γ·>γ·c, the α-relaxation process appears to merge with the β-relaxation process when τα(γ·) obtains a value on the order of the “fast” relaxation time, τf, a time scale on the order of a picosecond in laboratory units.

We now show that the changes in τα(γ·) and β(γ·) follow an apparently universal scaling behavior, and these quantities can be scaled onto master curves involving a reduced shear rate covering the full range of T and γ· where there are no fitted adjustable parameters. The behavior of the relaxation time is similar to that of the effective viscosity, with an initial T-dependent plateau followed by a decrease with increasing γ· up to the point where the T dependence vanishes and the τα value is comparable with the fast β-relaxation time.

We draw upon basic scaling theory, symmetry arguments, and previous modeling of the shear rate dependence of the structural relaxation in sheared near-critical fluids to develop an approximant for the shear rate dependence of the structural relaxation time τα(γ·) of glass-forming liquids. First, we note τα is invariant to the choice of the direction of shear so that τα(γ·) must be an even function of γ·. A natural first guess would be to assume that τα(γ·) is an analytic function of γ·, so that τα(γ·) would then be described by a series of even powers of γ·. Previous space shuttle measurements on the shear thinning of Xe near its critical point (4) provide some additional guidance regarding the form of this expansion, and the theory of shear thinning in fluids near their critical point gives us some physical and mathematical insights into how to model the shear thinning in glass-forming liquids (10). In particular, measurements on sheared near-critical Xe have indicated that the shear viscosity under steady shear is better described as a function of γ·τα(γ·=0) rather than even powers of the dimensionless shear rate, γ·τα(γ·=0). Correspondingly, we adopt the assumption that τα(γ·) is a function of γ·τα(γ·=0), on the basis of a presumed physical analogy between shear thinning in critical fluids and glass-forming fluids discussed below.

The dynamical renormalization group (RG) theory of shear thinning in critical fluids indicates that the reduction in the “order parameter relaxation time” τ, the viscoelastic relaxation time describing the lifetime of dynamic particle clusters associated with incipient phase separation, can be understood from the progressive breakdown of these clusters under steady shear (10, 11). In the absence of steady shear deformation, the shear viscosity and normal stresses of the near-critical fluid diverge as a power of the correlation length ξ describing the average size of the dynamic particle clusters, where the precise power relating the zero-shear viscosity η(γ·) to ξ, and also the shear thinning exponent describing reduction in τ under shear are both precisely predicted by the dynamical RG theory. Physically, the infinite order ϵ-expansion shear thinning exponent for the shear viscosity can be exactly understood to arise physically from the breakdown of the critical clusters under flow (11), so that the shear viscosity of critical fluids is completely independent of shear rate at high shear rates (10), i.e., the critical clusters are completely disintegrated at high shear rates. Below, we estimate the shear thinning exponent describing the dependence of τα(γ·) on shear rate based on the assumption that the dynamic clusters of immobile particles in cooled liquids (12, 13) are likewise progressively broken down with increasing shear rate, so that τα(γ·) is likewise independent of τα(γ·=0) in the limit of high shear rate. Our model of shear thinning in glass-forming liquids is then predicated on the hypothesis that the dynamic clusters of immobile clusters of chain segments (12, 13) are the fundamental structural elements in glass-forming liquids responsible for the growth of the viscosity upon cooling.

The next step in our scaling argument requires a specification of the general functional form τα(γ·) that is applicable for any steady shear rate between the low shear rate limit of the equilibrium glass-forming liquid and the high shear limit, where the dynamic heterogeneity is presumably suppressed, and relaxation correspondingly accelerated. The RG theory takes a general “crossover” form in the leading-order calculation of many properties of condensed materials transitioning between limiting behaviors (or fixed points of the renormalization scheme). We adopt this crossover form as natural functional form where the limiting behaviors are the homogeneous fluid and the dynamically heterogeneous fluid exhibiting large transient clusters of immobile particle tending to increase the fluid viscosity as in near-critical fluids. In particular, we adopt the functional form (19)τα(γ·)=τα0{1[1+(γ·τα0φ/γ*)]α}{1+aλ+O(λ2)}(3a)λ=[γ·τα0φ/γ*]/[1+(γ·τα0φ/γ*)](3b)where a is a small constant calculated from the RG theory, τα0τα(γ·=0), and γ* is a constant describing the rate of crossover between the zero shear and high shear rate limiting scaling behavior. On the basis of the physical arguments given above, we also specify the “crossover exponent” ϕ to equal 1. Further, in our initial work, we are interested in just the leading-order scaling behavior of τα(γ·), so we adopt the simplified crossover equationτα(γ·)τα0{1[1+(γ·τα0/γ*)]α}(4)in an attempt to reduce our simulation estimates of τα(γ·) to a nearly universal reduced variable description. The prediction of the dependence of τα(γ·) on γ· requires some means of specifying the shear thinning exponent α and γ*, and we invoke a physical model of the influence of shear on glass-forming liquids based on the well-studied example of shear thinning in near critical liquids. In particular, we observe that there is only one choice of the exponent α (i.e., α = 1) that is compatible with τα(γ·) in Eq. 4 becoming independent of τα0 in the limit of high shear rates. The liquid in this high shear rate limit is presumed to be dynamically homogeneous, in the sense that the immobile clusters are broken down by shear. Note that both Eqs. 3 and 4 correspond to the Carreau-Yasuda class of functions that are often used phenomenologically to describe the shear dependence of the viscosity in complex fluids (20), but our motivation for invoking this relation based on RG theory crossover scaling is different from conventional arguments for this equation.

An acceptable scaling expression for τα(γ·) must also address the fact that the relaxation time cannot become arbitrarily small at high shear rates, and this requires an additional modification of our scaling form for τα(γ·). In particular, we see in our simulations that τα(γ·) decreases progressively with shear until the α- and fast β-relaxation processes merge at some “critical” shear rate. This means that τα(γ·) must approach τβ at high shear rates, τα(γ·)=τα=τβ, where τβ is a temperature-insensitive time on the order of 1 ps in laboratory units. To account for this limiting behavior, we then revise our expression to the simple formτα(γ·)τατα0τα=1[1+(γ·τα0/γ*)](5)

We find that our data best fit Eq. 5 by taking τα=0.5. Since our dimensionless time unit is on the order of 1 ps, and τβ for our model glass-forming liquids and glass-forming liquids in general is on the order of 1 ps (21), then τα has the same order of magnitude as τβ. More generally, this characteristic time probably has a weak temperature dependence (as in the case of τβ), but this small variation is unessential in our data reduction below, as this only affects the data points at the extremely high shear rates. An analysis of the robustness of the data collapse as a function of τα=0.5±0.1 is shown in the Supplementary Materials.

Many previous theoretical and experimental studies have fit shear relaxation data to the functional form of Eq. 4 or 5 or the more general Carreau-Yasuda functional form where ϕ is commonly taken to equal 2 and the constant of proportionality in the reduced shear rate is taken as an adjustable “nonuniversal” constant [see table 1.7 of Tanner’s Engineering Rheology (22)]. Correspondingly, we have found that we can fit our τα(γ·) data rather well following this procedure, where we find that γ* is an empirical constant whose significance was at first obscure to us. Previous experimental studies have introduced a reduced variable of the observed shear thinning that we have found to provide insight into γ*. In particular, previous experimental investigations (23, 24) of η(γ·) of some glass-forming polymeric liquids indicated that η(γ·) could be expressed in terms of a nearly universal function of a dimensionless shear rate defined by the ratio of the shear stress of the fluid in its reference zero shear Newtonian state, σN=η(γ·=0)γ·, to a critical stress σc signaling the onset of “elastic turbulence” where a sharp change in the “fluidity” of the melt and flow instabilities in polymer melt extrusion measurements are observed. The utilization of a physical condition for defining a reduced shear rate offers an attractive method for making our expression for τα(γ·) more predictive, since if we could determine γ*, there would be no free parameters in our expression for τα(γ·). We next consider a reduced stress description of our structural relaxation time data to specify the undetermined parameter γ* in Eq. 5.

Within the RG scaling scheme discussed above, it is natural to define a dimensionless shear stress by a quantity that delineates an approach to a strong coupling scaling regime, and we follow former experimental practice in how we identify our reduced stress variable to obtain a unified description of the shear dependence of τα(γ·) at different T and γ·. In particular, (23, 24) have identified a critical shear stress based on the observation of a transition to unstable fluid flow, and this stress was used to define a reduced shear stress. Correspondingly, we identify the condition at which τα merges with τβ as corresponding to this flow instability condition, because this condition implies that the relaxation of the fluid by momentum diffusion (shear viscosity) has been essentially suppressed, so that the fluid can be viewed as roughly equivalent to an inviscid fluid from a modeling standpoint. “Turbulent” or “chaotic” flow is a common phenomenon in the flow of fluids having a low viscosity, and even if there is no exact theory of this ubiquitous phenomenon, this general tendency can be qualitatively understood from recent computational investigations of particle motion in idealized inviscid fluids, where convincing numerical evidence is presented that the motion of a solid body through an incompressible, inviscid fluid, moving irrotationally and otherwise at rest, is inherently chaotic (25). On the basis of previous experimental practice and physical argument just stated for identifying the approach of τα to τβ with the critical condition for unstable flow, we define our reduced stress σR by σNc.This critical stress can be estimated both by simulation and experiment, at least through extrapolation. The β-relaxation time of the unsheared material is on the order of 1 ps in the Fself(q0, t) data discussed above [τβ is on the general order of 1 ps for molecular liquids (21)] and in measurements (26), so τα0 also becomes on this order of magnitude when T approaches the onset temperature TA for non-Arrhenius relaxation (18). If we define the critical shear rate γ·c by the condition at which τα(T,γ·) equals τα(TA,γ·=0)O(1012 s), we can estimate a critical shear stress σc=σ(γ·c) from either simulated and measured σ(γ·) curves (see inset of Fig. 5 in Materials and Methods). Consistent with this criterion, previous experimental studies have identified a characteristic temperature (“pole temperature”) where polymer fluids appear to lose their cohesive nature and then flows with little resistance, a temperature that has been estimated for many polymers (27). To further specify an appropriate stress variable that is defined in terms of experimentally and computationally measurable properties, we also need an estimate of the relation between η(γ·=0,T)η0 and the structural relaxation time τα. For this estimate, we assume the validity of the Maxwell relation η0 = Gτ for the fluid at equilibrium, where G is the experimentally measured shear modulus Gg in a high frequency shear relaxation measurement, and τ is the shear stress relaxation time. In a glass-forming liquid at low temperatures, G corresponds to the “glassy modulus” and τ to the relaxation time, but in the simple fluid regime, where appreciable chain diffusion occurs on the relaxation time of the fluid, G corresponds to the Rouse or entanglement modulus and the relaxation time corresponds to the Rouse or “terminal” relaxation time of the polymer melt, respectively. Our focus on τα in the present paper allows us to focus on the influence of shear on the segmental relaxation time, which is not sensitive to either molecular mass or large-scale polymer relaxation associated with large scale chain diffusion. We may then reasonably take G in the Maxwell relation to correspond to the glassy modulus and the relaxation time τ to correspond to the segmental relaxation time, which is of the order of τα. Recent work has shown that Gg of the glass-forming material should not be identified with the infinite frequency shear modulus G, but rather with an intermediate time “plateau” in the shear relaxation function corresponding to the glassy relaxation process, in analogy with Aα in the decay of the intermediate scattering function (see Eq. 2) (28). Consistent with our proposed data reduction in our τα(γ·) data, the observed (28) high frequency modulus Gg has previously been used for the effective data reduction in shear thinning measurements on the shear viscosity of soda-lime-silica glasses (29). We note that the shear viscosity of polymeric glass-forming liquids inherently involves diffusion and relaxation processes at the scale of the polymer, so the viscosity of the polymer fluid should be much more sensitive to the polymer massthan τα and it is not clear that the reduced shear rate deduced above for τα should apply to the shear viscosity of polymer melts under steady shear. The α-relaxation time is a probe of dynamics at a segmental scale and is an inherently simpler quantity than the shear viscosity, diffusion coefficient, normal stresses, and other polymer melt transport properties.

Now, we have defined our reduced stress variableσR=σN/σcGgγ·τα0/σc(6)assuming the validity of the Maxwell equation for the equilibrium fluid (as noted above), the identification of the shear stress relaxation time τ with the segmental relaxation time τα obtained from the intermediate scattering function (6, 12, 24), and the definition noted above of the critical stress σc in terms of the magnitude of the critical strain rate γ·c at which the α- and β-relaxations merge. The parameter γ* of Eq. 5 is then identified with the ratio σc/Gg. Eq. 5, in conjunction with the approximate relation of Eq. 6, then provides a predictive relation for τα(γ·,T) without any adjustable parameters. We show the comparison of Eq. 5 to our simulation data in Fig. 3. Here, we have estimated Gg(T) from τα0 based on the former calculations of Puosi and Leporini (28). We observe that Eq. 5 describes our τα(γ·) data very well over the entire accessible range of γ· and T.

Fig. 3 Reduced variable description of α-relaxation time and stretching exponent β as a function of dimensionless shear rate.

Main panel: Relaxation times τα(γ·) in the inset of Fig. 1 collapse onto a near reduced variable description after normalizing the y axis by the limiting values τα and τα0 and using the reduced variable σNc for the x axis. Details of this data reduction are described in the text. The data then comply with the crossover function of Eq. 5. σNc involves the ratio Ggc, which ranges between 6.75 and 4.5 in the temperature range investigated. Note that the exponent of the power law −1 is robust, as even a change of the exponent by a few percent worsens the agreement with the data. Inset: Collapse of the β(γ·) stretching exponents onto an approximate master curve after normalizing by a T-dependent shear rate γ·β.

It is also notable that the shear thinning that we observe, and its physical interpretation of the breakdown of self-assembled structures under flow, seem to be broadly consistent with observations on naturally occurring complex fluids exhibiting thermally reversible self-assembly (actin, tubulin, etc.) (15). All these fluids exhibit a similar “shear thinning” at high shear rates with a shear-thinning exponent near −1 in the limit of high shear rates, defining a hydrodynamic limit of flow instability, i.e., the flow shear stress becomes essentially indeterminant (15). We note that the shear thinning that we observe apparently has little to do with the polymeric nature of fluid but rather derives from the emergent dynamic heterogeneity of the cooled glass-forming liquid as these fluids are cooled. This is evidenced by the essential disappearance of the shear thinning upon raising the temperature toward the onset for non-Arrhenius relaxation TA and by the observation of a very similar shear thinning in binary Lennard-Jones (LJ) fluids introduced to model metallic glass materials (see the Supplementary Materials). We view the polymeric nature of our glass-forming liquid as a convenient way to suppress any tendency toward crystallization so that dynamic heterogeneity occurs as an equilibrium phenomenon. The dynamics of crystallizing glass-forming materials, such as the binary LJ liquid, is complicated by equilibrium “aging” phenomena associated with the metastable nature of the supercooled state, yet this class of glass-forming liquids seems to exhibit a similar dynamics to our model polymeric glass-forming liquid.

By a similar reasoning to Eqs. 5 and 6, we may also deduce an approximate reduced variable description for the stretching exponent β(γ·) of the intermediate scattering function based on the same physical picture of a breakdown of immobile clusters under shear. Stukalin and co-workers (30) have shown that the β exponent for stress relaxation in solutions in which there are growing dynamic clusters can be exactly described by a sigmoidal variation between 1 at high temperatures to a value that saturates near 1/3 at lower temperatures, when the clusters have fully grown to large dimensions and persist for very long times. Since increasing shear has an effect on the immobile clusters that is similar to increasing temperature, we may generally anticipate a sigmoidal variation of β with shear rate as well. This expectation is again based on assumption that the decay of the intermediate scattering function and the structural relaxation time are physically closely related, a common, but unproven, assumption that explains why τα is designated the “structural relaxation time.” At any rate, simulations by Riggleman and co-workers (3, 31), although for a different mode of deformation, have indicated a general tendency of β to increase under deformation, consistent with our cluster breakdown interpretation of the variation of β with applied shear. For our system, β increases monotonically with increasing γ· from a low shear rate value near β0 near 2/3 at γ·=0 up to a high shear rate value β that is larger than unity. β apparently approaches a limiting value around 1.45 for large γ·, a value that is close to βf for the fast β-relaxation process. Here, we have further evidence that the α-relaxation process is being suppressed at high shear rates, where the overall relaxation progressively becomes more like the β-relaxation process seen only at short times in the quiescent fluid. A shortcoming of the equilibrium polymer model of stress relaxation for glass-forming liquids is that this model does not account for the fact that the immobile particle clusters form randomly branched polymers, which makes the low-temperature value of the β exponent somewhat larger in this type of dynamically heterogeneous fluid (32). Moreover, this model of a fluid exhibiting dynamic clustering corresponds to particles in two states at any given time: particles in polymeric clusters or undergoing Brownian motion in their “mobile” or “unassociated” particle state, so the short time dynamics in the equilibrium polymerization model is Brownian rather than inertial. Despite the quantitative disparities of the solution chain assembly model and the immobile clustering process observed in simulations of glass-forming liquids, the overall trend that we see in the shear dependence of β is qualitatively consistent with the expected behavior arising from immobile particle cluster breakdown under shear in glass-forming liquids, along with the contributing factor of the emergence of relaxation dominated by inertia at high shear rates. In particular, the transition from β0 to β is consistent with a transition between diffusion-dominated relaxation to a predominantly inertially dominated relaxation regime where the α-relaxation is entirely suppressed. For each temperature, we then consider the characteristic shear γ·β for which β(γ·)β0ββ0=0.5, the midpoint of the transition between the two regimes. γ·β decreases with decreasing T, going from 0.077 at T = 0.55 to 0.013 at T = 0.43. Introducing the reduced variable γ·/γ·β described above, we find that we can achieve a collapse of all of our β data, as shown in the inset of Fig. 3. In particular, β(γ·) can be described phenomenologically by the simple sigmoidal function corresponding to the function λ in Eq. 3bβ(γ·)β0ββ0=γ·/γ·β1+γ·/γ·β(7)

Our scaling expression for τα(γ·,T) is based on the hypothesis that the growth of τα0 upon cooling is associated with the growth of stress-bearing and persistent clusters of immobile polymer segments in cooled liquids. We are not as confident about the universality of this data reduction in β(γ·) as our τα(γ·) data because the determination of the crossover shear rate γ·β relies on curve fitting, rather than on a predicted criterion based on a physical argument. We also do not presently see a fundamental reason why the magnitude of β should generally equal 1.5 for all liquids. Nonetheless, this cross-function should provide a good means of quantifying this crossover from a purely correlative perspective.

Previous simulation studies by Starr et al. (12) have shown that fractal dynamic particle clusters of relatively immobile polymer segments exist in our coarse-grained polymer model and that the lifetime of those clusters is proportional to τα0 to an excellent approximation. Weitz and co-workers and Tracht and co-workers have also emphasized the importance of such immobile particle clusters for understanding stress relaxation in glass-forming liquids (33, 34). Riggleman and co-workers have further shown that immobile clusters exist in the glassy state, where they persist for extremely long time scales (35). Evidence for immobile particle clustering can also be inferred from nanoscale imaging of the interface of polymeric and metallic glass-forming materials (36). We note that our interpretation of the reduction in the relaxation time through the breakdown of immobile particle clusters is consistent with the recent comprehensive experimental studies by Ediger and co-workers (2), where it is emphasized that a wide range of measurements indicate that the application of stress accelerates mobility in immobile regions in the material much more than the mobile regions, as evidenced by the increase in β under deformation sufficient to cause flow. Correspondingly, we next check our hypothesis that the variation of τα(γ·) and the β exponent is linked to immobile particle clusters by directly calculating how the cluster size is altered under flow.

Starr et al. (12) define an autocorrelation function defining the persistence of the immobile particles in their “caged” state, and this method was also used in a recent work on metallic glass-forming liquids (13), where similar results were obtained. We follow these previous works in our definition of the immobile particles and then determine how the application of shear alters the persistence time of the immobile particle clusters and how this time relates to τα(γ·,T). Figure 4 shows the fraction of caged particles as a function of γ·. See Starr et al. (12) for a detailed discussion. In the present work, we confine our attention to the yz plane to be consistent with our analysis of Fself(q0, t) to consider how γ· alters the caged fraction C(t) at the relatively low temperature T = 0.43.

Fig. 4 Time dependence of the fraction C(t) of caged particles confined by neighbors for T=0.43 and variable γ·.

The decay time for the fraction of caged particles is comparable to the shear-dependent relaxation time shown in Fig. 2. The inversion point at t = 1 is due to the definition of the caged fraction and the decreasing dynamical heterogeneity for higher γ·. The two snapshots use larger, colored spheres to highlight the less mobile particles measured after t = τα at equilibrium and for the highest γ·, respectively. The colors correspond to the size of the clusters of neighbor immobile particles. The large clusters in the equilibrium system (blue beads in the top right snapshot) are broken down and spatially dispersed under high γ·. More quantitatively, the average cluster sizes for the two snapshots are 4.0 and 2.4, respectively.

Consistent with the equilibrium analysis of the lifetime of the immobile clusters under variable temperature conditions (12), the decay time is comparable to the structural relaxation time τα derived from the intermediate scattering function, becoming shorter with increasing γ·. In Fig. 4, we also show two snapshots representative of the behavior of the clusters of immobile particles. The 10% less mobile particles of the system [a rough estimate of the fraction of caged particles, see Figure 23 of (12)] in a time τα(γ·) are highlighted for the system at equilibrium and at the highest γ·. Different clusters of neighbor particles (chosen with a cutoff radius of rc = 1.46) have a different color. It is apparent that while the equilibrium system shows larger immobile particle clusters, they are indeed progressively broken down into smaller and more spatially heterogeneous clusters with increasing γ·, an effect comparable with an increasing T at equilibrium.

Of course, the mobile particles also play an important role in glass-forming liquids at equilibrium (12) as these particle determine the rate at which particle exchange between their mobile and immobile states, and string-like motion has been observed in experimental systems under shear (15, 14). A more detailed analysis of the string-like motion associated with the most mobile particles under shear will be presented in future works (see Discussion).


Given the practical importance of material deformation and flow on the processing and performance of glass-forming polymer materials, there have naturally been numerous experimental and computational studies of how relaxation and transport processes in glass-forming materials are influenced by molecular and thermodynamic properties and flow conditions. There is an extended scientific literature on this topic, recently reviewed by Hebert et al. (37), who summarize the findings of numerous experimental and computational studies. The measurements summarized in this work indicate, under constant strain rate conditions and for temperatures 10 K to 20 K below Tg, that the relaxation time in the postyield regime decreases with an apparent power law behavior, with an exponent close, but slightly below, −1. A corresponding reduction in the relaxation exponent β with strain rate was also observed. For example, Hebert and co-workers indicate a shear thinning exponent of about −0.8 to −0.9 for cross-linked polymers slightly below Tg, which compares well with our simulations at a much higher shear rate and at a temperature above Tg. For example, a log-log plot of our raw data of τα as a function of γ· at the lowest temperature simulated (T = 0.43) indicates a power law dependence with a slope of −0.79 ± 0.03 (see the Supplementary Materials for graphical illustration of this finding). The dynamic mean field theory of Schweizer and co-workers (38) also indicates an effective power law having an exponent with similar magnitude (≈0.86) for shear thinning flow when the results are analyzed in a similar fashion. We note that the model of Schweizer and co-workers (38) and its antecedent mean-field model by Eyring, emphasizing changes on the activation energy for relaxation under material deformation, cannot address changes in β since this class of models assumes that glassy materials are perfectly homogeneous.

There has apparently been no previous attempt to obtain a universal reduced variable description of τα(γ·,T) in simulated glass-forming liquids, as we have obtained in the present paper. Previous experimental studies on inorganic (4) and polymeric liquids (23, 24) have introduced reduced variable descriptions of the shear viscosity that are close in spirit to our simulation findings, however. Our τα(γ·) data reduction is based on fluid properties that may be determined either by experiment or simulation, so it should be possible to test our crossover expression directly with experimental measurements of intermediate scattering function or other measurement of segmental relaxation in polymeric or small-molecule liquids. We also gain qualitative insight into previous experiments into the shear rate dependence of segmental relaxation in glassy polymer materials (1, 2). The small variability in the apparent inverse power law scaling with shear rate probably derives from the fact that power law scaling in Eq. 5 is only approached asymptotically at extremely large strain rates. We then suggest that the somewhat smaller exponents, sometimes reported in glass-forming liquids and also seen in measurements on sheared self-assembled protein solutions (15), are probably a consequence of this type of crossover effect (see the Supplementary Materials).

We briefly comment on modeling changes of τα based on the assumption of high shear limit corresponding to a quasi-static deformation of the material in the limit of zero temperature. Our simulations indicate that the dynamics in the high shear limit approaches a state that is more similar to a gas than a zero-temperature solid. For our highest shear rates, the structural relaxation time is reduced by steady shear to a point that the material can no longer support any shear stresses at any frequency and particle motions become nearly ballistic rather than caged by surrounding particles (see the Supplementary Materials where the near ballistic nature of molecular displacements are indicated at high shear rates). At extremely high shear rates, but lower than the “critical” value at which τα(γ·)τβ, thermally activated transport remains important and an increasing shear rate plays a role in reducing the structural relaxation time that is similar to increasing temperature. There is though no perfect equivalence between varying shear rate and temperature, since shear rate and temperature have different influence on the modes amplitude Aα (see the Supplementary Materials). Moreover, the potential energy landscape becomes fundamentally altered when shear is imposed so that drawing an equivalence of the material under steady shear to a material at equilibrium cannot be possible beyond the weak deformation case.

This lack of equivalence between strain rate and temperature has also been demonstrated in the four-point density correlation analysis where anisotropy in transport is indicated in the steadily strained material (6). The two-point density correlation as well as the static and dynamic structure factors are insensitive to this shear-induced anisotropy, and we view this as roughly supporting the idea of a shear rate–dependent effective temperature in a mean field sense. At any rate, relaxation in sheared glass-forming liquids remains an activated process, and increasing the shear rate progressively suppresses α-relaxation processes mediated by diffusion in favor of inertial relaxation processes and the same general trends arise when we raise the temperature of the fluid toward its liquid-vapor critical point.

In summary, we find a predictive relation for the structural relaxation time under shear τα(γ·,T) with no adjustable parameters that takes into account the τα0 and viscoelastic modulus Gg of the liquid at equilibrium. The theoretical framework adopted is consistent with the idea that the clusters of immobile particles relevant for the dynamics of supercooled liquids at low temperatures are broken down under shear (37). It should be appreciated that clusters of highly immobile particles that also arise in cooled liquids are implicated in understanding the rate of molecular diffusion and relaxation processes, and those clusters have been emphasized by Onuki and co-workers (39) in understanding shear thinning of model glass-forming liquids.

Mobile particles involving collective particle exchange have been shown to be highly correlated with changes of the activation free energy of diffusion and relaxation so that these clusters must also influence τα through their effect on the rate at which the immobile particle exchange with the surrounding fluid particles. It is this type of change in the activation free energy under deformation that is envisioned in the Eyring model, and evidence indicates that this type of model has merits at low degrees of deformation before the material exhibits “yield” and a transition to material flow. It is anticipated that this distinct type of dynamic heterogeneity (12) plays an important role in understanding shear banding, plastic instabilities, and other instabilities, brought about by the progressive conversion of immobile particles into mobile particles through the application of shear. This idea is supported by experimental reports of increases in mobility by many orders of magnitude in shear bands (40) and modeling indicating that shear banding occurs from the deformation-induced emergence of regions of transiently high mobility in the glass-forming liquid (41). Simulation observations of large-scale rearrangements of relatively mobile particles undergoing cooperative rearrangements in metallic glass materials undergoing plastic deformation (42), and an observed quantitative relation in molecular dynamics simulations between the extent of string-like collective motion and the activation energy of Johari-Goldstein β-relaxation process in metallic glasses, a relaxation process highly correlated with the occurrence of plastic deformation and brittleness of metallic glasses (43), further support this possibility. The emergence and significance of the mobile particle clusters under strong shear conditions have been emphasized by Tanaka and co-workers (6), on the basis of a four-point density correlation analysis in which the corresponding cooperative motions are found to be more prevalent than for Fself(q0, t). There have also been experimental (14) and simulation (44) reports of the appearance of string-like organization of particles in colloids and molecular fluids at high rates of γ·. Preliminary results for our sheared fluid indicate that string-like collective motion is also suppressed under steady shear up to critical shear rates close to γ·c, comparable to an increase in temperature (18). For γ·γ·c, however, we find that the string-like collective motion is enhanced again. We plan to investigate the impact of the mobile particle clusters on the evolution of stress in sheared fluids and under physical aging conditions in future work.

As a final point, we note that care must be taken in inferring change of the shear viscosity based on estimates of the relaxation time τα(γ·,T). Although we may expect η(γ·) to be reasonably estimated by a Maxwell model of a viscoelastic fluid, η(γ·)Ggτ(γ·,T), the effective high-frequency modulus Gg can also be expected to depend on γ·, since the application of shear alters the cohesive interaction strength of the fluid. Consistent with this expectation, we find that η(γ·) can be described by a reduced shear rate expression, but where the shear thinning exponent is closer in magnitude to 2/3 rather than 1. We will report on this phenomenon, along with the dependence of variation of the first and second normal stresses (N1 and N2) under shear deformation, in a separate publication.


We model polymers as fully flexible linear chains of beads linked by harmonic springs. Each chain has M = 20 monomer, far below the entaglement regime. Nc = 500 chains are simulated for each system, for a total of N = 10,000 monomers. Nonbonded monomers belonging to the same or different chains interact with a truncated LJ potentialULJ(r)=ε[(σ*r)122(σ*r)6]+Ucut(8)up to rcut = 2.5σ, σ* = 21/6σ is the position of the potential minimum with depth ε. The value of the constant Ucut is chosen to ensure ULJ(r) = 0 at rrc = 2.5 σ. Bonded monomers interact with a harmonic potential Ub(r) = k(rr0)2 with k = 555.5 ε/σ2 and r0 = 0.97 σ. Henceforth, all quantities are expressed in terms of reduced LJ units, i.e., ε = 1, σ = 1, with unit monomer mass and Boltzmann constant. The reduced units can be mapped onto physical units relevant to generic nonequilibrium fluids, by taking molecular dynamics (MD) time, length, and energy units as corresponding roughly to about 2 ps, 0.5 nm, and 3.7 kJ/mol, respectively. MD simulations were carried out with the LAMMPS (large-scale atomic/molecular massively parallel simulator) code ( The systems were initially equilibrated in the NPT ensemble (constant number of particles N, pressure P and temperature T) using a Nose-Hoover thermostat and barostat with 〈P〉 = 0 to allow full correlation loss of the end-to-end vector of the polymer chains. After equilibration, shear deformation was applied to the simulation box at fixed shear rate in the xy plane. Periodic boundary conditions are used during the NPT equilibration. A range of temperatures from T = 0.43 to T = 0.55 and a range of shear rates from γ·=105 to γ·=100 were explored. Twelve independent replicas of each state were considered to ensure suitable statistical average.

Deformation protocol

Simple shear is performed in the xy plane at fixed rate γ·. The deformation is performed at constant volume. The SLLOD equations of motion are used along with the thermostat, and Lees-Edwards boundary conditions are applied. The microscopic stress σxy is considered as a function of the total strain γ=γ·t and defined as the average value of the per-monomer stress σxyiσxy=1Ni=1Nσxyi(9)where the per-monomer stress in the atomic representation isσxyi=12vjirxijFyij(10)where Fxkl and rxkl are the x components of the force between the kth and the lth monomer and their separation, respectively, and v is the average per-monomer volume, i.e., v = L3/N, where L is the box size.

The main panel of Fig. 5 reports the typical stress loading curve for our systems during shear. The stress along the plane of shearing σxy initially increases with increasing strain for small deformations. After reaching a maximum value (the yield strength), it then drops to a lower value (strain softening) and remains constant. The properties of the stress overshoot are widely studied and present a rich phenomenology (38), but this interesting transient effect is not the focus of the present paper. We focus on the steady-state flow regime that develops after the overshoot. Every system under study is initially deformed up to total strain γ = 102, beyond the overshoot that is observed around γ ∼ O(1) for all temperatures and shear rates studied. The inset of Fig. 5 shows the steady-state stress value for all the T and γ· investigated. All results shown in the Results are taken in the steady-state regime, starting from the end of this initial deformation.

Fig. 5 Shear stress versus total strain should approach steady flow limiting behavior.

Main panel: Typical loading curve of the shear stress as a function of total strain for the system at T = 0.55, γ·=101. An overshoot of the stress is observed after an initial linear increase, with a maximum around γ = 2.6. After the stress drops, a steady flow regime is observed where the system is stationary and a quasi-equilibrium condition is established. Inset: Value of the steady-state stress for all T and γ· considered. A strong T dependence is observed for low γ·, while it disappears for shear rates higher than γ·=101. Lines are a guide for the eye. Similar results were found for an LJ binary glass-forming mixture (7).


Supplementary materials for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: This work was supported in part by National Institute of Standards and Technology (NIST) Award 70NANB15H282. A grant of computing time from IT Center, University of Pisa and Dell EMC Italia is also acknowledged. Author contributions: A.G., J.F.D., F.P., and D.L. conceived the research. A.G. and J.F.D. designed the simulations and analyses. A.G. and J.F.D. conducted the analyses. All authors contributed to the writing of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article