Tension-Induced Nematic Phase Separation in Bidisperse Homopolymer Melts.

We use an analytical mean-field theory and all-atom molecular dynamics (MD) simulations to predict that external tension, together with the nematic coupling interactions, can drive phase separation of long chains from short ones in bidisperse homopolymer melts. The nematic coupling parameter α for polyethylene (PE) oligomers under applied tension is extracted from the MD simulations and used in the mean-field free energy to predict the phase boundary for bidisperse melts in which the longer chains are stretched by uniaxial tension. The predicted phase diagram is validated by direct MD simulations. We also show that extensional flow, and possibly even shear flow, may lead to nematic phase separation in molten PE oligomers, because the flow can impose a stronger tension on the longer chains than the short ones.

Introduction

Polymer blends and block copolymers can undergo macro- and microphase separations, respectively, driven by the incompatibility of different species, which is characterized by the Flory–Huggins χ parameter.[1] The value of χ depends on the mismatch in enthalpic interactions and local packing of different monomers.[2] When χ exceeds a critical value beyond which the excess free energy of mixing compensates the loss in translational entropy, inhomogeneous polymers demix.

For chemically similar polymers with negligible χ, orientational interactions can also lead to phase separation. For example, stiff polymers that exhibit a nematic phase can undergo macroscopic phase separation from more flexible ones that form an isotropic phase.[3] Mixing stiff and flexible chains dilutes the local nematic order, resulting in a smaller gain in orientational free energy, which in turn drives phase separation. In this case, the nematic interactions must be strong enough that one of the polymers can spontaneously undergo the isotropic-to-nematic (IN) phase transition.

Away from the IN phase boundary, the nematic interactions can still drive phase separation when external aligning fields are applied. Using a mean-field model, Olmsted and Milner predicted that nematic interactions and external strain fields that act only on the longer chains can drive phase separation in bidisperse homopolymer melts and gels.[4] The orientational energy of the strained system includes the Maier–Saupe interactions and the entropic orientational free energy that is calculated perturbatively as a function of nematic order parameters using the Landau theory. This phase separation arises form the competition between mixing entropy and gain in orientational free energy. The critical strain for phase separation decreases with increasing Maier–Saupe parameter α, which quantifies the interactions between monomer tangent vectors and the nematic field. The OM theory, however, is only strictly valid for systems with low applied strain and weak nematic order so that a large α is required for phase separation.

Here we combine an analytical theory and all-atom molecular dynamics (MD) simulations to demonstrate that external tension and nematic interactions can indeed induce phase separation in bidisperse polyethylene (PE) oligomers. PE is a rather flexible polymer and does not exhibit a spontaneous nematic phase because the nematic coupling parameter α in PE is weak under quiescent conditions. Even if α is too small on its own to drive a nematic transition, we show that phase separation can be observed when sufficient tension is applied to the chains. Increasing external tension can lead to increased anisotropy in chain conformation, which in turn increases the nematic coupling strength.

Because external fields such as flows can stretch and align chains to an extent that depends on their length, a tension-induced nematic phase separation may occur during polymer processing. While previous authors have developed theories for elastic stress-concentration coupling in sheared polymer solutions with negligible nematic interactions, these theories lead only to enhanced concentration fluctuations and not phase separation.[5,6] When sufficiently strong nematic interactions are present, however, flows may drive phase separation in polymer melts.[7,8] Using nonequilibrium MD (NEMD) simulations, we demonstrate, for the first time, that uniaxial extensional and shear flows may lead to phase segregation in bidisperse homopolymers. This flow-induced nematic phase separation could be critical for understanding flow-induced crystallization in polymers, including the role of long chains and the origin of critical flow rate and stress in the formation of precursors to shish-kebab structures.[9,10] Even noncrystallizable polymers might be strongly affected by flow-induced local segregation of polymers of differing molecular weight.

Results and Discussion

To estimate the nematic coupling parameter α for PE oligomers, we combine MD simulations with self-consistent field theory (SCFT).[11] We first perform NPT simulations at 550 K and 1 bar for 20 ns using the GROMACS package[12] and the OPLS-AA force field to generate equilibrated melts of C20, C40, and C80 alkanes. The total number of carbon atoms in each simulation is 16 000. Without applied tension, PE chains are isotropic. Using these isotropic chains, we compute the intrachain tangent–tangent correlation function ⟨t0·t⟩, where s is a monomer index. The persistence length measured in units of numbers of monomers Np = 4 is then obtained using ⟨t0·t⟩ = exp(−s/Np) (Figure ).

Figure 1

Tangent–tangent correlation function for isotropic C40. Curve: ⟨t0·t⟩ = exp(−s/4).

Starting with isotropic melts, we apply uniaxial tension to all chains at 550 K in the simulations by stretching the head and end monomers with forces of magnitude f in opposite directions along the x-axis (x̂) for 10 ns. This simulation time is greater than the Rouse relaxation time of the end-to-end vector of C80 (τC80 = 1.4 ns). Pressure coupling is only applied in ŷ and ẑ directions because dynamic scaling of the box dimension is not allowed in the tension direction in GROMACS. To quantify the induced order, we use a nematic order tensor Q = ⟨tt⟩ – δ/3, where t is the backbone tangent, δ is the Kronecker δ function, and the indices are i,j = (x, y, z). Because the alignment is along the x-axis, we simplify Q using a scalar order parameter q = ⟨P2(μ)⟩, where μ = t·x̂, and P2 is the second order Legendre polynomial. The value of q then equals 1.5λ, where λ is the largest eigenvalue of Q. After short equilibration periods, q fluctuates around the equilibrium value in our simulations. The value of q increases with increasing applied dimensionless tension βfa, where β = 1/KT, f is the applied force, and a is the monomer size (Figure a).

Figure 2

(a) Nematic order q vs time for C40 under various dimensionless tensions. (b) SCFT order q vs dimensionless coupling parameter βα for C40 with reduced tension βfa = 0.25 (solid curve). Fitting to MD (dotted lines). (c) Tangent distribution function for C40 with βfa = 0.25. MD results (symbol), SCFT (curve). (d) βα vs reduced tension for PE oligomers. Linear fits (solid lines). Inset: βα vs reduced work done by tension. Linear dashed line is to guide the eye.

Together with α, the applied tension f governs the nematic order of PE chains. By treating PE oligomers as wormlike chains, we write the single-chain Hamiltonian for a chain of length N as

Using this Hamiltonian, we write the propagator z(μ;S) asin which z0 normalizes the propagator with respect to a free chain. The value of z(μ;S) gives the statistical weight of a chain of length S, starting with backbone tangent orientation μ. The evolution of z(μ;s) follows a biased diffusion equation:in which ∇⊥2 is an angular Laplacian. The diffusion equation can be solved by expanding z(μ;s) in the eigenfunctions ψ of the right-hand side operator, and then applying the Legendre polynomial expansion to ψ.[11]

The backbone tangent distribution function is thus proportional to the product of two propagators, integrated over the whole chain contour N:By integrating P2(μ) over the distribution p(μ), we can self-consistently obtain the nematic order q as a function of α and f.

By fitting the calculated q(α, f) to the nematic order in the simulations, we obtain α(f) for each chain length N (Figure b). Using the fitted α(f), we also compute the tangent distribution function, which gives a more detailed description of chain orientation than does q alone. Our calculations agree well with the MD simulations (Figure c). In the absence of the nematic interactions, our model yields the force–extension relation of an isolated wormlike chain,[13] which underestimates the nematic order of the stretched molten chains. This is because the induced order is enhanced by the nematic interactions.

The value of α determined from the above procedure increases linearly with f for chains with different lengths (Figure d). The work done by tension fR (R is the mean end-to-end distance of chain along x̂) leads to a higher effective bending stiffness for the wormlike chain, which gives rise to an increased nematic coupling constant α.[14−16] We expect α to plateau eventually at large applied tension when the backbone becomes completely rodlike. By extrapolation in Figure d, we also find that the values of α for C20, C40, and C80 without applied tension are the same, about 1.1 kT. Such a small α is expected because PE oligomers do not exhibit a spontaneous nematic phase.

Unequal tension may be imposed on chains in a bidisperse melt so that different degrees of alignment qs and ql can be observed for the short and long chains, respectively. In this case, we write the mean nematic field qeff = ϕlql + ϕsqs, in which ϕ is the volume fraction of each type of chain, and ϕl + ϕs = 1. The presence of less ordered short chains lowers the mean nematic field relative to that for all long, highly ordered chains. Using fl and fs and the resulting nematic coupling parameters αl(fl) and αs(fs), we simultaneously solve the propagator equations (q is replaced by qeff in eq ) to obtain ql and qs for the long and short chains. Because the composition fluctuations are neglected, slight discrepancies between the mean-field theory and MD simulations may be observed for melts near the binodal (Figure ).

Figure 3

Tangent distributions p(μ) of C40 (blue) and C80 (red) in one-phase melts near the binodal. C80 volume fraction: (a) ϕl = 0.2 and (b) ϕl = 0.5. Theory (curve). Simulations (symbols). Reduced tension βfa = 0.11 applied only to C80.

With high enough coupling strength and applied tension on the long chains, the bidisperse chains can phase separate, driven by the excess nematic free energy upon mixing. To show this, we construct a mean-field free energy for the binary system asin which Zl and Zs are the single-chain partition functions for long and short chains, calculated using the propagator: Z = ∫–11 dμ ∫0 dsz(μ;s) z(μ;N – s), and assuming a chain is in an external field qeff. For a multichain system, the mean alignment of all chains gives rise to qeff. Thus, the third and fourth terms in eq double-count the gain in nematic energy per monomer, which is corrected by the last term in eq . The mixing free energy is concave up as a function of ϕl, and the blend forms two phases when sufficient tension is applied (Figure a).

Figure 4

(a) Mixing free energy vs composition for bidisperse C40/C80 with only C80 under tension. Common tangent construction (dashed line). (b) Phase diagram for bidisperse C40/C80. MD results (symbols). One-phase (blue). Two-phase (red). Two-phase with hexagonal packing of C80 (purple). (c) Snapshots for bidisperse C40/C80 (ϕl = 0.2). Only carbons on C80 shown for clarity. (d) Partial interchain pair correlation function for C80 (ϕl = 0.2).

From the mixing free energy, we obtain the binodal using a common tangent construction and the spinodal by solving ∂2βΔFmix/∂ϕl2 = 0 (Figure b). To validate the phase boundaries, we also perform 20 ns equilibration MD simulations for C40/C80 with different compositions and applied tensions. For systems with βfa = 0.44, the initial configurations are one-phase melts, obtained after 10 ns of NPT simulations at 550 K and 1 bar. For other systems, the initial configuration is the last configuration of the previous simulation to which the next higher uniaxial tension is applied. Snapshots for systems with ϕl = 0.2 are shown in Figure c (other snapshots are included in the SI). The mean-field phase boundaries agree well with MD simulations, in which phase separation, indicated by segregation of C80 chains, can be observed (in Figure c, for example) when f is large enough with all observations from the snapshots summarized in Figure b. For the system with ϕl = 0.6, external tension can even result in hexagonal order in the C80-rich domains (see the SI). The “tension-induced crystallization”, however, is beyond the scope of this paper.

The evidence of phase separation can be also found in the interchain partial pair correlation function g(r) of C80 carbon atoms, defined asin which N is the total number of carbon atoms on all C80 chains, N* = N – 80 is the slightly smaller number of interchain carbons, and the partial C80 density ρ* = N*/V where V is the system volume. The last 5 ns of simulation data of each run is used to generate g(r), given that the system needs to re-equilibrate after a new tension is applied. When the external tension is weak (for example, βfa = 0.11), the system is homogeneous with g(r) almost identical to the g(r) of a quiescent melt (Figure d). When phase separation occurs, the intensity of g(r) is enhanced, arising from the formation of C80-rich domains.

To demonstrate the flow-induced nematic phase separation, we perform NEMD simulations of bidisperse melts composed of 20% C80 and 80% C40 (ϕl = 0.2) in uniaxial extensional and shear flows at 550 K using LAMMPS.[17] The initial configurations are one-phase melts obtained using standard NPT simulations. To simulate chains in steady extensional flows, a generalized Kraynik–Reinelt boundary condition,[18,19] implemented by the Rutledge group in LAMMPS,[20] is applied together with the SLLOD equations of motion.[21] During our simulations, the simulation box deforms and is remapped to conserve image locations and avoid too close distances between particles and their periodic images. The initial simulation box is cubic (10.43 nm3) and comprises 32 000 carbon atoms. The rather large box is used for extensional flow simulations so that the box dimension is greater than the contour length of C80. To simulate chains in steady shear flows, the Lees–Edwards boundary condition is applied in the velocity gradient direction.[22] The initial dimensions of the shear simulation box (contains 16 000 carbons) are 10.4 nm in the flow direction (greater than the chain length of C80), and 7.4 nm in the gradient and vorticity directions.

We perform three independent simulations for 8 ns using different initial configurations for each flow condition. The simulation time is rather short so that the total computational cost is kept manageable. In our flow simulations, q for C80 and C40 start to fluctuate about steady values within about 2 ns. The nematic order of both C80 and C40 at steady state increases with increasing flow rate ϵ̇ (Figure a). Extensional flows, however, can orientate chains more effectively than can shear flows. The ratio of nematic order for C40 to that of C80 also increases with ϵ̇. This is because flow can stretch the short chains and enhance the nematic coupling of C40 to itself and to C80.

Figure 5

(a) Nematic order for C80 in extensional and shear flows (upper) and ratio of nematic order for C40 to C80 (lower). Weissenberg number Wi = ϵ̇τC80. (b) Partial interchain radial distribution function g(r) for C80 under different flow conditions. g(r) of homogeneous (quiescent) and weakly phase-separated C80s under applied tension (βfa = 0.19) shown as references (ϕl = 0.2). Inset: snapshots of bidisperse C40/C80 in quiescent condition and extensional and shear flows. Only carbons on C80 are shown. (c) Nnear (as defined in the text) vs tension and flow rate in shear and extensional flows. Nnearmelt obtained from quiescent melt. One-phase melts (red), two-phase melts (green), and melts near the binodal (blue).

At sufficiently high flow rates, we observe weak segregation of C80 chains in NEMD simulations (Figure b inset). To quantify this segregation of C80 from C40 chains, we compute g(r) for carbons on C80 chains from the last 5 ns of simulation trajectories (Figure b). For each flow condition, g(r) is averaged over three independent simulations. Figure c shows that the peak intensities of g(r) for C80 in extensional flow are higher than that obtained from equilibrium MD simulations for a weakly phase-separated C40/C80 melt, in which βfa = 0.19 is applied to 20% C80 (see Figure c). The enhanced g(r) indicates that extensional flow may induce phase separation in the bidisperse C40/C80 melts. Shear flows, however, are less effective in driving the nematic phase separation than are extensional flows, given that the enhancements in g(r) are not as significant. The nematic order induced by shear flow is rather weak so that the excess nematic energy of mixing is low.

We also compute the number of nearest interchain neighbors for C80s using Nnear = ∫0 dr 4πr2ρ*g(r), where rmin is the location of the first minimum in g(r) (0.75 nm). In MD simulations, the value of Nnear increases with increasing applied tension (Figure c). Using Nnear that we obtained in MD simulations for homogeneous (βfa = 0.11) and phase-separated (βfa = 0.19) C40/C80 (ϕl = 0.2), we define the boundaries of one-phase and two-phase melts. We observe that the bidisperse melts may phase separate even in rather slow extensional flows. Chains in shear flows, however, at best barely reach the binodal.

For chains in either extensional or shear flows, the intensities of g(r) and the values of Nnear tend to plateau at intermediate flow rates and drop at high flow rates. This is different from the MD simulations in which g(r) and Nnear increase monotonically with tension applied to C80 only. In fast flows, both long and short chains may become quite stretched, thus reducing the driving force for them to phase separate from each other. In addition, convection, together with the low surface tension between the two chemically identical phases, likely prevents the formation of large domains. To observe clear nematic phase separations in slow flows, and especially in shear flows, we may need longer bidisperse chains with a greater contrast in molecular weight. Simulating such chains, however, is expensive because large systems are required.

We expect that the flow-induced nematic phase separation may be also detected using in situ neutron or X-ray scattering. Although the high flow rates required to drive phase separation in short alkanes (ϵ̇ > τR,l, where τR,l is the Rouse time of the long chain) are difficult to achieve in experiments, observing the nematic phase separation in entangled polymer samples should be feasible with flow rates in the range τR,l–1 < ϵ̇ < τd,s–1, where τd,s is the reptation time of the short chains. The value of ϵ̇ is bounded by τd,s–1 because fast flows may induce orientational order to the entanglement strands of the short chains, and in turn reduce the driving force for phase separation. The required flow rate also implies that well controlled molecular weight distribution in the sample is critical. A large contrast in molecular weight Ml2 ≫ 3Ms3/Me is required so that τd,s ≪ τR,l, where Ml, Ms, and Me are the molecular weights of the long and short chains, and the entanglement strand, respectively. For example, to observe the flow-induced nematic phase separation for a bidisperse PE sample, in which Ms is about 10 kg/mol, the long chain molecular weight Ml needs to be at least 1700 kg/mol (Me is about 1 kg/mol for PE[23]).

Nonetheless, the flow-induced nematic phase separation may play an essential role in polymer processing because the segregation of more ordered long chains should affect the crystallization behavior once the melt is cooled. When the effective tension imposed by flow on long chains is not sufficient, the ordered long chains remain dispersed, and may in turn facilitate the crystallization of surrounding short chains. This could lead to the formation of oriented “rice grain” crystallites.[10] When phase separation between long and short chains occurs prior to crystallization, the domains that are rich in aligned long chains may act as the precursors for “shish” formation. This may explain why shish forms when a critical stress is achieved.[10] Of course, these domains may be small due to convection so that the phase-separated long chains only facilitate the formation of an initial shish. Less ordered short chains, especially those in the interfacial regions, can join the growing shish easily because they are aligned along the shish axis by the nematic interactions and flow. As such, the final shish structures are not dominated by the long chains, as demonstrated experimentally using small-angle neutron scattering.[9]

Our current work may also shed some light on the formation of precursor state for polymer crystallization. For example, spinodal scattering can be observed in crystallizing polymer melts in shear flow,[24,25] indicating the formation of a phase-separated precursor state. Olmsted and co-workers constructed a phenomenological model and suggested that the phase separation is driven by the density-conformation coupling of polymers.[26] Shear flows only shift the phase boundaries to higher temperatures. Our predicted tension-induced nematic phase separation may provide another possible origin of the observed spinodal scattering. In addition, we expect the nematic coupling parameter α to increase as temperature decreases,[11] and may in turn lead to the formation of nematic domains (i.e., nucleation of the nematic order) at low temperatures before crystal nucleation occurs in the quiescent condition. In fact, recent united-atom MD simulations suggested that polymers may first form uniaxially aligned domains after quench, in which crystal nuclei grow.[27]

Overall, we expect the nematic interactions and the nematic phase separation discussed in this paper eventually to lead to better understanding of polymer crystallization, especially flow-induced crystallization. While our equilibrium theory cannot be directly applied to polymers in flows, we have provided strong evidence that flow-enhanced nematic coupling can drive local phase separation in homopolymer melts. This local phase separation may play an important role in partial segregation of chains by length during subsequent crystallization and may be important also in flowing polymers that do not crystallize. To predict the behaviors of polymers in realistic processing conditions, a dynamic theory which includes the chain conformation-dependent nematic interactions is necessary, but left to future work.