Interlude of metastability in the melting of polymer crystals

Consider a crystal made of low molar mass simple molecules. Upon heating this crystal to a temperature above the equilibrium melting temperature, the crystal melts spontaneously. This spontaneous melting behavior is not common for semicrystalline polymers. The nonspontaneity of melting of polymers is widely recognized to arise from free energy barriers for melting.^1–3 It is reasonable to expect free energy barriers for melting of polymers due to conformational correlations associated with the loss of crystalline order for temperatures below the equilibrium melting temperature. However, for temperatures higher than the equilibrium melting temperature, the origin of nonspontaneity of melting remains to be fully understood. In general, the melting of semicrystalline polymers is a process that has been investigated far less than its opposite process of crystallization. Along with its relevance to technological applications, understanding the process of melting of a tightly bound crystal is key to uncovering several long-standing mysteries in the field of polymer crystallization, such as the memory effect of the melt temperature on the crystallization shown experimentally.^4–19 The melting of a homopolymer crystal represents the unraveling of the simplest macromolecule, which fundamentally represents the competition between thermal forces and the conformational entropy of polymers. Understanding how such crystals melt could also significantly impact our knowledge of several diseases such as Alzheimer’s, which are commonly known to be caused due to misfolding of proteins,^20,21 as they are the complicated counterparts of homopolymers.

The process of melting was historically studied using inferences drawn from X-ray scattering,^22,23 which seemed to indicate that melting occurred from the surfaces of crystalline lamellae in both melt-crystallized as well as solution crystallized polymer crystals. More recently, with the availability of several more advanced techniques, attention has shifted to the processes that occur during the melting itself. Barham and Sadler²⁴ first showed through neutron scattering that the melting of polyethylene single crystal mats results in an almost “explosive” change to a random coil structure, with the radius of gyration changing by 5 nm/s. The time taken for this change in size was observed to be on the order of seconds and independent of molecular weight. In a later publication,²⁵ Barham showed that this melting time is uncorrelated to the relaxation time of individual molecules.

With the development of differential scanning calorimetry (DSC), several experimentalists then focused on the energetics and kinetics of melting.²⁶ The use of sinusoidal temperature modulation²⁷ in several experiments indicates that partial melting is reversible up to a certain degree.^28–30 This indicated a preference for the melting to occur, along with keeping a nucleus intact onto which the partially molten parts could recrystallize back, without needing to overcome a free energy barrier. Further, Hu and co-workers³¹ extended this idea to show that this reversibility in melting is dependent on the ability of the polymer to carry out a sliding diffusion. Later, it was believed that there could be an activation barrier in the melting kinetics of polymer crystals,^1,32 which was attributed to the entropic barrier due to pinning.² 

Several authors have also tried to observe the melting process by annealing different polymers using atomic force microscopy. An overwhelming majority of works indicate that melting begins from the edge of the single polymer crystal.³³ Some other authors report stepwise unfolding³⁴ and crystal thickening when observing the melting of oligomers³⁵ and ultrahigh molecular weight polyethylene.³⁶ Several experiments reported that the melting of single polymers caused several cavities to appear.^37–41 Later, Hobbs showed that the parts of the crystal which were crystallized at lower temperatures melt before the material in confined geometries, followed by the rest of the lamellae from the edges inward.⁴² This is the reverse of the sequence during crystallization, where thick crystals form first before thin crystals. This was observed separately by Dubreuil and co-workers as well.⁴³ Sequential melting in different sectors of polyethylene single crystals has been observed as well.^44,45

While most of the experiments have focused on single crystals, the study of highly concentrated systems has also yielded several insights into the mechanism of melting. In their work, Rastogi et al.⁴⁶ showed that there is a long-lived heterogeneous state that is formed when ultrahigh molecular weight polyethylene is melted, which they argued was formed because of restricted reptation due to a heterogeneous distribution of entanglements. In a subsequent work,⁴⁷ they showed that there are two possible mechanisms of melting depending on the presence or absence of entanglements, corresponding to different bundle sizes of chains that detach from the crystal. The analysis, however, was complicated due to an annealing process below the true melting point of 141 °C. Later, Lippits and co-workers⁴⁸ showed evidence for important processes occurring before the melting point is reached. Rastogi and co-workers⁴⁹ then followed the melting of chains using solid state NMR experiments and argued that adjacently re-entrant domains melt at lower temperatures and clusters of chains melt at higher temperatures, resulting in two types of structures depending on whether melting happens fast or slow. Yao et al.⁵⁰ then showed that the entropic differences between the crystalline and amorphous parts of semicrystalline polyethylene could decide the diffusion of chains during melting. Pandey and co-workers⁵¹ then observed the melting temperatures (T_m) of entangled and unentangled polyethylene and clearly showed increasing melting temperature with heating rate. This was attributed to a lag arising out of instrumental measurement constraints. However, as mentioned earlier, Toda et al.^1,2 found an increase in melting temperature with heating rate. This observation was after accounting for instrumental deficiencies. This indicates the possibility of activation barriers being present during melting.

Recently, with the development of fast-scan Differential Scanning Calorimetry (DSC),⁵² crystals have been melted at very high heating rates (∼10⁵ K/s). This meant that the crystals which were formed at lower temperatures and were susceptible to reorganization at lower melting temperatures no longer reorganize before the other parts. The artificial increase in the melting point due to this reorganization can be avoided, and the true nature of melting can be elucidated.⁵³ In one such work, Toda and co-workers suggested the existence of activation barriers during melting.⁵⁴ 

Several of these experimental findings have been accompanied by theoretical and simulation calculations as well. Using kinetic Monte Carlo simulation of melting of a single molecule crystal in a simple cubic lattice, Hu and co-workers^55,56 argued that the free energy barrier for melting is the reverse of that for crystallization. In their simple model, only parallel polymer bonds are subjected to an attractive interaction (with effective energy E_p) and all the other significant features such as the trans-gauche conformational states and polymer-polymer, polymer-solvent, and solvent-solvent interactions are ignored. Based on this model, they showed that a first-order phase transition occurred at the melting temperature of 3.125 E_p/k_B, where k_B is the Boltzmann constant. They then addressed the issue of free energy barriers for three temperatures below this phase transition temperature of 3.125, namely, 2.941, 2.967, and 3.030 in units of E_p/k_B. The barrier at the temperature of 3.030 is about 25 k_BT. Such high barriers would result in enormously long times for melting, which is unseen in experiments. The key conclusion from this simulation work is that there is a high barrier for melting. Most importantly, the focus of the work has been at temperatures below the equilibrium melting temperature.⁵⁶ The analysis at the different temperatures has been complicated because of discrepancies in the determination of the equilibrium melting temperature between their actual simulation data and the post analysis of their data in constructing the free energy profile. Furthermore, such large barriers are estimated using the total system without explicit considerations of an order parameter, such as the global order parameter of Liu et al.⁵⁷ or Welch et al.⁵⁸ This has hidden the nature of melting much above the melting temperature of the polymer. Furthermore, in the above work and a subsequent publication,⁵⁶ Hu and co-workers treat the free energy landscapes of melting and crystallization to be the reverse of one another, with a single barrier in the free energy landscape. On the other hand, the pathway of crystallization into the final semicrystalline state of global free energy minimum consists of multiple entropic barriers,^59,60 and the reverse process is likely to follow a different set of entropic barriers. Additionally, it is important to note that the free energy barriers of ∼25k_BT are at temperatures below the melting point of the polymers. Moreover, a recent theory from the recrystallization of once-molten crystals invokes the presence of long-lived metastable states during the melting of a polymer³ and explains several experimental findings^7–11,18,61 using such states.

The goal of the present paper is to gain insight into the molecular basis of the nonspontaneity of polymer melting even at temperatures above the equilibrium melting temperature. To understand this, we study the melting process by use of Langevin dynamics simulations. Using our molecular model,⁶² which is derived after making modifications to another model that essentially reproduces several thermophysical properties of polyethylene,⁶³ we estimate the free energy landscapes using parallel tempering Langevin dynamics simulations along a new local order parameter, which serves as a reaction coordinate. We study two separate systems, single chain and multichain crystal, and show the kinetics of melting in both systems. We show that both single crystals and multichain crystals pass through a globular metastable state at lower melting temperatures, before forming more expanded random coil-like structures at higher melting temperatures, which for the latter system results in individual chains peeling off the crystalline core. For the multichain crystals, we also study melting at temperatures close to the equilibrium melting temperature and show the existence of long-lived metastable states of partial order. We estimate free energy landscapes using parallel tempering simulations and show the existence of a metastable state during the melting of both systems. For single chain melting, we show one free energy barrier, while for multichain melting, we show two free energy barriers. Our work highlights the melting behavior of polymers in an ideal case, where all the monomers are subjected to the same temperature upon melting, within statistical fluctuations of the Langevin thermostat. The paper is organized as follows: in Sec. II, we describe the simulation model and the sampling procedure; in Sec. III, we discuss our findings from the simulations.

Our simulation model⁶² incorporates just enough detail to obtain chain-folded structures without being very computationally intensive. This model is designed to mimic a polyethylene chain, with each methylene(CH₂) group taken to be a bead. We therefore have long chains of a degree of polymerization upto 700. All the interaction parameters have been modeled after the work of Paul, Yoon, and Smith,⁶³ which has successfully reproduced several thermophysical properties of polyethylene. Although they treat the terminal methyl (CH₃) group separately, we treat all the beads of the chain to be the same. We make another small modification to the model, which we discuss later.

Our model incorporates force fields to represent chain connectivity, chain rigidity, and nonbonded interactions. Therefore, the total potential energy arises due to contributions from the chemical bonds, bond angles, dihedral angles, and the van der Waals-like short ranged interaction. The bond angle θ is the angle formed between the bond vectors of any two successive bonds $r→12$ and $r→23$⁠. The dihedral angle ϕ is the angle between the two planes, which are in turn formed by three successive bond vectors $r→12$⁠, $r→23$⁠, and $r→34$⁠. The potential energy associated with the chemical bonds is taken to be harmonic

where r is the bond length and r₀ is the equilibrium bond length. In this potential, k is taken to be 115 kcal/mol Ų and r₀ is taken to be 1.54 Å. The potential energy associated with the bond angle is taken to be of cosine-squared form

where θ₀ is 109° and k_θ is 60 kcal/mol.

We also model chain stiffness by invoking a potential energy arising due to the dihedral angle ϕ and this is of the multiharmonic type given by

where k₁ is 3.02 kcal/mol, k₂ is −0.560 kcal/mol, and k₃ is 2.58 kcal/mol.

We model the nonbonded interactions by means of a modified Lennard-Jones potential, and this is given by

where the interaction strength ϵ is set to 0.112 kcal/mol. The equilibrium distance σ is 4.53 Å for beads further than five repeat units apart along the chain backbone. Beads closer than five repeat units interact with a reduced value of σ equal to 1.54 Å. This is expected to have minute effects on the behavior of the chain other than slightly increasing the chain’s local flexibility. The form of the Lennard-Jones is slightly different from the usual form because here the Lennard-Jones has its minimum value at σ rather than the expected $216σ$⁠.

For the simulation, we use reduced units, and all the results that we present are in these units. All quantities are derived from the basic units of mass m (130 g/mol), equilibrium bond length l₀ (1.54 Å), and Lennard-Jones energy ϵ (0.112 kcal/mol). Hence, the united atom quantities of mass, equilibrium bond length, and Lennard-Jones interaction energy are all taken as 1. The reduced time t^* is then given in units of $ml02ϵ$ (2.56 ps), and the reduced temperature T^* is given in units of $ϵkBNA$ (56.38 K), where k_B is the Boltzmann constant and N_A is the Avogadro number.

We have calculated an approximate value of the persistence length l_p according to the relation⁶⁴ 

where l_i represents the ith bond vector, l represents the equilibrium bond distance, and l_p represents the persistence length. By taking an average over 1000 snapshots of a linear polymer consisting of 700 monomers at a reduced temperature of T^* = 12, we find the value of the persistence length to be 7.52 (11.58 Å), after all the modifications to the model.

The interaction potentials are combined and integrated⁶⁵ with respect to time by using the methodology of Langevin dynamics.^66,67 In Langevin dynamics, the motion of the particles is described by the Langevin equation of motion, wherein the particles are subjected to the Langevin thermostat along with the interaction potentials. This is collectively given by the equation

Here, U_i is the net potential obtained by adding all the aforementioned potentials together. The $ζṙi$ term represents the frictional drag, and ζ is the friction coefficient, while Γ represents the random noise.

The Langevin dynamics method simulates the solvent by assuming that the solvent particles collide with the solute particles, but their negligible size in comparison with the solute molecules means that the entire collision can be parametrized in terms of a Gaussian noise Γ. The Gaussian noise is related to the friction coefficient by means of the fluctuation-dissipation theorem

Additionally, this Gaussian noise has the property of being time-averaged to 0,

We set the friction coefficient ζ to be 1 and integrate the Langevin equation using the velocity Verlet finite difference algorithm.^63,67 In this algorithm, the velocities of the particles are calculated at every half time step, before computing the new positions using those half time step velocities, leading to greater accuracy. The time step size that we use in these simulations is 0.001.

In an effort to monitor the ordering of the polymer on a local level, we introduce a new local order parameter in the spirit of Hu et al.⁵⁵ and Nicholson et al.⁶⁸ We start from the global order parameter,^57,58,60,69

where θ is the angle between all pairs of orientational vectors $r→i+1−r→i−1$ for any bead i. Then, we modify this global order parameter to define it for each monomer by identifying all the monomers which fall within a cut-off radius. The cut-off radius is set at 6.5 (10.01 Å), in keeping with the distance at which the Lennard-Jones potential U_LJ is truncated. Once the monomers are identified, the angles between the orientation vectors of these monomers are then averaged according to Eq. (9) to give a local value of the order parameter. If this local order parameter exceeds a value of 0.5, then the monomer is designated as a crystalline monomer and amorphous otherwise. The threshold value has been chosen as 0.5 after making inferences from the structural order parameter values obtained first by Liu and Muthukumar.⁵⁷ For the multichain crystals, the monomers from all the chains are utilized to generate the orientational vectors, as long as the monomers fall within the cut-off radius. In this manner, we develop a count of the number of amorphous monomers m_amorphous for every time step and use this as our order parameter for melting.

To calculate the free energy landscapes for the melting of single chains and multichain crystals, we use the method of parallel tempering along with Langevin dynamics^70,71 as done earlier by Mahalik and Muthukumar.⁷² In the parallel tempering method, there are a certain number of replicas of the system. We run Langevin dynamics with each of these replicas at their prescribed temperatures. Periodically, the information on the configurational details in a replica is exchanged into another replica chosen at random, with an acceptance probability, and the system would evolve in the new replica at the new temperature. After an elapse of the swapping time period, the swapping of configurational information is performed into another replica, and this process is continued. We have shown the swapping time period for all simulations in the supplementary material. If a state α is chosen to swap with a state β, then the acceptance probability P_α→β for this swap is given by

where U_α is the total potential energy of the α^th replica and T_α is the temperature of the α^th replica. To ensure that the average kinetic energy in any replica stays at $32kBT$⁠, we re-scale the momentum p of each monomer i of the system following a configuration swap as given by

In order to make sure that the energies of replicas overlap sufficiently for Eq. (10) to be used, we use an evenly spaced temperature grid containing the temperatures T_c = 9, T_m − 0.4, T_m − 0.2, T_m, T_m + 0.2, and T_m + 0.4, where the melting process has been studied for different melting temperatures T_m of 12, 14, 16, 18, and 20. Here, T_c and T_m are crystallization and melting temperatures, respectively. Additionally, we have also studied the melting of multichain crystals at temperatures of 11 and 11.6 to understand the kinetics close to the equilibrium melting temperature. Mahalik and Muthukumar⁷² have suggested other methods to improve convergence. But as the results of the parallel tempering method will show in Sec. III, convergence was attained using the chosen temperature grid. Therefore, no other temperature grids were utilized. Several crystalline configurations were maintained at T_c = 9 and were also part of the sampling scheme so that the system could attempt swaps with crystalline configurations as well. The total number of steps for which the simulation was run for each melting temperature and the frequencies of the swap step are provided in Table S1 in the supplementary material. The relative free energy contribution of each conformation to the free energy landscape along the melting reaction coordinate m_amorphous is given by

where $nmamorphous$ represents the number of conformations found that contain m_amorphous number of amorphous monomers and $N$ represents the total number of conformations sampled.

We first studied the melting of a single polymer crystal consisting of 700 monomers. The crystals were obtained from the crystallization process of linear polymers of N = 700 as described in our previous work.⁶⁹ Then, these crystals were melted isothermally by allowing them to equilibrate at their melting temperatures T_m for several thousands of simulation time units through the use of Langevin dynamics.

The melting process was studied at several different values of T_m. The equilibrium melting temperature $Tm0$ was determined in our previous work as 10.74 ± 0.20.⁶⁹ To study melting, we chose five different values of T_m—12, 14, 16, 18, and 20. As a reference for our model polymer, the difference between the temperatures of 12 and 10.74 in real units is 71.04 K. The choice of temperatures was made so as to explore the response of the polymer crystals to increasing driving force for melting.

Figure 1 (Multimedia view) shows typical trajectories of melting at two different temperatures, rendered using the Visual Molecular Dynamics (VMD) package.⁷³ Figure 1(a) (Multimedia view) shows a trajectory when the polymer crystal was allowed to equilibrate at T^* = 12, which represents a case of weak driving force, while Fig. 1(b) (Multimedia view) shows a trajectory when the polymer crystal was allowed to equilibrate at T^* = 18, which represents a case of strong driving force. At low melting temperatures like the one shown in Fig. 1(a) (Multimedia view), the polymer is found to progressively and uniformly lose its crystalline order throughout the crystal. Once the polymer loses its order, it equilibrates to a dense globulelike structure, until several thousands of Langevin dynamics time units, at which point we end the simulations. However, at higher melting temperatures, it is found that the polymer is able to escape the excluded volume attractions (Lennard-Jones attractions) and expand to a coil-like structure, after passing through the dense globulelike state as an intermediate. In essence, the melting process comprises of a step where the chain loses order and the other step where the chain expands. The accessibility of the latter step depends on the temperature.

In order to show the change in size of the polymer with an increase in the melting time, we show the ensemble-averaged radius of gyration with respect to simulation time. We have averaged the radius of gyration over 70 simulations. Figure 2 shows the radii of gyration at different values of T_m. The figure shows that at the temperatures of 12 and 14, the polymer does not expand beyond a certain radius of gyration. In fact, the radius of gyration becomes smaller at the melting temperature of T^* = 12 as the polymer changes from an aligned lamella to a globular state. However, at higher melting temperatures, the polymer chain undergoes a transition similar to that discovered by Barham and Sadler²⁴ using neutron scattering. This is seen in Fig. 2, where we can observe a sudden transition in the radius of gyration within a few hundred time units of the beginning of the simulation at the melting temperatures of 16, 18, and 20. It should be noted that Barham and Sadler studied different molecular weights of polyethylene lamellae at 145 °C and found the melting kinetics to be independent of molecular weight. An exact mapping of the melting behavior to that of Barham and Sadler²⁴ is beyond the scope of this work. Nevertheless, our simulation results are in qualitative agreement with the experimental results.

In addition to the change in shape of the polymers upon melting, we also describe the evolution of the crystalline order of the polymers. We describe the crystalline order of the polymer using the global average of the second order Legendre polynomial, as described by Eq. (9). Figure 3 shows the evolution of the ensemble-averaged value of P₂ at different T_m. The polymer is seen to lose order rapidly, and an associated time scale is evident from the figure. We assign the time at which P₂ first passes a value of 0.1 as the melting time. This choice has been made after the global order parameter values first provided by Liu and Muthukumar.⁵⁷ Based on this melting time, we decide the swapping frequency for the parallel tempering with Langevin dynamics simulations to describe the free energy barrier of this melting step. Table I in the supplementary material describes the melting time for each of these melting temperatures. Our results, however, are independent of this choice of 0.1 for the order parameter. It can be observed from the figure that the rate with which the crystal loses order increases with temperature. We then proceed to study the kinetics of the melting process.

To quantify the free energy landscape of the melting process, we have performed parallel tempering Langevin dynamics simulations. Parallel tempering ensures that the sampling is not biased toward any particular region in the landscape and that detailed balance is maintained.^70,71

To perform parallel tempering, we run 6 different Langevin dynamics simulations at 6 different temperatures, as described in Sec. II. This ensures that the sampling has access to all the sections of m_amorphous. The number of simulations, the sampling frequency, and the other details of the parallel tempering are shown in Table I of the supplementary material.

Using the parallel tempering in such a manner, we have then computed the free energy landscape along the number of amorphous monomers. A value of 0 amorphous monomers corresponds to a completely crystalline system, while a value of 700 corresponds to completely amorphous monomers. As examples, we show structures corresponding to m_amorphous = 83 and m_amorphous = 698 in Fig. 4(a). The crystalline monomers are marked in blue, while the amorphous monomers are marked in red. m_amorphous can then be thought of as a reaction coordinate. We have computed a free energy landscape using Eq. (12), by discretizing the reaction coordinate using a bin size of 50 for the number of amorphous monomers, and using a histogram method, like the one used by Welch and Muthukumar.⁵⁸ Our free energy estimates are shown in Fig. 4.

In Fig. 4(b), we first discuss the crystalline region of the landscape. Our initial crystalline structures contain m_amorphous = 83, as shown in Fig. 4(a). Thus, the contribution to the landscape below the m_amorphous = 100 region comes from purely crystalline structures. Even though the free energy values for m_amorphous = 50 appear to have a very high value in comparison with that of the intermediate regions of the landscape, they correspond to highly crystalline configurations.

Next, we discuss the intermediate regions of the landscapes. The free energy landscapes reveal clear barriers in going from the crystalline to the molten state, as seen from the free energy values when going from m_amorphous = 350 to m_amorphous = 400. Although on a relative energetic scale, these barriers seem to be nearly similar, we have also proceeded to compute the exact heights of these barriers F^* and have given those in Table I. These free energy barriers are of the order of 1.5 k_BT. The heights of these barriers are considerably less than the heights of the barriers shown by Hu and co-workers,⁵⁵ whose predictions were of a barrier height of 25 k_BT. Such high barriers would mean that the system has a very low probability of transitioning between the two states. Such a large barrier could arise due to several reasons. First, Hu et al. use a rigid rod-cylinder model to represent their polymer on a simple cubic lattice. This model only has bond-bond interaction energy without any consideration of trans-gauche conformations and polymer-polymer, polymer-solvent, and solvent-solvent interactions. Second, their estimate of the barrier height is at a temperature below $Tm0$⁠. In contrast, our simulation is an off-lattice simulation, accounting for conformational flexibility and Lennard-Jones interaction between the united atoms. This provides the polymer with far more freedom to explore its conformations and thereby transform from a packed crystal to a molten state. This could explain why we see barriers of considerably smaller heights than those of Hu et al. Even more significantly, our melting kinetics is followed at temperatures higher than the equilibrium melting temperature. In general, according to the classical nucleation theory, when a bulk crystal is superheated above its melting point, there is a nucleation barrier for the birth of droplets of the liquid phase. The stability of the nuclei is given by the melting spinodal, at which the nucleation barrier vanishes. For the single-chain crystal studied here, we are unable to record any such intrachain nuclei of the liquid phase, as evident from Fig. 1 (Multimedia view). The barrier in Fig. 4(b) is distinct and not thermal noise, and thus, it cannot be attributed to melting spinodal. It is also to be noted that the temperature range in Fig. 4(b) is very wide (450 K in experimental units!), and the relative insensitivity of F^*/k_BT to temperature (Table I) emphasizes that the barrier for melting is entropic in origin.

We now address the melting behavior of multichain crystals. The multichain crystal was grown by sequentially growing the crystal one chain after another as is explained in detail in our previous work.⁶⁹ Upon completion of the growth procedure, the multichain crystal contained 21 polymer chains, each containing N = 200 monomers. The crystals were melted isothermally by allowing them to equilibrate at their melting temperatures T_m for several thousands of time units through the use of Langevin dynamics similar to that for single crystal melting.

We choose seven different melting temperatures T_m—11, 11.6, 12, 14, 16, 18, and 20 to judge the change in the conformation of the multichain crystal and the kinetics of disassembly with respect to increasing driving force for melting. As an example, we show the trajectories of two such melting temperatures T^* = 12 and T^* = 18 in Fig. 5 (Multimedia view). The visualization has been rendered using the Visual Molecular Dynamics (VMD) package.⁷³ The 21st chain has been marked in orange in Fig. 5(a) (Multimedia view) and in blue in Fig. 5(b) (Multimedia view) to understand the effect that a rise in temperature has on the outermost chain in the aggregate.

Analogous to the melting of the single crystal, we observe differences in the melting of the multichain crystal with increasing temperature. At melting temperatures close to the equilibrium melting temperature $Tm0$⁠, like the one shown in Fig. 5(a) (Multimedia view), it was found that even after running the Langevin dynamics simulation for several thousands of time units, the multichain crystal does not separate into its individual chains but remains closely packed together in a dense disordered state. This is akin to the mesomorphic state that Hafele et al.¹⁸ seem to observe in their recrystallization experiments on random polyethylene-co-octene and the inhomogeneous intermediate state that Muthukumar³ invokes to explain the same phenomenon. Additionally, the entire crystal seemed to lose its order uniformly.

When the temperature was increased, it was found that the system escaped out of its dense compact shape into a free state, where the individual chains peel off from the core and behave as independent chains. An example of this is shown in Fig. 5(b) (Multimedia view), where the chains begin to disengage from the crystal toward the end of the simulation. In essence, the melting process seems to comprise of a step where the crystal loses its order and then expands in the next step. The accessibility of the latter step is dictated by the melting temperature, analogous to the melting of the single crystals.

We studied melting at temperatures even closer to the equilibrium melting temperature $Tm0$ to further the longevity of the metastable state. As described previously, the multichain crystals were allowed to equilibrate at temperatures of T^* = 11 and T^* = 11.6. At a temperature of T^* = 11, we observed that the multichain crystal continues to remain in its crystalline state, while at the increased temperature of T^* = 11.6, we observed the aggregate to be in a partially ordered state. We show representative trajectories in Fig. 6 (Multimedia view) and describe the evolution of disordering in the system next.

To explore the time scales involved in how the crystalline order of the polymer changes with temperature, we show the time evolution of the global order parameter averaged over all of the chains in Fig. 7. We show the global average of the orientational order parameter at different melting temperatures. Figure 7(a) shows the change in crystalline order for temperatures close to the equilibrium melting temperature, while Fig. 7(b) shows the evolution of crystalline order for temperatures well above the equilibrium melting temperature. The insets in these figures represent the temperatures for which simulations were run for much longer times [T^* = 11 for Fig. 7(a) and T^* = 12 for Fig. 7(b)]. They are shown separately to clearly distinguish the time scales involved. It can be seen from the figure that the crystalline order drops to negligible values very quickly for T_m = 14, 16, 18, and 20, while the inset [Fig. 7(b)] shows that a long time is required at T^* = 12 for the aggregate to become completely disordered. The rate at which the crystal loses order can be seen to increase with melting temperature. Moreover, Fig. 7(a) shows the evolution of the crystalline order much closer to the equilibrium melting point. At the temperature of T^* = 11 as seen in the inset of Fig. 7(a), the system stays in a state of order that is very close to its original state of order. More interestingly, at a temperature of T^* = 11.6 as seen in Fig. 7(a), the system equilibrates to a state of partial order denoted by the value of P₂ at about 0.35. To quantify such metastable states, we then proceed to study the kinetics of melting of the multichain crystals.

The time evolutions of the orientational order parameter P₂ as given in Figs. 7(a) and 7(b) are presented in Fig. 7(c) on a semilogarithmic plot. For high melting temperatures (T^* > 14), P₂ decays with time essentially as a single exponential. On the other hand, at temperatures of 11, 11.6, and 12, when the system stays in the metastable intermediate state, the kinetics is quite complex, representing a hierarchy of local back-and-forth segmental kinetics. For such temperatures, which are close to but above $Tm0$⁠, the time evolutions appear qualitatively as a stretched exponential (⁠$e−tβ$⁠), with the exponent β roughly $13$⁠.

To quantify the free energy of the melting of multichain crystals, we have performed parallel tempering Langevin dynamics simulations similar to those for single crystal melting. The details of the parallel tempering simulations are shown in Table I of the supplementary material. Each melting temperature was explored by running a six-temperature tempering at T = T_c, T_m − 0.4, T_m − 0.2, T_m, T_m + 0.2, T_m + 0.4. The crystalline conformations were allowed into the simulation so as to ensure that swaps with crystalline regions were also possible. We computed the free energy landscape along the m_amorphous reaction coordinate using Eq. (12). The number of amorphous monomers were discretized using bins of size 100, and the bins were then transformed to free energy using a histogram method similar to the one used by Welch and Muthukumar.⁵⁸ The free energy estimates are shown in Fig. 8. The crystalline configurations at the beginning of the simulation have m_amorphous in the range of 500–800. Therefore, the free energy well in the region around 500 in Fig. 8 corresponds to the starting crystalline conformations. The points with m_amorphous < 500 correspond to the low probability that annealing of the crystal at T_c resulted in crystal thickening, causing the total number of amorphous monomers to reduce. The low probability of this step is reflected in the very high free energy of the points below m_amorphous = 500.

The free energy landscapes shown in Fig. 8 show the presence of two barriers, and they have been tabulated in Table II. First, all the free energy landscapes show a free energy barrier at about m_amorphous = 1000 (⁠$F1*$⁠). The free energy landscape at T^* = 11 shows that the slightly molten crystals of Figs. 6(a) (Multimedia view) and 7(a) correspond to the valley in free energy at about m_amorphous = 500 in Fig. 8. The system cannot escape the barrier at this temperature. As the temperature is increased, we find that the system has to navigate through a barrier of finite height. The heights of these barriers as tabulated in Table II show that these barriers are of nearly equal height, relative to k_BT.

Next, the free energy landscapes of Fig. 8 also show a second barrier (⁠$F2*$⁠) near the middle of the landscape (around m_amorphous = 2500). The heights of these barriers have been tabulated in Table II. From these values and from Fig. 8, we can clearly see that the barrier height decreases with increasing melting temperature. At a temperature of T^* = 11, this barrier is too high for the system to surmount. When the temperature is increased to T^* = 11.6, the barrier reaches a finite size of ∼2k_BT. When the temperature is increased further, the barrier height drops even more, before disappearing at temperatures above T^* = 14. The positions of these barriers (near the middle of the reaction coordinate) prove the existence of the partially molten long-lived metastable states, which were shown earlier in Figs. 6 (Multimedia view) and 7(a). The overall mechanism of the crossing of the two barriers at the temperatures of 11, 11.6, and 12 is quite rich and involves a hierarchy of local dynamics, as evident in Fig. 7(c). On the other hand, the melting mechanism at melting temperatures higher than 14 essentially follows an exponential kinetics, as seen in Fig. 7(c).

We have studied the melting of single molecule crystals and multichain crystals in solutions by the use of Langevin dynamics simulations. We have performed this using our coarse-grained united atom model, which was generated after making modifications to a model that represented polyethylene.⁶³ Using this system, we have observed the responses of the semicrystalline polymer systems to melting by allowing preformed crystals to equilibrate at temperatures above the equilibrium melting point, which was determined in our previous work.⁶⁹ To analyze the changes in the morphology of single crystals during melting, we have monitored the ensemble-averaged radius of gyration and the global order parameter. Additionally, we have developed a new order parameter by locally averaging the second order Legendre polynomial P₂ and determining if a particular monomer can be categorized as a crystalline or amorphous monomer. This idea has been adapted from some previous works^55,68 and modified to suit our system. Using this new order parameter, we have computed free energy landscapes using parallel tempering Langevin dynamics simulations.

The melting of single molecule crystals is found to be qualitatively similar to the experimental results of melting of single polyethylene mats by Barham and Sadler,²⁴ when the crystals undergo a transition to a coil-like state. However, we find that the melting process is comprised of two steps: one where the crystal loses its order and the other where the globulelike chain expands into a coil-like one. We find the former step to occur at every temperature studied, but the latter is accessible only when the temperature is increased relative to the equilibrium melting temperature. However, we did not find any evidence to suggest preferential directions for melting, as seen in several experiments.^{33–36,42,43} We suspect that this might be because of constant temperature provided to all monomers by the Gaussian random noise in our simulations. The kinetics of the melting process was quantified after sampling the free energy from parallel tempering Langevin dynamics simulation over 6 different temperatures (one of which was constrained to sample the crystalline configurations alone), and these revealed the presence of long-lived partially crystalline metastable states during the melting process. This has been observed in experiments as well.^18,61 We quantified the free energy barrier and found that there was one significant barrier in this process for single chain crystals. The position of the barrier, interestingly, was at the middle of the reaction coordinate, and it had a height of approximately 2 k_BT.

For the melting of multichain crystals, we find that the crystal initially loses order uniformly, without any preferred direction. This is similar in feature to the melting of the single molecule crystal. The melting process is comprised of a step where there is a transition to a coil-like state formed before transforming to a free state where chains begin to peel off the crystal. The accessibility of the second step was dependent on the temperature, with higher temperature providing easier access. These simulations revealed that the system escapes its native state to a state of intermediate order at temperatures close to the equilibrium melting temperature, as addressed earlier by Muthukumar in previous work,³ consistent with the hypothesis of mesomorphic state by Strobl.^18,61 The extent of this intermediate ordering decreases with increasing temperature. We quantified the free energy landscape using parallel tempering Langevin dynamics simulations. Simulations at all of the melting temperatures revealed a free energy barrier closer to the beginning of the reaction coordinate. At temperatures closer to the equilibrium melting temperature, the simulations revealed a second barrier near the middle of the landscape. These indicate that the partially molten metastable states are preferred by the system above the equilibrium melting temperature, with the preference for such states decreasing with increasing temperature consistent with the earlier theory by Muthukumar on melt memory.³ At a certain temperature above the equilibrium temperature, the barrier near the middle of the reaction coordinate vanishes and the system prefers to be in a completely molten state.

See supplementary material for movies of melting of single molecule crystals at T^* = 12 and T^* = 18 and for videos of melting of multichain crystals at T^* = 11, T^* = 11.6, T^* = 12, and T^* = 18.

We have also explained the details of the parallel tempering procedure used, followed by details of how we have calculated the number of amorphous monomers m_amorphous.

This research was supported by the National Science Foundation (Grant No. DMR-1713696) and AFOSR (Grant No. FA 9550-17-1-0160).