Flow-enhanced nucleation (FEN) of n-pentacontahectane (C150) under biaxial extensional flows of varying strain rate ratios is studied using nonequilibrium molecular dynamics simulation. The nucleation rates thus calculated are used to test previously published FEN models based on invariants of the conformation tensor of Kuhn segments and the extra stress tensor. Models based on the conformation tensor provide a more accurate description of FEN observed in biaxial flow simulations than those based on the extra stress tensor. In addition, the formation of nematic domains previously reported to be stabilized by shear or extensional flow is absent in equibiaxial flows. However, such domains do form in non-equibiaxial flows, and nucleation occurs in these domains preferentially. The shape and orientation of nuclei formed under biaxial flows of various strengths and strain rate ratios are also reported.

## I. INTRODUCTION

Semicrystalline polymers constitute the largest fraction of industrial plastics. They are used for fabrication of fibers, films, packaging, and composites.^{1} Production processes for many of these products, such as injection molding, blown-film extrusion, and heat-sealing, involve complex flows coupled with solidification from the melt. The effect of flow on the final structure and properties, especially crystallinity, is complicated. Crystallization involves rearrangement of molecules into ordered regions; flow aids in this ordering through reorientation of polymer chain segments. The amount of stretching and orientation that the flow introduces dramatically affects the crystallization rate and consequent morphology.^{2} The evidence for rate enhancement in blown-film processes^{3} and changes in morphology in injection molding^{4} exemplify this effect of flow, commonly referred to as Flow-Induced Crystallization (FIC). The earliest stage of FIC is Flow-Enhanced Nucleation (FEN), which involves small molecular clusters that arise suddenly in the solution or melt^{5} but plays an important role in determining the final morphology of a material.

Many experimental studies have probed the phenomena of FIC and FEN.^{6–13} While experiments have provided valuable insights into these phenomena, their observations are typically indirect. By contrast, computational techniques, such as molecular dynamics simulations, are well suited to resolve the fine spatiotemporal scale of a nucleation event.^{14–22} However, to date, only the effects on nucleation due to shear or uniaxial extensional flow fields have been studied. Biaxial extensional flow (BEF) fields are industrially relevant for processes such as blown-film extrusion (BFE), but their effects on polymer crystallization have received only sparse attention.^{23–30} Bullwinkle *et al.*^{29} studied BFE using small angle light scattering. They modeled FIC but did not specifically study the effects of biaxial flow. Van Drongelen *et al.*^{30} measured crystallinity and molecular orientation in BFE using wide-angle x-ray diffraction. They found that the process conditions that lead to high molecular orientation also lead to higher crystallinity. Zhang *et al.*^{26} and Chen *et al.*^{25} studied experimentally the effect of transient biaxial stretching on the crystallization of polyethylene terephthalate (PET) and natural rubber, respectively. Both reports highlighted the importance of stretch ratio and orientation of amorphous chains. Nie *et al.*^{27} modeled nucleation kinetics under biaxial stretching using nonequilibrium molecular dynamics (NEMD) simulations and a coarse-grained model of polyvinyl alcohol; they applied a start-up flow followed by quenching below the melting point with cessation of flow and showed that the induction time for nucleation correlates with the average local orientation parameter at the final state of strain. However, the simulation protocol employed there convolves the effects of nucleation with relaxation after cessation of flow.

Nicholson and Rutledge^{21} examined nucleation rates for a united atom model of n-pentacontane (C150) in steady state shear or uniaxial extensional flow using NEMD simulations with lattice reduction techniques for generalized homogeneous flows. They simulated a range of strain rates and employed an iso-Weissenberg number quench to reduce the effect of the quench on chain conformation. They used these results to test several nucleation rate models proposed in the literature. Each model was classified according to the nature of the independent variable in the model: (a) strain rate, (b) stress, (c) conformation, or (d) free energy. Models based on measures of stress or conformation were found to perform significantly better than the others; in some cases, the two were indistinguishable. They argued that the similarity in performance is a consequence of the stress-optical rule and arises from coaxiality of stress and conformation tensors in shear and uniaxial extensional flow.

In this work, we report the effect of steady state biaxial flow on the orientation of C150 chain segments in the melt state, and in turn on nucleation kinetics, using the methods of Nicholson and Rutledge.^{21} Equibiaxial and non-equibiaxial flows are characterized by different ratios of the strain rate in orthogonal directions within the plane of biaxial extension. The orientation of chain segments under these flows is quantified using the conformation tensor and nematic order parameters. Nucleation kinetics under BEF are compared and contrasted with those published previously for uniaxial and shear flows. A model in which FEN is correlated with the P_{2} order parameter tensor for Kuhn segments is found to work well under all of the flow conditions studied.

## II. METHODS

### A. Flow simulation

The simulated system consisted of a melt of 100 C150 chains in which each CH_{2} or CH_{3} unit was modeled as a united atom (UA). The) force field of Paul, Yoon, and Smith (PYS)^{31} with modifications by Waheed *et al.*^{32,33} was used for the simulations; it has been shown previously to provide a realistic description of polyethylene melts and crystal nucleation.^{34–37}

^{38}was used to perform all the simulations. The USER-UEF

^{39}package allows for simulation of steady-state biaxial extensional flows with lattice reduction and remapping of compatible boundary conditions. The package employs the SLLOD equations of motion with boundary conditions as proposed by Dobson

^{40}and modified by Hunt.

^{41}A rate-of-deformation tensor as shown in the following equation was imposed in the simulation for biaxial flow:

*ξ*≤ 1. Evidently,

*x*and

*y*form the directions of extension and

*z*is the direction of compression.

*ξ*equal to 1 or 0 corresponds to planar extensional flow in the

*x*or

*y*direction, respectively; ξ = 0.5 corresponds to equibiaxial flow in the

*xy*plane. For convenience, we define a strain rate ratio Λ in the plane of extension as

*ξ*to the range 0.5 ≤

*ξ*< 1. Under these conditions, Λ = 1 corresponds to equibiaxial flow, Λ > 1 corresponds to biaxial flow with the major axis of extension in the

*x*direction, and planar extensional flow in the

*x*direction is obtained in the limit that Λ becomes very large. An effective strain rate ($\gamma \u0307$) was computed from the deformation rate tensor according to the following equation:

**D**= ∇

**u**+ ∇

**u**

^{T}is the strain rate tensor and $J2(D)=trD2$ is its second invariant. This deformation rate tensor is volume-conserving. However, to accommodate density changes associated with crystallization, the simulations were implemented with a constant stress boundary condition in the compression direction; the average stress in that direction was controlled to be 1 atm. In the absence of crystallization, the change in volume at a constant

*z*-component of stress averages to zero. For all simulations, the LAMMPS implementation of the Nosé–Hoover thermostat and barostat was used, with time constants of 0.4 and 4 ps, respectively. The equations of motion were integrated using the reversible reference system propagator algorithm (rRESPA)

^{42}with a 2 fs time step for bonded interactions and a 4 fs time step for nonbonded interactions. More details on the simulation procedure can be found in the work of Nicholson and Rutledge.

^{21}

For the simulation procedure, a steady state biaxial flow was imposed in the melt state at 440 K. This temperature is above the experimental melting point of 396.4 K for C150.^{43} The PYS force field has been shown elsewhere to reproduce experimental melting points for n-alkanes within about ±2%.^{35,44,45} The melt state flow simulation was run for 90 ns to reach the steady state, as confirmed by fluctuation of the Kuhn segment order parameter (defined below) about a steady average. The simulation was run for an additional 750 ns to analyze the melt flow state. 25 snapshots were sampled at intervals of 30 ns for subsequent nucleation studies. This time interval is slightly greater than the Rouse time (*τ*_{R}) of 26 ns^{36} for C150 at 440 K, ensuring independent, uncorrelated starting configurations. These configurations were then quenched to a crystallization temperature of 280 K, corresponding to a 29% supercooling ($Tm\u2212T/Tm$), to initiate crystal nucleation trajectories. *τ*_{R} for C150 at 280 K is estimated to be 340 ns.^{36} The strain rate of the simulation after quenching was chosen such that the Rouse Weissenberg number, defined as $Wi=\gamma \u0307\tau R$, was the same before and after the quench (an “iso-Weissenberg number quench”). Example strain rates before and after quench are presented in the supplementary material, Table S1. The steady-state NEMD simulations were conducted for strain rate ratios corresponding to Λ = 1, 1.2, 1.5, 3, and 10. Out of these simulations, configurations from simulations having Λ = 1, 1.5, and 3 were used for nucleation simulations, as these were sufficient to draw the conclusions of this work; with further increases in Λ, the deformation approaches planar extensional flow, which is similar to uniaxial extensional flow, for which data already exist.^{21}

### B. Analysis methods

**R**

_{i}and a Kuhn segment vector

**R**

_{K,i}for each united atom

*i*in the system and their unit vectors

**d**

_{i}and

**d**

_{K,i}, respectively. A Kuhn segment is taken to have

*n*

_{K}= 12 carbon atoms,

^{46–48}with the vector

**R**

_{K}going from the first to the last carbon of the Kuhn segment. The chain ends are treated differently, as shown in the definition in Eq. (5),

**r**

_{i}is the Cartesian coordinate of site

*i*and n

_{mol}is the total number of UAs in the chain, i.e., n

_{mol}= 150.

#### 1. Crystalline order parameter, clustering, and nucleation rate

*P*

_{2,i}, by which crystalline UAs are distinguished from noncrystalline UAs,

_{i}is the local neighborhood of the

*i*th UA, comprising all UAs

*j*within a cutoff distance of

*r*

_{local}= 1.5σ

_{LJ}≈ 0.6 nm from the

*i*th UA. Any UA that had P

_{2,i}greater than a threshold value P

_{2,th}= 0.7 was considered crystalline. These parameter values are consistent with those used by Nicholson and Rutledge for crystallization in flowing systems.

^{21}Two crystalline UAs within a threshold distance r

_{th}= 1.3σ

_{LJ}≈ 0.52 nm of each other were considered to be a part of the same crystalline cluster. To calculate a nucleation rate, the size of the largest cluster was tracked throughout the simulation. The first-passage time of the largest cluster FPT(

*n*) (that is, the time after quench at which the largest cluster first reaches a cluster of size

*n*UAs) was calculated for each of the 25 nucleation runs. Both the mean, τ

_{FPT}(

*n*), and standard deviation, σ

_{FPT}(

*n*), of the first passage time (FPT) distribution for the 25 trajectories were calculated. For a Poisson process, the two functions should be equal; however, in simulations where the growth rate of the clusters after nucleation is non-stochastic and measurable on the nucleation time scale, τ

_{FPT}(

*n*) continues to increase with increasing

*n*, whereas σ

_{FPT}(

*n*) levels off for values of

*n*several times greater than the critical nucleus size. For this reason, induction time was estimated in this work from the plateau value σ

_{plateau}of the standard deviation of the FPT distribution at large

*n*, as proposed by Shneidman.

^{49}The nucleation rate was then simply calculated as $Is=(\sigma plateauV)\u22121$, where V = 420.54 nm

^{3}is the volume of the system (corresponding to a density of 0.83 g/cm

^{3}). As discussed later, at sufficiently large Wi, σ

_{FPT}(

*n*

_{max}) was found not to level off in accordance with classical nucleation. In such cases, the plateau value was estimated as $\sigma plateau=\sigma FPTn150<n<300$. To estimate the error in nucleation rate, a jackknife or “leave-one-out” approach was used. In this method, nucleation rates were calculated based on all possible combinations of 24 simulations out of the 25. The standard deviation of these nucleation rates was used as the measure of uncertainty for the mean nucleation rate.

#### 2. Cluster fractal dimension, shape anisotropy, and orientation

For supercritical crystalline clusters, we characterized the compactness of the cluster using its fractal dimension, D_{f}. For this purpose, we followed the box counting algorithm described by Gagnepain and Roques-Carmes^{50} and also employed by Nicholson and Rutledge.^{39} Briefly, starting with a cubic box of side length L that just encloses the cluster, the box was further divided into voxels of side length δ = L/2, L/3, etc. For a given δ, the number of voxels Q(δ) that intersect with any part of a UA in the cluster was counted. D_{f} was then determined from the slope of −log(Q(δ)) vs log(δ) in the range 1 nm < δ < σ_{LJ}, where σ_{LJ} is the Lennard-Jones radius of a UA. D_{f} was found to be not very sensitive to the range over which δ was averaged (the 1 nm cutoff is chosen for consistency). For a three-dimensional object (such as our nuclei), D_{f} = 3 for an object that is compact and completely fills space. D_{f} < 3 for clusters that are not compact and space-filling.

^{2}) of a supercritical crystalline cluster was computed from the eigenvalues (λ

_{1}> λ

_{2}> λ

_{3}) of its gyration tensor as follows:

^{2}= 1 for a very thin cylinder, and κ

^{2}= 0 for a perfect sphere.

For each simulation, D_{f}, κ^{2}, and *n*_{max}, the size of the largest cluster, were each averaged over five consecutive 3 ns windows, and then, <D_{f}> and <κ^{2}> were plotted against *n*_{max}. For each biaxial flow mode (Wi, Λ), the results for 25 simulations were binned, and the mean and standard deviation of each bin were plotted vs <*n*_{max} > for that bin.

**S**

_{cluster}; the director is then the eigenvector

**e**

_{1}associated with the largest eigenvalue of

**S**

_{cluster}. The vector

**e**

_{1}was calculated for the largest cluster the first time it crossed

*n*= 200. The choice of

*n*= 200 is somewhat arbitrary so long as it is large enough to ensure the clusters are supercritical; in any case, the results were not sensitive to

*n*for any reasonable choice.

#### 3. Melt state structure

**C**

_{K}), defined as

**R**

_{K,i}for all UA sites in the system. $RK,02=trRKRK0$ is the mean squared Kuhn segment length at equilibrium in the absence of flow. The conformation tensor can be decomposed as $CK=\lambda K2SK$, where λ

_{K}is a scalar stretch factor, $\lambda 2K=trCK$, and

**S**

_{K}is a tensorial orientation factor. The global nematic Kuhn segment order tensor,

**P**

_{2,K}, is then given by

**I**is the identity matrix. The second invariant of this nematic order tensor

*J*

_{2}(

**P**

_{2,K}) has been shown previously to provide a good measure for the effect of flow-enhanced nucleation in shear and uniaxial extension.

^{21}

*P*

_{2,K,i}for nematic UAs in the melt,

_{i}is the local neighborhood of the

*i*th UA as defined in Eq. (8).

^{22}who introduced a co-alignment parameter to gauge the degree to which Kuhn segments associated with sites

*j*in the local neighborhood of site

*i*are co-aligned with the nematic director of that neighborhood (the directrix). Analogous to the cluster orientation analysis, a local orientation tensor for Kuhn segments is defined as

_{i}is the local neighborhood of the

*i*th UA as defined in Eq. (8) and the same

*r*

_{local}. The director $d\u0304k,i$ for the neighborhood of the

*i*th UA is the eigenvector (

**e**

_{1,K}) of

**S**

_{K,i}associated with its largest eigenvalue. The co-alignment parameter is then 1 minus the variance of Kuhn segment vectors with respect to the director for the local neighborhood of the

*i*th UA, that is,

Nicholson and Rutledge^{22} showed that two-dimensional distributions of P_{2,K,i} and ν_{K,i} split into two distinct peaks above a certain Wi for both uniaxial and shear flows. The peak associated with the lower values of P_{2,K,i} and *ν*_{K,i} was dubbed “flow-oriented,” whereas the peak associated with higher values of P_{2,K,i} and *ν*_{K,i} was dubbed “nematic.” Nematic UAs in the melt are those belonging to the nematic peak in this two-dimensional distribution.

*t*originated. Following Nicholson and Rutledge,

^{22}we define a memory function

*Z*

^{L}(

*t*),

_{i,0}is a measure of order for the

*i*th UA in the melt at time zero. For example, L

_{i,0}is some function of P

_{2,K,i,0}and ν

_{K,i,0}used to identify UAs within nematic domains in the melt at the time of quench (t = 0); in this work, we use L

_{i,0}= ν

_{K,i,0}, as explained later. $Z\nu K(t)$ is then the average of ν

_{Κ,i,0}for all UAs that are crystalline at time t. Again, P

_{2,th}= 0.7 is used to distinguish sites that are crystalline. Basically, $Z\nu K(t)$ averages the value of ν

_{K}at time 0 for those sites that are crystalline at time

*t*. A function $W\nu K(n)$ can then be defined based on $Z\nu K(t)$ as

*x*

_{i}is 1 if the

*i*th UA is crystalline and 0 otherwise. Thus, $W\nu K(n)$ is the average value of $Z\nu K(t)$ for all simulation snapshots (i.e., all times) in which the total number of crystalline UAs is

*n*, or crystallinity χ =

*n*/

*N*. In the limit that χ goes to 1, $W\nu K(n)$ necessarily converges to $\nu K,i,0i$, the system average at the time of quench,

*t*= 0. However, if nucleation occurs primarily within nematic domains (i.e., domains with higher value of ν

_{K}), $W\nu K(n)$ takes values of ν

_{K}higher than the system average for small

*n*, or vice versa if nucleation occurs primarily in flow-oriented domains. Thus, this definition of $W\nu K(n)$ readily shows whether nucleation occurs preferentially in nematic domains or flow-oriented domains. More details on this function can be found in the work of Nicholson and Rutledge.

^{22}

## III. RESULTS AND DISCUSSION

### A. Kuhn segment orientation in the melt

The imposition of flow induces orientation of the polymer chain segments to varying degrees in the melt state at 440 K. The degree of orientation was characterized by the second invariant of the global nematic Kuhn segment order tensor, *J*_{2}(**P**_{2,K}), which can range between 0 and 1.5. This second invariant is plotted vs Wi in Fig. 1 for the different flow conditions. The uniaxial extensional flow imposes the highest degree of orientation among all flows studied. *J*_{2}(**P**_{2,K}) is approximately zero at low Wi < 1, indicative of isotropic melts at small strain rates. *J*_{2}(**P**_{2,K}) begins to increase around Wi = 1, indicative of oriented melts, and then to level off above Wi = 10, as orientation becomes saturated. The same qualitative behavior is observed for the biaxial flows, but with *J*_{2}(**P**_{2,K}) being lower for the lower strain rate ratio in the range of Λ = 10 to 1, with the last corresponding to equibiaxial flow. The shear flow curve overlaps approximately the curve for biaxial flow at Λ = 1.2. These behaviors show that orientation can vary over a wide range under biaxial flow conditions, from orientations weaker than in shear flow for equibiaxial flow to orientations approaching that in uniaxial extension for strongly non-equibiaxial flows.

### B. Nucleation rate

To study nucleation kinetics, 25 independent configurations from the melt flow state were quenched to a crystallization temperature T_{c} = 280 K at constant Wi (Table S1 in the supplementary material shows typical strain rates before and after the quench). Figure 2 shows the size of the largest crystalline cluster as a function of time after quench for a typical trajectory. The first passage time (FPT) of the largest cluster in the system was tracked for each of the 25 nucleation simulations to create a first-passage time distribution, FPTD. Figure 3 shows the mean and standard deviation, τ_{FPT} and σ_{FPT}, of the FPT distribution as a function of the size, *n*, of the largest cluster for both small and large Wi. At small Wi [Fig. 3(a)], σ_{FPT} levels off at the end of the stochastic nucleation process, whereas τ_{FPT} continues to increase as growth processes become important for supercritical clusters. For sufficiently large Wi [Fig. 3(b)], neither τ_{FPT} nor σ_{FPT} levels off, indicative of a breakdown of classical nucleation theory, which is consistent with previously reported observations.^{22} For this reason, σ_{FPT} was evaluated for a fixed range of cluster sizes that was deemed to be only slightly supercritical (150 < *n* < 300), as described in Sec. II. This range was also chosen for consistency with prior work.^{21}

Figure 4 shows the nucleation rate as a function of Wi for different flow conditions. The nucleation rate, *I*_{S}, increases with increasing Wi for all flows, indicative of FEN. At low Wi < 1, the increase in the nucleation rate is approximately independent of the flow conditions, whereas for Wi > 1, the nucleation rate is strongly dependent on the flow conditions. As was the case for *J*_{2}(P_{2,K}) in Fig. 1, *I*_{S} decreases with the decreasing strain rate ratio between 3 and 1 in BEF at any given Wi, from the highest values for strongly non-equibiaxial flow to the lowest values for equibiaxial flow. For the case of equibiaxial flow, Λ = 1, the nucleation rate appears to level off around Wi = 10. This behavior is a preliminary indication that nucleation rate correlates with chain segment orientation. This correlation is explored quantitatively in the discussion of FEN models that follows.

### C. Comparison of FEN models

Figure 5 shows the correlation between the second invariants of the extra stress tensor **τ** = **σ** − *P*_{ext}**I**, where **σ** is the stress tensor and P_{ext} is the externally applied pressure (1 atm in our case), for shear, uniaxial, and biaxial flows over a wide range of Wi. The superposition of data for shear and uniaxial flow is indicative of the stress-optical rule and the coaxiality mentioned in the Introduction. Also shown in Fig. 5 are the dependence of these second invariants for the more general case of biaxial extensional flows for three different strain rate ratios. The biaxial flow data do not superpose. At low Wi (where the second invariants are small), the biaxial data are shifted upward relative to the shear and uniaxial data, indicative of elevated extra stresses in biaxial flows. At high Wi, the results for non-equibiaxial flows parallel the uniaxial flow case, whereas the equibiaxial data deviate toward still higher levels of extra stress. Evidently, significant deviations from the stress-optical rule are manifested in biaxial flows at large Wi.

^{21}using nucleation rates obtained by NEMD simulations of C150 in the more general case of biaxial flows. For conformation, we choose a model based on the second invariant of the

**P**

_{2,K}tensor,

_{K}is a model parameter. The quiescent nucleation rate I

_{S,q}can be treated as a second model parameter, or it can be determined independently by simulation of nucleation in quiescent systems, as done by Yi

*et al*.

^{36}Here, we treat I

_{S,q}as a model parameter, but the value obtained is of the same order of magnitude as the independently determined value. For stress, we choose a model based on the extra stress tensor,

Figure 6(a) shows the nucleation rate for shear, uniaxial, and biaxial flow conditions at three different strain rate ratios, plotted as a function of J_{2}(**P**_{2,K}). Also shown is the conformation-based model for nucleation rate fitted to all the data points shown in the plot. The model parameters found in this study, I_{s,q} = 32.46 × 10^{30} m^{−3} s^{−1} and θ_{K} = 6.01, are similar to the values reported earlier based on uniaxial and shear data alone (I_{s,q} = 38.5 × 10^{30} m^{−3} s^{−1} and θ_{K} = 5.79).^{21} The model based on the nematic order tensor, Eq. (18), provides a good description of all the data from NEMD simulations for shear, uniaxial, and biaxial flow conditions.

Figure 6(b) shows the same data for the nucleation rate plotted as a function of J_{2}(**τ**). Clearly, the data do not collapse onto a single curve, and no single parameterization can describe all the data. The model based on the extra stress tensor, Eq. (19), using parameters reported by Nicholson and Rutledge,^{21} is also shown. The nucleation rates for the three biaxial flows follow a similar trend as the model at low Wi but deviate at high Wi, indicating a reduced sensitivity of the nucleation rate to extra stress at large applied stresses. This deviation reveals that the extra stress tensor is not as robust as the Kuhn segment order parameter to correlate FEN.

In their experimental study of strain-induced crystallization of natural rubber, Chen *et al.*^{25} observed that crystallization did not occur under biaxial stretching conditions for stretch ratios less than 1.6. This suppression of crystallization was attributed to frustration in the alignment of amorphous segments under nearly equibiaxial stretch, which is consistent with the low levels of Kuhn segment orientational order observed under equibiaxial flow conditions for all strain rates (Wi) in the current work (Fig. 1). The low levels of orientational order achieved, despite the application of relatively high Wi (and correspondingly high applied stress, as shown in Fig. S2 of the supplementary material) seems to be indicative of this frustration. As the strain rate ratio increases, a dominant direction of extension emerges and the strain rate required to realize a given amount of orientation decreases.

### D. Cluster shape

The post-critical nuclei can be observed to understand the effect of flow on the shape and morphology of these nuclei. However, calculating the shape metrics from a single snapshot can be misleading since short-lived appendages appear and disappear with some frequency. Figure 7 shows examples of nuclei showing all UAs appearing as part of the largest cluster in a 3 ns window. The UAs that have low persistence are more transparent and yellow. Within 3 ns, these nuclei can fluctuate about the same shape [Figs. 7(a) and 7(b)], have short-lived appendages [Figs. 7(c) and 7(d)], or consider two nuclei that are close to each other [Figs. 7(e) and 7(f)]. Thus, averaging shape metrics (as detailed in Sec. II) gives us a better estimate of the persistent part of the nuclei. More details on the averaging and rendering of Fig. 7 are given in the supplementary material.

Fractal dimensions and relative shape anisotropy are plotted as functions of cluster size in Figs. 8 and 9, respectively. The fractal dimension values are all above 2.5, which for a three-dimensional object means that they are largely space filling. The shape anisotropy values are also below 0.5, with most of them being in the 0.2–0.3 range. Thus, the nuclei are mostly spherical and not cylindrical or rod-like. For a given Wi, increasing the size of nucleus decreases the fractal dimension, but the anisotropy is largely unaffected. Thus, the larger nuclei are not as completely filled as the smaller ones and become slightly more spherical as they grow. On the other hand, increasing Wi decreases the fractal dimension of nuclei and increases their shape anisotropy. Thus, slower flows allow nuclei to be more completely filled and more spherical in shape. Slower flows lead to slower formation of nuclei, which gives polymer chains more time to organize themselves in more filled and spherical nuclei. Comparing between Λ = 1 and Λ = 3, there are no significant differences in shape metrics. Thus, we can say that the type of biaxial flow does not affect the shape of nuclei for the flow conditions examined in this work.

### E. Cluster orientation

Flow biases the orientation of Kuhn segments in the melt, so it may be expected that the orientation of crystallites is similarly biased. The crystallite orientation can be conveniently characterized by the direction cosines of the director of **S**_{cluster} with respect to the axes of flow. Figure 10(a) shows the relevant section of a triangle plot of the square of the direction cosines of the cluster director, averaged over clusters of size (or n) approximately equal to 200 from 25 trajectories. To show the variation in orientations of these 25 nuclei, the director is projected onto the plane of extension (xy plane). The angle this projection makes with the x axis is called θ_{proj} and the mean and standard deviation of θ_{proj} is shown in Fig. 10(b). The data point for crystallite orientation in the absence of flow (quiescent crystallization) lies close to the centroid of the triangle in Fig. 10(a), showing that crystal stems have no orientation preference in this case. With increasing Wi, all data points for biaxial flows are shifted toward the xy plane. The data points for equibiaxial flow (in blue) are all approximately equidistant from the x and y vertices, but with a very broad distribution of values in the angle θ_{proj}. Figure 10(b) confirms that there is a large variation in orientation in the xy plane for the case of Λ = 1. Taken together, these observations indicate that the clusters are randomly oriented in the xy plane for all Wi when Λ = 1. For strain rate ratios Λ > 1, the crystallites orientation shifts toward the x axis with increasing Wi, and the distribution of orientation becomes rather narrow. Thus, crystallites in strong biaxial flows with Λ > 1 are oriented along the x axis (i.e., the axis of strongest extension). Some example angle distributions are shown in the supplementary material, Fig. S3.

### F. Nematic ordering in melt phase

Following the method of Nicholson and Rutledge,^{22} the formation of nematic domains in the melts under biaxial flow is examined using two-dimensional distributions of the populations of UAs as functions of the Kuhn segment order parameters ν_{K} and P_{2,K}. For Λ > 1, these distributions exhibit two peaks above a critical Wi. This behavior is shown in Fig. 11.

Although the individual order parameters are clearly correlated, projections of the distributions onto the individual order parameter planes reveal that the two distinct peaks are more evident as a function of ν_{K} than P_{2,K}. Thus, Fig. 12 shows projections of the UA populations on the ν_{K} plane for a range of Wi at three different strain rate ratios, Λ = 1.0, 1.5, and 3.0. Figure 12(a) and the first column of Fig. 11 show that equibiaxial flow does not give rise to the nematic order peak, and only a shoulder can be seen at the highest Wi. On the other hand, higher strain rate ratios [Figs. 12(b) and 12(c) and the second and third columns of Fig. 11] show bimodal distributions. The lower valued peak is designated the “flow-oriented” population, and the higher valued peak is designated the “nematic” population.

Figure 13 shows the function $W\nu K(n)$ as a function of *n* for several combinations of Wi and Λ. In all cases, one observes $W\nu K(n)$ > $\nu K,i,0i$ for all *n*, with $W\nu K(n)$ approaching$\nu K,i,0i$ in the limit that *n* approaches *N*, the total number of UAs in the system. This observation confirms that nucleation occurs preferentially in nematic domains as measured by ν_{K}.

Nie *et al.*^{27} performed NEMD simulations of a polyvinyl alcohol (PVA) undergoing biaxial stretch; they also reported a spatial correlation between the more ordered sites in the melt and the sites found in crystal nuclei. Using all-atom MD simulations of polyethylene (C1000) under shear flow, Tang *et al.*^{51} concluded that FEN is a two-step process in which nucleation occurs in regions of higher density. The higher density of these regions was attributed to their being more ordered. Thus, similar to this work, they observe nucleation occurring preferentially from more ordered regions.

Figure 14 shows a phase diagram based on whether a nematic phase was observed in the simulations at 440 K. Each of the biaxial flows (except Λ = 1) shows a minimum Wi for the formation of the nematic phase. Wi at which nematic domains emerge decreases as the strain rate ratio increases, suggesting that such nematic domains are stabilized by more anisotropic, or non-equibiaxial, flow conditions. This idea is supported by looking at extremes of anisotropy in flow. Equibiaxial flow, which is the most isotropic, fails to show nematic phase formation even at the highest Wi studied. On the other hand, uniaxial flow shows nematic phase formation at Wi = 2.31,^{22} which is lower than for all biaxial flows studied here. Thus, flows with a higher strain rate ratio not only give a higher average Kuhn segment order but also are more efficient in the formation of nematic order domains within the melt.

## IV. CONCLUSION

We performed NEMD simulations to study FEN under biaxial flows. Varying the strain rate ratio (Λ), we observed a variety of orientations in the melt. These orientations range from low orientation in equibiaxial flows to high orientation in anisotropic biaxially flows, asymptotically approaching uniaxial flow for large values of the strain rate ratio. Thus, biaxial extensional flows (BEFs) provide a good candidate to test FEN models. It is found that a model based on the Kuhn segment conformation tensor best explains FEN for all flow conditions. The stress-optical rule fails for BEF, and consequently, the model based on the extra stress tensor does not fit the data. Segmental orientation at the Kuhn segment scale can be obtained by melt birefringence studies, allowing for experimental validation of this model.

An NEMD study allows us to study the post-critical nuclei in detail and observe orientation and shape. We show that flow greatly influences the direction of crystal stems. Equibiaxial flows introduce randomness in orientation within the plane of extension, whereas asymmetric flows tend to orient all nuclei in one direction. We also observe that nuclei are largely spherical even for strong and anisotropic flows. Thus, nuclei do not show evidence of change in shape up to the sizes studied here. Hence, any change in morphology observed in experiments must occur at larger length scales than the ones studied here.

We are also able to show nematic phase formation in melt for biaxial flows that are asymmetric enough. Equibiaxial flows behave differently and fail to show the formation of a nematic phase. This nematic phase is also shown to nucleate preferentially. Thus, we can say that FEN occurs through the formation of more ordered parts in the melt phase. Our observations agree with other FEN and FIC studies conducted in the literature. From these observations, one can say that biaxial flows (especially equibiaxial) behave differently from uniaxial or shear flows previously studied. They orient melt and crystal stems differently and may not form a nematic phase. Thus, these flows must be treated appropriately in FEN experiments, simulations, and models.

## SUPPLEMENTARY MATERIALS

Additional information and images are provided as the supplementary material in a separate file.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge support from Dow through the University Partnership Initiative. The authors also acknowledge the Lammot du Pont Chair in Chemical Engineering to G.C.R. for funding of this work.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Chinmay S. Gangal**: Formal analysis (lead); Investigation (lead); Software (lead). **Gregory C. Rutledge**: Conceptualization (lead); Supervision (lead).

## DATA AVAILABILITY

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

## REFERENCES

*Polymer Science: A Comprehensive Reference, 10 Volume Set*

*Crystallization Modalities in Polymer Melt Processing*

*Polymer Crystallization II. Advances in Polymer Science,*

*n*-alkane melt by molecular simulation

*n*-eicosane

*n*-eicosane melts

*n*-eicosane melts under steady shear and uniaxial extension

*n*-octane melts

*Polymer Physics*

*Physical Properties of Polymers Handbook*