Clustering of Entanglement Points in Highly Strained Polymer Melts

Polymer melts undergoing large deformation by elongation are studied by molecular dynamics simulations of bead–spring chains in melts. By applying a primitive path analysis to strongly deformed polymer melts, the role of topological constraints in highly entangled polymer melts is investigated and quantified. We show that the overall, large scale conformations of the primitive paths (PPs) of stretched chains follow affine deformation while the number and the distribution of entanglement points along the PPs do not. Right after deformation, PPs of chains retract in both directions parallel and perpendicular to the elongation. Upon further relaxation we observe a long-lived clustering of entanglement points. Together with the delayed relaxation time this leads to a metastable inhomogeneous distribution of topological constraints in the melts.


■ INTRODUCTION
Materials based on polymeric compounds are essential for many areas of modern technology. Their outstanding properties originate from chemical details determining local interactions and the fact that polymers typically are long chain molecules. The latter, generic aspect is of interest for the present work, where the fact that chains cannot cross through each other plays the central role. Such complex topological constraints resulting in entanglements play an essential role for dynamical and rheological properties of polymer melts. 1−8 In the linear viscoelastic regime, these are well described by reptation theory (the tube model), 1,2,9 and have been confirmed by many simulations 10−21 and experiments. 22−26 However, despite this remarkable achievement, a precise definition of an entanglement within the reptation concept still is lacking, and attempts to include multichain effects analytically were of very limited success only. 27 Everaers et al. 4 introduced the primitive path analysis (PPA) based on the concept of Edwards' tube model 28 to identify the backbone of the tube, i.e., the primitive path (PP) of each polymer chain in a melt, and applied it to bead−spring chains. 11 A detailed discussion regarding self-entanglements, local self-knot effect, and finite-size effect is given in refs 5 and 29. Kroger et al. 30,31 subsequently developed the Z-code and its updated version, Z1-code, where PPs are treated as infinitely thin and tensionless lines using geometrical operations rather than multibead chains. Once all PPs, represented by the shortest (optimal) paths (SPs), are obtained, the entanglement molecular weight N e is simply obtained from the length of the path. Alternatively one can also count the kinks of adjacent segments along the chain to get a first rough estimate. Tzoumanekas et al. 32 implemented another algorithm, CReTA, which is capable of reducing the atomistic configuration of a computational polymer sample to a network of corresponding PPs where the topological constraints are conserved. Other similar methods dealing with entanglements in polymer melts are given in refs 33−35. Thus, there are a number of methods available, which analyze the role of topological constraints, i.e., the role of entanglements in a polymer melt. Some of them, such as PPA, are reproducing the entanglement length quantitatively correctly. 4,20,36 Beyond the analysis of constraints themselves, there are many discussions in the literature 8 dealing with the motion of polymers due to entanglement constraints and the time such constraints last. For entangled chains, processes of contour length fluctuation (CLF) 37−39 of the PPs and constraint release (CR) 38,40−42 also contribute in addition to pure reptation, i.e., the motion of the polymer along the hypothetical reptation tube, and have to be considered. However, it is not yet clear in which way all the abovementioned concepts can be employed for describing the behavior of polymer melts in the nonlinear viscoelastic regime. To shed light on this problem, we recently have started to approach this problem from the computational side. 43,44 We have compared the predictions of chain conformations given by the Doi−Edwards tube model 2 and its extensions based on the Graham−Likhtman−McLeish−Milner (GLaMM) tube model 45 to extensive simulations of polymer melts in the nonlinear viscoelastic regime. For this we decided to concentrate on the isochoric elongation of polymer melts. The chain retraction mechanism, as predicted by the GLaMM concept, which sets in right after deformation, has been investigated with contradicting results, based on both experimental and simulation data 44,46−50 while our recent results for longer chains support this general scheme. 44 The GLaMM model also correctly accounts for stress undershoots in transient shear viscosity. 45,51,52 Although primitive path network models can account for viscosity changes in entangled polymers upon elongational and shear flow 53−56 in the linear viscoelastic regime, they seem to fail for strain rates ε̇faster than the inverse Rouse time τ R,N of the whole chains. For example, the monotonic elongational thinning behavior of entangled polymer melts in elongational flow observed from the experiment 57 at ε̇> 1/τ R,N cannot be described by simulations using the primitive path network models. 55 The present work intends to make a contribution to filling a gap in this research field.
Here we start from well-equilibrated and highly entangled polymer melts composed of weakly semiflexible bead−spring chains at a monomer density ρ = 0.85σ −3 , prepared by a new, efficient hierarchical methodology. 20,29,58 These melts are subject to strong deformation by isochoric elongation in the nonlinear rheological regime. Following this deformation, we investigate in detail the subsequent relaxation. Applying the primitive path analysis (PPA) 4,5 to strongly deformed polymer melts in the nonlinear viscoelastic regime, we recently have shown that the force pattern along the primitive paths (PPs) qualitatively matches that of the corresponding original paths (OPs). 43 This indicates that the conformations of OPs of chains within fuzzy tubelike regimes are well represented by their corresponding PPs. We directly can relate sign switches of the tension force to kinks, "effective entanglement points", of high curvature along the PPs. Based on these findings, we use the relaxation of PPs of deformed chains in a melt in order to shed some light on the role of topological constraints for the relaxation of highly entangled deformed polymer melts.
The outline of this paper is as follows: in the next section, we summarize the main features of our model and the simulation techniques. In the third section we describe the conformational changes, the characteristics of topological constraints, and the stress relaxation of the deformed polymer melts followed by our conclusions in section four.

■ MODEL AND SIMULATION METHODS
Melts of Bead−Spring Chains with a Weak Bending Stiffness. For our simulations, a polymer melt consisting of n c polymer chains of chain size N, i.e., the number of monomers, is described by a standard bead−spring model 11 at a monomer density ρ = 0.85σ −3 = n c N/V, where σ = 1 is the unit of length and the size of a monomer. V = L x L y L z is the volume of the simulation box with three orthogonal linear dimensions, L x , L y , and L z . Any pair of bonded and nonbonded monomers located at a distance r apart interact via a shifted, purely repulsive Lennard-Jones (LJ) potential 59−62 U LJ (r) where ϵ is the energy unit of the pairwise interaction and r cut = 2 1/6 σ is the cutoff in the minimum of the potential such that force and potential are zero at r cut . Any pair of bonded monomers interacts via the finitely extensible nonlinear elastic (FENE) binding potential 59 U FENE (r) where k = 30ϵ/σ 2 is the force constant and R 0 = 1.5σ is the maximum value of bond length. These Lennard-Jones units σ and ϵ also provide a natural time definition via τ σ = ϵ m/ where m = 1 is the mass of the particles. In addition, a weak bond-bending potential 4 U BEND (θ) with a chain stiffness parameter k θ is introduced where θ is the angle between two subsequent bonds, i.e., , and b j = r j − r j−1 is the bond vector between monomers j and j−1 along the chain. Choosing k θ = 1.5ϵ, where chains become weekly semiflexible, the mean-square radius of gyration for unperturbed (i.e., fully equilibrated) chains in a melt is where ⟨R e 2 ⟩ 0 is the mean square end-toend distance and l b = ⟨b 2 ⟩ 0 1/2 ≈ 0.964σ is the root-mean-square (rms) bond length. This gives ⟨R e 2 ⟩ 0 , ⟨R g 2 ⟩ 0 ∝ N 2ν with the Flory exponent ν = 1/2. Here the average ⟨···⟩ 0 denotes the average over all chains and over all independent configurations of unperturbed melts. Polymer chains in a melt behave nearly as Gaussian chains. 20 For the above parameters the corresponding entanglement length N e = N e,PPA (0) ≈ 28, estimated both from the plateau modulus G N 0 = (4/5)ρk B T/N e as well as from the primitive path analysis (PPA). 4,5,20,29 The Rouse relaxation time of a subchain of entanglement length N e and of the overall chain of size N is τ e = τ 0 N e 2 and τ R,N = τ 0 N 2 , respectively. Here the characteristic time prefactor τ 0 ≈ 2.89τ is determined from the estimate of the mean-square displacement of inner monomers 20 for n c = 1000 chains containing 18 ≤ Z ≡ N/N e ≤ 72 entanglements, g 1 (t). We here focus on the above-mentioned cases and use the molecular dynamics simulations package 63 ESPResSo++ for all runs with simulation time step set to Δt = 0.01τ. The temperature is kept constant (T = 1ϵ/k B , k B being the Boltzmann factor) through a Langevin thermostat with a weak friction constant Γ = 0.5τ −1 .
Deformation Mechanism: Isochoric Elongation. We start from fully equilibrated, highly entangled polymer melts (unperturbed polymer melts) originally in a cubic simulation box, i.e., L x = L y = L z = L 0 , with periodic boundary conditions along the three orthogonal directions. We apply a simple "isochoric elongation" deformation mechanism. At each elongation step the whole simulation box is instantaneously stretched by a factor of 1.02 along the x-direction, contracted in the y-, z-directions by a factor of 1/ 1.02 such that the box size V = L x L y L z = L 0 3 is kept as a constant. This deformation step is so small that it does not induce any instabilities in the simulation. Then the system is given a short time to relax. By that an average fixed strain rate ε̇in the range τ R,N −1 <ε̇< τ e −1 is introduced, namely i.e., L = L 0 exp(εṫ). We set the strain rate ε̇to ετ R,N = 77, that is, ετ e = 77(N e /N) 2 (e.g., ετ e ≈ 0.015 for N = 2000) at each elongation step.
To obtain this deformation rate, the system can relax for (0.02τ R,N / 77)τ between two such steps. Thus, one could expect that locally, i.e., below a few N e ≈ ≈ N N ( 1/0.015 8.2 e e ), the chain can fully relax during the elongation, while globally, the chain follows the affine deformation implying that the macroscopic chain deformation follows the macroscopic strain. Altogether 81 elongation steps are performed, leading to a total strain of 43,44 λ ≈ 5.0. As a result the deformed polymer melts are deep in the nonlinear viscoelastic regime.
Primitive Path Analysis (PPA) and Definition of Significant Kinks. According to the original PPA procedure, 4,5 the two ends of all the chains in the polymer melt are fixed at their actual position in space. Then the intrachain excluded volume as well as the bondbending interaction are switched off, i.e., U LJ (intra) (r) = U BEND (θ) = 0, while interchain excluded volume interactions remain in order to prevent bond crossing and to preserve interchain topological constraints. Intrachain topological constraints were shown to be insignificant 5 within the error bars achievable here. Then the temperature is set to zero, so that the chains contract to the shortest paths between the two ends, observing all interchain topological constraints. Practically, the temperature is set to T = 0.001ϵ/k B (close to zero), and the basic time step Δt is reduced to 0.006τ. The friction constant is set to Γ = 20τ −1 during the first 10 3 MD steps, and Γ = 0.5τ −1 afterward. 5,20,29 Thus, chains straighten out when the bond springs try to reduce the average bond length in order to minimize the energy from l b ≈ 0.964σ of the OPs to ⟨b PP ⟩ 0 = 0.31σ of the PPs for unperturbed polymer melts and ⟨b PP ⟩ λ ≈ 0.57σ for strongly deformed polymer melts at λ ≈ 5.0, respectively.
A typical snapshot of the PP of one selected chain of size N = 2000 in a deformed melt is shown in Figure 1a. The PP consists of straight pieces with relatively sharp kinks created by the excluded volume interactions with other chains. The distribution of kinks along the selected chain strongly depends on the surrounding chains (see Figure   Figure 1. ) 2 ⟩ 0 ] 1/2 for the PPs of chains (b), plotted versus the strain λ. Rescaled mean-square internal distance, C⟨R 2 (n)⟩ λ /⟨R 2 (n)⟩ 0 for the OPs of chains (c) and C⟨(R (PP) (n)) 2 ⟩ λ /⟨R (PP) (n) 2 ⟩ 0 for the PPs of chains (d), plotted versus the rescaled chemical distance n/N e . Two components in the directions perpendicular (⊥) and parallel (||) to the direction of stretching are shown, as indicated. Data are for elongated polymer melts of chain sizes N = 500, 1000, and 2000, as indicated in (a), (b), but only of N = 2000 in (c), (d). In (a), (b), the expected scaling laws for affine deformation are shown by straight lines. In (c), (d), five strain values of λ are chosen, as indicated, and n = 8.2N e is pointed out by an arrow (cf. text). Data for the elongated PP mesh are also included in (d) by a red curve.

Macromolecules
Article 1b where only short segments of surrounding chains near the test chain are shown). Therefore, a direct way of recognizing "entanglement points" is to analyze the curvature along the PP to identify these kinks. Since the sharp kinks are still rounded off due to the resulting short bonds and the remaining interchain excluded volume, we have chosen the bond angle θ j,j+5 between bonds b j and b j+5 for j = 2, 3, ..., N − 5 along the chain as shown in Figure 1c. For identifying significant kinks (entanglement points), all local maxima marked by symbols along the path denoted by j are sorted into four categories, θ j,j+5 > 90°, 90°≥ θ j,j+5 > 60°, 60°≥ θ j,j+5 > 30°, and 30°≥ θ j,j+5 > 15° (Figure 1d). Comparing to the force pattern along the same PP, we define an entanglement point to be located at r j+2 if the curvature of a kink corresponding to 43 θ j,j+5 ≥ 60°. There is, however, an ambiguity related to this analysis. While the number of topological constraints up to some end effects is strictly conserved, their physically relevant number will vary with time. When two kinks come very close along the backbone of the constraining chain, they probably act like one. Because of that we always present results, where we count them as one, if the distance along the same chain is less than 1σ.

■ SIMULATION RESULTS
We study three polymer systems each containing n c = 1000 chains of sizes N = 500, 1000, and 2000, and the corresponding lengths of the unperturbed simulation box being L 0 /σ = 83.79, 105.57, 133.01, respectively. All estimates of physical quantities are taken as averages over all n c = 1000 chains except if a selected chain in a melt is explicitly mentioned and discussed. All n c = 1000 chains are labeled by i = 1, 2, ..., 1000. For the fully equilibrated melts, chain conformations are characterized by the mean square radii of gyration ⟨R g 2 (N)⟩ 0 /σ 2 of approximately 221, 438, and 909 for the OPs of chains, ⟨(R g (PP) (N)) 2 ⟩ 0 /σ 2 ≈ 209, 424, and 895 for the PPs of chains. Correspondingly, the mean-square end-toend distances are ⟨R e 2 (N)⟩ 0 /σ 2 = ⟨(R e (PP) (N)) 2 ⟩ 0 /σ 2 = 1329, 2575, and 5354, for N = 500, 1000, and 2000, respectively. Results of the two components of the rms radius of gyration of the OPs (PPs) in the directions parallel (||) and perpendicular (⊥) to the stretching, rescaled to 1/3 and 2/3 of the estimates of ⟨R g 2 ⟩ 0 1/2 (⟨(R g (PP) ) 2 ⟩ 0 1/2 ), respectively, right after deformation are shown in Figure 2. Under elongation, we see that the overall conformations of the OPs as well as of the PPs deform affinely. Namely, [⟨R g,|| 2 ⟩ λ /⟨R g,|| ) 2 ⟩ 0 ] 1/2 = λ −1/2 . The corresponding scaling laws also hold for ⟨R e,|| 2 ⟩ and ⟨R e,⊥ 2 ⟩, respectively (not shown). In order to estimate the length scale at which the deformation of chains becomes affine, we also include the results of the rescaled internal mean-square distances of the OPs (PPs) right after deformation, C⟨R 2 (n)⟩ λ /⟨R 2 (n)⟩ 0 (C⟨R ( P P ) (n) 2 ⟩ λ / ⟨R (PP) (n) 2 ⟩ 0 ), in Figure 2. Here n is the chemical distance between two bonds along the OP (PP) of the same chain, and C is the scale parameter. According to the scaling law for affine deformation, C = λ and C = 1/λ 2 for the two components in the directions perpendicular and parallel to the stretching, respectively. Our data demonstrate that n ≈ 8.2N e related to the chosen strain rate ε̇represents the characteristic length scale along the chains and above which both OPs and PPs of chains deform affinely. However, chains obviously cannot fully relax on short length scales since the connectivity and the constraints lead to deviations from this picture. This qualitatively compares to the expectation based on the strain rate chosen, which allows for local but not global relaxation during the deformation. For comparison, we also include data for the affinely (instantly) elongated original primitive path mesh of the unperturbed melt named the elongated PP mesh thereafter in Figure 2d. Our results of the mean-square internal distances show that PPs of chains in the PP mesh generated from the corresponding elongated OPs of chains in a melt through the PPA 4,5 display a deformation pattern very similar to the one of the OPs.
Conformational Change of Single Chains in Deformed Melts. In the nonlinear viscoelastic regime, the ) 2 (t)⟩/⟨(R g (PP) ) 2 ⟩ 0 ] 1/2 , plotted versus subsequent relaxation times, t/τ R,N , at λ ≈ 5.0 on a log−log scale. Data are for n c = 1000 primitive paths of chains of sizes N = 500, 1000, and 2000, as indicated. Data for the OPs of chains are also shown for comparison. In (b), four typical configuration snapshots of chain i = 300 of N = 2000 are included for better illustration, as indicated. At a chosen relaxation time t, the PP of the selected chain is shown by a curve in red, the OP is shown by beads in green, and the other 6 OPs of the same chain at t ± 1.32τ e , t ± 0.88τ e , and t ± 0.44τ e with τ e ≈ 2266τ ≈ 0.0002τ R,N=2000 are also shown by beads in gray. The minimum occur at t/τ R,N ≈ 0.07, 0.17, and 0.32 for the PPs of N = 500, 1000, and 2000, respectively, are pointed out by arrows. Theoretical predictions from the GLaMM model 44,45 for Z = 72, 36, and 18 are also shown by black curves from top to bottom in (a) and bottom to top in (b) for comparison.

Macromolecules
Article tube model predicts for the initial relaxation right after deformation an overdamped initial retraction process of the individual chains in both directions parallel and perpendicular to the stretching direction. 2,45 Such an immediate chain retraction mechanism has been observed for the time evolution of the rescaled two components of the radius of gyration for the OPs of chains 43,44 of sizes N = 1000 and 2000. One should also expect a similar behavior for the end-to-end distance of the OPs of chains, but with larger fluctuation due to the end effect. Therefore, we still only focus on the chain conformations described by the radius of gyration here. For the PPs the same holds as shown in Figure 3. Typical snapshots of a selected chain of size N = 2000 in the melt before and after deformation and after different relaxation times t/τ R,2000 = 0.004, 0.26, and 1.0 are also presented in Figure 3b for better illustration of the conformational variations. In both parallel and perpendicular directions to the stretching direction, ⟨R g,|| (PP) (t) 2 ⟩ and ⟨R g,⊥ (PP) (t) 2 ⟩ decrease initially with increasing time showing the signature of chain retraction. [⟨(R g,⊥ (PP) ) 2 ⟩ λ /⟨(R g,⊥ (PP) ) 2 ⟩ 0 ⟩] 1/2 reaches a deeper minimum compared to the OPs, while parallel to stretching, data for the PPs and OPs, respectively, coincide. The minima of (⟨R g,⊥ (PP) (t) 2 ⟩/⟨R g,⊥ (PP) (t) 2 ⟩ 0 ) 1/2 for PPs also occur slightly later than for OPs, 43 but the difference for the occurrence times becomes negligible within fluctuation as the chain size increases (at t/τ R,N ≈ 0.30(4) (OPs), 0.32(4) (PPs) for N = 2000, at t/τ R,N ≈ 0.09(5) (OPs), 0.17(4) (PPs) for N = 1000, and at no time (OPs), 0.07(4) (PPs) for N = 500). The minimum becomes deeper with increasing N. However, the duration of this global retraction process is still below the predicted longest time scale, 2,45 i.e., the Rouse time of the whole chains in unperturbed melts, τ R,N . This indicates a rather strong contribution from tension along the PPs. Setting the two parameters c ν = 0.1 and R s = 2.0 in the GLaMM model 45 (see the Supporting Information in ref 44), results of the time evolution of the radius gyration perpendicular and parallel to the stretching direction have been obtained by solving the constitutive equation iteratively in our previous work. 44 Obviously, our results for both the OPs and PPs qualitatively capture the signature of the initial chain retraction mechanism 2,45 after a large step elongation, while quantitatively the GLaMM model predicts a significantly stronger signature of retraction for the same chain size N, i.e., the same number of entanglements Z = N/N e . Furthermore, the rate of retraction  (t) 2 ⟩ 0 ) 1/2 indicates a turn to an intermediate plateau (which is more pronounced for N = 2000 than for N = 1000) showing a substantially delayed conformational relaxation well above and significantly earlier than the regime predicted by the GLaMM model. This relaxation retardation of deformed chains, not accounted for in current theoretical models, has been attributed to an inhomogeneous distribution of entanglement points along the PPs 43 and will be discussed later. A similar delay has also been observed in the context of rheological experiments of very long, highly entangled polymer chains by several other authors. 64−66 Transient Entanglement Effects. In the following we focus on understanding the relaxation of effective entanglement points along the PPs. Figure 4 shows the comparison between the curvature characterized by the local bond angle maxima θ j,j+5 ≥ 60°(significant kinks) and the tension force F j,|| (PP) = −▽U FENE (PP) ·ex in the direction parallel to the stretching direction along the PP of chain i = 2 of size N = 2000 (a typical case). Before the system is elongated, significant kinks are roughly equally distributed along the PP of the chain. The actual number of kinks is N kink (i = 2) = 99, while it varies between 40 and 100 for the whole system. The sign or magnitude switches of the tension force pattern and the kinks are closely correlated (see Figure 4a). To set the stage for comparison, we also take the original PP mesh of the unperturbed melt and deform this affinely up to λ ≈ 5. The distribution of kinks along the PP of chain i = 2 in this case remains very similar (see Figure 4b) to the original one; however, some kinks become sharper, some less sharp, and the number of kinks is slightly smaller, N kink (i=2) = 83. Differently from that, kinks along the corresponding PP of chain i = 2 in the elongated melt not only become sharper but also the number of kinks is reduced, N kink (i=2) = 68 (see Figure 4c). At the same time, the correlation between curvature and force pattern becomes even more pronounced. The linear correlation between the kinks of high curvature and sign switches of the tension force along the PP has been shown in the Supporting Information of ref 43. Our results indicate that the entanglement points, i.e., significant kinks, already in the very beginning do not follow the affine deformation if compared to the results for the elongated PP mesh (see Figure 5), while the average conformation of chains does, cf. Figure 2. This inhomogeneous distribution even becomes more pronounced upon relaxation of the deformed system up to 0.5τ R,N=2000 (N kink (i=2) = 52) or even up to 1.0τ R,N=2000 (N kink (i=2) = 47). Obviously, the length of the regions without force sign change is increasing instead of decreasing, and kinks become less sharp. At t/τ R,N = 2.0, the inhomogeneous pattern still persists although the number of kinks (N kink (i=2) =66) again increases. It should be kept in mind that these processes are subject to (up to chain end effects) the (approximate) conservation of topological constraints, but N kink can vary. This implies that entanglement points initially do not redistribute along the chain backbone as the chain tries to retract within a still globally affinely deformed tube. Qualitatively all chains display very similar patterns. Taking the clustering of entanglement points, i.e., the most confining topological constraints, one expects more conformational freedom between these regions. This is reminiscent of knotted polymers, where entropic forces tend to pull knots tight 67,68 (e.g., jamming knots). In this "jammed" regime, the knot's diffusivity decays exponentially above a critical tension force such that monomers can only move slowly due to high monomeric friction. 68 This is consistent with our finding that the state of the system is stabilized by the significant clustering of original kinks along the PP.
The theoretical predictions of the GLaMM model are based on the assumption of a constant Kuhn length of the PP. The authors assume that for deformation rates as in the case of the present study the melt structure on the scale of the tube diameter d T is only weakly perturbed from that of an equilibrium melt, d T ≈ d T (0) (d T (0) being the tube diameter in an equilibrium melt) since entanglements should be considered as mutual, delocalized topological interactions acting on a length scale of d T . To facilitate this either the number of entanglement points would have to increase while the Kuhn length of the PP is not changing or the other way around. Another possibility would be that a significant amount of topological constraints does not lead to kinks, e.g., like in a mesh where many chains are somewhat aligned. Nevertheless, this is not the case here. In Figure 5a, we see that differently from such a naive extension of the original tube model, the average number of kinks for elongated melts, ⟨N kink ⟩ λ , neither increases nor remains constant. This might be partially due to the reason that the cross section of the tube perpendicular to the tube axis is not always a perfect circle but can be rather elliptical under deformation (see Figure 2b). Comparing to the estimates of ⟨N kink ⟩ λ along all PPs in elongated PP meshes, we see that both sets of data are quantitatively the same within error bars under small perturbation while the deviation between these two sets of data becomes stronger with the

Macromolecules
Article increase of strain λ. Thus, topological constraints in polymer melts under large deformation do not follow affine deformation, but polymer melts under small deformation do (see Figure. 2d). The results of the affinely elongated PP mesh also indicate that a significant amount of topological constraints does not lead to kinks. Apparently, our results indicate a limitation of the GLaMM model, namely the assumption of the unperturbed Kuhn length and the homogeneous distribution of kinks breaks down. During subsequent relaxation, the characteristic behavior of ⟨N kink (t)⟩ shown in Figure 5b is very similar to the retraction of the overall chain (see Figure 3b). ⟨N kink (t)⟩ initially decreases, reaches a minimum at t/τ R,N ≈ 0.23(3), 0.40(3), and 0.57 (12) for N = 500, 1000, and 2000, respectively, and then slowly begins to increase to eventually approach ⟨N kink ⟩ 0 . Taking the data for N = 500 this final relaxation seems to occur on the time scale of the disentanglement time τ d ≅ τ R,N (N/N e ) x , x = 1 in the original reptation scheme 1,2 and x ≅ 1.4 experimentally. 69 PP strands between two neighboring entanglement points are almost straight lines (see Figure 1). Thus, a PP consists of straight segments of fluctuating length, and the average length of straight segments gives the length of an "effective entanglement length", which in the case of unperturbed , where the two ends of the chain are treated as kinks. Thus, results of ⟨l str ⟩ are expected to be proportional to 1/⟨N kink ⟩. The distributions of l str and l str,k , P(l str , t) and P(l str,k , t), respectively, right after deformation, and upon relaxation to t/τ R,N = 1.5 for N = 2000 are presented in Figure  6. For comparison, results for unperturbed polymer melts and for elongated polymer meshes are also included in Figure 6. There the inverse of the estimated slope in Figure 6b roughly corresponds to ⟨l str,k ⟩ ≈ 22 which is close to N e,PPA (0) = 28. We see that the profiles of P(l str , t) are quite well described by shifted Gaussian distributions. The peak first shifts to significantly larger values of l str , while only after t/τ R,N > 0.5 the peak of P(l str , t) slowly begins to move toward smaller values of l str , still being far away from the unperturbed case, even for t ≈ 1.5τ R,N . Furthermore, the width of the distribution (⟨l str 2 ⟩ − ⟨l str ⟩ 2 ) seems to become even wider with time up to t ≈ 0.5τ R,N . Within the chains the distribution of the lengths of straight PP segments P(l str,k , t), Figure 6b, shows long tail distributions of l str,k . Right after deformation, the tail of P(l str,k , t) becomes significantly broader than for an unperturbed polymer melt. At subsequent relaxation times, instead of moving back to the distribution for unperturbed polymer melt, we see that the tail initially becomes even significantly longer up to t ≈ 0.5τ R,N , before it very slowly turns toward the distribution for the unperturbed melt as shown in Figure 4, which it eventually has to reach. The apparent slope for l str,k > 120 would indicate an intermediate effective N e of about 100 based on a small fraction of the total straight PP segments, which should not be confused with the entanglement length based on the theoretical considerations for the PPA, 4,5 discussed later. For the elongated PP mesh at λ ≈ 5.0, we see only a slight shift to larger values of l str for P(l str , t) compared to that of the unperturbed polymer melt since ⟨N kink ⟩ is only about 10% smaller than ⟨N kink ⟩ 0 right after deformation (Figure 5a). However, such effects are too weak to analyze them quantitatively, as they are within the error bars of the profiles of P(l str,k , t).
Our results clearly show that the number of kinks along all PPs in polymer melts under large deformation does not follow affine deformation, and the kinks along the individual PPs of chains distribute unequally right after deformation. The delayed relaxation observed in the profiles of ⟨N kink (t)⟩, P(l str , t), and P(l str,k , t) indicate that topological constraints play an essential role in the relaxation retardation of deformed chains.
Entanglement Point Distribution in Space. So far we have been focusing on conformations of individual chains and their respective primitive paths. We now turn to the distribution of topological constraints, as represented by entanglement points, in space. Figure 7 shows snapshots of such distributions as a function of time. Quantitatively, the distribution can be described by the corresponding density of entanglement points projected onto the x−y plane, ρ kink (x, y). ρ kink (x, y) is estimated by simply counting the number of entanglement points located at (x ,y) and normalized such that . For this the grid spacing is set to 2.0σ. For a large disorder system, structural inhomogeneities presented in such a projection would be smeared out. Thus, one can view Figure 7 as a rephrase native slice of a large system. Before and right after deformation, the entanglement points distribute homogeneously while the kinks become sharper. For the unperturbed melt (λ = 1.0), N kink = 88 128 for 60°< θ j,j+5 < 90°and N kink = 7617 for θ j,j+5 > 90°while N kink = 7920 for 60°< θ j,j+5 < 90°and N kink = 74 344 for θ j,j+5 > 90°for the polymer melt right after deformation (λ ≈ 5.0). Upon relaxation, we observe a clustering of entanglement points, and the clustering pattern becomes more distinct and seems to reach a maximum of around t = 0.5τ R,N . Even up to t = 2.0τ R,N , these clusters persist although the kinks along deformed chains become less sharp (N kink = 46 917 for 60°< θ j,j+5 < 90°and N kink = 8888 for θ j,j+5 > 90°). Since the larger clusters of entanglement points are rather fuzzy and since there are still many isolated entanglement points, it is difficult to identify a characteristic length scale which leads to the instability in the homogeneous distribution of entanglement points. Based on the original tube concept one would expect that the onset of the separation into these jammed areas should be of the order of at most a few tube diameters (times the elongation amplitude) since this is the only length originating from the topological constraints. The one-dimensional scattering function S (kink) (q || ) of entanglement points in the melt gives an impression of the correlation, where q || is the component of the vector q in the direction parallel to the stretching. The discretization in q-space along the stretching direction is given by q || = 2πm x /L x with m x = 0, 1, 2, ..., where L x = L 0 and L x ≈ 5.0L 0 before and after deformation, respectively. Results of S (kink) (q || ) plotted versus q || are shown in Figure 8. Here we only focus on the regime q || ≤ 1.0σ −1 . At first glance, the distributions of entanglement points seem to be more structured during initial relaxation after polymer melts are deformed compared to the structure for the unperturbed polymer melt. This is consistent with the observed clustering of entanglement points shown in Figure 7. If we only focus on clusters with high density of kinks, the characteristic distance d cluster can be roughly determined by d cluster = 2π/q || max , i.e., at q || = q || max , S (kink) (q || ) reaches a maximum along the stretching direction. For t = 0.25τ R,N , we see one broad maximum around q || max ≈ 0.08σ −1 which gives d cluster ≈ 78σ ∼ 16d T (0) , consistent with the above argument. At subsequent relaxation times, this broad maximum slowly moves to a larger value of q and becomes even broader. This clustering structure remains even up to t = 2.0τ R,N , and d cluster ≈ 40σ ∼ 8d T (0) (q || max ≈ 0.16σ −1 ). Thus, the clustering of the entanglement points is related to the relaxation retardation of the deformed polymer melt.
Mobility of Monomers in Nonequilibrium States. Polymer chain dynamics is usually characterized by different regimes of the mean-square displacement (MSD) of monomers. Reptation theory 1,2 predicts the crossover between these regimes at the characteristic time τ 0 for local fluctuations, the entanglement time τ e ∝ N e 2 , the Rouse time τ R,N ∝ N 2 , and the disentanglement time τ d,N ∝ N 3.4 (when contour length fluctuations, constraint release, and correlation hole effects are taken into account 37,38,70,71 ). The crossover scaling predictions of the MSD of monomers, g 1 (t), of monomer with respect to the corresponding center of mass, g 2 (t), and of the center of mass g 3 (t) have also been verified by our large time scales MD simulations for highly entangled and fully equilibrated polymer melts, 20 i.e., the unperturbed systems studied here.
In Figure 9, we compare the motions of inner (N/2 + 1) monomers (i.e., eliminating the strong fluctuations caused by chain ends 11,14,20 ) parallel and perpendicular to the stretching direction, g 1,|| (t) and g 1,⊥ (t), respectively, and of the center of mass, g 3 (t), in the polymer melt of size N = 2000 after deformation to that in the fully equilibrated polymer melt. Note that we here do not average over starting times as we begin to measure g 1 (t) and g 3 (t) right after deformation. We see that along the stretching direction, g 1,|| (t) for the deformed polymer melt initially follow the same scaling behavior up to the original tube diameter d T (0) as that for the unperturbed polymer melts (g 1,|| (t) ∼ t 1 for t < τ 0 , and g 1,|| (t) ∼ t 1/2 for τ 0 < t < τ e = τ 0 (N PPA (0) ) 2 ). For t > τ e , monomers move faster, g 1,|| (t) ∼ t 0.69 , until reaching t ≈ 0.30τ R,2000 , i.e., the duration of the chain retraction process. Within this time frame, the change of g 1,|| (t), 3 is of the same order as the change of ⟨R g,|| 2 (t)⟩, while the center of mass displays unperturbed diffusion g 3 (t) ∼ t up to 0.30τ R,2000 , but moves a little more slowly around t = τ e . For t > 0.3τ R,2000 , monomers' motion is slowed down again due to entanglement effects and the corresponding relaxation retardation, similar as in a tubelike regime created by the surrounding chains, g 1,|| (t) ∼ t 1/4 . However, the effect of entanglements varies in the process of equilibration. A strong relaxation retardation is also observed for g 3 (t) that g 3 (t) ∼ t 1/4 . Perpendicular to the stretching, the situation is somewhat different. Apparently, monomers move more slowly that a gradual deviation from g 1,⊥ (t) ∼ t 1/2 to g 1,⊥ (t) ∼ t 1/4 develops as t approaches τ e . After deformation, chains are somewhat aligned along the stretching direction such that monomers presumably have less freedom to move in the crowded space along the perpendicular direction.
To test our finding that the deformed chain conformations are stabilized state by the clustering/jamming of entanglement points, we estimate the MSD of groups of 20 inner monomers in the clustering (jammed) region and 20 monomers in the free region for all n c = 1000 chains. Here the clustering and free regions are identified according to the curvature the PPs of chains at t = τ R,N . For example, in the case of i = 2 (see Figure  4e), monomers j = 960 to j = 979 in the clustering region and monomers j = 1480 to j = 1499 in the free regime are considered. This is compared to MSD of inner (N/2 + 1) monomer during the relaxation process. Results of g 1 (t) plotted as a function of t/τ R,2000 are shown in Figure 9d. We see that monomers in the constrained regime move much slower compared to those in the free regime. Indeed, the MSD of inner monomers is dominated by the motion of inner monomers in the constrained regime. In Figure 9d, we have also included the MSD of 20 monomers in the free regime near one end of all chains. Monomers apparently move much faster Figure 10. Average bond length ⟨b PP ⟩ λ (a) and orientational order parameter ⟨Q PP (ϕ)⟩ λ (c), plotted versus λ for elongated polymer melts. Average bond length ⟨b PP (t)⟩ (b) and orientational order parameter ⟨Q PP (ϕ,t)⟩ (d), plotted versus the subsequent rescaled relaxation time t/τ R,N after elongation. In (a), (b), data are rescaled to ⟨b PP ⟩ 0 ≈ 0.31σ for the PPs of unperturbed chains in melts. In (a), (c) data for elongated PP mesh, and the theoretical predictions for isotropic PP, eq 6 and eq 8, respectively, are also shown for comparison.

Macromolecules
Article as we have expected due to the weak topological constraints at the end and the retraction process.
Intrinsic Properties of Primitive Paths and the Tube Picture. Based on the tube picture, 1,2,9 one would expect that entanglement effects appear at t ≈ τ e and monomers are restricted to move along the contour of an imaginary tube of diameter d T created by surrounding chains. The primitive paths represent the backbone of the tube that can be constructed following this picture. In this subsection we investigate the time-dependent intrinsic properties of PPs created by PPA 4,5 for polymer melts under strain and link them to the tube concept.
The where r i,j (PP) is the coordinate of the jth monomer of the PP of chain i. If we assume the distribution of the original bonds to be isotropic, taking the integral of the deformed bonds over the surface of a unit sphere, normalized by 4π, the average bond length as a function of the strain λ is given by where E is a diagonal deformation tensor with elements {λ, λ −1/2 , λ −1/2 }, u is a unit vector in the spherical coordinate, u = (cosϕ, sinϕcos θ, sinϕsin θ), ϕ is the polar angle between u and the x-axis, and θ is the azimuthal angle. One gets ⟨b PP ⟩ λ ≈ 2.56⟨b PP ⟩ 0 at λ = 5.0. The deviations between different approximations discussed in our previous work 43 are a consequence of the fact that there is a lower cutoff, given by N e . In general one should expect that ⟨b PP ⟩ λ /⟨b PP ⟩ 0 for the affinely elongated PP mesh follow eq 6, which indeed is seen in Figure 10a. At the same time the deviation between the elongated original PP mesh and the PP meshes of deformed polymer melts becomes more pronounced with increasing strain λ. At λ = 5.0, ⟨b PP ⟩ λ ≈ 1.84⟨b PP ⟩ 0 for the deformed polymer melt, which is about 28% below that of the deformed original PP mesh. At the subsequent relaxation after stretching, we see that in Figure 10b, ⟨b PP (t)⟩ reaches ⟨b PP ⟩ 0 for N = 500 and N = 1000 for t > τ R,N while for N = 2000, ⟨b PP (t)⟩ seems to settle at a slightly larger value than ⟨b PP ⟩ 0 .
Choosing the x-axis (along the stretching direction) as a reference, the orientational order parameter of bond vectors along PPs is defined by Q PP (ϕ) = (3cos 2 ϕ − 1)/2 where ϕ is the angle between the bond vector b PP and the x-axis. Assuming the distribution of the original bonds to be isotropic, ⟨Q PP (ϕ)⟩ λ for the PP meshes of polymer melts under elongation is given by Under elongation Figure 10c shows that the orientation distribution of bonds changes from isotropic, ⟨Q PP (ϕ)⟩ λ = 0.0 for λ = 1.0, to anisotropic, ⟨Q PP (ϕ)⟩ λ ≈ 0.8 for λ = 5.0, where most bond vectors align along the stretching direction. The subsequent relaxation of bond orientation (Figure 10d) is apparently significantly delayed compared to that of the bond length. Furthermore, the relaxation rate decreases as the chain size N increases, Despite quantitative discrepancy for N = 2000 both ⟨b PP (t)⟩ and ⟨Q PP (ϕ,t)⟩ show relaxation retardation starting at t ≈ τ R,N .
Finally, we turn to the entanglement length N e,PPA as estimated from the PPA for unperturbed and strongly deformed polymer melts, cf. Figures 11−13. For strongly deformed polymer melts, the estimate of N e,PPA becomes less obvious since the basic original concept assumes a Gaussian chain conformation. The Gaussian chain assumption still holds for each component individually; however, it is different in different directions since the contours of PPs globally deform affinely. Keeping this in mind, we nevertheless apply the standard formula ⟨R e 2 ⟩ = ⟨(R e (PP) ) 2 ⟩ = L PP l K (PP) = (N − 1) N e,PPA ⟨b PP ⟩ 2 , l K = N e,PPA ⟨b PP ⟩ being the Kuhn length of the PPs of chains. 4,5 For the overall affine deformation of the chains along the three orthogonal directions (see Figure 2), one would expect ⟨R e 2 ⟩ λ = λ 2 ⟨R e,x 2 ⟩ 0 + ⟨R e,y 2 +R e,z 2 ⟩ 0 /λ = (λ 2 + 2/ λ)⟨R e 2 ⟩ 0 /3. Following eq 6, we obtain for the elongated PP mesh follows affine deformation up to λ ≈ 5.0 while for the PPs of the deformed polymer melts it does not. This is consistent with the estimates of "effective entanglement length" ⟨l str ⟩ λ ∝ 1/ ⟨N kink ⟩ λ (Figure 5a). At subsequent relaxation time t ≈ 0.5τ R,N , N e,PPA (t) reaches a maximum for N = 500 and 1000 while for N = 2000, it seems to reach a plateau value without any further decay within the time window studied here, and the plateau value is very close to the effective entanglement length (100) extracted from the long tail probability distribution at t = τ R,N . This is another direct piece of evidence for a significantly delayed relaxation.
Along the original path the monomers can fluctuate in space, confined by fuzzy tube boundaries. If the average monomer positions and thus the tube are deformed affinely, one would expect that the average number of monomers located inside the tube of tube diameter , ⟨N tube ⟩, is approximately kept as a constant. The GLaMM model even assumes that d T = d T (0) remains unchanged, which suggests that even more monomers would be located inside the tube since the contour length of the tube becomes larger. From our results shown in previous sections, we find that this cannot be expected here. Still, following the assumption of the GLaMM model and fixing the tube diameter to that of t h e f u l l y e q u i l i b r a t e d m e l t , n a m e l y t o 2 0 for each chain i and ⟨N tube ⟩ upon deformation and subsequent relaxation. The way of counting N tube (i) for each chain i is as follows: To determine whether monomer k is located within its actual reptation tube, we assume the tube to be constructed of piecewise cylinders with diameter d T (0) and length d T (0) . Monomer k belongs to the tube if

Macromolecules
Article the shortest distance to any PP bond j with k − δ ≤ j ≤ k + δ is less or equal to d T (0) /2 with δ = 8, being determined by d T Note that the precise number of N tube depends on details of the algorithm and the chosen value for δ, but results of ⟨N tube ⟩ should be qualitatively the same. To provide some insight into the relation between the location of kinks and the fluctuation of monomers of OP along the PP for the same chain, we choose the same chain i = 2 of size N = 2000 as it is selected in Figure 4 for comparison. Indeed, the monomer k of the OP belonging to the bond j along the PP of chain i = 2 in a melt presented in Figure 12 shows approximate agreement with the data before deformation and at several selected subsequent relaxation times, but not right after deformation. This makes sense since monomers along the OP have less freedom moving away from the tubelike regime along its corresponding PP as the chain is strongly stretched, and thus N tube (i) increases. Once the deformed chain starts to relax, N tube (i) immediately decreases due to the dramatic changes of the entangled surrounding, and then monomers of OP between two neighboring kinks can start to fluctuate.
⟨N tube ⟩, cf. Figure 13a, b, fits to the overall scheme observed so far. In the linear regime (λ < 1.5), ⟨N tube ⟩ λ ≈ ⟨N tube ⟩ 0 . In the nonlinear regime (λ > 1.5), upon deformation ⟨N tube ⟩ λ increases with a rate larger with larger N. Right after deformation, at the subsequent initial relaxation time, ⟨N tube (t)⟩ decreases, if normalized by ⟨N tube ⟩ 0 . Eventually, it even drops below the equilibrium value, slows down for t/τ R,N > 0.2, and reaches a minimum around the Rouse time. However, the minimum is quite shallow and shows the delayed relaxation around the Rouse time. The relaxation back to the unperturbed tube occupancy at least takes several Rouse times of the chains, just as for other tube related quantities studied here. Choosing d T = d T (0) (N e,PPA (λ≈5.0) /N e,PPA (0) ) 1/2 ≈ 7.61σ ≈ 1.5d T (0) , the profiles of ⟨N tube (t)⟩/⟨N tube ⟩ 0 having the minimum around t ≈ τ R,N within fluctuation are similar to the results shown in Figure 13b, while quantitatively, the estimates of ⟨N tube (t)⟩ are somewhat larger. Again this agrees with the formation and growth of topologically highly congested areas along the chains. Regions of the low density of entanglement points, where configurations of monomers along the OPs can fluctuate significantly in space, seem to stabilize regions with the high density of entanglement points. The probability distributions P(N tube /N) as a function of (N tube /N) for deformed polymer chains of size N = 2000 in a melt at several selected strain values of λ and at several selected rescaled relaxation times t/ τ R,N after deformation are shown in Figure 13c, d. We see that the P(N tube ) is simply a shifted Gaussian distribution in terms of the mean value ⟨N tube ⟩ and the standard deviation 2 . The distribution P(N tube /N) remains the same for λ = 1.0 and λ ≈ 1.5, and then the profile of P(N tube /N) shifts to a larger value of N tube /N as λ increases. After deformation, the profile shifts to a smaller value of N tube / N with increasing t, and even moves to the left-hand side of the profile for unperturbed chains similar as the behavior of ⟨N tube (t)⟩ shown in Figure 13b. We see that the intrinsic properties of PPs, analyzed according to the tube concept, provide profound insights into the relaxation paths of entangled chains in deformed polymer melts although the "effective entanglement length" for deformed chains in a melt can no longer be extracted from the theoretical considerations for the PPA. 4,5 Stress Relaxation. After having discussed details of individual and collective conformational relaxation, we turn to the related stresses in the systems, a quantity that would be experimentally accessible more easily. Since the entanglement structure of melts is closely connected to their viscous and elastic properties, we expect that the previously described primitive path relaxation also leads to characteristic signals in the viscoelasticity of the strained melts. The viscoelasticity of polymer melts is normally characterized by the time-dependent stress relaxation modulus G(t). For fully equilibrated, entangled polymer melts, the short-time dynamics of chains is described by the Rouse model. This leads to G(t) ∼ t −1/2 for t < τ e , while for τ e < t ≪ τ d,N (τ d,N = τ R,N (N/N e ) 1.4 being the disentanglement time of chains of size N) where chains are assumed to move in a tubelike regime due to entanglements, G(t) reaches a plateau value G N 0 = (4/5)(ρk B T/N e ), which depends on the entanglement length or the molecular weight between entanglements as predicted by the Doi−Edwards tube model. 2,9 In the linear viscoelastic regime, G(t) and its approach to G N 0 for intermediate time scales are well understood. In contrast, our understanding of the stress relaxation scenario for strongly deformed polymer melts in the nonlinear viscoelastic regime and the relationship between the strain rate and stress relaxation is significantly less developed. 72 Results of the stress relaxation modulus G(λ, t) that characterize the viscoelasticity of entangled polymer melts in both linear (unperturbed and λ ≈ 1.2) and nonlinear regimes (λ ≈ 5.0) are shown in Figure 14 (parts of data have been shown in refs 20 and 43). For clarity only the case of N = 2000 (≈72N e ) is discussed here. For equilibrium polymer melts G(λ = 1.0, t) is calculated from the stress autocorrelation function of off-diagonal elements of the stress tensor, 20,73,74 using the Green−Kubo relationship. G(λ = 1.0,t) reaches a plateau value G N 0 with N e ≈ 28 for τ d,N ≫ t > τ e ≈ 2 × 10 −4 τ R,N . Alternatively G(λ, t) can also be given by the normal stress difference after the simulation box is deformed, i.e., stretched along the xdirection, nomalized by the damping function h(λ) with the neo-Hookean prediction and h(λ)=λ 2 −1/λ. For small and fast deformation (λ ≈ 1.2 and ετ R,N = 32 000), G(λ ≈ 1.2,t) follows the data 20 calculated from the Green−Kubo relation for t > τ e and is in good agreement with experiments 75,76 and the tube model. 1,2,9 For the large but relatively slow deformation (λ ≈ 5.0 and ετ R,N = 77), G(λ ≈ 5.0, t) already from the very beginning displays a plateau below the plateau modulus G N 0 of the linear regime. This is consistent with the experimental finding that one observes elongational thinning in melt under such conditions. 51,77 Only after a time corresponding to the inverse strain rate, i.e., t = ε̇− 1 = τ R,N /77 ≈ 10 −2 τ R,N a further decay of G(t) is observed. Keeping λ ≈ 5.0 but choosing the same very fast strain rate ετ R,N = 32 000 as for λ ≈ 1.2, we see that G(λ ≈ 5.0, t) initially follows the data of G(λ ≈ 1.2, t) up to ε̇− 1 = τ R,N /32 000 ≈ 10 −5 τ R,N . Then it deviates and after the relaxation time reaches t = ε̇− 1 =τ R,N /77, both curves for λ ≈ 5.0 follow the same softening pattern. This indicates that relevant chain chain interpenetration did not change significantly during the slow stretching process (ε̇= 77/τ R,N which is still quite fast compared to τ d,N −1 ), that is, topological constraints are not released up to some chain end effects. Furthermore, this supports the concept that the longtime behavior of the stress relaxation also in the nonlinear regime seems to be just a function of the ultimate deformation and is independent of the original strain rate τ R,N −1 < ε̇< τ e −1 . From the above we expect that estimates of ⟨N kink (t)⟩ for deformed polymer melts with two different strain rates at λ ≈ 5.0 in Figure 14b should eventually coincide. Indeed, as expected for t/τ R,N > 1/77 ≈ 0.013, the two sets of data perfectly coincide with each other, i.e., entanglement effects in both deformed polymer melts remain the same.

■ CONCLUSIONS
In this paper, we have employed extensive molecular dynamics simulations to study highly entangled polymer melts subject to an isochoric elongation in the nonlinear viscoelastic regime. Focusing on the analysis of topological constraints, as identified through a primitive path analysis, 4,5,43 our simulation results can serve as benchmarks and starting points for further experimental and theoretical studies of monodisperse polymer melts under a large step elongation. An affine deformation of the overall conformations of polymer chains in a melt is to be expected by setting the strain rate τ R,N −1 < ε̇< τ e −1 (Figure 2). The chain retraction mechanism predicted by the Doi− Edward tube model 2 and its refined GLaMM tube model 45 after large step deformation are demonstrated to hold qualitatively ( Figure 3). Perpendicular to stretching, the signature of predicted chain retraction enhances with the increase of molecular weight, i.e., number of entanglements. However, the onset of delayed relaxation in the conformation of chains occurs earlier, and the duration of chain retraction process is shorter than the predicted Rouse time of chains in unperturbed melts. 2,45 Since the central assumption of an isotropic Gaussian chain conformation is not valid anymore here, entanglement points are identified as significant kinks along the PP by comparing the curvature and the tension force pattern along the PP (Figure 4). The resulting average number of entanglement points (kinks) ( Figure 5) shows that the distribution of the entanglement points does not deform

Macromolecules
Article affinely upon elongation in the nonlinear viscoelastic regime (λ > 1.5). At subsequent relaxation times, the average number of kinks first decreases and reaches its minimum value before it moves toward the value for unperturbed polymer melts again implying delayed relaxation. The distributions of "effective entanglement length" between two neighboring kinks along the same primitive path for the whole deformed polymer melt also show the similar behavior ( Figure 6). Furthermore, we have observed a clustering and inhomogeneous distribution of entanglement points not only along the individual PPs of chains ( Figure 4) but also for the whole melts (Figure 7). The distance between these jammed areas is found to be of the order of a few tube diameters of unperturbed melts. Tracking the mean-square displacement of inner monomers and of the center of mass right after polymer melts are deformed, significant deviations from that for unperturbed polymer melts are discovered for t > τ e (Figure 9). The relaxation retardation picture has also been confirmed by investigating the time-dependent intrinsic properties of PPs of stretched polymer melts characterized by the bond length, the orientational order parameter, and the entanglement length estimated using the original recipe of the PPA 4,5 ( Figures  10−11). However, all these quantities do not follow affine deformation under strain. Our results also show that the average internal conformations of chains are straightened for n > N e (Figures 2c, d), and most bonds are weakly aligned along the stretching direction (+x-axis) (Figures 10c, d) with the strain rate ετ R,N = 77, however, not enough for nematic order. The similar behavior has been previously observed in the study of nonlinear extensional flows for Rouse−Weissenberg number 78 Wi R > 1, there the chains are much shorter (z < 18). Following the tube picture assumption, we also count the number of monomers confined in a tubelike regime of fixed tube diameter d T = d T (0) for the unperturbed polymer melt to measure the transverse fluctuation of monomers during the elongation and relaxation processes (Figures 12, 13). Our results again show a significantly delayed relaxation compared to the current theoretical considerations. 45 Finally, we see that the time evolutions of the stress relaxation modulus for deformed polymer melt only initially depend on the strain rates while they follow the same softening patterns ultimately in the nonlinear viscoelastic regime (Figure 14a). The same scenarios are also seen from the time evolutions of the number of effective entanglement points (Figure 14b). All our results indicate that in the nonlinear viscoelastic regime, the topological constraints in highly strained large polymer melts are better described by the PP meshes than the corresponding elongated OPs of chains. The resulting entanglement effects play an important role in the relaxation retardation. However, it is still unclear how long the relaxation retardation would last and what could be the time for equilibrating the deformed polymer melts if it would happen. For all our data, where we observe a characteristic change from initial to later time relaxation, we observe deviations from Rouse scaling. While the equilibration time of deformed chains is expected to follow the same power law as that of τ d,N , the positions of the minima in, e.g., Figure 5 suggest a different time scale τ new,N = Aτ R,N N x in between the Rouse and the disentanglement time. Assuming a constant prefactor A the minima are located at a time proportional to N 2+x with an exponent x somewhere in between 0.5 and 0.8. Any definite statement here would, however, require more systematic studies with different chain lengths and elongation schemes. These findings so far have not been considered in the current theoretical framework, and similarities to knotted polymers are worth further considerations. The clustering and inhomogeneous distribution of entanglement points offer interesting options for both experimental and simulation investigations, especially in the vicinity of the glass transition point. 79,80 Especially, it is important to understand whether the complex topological constraints in deformed polymer melts are the crucial polymer characteristic for understanding the dynamical and thermodynamic properties of polymer chains in a deformed melt.