Molecular models for creep in oriented polyethylene fibers

High-performance polyethylene (PE) fibers are highly oriented polymer materials with stiffness comparable to steel while having one eighth the density.^1–3 Polyethylene’s simple chemical structure and controllable semicrystalline morphology makes it an affordable and versatile material for applications including ship sails, mooring lines, vehicular chassis, prosthetic implants, composite armors, and fabrics.^4,5

Fibers are typically formed by plastically drawing a melt or solution of entangled linear chains. Drawing orients and extends chains along the fiber axis, creating extended chain crystals that compose more than 90% of the final fiber.^6,7 The crystal domains are oriented with chain backbones along the fiber axis and interconnected to form a mesoscale network spanning the fiber.

Chains in PE fibers are not permanently cross-linked as in rubbers or elastomers. Thus, most fibers creep during sustained loading.⁸ Creep causes fibers to fail at stresses well below their ultimate strength, severely limiting their potential applications. Understanding the microscopic mechanisms mediating fiber creep is important for developing new processing protocols to make better fibers. However, creep processes usually occur at the nanoscale where direct observation is difficult. Instead, experiments characterize creep mechanisms indirectly by observing how changes in fiber length depend upon temperature and tensile load.^9–16 

In a typical creep experiment, fibers are subjected to a constant tensile load σ at a constant temperature T. The time-dependent creep rate is defined as $ε̇≡d⁡ln⁡L/dt$⁠, where L is the fiber length. Measurements are usually interpreted using an Eyring model for stress-biased, thermally activated hopping.^8,11,17 Creep is assumed to result from the activated hops of some degree of freedom over a barrier $Ea(σ)=Ea0−V*σ$ that decreases linearly with the applied tensile stress. The proportionality constant V^* is called the activation volume. V^* represents the sensitivity of the barrier to stress rather than a literal volume of hopping regions, but V^* is often comparable to the monomer volume.⁸ If the attempt frequency for activation is ω₀, Eyring theory gives a creep rate¹⁷ 

where the hyperbolic sine results from a competition between forward and backward hops that ensures the rate vanishes at zero stress.

PE’s hierarchical semicrystalline structure^6,7 produces multiple relaxation mechanisms during tensile loading.^11,14,18,19 Data are often fit to generalized Eyring models with several barriers that represent different activated relaxation processes operating in parallel. Wilding and Ward showed that creep data for a wide variety of PE fibers could be fit by assuming just two Eyring processes.^11,20 One process in their model accounts for the rapid creep that occurs as the mesoscale crystal network initially responds to the load. This “network relaxation” dominates the first 2%–4% of strain and is highly sensitive to the composition and processing history of a fiber. This suggests that it may reflect processes in the noncrystalline regions. At larger strains, $ε̇(t)$ plateaus to a steady creep rate that is controlled by a second Eyring process. Experiments suggest that this steady creep process is localized in the crystalline phase.²⁰ Unlike the first process, its contribution and activation volume are largely insensitive to processing history and only depend upon the percent crystallinity of a fiber. Of the two processes, the steady creep rate dominates the lifetime and rupture of loaded fibers because it is largest at long times. It will be the focus of this paper.

Many experiments have measured E_a and V^* for the steady, crystal creep process. Regel’ et al.⁹ and Dijkstra et al.¹⁰ studied how fiber rupture times vary with load and temperature. Later Wilding and Ward ^11,20 and Govaert et al.^12–14 tracked and analyzed the time-dependent total strain ε(t). Most recently, Jenket systematically studied how temperature and rate affected the failure of hundreds of high-performance PE fibers.¹⁶ All these studies suggest a range^10,20 V^* = 10 ų–100 ų and E_a ≈ 1 eV–1.7 eV^{12,13,16,20,21} for steady creep. The value of V^* is comparable to the volume of a single monomer in the crystal phase, which is about 47 ų.^22,23 The value of E_a was initially attributed to scission of chain backbones within the PE crystal,^9,20 but failure by scission alone is inconsistent with the experimental measurements of fiber strength.^3,21,24,25

The strongest fibers reported in experiments fail at stresses ∼6 GPa–7 GPa.^26–28 This is much lower than the theoretical strength of ∼20 GPa–40 GPa calculated for C–C bonds near room temperature.^29,30 If fibers deform plastically by chain scission, then their ultimate strengths should be substantially larger than 7 GPa. However, the past 40 years of process development have not produced substantial increases in fiber strengths. The energy barrier for C–C bond scission (∼4 eV) is also significantly larger than the ∼1.5 eV measured in creep experiments. This would lead to a dramatic difference in the rate of creep, since E_a enters in the exponential prefactor of Eq. (1).

An alternative is that chains slip past each other without breaking.^3,8,31–33 This disrupts the weaker van der Waals bonds between chains rather than the covalent bonds involved in scission. However, coherent motion of an entire chain would still be prohibitive because chains contain ∼10⁵ monomers. Instead, chains can slip by translation of a localized defect,^34–37 such as a dislocation, along the chain [Fig. 1(a)].

In a recent paper,³² we applied molecular dynamics (MD) simulations to study the effect of chain end defects on the strength of otherwise perfect, crystalline fibers. Chain ends must be present in any macroscopic fiber, and studies of this ideal case provide an upper bound for the strength of semicrystalline fibers. As illustrated in Fig. 2(a), the tension carried along the covalent backbone must be replaced by shear forces from neighboring chains within a characteristic length ξ from each chain end. Similar changes in stress and displacement occur near a dislocation between the chain and neighboring chains [Fig. 1(a)]. Our simulations found that dislocations would nucleate from chain ends above a limiting stress, allowing the chain to slip relative to its neighbors. The limiting stress at room temperature ∼6 GPa–7 GPa is comparable to the strongest macroscopic fibers and could be understood from a simple one-dimensional Frenkel–Kontorova (FK) model for the chain. These results suggest that chain slip limits the strength of PE fibers and prevents the stress from rising high enough to cause chain scission.

In this paper, we extend our studies to activated creep. Molecular dynamics simulations are used to parameterize a one-dimensional FK model for chain ends and dislocations. Four different interaction potentials are used to show that the behavior is robust and to estimate quantitative uncertainties. From the FK model, we obtain the stress-dependent activation barrier E_a for dislocation nucleation at chain ends. Once nucleated, the dislocation can move nearly freely away from the end, allowing the chain to slip by a lattice constant. The activation energy for this creep mechanism drops nearly linearly with stress, and the values of E₀ and V^* are in the experimental range. The model shows that nucleation of dislocations away from chain ends is too costly for it to contribute significantly to creep. Other possible defect mechanisms are also discussed.

Section II describes the one-dimensional FK model and its predictions for activation energies and volumes for Eyring models of creep. Section III describes the molecular dynamics interaction potentials and procedures for fitting the FK model and gives the results for four potentials. These results are compared to the experimental data for failure and creep in Sec. IV, and the possible role of other defects is discussed. Section V presents a summary and conclusions.

The conformation of chains near end defects and in dislocations is determined by the competition between tension in the chain backbone and lateral forces from interactions with neighboring chains. The FK model captures both of these effects with the minimal model illustrated in Fig. 2(b). Each C₂H₄ monomer is represented by a single bead. As discussed below, the internal degrees of freedom in the monomer have only a limited effect on chain translation. The interactions along the chain are represented by ideal harmonic springs between beads with stiffness k and equilibrium spacing b₀ = 0.254 nm determined by the equilibrium crystal structure.²² The surrounding crystal imposes a periodic potential U that is approximated by a sine wave with amplitude U₀. If the crystal is under an elastic tensile strain ε_e, the period of the potential becomes a = (1 + ε_e)b₀. The macroscopic strain of the fiber will reflect both the local elastic strain from tensile stresses and the plastic strain ε_p from chain-end slip or other effects: ε = ε_e + ε_p.

The free energy F of the FK model can be expressed in terms of u_n = r_n − na, the displacement of the position r_n of the n-th monomer relative to the minimum of the n-th potential well. Stable solutions minimize^38,39

The first term reflects the covalent energy cost of deviations from the equilibrium length b₀, and the second one gives the free energy of monomers in the corrugated potential from the surrounding crystal. The competition between the two terms introduces a characteristic length scale that corresponds to the core size of dislocations,

where τ = πU₀/a is the peak traction force from the periodic potential.

Frank and van der Merwe showed that the discrete model of Eq. (2) can be approximated by its continuum limit when ξ is larger than several lattice spacings.³⁹ This condition is satisfied for PE,³² and it allows us to use analytic solutions for many quantities. Replacing n by the undeformed coordinate x = na, stable solutions of the continuum equations satisfy³⁹ 

Equation (4) supports a commensurate solution u(x) = 0 that corresponds to all bonds stretching in registry with the crystal potential. It also supports dislocation and anti-dislocation solutions that translate the registry of the chain by one site in the crystal potential. Figure 1 shows a compressional dislocation with an extra interstitial monomer. The analytic solution of Eq. (4) for an isolated dislocation centered at the origin is given by³⁸ 

As the dislocation travels rightward, it translates monomers by one lattice spacing to the right. Bonds in the core of the dislocation are compressed from the uniformly stretched state b(x) = a. The fractional change in the length of bonds is

An anti-dislocation would produce a monomer vacancy, and u and du/dx would have the opposite sign.

In the continuum limit, the energy of an isolated dislocation relative to the uniformly stretched state is given by³⁹ 

with ε_m ≡ 2/(π²ξ). Dislocations are energetically unfavorable in an unstrained crystal, but applying a tensile strain decreases E_d for compressional dislocations because they relax tension along the backbone. For ε_e > ε_m, dislocations become energetically favorable, but they cannot form in the center of the chain without creating an anti-dislocation. Tension increases the energy of anti-dislocations, and the net cost of creating a dislocation/anti-dislocation pair is 8/πξτa. Since this barrier is constant, creep through the creation of pairs of defects would not follow the Eyring model. Moreover, we shall see that the energy barrier is quite high (120 k_BT), so double defect nucleation is extremely unlikely under experimental conditions. However, individual dislocations can nucleate from chain ends and travel into the chain.

Figure 3(a) plots u(x) for a finite chain under tension (ε_e > 0). The solution looks like part of a dislocation is entering the chain from each end. Since PE chains are much longer than ξ, the behavior at the two ends can be decoupled. In the center of the chain, u(x) = 0 and bonds are stretched by a − b₀ = ε_eb₀ to remain in registry with the stretched crystal potential from their neighbors. This stretch requires a constant tension along the chain backbone, which must vanish near the chain ends because there is a missing covalent bond. This requires that the bond length equals b₀ or du/dx → −ε_e at chain ends. As ε_e increases, more of the dislocation moves in from each end. For ε_e > ε_u = a/πξ, there is no stable solution, and the dislocation is free to move into the chain and relax stress.

For ε_e < ε_u, there is a finite energy barrier for a dislocation to enter the center of the chain. Figure 3(b) shows the transition state corresponding to a linearly unstable solution to the equations. Further motion of the dislocation to the right lowers the energy, and the dislocation will move rapidly into the chain causing the region behind it to shift by a lattice constant. For N ≫ ξ/b₀, the energy barrier between stable and unstable solutions is^38,39

where ε_* ≡ ε_e/ε_u. The dependence of E_a and E_d is shown for one model of PE in Fig. 4. The barrier for dislocation entry, E_a, decreases with increasing ε_e because the equilibrium configuration of dislocations at the ends moves further into the chain to accommodate the larger central strain. E_a vanishes once ε_e increases to ε_u. This defines the strain where the background tension is equal to the maximum compression produced by the dislocation. Chain ends cannot be stable at larger elastic strains and will spontaneously nucleate dislocations to release tension via plastic strain until stability is re-established.

An important subtlety of the above equations is that we have used the strained state a = (1 + ε_e)b₀ as the reference. Relative to their equilibrium values, denoted by subscript 0, $ξ=ξ0(1+εe)2$ and $τ=τ0(1+εe)−1$⁠. This dependence introduces additional terms that add small corrections to the energies predicted by Eqs. (8) and (7). These corrections are only of order ε_u, which we find to be only a few percent for PE. Thus, we calculate and report equilibrium values in Sec. III.

To describe elongational creep experiments, we note that the elastic strain will be proportional to the applied tensile stress: ε_e = σ/Y, where Y is Young’s modulus. In experiments on PE fibers,^6,8,20 creep is usually measured at stresses below 1 GPa and thus at elastic strains much less than 1%. In this limit, E_a and E_d are nearly the same (Fig. 4), and we can use a ≈ b₀ as discussed in the last paragraph. Then, from Eqs. (7) and (8), the activation energy for dislocations fits the Eyring model with

Nonlinear effects become important at larger stresses, and the fiber must fail catastrophically when the stress exceeds σ_max = Yε_u and chain ends slip with no activation barrier.

The above results imply that activation of dislocations at chain ends will produce plastic creep that follows Eyring theory and may explain experiments. Since the energy barrier for a dislocation to move in from chain ends drops linearly with the applied tension at small strains, the rate of dislocation nucleation will follow Eq. (1). Once nucleated, each dislocation will move along the chain in whatever direction lowers the stress. Their motion allows chains to slip relative to their neighbors and contribute to the plastic strain in proportion to their rate of nucleation. In principle, dislocation motion is hindered by a Peierls barrier from the periodicity of the chain,⁴⁰ but this is negligible for ξ/a ≫ 1, which is the same condition needed to use the continuum approximation above. The barrier for the dislocation leaving either end of the chain is small because E_d ∼ E_a for small strains.

The only parameters needed to parameterize the FK model are k, b₀, A, and τ. Section III describes how they can be determined from molecular dynamics simulations. One can make a simple estimate for V^* using the fact that covalent bonds along the chain backbone are much stronger than the van der Waals bonds between molecules. In this limit, the elastic coupling between tensile and perpendicular strains is small, and Y is approximately given by the modulus C_zz at a constant cross-sectional area. The latter is simply related to the spring stiffness k and the number of chains per unit area. This gives Y ≈ knb₀/A, where A is the area per unit cell perpendicular to the tensile axis and n is the number of chains per unit cell. Inserting this in toEq. (9) gives V^* equal to the volume per monomer within the crystal, which is roughly consistent with the experiment.

Following Ref. 32, we fit FK theory using small molecular dynamics simulations of perfect polymer crystals. This prior study showed that results from explicit MD simulations of fibers with chain-end defects were accurately reproduced by the FK models parameterized from simulations of perfect crystals. The major source of uncertainty is the appropriate choice of interaction potential. Potentials are normally fit to specific quantities over a restricted range of temperatures and pressures and may not be accurate when applied outside this range or to different quantities. To estimate the magnitude of uncertainty, we present results for four widely used interatomic potentials for hydrocarbons: AIREBO-M, AIREBO, L-OPLS, and PCFF. The Adaptive Intermolecular Reactive Empirical Bond-Order Potential (AIREBO) and its variant AIREBO-M are reactive potentials for hydrocarbons, and AIREBO-M has been shown to accurately reproduce the static and dynamic behavior of crystalline PE over a wide range of pressures.^25,32,41 The Long-Chain Optimized Potential for Liquid Systems (L-OPLS) and the Polymer Consistent Force Field (PCFF) are all-atom potentials that are commonly used for modeling polymer melt mechanics and thermodynamics.^42,43 All potentials give relatively similar results for the strong covalent backbone interactions, but they have different low energy van der Waals interactions. van der Waals interactions are harder to parameterize and are crucial in determining interchain interactions and the crystal structure. This leads to a significant variability in k and τ for different potentials.

For all potentials, we simulate 5 × 8 × 12 unit cells of a perfect and periodic PE crystal with the a, b, and c lattice vectors aligned along the x, y, and z Cartesian axes, respectively. We use the LAMMPS software package⁴⁴ to numerically solve the equations of motion with a velocity-Verlet integrator and a 0.5 fs time-step. Crystals are equilibrated to T = 300 K or T ≈ 0 K and P = 1 bar with a Nose–Hoover thermostat and barostat. The thermostat and barostat relaxation times are 0.25 ps and 1.0 ps, respectively. We measure b₀ by taking the average monomer spacing along the chain backbone at T = 300 K.

Chain stiffnesses k are calculated from uniaxial tension tests on crystals at T = 300 K. Crystals are loaded in tension along the fiber axis with a constant engineering strain rate $ε̇=108s−1$⁠, while a Nose–Hoover barostat maintains a constant 1 atm pressure in the lateral directions. Young’s modulus Y is measured by linearly fitting the axial stress σ_zz as a function of strain for the first 1% of tensile strain. Values for Y are recorded in column 3 of Table I.

While Y is most easily compared to the experiment, the stiffness k is more directly related to C_zz, the elastic modulus at a constant cross-sectional area. In Ref. 45, we found the elastic moduli for PE with the AIREBO-M potential. The value of C_zz = 270 GPa is only slightly larger than Y because the off diagonal terms from van der Waals interactions, C_xz and C_yz, are only 4 GPa–5 GPa. We, thus, take

with the understanding that this may introduce errors of order 3% that are smaller than variations between potentials.

We measure the amplitude of the intermolecular friction τ by translating one chain along the fiber axis and measuring the periodic force exerted on the chain by the surrounding crystal. Due to our small simulation sizes, collecting accurate statistics for the periodic potential at T = 300 K can be challenging. The van der Waals energies between chains are of similar order to the thermal energy scale, and this leads to significant fluctuations of the intermolecular forces about their average values. In order to make accurate measurements, we follow our previous work³² and measure the periodic potential at T ≈ 0 K where thermal fluctuations are absent. This prior study³² found that both k and τ show only a weak temperature dependence for T < 300 K, and using low temperature measurements of τ produced FK model predictions that agreed well with simulations at T = 300 K.

The z-component of the translating chain relative to the center of mass of the surrounding crystal is maintained by a stiff spring (4.34 eV/Ų). The equilibrium length of the spring is repeatedly increased to shift the axial registry of the translating chain in 0.1 Šincrements. At each increment, the average intermolecular force exerted on the chain is measured over 3 ps of simulation time, which is much longer than the time required for the crystal to accommodate the displacement.

Simulation measurements of FK parameters for the four PE models are given in the first four rows of Table I. The values for b₀ are shown in column 2 of Table I and are 2.54 Å for all four potentials, which is consistent with WAXS experiments.²² The four models show good agreement since b₀ is set by the length and angle of the C–C covalent bonds in the chain backbone, which is well understood chemistry. The backbone bond stiffness k is also dominated by covalent chemistry, and the values for the four models (column 5 of Table I) show relatively good agreement, varying from 8.7 eV/Ų to 12 eV/Ų for PCFF and AIREBO, respectively.

Values for Young’s modulus can also be compared to past experiments and theoretical predictions. WAXS measurements give values of 235 GPa–260 GPa^46–48 for crystalline regions, and lattice dynamics theories give 257 GPa.⁴⁹ The latter value is very close to the results from AIREBO and AIREBO-M. In contrast, L-OPLS and PCFF give moduli that are substantially below the experimental range.

The spring constant k depends on both Y and the monomeric volume V_mon = Ab₀/2. The fourth column in Table I shows that the models give a range of volumes from 45 ų to 50 ų. AIREBO-M was parameterized for crystal systems and gives a volume only about 1% lower than experiments. AIREBO and PCFF give volumes that are too large, and L-OPLS gives too small a volume. Note that a modified version of AIREBO that gives more accurate densities has been developed but is not part of the current LAMMPS distribution.⁵⁰ 

The fifth column of Table I gives values for k for each model and experiment using Eq. (10). Once again, AIREBO-M is closest to the experiment. AIREBO gives too large a value because of the low density, and the other models give values that are too soft by more than 10%.

The periodic force traces for the four PE models are shown in Fig. 5. The intermolecular force per monomer is plotted against the translating chain’s registry relative to equilibrium for one complete lattice period. The four PE models produce periodic potentials that are qualitatively similar to each other and to the sine potential of the FK model. All periodic potentials exhibit a small inflection at a displacement of 0.5a where the chain has a maximum energy and is linearly unstable. This inflection results from a small rotation of the chain backbone (∼5°) that allows hydrogen atoms on adjacent chains to reduce their overlap during translation. This internal degree of freedom in the monomer is not included in our FK theory but is only important near the energy maximum. The predictions for dislocation shapes and energies are mostly sensitive to displacements from equilibrium to the peak force, which is reasonably approximated by a sine wave. MD simulations of chain ends are quantitatively consistent with the FK parameterization, as demonstrated in Ref. 32.

All four curves show an increasing restoring force for small shifts in registry that reaches a maximum amplitude τ at about ±b₀/4. However, the values of τ vary significantly for the four models. PCFF and AIREBO-M have similar values of τ = 39 meV/Å and 46.5 meV/Å, respectively, whereas L-OPLS and AIREBO give values that are 50% larger. This variability is typical of different atomistic potentials since the measurements of vdW forces with either experiments or with ab initio methods are subject to larger uncertainties.⁵¹ 

The specific functional form used for vdW interactions may also contribute to the variability in τ. L-OPLS and AIREBO give the largest values of τ, and both model vdW forces with a Lennard-Jones (LJ) potential, while PCFF and AIREBO-M both use a softer exponential repulsion and give similar values for τ.^42,52 As mentioned previously, the LJ potential is known to overestimate the intermolecular interactions between chains in PE crystals and, thus, may overpredict the value of τ.^41,53 By again comparing AIREBO and AIREBO-M, which only differ in the functional form of the vdW potential, we see AIREBO gives τ = 64.5 meV/Å, which is 40% larger than the 46.5 meV/Å of AIREBO-M. Thus, the larger values of τ given by AIREBO and L-OPLS may be from overestimating steric repulsion from the LJ potential. We expect that PCFF and AIREBO-M provide more accurate measures for τ than AIREBO and L-OPLS. Note, however, that PCFF gives too small a value for k.

The measured values for b₀, k, and τ completely parameterize the FK model for each of the four potentials. The calculated predictions for the core size ξ are given in column 7 of Table I. AIREBO-M and PCFF have similar dislocation core sizes ξ ≈ 10. Even though k and τ are smaller for PCFF than for AIREBO-M, the core size and critical strains are sensitive to the ratio kb₀/τ that is similar for the two models. Due to their much larger interchain friction, AIREBO and L-OPLS give smaller core sizes equal to 8.7 and 6.3 monomers, respectively.

When comparing to experiments, it is important to remember that our model considers an ideal and perfectly crystalline fiber that cannot be achieved in the laboratory. Real fibers have additional plastic relaxation processes associated with disorder and other types of defects. These mechanisms can produce additional creep that our model does not account for. We have also assumed that slip events are independent of each other, while the plasticity of real fibers may be correlated. For instance, when a chain slips, it transfers the tensile load it was carrying to neighboring chains, which may increase their tendency to slip and disorder the lattice.³³ In addition, disorder in real fibers leads to variations in the local stress that activates creep, reducing the effective activation energy. Thus, we expect our theoretical predictions to set an upper bound on fiber strength and a lower bound on creep observed in experiments.

The FK theory predicts ε_u ≈ b₀/(πξ₀) is the maximum strain a crystal can support before yielding. Above this strain, the crystal is unstable and chains can spontaneously slip at their chain ends. This sets an upper bound on the high-rate/low-temperature strength of a PE fiber. Assumming Young’s moduli do not change much over this strain interval, the maximum load a fiber can sustain is

The predicted ε_u and σ_u for the four FK models are recorded in column 2 of Table II. The highest experimental values from the literature for ε_u^16,21 and σ_u^26–28 are also provided.

The four models predict ultimate strengths ranging from 6.02 GPa for PCFF to 10.60 GPa for L-OPLS with an average value of σ_u = 8.5 GPa. Note that all values are much lower than the failure stress of 20 GPa–40 GPa that is predicted to be necessary for C–C scission.^29,30 The highest published strengths are 6 GPa–7 GPa and tend to be found for fibers with the highest crystal concentration. This is consistent with the average value of σ_u from the four models. Only the PCFF model gives σ_u in the achieved experimental range, suggesting that it underpredicts τ and k.

FK model predictions for the maximum strain ε_u also agree well with the recent experimental analyses by Jenket and by Wang et al.^16,21 Jenket studied the tensile failure of hundreds of PE fibers while systematically varying the temperature and strain rate. For all rates, he observed a limiting strain at failure of ε ≈ 0.03 as fiber temperatures were decreased.¹⁶ Through extrapolation of their own data, Wang and Smith predicted a maximum strain at failure of 0.034 for a perfectly crystalline fiber.¹⁶ The value of ε_u from L-OPLS is about 50% higher than this value, but all other models give values within 10% of this extrapolated value.

Note that the strength of a fiber can exceed σ_u if the chain ends are removed from the fiber interior. This may be impractical for conventional processing of high-performance fibers, but it was recently demonstrated at the nanoscale by Shrestha et al.⁵⁴ They used thermal annealing to create a highly ordered, nanoscale PE filament that is much smaller in length and diameter than the extended length of the chains. This nanofiber is expected to contain few chain ends, and Shrestha et al. reported failure strengths approaching 12 GPa, which are the largest ever reported for PE fibers.

Creep is controlled by the energy barriers for activated motion. Figure 4 shows the barrier for a dislocation to enter from a chain end, E_a, as a function of elastic strain for the FK model parameterized by AIREBO-M simulations. The barrier vanishes at ε_u, leading to the catastrophic failure described in Sec. IV A. Creep experiments are typically performed at much lower elastic strains. The largest loads applied in the creep experiments by Ward et al. were only about 0.5 GPa,^11,20 corresponding to ε_e ∼ 0.002. In this regime, dislocations should remain dilute, and their activation energy decreases nearly linearly with stress or elastic strain. Thus, their contribution to creep can be described by Eyring theory. Note that high rate experiments see failure at total strains that may be outside this linear regime (ε > 0.02).⁵⁵ 

In Fig. 6, we directly compare model predictions for the stress dependence of E_a to the experimental data. We plot Eq. (8) using values for $Ea0$ and V^* averaged over all four models (solid line). The blue region gives a rough measure of the uncertainty in $Ea0$ derived from the standard deviation of the four models. Experimental data of Jenket,¹⁶ Govaert et al.¹² and Regel’ et al.⁹ are shown as points. For the data of Jenket and Govaert, only an activation energy at 0 load is given. Broken lines show Eyring fits to the experimental data as reported by Ward and Wilding (dashed lines)²⁰ and Dijkstra (dashed-dotted line).¹⁰ Ward and Wilding fit many fiber types by first finding an optimal $Ea0≈1.3$ eV and then varying V^* for each fiber. The span of their fits is indicated by the gray shaded region. Except for the data of Dijkstra et al.,¹⁰ the average model prediction sits above experimental values for all tested loads. We expect this to be the case since our perfect fibers at zero temperature should set an upper bound on $Ea0$ for creep in actual fibers. Indeed, even in crystals, there are variations in the distance between chain ends and other defects that will produce lower barriers for some chain ends, and fibers have non-uniform local stresses that will facilitate local nucleation.³² 

The FK model predictions for $Ea0$ are recorded in column 4 of Table II. Models give values ranging from 1.2 eV to 1.9 eV with an average of ≈1.5 eV. Model predictions are consistent with the experimental range of 1.1 eV–1.7 eV and well below the energy for chain scission of C–C bonds ∼4 eV. AIREBO-M gives $Ea0=1.4$ eV ≈ 60 k_BT, which is closest to the average for both models and experiments.

Model predictions for V^* are also consistent with the range of experimental values. The four FK models give similar values for V^* between 45 ų and 50 ų that are recorded in column 5 of Table II. As discussed previously, FK theory predicts that V^* for chain slip is equal to the monomer volume within the crystal, which x-ray scattering experiments find to be about 45 ų.²² Few experiments have measured V^* during creep, but available data suggest that it lies between 10 ų and 100 ų.^9,20 While the experimental variation in V^* is high, the range is compatible with the FK model values for chain slip.

PE hosts a variety of crystal defects including braids, hair-pins, crank-shafts, and twistons.^{8,34–36,56,57} Twistons, in particular, have been considered as a source for creep. They are soliton defects that smoothly rotate a single chain backbone by 180° while also stretching or compressing the chain backbone by half a monomer spacing.^34–36 These are called twist-extension (TE) and twist-compression (TC) twistons, respectively. Twistons have been observed in equilibrium crystals with nuclear magnetic resonance (NMR) and are known to play an important role in the disordering of the crystal phase during melting.^36,37,58

In order to explain the success of Eyring theory, twistons would have to produce plastic flow via an activated process whose activation energy had the right magnitude and a linear decrease with stress. TC and TE twistons must form in pairs within a chain.⁵⁸ TE twistons are extended defects such as the dislocations considered here, and the calculated activation energies for PE^34–36 range from ∼0.7 eV at 0 K to ∼0.43 eV at T = 350 K in unstrained crystals. Thus, the creation of a pair of defects has too low an energy to explain creep. Moreover, tension would increase the energy of TE twistons, while Eyring theory requires a mechanism whose energy decreases. The energy of TC twistons is lowered by tensile stress, but they collapse into a localized sequence of highly rotated dihedral angles.^34,35 In contrast to dislocations, Mowry and Rutledge showed that these collapsed TC defects are relatively immobile with large energy barriers (0.14 eV–0.65 eV) for translation along the chain.⁵⁶ TC twistons could in principle nucleate from chain ends like dislocations, but our own repeated attempts to force nucleation of such defects at chain ends in MD simulations failed–only producing localized gauche configurations of the terminal bonds–and they were never observed spontaneously.

A simple analytic model of thermally activated chain slip was used to predict the steady creep behavior of PE fibers. Our theory is informed by prior work where MD simulations identified chain slip as the strength limiting failure mechanism of PE fibers. Chain slip is mediated by the nucleation of 1D dislocations at chain ends, which can be modeled with a simple Frenkel–Kontorova model. The FK model has three independent parameters that are parameterized by MD simulations of perfect PE crystals.

We parameterized FK models for four widely used atomic potentials: L-OPLS, PCFF, AIREBO, and AIREBO-M. Of the four, AIREBO-M best represented the average response and produced an equilibrium crystal structure and elastic moduli in best agreement with the available literature.^23,45–48 Averaged over the four models, the FK theory predicts an average activation barrier $Ea0≈1.5$ eV and activation volume V^* ≈ 47 ų for chain slip in a perfect crystal. These effective Eyring parameters should set an upper bound on the energy barrier for creep in real fibers. Additional defects and disorder within a real fiber should tend to decrease $Ea0$ and increase V^*. The effective Eyring model for chain slip gives a stress-dependent activation barrier E_a(σ) that remains near but above all available experimental data, except for the data of Dijkstra et al.¹⁰ The average $Ea0$ from FK theory lies above all considered experimental data except for that of Dijkstra et al.¹⁰ 

The FK theory for chain slip is relatively simple, but it produces an activation energy barrier that is consistent with steady creep experiments and is unambiguously distinct from the barriers for twiston formation or C–C bond scission.^34–36,56 The chain slip mechanism also sets upper bounds on the ultimate strain and strength of fibers that is consistent with experiments.^16,26–28 We believe these results make a strong case for 1D dislocations being the microscopic mechanism mediating steady creep within PE crystals. At present, we know of no other microscopic model that gives comparable predictions for fiber creep.

Our theory only considers the nucleation of dislocations at chain ends in an otherwise undefected fiber. Real fibers contain amorphous domains and crystal–amorphous interfaces. These disordered interfaces may serve as additional locations for chain slip to occur and could lead to lower experimental values of E_a and V^*.^18,19 Thus, our Eyring parameters should be seen as setting an upper bound on thermally activated creep in PE fibers.

Even as an upper bound, the $Ea0$ predicted for chain slip is significantly lower than the 4 eV binding energy needed for C–C bond scission. This suggests chain slip will be active and dominate crystal creep well before chain scission. While we do not expect chain scission to control the creep rate, EPR experiments have shown that scission does occur during tensile loading. In order to be consistent with our predictions, these chains should not be breaking within the bulk crystal phase. Instead, we expect chain scission to be occurring primarily on “tie-chains” bridging separate crystallites within the fiber, which may be highly stretched due to the low stiffness amorphous regions. Recent AFM experiments characterized the internal structure of high-performance PE fibers and observed numerous tie-chains between adjacent crystal domains.⁷ These can be pulled taut and break as adjacent crystallites shear past each other during deformation.^59,60

T.C.O. dedicates this last collaboration to his cherished mentor and friend, Professor Mark O. Robbins, who passed away during the publication of this manuscript.

This material is based on the work supported by the National Science Foundation (Grant No. CMMI-1628974) and by the Army Research Laboratory under the MEDE Collaborative Research Alliance (Grant No. W911NF-12-2-0022). Simulations were performed at the Maryland Advanced Research Computing Center. This work was performed, in part, at the Center for Integrated Nanotechnologies, an Office of Science User Facility operated for the U.S. Department of Energy (DOE) Office of Science. Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. DOE’s National Nuclear Security Administration under Contract No. DE-NA-0003525. The views expressed in the article do not necessarily represent the views of the U.S. DOE or the United States Government.

The data that support the findings of this study are available from the corresponding author upon reasonable request.