Clustering of entanglement points in highly strained polymer melts

TL;DR

This study uses molecular dynamics simulations to analyze how entanglement points in highly strained polymer melts cluster and relax, revealing persistent inhomogeneities and the dynamics of topological constraints during deformation.

Contribution

It introduces a detailed analysis of entanglement point clustering and relaxation in highly deformed polymer melts using primitive path analysis.

Findings

Entanglement points cluster persistently after deformation.

Primitive path conformations follow affine deformation.

Relaxation involves long-lived inhomogeneous distribution of constraints.

Abstract

Polymer melts undergoing large deformation by uniaxial elongation are studied by molecular dynamics simulations of bead-spring chains in melts. Applying a primitive path analysis to strongly deformed polymer melts, the role of topological constrains in highly entangled polymer melts is investigated and quantified. We show that the over-all, 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.

Equations

$$U_{\mathrm{LJ}}(r)$$

$$U_{\mathrm{FENE}}(r)$$

$$U_{\rm BEND}(\theta)=k_{\theta}(1-\cos\theta)\$$

$$\dot{\epsilon }=\frac{d\epsilon }{dt}=\frac{1}{L}\frac{dL}{dt}=\frac{d\ln (L/L_0)}{dt},$$

$$⟨b_{PP}⟩=\frac{1}{n_c(N-1)}\sum _{i=1}^{n_c}\sum _{j=1}^{N-1}\mid {\mathbf{r}}_{i,j+1}^{(PP)}-{\mathbf{r}}_{i,j}^{(PP)}\mid$$

$$⟨b_{PP}{⟩}_{\lambda }$$

$$h(\lambda )={\lambda }^2-1/\lambda$$

$$⟨Q_{PP}(\varphi ){⟩}_{\lambda }$$

$$\frac{N_{e,PPA}^{(\lambda )}}{N_{e,PPA}^{(0)}}=\frac{({\lambda }^2+2/\lambda )}{3}{\left\{\frac{\lambda }{2}+\frac{1}{4\lambda h(\lambda {)}^{1/2}}\ln \left|\frac{h(\lambda )+\lambda \sqrt{h(\lambda )}}{-h(\lambda )+\lambda \sqrt{h(\lambda )}}\right|\right\}}^{-2}.$$

$$G(\lambda ,t)={\sigma }_{norm}(\lambda ,t)/({\lambda }^2-1/\lambda )$$

Full text

1 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, 2, 3, 4, 5, 6, 7, 8. In the linear viscoelastic regime these are well described by reptation theory (the tube model) 1, 9, 2, and have been confirmed by many simulations 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20 and experiments 21, 22, 23, 24, 25. However, despite this remarkable achievement, a precise definition of an entanglement within the reptation concept still is lacking and attempts to include multi-chain effects analytically were of very limited success only 26. Everaers et al 4 introduced the primitive path analysis (PPA) based on the concept of Edwards’ tube model 27 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, 28. Kröger et al. 29, 30 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 multi-bead 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. 31 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. 32, 33, 34. Thus, there is 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, 35, 19.

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) 36, 37, 38 of the PPs, and constraint release (CR) 37, 39, 40, 41 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 above mentioned concepts can be employed for describing the behavior of polymer melts in the non-linear viscoelastic regime. To shed light on this problem, we recently have started to approach this problem from the computational side 42, 43. 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 44 to extensive simulations of polymer melts in the non-linear viscoelastic regime. For this we decided to concentrate on 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 45, 46, 47, 48, 49, 43 while our recent results for longer chains support this general scheme 43. Although primitive path network models can account for viscosity changes in entangled polymers upon elongational and shear flow 50, 51, 52, 53 in the linear viscoelastic regime, they seem to fail for strain rates $\dot{\varepsilon}$ faster than the inverse Rouse time $\tau_{R,N}$ of the whole chains. For example, the monotonic shear thinning behavior of entangled polymer melts in elongational flow observed from the experiment 54 at $\dot{\varepsilon}>1/\tau_{R,N}$ cannot be described by simulations using the primitive path network models 52. 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 $\rho=0.85\sigma^{-3}$, prepared by a new, efficient hierarchical methodology 55, 28, 19. These melts are subject to strong deformation by isochoric elongation in the non-linear 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 non-linear viscoelastic regime, we recently have shown that the force pattern along the primitive paths (PPs) qualitatively match that of the corresponding original paths (OPs) 42. This indicates that the conformations of OPs of chains within fuzzy tube-like 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 about 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.

2 Model and simulation methods

2.1 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 $\rho=0.85\sigma^{-3}=n_{c}N/V$ where $\sigma=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 non-bonded monomers located at a distance $r$ apart interact via a shifted, purely repulsive Lennard-Jones (LJ) potential 56, 57, 58, 59 $U_{\rm LJ}(r)$,

[TABLE]

where $\epsilon$ is the energy unit of the pairwise interaction, and $r_{\rm cut}=2^{1/6}\sigma$ is the cutoff in the minimum of the potential such that force and potential are zero at $r_{\rm cut}$. Any pair of bonded monomers interacts via the finitely extensible nonlinear elastic (FENE) binding potential 56 $U_{\rm FENE}(r)$,

[TABLE]

where $k=30\epsilon/\sigma^{2}$ is the force constant and $R_{0}=1.5\sigma$ is the maximum value of bond length. These Lennard-Jones units $\sigma$ and $\epsilon$ also provide a natural time definition via $\tau=\sigma\sqrt{m/\epsilon}$ where $m=1$ is the mass of the particles. In addition, a weak bond-bending potential 4 $U_{\rm BEND}(\theta)$ with a chain stiffness parameter $k_{\theta}$ is introduced,

[TABLE]

where $\theta$ is the angle between two subsequent bonds, i.e., $\theta=\cos^{-1}\left(\frac{{\bf b}_{j}\cdot{\bf b}_{j+1}}{\mid{\bf b}_{j}\mid\mid{\bf b}_{j+1}\mid}\right)$ and ${\bf b}_{j}={\bf r}_{j}-{\bf r}_{j-1}$ being the bond vector between monomers $j$ and $j-1$ along the chain. Choosing $k_{\theta}=1.5\epsilon$, where chains become weekly semiflexible, the mean square radius of gyration for unperturbed (i.e. fully equilibrated) chains in a melt is $\langle R_{g}^{2}\rangle_{0}\approx\langle R_{e}^{2}\rangle_{0}/6\approx 0.484N\ell_{b}^{2}$ where $\langle R_{e}^{2}\rangle_{0}$ is the mean square end-to-end distance and, $\ell_{b}=\langle{\bf b}^{2}\rangle_{0}^{1/2}\approx 0.964\sigma$ is the root-mean-square (rms) bond length. This gives $\langle R_{e}^{2}\rangle_{0}$, $\langle R_{g}^{2}\rangle_{0}\propto N^{2\nu}$ with the Flory exponent $\nu=1/2$. Here the average $\langle\ldots\rangle_{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 19. For the above parameters the corresponding entanglement length $N_{e}=N_{e,PPA}^{(0)}\approx 28$, estimated both from the plateau modulus $G_{N}^{0}=(4/5)\rho k_{B}T/N_{e}$ as well as from the primitive path analysis (PPA) 4, 5, 28, 19. The Rouse relaxation time of a subchain of entanglement length $N_{e}$, and of the overall chain of size $N$ is $\tau_{e}=\tau_{0}N_{e}^{2}$ and $\tau_{R,N}=\tau_{0}N^{2}$, respectively. Here the characteristic time prefactor $\tau_{0}\approx 2.89\tau$ is determined from the estimate of the mean square displacement of inner monomers 19 for $n_{c}=1000$ chains containing $18\leq Z\equiv N/N_{e}\leq 72$ entanglements, $g_{1}(t)$. We here focus on the above mentioned cases and use the molecular dynamics simulations package 60 ESPResSo++ for all runs with simulation time step set to $\Delta t=0.01\tau$. The temperature is kept constant ($T=1\epsilon/k_{B}$, $k_{B}$ being the Boltzmann factor) through a Langevin thermostat with a weak friction constant $\Gamma=0.5\tau^{-1}$.

2.2 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/\sqrt{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 $\dot{\varepsilon}$ in the range $\tau_{R,N}^{-1}<\dot{\varepsilon}<\tau_{e}^{-1}$ is introduced, namely

[TABLE]

i.e., $L=L_{0}\exp(\dot{\varepsilon}t)$. We set the strain rate $\dot{\varepsilon}$ to $\dot{\varepsilon}\tau_{R,N}=77$ that is $\dot{\varepsilon}\tau_{e}=77(N_{e}/N)^{2}$ (e.g. $\dot{\varepsilon}\tau_{e}\approx 0.015$ for $N=2000$) at each elongation step. To obtain this deformation rate the system can relax for $(0.02\tau_{R,N}/77)\tau$ between two such steps. Thus one could expect that locally, i.e. below a few $N_{e}$ $(\approx\sqrt{1/0.015}N_{e}\approx 8.2N_{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 42, 43 $\lambda\approx 5.0$. As a result the deformed polymer melts are deep in the non-linear viscoelastic regime.

2.3 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 bond bending interaction are switched off, i.e., $U_{\rm LJ}^{\rm(intra)}(r)=U_{\rm BEND}(\theta)=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\epsilon/k_{B}$ (close to zero), and the basic time step $\Delta t$ is reduced to $0.006\tau$. The friction constant is set to $\Gamma=20\tau^{-1}$ during the first $10^{3}$ MD steps, and $\Gamma=0.5\tau^{-1}$ afterwards 5, 28, 19. Thus, chains straighten out when the bond springs try to reduce the average bond length in order to minimize the energy from $\ell_{b}\approx 0.964\sigma$ of the OPs to $\langle b_{PP}\rangle_{0}=0.31\sigma$ of the PPs for unperturbed polymer melts and $\langle b_{PP}\rangle_{\lambda}\approx 0.57\sigma$ for strongly deformed polymer melts at $\lambda\approx 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 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 $\theta_{j,j+5}$ between bonds ${\bf b}_{j}$ and ${\bf b}_{j+5}$ for $j=2$, $3$, $\ldots$, $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, $\theta_{j,j+5}>90^{o}$, $90^{o}\geq\theta_{j,j+5}>60^{o}$, $60^{o}\geq\theta_{j,j+5}>30^{o}$, and $30^{o}\geq\theta_{j,j+5}>15^{o}$ (Figure 1d). Comparing to the force pattern along the same PP, we define an entanglement point to be located at ${\bf r}_{j+2}$ if the curvature of a kink corresponding to 42 $\theta_{j,j+5}\geq 60^{o}$. 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\sigma$.

3 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}/\sigma=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 labelled by $i=1$, $2$, $\ldots$, $1000$. For the fully equilibrated melts chain conformations are characterized by the mean square radii of gyration $\langle R_{g}^{2}(N)\rangle_{0}/\sigma^{2}$ of approximately $221$, $438$, and $909$ for the OPs of chains, $\langle(R_{g}^{(PP)}(N))^{2}\rangle_{0}/\sigma^{2}\approx 209$, $424$, and $895$ for the PPs of chains. Correspondingly, the mean square end-to-end distances are $\langle R_{e}^{2}(N)\rangle_{0}/\sigma^{2}=\langle(R_{e}^{(PP)}(N))^{2}\rangle_{0}/\sigma^{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 $(\perp)$ to the stretching, rescaled to $1/3$ and $2/3$ of the estimates of $\langle R_{g}^{2}\rangle_{0}^{1/2}$ ($\langle(R_{g}^{(PP)})^{2}\rangle_{0}^{1/2}$) , respectively, right after deformation are shown in Figure 2. Under uniaxial elongation, we see that the over-all conformations of the OPs as well as of the PPs deform affinely. Namely, $\left[\langle R_{g,||}^{2}\rangle_{\lambda}/\langle R_{g,||}^{2}\rangle_{0}\right]^{1/2}=\left[\langle(R_{g,||}^{(PP)})^{2}\rangle_{\lambda}/\langle(R_{g,||}^{(PP)})^{2}\rangle_{0}\right]^{1/2}=\lambda$ and $\left[\langle R_{g,\perp}^{2}\rangle_{\lambda}/\langle R_{g,\perp}^{2}\rangle_{0}\right]^{1/2}=\left[\langle(R_{g,\perp}^{(PP)})^{2}\rangle_{\lambda}/\langle(R_{g,\perp}^{(PP)})^{2}\rangle_{0}\right]^{1/2}=\lambda^{-1/2}$. The corresponding scaling laws also hold for $\langle R^{2}_{e,||}\rangle$ and $\langle R^{2}_{e,\perp}\rangle$, 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\langle R^{2}(n)\rangle_{\lambda}/\langle R^{2}(n)\rangle_{0}$ ($C\langle R^{(PP)}(n)^{2}\rangle_{\lambda}/\langle R^{(PP)}(n)^{2}\rangle_{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=\lambda$ and $C=1/\lambda^{2}$ for the two components in the directions perpendicular and parallel to the stretching, respectively. Our data demonstrate that $n\approx 8.2N_{e}$ related to the chosen strain rate $\dot{\varepsilon}$ 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.

3.1 Conformational change of single chains in deformed melts

In the non-linear viscoelastic regime, the tube model predicts for the initial relaxation right after deformation an over-damped initial retraction process of the individual chains in both directions parallel and perpendicular to the stretching direction 2, 44. Such an immediate chain retraction mechanism has been observed for the time evolution of the rescaled two components of radius of gyration for the OPs of chains 42, 43 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/\tau_{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, $\langle R_{g,||}^{(PP)}(t)^{2}\rangle$ and $\langle R_{g,\perp}^{(PP)}(t)^{2}\rangle$ decrease initially with increasing time showing the signature of chain retraction. $\left[\langle(R_{g,\perp}^{(PP)})^{2}\rangle_{\lambda}/\langle(R_{g,\perp}^{(PP)})^{2}\rangle_{0}\rangle\right]^{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 $(\langle R_{g,\perp}^{(PP)}(t)^{2}\rangle/\langle R_{g,\perp}^{(PP)}(t)^{2}\rangle_{0})^{1/2}$ for PPs also occur slightly later than for OPs 42, but the difference for the occurrence times becomes negligible within fluctuation as the chain size increases (at $t/\tau_{R,N}\approx 0.30(4)$ (OPs), $0.32(4)$ (PPs) for $N=2000$, at $t/\tau_{R,N}\approx 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, 44, i.e. the Rouse time of the whole chains in unperturbed melts, $\tau_{R,N}$. This indicates a rather strong contribution from tension along the PPs. Setting the two parameters $c_{\nu}=0.1$ and $R_{s}=2.0$ in the GLaMM model 44 (see the supplementary material in Ref. 43), 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 43. Obviously, our results for both the OPs and PPs qualitatively capture the signature of the initial chain retraction mechanism 2, 44 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 $(\langle R_{g,||}^{(PP)}(t)^{2}\rangle/\langle R_{g,||}^{(PP)}(t)^{2}\rangle_{0})^{1/2}$ becomes smaller with time for $N=1000$ and $2000$, and eventually $(\langle R_{g,||}^{(PP)}(t)^{2}\rangle/\langle R_{g,||}^{(PP)}(t)^{2}\rangle_{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 42, 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 61, 62, 63.

3.2 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 $\theta_{j,j+5}\geq 60^{o}$ (significant kinks) and the tension force ${\bf F}_{j,||}^{(PP)}=-{\bf\bigtriangledown}U_{\rm FENE}^{(PP)}\cdot\hat{e}_{x}$ 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_{\rm 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 $\lambda\approx 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_{\rm 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_{\rm 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 supplementary information of Ref. 42. 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\tau_{R,N=2000}$ ($N_{\rm kink}^{(i=2)}=52$) or even up to $1.0\tau_{R,N=2000}$ ($N_{\rm 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/\tau_{R,N}=2.0$, the inhomogeneous pattern still persists although the number of kinks ($N_{\rm 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_{\rm 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 64, 65 (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 65. 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}\approx 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, $\langle N_{\rm kink}\rangle_{\lambda}$, 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 perfect circle, but can be rather elliptical under deformation (see Figure 2b). Comparing to the estimates of $\langle N_{\rm kink}\rangle_{\lambda}$ 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 increase of strain $\lambda$. 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 $\langle N_{\rm kink}(t)\rangle$ shown in Figure. 5b is very similar to the retraction of the overall chain (see Figure 3b). $\langle N_{\rm kink}(t)\rangle$ initially decreases, reaches a minimum at $t/\tau_{R,N}\approx 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 $\langle N_{\rm kink}\rangle_{0}$. Taking the data for $N=500$ this final relaxation seems to occur on the time scale of the disentanglement time $\tau_{d}\cong\tau_{R,N}(N/N_{e})^{x}$, $x=1$ in the original reptation scheme 1, 2 and $x\cong 1.4$ experimentally 66.

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 melts is close to $N_{e}$. The number of bonds along such an almost straight segment of the PP between the $k$th and $(k-1)$th entanglement points of a chain is denoted by $l_{{\rm str},k}$,and the average length of the straight segments along the PP of each chain is given by $l_{\rm str}=\sum_{k}l_{{\rm str},k}/(N_{\rm kink}+1)$ with $\sum_{k}l_{{\rm str},k}=N-1$, where the two ends of the chain are treated as kinks. Thus, results of $\langle l_{{\rm str}}\rangle$ are expected to be proportional to $1/\langle N_{\rm kink}\rangle$. The distributions of $l_{\rm str}$ and $l_{{\rm str},k}$, $P(l_{\rm str},t)$ and $P(l_{{\rm str},k},t)$, respectively right after deformation, and upon relaxation to $t/\tau_{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 $\langle l_{{\rm str},k}\rangle\approx 22$ which is close to $N_{e,PPA}^{(0)}=28$. We see that the profiles of $P(l_{\rm str},t)$ are quite well described by shifted Gaussian distributions. The peak first shifts to significantly larger values of $l_{\rm str}$, while only after $t/\tau_{R,N}>0.5$ the peak of $P(l_{\rm str},t)$ slowly begins to move towards smaller values of $l_{\rm str}$, still being far away from the unperturbed case, even for $t\approx 1.5\tau_{R,N}$. Furthermore the width of the distribution $(\langle l_{\rm str}^{2}\rangle-\langle l_{\rm str}\rangle^{2})$ seems to become even wider with time up to $t\approx 0.5\tau_{R,N}$. Within the chains the distribution of the lengths of straight PP segments $P(l_{{\rm str},k},t)$, Figure 6b, shows long tail distributions of $l_{{\rm str},k}$. Right after deformation, the tail of $P(l_{{\rm 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\approx 0.5\tau_{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_{{\rm 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 $\lambda\approx 5.0$, we see only a slight shift to larger values of $l_{\rm str}$ for $P(l_{\rm str},t)$ compared to that of the unperturbed polymer melt since $\langle N_{\rm kink}\rangle$ is only about $10\%$ smaller than $\langle N_{\rm kink}\rangle_{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_{{\rm 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 $\langle N_{\rm kink}(t)\rangle$, $P(l_{\rm str},t)$, and $P(l_{{\rm str},k},t)$ indicate that topological constraints play an essential role in the relaxation retardation of deformed chains.

3.2.1 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, $\rho_{\rm kink}(x,y)$. $\rho_{\rm kink}(x,y)$ is estimated by simply counting the number of entanglement points located at $(x,y)$ and normalized such that $\int_{0}^{L_{x}}\int_{0}^{L_{y}}\rho_{\rm kink}(x,y)dxdy=1$. For this the grid spacing is set to $2.0\sigma$. For a large disorder system, structural inhomogeneities presented in such 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 ($\lambda=1.0$), $N_{\rm kink}=88128$ for $60^{o}<\theta_{j,j+5}<90^{o}$ and $N_{\rm kink}=7617$ for $\theta_{j,j+5}>90^{o}$ while $N_{\rm kink}=7920$ for $60^{o}<\theta_{j,j+5}<90^{o}$ and $N_{\rm kink}=74344$ for $\theta_{j,j+5}>90^{o}$ for the polymer melt right after deformation ($\lambda\approx 5.0$). Upon relaxation, we observe a clustering of entanglement points and the clustering pattern becomes more distinct and seems to reach a maximum around $t=0.5\tau_{R,N}$. Even up to $t=2.0\tau_{R,N}$, these clusters persist although the kinks along deformed chains become less sharp ($N_{\rm kink}=46917$ for $60^{o}<\theta_{j,j+5}<90^{o}$ and $N_{\rm kink}=8888$ for $\theta_{j,j+5}>90^{o}$). 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^{\rm(kink)}(q_{||})$ of entanglement points in the melt gives an impression of the correlation, where $q_{||}$ is the component of the vector ${\bf q}$ in the direction parallel to the stretching. The discretization in $q$-space along the stretching direction is given by $q_{||}=2\pi m_{x}/L_{x}$ with $m_{x}=0,1,2,\ldots$, where $L_{x}=L_{0}$ and $L_{x}\approx 5.0L_{0}$ before and after deformation, respectively. Results of $S^{\rm(kink)}(q_{||})$ plotted versus $q_{||}$ are shown in Figure 8. Here we only focus on the regime $q_{||}\leq 1.0\sigma^{-1}$. At first glance, the distributions of entanglement points seem to be more structured during initial relaxation after polymer melts are deformed comparing 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_{\rm cluster}$ can be roughly determined by $d_{\rm cluster}=2\pi/q_{||}^{\rm max}$, i.e., at $q_{||}=q_{||}^{\rm max}$, $S^{\rm(kink)}(q_{||})$ reaches a maximum along the stretching direction. For $t=0.25\tau_{R,N}$, we see one broad maximum around $q_{||}^{\rm max}\approx 0.08\sigma^{-1}$ which gives $d_{\rm cluster}\approx 78\sigma\sim 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\tau_{R,N}$, and $d_{\rm cluster}\approx 40\sigma\sim 8d_{T}^{(0)}$ ($q_{||}^{\rm max}\approx 0.16\sigma^{-1}$). Thus, the clustering of the entanglement points is related to the relaxation retardation of the deformed polymer melt.

3.3 Mobility of monomers in non-equilibrium 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 $\tau_{0}$ for local fluctuations, the entanglement time $\tau_{e}\propto N_{e}^{2}$, the Rouse time $\tau_{R,N}\propto N^{2}$, and the disentanglement time $\tau_{d,N}\propto N^{3.4}$ (when contour length fluctuations, constraint release, and correlation hole effects are taken into account 36, 67, 68, 37). 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 19, 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, 19) parallel and perpendicular to the stretching direction, $g_{1,||}(t)$ and $g_{1,\perp}(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)\sim t^{1}$ for $t<\tau_{0}$, and $g_{1,||}(t)\sim t^{1/2}$ for $\tau_{0}<t<\tau_{e}=\tau_{0}(N^{(0)}_{PPA})^{2}$) For $t>\tau_{e}$, monomers move faster, $g_{1,||}(t)\sim t^{0.69}$, until reaching $t\approx 0.30\tau_{R,2000}$, i.e., the duration of the chain retraction process. Within this time frame, the change of $g_{1,||}(t)$, $\Delta g_{1,||}(t)\approx{\cal O}(10^{3})$ is of the same order as the change of $\langle R_{g,||}^{2}(t)\rangle$, while the center of mass displays unperturbed diffusion $g_{3}(t)\sim t$ up to $0.30\tau_{R,2000}$, but moves a little more slowly around $t=\tau_{e}$. For $t>0.3\tau_{R,2000}$, monomers motion is slowed down again due entanglement effects and the corresponding relaxation retardation, similar as in a tube-like regime created by the surrounding chains, $g_{1,||}(t)\sim 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)\sim t^{1/4}$. Perpendicular to the stretching, the situation is somewhat different. Apparently, monomers move more slowly that a gradual deviation from $g_{1,\perp}\sim t^{1/2}$ to $g_{1,\perp}\sim t^{1/4}$ develops as $t$ approaches $\tau_{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=\tau_{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/\tau_{R,2000}$ are shown in Figure 9d. We see that monomers in the constrained regime move much slower comparing 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 as we have expected due to the weak topological constraints at the end and the retraction process.

3.4 Intrinsic properties of primitive paths and the tube picture

Based on the tube picture 1, 9, 2, one would expect that entanglement effects appear at $t\approx\tau_{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 average contour length of PPs is defined by $L_{PP}=(N-1)\langle b_{PP}\rangle$ with the average bond length of PPs

[TABLE]

where ${\bf r}_{i,j}^{(PP)}$ is the coordinate of the $j$th 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 unit sphere, normalized by $4\pi$, the average bond length as a function of the strain $\lambda$ is given by

[TABLE]

with

[TABLE]

where ${\bf E}$ is a diagonal deformation tensor with elements $\{\lambda,\lambda^{-1/2},\lambda^{-1/2}\}$, ${\bf u}$ is a unit vector in the spherical coordinate, ${\bf u}=(\cos\phi,\sin\phi\cos\theta,\sin\phi\sin\theta)$, $\phi$ is the polar angle between ${\bf u}$ and the $x$-axis, and $\theta$ is the azimuthal angle. One gets $\langle b_{PP}\rangle_{\lambda}\approx 2.56\langle b_{PP}\rangle_{0}$ at $\lambda=5.0$. The deviations between different approximations discussed in our previous work 42 are a consequence of the fact that there is a lower cutoff, given by $N_{e}$. In general one should expect that $\langle b_{PP}\rangle_{\lambda}/\langle b_{PP}\rangle_{0}$ for the affinely elongated PP mesh follow eq. 12, 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 $\lambda$. At $\lambda=5.0$, $\langle b_{PP}\rangle_{\lambda}\approx 1.84\langle b_{PP}\rangle_{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, $\langle b_{PP}(t)\rangle$ reaches $\langle b_{PP}\rangle_{0}$ for $N=500$ and $N=1000$ for $t>\tau_{R,N}$ while for $N=2000$, $\langle b_{PP}(t)\rangle$ seems to settle at a slightly larger value than $\langle b_{PP}\rangle_{0}$.

Choosing the $x$-axis (along the stretching direction) as reference, the orientational order parameter of bond vectors along PPs is defined by $Q_{PP}(\phi)=(3\cos^{2}\phi-1)/2$ where $\phi$ is the angle between the bond vector ${\bf b}_{PP}$ and the $x$-axis. Assuming the distribution of the original bonds to be isotropic, $\langle Q_{PP}(\phi)\rangle_{\lambda}$ for the PP meshes of polymer melts under elongation is given by

[TABLE]

Under elongation Figure 10c shows that the orientation distribution of bonds changes from isotropic, $\langle Q_{PP}(\phi)\rangle_{\lambda}=0.0$ for $\lambda=1.0$, to anisotropic, $\langle Q_{PP}(\phi)\rangle_{\lambda}\approx 0.8$ for $\lambda=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 $\langle b_{PP}(t)\rangle$ and $\langle Q_{PP}(\phi,t)\rangle$ show relaxation retardation starting at $t\approx\tau_{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, is different in different directions since the contours of PPs globally deform affinely. Keeping this in mind, we nevertheless apply the standard formula $\langle R_{e}^{2}\rangle=\langle(R_{e}^{(PP)})^{2}\rangle=L_{PP}\ell_{K}^{(PP)}=(N-1)N_{e,PPA}\langle b_{PP}\rangle^{2}$, $\ell_{K}=N_{e,PPA}\langle b_{PP}\rangle$ 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 $\langle R_{e}^{2}\rangle_{\lambda}=\lambda^{2}\langle R_{e,x}^{2}\rangle_{0}+\langle R^{2}_{e,y}+R^{2}_{e,z}\rangle_{0}/\lambda=(\lambda^{2}+2/\lambda)\langle R_{e}^{2}\rangle_{0}/3$. Following eq. 12, we obtain

[TABLE]

Figure 11a shows that $N_{e,PPA}^{(\lambda)}$ for the elongated PP mesh follows affine deformation up to $\lambda\approx 5.0$ while for the PPs of the deformed polymer melts it does not. This is consistent with the estimates of “effective entanglement length” $\langle l_{\rm str}\rangle_{\lambda}\propto 1/\langle N_{\rm kink}\rangle_{\lambda}$ (Figure 5a). At subsequent relaxation time $t\approx 0.5\tau_{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=\tau_{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 $d_{T}=d_{T}^{(0)}/\sqrt{\lambda}$ , $\langle N_{\rm tube}\rangle$ 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 the fully equilibrated melt, namely to 19 $d_{T}^{(0)}\approx\sqrt{2}\langle R_{g}^{2}(N_{e})\rangle_{0}^{1/2}\approx 5.02\sigma\propto\sigma\sqrt{N_{e,PPA}^{(0)}}$, we follow the variation of $N_{\rm tube}^{(i)}$ for each chain $i$ and $\langle N_{\rm tube}\rangle$ upon deformation and subsequent relaxation. The way of counting $N_{\rm 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 the shortest distance to any PP bond $j$ with $k-\delta\leq j\leq k+\delta$ is less or equal to $d_{T}^{(0)}/2$ with $\delta=8$, being determined by $d_{T}^{(0)}$ ($\delta\langle b_{PP}\rangle_{0}\sim d_{T}^{(0)}/2$). Note that the precise number of $N_{\rm tube}$ depends on details of the algorithm and the chosen value for $\delta$, but results of $\langle N_{\rm tube}\rangle$ 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 show approximate agreement with to 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 tube-like regime along its corresponding PP as the chain is strongly stretched, and thus $N_{\rm tube}^{(i)}$ increases. Once the deformed chain starts to relax, $N_{\rm 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.

$\langle N_{\rm tube}\rangle$, cf. Figure 13ab, fits to the overall scheme observed so far. In the linear regime ($\lambda<1.5$), $\langle N_{\rm tube}\rangle_{\lambda}\approx\langle N_{\rm tube}\rangle_{0}$. In the non-linear regime ($\lambda>1.5$), upon deformation $\langle N_{\rm tube}\rangle_{\lambda}$ increases with a rate larger with larger $N$. Right after deformation, at the subsequent initial relaxation time, $\langle N_{\rm tube}(t)\rangle$ decreases, if normalized by the Rouse time. Eventually, it even drops below the equilibrium value, slows down for $t/\tau_{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^{(\lambda\approx 5.0)}_{e,PPA}/N^{(0)}_{e,PPA})^{1/2}\approx 7.61\sigma\approx 1.5d_{T}^{(0)}$, the profiles of $\langle N_{\rm tube}(t)\rangle/\langle N_{\rm tube}\rangle_{0}$ having the minimum around $t\approx\tau_{R,N}$ within fluctuation are similar as the results shown in Figure 13b, while quantitatively, the estimates of $\langle N_{\rm tube}(t)\rangle$ are somewhat larger. Again this agrees with the formation and growth of topologically highly congested areas along the chains. Regions of low density of entanglement points, where configurations of monomers along the OPs can fluctuate significantly in space, seem to stabilize regions with high density of entanglement points. The probability distributions $P(N_{\rm tube}/N)$ as a function of $(N_{\rm tube}/N)$ for deformed polymer chains of size $N=2000$ in a melt at several selected strain values of $\lambda$ and at several selected rescaled relaxation times $t/\tau_{R,N}$ after deformation are shown in Figure 13c,d. We see that the $P(N_{\rm tube})$ is simply a shifted Gaussian distribution in terms of the mean value $\langle N_{\rm tube}\rangle$, and the standard deviation $\sigma(N_{\rm tube})=\sqrt{\langle N^{2}_{\rm tube}\rangle-\langle N_{\rm tube}\rangle^{2}}$. The distribution $P(N_{\rm tube}/N)$ remains the same for $\lambda=1.0$ and $\lambda\approx 1.5$, and then the profile of $P(N_{\rm tube}/N)$ shifts to a larger value of $N_{\rm tube}/N$ as $\lambda$ increases. After deformation, the profile shifts to smaller value of $N_{\rm tube}/N$ with increasing $t$, and even moves to the left-hand side of the profile for unperturbed chains similar as the behavior of $\langle N_{\rm tube}(t)\rangle$ shown in Fig. 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.

3.5 Stress relaxation

After having discussed details of individual and collective conformational relaxation we turn to the related stresses in the systems, a quantity which 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)\sim t^{-1/2}$ for $t<\tau_{e}$, while for $\tau_{e}<t\ll\tau_{d,N}$ ($\tau_{d,N}=\tau_{R,N}(N/N_{e})^{1.4}$ being the disentanglement time of chains of size $N$) where chains are assumed to move in a tube-like regime due to entanglements, $G(t)$ reaches a plateau value $G_{N}^{0}=(4/5)(\rho 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 9, 2. 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 non-linear viscoelastic regime and the relationship between the strain rate and stress relaxation is significantly less developed.

Results of the stress relaxation modulus $G(\lambda,t)$ that characterize the viscoelasticity of entangled polymer melts in both linear (unperturbed and $\lambda\approx 1.2$) and non-linear regimes ($\lambda\approx 5.0$) are shown in Figure 14 (parts of data have been shown in Refs. 19, 42). For clarity only the case of $N=2000(\approx 72N_{e})$ is discussed here. For equilibrium polymer melts $G(\lambda=1.0,t)$ is calculated from the stress autocorrelation function of off-diagonal elements of the stress tensor 69, 70, 19, using the Green-Kubo relationship. $G(\lambda=1.0,t)$ reaches a plateau value $G_{N}^{0}$ with $N_{e}\approx 28$ for $\tau_{d,N}\gg t>\tau_{e}\approx 2\times 10^{-4}\tau_{R,N}$. Alternatively $G(\lambda,t)$ can also be given by the normal stress difference after the simulation box is deformed, i.e. stretched along the $x$-direction,

[TABLE]

where $\sigma_{norm}(\lambda,t)=\sigma_{xx}-\frac{1}{2}(\sigma_{yy}+\sigma_{zz})$. For small and fast deformation ($\lambda\approx 1.2$ and $\dot{\varepsilon}\tau_{R,N}=32000$), $G(\lambda\approx 1.2,t)$ follows $t>\tau_{e}$ the data 19 calculated from the Green-Kubo relation, and in good agreement with experiment 71 and the tube model 1, 9, 2. For the large but relatively slow deformation ($\lambda\approx 5.0$ and $\dot{\varepsilon}\tau_{R,N}=77$), $G(\lambda\approx 5.0,t)$ already from the very beginning displays a plateau modulus $G_{N}^{0}$ below that of the linear regime. This is consistent with the experimental finding that one observes shear thinning in melt under such conditions 72, 73. Only after a time corresponding to the inverse strain rate, i.e. $t=\dot{\varepsilon}^{-1}=\tau_{R,N}/77\approx 10^{-2}\tau_{R,N}$ a further decay of $G(t)$ is observed. Keeping $\lambda\approx 5.0$, but choosing the same very fast strain rate $\dot{\varepsilon}\tau_{R,N}=32000$ as for $\lambda\approx 1.2$, we see that $G(\lambda\approx 5.0,t)$ initially follows the data of $G(\lambda\approx 1.2,t)$ up to $\dot{\varepsilon}^{-1}=\tau_{R,N}/32000\approx 10^{-5}\tau_{R,N}$. Then it deviates and after the relaxation time reaches $t=\dot{\varepsilon}^{-1}=\tau_{R,N}/77$, both curves for $\lambda\approx 5.0$ follow the same softening pattern. This indicates that relevant chain chain interpenetration did not change significantly during the slow stretching process ($\dot{\varepsilon}=77/\tau_{R,N}$ which is still quite fast comparing to $\tau_{d,N}^{-1}$), that is, topological constraints are not released up to some chain end effects. Furthermore this supports the concept that the long time behavior of the stress relaxation also in the non-linear regime seems to be just a function of the ultimate deformation and is independent of the original strain rate $\tau^{-1}_{R,N}<\dot{\varepsilon}<\tau^{-1}_{e}$. From the above we expect that estimates of $\langle N_{\rm kink}(t)\rangle$ for deformed polymer melts with two different strain rates at $\lambda\approx 5.0$ in Figure 14b should eventually coincide. Indeed, as expected for $t/\tau_{R,N}>1/77\approx 0.013$, the two sets of data perfectly coincide with each other, i.e. entanglement effects in both deformed polymer melts remain the same.

4 Conclusions

In this paper, we have employed extensive molecular dynamics simulations to study highly entangled polymer melts subject to an isochoric elongation in the non-linear viscoelastic regime. Focusing on the analysis of topological constrains, as identified through a primitive path analysis 4, 5, 42, 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 $\tau_{R,N}^{-1}<\dot{\varepsilon}<\tau_{e}^{-1}$ (Figures 2). The chain retraction mechanism predicted by the Doi-Edward tube model 2 and its refined GLaMM tube model 44 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, 44. 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 affinely upon an uniaxial elongation in the non-linear viscoelastic regime ($\lambda>1.5$). At subsequent relaxation times, the average number of kinks first decreases, and reaches its minimum value before it moves towards 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 center of mass right after polymer melts are deformed, significant deviations from that for unperturbed polymer melts are discovered for $t>\tau_{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. Following the tube picture assumption, we also count the number of monomers confined in a tube-like 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 44. Finally, we see that the time evolutions of the stress relaxation modulus for deformed polymer melt only initially depends on the strain rates while they follow the same softening patterns ultimately in the non-linear 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 non-linear 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 for 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. 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 offers interesting options for both experimental and simulation investigations, especially in the vicinity of the glass transition point. 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.

{acknowledgement}

We are grateful to M. Doi, G. S. Grest, D. Vlassopoulos, T. Ohkuma for stimulating and helpful discussions, and K. Ch. Daoulas for a critical reading of the manuscript. This work has been supported by European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC Grant Agreement No. 340906-MOLPROCOMP. We also gratefully acknowledge the computing time granted by the John von Neumann Institute for Computing (NIC) and provided on the supercomputer JUROPA at Jülich Supercomputing Centre (JSC), and the Max Planck Computing and Data Facility (MPCDF).