Molecular dynamics simulations of crystal nucleation in entangled polymer melts under start-up shear conditions

Understanding the ﬂow induced crystallisation (FIC) process is necessary due to its technological relevance to polymer processing. Polymer crystallisation controls the morphology of semi-crystalline polymers and hence the properties of the end product. We perform molecular dynamics simulations of polymer melts consisting of suﬃciently entangled linear chains under shear ﬂow. We determine the Rouse relaxation time ( τ R ) for linear polymer chains using an established rheological model at diﬀerent temperatures and ﬁt the simulation data with the Arrhenius and Williams-Landel-Ferry (WLF) equations. We simulate the crystallisation induction times for diﬀerent values of the Rouse Weissenberg number (W iR = ˙ γτ R ) at diﬀerent temperatures. We observe that the level of strain and stretch required to induce crystallisation increases with temperature. We ﬁnd that the induction times follow a power law in shear rate and observe a more pronounced eﬀect of ﬂow rate for higher temperatures than at lower temperatures. Moreover, we determine that nucleation events occur relatively early in the shear transient and at a stretch value that is smaller than its steady state value. We also report the values of strain at which the occurrence of a nucleation event is most likely to happen.


I. INTRODUCTION
The commercial importance of polymer processing means that polymer crystallisation remains an active and vital field even after extensive efforts over several decades from various fields of science, engineering and technology. This topic is not only interesting scientifically but its understanding is also necessary for industry to produce products of desired mechanical, thermal and optical properties. Polymer crystallisation defines the morphology of semi-crystalline polymers, which has a strong influence on many of the properties of the end product. Therefore, in order to achieve the desired product properties, we must have better control over the crystallisation process in order to select the morphology of these materials. This requires a thorough understanding of the crystallisation process.
Experimentally, it is extremely difficult to detect flow induced crystallisation at its earliest stage because the small spatio-temporal scales are below the resolution of available experimental techniques 33 . In such scenarios, computer simulations are a good candidate to investigate the problem. On a mesoscopic level, kinetic Monte Carlo 34 simulations have been performed and the effects of shear rate, shear time and shear strain on the crystallisation kinetics, morphology and rheology have been discussed. Graham and Olmsted 35,36 have developed a coarse-grained kinetic Monte Carlo model. It can capture many features of flow induced crystallisation and is able to reproduce the effect of flow rate on crystallisation kinetics, which is in a good agreement with experimental data. However, it is a phenomenological model which makes pre-assumptions about the nucleus shape and relies on unknown thermodynamic parameters. Therefore there is a strong need for detailed molecular simulations with full resolution down to the monomer scale. On a molecular level, Monte Carlo (MC) and molecular dynamics (MD) simulations have been performed to address polymer crystallisation under quiescent conditions  and under flow or large deformation [61][62][63][64][65][66][67][68][69][70] . More recently, Graham 71 has reviewed and discussed the recent advances and future directions for a better understanding of FIC, with an emphasis on the role of molecular simulations. This review covers advantages, disadvantages and limitation of experiments, theory and multiscale simulations of FIC. In the case of flow induced crystallisation, there are very few MD studies 70,72 that quantify flow in-duced nucleation. These studies focused on short chain alkanes of length 20 monomers and 150 monomers. They reported the nucleation mechanism and kinetics for these short chains. However, longer chains that are sufficiently entangled have not been investigated yet using molecular dynamics simulations. Therefore, in this work, we report flow induced crystallisation of linear polymer chains consisting of 1000 monomers (entanglements ∼ 13) under different degrees of super-cooling (4-14%) and at different values of the Rouse Weissenberg number . Another distinct feature of this work is that we use a rheological model to contextualise our data, because flow induced crystallisation and polymer rheology are strongly linked.

A. Model
We use a united atom model, which has been used extensively by Rutledge and co-workers [41][42][43]70 and Schilling and co-workers 45,72,73 to study crystallisation under quiescent and under flow conditions. This model was proposed by Paul et.al. 74 and then modified by Waheed et.al 56 . In this model, CH 2 and CH 3 groups are represented by beads or "united atoms". These beads interact with each other via bonded and nonbonded potentials. The nonbonded interaction consists of Lennard Jones interactions and can be expressed by following relation for a distance r ij between monomers i and j: The bonded potential, which acts between monomers along the chain, consists of a harmonic bond potential a harmonic bond angle potential where θ is the angle between two consecutive bonds, and a dihedral potential (4) where φ is the dihedral angle defined by three consecutive bonds.
We provide parameters for all potentials in Table I.  The model has been optimised to reproduce the dynamical and structural properties of the melt, the melting  point, and the rotator phase properties.  In order to identify the crystalline regions within a melt, we use a crystallinity order parameter that has been used recently by Schilling and co-workers 45,72,73 . This order parameter is based on the local alignment of segments of chains.
First of all, we identify the neighbouring particles j of a given particle i, which lie within a cut-off radius of r c = 1.5σ from particle i, with particle i and j not on the same chain.
We associate a unit vectorê to every particle i pointing from the particle centre of i − 1 to the centre of i + 1. Then, we determine the angle between particle i and every neighbour j.
Those particles whose number of "aligned" neighbours exceeds the threshold value (8 monomers in our study) are called crystalline particles. This threshold value is determined from the following analysis of the probability distribution of aligned neighbours in the bulk melt. We draw a probability distribution of aligned neighbours for an equilibrated melt and select the threshold number of aligned neighbours to be the value where the probability reaches zero on the right-hand side of the bell shape curve of the probability distribution. This shows that no melt particles in the system have more aligned neighbours than this threshold value. Hence, if any particle has more aligned neighbours than this threshold value, it is a crystalline particle.
Finally, the clusters of crystalline particles are identified using a standard clustering algorithm. The crystalline clusters in the system are identified by picking a random particle and checking if it is crystalline or not. If it is crystalline, we count it as the first particle of this cluster and examine its neighbours. If any neighbouring particle is also crystalline, it is counted as the second particle of the same cluster. Similarly, we move recursively from neighbour to neighbour to compute the cluster size. If a particle does not have any new crystalline neighbours, then we move to the next particle to identify the second cluster in the system and so on. At the end, all cluster sizes are compared and the largest cluster is identified.
We note here that, in our simulations whenever a nucleus that is significantly bigger than the critical nucleus (as define later) is seen, this always develops into a large crystalline structure.

C. Simulation details
We have performed molecular dynamics simulations of linear polymer chains to study nucleation under shear flow. The system consists of 300 chains of length 1000 monomers. We equilibrated the system at 550 K which is well above the equilibrium melting temperature (396.4 K) of polyethylene. This value of the melting temperature has been determined experimentally 76 and this force field reproduces the experimental melting temperature for short chain alkanes 42 and has been used for longer chains in previous studies 43 . We consider a system to be equilibrated once the chains have diffused a distance equal to their radius of gyration. Furthermore, we ran one very long equilibration trajectory of approximately three times this diffusion time. Our nucleation under flow results these two different lengths of equilibration are indistinguishable. We then run NPT simulations at temperatures at which we want to see nucleation and let it relax to the density that corresponds to 1 atm pressure. The system remained amorphous during these density relaxation simulations as the density relaxation at T c , which was performed quiescently, took significantly less time than the flow-induced induction times. After getting these relaxed configurations, we fix the density again and raise the temperature to 550 K and let the melt relax at this temperature. We did this because we wanted to quench a fully relaxed system to the desired crystallisation temperature at the correct density. The densities of the metastable melt at 1 atm pressure and the corresponding temperatures are given in Table II. All shear simulations have been performed under constant volume and constant temperature conditions. The shear flow has been generated using Lees Edwards boundary conditions 77 and the DPD thermostat 78 in the same way as Schilling and co-workers 72 . We took the friction coefficient for the thermostat to be 0.5 τ −1 and the cutoff radius to be R c = 1.3σ, where τ = kBT mσ 2 and m is the mass of the bead, k B is the Boltzmann constant, T is the temperature in Kelvin and σ is the size of the beads. We used the ESPResSo 79 molecular dynamics package for all simulations.
During nucleation simulations, we quenched the equilibrated configurations from 550 K to 380 K, 360 K and 340 K and applied the shear rate (γ). We run simulations at three different temperatures and at three different values of the Rouse Weissenberg number (W iR =γτ R ). The integration timestep used in the simulations at 340 K was 0.0066τ . The factor 0.0066 was changed at different temperatures to keep the time step the same in seconds. We choose the flow timescale (γ) to be slower than the Rouse time of an entanglement segment τ e (γτ e < 1) so that, on the timescale of the flow, monomers have enough time to experience their surrounding constraints and entanglements have a strong influence on how the flow deforms the polymer chains. This is because any substantial strain ( γ > 1) takes longer than the Rouse time of an entanglement segment and so the response to flow is always influenced by entanglements. In contrast, the timescale of a successful nucleation event is shorter than τ e and so nucleation, itself, is not directly influence by entanglements. We define the timescale of a successful nucleation event as follows: we begin with a well formed crystal and run time backwards through the saved trajectory until the crystal is the critical nucleus size and call this τ crit . We then continue to run time backwards, stopping when no monomers in the critical nucleus are crystalline and label this τ start . This is the time just before the successful nucleation event begins. The time τ crit − τ start measures the time from inception to stable nucleation for a single successful nucleation event. This observation is supported by quiescent nucleation simulations by Yi et. al 43 . They saw the same nucleation barrier for both entangled and unentangled chains, suggesting that nucleation, at least at this level of undercooling, is a local event occurring on a lengthscale shorter than the tube diameter. We used a total of approximately 2000 kCPU hours on the University of Nottingham and Midland Plus HPC facilities to carry out these simulations (equilibration, rheological characterisation and nucleation simulations).

A. Rouse relaxation time
Shear flow is expected to significantly change the nucleation rate above a shear rate that corresponds to a Rouse Weissenberg number of order 1. To determine the Weissenberg number, we must know the Rouse relaxation time of our system. There are many methods to determine the Rouse time from molecular dynamics simulations. Recently, Rutledge and co-workers 43 estimated the Rouse time from equilibrium molecular dynamics simulations by computing the end to end vector auto-correlation function. Cao and Likhtman 80 determined the Rouse time from equilibrium molecular dynamics simulations by calculating the stress relaxation auto-correlation function. We wish to calculate the Rouse relaxation time through stress simulation data to be consistent with rheological models and measurements. The tube model is widely used to study the rheology of polymers and it separates the molecular deformation into chain stretching and orientation. Chain stretching defines the Rouse time (τ R ) and chain orientation defines the reptation or terminal time (τ d ∼ 3Zτ R ), where Z is the number of entanglements. We found that the usual method to extract the relaxation time based on equilibrium simulations of the stress relaxation auto-correlation function was prohibitively expensive, requiring a simulation of ∼ 10 τ d to get good statistics for the terminal time 81 . Therefore, we devised a new way to extract the Rouse relaxation time based on non-equilibrium molecular dynamics simulations which requires a simulation of length ∼ 4 τ R . We describe this method in detail here.
In this method, we run non-equilibrium molecular dynamics simulations at 450 K to compute the transient shear viscosity, which we show in FIG. 1. We run simulations at shear rates ofγ = 0.000002τ −1 anḋ γ = 0.0001τ −1 . In this figure, the simulation data are shown with symbols and the GLaMM model 82 fitting is represented with lines. In the GLaMM model, three input parameters are used to obtain the best fit. These parameters are the entanglement modulus (G e ), the Rouse relaxation time of an entanglement segment (τ e ) and the number of monomers per entanglement (N e ). We used the standard value of 0.1 for the constraint release parameter (c ν ). The Rouse time can be calculated using the relation τ R = τ e Z 2 , where Z = N/N e and N is the number of monomers in a chain. We have also computed N e using the Z1 code [83][84][85][86] , which uses a primitive path analysis approach to estimate the number of chain entanglements. We find the entanglement length to be 72 monomers (Z ∼ 13.84), which is in very close agreement with the rheologically determined value.
Once we have the Rouse time at a single temperature, we can then estimate the Rouse time at all other temperatures by running short equilibrium molecular dy-  namics trajectories at the desired temperatures. From these trajectories, we calculate the stress relaxation modulus (G(t)) and then we shift all curves to the reference temperature curve based on the time-temperature superposition (TTS) concept 87,88 using the same approach as widely used for experimental data [88][89][90][91][92][93][94] . We apply TTS to G(t) in the following way: where T 0 is the reference temperature (T 0 = 450K in this work), t is the time, T is the temperature from which G(t) is to be shifted and b T and a T are the vertical and horizontal shift factors respectively. The value of the vertical shift b T is specified by the densities: where T 0 and ρ(T 0 ) are the temperature and density to which the data is shifted. All densities were determined by NPT simulations and are in Table II. The temperature dependence of the horizontal shift factor a T can be expressed using the Arrhenius equation, given as follows: where E a is the activation energy of flow and R is the universal gas constant. This relation gives a good fit of the data in the plateau and terminal zones for temper-atures well above the glass transition temperature T g , which has been reported as 223.04 ± 0.22K 43 for this force field. The temperature dependence of the horizontal shift factor a T for a temperature range between T g and T g + 100 K can be expressed by the Williams-Landel-Ferry (WLF) equation 91,95 : where c 1 and c 2 are material parameters. We determined the value of the horizontal shift factor a T by manually shifting the G(t) curves to the G(t) curve for the reference temperature T 0 = 450K, which generated a master curve for G(t).
The Rouse time (τ T R ) can then be estimated for different temperatures from the horizontal shift factor (a T ) and the Rouse time at the reference temperature (τ T0 R ) as follows: In FIG. 2(a), we show simulation data for the stress relaxation modulus (G(t)) at different temperatures for a short period of time. We select the curves from the highlighted region because the curves are statistically rich in this region when compared to the later time curves. We already have the Rouse time at 450 K, which we have calculated from non-equilibrium simulations and the GLaMM model fit. Now, we have to shift the curves vertically using EQU 7 and then move the curves horizontally to superpose the curves on the reference curve (450 K). In FIG. 2(b), we show the shifted and superposed curves from different temperatures and the reference curve at 450 K. From here we obtain different values of the horizontal shift factor (a T ) for different temperatures which we show in FIG. 2(c) and in tabulated form in Table II. We also show fits to the data of the horizontal shift factors (a T ) using the Arrhenius equation (EQU 8) and the WLF equation (EQU 9) with solid black and blue lines respectively. The value of the activation energy (E a ) was found to be 31 KJ/mole for the best fit. This value is very close to the experimentally reported value (26.77 KJ/mole) 96 and numerically reported value (21.46±3.30KJ/mole) 43 . A detailed discussion of the variation of the activation energy with different measurement techniques can be found in these works 93,94 . The Arrhenius equation fits well for temperatures well above the glass transition temperature, while the WLF equation fits better for the whole range of temperatures with material constants c1 = 2.05 and c2 = 300K.   We perform molecular dynamics simulations at three different temperatures for three different values of the Rouse Weissenberg number to calculate the induction time (τ * ). In nucleation studies, estimation of the induction time is an important step particularly from molecular dynamics simulations. This involves the appropriate selection of an order parameter to identify the largest crystal cluster in the melt from the trajectories and then selection of a method to estimate the induction time. We use the crystallinity order parameter that we described in section II B and we use a mean first passage time (MFPT) approach 97 , which is based on classical nucleation theory, to estimate the induction time. Mean first passage time analysis is performed on the evolution of the largest cluster in the system to define the average time of the first appearance of a cluster with size n max : where M is the total number of trajectories and τ (i) nmax is the time when a cluster with size n max first appears. As nucleation is followed by fast cluster growth, τ (n max ) has a sigmoidal shape and can be fitted by the equation: where τ * is the induction time, n * is the critical nucleus size, Z is the Zeldovich factor and the error function is erf = 2 √ π x 0 e −x 2 dx. This mean first passage time has been successfully used for this purpose in many recent studies [41][42][43]45,72 .
It is a common perception that flow enhances the nucleation rates by inducing local orientation of polymer chains. Under flow, temperature plays an important role in two ways. Firstly, at higher temperatures, the relaxation times are shorter, meaning lower orientation. Secondly, temperature affects the quiescent crystallisation kinetics. Therefore, the coupled effect of flow and temperature needs to be addressed. The effects of flow and temperature on the induction time have been studied experimentally 15,98 and theoretically 15 . These studies show that the effect of flow on the induction time can be divided into two distinct regions. For low Weissenberg numbers there is no effect of shear flow on the induction time, while for high Weissenberg numbers, the induction time decreases as a power law in shear rate. Furthermore, the slope of the curves in the second region (W iR ≥ 1) at higher temperatures is greater than the slope at lower temperatures. This means that at very high shear rates, the induction times for different temperatures start to converge. We show experimental observations of these effects in FIG. 4(b), which has been taken from the literature 15 . This same effect has been seen in the Graham and Olmsted 35,36 kinetic Monte Carlo model. In this model, the increased sensitivity to flow occurs because higher undercooling produces larger critical nuclei, which can incorporate more stretched chain segments 99 .
We computed the induction time, as show in FIG. 3. This figure shows that the induction times are lower at higher temperatures for the same Rouse Weissenberg number (W iR ) which is contrary to the observations made in experiments 15,98,100,101 and theory 15,34 . However, when we incorporated the monomer friction by dividing the induction times by Rouse times (τ R ) (see FIG. 4(a)), we too observe the same trend as reported in the previously published literature from experiments and Monte Carlo simulations, namely a more pronounced effect of flow at higher temperatures. We can explain need to divide the induction time by the Rouse time using classical nucleation theory. According to classical nucleation theory, the nucleation rate depends on the product of a kinetic factor and the Boltzmann factor of the nu- cleation barrier height. The former is proportional to the monomer friction. In experiments, usually the degree of supercooling is very low and experiments are run at temperatures that are far from T g and have only a few degrees difference. Therefore, this monomer friction factor does not change significantly between these temperature. In contrast, in our molecular dynamics simulations, the kinetic prefactor changes more strongly between the different temperatures. This is because the degree of supercooling is relatively large so we are much closer to T g and the differences between the temperatures at which simulations are run are high (20K difference in our case). Therefore, we divide induction time by the Rouse time in MD simulations to correct for the strong dependence of the induction time on the monomer friction. We fitted power laws to the data from each temperature to obtain the exponent with shear rate, as shown in FIG 4(a). The window of induction times we were able to access in MD was too limited to observe the cross-over to quiescent behaviour as in experiments. Longer simulations, perhaps utilising a fast nucleation algorithm 102,103 , will be needed to see this effect. Our induction time data can also be converted to an induction strain γ I , computed by multiplying the induction time by the shear rate. The induction strain can be considered as the level of shear strain required to start crystallisation. We show the induction shear strain against the Rouse Weissenberg number from our simulations in FIG. 5(a). As expected, the induction strain increases with temperature. However, the dependence on shear rate shows some counter-intuitive behaviour. At the highest temperature γ I decreases with shear rate, at the intermediate temperature γ I is almost constant and at the lowest temperature γ I increases with shear rate. Nevertheless, these qualitative features have been seen in experiments 100 as shown in FIG. 5(b). In section IV, we offer an interpretation of this observation, by considering the degree of chain stretch at nucleation.

IV. INTERPRETATION OF NUCLEATION DATA
The value of γ I represents the macroscopic strain at the time of nucleation. However, the microscopic strain on the chain level is likely to be different due to chain relaxation. The microscopic strain can be observed via the chain stretch ratio, defined as the fractional increase in the tube contour length (Z * (t)/Z) as computed from the GLaMM model 82 . In order to compute the stretch from the GLaMM model, we used the parameters that we extracted in the relaxation time calculations by fitting the GLaMM model to our MD simulations. We plot the chain stretch, from the GLaMM model, at the induction time, from MD simulation, in FIG. 7. This induction stretch, λ I , shows three key features: λ I increases with increasing temperature; λ I increases weakly with flow rate; and the increase in λ I with flow rate is steeper at a)   lower temperatures. We interpret these data by postulating that the induction time comprises of two contributions: the time required to achieve sufficient stretch that nucleation becomes highly likely, τ λ * , minus the contributions from the nuclei that manage to cross the barrier before λ * is reached, τ early . This explains the three key features of FIG. 7: λ I increases with temperature because the increased quiescent nucleation barrier requires greater stretch to sufficiently lower the barrier (i.e. λ * increases with temperature); λ I increases with flow rate because, at slower flow rates, the longer time to reach λ * permits more early nucleation events (i.e. τ early decreases with increasing flow rate); and the lower barriers at lower temperatures allow more nucleation prior to λ * , giving a larger role for τ early . Indeed FIG. 7 shows that λ I is virtually constant at 380K, suggesting an insignificant effect of τ early , whereas somewhat steeper increases in λ I are seen at lower temperatures. This interpretation explains the counter-intuitive results for the variation in induction strain with shear rate in figure 5 as a combination of the strain required to obtain the critical chain stretch and the probability of early nucleation events before this critical stretch is reached. We show histograms of the observed induction times (computed from MD simulations) and stretch (computed from the GLaMM model 82 ) with red bars and blue lines respectively in FIG. 8. The stretch curve shows that all nucleation events occur in transient conditions and the stretch value is smaller than the steady stretch value. Further data for the stretch and induction times are shown in appendix B. From figures (FIG. 6 and FIG. 8), we can deduce that nucleation happens early in the transient, when the molecular stretch is significantly below both the transient maximum and steady state value. This suggests that at this under-cooling, the flow very readily induces nucleation even at fairly modest molecular deformation. This illustrates the type of information that is difficult to extract from experiments but can be obtained from MD.

V. CONCLUSIONS
We carried out united-atom molecular dynamics simulations to characterise the nucleation behaviour of entangled linear polymers under shear flow. We studied well-entangled C1000 chains, for which we determined the Rouse relaxation time by fitting a tube model to shear stress data from simulations of non-linear start-up shear. We characterised the temperature dependence of this timescale via time-temperature superposition. From our crystallisation simulations we computed the induction time at a range of Weissenberg numbers and for several temperatures. We found a more pronounced effect of flow rate on the induction times at higher temperatures. Also we observed a weak dependence of the induction strain on shear rate, with the slope of this dependence decreasing and changing sign with increased temperature. Both of these effects are also seen experimentally and our results provide a direct confirmation that these phenomena can be predicted from molecular dynamics simulations. We presented results for the distribution of induction times, compared to the shear stress transient predicted by the GLaMM model. For a given temperature and shear rate the distribution of induction times is quite narrow. Furthermore, all nucleation events occur early in the start-up flow, before the shear stress overshoot.
By using the microscopic chain stretch at the time of nucleation, as computed by the GLaMM model, we are able to explain the counter-intuitive results for the variation of induction strain with shear rate in figure 5. We postulated that the induction strain is a combination of the strain required to obtain the critical chain stretch