Local Transient Jamming in Stress Relaxation of Bulk Amorphous Polymers

Bulk amorphous polymers become stretched and parallel-aligned under loading stress, and their intermolecular cooperation slows down the subsequent stress relaxation process. By means of dynamic Monte Carlo simulations, we employed the linear viscoelastic Maxwell model for stress relaxation of single polymers and investigated their intermolecular cooperation in the stress relaxation process of stretched and parallel-aligned bulk amorphous polymers. We carried out thermal fluctuation analysis on the reproduced Debye relaxation and Arrhenius fluid behaviors of bulk polymers. We found a transient state with stretch-coil coexistence among polymers in the stress relaxation process. Further structure analysis revealed a scenario of local jamming at the transient state, resulting in an entropy barrier for stretch-coil transition of partial polymers. The microscopic mechanism of intermolecular cooperation appears as unique to polymer stress relaxation, which interprets the hydrodynamic interactions as one of essential factors raising a high viscosity in bulk amorphous polymers. Our simulations set up a platform of molecular modeling in the study of polymer stress relaxation, which brought new insights into polymer dynamics and the related mechanical/rheological properties.


INTRODUCTION
Stress relaxation is a kinetic process of liquids that releases the internal tension without changing the global strain, appearing as a fundamental behavior in the thermal fluctuations and the mechanical properties of condensed matter. The most basic linear viscoelastic behavior of stress relaxation can be described by the Maxwell model consisting of a series connection of a viscous dashpot and an elastic spring, [1] as illustrated in Fig. 1. One of the significant features for chainlike polymers different from small molecules is the substantial slowing down of stress relaxation due to a high viscosity in bulk amorphous polymers. [2,3] Among those factors raising the high viscosity, the intramolecular chain-connectivity of monomers has been well recognized by the Rouse model of polymer dynamics. [4] In addition, the effect of long-range intermolecular entanglement has been considered by the reptation model of polymer dynamics. [5] However, so far, the effect of short-range intermolecular cooperation of flexible polymer chains has commonly been treated as a frictional constant averaged over all the monomers in the fully draining mode, appearing as the meanfield hydrodynamic interactions in polymer melt, [6] and its microscopic mechanism unique to polymer stress relaxation remains unclear.
The stress loading in bulk amorphous polymers commonly leads to parallel alignment of orientated polymer chains. In this sense, we employed the Maxwell model for stress relaxation of each polymer and investigated the intermolecular cooperation upon stress relaxation in stretched and parallelaligned bulk amorphous polymers, by means of dynamic Monte Carlo simulations of athermal lattice polymers. The simulations could reproduce the Debye relaxation and Arrhenius fluids behaviors of bulk polymers. This linear response condition allows us to make thermal fluctuation and structure analysis on the transient state of the stress relaxation process, which will reveal a stretch-coil coexistence among polymers, indicating a scenario of local transient jamming upon intermolecular cooperation that raises an entropy barrier for polymer stress relaxation.

MODELING STRESS RELAXATION Preparing the Initial State
We put an array of 1920 regularly folded polymer chains, each holding 128 monomers with 8 stems in the folding length 16 (occupying the consecutive lattice sites), in parallel stacking side-by-side along the YZ plane in the lattice space of cubic entities 16×128×128 (X×Y×Z), and let them relax into the bulk amorphous coils under athermal conditions with neutrally repulsive hard walls at the boundaries while their two chain ends were not allowed to leave the corresponding YZ planes at the boundaries of X-axis. Therefore, polymers remained in their less entangled and parallel-stacked states in order to avoid the possible influence of chain entanglement in the present case. During relaxation, polymer chains performed micro-relaxation steps in the lattice space via single monomer jumping into its vacancy neighbor, sometimes allowing local sliding diffusion along the chain. [7] Both the double monomer occupation and the bond crossing were not allowed in the micro-relaxation steps to mimic the volume exclusion of polymers. This lattice model of polymers belongs to the single-site bond-fluctuation model, which could reproduce the time scaling law of the Rouse model in bulk amorphous polymers as well as the chain-length scaling law of coil sizes in polymer solutions. [7] Then, keeping the two chain ends separately on the parallel boundary planes, 1920 bulk amorphous polymers were stretched step-by-step along X-axis up to 100% strain. In each step, we selected an Xsite randomly, let YZ-plane cut the body into two halves, moved the right body for one more X-site, and reconnected all the broken ends on that cutting YZ-plane via local sliding diffusion in the rest body. [8] We then gave 4000 MC cycles of athermal relaxation to remove the new vacancy sites, which recovered the total volume of bulk polymers upon stretching. Here, the MC cycle was defined as the time unit on one trial move of micro-relaxation per monomer averaged over the total amount of monomers. During stretching, the lateral free-standing surfaces were maintained by the repulsion of those new vacancy sites from bulk polymers. The scheme of homogeneous stretching has been broadly applied in the simulations of straininduced crystallization in homopolymer, [8,9] random copolymers, [10,11] polydisperse polymers, [12] solution polymers, [13] blend polymers [14] and short-branch polymers. [15] After the sample system has reached 100% strain upon the above step-by-step stretching, one million MC cycles were measured for additional athermal relaxation to ensure the equilibrated states in order to eliminate the kinetic effect of the stretching with a limited strain rate. As demonstrated in Fig. 2, we used the final state of 100% strain as the initial state for our present simulations of stress relaxation.
Concerning to the lateral surfaces that might affect our results for bulk polymers, we calculated X-fraction of meansquare radius of gyration of those polymer coils with their mass centers distributed on various radius distances away from the center of their YZ-planes, in the initial state of bulk polymers for stress relaxation, as shown in Fig. 3. One can clearly see the decaying of coil sizes near the lateral surface of the cylindrical space of bulk polymers, indicating the skin layer of our sample system. Those polymers with their mass centers locating in the skin layer were not included in our following observations of stress relaxation. In other words, we paid attention only to those polymers locating in the central bulk amorphous region to eliminate the influence of the lateral surfaces.

Stress Relaxation Model of Single Polymer
From the initial state of 100% strain, we removed the restrictions on both the chain-ends and launched the stress relaxation Snapshots to demonstrate a homogeneous deformation (100% strain) of bulk polymers along X-axis, stretched with the strain rate of 6.25%/4000 MCCs under athermal conditions. The blue cylinders in the box represent bonds, and the yellow spheres represent the chain ends that were restricted at two YZ planes at the boundaries of X-axis. The left snapshot shows the initial state of deformation, and the right one shows the state of 100% strain as the initial state of stress relaxation.
process. Since 100% strain was not large for our polymers, and the lateral free-standing boundary held strong repulsion of those new vacancy sites, the global strain could remain constant during the whole stress relaxation process. Therefore, the residual modulus M(t) could be used to monitor the residual stress. We assigned the initial modulus M 0 to the internal tension of each polymer coil, and let it decay in the time (t) evolution by following Debye relaxation [16] with a characteristic time τ as given by On the directions of the stretching X-axis, we assumed that the probability to find the monomers around the mass center of each polymer coil roughly followed the Gaussian distribution function, as illustrated in the upper left corner of Fig. 4. Every monomer intended to come back to its most probable position at the mass center due to the elasticity of polymer conformational entropy, but this tendency was balanced by the internal tension as delivered from the external loading stress. The internal tension slowed down the stress relaxation in each coil, rather than being driven solely by the restoring of the conformational entropy. Therefore, according to the Boltzmann's relationship, each monomer held an elastic energy 1/2 × M × x 2 , and its derivation to x gave the elastic force M × x, while M was the modulus of internal tension in balance with the entropy elasticity of polymer coils and x was the monomer distance from the coil mass center along the stretching X-axis. In practice, we assigned an internal tension F on each monomer with its strength depending upon its location and its distance relative to the mass center of the coil on the stretching X-axis, as illustrated in Fig. 4.
Next, we employed the Maxwell model to describe the microscopic stress relaxation of each polymer coil in our simula-tions. The Debye relaxation of the residual modulus was synchronized with the Brownian motions of each polymer coil. Therefore, in Eq. (1), the time t was replaced by the diffusion time of the coil mass center, according to the Einstein's relationship [17] as given by where [r(t)-r(0)] 2 was the square displacement of the mass center of that polymer coil performing random Brownian motions within the time period t, and D was the diffusion constant of that polymer coil. The characteristic relaxation time τ was replaced as well by the characteristic time of that polymer coil diffusing over the square radius of gyration of the coil R g 2 , [18] as given by In practice, we chose the time segment as 200 MC cycles for each step of microscopic stress relaxation in each polymer coil. The evidence for the optimum selection of 200 MC cycles will be provided in the next subsection. Upon an integration of time segments along the path of random Brownian motions, the residual modulus M was delivered to the next segment as the initial modulus, and the coil size R g 2 at the beginning of the segment was used in the calculation of the segment. Segment-by-segment, we traced the residual modulus averaged over all the polymer coils at the end of each segment in order to monitor the macroscopic stress relaxation X-fraction of mean square radius of gyration <R g 2 > of polymers with their mass centers distributed on various radius distances away from the YZ-plane center in the initial state of the sample system for stress relaxation. One could clearly see the skin layer holding the decaying coil sizes of polymers from the central bulk amorphous region. The inset snapshot demonstrates the initial state of 100% strain polymers with their chain-bonds in blue color and chain-ends in yellow color. 2 , a, b>0

Fig. 4
Assuming a roughly gaussian distribution of monomers around the mass center on the stretching axis. The internal tension in each stretched polymer coil was assigned to the force F or F' that drags each monomer away from the mass center to resist the entropy elasticity in the polymer coil, with the strength depending upon its relative location (for example, downside) and distance (for example, d or d', respectively) to the mass center on the stretching axis (up and down directions). process.
In the above stress relaxation model of each polymer coil, the Brownian motions of polymer coils in the bulk amorphous phase slowed down the stress relaxation process in each polymer, and the resulted characteristic time could reveal the effect of intermolecular cooperation in the stress relaxation process.

Modeling Stress Relaxation in Bulk Amorphous Polymers
The conventional Metropolis sampling algorithm [19] was employed with the total potential energy change in each trial move, as given by where ∑x was summed over all the monomers involved in each step of micro-relaxation (x was negative if the monomer left away from the mass center along the directions of X-axis), M was calculated from the above procedure, T was the temperature, and k was Boltzmann's constant. We used kT/M 0 as the reduced temperature below. The internal tension M×x made the work in one lattice unit along the stretching X-axis to let the Brownian motions slow down the spontaneous retraction raised by the conformational entropy elasticity. Here, except for the volume exclusion and chain-connectivity, any other inter-and intra-chain interactions of monomers were not introduced, in order to expose the entropy side of intermolecular cooperation in the stress relaxation process. The enthalpy side of stress relaxation will be studied in the near future. The residual modulus averaged over all the polymer coils followed well with the Debye relaxation curves according to Eq. (1), as summarized in the proper temperature range in Fig. 5. Within the temperature window from 30 to 100, linear regression on the early-stage stress relaxation at 3.0×10 5 MC cycles could reach 0.999 level of R-square values. We checked various lengths of time segments for the quality of linear regression in Fig. 5 at three typical temperatures 30, 40 and 50. As shown in Fig. 6, above 200 MC cycles, the quality of linear regression becomes saturated. The time segment could be too short to ruin the characteristic feature of Brownian motions, [20] or too long to lower the simulation efficiency. Thus, 200 MC cycles of time segment was the optimum selection.
According to the Maxwell model, the characteristic time τ α obtained from the slopes in Fig. 5 could be regarded as equal to the characteristic time of polymer diffusion. Fig. 7 shows that within the temperature window from 30 to 80, they follow the Arrhenius equation [21] well as given by where ε was the activation barrier for polymer diffusion (also called diffusion energy) in the bulk amorphous phase. Since MC simulations could reproduce the Brownian motions provided by the broad enough time segment (here, 200 MC cycles), [20] fast diffusion at very high temperatures will ruin the validity of MC simulations in the time segment of 200 MC cycles, resulting in the deviation from Arrhenius fluid behaviors near the high temperature end. On the other side, at the low temperature end, the deviation of Arrhenius fluid behaviors could be attributed to more complicated mechanism of intermolecular cooperation, which will be further studied in the near future.

RESULTS AND DISCUSSION
We focused on the Debye relaxation of Arrhenius fluids of bulk polymers in the temperature region from 30 to 80. In this linear response region where the global stress relaxation process of bulk amorphous polymers reproduces the microscopic input of Debye relaxation in each polymers, the fluctuation-dissipation theorem allows us to make a thermal fluctuation analysis on the time-correlation function of residual moduli among polymers during this non-equilibrium process of stress relaxation. [22] To this end, we monitored the time evolution of the  autocorrelation function Δ(t) of the residual moduli M(t)/kT in the stress relaxation process of polymers, as the equation [15] given by where < … > was the assemble average over three repeating simulations. The strong thermal fluctuation peaks mainly occur around 1.0×10 5 MC cycles for three temperatures 30, 40 and 50 during their stress relaxation processes, as shown in Fig. 8. The higher peaks at the lower temperatures indicate a higher kinetic barrier for stress relaxation. The peak time for the occurrence of the transient state at the peak top seems insensitive to the temperature. Apparently, the thermal fluctuations among the residual moduli of polymers are responsible for the observed characteristic relaxation time of stress relaxation.
In order to study the microscopic mechanism of the kinetic barrier for stress relaxation, we further analyzed the distributions of residual moduli among polymers at various stages of stress relaxation at temperature 30, as shown in Fig. 9. One can clearly see the double peak of modulus distributions at the transient state of 1.0×10 5 MC cycles. The transient bistable state indicates that spontaneous dynamic heterogeneity was raised among polymers in their stress relaxation processes. Some polymers are obviously relaxing faster than the others. Since the residual moduli reflect the relaxation extents of stretched polymers, the bistable state demonstrates the coexistence of stretches and coils among polymers. Similar to a nucleation barrier for first-order phase transition, the free energy barrier that causes the stretch-coil coexistence could be attributed to the conformational entropy loss of those stretched polymers due to intermolecular cooperation. Therefore, the first-order stretch-coil transition implies an en-tropy barrier for stress relaxation in bulk amorphous polymers.
de Gennes used the simple dumbbell model to interpret the abrupt stretch-coil transition in single flexible polymer under large velocity gradient of flows. [23] The free energy bar- rier for the stability of stretches was attributed to the recovery of polymer conformational entropy detained by the frictional interactions of the flow. [18] Such a transition retards the relaxation of polymer chains in the shear and elongational flows of dilute solutions, as evidenced by the observations of video fluorescence microscopy in Chu's group [24,25] as well as by Brownian dynamics simulations in Larson's group. [26] So far, in the concentrated and bulk amorphous polymer phases even under fast flows, the large-scale phenomenon of disrupt stretch-coil transition has not yet been observed by the scattering experiments.
Here similarly, the entropy barrier for stretch-coil transition is apparently raised by the intermolecular cooperation during the stress relaxation process of bulk polymers. We made the structural analysis on the transient state of stretch-coil coexistence to study its microscopic mechanism. To this end, we tentatively split all the polymers into two fractions separately corresponding to those in the coil states and those in the stretch states, as marked by the dashed vertical line at M/kT=0.01 on the curve of 1.0×10 5 MC cycles for the stress relaxation process at the temperature 30 shown in Fig. 9. We calculated the pair correlation functions of the monomers along three axes for the two fractions of stretches and coils, as for instance along X-axis defined by [22] where L y and L z were the boundary sizes on the YZ plane, Δx was the lattice unit distance, ρ was the occupation density, N was the total amount of monomers in the stretch or coil fractions of polymers, and x ij was the X-axis distance between monomer i and monomer j. As shown in Fig. 10, the overlapping of constant pair correlation functions for the monomers of the stretch polymers along Y-and Z-axes was observed clearly, implying that the monomers of the stretch polymers are homogeneously distributed over both Y-and Z-axes, with no sign of segregation between the stretches and the coils on the dimensions of the YZ plane, as the homogeneous blue-colored polymers shown in the left inset snapshot viewed along X-axis. This observation implies no long-range intermolecular correlation on the YZ plane for the aggregation of stretches, but just local events of stretch-coil coexistence during the stress relaxation process of polymers. The local events explain why the transient states of stretch-coil coexistence are difficult for the detection of scattering experiments upon the stress relaxation process of the concentrated or bulk amorphous polymers. In Fig. 10, both the pair correlation function curves for the monomers of the stretches and of the coils along X-axis indicated the structure formation due to their segregation at the transient state. The monomers of the coil polymers are concentrated in the central area of X-axis like spheres, while the monomers of the stretch polymers are rich near two borders like dumbbells, as demonstrated separately by the red-and blue-colored polymers shown in the right inset snapshot viewed along Y-axis. Apparently from the structural point of view, the transient state holds the events of spherical coils jamming in the central space of X-axis and the coils are surrounded with dumbbell stretches at two borders, which are local as homogeneously distributed over the YX plane.
We further traced back and forth on the history of the monomer distributions of the stretch and coil fractions on the stretching X-axis, as shown in Fig. 11. One can see that, at the very beginning of the stress relaxation process, those coil polymers hold slightly more monomers in the central space. One can imagine that at beginning, all the stretch polymers tend to relax and occupy the central space of X-axis upon stress relaxation, raising traffic jamming in the central space. Just by chance, those polymers randomly distributed over the YZ plane hold relatively more monomers in the central space, and exhibit a priority of stress relaxation because the rest monomers near two borders are relatively less and hence  easier to come into the central space; thus they become the dominant monomers filled in the limited central space along X-axis; in consequence, they temporarily block the stress relaxation of their neighboring stretch polymers, till some of them diffuse away from the central space via random Brownian motions. The inset schemes in Fig. 11 illustrate the spheredumbbell coexistence of coil (red) and stretch (blue) polymers at the transient state of stress relaxation from the stretched and parallel-aligned polymers. The scenario of local jamming on the kinetic paths from stretches to coils of parallel-aligned polymers demonstrates an intrinsic entropy barrier unique to the intermolecular cooperation of polymers upon stress relaxation in the bulk amorphous phase. Such a microscopic mechanism of intermolecular cooperation significantly slows down the stress relaxation process. Therefore, it becomes one of the essential factors to raise the high viscosity of bulk amorphous polymers.

CONCLUSIONS
We employed the microscopic Maxwell model for linear viscoelastic stress relaxation of single polymers to analyze our dynamic Monte Carlo simulations of stress relaxation of stretched and parallel-aligned bulk amorphous polymers. The average residual modulus could reproduce the Debye relaxation and Arrhenius fluid behaviors of bulk polymers. Under this linear response condition, we carried out fluctuation and structural analysis on the non-equilibrium process of stress relaxation. The results demonstrated a transient jamming state with local stretch-coil coexistence among the polymers. The microscopic scenario of stress relaxation appears to be imaginable: along with the loading stress, the strain of bulk amorphous polymers raises the orientation and parallel-alignment of polymer chains. Upon stress relaxation, every stretched polymer coil tends to relax and occupy the common but limited central space, which causes local traffic jamming and slows down the stress relaxation process in bulk amorphous polymers. In the beginning, some randomly distributed polymers occasionally hold slightly more monomers in the central space and thus relax in priority to generate spontaneous dynamic heterogeneity; as a result, they take over the limited central space and tentatively block the relaxation of their stretch neighbors. Such a mechanism of intermolecular cooperation makes a barrier originated from polymer conformational entropy for stress relaxation in bulk amorphous polymers. If without this jamming, polymer chains were supposed to gain a higher conformational entropy (or a lower free energy) for a period of time in the transient state.
Our observations revealed the polymer feature of intermolecular cooperation that joins together with other polymer features like chain connectivity and chain entanglement to raise the high viscosity of bulk amorphous polymers.
Our present work sets up a platform of molecular modeling on polymer stress relaxation and paves the way towards a better understanding of polymer dynamics and the related mechanical/rheological properties. In the next steps, we will study how those chain-unit interactions dominate polymer stress relaxation and how the Arrhenius fluids of bulk polymers lose their linear viscoelasticity at low temperatures. We even add stress relaxation into our previous MC simulations of strain-induced polymer crystallization to study the mechanisms of stress-induced polymer crystallization in mono-and bi-axial stretching processes.