Coarse-grained simulations of flow-induced nucleation in semi-crystalline polymers

We perform kinetic Monte Carlo simulations of flow-induced nucleation in polymer melts with an algorithm that is tractable even at low undercooling. The configuration of the non-crystallized chains under flow is computed with a recent non-linear tube model. Our simulations predict both enhanced nucleation and the growth of shish-like elongated nuclei for sufficiently fast flows. The simulations predict several experimental phenomena and theoretically justify a previously empirical result for the flow-enhanced nucleation rate. The simulations are highly pertinent to both the fundamental understanding and process modeling of flow-induced crystallization in polymer melts.

I.

INTRODUCTION
The nucleation of microscopic crystallites in polymer liquids is profoundly influenced by flow [1,2]. This flow-induced crystallization (FIC), is a fascinating example of an externally driven, non-equilibrium phase transition, controlled by kinetics. Furthermore, FIC is ubiquitous in industrial processing of semi-crystalline polymers, the largest group of commercially useful polymers. A fundamental understanding of FIC promises extensive control of polymer solid state properties, as virtually every property of practical interest is determined by the crystal morphology. Flow can drastically enhance nucleation and trigger the formation of highly aligned, elongated crystals, known as shish kebabs [1]. Recent experiments on entangled polymers, have studied, in detail, shish kebab formation [3,4,5] and the role of blend concentration [6], molecular architecture [7] and molecular relaxation time [8]. Often the most pronounced flow-induced effects occur near the melting point, where quiescent crystallization is immeasurably slow [2,8].
The widely postulated mechanism for FIC states that flow forces the polymer chains into elongated configurations, which lowers the entropic penalty for crystallization [1]. However, this hypothesis has yet to be developed into a quantitative molecular model. FIC is extremely sensitive to the flow-induced configurations of the non-crystalline chains, so an accurate molecular flow model is an essential prerequisite. Unfortunately, most polymer flow models predict only the macroscopic stress tensor and not the full molecular configuration.
Alternatively, detailed simulations of polymer crystallization have provided much useful information on the growth process [9,10,11], yet simulating primary nucleation has proven difficult, especially at low undercooling, because of the extremely long nucleation times. At a much higher level of coarse-graining, models based on differential equations either assume an empirical dependence of the nucleation rate on the flow conditions [12], the stress tensor [13], or the chain stretch [14]; or assert that free energy changes under flow can be directly subtracted from the nucleation barrier [15,16]. In either case the postulated FIC mechanism remains untested. An intermediate level of coarse-graining is required to surmount these difficulties.
This letter presents coarse-grained kinetic Monte Carlo (MC) [17] simulations of anisotropic nucleation in flowing polymers. We compute chain configurations using a recent molecular flow model [18] that reliably predicts both neutron scattering [19,20] and bulk stresses. Our simulations predict both enhanced nucleation and elongated shish nuclei.
Kinetic MC has previously been used to model quiescent crystal growth in dilute polymers [21]. However, it is particularly suited to nucleation and our algorithm is tractable even at low undercooling, providing an efficient and highly flexible framework to simulate general anisotropic nucleation under external fields.

II. MODEL
We compute the transient chain configuration of the uncrystallized chains under flow using the GLaMM model [18], with finite chain extensibility included using Cohen's approximation [22]. The chains are divided into Z sub-chains, each corresponding to an entanglement segment of N e Kuhn steps of length b. We take N e = 100 throughout this Letter. One deterministic run of the model provides the end-to-end vector f i (t) = r i r i , where the ensemble average is for sub-chains of type i. The data for an entire transient flow are used later in the nucleation simulations. All flow timescales are in units of the subchain Rouse time τ e and r is normalized by Deformation of the amorphous chains has two effects on the nucleation kinetics: stretching reduces the entropic penalty for crystallization; and monomer alignment modifies the probability of compatible alignment with the nucleus. The change in elastic free energy ∆F el for chains with ensemble average constraints f = rr , but locally at equilibrium, can be calculated by statistical mechanics [23]. Although an analytic calculation is not possible for finitely extensible chains, steep free energy gradients in highly stretched chains suppress fluctuations. Thus our numerical calculations for uniaxial deformations show that ∆F el ( rr ) can be accurately approximated by an expression that interpolates between Gaussian elasticity [23] for small Tr f and Cohen's [22] approximation with r 2 = Tr f at high stretching, Similarly, numerical calculation of the monomer orientation distribution w(θ) for chains with a constraint f are well approximated by using r 2 = Tr f in the expression for w(θ) derived from a direct constraint on r [24].
where L −1 is the inverse Langevin function and θ is the angle between the monomer and the principle axis of f. Each subchain has an individual f i so is treated as a separate species with concentration φ i = 1/Z. Our coarse-grained simulations use the minimal nucleus description required for anisotropic nucleation. The nucleus comprises N T crystallized "monomers" or Kuhn steps, arranged in stems, with each stem formed from a single chain (see fig 1). The total number of stems N s and the number of monomers in each stem is simulated. The arrangement of the monomers within the crystal is not resolved. The nucleus is assumed to be spheroidal with the polar radius L parallel to the stems. Assigning a crystalline volume of b 3 0 to each monomer and normalizing all lengths by b 0 gives the equatorial radius W = N s /π, the volume V = N T and thus the polar radius L = 3N T 4Ns . We also simulate the unit vector v parallel to the polar radius. As in classical nucleation theory the nucleus free energy comprises the free energy of crystallization, proportional to the nucleus volume, and a free energy penalty proportional to the surface area S. Thus the free energy in units of k B T is where ǫ B and µ s are the coefficients of the volume and surface area free energies, respectively.
We simulate two types of MC moves: addition of a new stem containing one monomer and lengthening of an existing stem by adding a new monomer, each with a corresponding reverse move. As in [24] we assume that to attach, a monomer must be oriented within a solid tolerance angle Ω of the nucleus orientation. For small Ω, the fraction of monomers within this angle is w(θ)Ω, where θ is the angle between the nucleus polar radius and the sub-chain principle strain axis. The stem attachment rate for species i is proportional to its melt concentration φ i . In contrast, for stem lengthening only the next monomer along the chain forming the stem can crystallize. The concentration of this monomer at the nucleus   surface where the lengthening event occurs is taken to be unity. Thus the attachment (+) and detachment (−) move rates for stem addition (k st ) and stem lengthening (k len ) are: where τ 0 is the time for a monomer attachment attempt, ω i = 4πw i (θ i ) is the fraction of correctly aligned monomers (normalized to unity in the quiescent limit), and the constant ln(Ω/4π) has been added to ǫ B . The free energy change of attaching a new stem of species where ∆F el i is the flow-induced free energy change in sub-chain i. Similar calculations give the free energy changes for the other move types. A stem can only be detached if it contains a single monomer.
The kinetic MC algorithm requires a sum over all possible move rates. The area available for stem addition moves is A(N s ), which is taken to be proportional to √ N s , and to give spherical nuclei in the quiescent limit. To obey detailed balance the rate of stem removal must be multiplied by A(N s − 1)/N s , the probability of a given stem being at the nucleus surface. Each stem can lengthen or shorten from either the top of bottom. Thus the total sum over all possible move rates is: At each kinetic MC timestep one move is performed at random, with the selection probability weighted by the move rate. Time is then incremented by a stochastically determined interval After each MC step the nucleus orientationv is incremented over ∆t by Brownian dynamics. Flow rotates the spheroid nucleus through the Jeffery algorithm [26]. The angular diffusion time τ rot scales with both the nucleus volume N T and aspect ratio ρ following the expressions in [26], giving

III. RESULTS
Each simulation evolves a single nucleus from a single monomer. The algorithm is especially effective at low undercooling, as for small nuclei the rate sum is small, leading to large time steps. The quiescent free energy landscape can be calculated analytically from ǫ B and µ S to give a dimensionless nucleation barrier ∆f * and a critical nucleus of n * monomers. The simulated nucleation time τ N is the first time the polar and equatorial radii simultaneously exceed the critical radius r * = 3 3n * /4π. Choosing a larger threshold size for nucleation has little effect on τ N . The results are accurately approximated by τ N = τ 0 exp(∆f * ).
We obtained good statistics for barriers up to 25k B T in ∼ 50hrs on one 2GHz processor, giving a nucleation time of ∼ 10 11 τ 0 . The nucleation times are Poisson distributed, so the nucleation rate can be defined asṄ 0 = 1/ τ N .
Under flow, the nucleation kinetics depend on the evolving chain deformation. We define the instantaneous nucleation rateṄ (t) ≈ 1 where n(t) is the cumulative fraction of runs nucleated at time t. For shear ratesγ that are slow compared to the critical nucleus rotation time (τ n * rotγ ≪ 1) alignment effects can be ignored, equivalent to taking α = 0. The GLaMM model predictions for start-up of constant shear atγτ e = 0.1 of a Z = 25 melt and the resulting instantaneous nucleation rate are in fig 2(a). The results are independent of the ratio of flow and crystallization timescales for S 1. In fact, fixing the non-crystalline chain configuration to that corresponding to flow time t for the entire simulation and plotting the resulting quasi-static nucleation rate against t, reproduces the transient results (fig 2a). Thus nucleation is fully controlled by the instantaneous configuration of the surrounding chains, if τ e > τ 0 .
i is the dominant factor determining the nucleation rate, despite variations in flow rate, molecular weight and flow geometry (Fig 2(b)). This universal result is somewhat surprising as the distribution of stretch along the chain varies considerably with flow conditions, which may be expected to influence nucleation, especially if molecular weight and flow geometry is varied. This result will be useful in deriving simple differential models of FIC [13,14] as the nucleation rate can be described the ex-pressionṄ =Ṅ 0 exp(η(λ 2 − 1)), where η is a fitting constant. Also in fig 2(b) is the curvė N =Ṅ 0 + β(λ 4 − 1), which empirically fits measured nucleation rates from flowing melts [14]. The agreement with our simulation data for λ 3.5 theoretically justifies this empirical expression. In the inset to fig 2(b) the sensitivity of nucleation to chain stretching increases with decreased undercooling (decreasing ǫ B ), as seen experimentally [27]. an orientational order parameter, the P2 orientation function, against shear rate for a 2% high molecular weight blend [8], highlighting similarities with our predictions. Also in Fig 3 reduced undercooling, from reducing ǫ B , increases the anisotropy, as seen experimentally [29]. When crystallization is less favorable, the disparity in the kinetics of stretched and unstretched chain segments increases.

Our efficient kinetic Monte Carlo algorithm for flow-induced nucleation in polymer melts
is tractable even at low undercooling. The flow-induced nucleation rate is a universal function of the chain stretch ratio, independent of flow rate, molecular weight and flow geometry, but with decreasing undercooling causing increased flow sensitivity. This universal curve is very similar to an empirical relationship that accurately describes flow induced nucleation rate measurements [14], justifying this empirical result. Our successful quantitative comparison with nucleation measurements [25] establishes that changes in chain free energy under flow can describe flow enhanced nucleation in semi-crystalline polymers. In our bimodal blend simulations, a few percent of high molecular weight chains optimizes shish formation and the degree of anisotropy increases with shear rate and decreased undercooling, all of which are seen experimentally [4,5,6,7,8,29]. The simulation can readily be generalized to fully polydisperse melts, an essential step to model industrial polymer processing. The efficiency and flexibility of our algorithm makes it suitable to simulate general anisotropic nucleation under external fields.

V. ACKNOWLEDGEMENTS
We thank the EPSRC for funding (GR/T11807/01) and L Balzano, R Steenbakkers, G

Peters, O Mykhaylyk, T Ryan and T McLeish for useful discussions.
APPENDIX A: COMPARISON WITH NUCLEATION MEASUREMENTS.

Introduction
Here we quantitatively compare our simulation results with directly observed steady-state nucleation rates from a polymer melt under shear [25]. These measurements, by Coccorullo et. al, were made on a commercial grade isotactic polypropylene (iPP), known as iPP T30G.
These data show an exponential dependence of the nucleation rate with shear rate, which our model correctly predicts. In order to quantitatively model this system, the required parameters fall into two classes: molecular relaxation times, which can be accurately calculated from literature values; and crystallization parameters, for which order of magnitude estimates are provided by the literature, but whose specific values we determine from the nucleation data. The high polydispersity of iPP T30G is a further complication. As the full molecular weight distribution (MWD) for this resin was not reported, we choose a realistic MWD for highly polydisperse materials and demonstrate that this distribution is consistent with both the reported molecular weight averages and linear rheological measurements. For non-linear flow modeling, we approximate this polydisperse melt as a bimodal blend, since no polydisperse equivalent of the GLaMM mode is currently available. We note that no fundamental modification of our simulation algorithm would be required to take advantage of the results of a polydisperse tube model. In this data comparison we are able to correctly describe these steady state nucleation rate measurements against shear rate, using crystallization parameters that are consistent with estimates provided by literature measurements.

Molecular relaxation times
The relaxation times of iPP T30G can be determined by the tube model and confirmed against linear oscillatory shear measurements. Taking the tube model parameters for iPP from the literature [30] and shifting τ e to the experimental temperature of 140 • C using the time-temperature superposition parameters for iPP T30G [31], gives τ e = 90ns and M e = 4.4kg/mol. Here τ e is the Rouse time of an entanglement segment and M e is the molecular weight between entanglements. The mass of a Kuhn step of iPP is 187.8g/mol [32], meaning the number of Kuhn steps per entanglement segment is N e ≈ 25. From these parameters the chain Rouse time can be calculated using τ R = (M/M e ) 2 τ e [33].
The form of the molecular weight distribution (MWD) for iPP T30G is essential to our modeling as small amounts of high molecular weight (HMW) material can control FIC (See, for example, ref [6]). Unfortunately, the full MWD of iPP T30G was not reported, with only the weight and number average molecular weights being given [25]. Furthermore, standard distribution functions, such as the log-normal and generalized exponential distributions, typically under-estimate the HMW tail of highly polydisperse melts [34]. Fortunately, linear rheological measurements are highly sensitive to this HMW tail at low frequencies and these data provide a clear upper-bound on the amount of HMW material. Thus we select a realistic form for the full MWD and show that this is consistent with both the reported complex viscosity and molecular weight measurements (M w = 481kg/mol and M n = 75kg/mol). We describe the MWD using a bimodal log-normal distribution, since the augmented tail of this distribution accounts more accurately for the HMW tail of highly polydisperse melts than standard monomodal distributions. This takes the form where φ h is volume fraction of the HMW peak, M wl/h and M nl/h are the weight and number average molecular weights for the low and high molecular weight chains, respectively, and w LN is the monomodal log-normal distribution. Taking φ h = 0.01%, M wh = 3 × 10 5 kg/mol, M nh = 2 × 10 4 kg/mol, M wl = 451kg/mol and M nl = 75kg/mol gives the correct overall weight and number average molecular weights and correctly predicts the complex viscosity measurements. These values were chosen as they correspond to the longest HMW tail that correctly predicts the complex viscosity measurements . Fig 4(a) shows the molecular weight distribution and fig 4(b) compares the predicted complex viscosity resulting from this distribution with measurements on iPP T30G. Here the linear oscillatory shear response was computed using a polydisperse tube model for linear response, which combines the tube escape formula of Likhtman and McLeish [35] with the Rubinstein and Colby algorithm for constraint release [36] and is implemented within our in-house software for molecular rheology, RepTate [41]. Using a monomodal log-normal distribution, with φ l = 0, produces slightly inferior agreement with the complex viscosity measurements and cannot account for the nucleation measurements. The shaded area is the region of the molecular weight distribution defined as the high molecular weight tail in our modelling. (b) Linear rheological measurements on iPP T30G [25] along with the predictions of polydisperse tube model using the molecular weight distribution discussed above and with G e = 1.1 × 10 6 Pa).
As discussed above, because there is no suitable polydisperse tube model for non-linear flow, we approximate iPP T30G as a bimodal blend when modelling its non-linear response .   Fig 4(a) shows that there is no clear molecular weight at which to define the division between the two blend species. We choose the high molecular weight fraction to be all material with a molecular weight above M * = 4200kg/mol, giving a HMW volume fraction of 1%. Taking a lower value of M * produces a Rouse time that is too short to account for the flow nucleation measurements, but we note that larger values of M * may be equally valid choices. We note also that this fraction contains discernible contributions from both of the monomodal terms in eqn (A1). We compute the expectation value of the Rouse time of the high molecular weight tail via,

Crystallization parameters
To model iPP T30G we require its quiescent critical nucleus size and barrier height, at the experimental temperature. An estimate for the number of Kuhn steps in a critical nucleus can be found from literature data. The critical nucleus diameter has been shown to be of the order of the lamella thickness by transmission electron microscopy [37,38] for iPP and confirmed for a different polymer by atomic force microscopy [39]. For iPP the lamella thickness is ∼ 10nm at 140 • C. Using the iPP crystal density of 0.9g/cm 3 [40] and the Kuhn step mass for iPP [32] and assuming a spherical critical nucleus gives the number of Kuhn steps in a critical nucleus as ∼ 1000.
The barrier height can be estimated from the nucleation data of Coccorullo et. al [25] (see fig 6). Projecting these data to zero shear rate gives a quiescent homogeneous nucleation rate of 3.4 × 10 −11 /µm 3 /sec. This can be converted to a rate per Kuhn step using the melt density of iPP (0.85g/cm 3 ) and the iPP Kuhn step mass. Thus the density of Kuhn steps is ρ K = 2.7 × 10 9 /µm 3 so the quiescent nucleation rate isṄ 0 = 1.2 × 10 −20 /sec per Kuhn step.
Taking the crystallization timescale of a Kuhn step τ 0 to be of the order of the Kuhn step relaxation time τ K = τ e /N 2 e = 0.144ns, means thatṄ 0 τ 0 ∼ 10 −30 . Clearly nucleation in this system is extremely rare and such a separation of timescales cannot feasibly be simulated even by our highly efficient kMC algorithm. Instead we systematically project our results to larger energy barriers to allow a comparison with these data.

Projection to large nucleation barriers
To perform this projection we simulated the average nucleation time over a range of quiescent nucleation barrier heights, varying from 10.6 to 26.3k B T, with the critical nucleus size held fixed at 1081 Kuhn steps. Nucleus alignment was neglected for all of these simulations (α = 0, throughout). The steady-state configuration for both blend components for various shear rates was calculated using a recent generalisation of the GLaMM model to bimodal blends [28]. The dimensionless steady-state nucleation timeτ N = τ N /τ 0 against shear ratė γ was simulated using the GLaMM model predictions, in steady state, for a range of shear rates and molecular weights. The results are independent of molecular weight when the shear rate is expressed as a Rouse Weissenberg numberγτ R , but do depend on the height of the quiescent nucleation barrier, ∆f * , as shown in figure 5(a). The simulated nucleation times vary exponentially with shear rate and can be described by the expressioñ barriers was used only to confirm the projection and was not involved in the fitting. This procedure allows the simulation results to be projected to nucleation barriers that are too high to be simulated.

Data comparison
We can use our projected results to model the data of Coccorullo et. al [25]. The value of the projected zero shear nucleation rate from these data, combined with the projection  In summary, our model correctly predicts the exponential dependence of nucleation rate on shear rate as observed by Coccorullo et al. [25]. Furthermore, with a characterization of the full molecular weight distribution and extrapolation of our simulation results to high nucleation barriers, we can make a direct quantitative comparison with these data. We use crystallization parameters that are estimated from the literature but whose final values are determined in response to the nucleation data. This results in very good agreement with directly measured nucleation rates under flow using parameters that are consistent with literature estimates.