Nonlinear Shear of Entangled Polymers from Nonequilibrium Molecular Dynamics

: This study aims to use molecular dynamics (MD) simulations of Kremer – Grest (KG) chains to inform future developments of modelsof entangled polymer dynamics. We per-form nonequilibrium MD simulations, under shear ﬂ ow, for well-entangled KG chains. We study chains of 512 and 1000 KG beads, corresponding to 8 and 15 entanglements, respectively. We compute the linear rheological properties from equilibrium simulations of the stress autocorrelation and obtain from these data the tube model parameters. Under nonlinear shear ﬂ ow, we compute the shear viscosity, the ﬁ rst and second normal stress differences, and chain contour length. For chains of 512 monomers, we obtain agreement with the results of Cao and Likhtman (ACS Macro. Lett. 2015, 4 , 1376). We also compare our nonlinear results with the Graham, Likhtman and Milner-McLeish(GLaMM) model.We iden-tify some systematic disagreement that becomes larger for the longer chains. We made a comparison of the transient shear stress maximum from our simulations, two nonlinear models and experiments on a wide range of melts and solutions, including polystyrene (PS), polybutadiene, and styrene – butadiene rubber. This comparison establishes that the PS melt data show markedly different behavior to all other melts and solutions and KG simulations reproduce the PS data more closely than either the GLaMM or Xie and Schweizer models. We discuss the performance of these models against the data and simulations. Finally, by impos-ing a rapid reversing ﬂ ow, we produce a method to extract the recoverable strain from MD simulations, valid for suf ﬁ ciently entangled monodisperse polymers. We explore how the resulting data can probe the melt state just before the reversing ﬂ ow. © 2019 The Authors. Journal of Polymer Science Part B: Polymer Physics published by Wiley Periodicals, Inc. J. Polym. Sci., Part B: Polym.Phys. 2019 , 57 , 1692 – 1704

INTRODUCTION Long-chain polymers molecules are at the heart of the multibillion pound plastics industry. They undergo large and rapid strains during processing, which strongly influence the properties of products. Entangled polymers exhibit complex viscoelastic behavior under flow due to involvement of a wide range of time and length scales, and topological constrains due to entanglements with surrounding chains. Understanding the dynamics of polymers under flow is central to controlling how polymers process and achieving desirable properties of end products. Relevant conditions to polymer processing are long chains under strongly nonlinear shear flow.
Theoretically, the dynamics of polymer chains are described by the tube theory, developed by Doi and Edwards 1 and based on de-Gennes' work. 2 This theory relies on the assumption that a single-chain model can be developed by treating the surrounding chains as a mean field that creates an effective tube around the test chain. Over time, the tube model has been improved to predict viscoelastic properties of entangled polymers by including important mechanisms such as contour length fluctuation, [3][4][5] constraint release, 6-10 and chain stretching. [11][12][13] An alternative theoretical picture for the dynamics of chain stretch has recently been presented by Xie and Schweizer,14,15 whose model includes a new interchain "grip force" that delays chain retraction until strands exceed a critical tension. Experimentally, the viscoelastic behavior of entangled polymers has been addressed for solutions at a range of concentrations [16][17][18][19][20][21] and for melts. 12,[22][23][24][25][26] Molecular dynamics (MD) is potentially an extraordinarily useful tool due to its enormously high spatial and temporal resolution. Specifically, it allows the user to track the individual and collective motion of chains and subchains in ways that are impossible by experiment. Moreover, it provides useful information about the relations between viscoelastic properties at the macroscopic level and the related molecular mechanisms at the microscopic level. 27 Recently, MD has been used to investigate the viscoelastic behavior of coarse-grained polyethylene chains, in which CX denotes a melt of linear chains containing X carbon atoms. This system has about 75 carbon atoms per entanglement. 28 Unentangled C50 and slightly entangled C178 chains have been simulated 29 in bulk and in confined systems. MD has been used to explore the structural, conformational, rheo-optical, and topological properties of an entangled C400 chains in the linear and highly nonlinear regimes 30,31 ; the dynamics and rheology of supercooled melts of C10 32 C20 and C150, 33 and C1000 28 chains; the steady 34 and transient 34 behavior of a C100 melt; the steady and transient 31,35 behavior of a C400 melt; and steady and startup shear of a C700 melt. 36,37 In polymer fluids, long-chain entanglement leads to the slow dynamics that makes polymer processing phenomena rich and difficult to predict empirically. These slow dynamics mean that MD of long polymers is extremely difficult without cutting-edge hardware and code. Consequently, a fast and highly universal polymer force field was developed by Kremer-Grest (KG), 38 which aims to describe a generic entangled melt rather than specific polymer chemistry. For this KG force field, we use NX to denote a melt of linear chains containing X beads. The KG model has been investigated 39,40 over a period of nearly 30 years but has only recently been applied to start-up of nonlinear shear of entangled melts. Initial results from Brownian dynamics simulations of mildly entangled N500 chains by Lu et al. 41,42 contradicted the tube model in several significant qualitative ways. Specifically, they reported chain stretching at flow rates that are expected to be too low. Furthermore, the results of Lu et al. 41,42 imply a violation of the stress-optical rule. Attempts to verify these results from other groups did not reproduce these findings. MD simulations of N200 chains by Masubuchi 45 There has also been some recent intense examination of the tube model's retraction mechanism, involving the interplay of experiments, simulation, and theory. This has focused on signatures of chain retraction predicted by the tube model that should be detectable by small-angle neutron scattering (SANS). Specifically, a nonmonotonic transient response in the radius of gyration perpendicular to flow, following a uniaxial extension. Recent SANS experiments by Wang et al. 46 did not observe this signature for a stretch of λ = 1.8 for polystyrene (PS) chains with Z = 34 entanglements, despite earlier SANS experiments confirming the signature for a λ = 1.7 stretch of Z = 54 polyisoprene chains. 47 MD simulations by Xu et al. 48 for λ = 1.8, Z = 34 KG chains confirmed the lack of retraction signature. However, the signature was observed in MD simulations of KG chains by Hsu and Kremer 49 when either the strain or the number of entanglements was increased. These simulations also showed that, during chain retraction, the leading anisotropic term of the single-chain structure factor is well predicted by the GLaMM model, particularly for moderately entangled chains. In summary, the GLaMM model predicts that the SANS signatures of retraction occur rather more sharply than is observed in PS experiments and MD simulations. Consequently, the GLaMM model predicts signatures of retraction for moderately entangled chains (Z = 34) at small strain (λ = 1.8) but these are not observed in PS experiments and KG simulations until longer chains 47,49 or larger strains 49 are employed. Finally, in the slip-link simulations of Masubuchi 50 the signatures are absent under milder flow conditions, in agreement with Refs, 46,48 even though chain retraction is present. As slip-link simulations are based on the tube model, but avoid some of the averaging approximations of the GLaMM model, this suggests that the overprediction of the retraction signature is due to mathematical approximations in the GLaMM model rather than a fundamental issue with the tube model. Indeed, Masubuchi tentatively ascribes the difference to entanglement fluctuations inducing inhomogeneities in the entanglement density along the chain contour.
In this paper, we aim to improve the use of MD simulations of KG chains to inform future developments of models of polymer dynamics. The above controversy indicates the challenging nature MD of simulations for nonlinear flows of entangled polymers and the need for careful verification of the underlying code base and of emerging simulation results. The studies of Masubuchi and Watanabe 43 and Cao and Likhtman 44 both used in-house codes, which are not publicly available. There is a clear need for a corresponding open-source code to widen access to this powerful technique, to benefit from more widespread scrutiny and testing of the underlying code, and to produce results for more strongly entangled chains. Recently, Schilling's group added the necessary shear boundary conditions 33 to the widely used MD code ESPResSo. 51 This enables nonlinear shear of polymers to be simulated within this welltested and highly parallel open-source software. In this study, we use this code to run nonequilibrium MD (NEMD) simulations for KG chains. We reproduce the results of Cao and Likhtman 44 for N512 chains and produce new data for chains of approximately double the number of entanglements. We also compare the peak stress from our simulations with a range of experiments and models and interpret the picture emerging from this comparison. Finally, we present a method to extract the recoverable strain from MD simulations on entangled polymers and test modeling ideas that link the recoverable strain to the polymer stress in the velocity gradient direction. 52 To assist with future work we have made examples of our input scripts available in the supporting information for this manuscript.

MODEL AND SIMULATION DETAILS
As in Cao and Likhtman, 44 we represent polymer chains by the conventional KG polymer model. 38,53 We used the standard KG model to equilibrate our systems, rather than using a soft nonbonded potential along with harmonic bonds as in Cao and Likhtman. 44 Then, we ran simulations under shear conditions and computed the transient shear viscosity using the following relation: Here η is the shear dependent viscosity, P xy is the xy component of the pressure tensor, and _ γ is the imposed system average shear rate. The pressure tensor for a molecular system can be written as On the right-hand side of eq 2, the first summation term is the kinetic contribution, and the second term is the potential contribution to the pressure tensor. Here i represents the particle index, α and β represent Cartesian components, m i is the mass of particle i, v iα , and v iβ are the velocity components of particle i in the α and β directions, respectively. In the second summation, r ijα represents the α component of the distance vector between particle i and j, and f ijβ is the β component of the force exerted on particle i by particle j.
In our NEMD simulations, we generated shear flow using Lees Edwards boundary conditions 54 with a dissipative particle dynamics thermostat. 55,56 All simulations were performed under constant volume and constant temperature conditions. We used a temperature of ϵ/k B , where ϵ is the Lennard-Jones interaction strength. We took the friction coefficient for the thermostat to be 0.5 τ −1 and the cut-off radius to be R c = 1.3σ, and m is the mass of the bead, k B is the Boltzmann constant, T is the temperature, and σ is the Lennard-Jones size of the beads. We used the ESPResSo 51 MD package for all simulations. We simulated systems consisting of 300 and 420 chains of lengths 512 and 1000, respectively. We ran simulations at shear rates _ γ ð Þ of 0.000005, 0.00001, and 0.00003 τ −1 . The integration timestep used in our simulations was 0.012 τ. We ran our equilibrium simulations on 8 and 12 nodes for N512 and N1000, respectively, and for the nonequilibrium simulations we used 5 and 7 nodes, for N512 and N1000, respectively. In all cases each node had 28 CPUs.

RESULTS AND DISCUSSION
We present results from nonequilibrium MD simulations. We first discuss the rheological properties of KG chains in equilibrium, then we focus on the transient behavior of these chains under nonlinear shear flow.

Linear Rheology
In Figure 1(a), we show the stress autocorrelation function G(t) against time. We computed G(t) using a recently developed efficient method. 57 We reproduce Cao and Likhtman's results 44 for N512 chains to confirm consistency of the two codes in linear response. The two N512 curves are in close agreement. We ran simulations for the much more entangled N1000 chains (solid squares in Fig. 1(a)), and are able to resolve well the terminal time. This required a long simulation involving 9 × 10 9 MD simulation time steps, which corresponds to 1.08 × 10 8 τ or approximately 16 terminal times, if we use double reptation to approximate the terminal time as half the Likhtman-McLeish 5 reptation time (see Table 1). Our simulation parameters are compared with those of Cao and Likhtman in Table 1.
We computed the storage G 0 (ω) and loss G 00 (ω) moduli by fitting G(t) with Maxwell modes, as in Likhtman and Cao. 44 Tables 2 and 3 show the resulting Maxwell mode parameters and Figure 1(b) plots G 0 (ω) and G 00 (ω). We then fitted these curves with the Likhtman and McLeish model, 5 which is a quantitative tube model for linear response. For this fitting, we fixed the constraint release parameter at c ν = 0.1 and fitted the remaining three tube model parameters, namely, the (a) (b) entanglement modulus, G e , the Rouse relaxation time of an entangled segment, τ e , and the number of monomers per entanglement, N e , which should be a constant for a given chemistry and volume fraction. We fitted our N512 and N1000 data simultaneously with the same parameter values.
Our resulting parameters are shown in Table 1 and compared with Cao and Likhtman's values, obtained by fitting their N350 and N512 data simultaneously. Table 1 also contains the number of entanglements, Z = N/N e , the chain Rouse time τ R = Z 2 τ e and the reptation time (computed using eq 12 from reference 5 ). We expect our parameters to be superior to Cao and Likhtman's because we fit to a wider range of molecular weights (almost a factor of 2 in our case), we exclude the weakly entangled N350 data and because our set includes a well-entangled system (N1000). Table 1 shows that we obtain the same values for τ e and G e as Cao and Likhtman. This is expected as these parameters can be determined accurately from moderately entangled chains and this requires only high frequency data, which can be obtained with good statistical resolution from a comparatively short simulation. That is, adding a highly entangled sample does not update these parameters from Cao and Likhtman's results. In contrast, we see a noticeable change in N e , because this value is sensitive to the terminal region and is best determined from data on well-entangled chains, as we provide here. The value of N e is important for accurate tube model predictions for fast nonlinear flow as N e is required to compute the overall chain Rouse time, τ R . Indeed, Table 1 shows that we determine a somewhat smaller τ R for N512 than Cao and Likhtman. Later, this improved τ R will be important when comparing our MD results to experiments. Figure 1(b) shows that the Likhtman and McLeish model generally predicts G 0 (ω) and G 00 (ω) very well but that we have enough statistical resolution in our simulation data to expose a systematic disagreement for the shape of the peak in G 00 (ω).

Nonlinear Shear Stresses
We ran MD simulations for start-up of nonlinear shear at a range of shear rates. We computed the transient shear stress using eq 1 and we used logarithmic bins in time for averaging as described by Cao in his thesis. 58 In this technique, the The tube model parameters were obtained by fitting the Likhtman-McLeish model 5 to the storage and loss modulus data in Figure 1(b).  quantity C(t) at any time (t) after start of shear can be computed using the following relation: Here t i are averaging time intervals and t i = t 0 + P i−1 k = 0 t 0 k , t 0 is the shear start time, t 0 i = DM i , D is constant for the first time interval, and M is a multiplication factor. In this work, D and M are set to 8Δt and 1.1, respectively, where Δt is the simulation time step.
Our results for N512 are shown in Figure 2, along with a comparison of the results from Cao and Likhtman 44 and this work, using open and solid symbols, respectively. We do not see any significant difference between our results and Cao and Likhtman's results. 44 At these shear rates, the shear stress exhibits a transient overshoot. Also included in Figure 2 is the linear viscoelastic envelope, η 0 (t), which we calculated from using the Maxwell modes we obtained by fitting our linear simulations of G(t) as in Figure 1 and Table 2. Both sets of nonlinear simulation data conform to η 0 (t) at low strain. Finally, we also show the predictions of the GLaMM model 13 using the tube model parameters obtained from our linear simulations (Table 1) and including the minor modifications proposed by Auhl et al. 59 Our revised tube model parameters are important to the comparison of our MD to experiments because our improved value of N e reduces the uncertainty in the Rouse time of the simulations. The GLaMM model shows a systematic overprediction of the time and height of the shear stress overshoot relative to the simulation data.
In Figure 2, we show the shear viscosity as a function of time at different shear rates for the longer N1000 chains. As with Figure 3, we also include the linear viscoelastic envelope and the predictions of the GLaMM model using the parameters from linear response ( Table 1). The systematic disagreement is somewhat more pronounced for these longer chains. The model overpredicts the peak and steady-state stress, with the effect becoming stronger at higher rates. Some of this greater disagreement, relative to N512, can be explained by the greater Rouse Weissenberg numbers for N1000.
In Figure 4, we show the first normal stress difference coefficient as a function of time for N1000. The first normal stress difference overshoots, like the shear viscosity, but the peak value is shifted to later times. It was also observed by Jeong et al. 35 that the peak in the first normal stress difference occurs at t ψ max The dashed line in Figure 5 is the linear viscoelastic envelope, obtained from which describes the nonlinear data at small strains. Also shown are the GLaMM model predictions, which are somewhat above the MD data for the highest shear rate. Figure 5 shows the transient second normal stress difference, N 2 , for N1000 chains, plotted as −N 2 . We observe an undershoot in N 2 , manifest as a maximum in the negative N 2 .The minimum at In addition to the undershoot, our simulated N 2 is small, negative, and shear thinning, all of which are seen in experiments. 60 Despite its small magnitude the second normal stress difference is significant as it has been linked to the onset of edge fracture. 60,61 The GLaMM model, like many constitutive models, predicts a zero second normal stress difference under shear.

Transient Chain Contour Under Shear
In order to examine the chain deformation directly from our MD simulations, we computed the chain contour length, as follows. We divided each chain into Z = N/N e subchains, with each segment having 65 (N e ) monomers. We computed the square end-to-end vector of each subchain R i and took the ensemble average. We then computed the contour length of the whole chain by summing over the square root of each mean square subchain length. That is, We chose this method of averaging as it is a measure of the chain contour length that is readily computed from the GLaMM model (see eq 20 of Ref. 13 ). In Fig. 6, we show the ratio between the average contour length L and the equilibrium length L 0 against time for N512 chains under shear, computed from both MD and the GLaMM model. A corresponding plot for N1000 chain is shown in Figure 7. The worst agreement with the GLaMM model occurs where the MD stresses are overpredicted, namely, in steady state for both N512 and N1000 and at the highest shear rate for N1000. This shows that an overprediction of the chain stretching in the GLaMM model is at least partially responsible for the overprediction of the stresses in MD. Generally, the stretch evolution in the simulations occurs later than the GLaMM model predictions, which likely contributes to the GLaMM model predicting stress peaks that are earlier than those in the simulations.
There has been some discussion in the literature of how reliable coarse-grained forcefields are for high flow rates. Jeong et al. 35 concluded that the chain orientation is the primary cause for the macroscopic stress overshoot of entangled polymeric systems under start-up shear. Baig and Harmandaris 62 reported larger chain deformation in their course-grained simulations of PS melts than their atomistic simulations, for Weissenberg numbers above 5. We confirmed that our results  show no unphysically large stretching. Finally, we will show in Section 3.4 that our simulated peak stresses follow the same scaling law with flow rate as experiments for all shear rates in our study, suggesting no onset of unphysically large deformation in our simulations.

The Shear Stress Maximum and Comparison with Experiments
The shear stress overshoot is a key transient feature at high shear rates. To allow a comparison between our simulations, experiments, and models, we have plotted, in Figure 8, the features of the shear stress maximum for varying shear rate for several difference chemical species. These include PS, polybutadiene (PBD), and styrene-butadiene rubber (SBR), encompassing both melts and entangled solutions. These plots require τ R and the entanglement modulus, G e , for each experimental system. For some materials, 18 these values have previously been determined 13 from linear rheology using the Likhtman and McLeish model. 5 For all other materials τ R was reported in the original experimental papers and we have used these values. In addition, where the original experimental papers report a plateau modulus, G N , we convert this value to G e by multiplying by the standard factor 63 of 5/4. We have also applied this 5/4 factor as a downward shift to the Xie and Schweizer model in Figures 8(b-d) The data in Figure 8 fall into three distinct groups: low Z, high Z, and the PS data. To illustrate this behavior, we separate the peak strain plots into two parts, low Z and PS data in Figure 8 (a) and high Z data in Figure 8(b). We employ a similar separation for the peak stress data in Figure 8(d,e). The PS curves are clearly different to all the other materials and this distinct behavior for PS is somewhat unexpected. We used PS data from two different groups 24,25 and confirmed that both lie on the same mastercurve. We also confirmed the values of τ R by fitting the Likhtman and McLeish model to the linear rheological experiments and obtained essentially the same τ R as reported in the original experimental papers. Thus, the separate behavior for PS is clearly established. Interestingly, our MD simulations follow the PS mastercurve and not that of the other materials.
The peak strain plots in Figures 8(a,b) show three types of behavior: all of the PS data collapse to a single master curve, which is captured well by the Xie and Schweizer model; all of the remaining low Z data are well described by the GLaMM model for _ γτ R ≤ 10; and the high Z data are also captured by the GLaMM model for _ γτ R ≤ 10, but with some disagreement for some of the lowest Weissenberg number measurements. Specifically, at very low Rouse Weissenberg numbers ( _ γτ R < 0:3), the GLaMM model fails to capture the SBR melt with Z = 76, which is the most entangled sample for which there is low _ γτ R data. Here, the model slightly overpredicts the maximum strain. This flow regime is expected to be dominated by convective-constraint release (CCR), and so this may be experimental evidence of some deviation from the CCR implementation in the GLaMM model. Finally, we note that the high Z experiments and the simulations show γ max / _ γ 1=3 scaling, which is captured by the Xie and Schweizer model but not the GLaMM model. Hence, the Xie and Schweizer model captures the shape of these data but gives numerical values that are somewhat too large.
Further details of the shear stress maximum are revealed by plots of the shear stress maximum. Figure 8(c), which plots the maximum stress against maximum strain, again shows two mastercurves: the PS experiments collapse to one curve and the non-PS experiments collapse to a separate curve. The collapse is very clear at moderate and high γ max but somewhat less clean at lower γ max , where all of the low Z data show some dependence on degree of entanglement, Z. Nevertheless, the PS mastercurve is noticeably below the non-PS mastercurve. As with the previous plots ( Figs. 8(a,b)), the GLaMM model captures well the non-PS experiments, for _ γτ R < 10. In contrast, the Xie and Schweizer model predicts a σ max xy that is above all of the experiments. In particular, this model deviates mostly strongly from the PS data, despite successfully predicting γ max for PS. Our MD simulations capture both γ max and σ max xy for the PS materials. A similar picture is apparent in Figs 8(d,e), namely, that the GLaMM model captures the non-PS data for _ γτ R < 10; the Xie and Schweizer model is higher than all data; and the MD simulations are close to the PS data. However, again the Xie and Schweizer model has the correct shape over a wide range of flow rates for the high Z data, but with predicted value of a   _ γτ R > 1, similarly to previous studies. 44 Future MD studies exploring lower Weissenberg numbers, although computationally expensive, might elucidate the physics behind this lower flow rate behavior in the PS experiments.
We now consider reasons for the divergence between the GLaMM model and some of the experiments. The GLaMM model assumes Gaussian chain statistics and constant monomer friction under flow. Both of these assumptions are most accurate for materials where N eK , the number of Kuhn steps between entanglements, is high. Thus we expect the best agreement for highly diluted, entangled solutions. Our data comparison shows that the SBR melts conform to the same master curve as the solutions and so appear to have a sufficiently large N eK to agree with the GLaMM model. Auhl et al. 59 have also shown good agreement with polyisoprene melts, which is consistent with the comparatively large melt N eK for this chemical species. Regarding Z, early versions of the tube model generally neglected some mechanisms that are important for shorter chains, such as constraint release and contour length fluctuations and so were expected to work best for large Z. However, the GLaMM model includes these corrections so can predict melt and solution data as low as Z = 7. 59,64 Below this threshold linear tube models begin to breakdown, as discussed in Ref. 59.
Validation of the GLaMM model for very large Z (Z > 50) is less comprehensive due to the rarity of nonlinear data for such highly entangled chains. Figure 8 includes the most entangled samples for which there is available nonlinear shear data. The agreement of the non-PS data is equally good for all values of Z, provided _ γτ R < 10, apart from a small number of minor deviations detailed above.
Some recent modeling work 65 suggests a possible reason for the different behavior of PS compared to other materials, N eK , which is lower for PS than for other melts. 66 It is expected that the monomer friction will reduce with local alignment of monomers 65 and this will become important for nonlinear flows of low N eK fluids. Hence PS is more susceptible to this effect than other melts. Experimental evidence on PS melts qualitatively supports this. Specifically, Ianniruberto et al. 65 observed that modeling alignment-dependent friction is more important to capture PS extensional rheology, compared to melts of other chemistries, due to its small N eK . Furthermore, Robertson et al. 67 showed that including alignment-dependent friction in the tube model produces much-improved predictions of extrudate swell in PS, particularly for near-monodisperse polymers. In solutions reduced friction is expected to be absent because the local environment is dominated by the solvent. 68 Hence, because the non-PS melts and the solutions in Figure 8 shared a single master curve, alignmentdependent friction seems negligible here. Furthermore, this conclusion is also supported by the agreement of these experiments with the GLaMM model (for _ γτ R < 10), which neglects changes in friction under flow.
It is less clear why our KG chain simulations are close to the PS experiments, rather than the other materials. Our KG chains have approximately two beads per Kuhn step, 69 meaning N eK = N e /2 = 32.5 which is small, but not as small as PS. There may also be a role for bead density and temperature in KG simulations. The standard choices for these, taken herein, give a relatively high density compared to the Lennard-Jones interaction distance (ρ = 0.85σ 3 ) and a low temperature relative to the interaction energy (T = ϵ/k B ). Both of these factors mean that monomer-monomer interactions will involve significant excluded volume in addition to noncrossability of chains. Hence these factors may increase the role of alignmentdependent friction relative to non-PS melts and solutions. There may also be a role for fluctuations in the entanglement density, 50 absent in the GLaMM model, to explain the PS and MD data.
We now examine support from our simulations for the intermolecular gripping force in the model of Xie and Schweizer. 14,15 The Xie and Schweizer model predicts well the peak strain from our simulations (see Fig. 8a) as well as the corresponding data for PS melts. In particular, the model has the correct scaling exponent, as reported in Refs 14,15. The model achieves this because the intermolecular gripping force delays chain retraction relative to the GLaMM model for moderate values of _ γτ R by requiring strands to exceed a critical tension. Consequently, the Xie and Schweizer model predicts a noticeably higher σ max xy than the MD and PS data in Figure 8(b, c). Here, the MD and PS data fall below the non-PS data and the GLaMM model, suggesting less chain stretching in these data. This picture of suppressed stretching in the MD data, relative to the GLaMM model, is also observed directly in our simulation results for transient chain stretching in Figures 6  and 7.

Strain Recovery from a Rapid Reversing Shear Flow
There has been recent experimental 70,71 and theoretical 52,72 interest in recoverable strain in polymers. Experiments by Lee et al. 71 have linked recoverable strain to microscopic structure and normal stress in two very different polymeric systems. In entangled polymers, the bulk of the strain is typically recovered in a very fast, time-dependent, reversing shear flow, as the slow-relaxing polymeric stress far exceeds the viscous stresses. 52,72 In such a fast reversing flow, the stress-strain behavior is independent of the strain rate as there is insignificant relaxation during the strain recovery. Thus we model the recovery period in our simulations as a rapid, constant rate reversing shear flow. The recoverable strain corresponds to the strain at which the shear stress returns to zero. Our MD simulations to investigate recoverable strain are shown in Figure 9. We picked a point from the shear stress transient in Figure 2, corresponding to an imposed strain of γ imposed and applied a reversing shear at a rate of _ γ r . The shear stress transients for these reversing flows are shown in the insets of Figure 9, in which γ rec is the shear strain applied in the reversing direction and σ 0 xy is the shear stress at the point the shear is reversed. For γ imposed = 16, we imposed three different reversing shear rates in the range j _ γ r j τ R = 23−69. Figure 9 shows that, as expected, the stress-strain behavior is independent of the flow rate. We repeated the test for two rates for γ imposed = 9, obtaining the same collapse of the stress-strain results. Thus using this shear rate window, we could determine the recoverable strain from the reversing strain required to return the shear stress to zero. Figure 9 shows the recoverable strain extracted from our simulations. In these simulations, the imposed strain is almost entirely recovered for imposed strains up to the shear stress maximum, with the recovery dropping off rapidly beyond γ imposed = γ max as in experiments (see Figure 10). 70 It is clear from the insets of Figure 9 that the reversing shear stress transients show considerable curvature. In particular, there is a rapid initial drop-off in the shear stress shortly after the reversing flow is applied. To examine this more closely, we plotted in Figure 11(a) the reversing viscosity, defined as where t = 0 corresponds to the start of the reversing flow. Our simulation results lie between η 0 (t) (the linear viscoelastic envelope for start-up shear, using the spectrum from Table 3) and η fast (t), the same spectrum but only including modes for which τ j _ γ rev j < 1. Between these two limiting curves the MD results follow an apparent power law with the exponent weakening with increasing γ imposed . In contrast, for sufficiently entangled polymers, the recoverable shear modeling of Holroyd et al. 52 predicts η rev (t) = G eff t, where G eff is an effective modulus. The apparent power-law behavior in our simulation results is likely due to the linear viscoelastic effect of fast relaxing modes that are nonnegligible as the chains are only moderately entangled. To correct for this we subtracted η fast (t) from η rev (t) in Figure 11(b).
This corrected quantity has the expected linear scaling with t and from this we could extract a G eff . This suggests that the response of the fast modes to the change in shear rate is the same as for a start-up flow, whereas the slow modes retain a memory of the earlier flow. Our results for G eff are shown in Figure 12(a), normalized by the effective modulus from the linear viscoelastic envelope, G 0 . The recoverable shear modeling of Holroyd et al. 52 predicts G eff /G 0 = σ yy /G N , where σ 0 yy is the polymer stress in the yy direction at the point that the shear is reversed and G N = 4/5G e is the plateau modulus. The GLaMM model predictions for σ yy /G N are shown in Figure 11(b). They agree well with the G eff /G 0 values from our MD simulations. The key area of disagreement is that the GLaMM model predicts a small rise in σ yy at large strains, as σ yy undershoots, but this is not reflected in G eff . Our results suggest an experimental protocol, involving only shear stress measurements, that could provide information about the polymer stress in the velocity gradient direction.

CONCLUSIONS
We performed MD simulations to investigate the nonlinear behavior of entangled KG chains. We report results for two chain lengths of 512 and 1000 beads, which have 8 and 15 entanglements, respectively. These results provide benchmarking for the weakly entangled system and novel data for a more strongly entangled system. We also present novel data for the normal stresses and transient chain stretching for both chain lengths. For both chain lengths we computed, from our simulations, the linear stress relaxation function G(t) and, under nonlinear shear, we computed the transient shear stress and chain contour length. For N1000 we also generated MD data for the first and second normal stress differences. Our shear stress data for N512 are entirely consistent with those of Cao and Likhtman, 44 providing an independent verification of their results using a different thermostat. This is significant given the contradictory results from Lu et al. 41,42 Our stress results for N1000 have many qualitative features in common with experiments. We observe an overshoot in the shear viscosity and the first normal stress difference, with the peak in the first normal stress difference coefficient, occurring at t ψ max The second normal stress difference is small, negative, shear thinning, and shows an undershoot, occurring at t ψ min 2 À Á ≈t η max ð Þ. We ran our simulations within the widely used MD code ESPResSo, 51 well-tested and highly parallel open-source software. The multinode parallel functionality of this code was essential to access the very long simulation times that are needed for highly entangled chains. Furthermore, as ESPResSo is open source, it is available for use, scrutiny and modification by the whole research community.
We analyzed our MD data using standard linear 5 and nonlinear 13 implementations of the tube model. For the linear rheological simulations, we build on the results of Cao and Likhtman 44 by adding a significantly more entangled chain (N1000). This wider range of molecular weights enables more effective fitting with  Table 3 and η fast (t) is the same spectrum but only with modes for which τ j _ γ rev j < 1. the Likhtman-McLeish 5 model, which produces an improved value for the number of monomers per entanglement, N e , for KG chains. This improved parameterisation leads to improved values for the Rouse time of KG chains, which is a key element in comparing these simulations to experimental data. Many of the features in our simulated stresses are qualitatively similar to experiments. Indeed, our comparison with experiments shows that our KG results capture well PS melt experiments. The N 2 data that we present may be useful in informing new constitutive models, particularly those aimed at capturing edge fracture and similar instabilities.
We made a systematic comparison of the transient shear stress maximum for our simulations, experiments on a wide range of both melts and solutions, including PS, SBR, and PBD. We also compared with predictions from the GLaMM and Xie and Schweizer models. This comparison establishes that PS data and non-PS data show markedly different behavior and that KG chains reproduce the PS data more closely than the GLaMM or Xie and Schweizer models. We ascribe the PS behavior to the low N eK value of PS melts and KG chains. The model of Xie and Schweizer 14,15 captures the strain of peak stress for PS melts and KG simulations, but overpredicts the stress peak in experiments and simulations. However, the Xie and Schweizer model has the correct shape over a wide range of flow rates for the high Z data despite its overprediction relative to these experiments. The GLaMM model 13 captures both the peak stress and strain for the non-PS experiments for shear rates up the _ γτ R1 0, beyond which the peak height and strain are over predicted. The GLaMM model overpredicts the chain stretch in KG chains and this is a likely reason for its difficulty in predicting the features of the shear stress maximum in PS melts. The Xie-Schweizer model also overpredicts these stretch data. We identified alignment dependent monomer friction as a potential low N eK correction for the GLaMM model. There may also be a contribution from inhomogeneities in the entanglement density along the chain contour.
We developed a method to extract the recoverable strain from MD simulations, via a rapid reversing flow, which is valid for sufficiently entangled monodisperse polymers. We generated recoverable strain data under rapid reversing flow for N1000 chains. We also analyzed the reversing flow stress transients and identified fast relaxing parts of the melt spectrum that behave as though fully relaxed when the flow is reversed and a slow relaxing portion of the spectrum that retains a memory of the melt state just before the reverse flow. From this separation of the reversing flow transient we could extract from our MD data, an effective modulus G eff for the reversing flow. Our data for G eff show that the initial forward shear reduces the modulus that opposes the reversing flow, an effect that increases with the strain imposed during the forward shear. The modeling by Holroyd et al. 52 predicts that G eff / G 0 = σ yy G N and we verify this prediction by showing that the GLaMM model captures the shear induced changes in G eff . Our reversing flow protocol requires only shear stress data from a rapid reversing flow and the linear viscoelastic spectrum and so should be achievable in experiments, to provide an experimental method to probe the polymer stress in the velocity gradient direction from shear stress measurements alone. Indeed, recent experiments by Lee et al. 71 also linked recoverable strain to normal stress measurements.