Shear induced martensitic transformations in crystalline polyethylene: direct MD simulation

TL;DR

This study uses molecular dynamics simulations to explore shear-induced martensitic phase transitions in crystalline polyethylene, revealing how shear rate influences the transition pathways and mechanisms at the atomic level.

Contribution

First MD simulation of shear-induced martensitic transformations in polyethylene, demonstrating different transition modes depending on shear rate and deformation speed.

Findings

Rapid shear induces direct transition to triclinic phase.

Slow shear leads to intermediate monoclinic phase before triclinic.

Chain rotations involve short twist defects diffusing along chains.

Abstract

For the first time, we carry out molecular dynamics (MD) simulation of shear-induced martensitic phase transitions between the orthorhombic and non-orthorhombic (triclinic and monoclinic) phases of crystalline polyethylene (PE) in the framework of a realistic all atom model of the polymer. We show that the variation of the shear rate allows observing on a nano-sample both a strongly nonequilibrium phase transition occurring by random nucleation and irregular growth of a new phase ('civilian' way, for rapid deformations) and the coherent, or 'military', kinetics (generally considered as usual for martensitic transformations). We induce transitions from the orthorhombic to the triclinic phase according to two transformation modes observed in experiment on PE single crystals. Rapid deformation favors the transition directly to the triclinic phase, slow deformation - first to the intermediate monoclinic, and only then - to the triclinic phase. The second way corresponds to the experiment on extended chain PE. We explain this result and analyze the competition between different transformation and plastic deformation modes. Rotations of PE chains around their axes necessary for the transition between the orthorhombic and non-orthorhombic phases are executed by short twist defects diffusing along the chains. The transition between the monoclinic and triclinic phases occurs through half-chain-period translations of the chains along their axes, mostly collectively, as crystallographic slips.

Tables

Table 1. Table 1: Potentials of interactions between the atoms (C - carbon, H - hydrogen) in the used all atom PE model - the AMBER parm99 force field Cheatham  et al. ( 1999 ) .

Valence bond potential: $U\left(L\right)=K_L{\left(L-L_0\right)}^2$
C-C: L0=1.526Å, KL=310 $\mathrm{kcal}{\mathrm{mol}}^{-1}$Å-2
C-H: L0=1.090Å, KL=340 $\mathrm{kcal}{\mathrm{mol}}^{-1}$Å-2
Valence angle potential: $U\left(\theta \right)=K_{\theta }{\left(\theta -{\theta }_0\right)}^2$
C-C-C: ${\theta }_0={109.5}^o$, $K_{\theta }=40$ $\mathrm{kcal}{\mathrm{mol}}^{-1}{\mathrm{rad}}^{-2}$
H-C-H: ${\theta }_0={109.5}^o$, $K_{\theta }=35$ $\mathrm{kcal}{\mathrm{mol}}^{-1}{\mathrm{rad}}^{-2}$
H-C-C: ${\theta }_0={109.5}^o$, $K_{\theta }=50$ $\mathrm{kcal}{\mathrm{mol}}^{-1}{\mathrm{rad}}^{-2}$
Torsion angle potential:
$U\left(\phi \right)={\sum }_{n=1}^3{K}_{\phi n}\left(1+\cos \left(n\phi +{\phi }_{0,n}\right)\right)$
C-C-C-C:	$K_{\phi 1}=0.2$ $\mathrm{kcal}/\mathrm{mol},{\phi }_{0,1}={180}^0$
	$K_{\phi 2}=0.25$ $\mathrm{kcal}/\mathrm{mol},{\phi }_{0,2}={180}^0$
	$K_{\phi 3}=0.18$ $\mathrm{kcal}/\mathrm{mol},{\phi }_{0,3}=0^0$
C-C-C-H:	$K_{\phi 1}=K_{\phi 2}=0$
	$K_{\phi 3}=0.04$ $\mathrm{kcal}/\mathrm{mol},{\phi }_{0,3}=0^0$
H-C-C-H:	$K_{\phi 1}=K_{\phi 2}=0$
	$K_{\phi 3}=0.0375$ $\mathrm{kcal}/\mathrm{mol},{\phi }_{0,3}=0^0$
van der Waals pair interactions between atoms separated by more than 3 bonds or belonging to different molecules:
$U\left(r\right)=U_{LJ}\left(r\right)S\left(r\right)$
$U_{LJ}\left(r\right)=\epsilon \left({\left(R_{min}/r\right)}^{12}-2{\left(R_{min}/r\right)}^6\right)$
${\epsilon }_{CC}=0.1094$ $\mathrm{kcal}/\mathrm{mol}$, $R_{min,CC}=3.816$Å
${\epsilon }_{HH}=0.0157$ $\mathrm{kcal}/\mathrm{mol}$, $R_{min,HH}=2.974$Å
${\epsilon }_{CH}={\left({\epsilon }_{\mathrm{CC}}\cdot {\epsilon }_{HH}\right)}^{1/2}$;
$R_{min,CH}=\left(R_{min,CC}+R_{min,HH}\right)/2$
$r_{in}=12$Å, $r_{out}=13.5$Å
For the atom pairs belonging to the same molecule and separated by 3 bonds, $S\left(r\right)$ was multiplied by 0.5

Table 2. Table 2: Results of MD simulations for the large samples, all but one consisting of (short) chains of 50 carbon atoms. The one of 100 carbon atoms is marked as ’long’. For all the velocities of the lid, we listed the observed modes, the resulting phase, the maximal (or minimal for negative stresses) value of S xy component of the stress tensor, the values of engineering shear strain at the beginning (g 1 ) and at the end (g 2 ) of the transition and the duration of the transition.

shear for:	vel.,	observed modes	result	Sxy,	g1,	g2,	time,
	$m/s$	in order of appearance		MPa	$deg.$	$deg.$	$ns$
O$\to$T,
$T{2}_1$	1	O-T $T{2}_1$(0.35); O-M $T{3}_1$(0.25); T-M $T{2}_1$(0.215): all civilian	T+M(long)	185	-3.61	-12.75	1.42
			M	162	-3.07	-8.24	0.744
	1/4	O-M $T{1}_1$(0.04), civil.; O-T $T{2}_1$(0.35), military; M-T $T{3}_1$(0.31), milit.	T	150	-2.23	-16.48	9.18
	1/16	O-M $T{1}_1$(0.04), civil.; twinning of M, military; M-T $T{1}_1$(0.17), civil.	T	193	-4.29	-20.17	31.2
O$\to$T,
$T{1}_1$	1	O-M $T{1}_1$(0.04); O-T $T{1}_1$(0.18); M-T $T{1}_2$(0.17): all civilian	T	-247	2.85	9.69	1.26
	1/4	O-M $T{1}_1$(0.04), civilian; M-T $T{1}_2$(0.17), military	T	-212	2.45	9.11	4.84
T$\to$O,
$T{2}_1$	1/4	O-M $T{1}_1$(0.04), civil.; T-O $T{2}_1$(0.35), military; T-M $T{4}_1$(0.38), milit.	O+M	-63	0	13.1	8.5
with nucl.

Full text

I Introduction

Deformation induced martensitic transformations can be accompanied by twinning and generation of dislocations in the parent phase. This is the case in steel and other alloys, such as CoNi. Then the defects can prevent the reverse transition. In shape memory alloys the energy barrier for the phase transition is, on the contrary, much lower than for the plastic deformation modes. Under load, they are not actuated both in parent and product phases, and the transition is almost perfectly reversible. In some materials, the heights of these two barriers are comparable. In this case, the initiation of the transition can be impeded, and the behavior of the system under load can depend on many factors in a complicated manner. An example of such a transition is the transition from the orthorhombic (O) to the triclinic (T) phase in polyethylene (PE) (see Fig. 1).

When solidified from a melt under ordinary conditions, polyethylene (PE) forms a semicrystalline sample (about 70% crystallinity), in which tiny O (high symmetry and high temperature phase) crystallites (lamellas) consisting of folded chains are embedded by ’tie’ molecules into the amorphous matrix Sperling (2006). The ground state of the system (an extended chain crystal) does not form because of kinetic reasons: long polymer chains have no time for it. In the lamellas, a load applied to the semicrystalline sample causes (in order of frequency) sliding on some crystallographic planes, twinning, and transitions of O crystallites into non-O phases.

The non-O phases are always present in commercial PE. A large amount of them was found in PE crystals formed on substrates Wittmann and Lotz (1989); Furuheim et al. (2004), after polymerization inside nanochannels Nair et al. (2005) and in reactor powders Aulov et al. (2004); Marikhin et al. (2018). These phases also appear under pressure, tension or shear, but the complete or almost complete transitions from the O to the non-O phases were observed only on specially prepared samples: on single crystals grown from a solution Kiho et al. (1964); Allan and Bevis (1974a), on extended chain (poly)crystals Miyaji et al. (1987); Fontana et al. (2007), or on (poly)crystals consisting of stacks of chain-folded lamellas Seto et al. (1968); Young and Bowden (1974). In all these cases, the chains in all the crystallites were parallel.

High-symmetry O phase is preferable at high temperatures. If one deforms an O single crystal at 110∘C, the transition to the T phase does not occur Kiho et al. (1965a). Annealing for 5 minutes at 100∘C of a semicrystalline sample containing T crystallites leads to their complete transition back to the O phase Seto et al. (1968). The single crystals that have transformed into the T phase, return to the O phase after removing the applied stress if they relax freely Kiho et al. (1965b).

On the other hand, when a sample made of chain-folded lamellas is released after compression, only a part (which appeared at the late stage of deformation) of the formed T phase returns to the O phase, another part just twins Seto et al. (1968). In extended chain crystals, after a compression-decompression cycle even at temperatures more than 2000C, the coexistence of the O and non-O phases is observed down to 2.3 GPa (under compression, non-O phases appear only at 6 GPa) Fontana et al. (2010). The very high temperature used by Fontana et al Fontana et al. (2010) seems to be needed for the full recovery of the O phase under conditions of the experiment. From these facts, one has to conclude that there should be high barriers dividing the free energy minima corresponding to the O and non-O phases.

In addition, the picture is complicated by the existence of two non-O phases having higher free energies than the O phase: monoclinic (M) and T (Fig. 1). They both appear at compression of extended chain crystals at 280${}^{0}C$ (M at 6 GPa, T at 15.5 GPa), they coexist at high pressures (40 GPa), and they both are present at decompression, the share of the T phase increasing with pressure decrease Fontana et al. (2007, 2010). It is possible that the M phase also appeared in the other experimental works cited above, but it was not detected because of similarity between the lattice parameters of the M and T phases.

The densities of the O, T and M phases are very close in the wide range of pressures Fontana et al. (2007), therefore the transitions between these phases are considered to be martensitic. More exactly, the deformation needed for transitions between the phases is an invariant plane strain, almost a simple shear. Based on the observed diffraction patterns, Seto Seto et al. (1968) determined the approximate relative orientation of the product T and parent O phases. The crystallographic analysis of the phases shows that there is a whole series of possible transformation modes bringing one lattice to another Bevis and Crellin (1971). The smallest displacements of molecules in $(xy)$ plane are provided by modes $T1_{1}$, $T1_{2}$ (relative locations of the molecules in these two modes are very close to each other and consistent with the conclusion of Seto’s analysis Seto et al. (1968), but the invariant planes are different) and $T2_{1}$, $T2_{2}$ (see Fig. 2).

In tensile deformations of monocrystals, modes $T1_{1}$ and $T2_{1}$ are operative Allan et al. (1973); Allan and Bevis (1974b), while only mode $T1_{2}$ was observed in compression experiment on single crystal textured (bulk) PE Young and Bowden (1974). The further experimental investigation of the transitions between O and T phases in strained PE monocrystals was hindered by two unexpected problems. First, the mode of transformation proved to be dependent on fold sector (the pattern of returning of folded chains into the crystal) Allan and Bevis (1974b), i.e. on boundary conditions. Secondly, as the martensitic transitions are caused by shear, the tension of a crystal induces too many competing modes of deformation, such as twinning and crystallographic slip Allan and Bevis (1975). The resulting picture was too complicated for productive theoretical analysis. Both the obstacles are also present for compression experiments on samples consisting of stacks of lamellae. They rotate and twin both in parent and product phases Seto et al. (1968); Young and Bowden (1974). A successful shear experiment on any proper sample has never been reported. The reverse T$\rightarrow$O transition (induced by annealing) has been purposely studied only in one work Takahashi et al. (1988).

A molecular dynamics (MD) simulation of the transition can become a guide for the halted experimental studies. But, although polymorphism is the most prominent feature of polymer crystals, its theoretical investigation is hindered by the lack of adequate realistic models reproducing the phase transitions under study. Besides, the adequate statement of numerical experiments is not evident as regards boundary conditions, shear directions, and sizes of samples.

Luckily, for PE, the simplest polymer, there is a phenomenological all atom force field, AMBER parm99 Cheatham et al. (1999), which reproduces all the three phases of PE, their energies being in the required order (T is the ground state, O possesses the highest energy at 0K, the difference in energies and densities is very small, the barriers to transitions are rather high). All the phases are (meta)stable at normal conditions. Thus, it is possible to study the transitions between the phases in crystalline PE by direct molecular dynamics simulations, and to observe the kinetics of the transitions at the level of individual molecules.

In this work, we present a series of numerical experiments allowing to answer the most intriguing question on the O-T transition in PE: why is there the M phase between the parent O and the product T phases in the O-T transition Fontana et al. (2007)? What prevents the direct transformation? Is there a way to actuate it?

II Details of MD simulations

We wanted to simulate ’ideal’ martensitic transitions between the lattices of crystalline PE. Namely, the transition should be induced by shear, not by tension or compression (as any PE lattice may be transformed to another by (almost) simple shear), and from the edges of samples. The samples should be small enough so that the stress field during the deformation at finite temperature was as homogeneous as possible, but the samples should be large enough so that the influence of the edges was minimal. We also wanted to avoid the dependence of the transition mode on the fold geometry of the polymer crystal Allan and Bevis (1974a): chain folds can prevent some modes from being operative.

Therefore, the samples consisted of finite chains (with their axes oriented along $z$ direction): 48 CH2-groups (plus two terminal CH3-groups) or 98 CH2-groups (plus two terminal CH3-groups). In one direction ($y$), the sample was restricted by two planes (orthogonal to the $y$ axis) approximately parallel to the invariant plane (shear plane K1 Bevis and Crellin (1971)) of the transformation mode under study. The shear direction in this plane $\eta_{1}$ coincided with the $x$ axis. To stabilize the samples under shear, we imposed periodic boundary conditions along the $x$ axis. One of the samples is shown in Fig. 3.

The molecules belonging to the upper and the bottom edges of the sample formed two ’lids’ by which the deformation was carried out. The molecules of both the lids were rigid zigzags with axes parallel to the $z$ axis. The molecules could shift along $z$ axis and rotate around it, but not around $x$ or $y$ axes. Such lids kept the sample intact.

In the process of shear, the axes of the molecules of one of the lids were shifted along the $x$ axis with a constant speed of 1, 1/4, or 1/16 m/s. They also could move along the $y$ axis independently of each other. The $x$ and $y$ coordinates of the mass centers of the molecules of the second lid were fixed, except for the action of a barostat which could change the $x$ coordinates of the mass centers of both the lids. Thus, the deformation we carried out was slightly different from a simple shear, but at this cost the sample had sufficient space to accomplish the phase transformations which requires, as we have mentioned in Introduction, not exactly the simple shear but an invariant plane strain.

To perform the MD simulations, we used the LAMMPS Plimpton (1995) package LAM (a) with the AMBER parm99 potential set Cheatham et al. (1999) (see the force constants in Table 1). The time step was equal to 0.5 fs. The samples were kept at 300K with the use of the Langevin thermostat (the damping parameter of 1ps) applied to non-rigid chains. The Nose-Hoover barostat LAM (b); Shinoda et al. (2004) with the same damping parameter of 1ps kept zero stress along the x axis (Sxx component of the stress tensor). The visualization was made with VMD VMD support.

III Direct O to T transition

For direct O$\rightarrow$T transitions, only three orientation relations were detected in experiment. They correspond to transformation modes Bevis and Crellin (1971) $T1_{1}$, $T2_{1}$, and $T1_{2}$ (see Fig.2). We intended to implement modes $T1_{1}$ and $T2_{1}$ known to be operative in single crystals of PE Allan et al. (1973); Allan and Bevis (1974b). The results of the MD simulations are summarized in table 2.

We compared the transitions in large and small (in cross section) samples consisting of short and long chains. For the large samples, the rigid molecules in the lids constitute 10% of the sample, and for the small samples - 20%. One may expect that for the small samples the picture of the phase transition will be distorted (which turned out to be not the case), while in the large samples the stress field and the deformation will be inhomogeneous (which proved to be correct in some cases). We stopped the simulation when a plastic deformation of the formed composite samples started. The residual stress almost instantly relaxed upon release of the sample, without causing the reverse transition. We shifted the lid at four speeds - 1, 1/4, 1/8, and 1/16 m/s (the first speed corresponds to shear rate of 108s*-1* for the large samples). The transitions took place within 1-32 nanoseconds.

III.1 Transformation mode O-T T2

Mode $T1$ requires the smallest displacements of the molecules, while the transition according to mode O-T $T2$ should noticeably bend the crystal (see Fig. 2), and, indeed, the vertical rows in the samples visibly curve during the shear. We start with mode T2, as it shows more variety in behavior depending on the deformation rate and the size of the samples.

III.1.1 O-T T2, rapid deformation (1 m/s): random nucleation and irregular growth of the T phase, strong competitor: O-M T3

Transformation begins with a turn of one chain through 90∘ (see Fig. 3). The chain becomes a nucleus of the product T phase which grows with the deformation. As expected, the phase is oriented relative to the parent O phase according to transformation mode $T2_{1}$. During deformation, several nuclei appear and grow rather irregularly (in ’civilian’ way Christian (2002)): in different directions, but when they are small - preferentially along the $x$ axis, and after formation of a horizontal layer - in the $y$ direction. Rotations of the chains of the O phase happen more often near the T nuclei. Here arise M domains oriented according to mode O-M $T3_{1}$ (see Fig. 2). In the small samples, it happens rather late, and the M phase does not noticeably spread. The result of deformation of the samples consisting of long chains and the samples consisting of short chains is more or less the same (Fig. 3). Namely, the T phase (mode O-T $T2_{1}$) dominates. At the boundary of this phase, there are a domain of the M phase (O-M $T3_{1}$) and a remnant of the parent O phase. In the two small samples, the location of the non-T phases differed. It was near the upper lid for the sample with short chains (50 carbons), and near the bottom lid - for the sample with long chains (in both the cases, we shifted the bottom lid to the right). In the large samples, the domains of the M phase (mode O-M $T3_{1}$) grow especially easily in the areas between the formed T nuclei (Fig. 4(1)).

In the large sample with long chains, the phases T and M ($T3_{1}$) spread out at similar rates, and, at the end of the deformation, share the sample approximately equally (see file O-T_T2_1_(100CH2).mp4 in Supplemental Material mov ), and this ratio does not change if one continues the shear. In the sample with short chains the M phase ($T3_{1}$) emerges simultaneously with the T ($T2_{1}$) phase, propagates much sooner, and, in addition, easier consumes small T domains (by mode T-M $T2_{1}$). At the end of the deformation, the M phase ($T3_{1}$) occupies the whole sample (see file O-T_T2_1.mp4 in Supplemental Material mov ).

III.1.2 O-T T2, ’regular’ deformation (1/4 m/s): ’military’ propagation of the T phase, weak competitor: O-M T1

If we decrease the rate of the lid in the large sample with short chains to 1/4 m/s, the result of the deformation and the evolution of the sample are quite different (file O-T_T2_0.25.mp4 in Supplemental Material mov ). First, in the whole sample, only one nuclei of the T phase (mode $T2_{1}$) emerges and grows. Before its formation, a small domain of the M phase (mode $T1_{1}$) appears. Later, it is being gradually occupied by the expanding T phase. Secondly, the growth of the T phase is not random, but very close to what can be described as a ’military’ way Christian (2002). Formed in the center of the sample, the horizontal layer of the T phase expands to the lids, its boundaries being almost parallel to the invariant plane, and the sample deforming correspondingly (Fig. 4(2)). The result of the ’regular’ (0.25 m/s) deformation is the expected pure T phase (mode $T2_{1}$), contrary to the rapid (1 m/s) shear of the same sample that has led to pure M phase (mode $T3_{1}$), which does not even appear at ’regular’ rate.

III.1.3 O-T T2, slow deformation (1/16 m/s): two intermediate M phases, twinning and crystallographic slips

If the shear rate is further reduced by a factor of 4 (down to 1/16 m/s), the deformation results, as at velocity of 1/4 m/s, in a pure T phase (mode $T2_{1}$), but the transition flows differently (file O-T_T2_0.0625.mp4 in Supplemental Material mov ). The sample first transforms into the M phase oriented according to mode O-M $T1$. We saw one row of this phase at the velocity of 1/4 m/s. This phase grows not through a front propagation. The nuclei of the M phase arise mainly at the lids of the sample and represent one molecule turned through 90∘. The growth of the M phase occurs only along the vertical row in which the nucleus appears. Molecules flip one after another, starting from the lid. In the following deformation, the M phase changes its orientation to mode O-M $T2_{1}$ (relative to the parent O phase). The transformation goes not through gradual distortion of the cell geometry and gradual change of setting angles of the molecules. Instead, we saw the formation and development (in the ’military’ way) of a twin of the first intermediate (O-M $T1$) M phase. The twin quickly caught all the sample. The following deformation of the second (O-M $T2_{1}$) M phase went through crystallographic slips along directions in ($xy$) plane until, in the end, this phase transformed to the T phase according to mode M-T $T1_{1}$ (see Fig. 2) by crystallographic slips parallel to the chain axis.

III.2 Transformation mode O-T T1

III.2.1 O-T T1, rapid deformation (1 m/s): intermediate M phase, random nucleation and irregular growth of the T phase

When the lids are parallel to the invariant plane of transformation mode O-T $T1_{1}$, the intermediate M phase (mode O-M $T1_{1}$) appears during the transition even at the high strain rate. The nucleation and growth of this M phase are similar to that of the intermediate M phase emerging in the slow deformation (1/16 m/s) according to mode O-T $T2_{1}$. After a noticeable spread of the domains of the M phase, it begins to transform to the T phase through the slip of half the chains along their axes. T domains emerge at random places and grow independently of each other, in ’civilian’ way (see video file O-T_T1_1.mp4 in Supplemental Material mov , till 0.89ns), and their lattices may not fit together. But when two ’incoherent’ domains come into contact, one of them reorganizes (through the transition to the M phase) so that they form a whole (0.94ns). Thus, during the transition, the domains of all the three phases, O, M, and T, coexist, but it does not cause a noticeable deformation of the sample (Fig. 4(3)) because the lattice parameters for the properly chosen cells are very close at this orientation relation (according to mode $T1_{1}$). The result of the transition is an ideal T phase. The kinetics of the transition in the samples composed of long and short chains is approximately the same. The only feature of the transition in the sample with short chains is that in it the nuclei of not M, but immediately T phase also occasionally emerge.

III.2.2 O-T T1, ’regular’ deformation (1/4 m/s): intermediate M phase, almost ’military’ propagation of the T phase

At a lower strain rate (1/4 m/s) in the large sample with short chains, the transition begins in the same way as at the high rate, with the growth of the M phase nuclei (according to the same O-M $T1_{1}$ mode) starting from the edges. However, the T phase begins to appear much later, when most of the sample is already in the M phase. The transition to the T phase starts mainly from one of the edges, the M phase and the remnants of the O phase being captured rather in ’military’ than in ’civilian’ way. Although the boundary of the T phase is noticeably different from a plane, but, as in the case of ’regular’ deformation according to mode O-T $T2_{1}$, this boundary is marked by a slight but still visible curvature of the vertical rows in the sample. Thus, the slowing down of the shear leads, as for mode O-T $T2_{1}$, to more cooperative kinetics of the transition (file O-T_T1_0.25.mp4 in Supplemental Material mov ) and to the transition first to the intermediate (O-M $T1_{1}$) M phase, and only then to the T phase.

IV Reverse T to O transition

To study the reverse T to O transition, we stopped the deformation in O-T $T2_{1}$ ’regular’ (1/4 m/s) transition (the large sample with short chains) at the moment close to the finish and changed the velocity of the moving lid to the opposite (see file T-O_T2_0.25.mp4 in Supplemental Material mov ). As expected, the boundaries of the T domain went back to the lids. However, in the remainder of O phase, there appeared the common ’parasitic’ M phase (O-M $T1_{1}$) that spread together with the O phase catching the T domain (there is the proper mode T-M $T_{4}$, see Fig. 4(4)). And we remember from the direct transitions that this direction of shear favors the transition of the O phase to this M phase on their boundaries parallel to the $y$ axis. Therefore, the sample after the shear was a mosaic of the O and M phases.

V Discussion.

V.1 Mode O-M T1: intermediate and parasitic.

After the described simulations, we can answer the question posed at the end of the Introduction about the reason for the appearance of the M phase between the O and T phases in the compressed sample Fontana et al. (2007). Geometrical parameters of the M lattice are much closer to those of the O lattice than the geometrical parameters of the T lattice. Namely (see Fig. 2), the M phase can be obtained from the O phase by a simple shear with a very small magnitude of 0.04 (mode O-M $T1$). An analogous (with the close orientation relation, $T1$) mode for the O-T transition has a magnitude of 0.18. Because the corresponding shear directions are fairly close, it is not unexpected that we have seen that the shear in the direction of mode O-T $T1$ causes a transition to the M phase according to mode O-M $T1$ at any shear rate.

Moreover, mode O-M $T1$ may, due to its tiny magnitude, lead to the growth of the M phase also under shear in any direction. Indeed, there are two modes of transformation O-M $T1$: $T1_{1}$ and $T1_{2}$, their invariant planes being almost perpendicular (see Fig. 2). A simple shear in any direction along any invariant plane will have a positive projection on the shear direction of one of the modes $T1_{1}$ or $T1_{2}$. Therefore, under quasistatic shear in any direction in the O lattice, mode O-M $T1$ should be actuated. And the shear directions for two modes O-T with the smallest amplitudes $T1$(0.18) and $T2$(0.35) are very close to the shear direction for mode O-M $T1$.

Only in the rapid shear in the direction far from the direction of O-M $T1$, the M phase did not have time to be born before the deformation becomes large enough to actuate another mode with larger magnitude, but in a direction closer to the shear. In our simulations, we did not observe the mode O-M $T1$ only in the transition O-T $T2$ at a speed of 1 m/s (then the modes O-T $T2_{1}$(0.35) and O-M $T3_{1}$(0.25) appeared, see Fig. 4(1)). At both lower speeds 1/4 and 1/16 m/s this mode was present. It was even actuated as a ’parasitic’ in the remnants (nuclei) of the O phase in the reverse T-O transition (Fig. 4(4)). This mode can be operative also in rapid shear, if the shear direction is closer to the needed direction: in the transition O-T $T1$ at the speed of 1 m/s, the mode O-M $T1$ was the first to appear.

It is unlikely that in hydrostatic compression experiment by Fontana et alFontana et al. (2007), local shear rates could be higher than in our simulations (see the estimates in section V.2). Therefore, it is not surprising that at the beginning of the experiment the M phase was observed. It is equally natural that after the transition to the M phase, the crystallites under the action of the same stresses passed into the T phase: the direction of shear for the transition O-M $T1$ is close to the direction of M-T $T1$ (see Fig. 4(3)). Thus, it can be said that the path from the O to the T phase naturally lies through the M phase.

V.2 Rate of shear and kinetics of transitions.

At first glance, the used velocities of the molecules of the lid of our nano-sample –1/16-1 m/s - and the corresponding shear rates – ($0.6\cdot 10^{7}$-$10^{8}$)s*-1* - seem unrealistically large. Even under impact deformation, shear rate reaches only 105s*-1* at striker velocities of about tens of m/s. On the other hand, the speed of sound along PE chains (in our model) is 12, 2.6 or 1.9 km/s (depending on polarization), and in transverse directions – 2-4 km/s. The velocities of the front propagation were about 0.55 m/s for mode $T2_{1}$ and 2.2 m/s for mode $T1_{1}$, which is three orders of magnitude lower than the speed of sound in this direction. The experimental estimate of the rate of growth of the martensite phase in iron under rapid heating is about 1 km/s (sound velocity is about 6 km/s) Nishiyama (1978).

There are several competing deformation modes in the O lattice of crystalline PE. Shifts in crystallographic planes (especially parallel to $z$ axis), twinning, and the O-M $T1$ mode (0.04) have minimal magnitudes. The case of the transition O-T $T2$ (see Table 2) can illustrate the possible dependence of the transition kinetics and even of the product phase on the strain rate when several transformation modes can be actuated. One can assume that (if the boundary conditions allow it), under quasistatic shear, the first mode to be actuated will be the one for which the projection of this shear on its direction will first reach the needed magnitude.

At the lowest speed (1/16 m/s) of shear in the direction of mode O-T $T2$(0.35), we observed sequentially: a transition to the M phase according to the O-M $T1$ mode, twinning of the M phase, and a transition to the T phase according to M-T mode (mainly through shifts in crystallographic planes parallel to the $z$ axis). However, the end result was the expected T lattice with orientation relation corresponding to the mode O-T $T2$. When we increased the speed up to 1/4 m/s, the transition to the M phase according to the O-M $T1$ mode happened only in one vertical row, and, accordingly, the entire subsequent chain of transitions did not follow. Instead, the expected O-T $T2$ mode was actuated, and the transition went according to it (in military way).

With a further increase in speed up to 1 m/s, in addition to the mode O-T $T2$(0.35), the mode O-M $T3$(0.25) is also actuated, predominantly at the boundaries of the T domains. The T and M phases develop simultaneously. In this case (as opposed to all the other simulations), the result of the transition depends on the length of the chains in the sample. In the sample of short chains (50 carbon atoms), the mode T-M $T2$(0.215) is also operative, and the transition results in the M phase. In the sample of long chains (100 carbon atoms), the T and M phases share the sample. Thus, the result of this highly nonequilibrium transition is critically dependent on boundary conditions.

The shear rate (together with boundary conditions) also determines what kind of kinetics of the transition to expect: military or civil. For the transition O-T $T2$(0.35), a speed of 1/4 m/s was ’regular’: it allowed to carry out this transition in standard military way without actuating modes with smaller magnitudes. The higher speed of 1 m/s led to the civil kinetics of the nucleation and growth of the T phase and the parasitic mode O-M $T3$. A lower speed of 1/16 m/s already gave time for the nucleation of the intermediate M phase according to the (low-amplitude) mode O-M $T1$.

For the transition O-T $T1$(0.18), which is closer to mode O-M $T1$ both in orientation relation and in magnitude, this intermediate M phase appears even at a speed of 1 m/s. The speed providing the military transition to the T phase has to be higher. The speed at which the direct nucleation of this T phase would go in civilian way is still higher. For this mode, we did not seek to find both the speeds because they are not currently available in experiment. The decrease in speed for the O-T $T1$ transition leads to the change in the kinetics of the second transition M-T $T1$(0.17) from the intermediate M to the final T phase: civilian kinetics at a speed of 1 m/s becomes almost military at a speed of 1/4 m/s.

The observed civilian kinetics is an example of a strongly non-equilibrium first order phase transition. The nanoscale domains of two or three phases coexist, separated by lengthy boundaries of complex shape (Fig. 4(1) and (3)). Formation of these boundaries requires more energy, which is provided by greater stresses (see Table. 2), accompanying the rapid transitions. The excess energy can be spent on formation of the boundaries between many nuclei of the new phase and the matrix.

V.3 Boundaries between the phases: twist defects and tension defects.

To accomplish a transition between the O and non-O phases, in addition to the displacements of molecules in the plane $(xy)$, their rotations around their axes are required. In transformation modes O-T and O-M $T1$ and $T2$, half of the chains in the sample have to turn through approximately 900. Both the transitions start with the generation of a T or M defect in the O lattice: a chain rotated counter-clockwise or clockwise through 90 degrees. The energies of these defects proved to be approximately equal to 0.2 kcal/mol per one CH2-group. Both the defects correspond to local minima.

Using the nudged elastic band (NEB) method JÓNSSON et al. (1998), we found that there is much lower barrier on the transition path to these minima through generation of a twiston (see Fig. 5) or, under periodic boundary conditions, a pair of twistons, than through the rotation of the chain as a whole (see Fig. 6, on the left).

Similar twist defects were observed in MD simulation of premelting in a PE crystal Zubova et al. (2012).

The generation of one twiston turning a chain of O lattice to the position as in T defect costs 4 kcal/mol, and as in M defect - 2 kcal/mol. The barrier between the O lattice and the T defect is 20 kcal/mol for a finite chain with 100 carbons (Fig. 6, on the right). But in a deformed sample, the energy of the T defect lowers with the shear. For example, under 31% of the shear required for transition O-T T11, the energy of the defect is the same as the energy in the original position in O lattice. To turn a chain under this or greater shear, one needs to create only the twiston (see Fig. 6, on the right). Therefore, in the described phase transitions, the upper boundary estimate for activation energy under shear equals to the energy of the twistons, 4 kcal/mol (for the direct O-T transition) or 2 kcal/mol (for the O-M transition).

The twistons emerge at the ends and sometimes (in pairs) in the middle of a chain. On long chains, we observed three twistons at the same time: a couple near the middle, and one near an end. The average overall width of a kink is 8-14 CH2 groups (the core of the defect - about 6), and the average speed - 250-650 m/s (which allows to turn a chain containing 100 carbon atoms in 20-50 picoseconds).

Twistons are often generated at close time moments on neighboring chains. It should be noted that the speeds of the sound waves propagating along PE chains and polarized parallel and perpendicular to the axes of the chains are 12 and 2.6 (or 1.9) km/s, respectively. We observed deceleration of the kinks and their movement in the opposite direction. In this case, the turn of a chain could not be completed, and the chain remained in the old phase. So, the motion of the twistons along the chains is a diffusion process, but it is a very quick diffusion.

Analogously to rotation around the axis, a long polymer chain in a PE crystal does not shift along the axis as a whole. A shift by a full chain period may be effected by motion of a vacancy (localized defect of tension) along the chain from one of its ends to another Zubova et al. (1999). The chain diffusion between amorphous and crystalline fractions in semicrystalline samples is implemented by motion of cheaper twist-tension defects (shift by half of a chain period and twist through 1800 brings a CH2 group into the position of the nearest neighbor) Zubova et al. (2007).

In the process of transition between the intermediate M and the product T phases, we observed shifts of chains by half-chain-period. Namely, in the transition from the M to the T phase according to mode M-T $T1$(0.17) half of the chains in M lattice shifts by half-chain-period (which is equivalent to a turn through 1800). These shifts are implemented by defects of tension. But the rigidity of a PE chain in tension (along its axis) is much more than the rigidity in torsion (around the axis). Therefore, the defects of tension are much wider than twist defects. For example, for the twist-tension defect responsible for the chain diffusion, the width of the area of tensile deformation is about 80 CH2 groups, while the area of twist deformation - 25 CH2 groups Zubova et al. (2007). Therefore, when a chain is 50 CH2 groups long, its shift looks as if the chain moves as a whole. The chain shifts very quickly, within a few picoseconds. We have often observed simultaneous displacements of several chains situated in one crystallographic plane. In this case, their cooperative shifts corresponded to crystallographic slips.

VI Conclusion

In the simplest polymer, PE, there are three crystalline phases (meta)stable under normal conditions, and able to transform one to another by shear-induced (first-order) phase transitions. There are many transformation modes ((almost) simple shears in different directions with different magnitudes) realizing these phase transitions. The height of the barriers for actuation of these modes is comparable to the height of the barriers for plastic deformation modes (crystallographic slips or twinning) Young and Bowden (1974). Thus, when a sample is deformed, the modes of deformation and of the phase transitions compete, and it is known that the result of the deformation depends not only on its type, but also on the structure of the sample. Another feature of phase transitions in polymer crystals is the micro-dimensions of single crystals or crystallites in semicrystalline samples and, accordingly, the strong influence of boundary conditions on these transitions.

An experimental study of the transition between the O and T phases of PE is hindered because of the need for special preparation of samples in which the competition with the deformation modes would be minimal. But even for such samples, the influence of the boundary conditions makes the productive theoretical analysis of the experiments very difficult. It is also not known whether there are another experimental conditions affecting the process of the transition. Many experimental results are unclear, in particular, the appearance of an intermediate M phase between the parent O and the product T phases in an extended chain sample under hydrostatic compressionFontana et al. (2007).

An MD study could serve as a guide for the halted experimental studies. Although for most polymers there are no force fields adequately modeling phase transitions between crystalline phases, the AMBER’s force field Cheatham et al. (1999) reproduces all crystalline phases of PE with correct relative energies.

In the present work, for the first time, in the framework of this realistic force field, we carried out an MD investigation of the phase transition between the O and T crystalline phases of polyethylene. We studied the first two (with minimal magnitudes) transformation modes which were experimentally observed. In our MD experiment, the transition was initiated by a shear along the invariant plane for the mode under study. The periodic boundary conditions were imposed only in one direction, and the chains had a finite length.

The study showed that the course and (to a lesser extent) the result of the transition depends on the rate of the shear deformation. A shear in the direction of mode O-T $T2$ (observed in single crystals) at high strain rate (velocity of the lid is 1 m/s) actuated not only the expected mode but also the transformation mode to the M phase O-M $T3$, and the transition went in civilian way: through random nucleation and irregular growth of the new phases. The result was (depending on the length of the chains in the sample) a combination of the T and M phases, or the M phase. At a lower speed (1/4 m/s), the same shear caused the coherent (military) transition to the expected T phase. The slow deformation (1/16 m/s) sequentially actuated the modes with smaller magnitudes and shear directions close to the shear direction of the deformation: O-M $T1$ (result: the M phase), twinning of this M phase (result: the M phase, oriented already according to mode O-M $T2$, similarly to the expected mode O-T $T2$), and finally, M-T $T1$ (result: the T phase, oriented according the expected mode O-T $T2$). So, the result of the slow deformation was the expected T phase oriented according to the mode we wanted to actuate, but the transition went through two differently oriented intermediate M phases.

The transformation mode to the first M phase O-M $T1$ has a very small magnitude of 0.04 (O-T $T1$ and $T2$ modes have magnitudes of 0.18 and 0.35, respectively). Because there are two nearly perpendicular directions of shear causing this mode, for any sufficiently slow simple shear, this mode will be actuated. For the shear in the direction of mode O-T $T2$, the critical speed was 1/4 m/s (one row of the M phase had already emerged, and was then absorbed by the T phase). Under the shear in the direction of mode O-T $T1$ (observed both in single crystals of PE and in bulk PE), the O-M $T1$ mode is being actuated even at a speed of 1 m/s. The resulting M phase transforms into the T phase by the same mode (M-T $T1$) as in the slow transition according to mode O-T $T2$. Thus, we see that the M phase (mode O-M $T1$) naturally emerges as an intermediate in the transition from the O to the T phase according to both modes O-T $T1$ and $T2$ observed in experiment.

The MD simulation also allows studying the transition process at the level of individual molecules, unavailable in experiment. In particular, we could observe the process of rotation of half of the chains through 90 degrees around their axes during the transition from the O phase to T or M. The turn proved to occur generally due to generation of a localized twist defect (twiston) on one end of a chain and the subsequent diffusion of the twiston to another end of the chain. The twiston is short (about 6 CH2 groups), moves at speed of 250-650 m/s and rotates a chain containing 100 carbon atoms in 20-50 picoseconds. The energy of its generation is 4 kcal/mol for the O-T transition and 2 kcal/mol for the O-M transition. In this connection, one can notice that the close heights of the barriers for the studied phase transitions and crystallographic slips in polymers with chain diffusion between amorphous and crystalline fractions is not an accident. The reason is that the molecular mechanism of dislocation generation and motion also includes generation and motion of localized twist-tension defects Nikolov and Raabe (2006); Spieckermann et al. (2017).

In transition from the M to the T phase (or vice versa) according to mode M-T (or T-M) $T1$, the crystal needs not only a (small) deformation of its cell, but also a shift of half of its chains along their axes by half a period (or the rotation of the chains around their axes through 180 degrees). In our samples, consisting of finite chains 50 or 100 CH2-groups long, we observed mainly collective half-chain-period translations corresponding to slip in a crystallographic plane. The shift of one chain occurs within a few picoseconds.

In the present work, only the first step is done in MD investigation of the phase transitions in PE. In particular, we did not study the reasons for the experimentally observed difficulty of the reverse transition to the O phase under normal conditions. We also did not analyze the initiation of the shear transformation modes, actuated along with the phase transitions in commonly used in experiment deformation processes: hydrostatic or uniaxial compression, or stretching. We also did not study the effect of boundary conditions (the length of chains and the direction of chain loops in single crystals and lamellae) on the actuation of a particular mode of phase transitions or deformations. We hope that our work will stimulate the research in these directions.

Acknowledgements.

The work was supported by the Russian Science Foundation (award 16-13-10302). The simulations were carried out in the Joint Supercomputer Center of Russian Academy of Sciences.