Predicting Nematic Phases of Semiflexible Polymers - American

Feb 27, 2015 - transition temperature TIN can be calculated. We apply our method to obtain α and TIN of a commonly studied semiflexible conjugated po...
1 downloads 0 Views 4MB Size
Download PDF
Article pubs.acs.org/Macromolecules

Predicting Nematic Phases of Semiflexible Polymers Wenlin Zhang,† Enrique D. Gomez,*,†,‡ and Scott T. Milner*,† †

Department of Chemical Engineering and ‡Materials Research Institute, The Pennsylvania State University, University Park, Pennsylvania 16802, United States ABSTRACT: The nematic coupling constant α, together with the chain stiffness κ, governs chain alignment and the isotropic-to-nematic (IN) transition temperature TIN for semiflexible polymers. We combine self-consistent field theory (SCFT) with atomistic molecular dynamics (MD) simulations of semiflexible chains under external tension in the isotropic phase to determine the nematic coupling constant α. Using α, we obtain the variational free energy of a multichain system, from which the IN transition temperature TIN can be calculated. We apply our method to obtain α and TIN of a commonly studied semiflexible conjugated polymer, poly(3-hexylthiophene) (P3HT). We predict TIN to be above the crystal melting temperature Tm for P3HT and to follow TIN(S) = 535K(1 − 1.64/S), in which S is the number of monomers.

■

INTRODUCTION Semiflexible polymers can exhibit nematic phases, resulting from the anisotropic shape of backbone segments. Semiflexible chains in a concentrated solution can be represented as a sequence of nearly rigid backbone segments dispersed in solvent. In the melt, flexible side chains act like “bound solvent”, tending to disperse the stiff backbones. Rigid rods undergo an isotropic to nematic transition as a function of concentration, as demonstrated by Onsager.1 The transition occurs when the volume fraction of rods ϕ is high enough that randomly placing and orienting the rods would lead to about one collision per rod. At high concentrations, the isotropic phase is no longer entropically viable; the rods tend to align with each other, resulting in a nematic phase. The critical volume fraction scales as the aspect ratio d/l, for rods of diameter d and length l. Likewise, semiflexible rods have been shown theoretically to undergo an isotropic to nematic transition as a function of concentration. For semiflexible chains, the critical volume fraction is governed by d/Lp, where Lp is the persistence length of the chain.2,3 The nematic phase is important for many semiflexible polymers. As liquid crystalline polymers, semiflexible polymers with nematic phases can be used in numerous applications, including displays, high strength fibers, microelectromechanical systems, and biomedical devices.4−7 Moreover, the existence of nematic phases enables better processing of functional semiflexible polymers. For example, crystallization from the nematic phase or other liquid crystal phases may promote the formation of large crystalline domains and ordered amorphous regions, which can in turn enhance the electrical properties of semiconducting polymers.8−11 The nematic coupling constant α quantifies the orientational coupling between backbone segments. Together with the chain stiffness κ, α governs the IN transition temperature TIN. For chains with stiffer backbones and stronger nematic coupling, the IN transition temperatures are higher. The value of α also © XXXX American Chemical Society

affects molecular ordering in mesophases and at interfaces. For example, Morse and Fredrickson predict that semiflexible chains tend to align parallel to the plane of the interface in the “stiff” regime (βκχ ≫ 1 where χ is the Flory−Huggins parameter).12 Such orientational ordering at the interface would be enhanced by nematic coupling between chain backbones. For polymers with accessible IN transitions, α can be estimated from the observed transition temperature. Olsen and co-workers estimated the nematic coupling constant from the IN transition temperature TIN for rod-like poly(2,5-di(2′ethylhexyloxy)-1,4-phenylenevinylene) (DEH-PPV).13 Ho et al. applied the same method to estimate the nematic coupling constants of two semiflexible polymers: poly(3-(2′-ethyl)hexylthiophene) (P3EHT) and poly(3-dodecylthiophene) (P3DDT).14 For many semiflexible chains, however, the IN transition is precluded by either crystallization or thermal degradation so that its location cannot be used for estimating α. In this paper, we report a method that combines selfconsistent field theory (SCFT) with atomistic molecular dynamics (MD) simulations of semiflexible chains under external tension in the isotropic phase to predict α from molecular structure. In MD simulations, we induce weak nematic order in semiflexible chains by applying uniaxial tension. For the same semiflexible chains, we can compute the nematic order parameter analytically for a given applied tension and α using SCFT. The value of α is determined by fitting the order parameter q(α) calculated by SCFT to the MD simulation results. Using the predicted α, we calculate the variational free energy for a multichain system. The mean-field free energy exhibits a first-order phase transition at the IN transition temperature TIN. Our method for predicting nematic phases is more convenient than directly observing the IN Received: January 3, 2015 Revised: February 23, 2015

A

DOI: 10.1021/acs.macromol.5b00013 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Using the aligning potential V(t), we write the Hamiltonian of an aligned semiflexible chain containing S monomers of length a as

transitions in simulations, for which equilibrium is challenging to achieve and the finite system size smears the transition. We apply our method to predict the nematic coupling constant α of poly(3-hexylthiophene) (P3HT), a commonly studied conjugated polymer. Using the predicted α, we obtained the mean-field free energy for P3HT, from which the IN transition temperature TIN(S) is calculated as a function of chain length S. We predict that TIN(S) lies above the experimental crystal melting temperature Tm(S) of P3HT as a function of chain length, indicating that the nematic phase is stable over a range of temperatures. Our predictions are consistent with experiments, in which a nematic phase has been reported for P3HT based on polarized optical microscopy (POM).14 Liquid crystal order has also been observed during solvent evaporation for P3HT in 1,2,4-trichlorobenzene using polarized Raman spectroscopy.15 The value of TIN is not available, however, because thermal degradation may preclude the IN transition in experiments.

H0 =

⎛ κ dt s ds⎜⎜ 2 d s ⎝

2

⎞ + V (ts)⎟⎟ ⎠

(2)

Here s is the monomer index, and ts is the unit tangent vector of monomer s. The bending modulus is denoted by κ, given by κ = kT Np, where Np is the persistence length measured in repeat units. The first term in eq 2 is the bending energy of a worm-like chain. We expand the uniaxial aligning potential V (eq 1) in Legendre polynomials Pn(cos θ). Because nematic ordering has azimuthal symmetry, we only consider the dependence of V on μ = cos θ, in which θ is the deflection angle between a backbone tangent and the nematic director (which coincides with the axis of stretching). The quadrupolar order parameter q in the tensor Q can be then written with the second-order Legendre polynomial P2:

■

METHODS A Single Semiflexible Chain under Tension. We determine the nematic coupling constant α from the tensioninduced chain ordering of semiflexible chains in the simulated isotropic phase. When tension is applied to a chain in a melt, it stretches along the direction of force, and its backbone segments tend to align in that direction. When tension in a given direction is applied to all chains in a melt, an average segmental alignment results. If the segments interact with a nematic coupling α, the segmental alignment is amplified, beyond what would be expected for a single chain under tension. By comparing the observed segmental alignment in MD simulations of chains under modest tension in the isotropic phase to predictions using SCFT, we can infer the nematic coupling. We perform MD simulations of the chains under tension because simulations at elevated temperatures equilibrate more rapidly. For atomistic MD simulations of even rather short semiflexible chains, it is challenging to carry out properly equilibrated simulations of the nematic phase directly. Instead, we observe chain alignment under tension to determine α, extrapolate its value to lower temperatures, and predict the nematic phase behavior using SCFT. To carry out this procedure for determining the nematic coupling constant α, we construct a self-consistent field theory for a semiflexible chain under external tension, immersed in a sea of neighboring chains under the same applied tension. The applied tension induces segmental alignment at second order, creating a mean aligning quadrupolar field Q. For a given α, we can determine the mean-field Q self-consistently using the single-chain model. The nematic coupling constant α is determined by fitting the calculated Q(α) to the observed Q in our MD simulation of many stretched chains. The aligning potential acting on a backbone tangent t of size a on a semiflexible chain with applied tension f is V (t) = −α t·Q·t − a f·t

∫0

S

1

q=

∫−1 dμ Π(μ)P2(μ)

(3)

in which Π(μ) is the probability distribution of backbone tangents. Using the order parameter q, we write V(μ) as V (μ) = −αqP2(μ) − fP1(μ)

(4)

To self-consistently compute the tangent distribution Π(μ) for a given α, we write a propagator Z(μ;S) to describe the statistical weight for a chain with length S, starting with tangent orientation μ: Z(μ; S) = AS

∫ dμ1 dμ2 ...dμS e−βH δ(μ1 − μ) 0

(5)

in which β = 1/kT. The factor A normalizes the propagator with respect to the free chain (V = 0). Our Z(μ;S) is analogous to the propagator Z(r,μ;S) used by other authors to describe semiflexible homopolymer blends or block copolymers, for which not only the tangent orientation μ is important but also the spatial position r.12,16,17 The tangent distribution Π(μ) is proportional to the product of two propagators, integrated over the whole chain length S: S

Π(μ) =

∫0 dn Z(μ; n)Z(μ; S − n) Z0

(6)

Here Z0 is the single chain model partition function, written as 1

Z0 =

∫−1 dμ Z(μ; S)

(7)

Using the propagator, we can determine the nematic order parameter q for a given nematic coupling constant α selfconsistently. Self-consistency is satisfied when the calculated q (eq 3) agrees with the input q in the single chain Hamiltonian H0 (eq 2). The propagator Z(μ;s) satisfies a diffusion equation, in which the backbone tangent vector ts executes a biased random walk on the unit sphere:

(1)

The mean quadrupolar field Q is a symmetric and traceless second-order tensor with the form Qij = ⟨titj⟩ − δij/3. With tension applied along x̂, the mean quadrupolar field can be written as Qxx = q, Qyy = −q/2, Qzz = −q/2, in which q is a quadrupolar order parameter.

∂Z(μ; s) = D∇⊥2 Z(μ; s) − VZ(μ; s) ∂s

(8)

The diffusion coefficient D is given by D = 1/(2βκ); tangent− tangent correlations decay more slowly for stiffer chains. Here B

DOI: 10.1021/acs.macromol.5b00013 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules ∇⊥2 is the angular Laplacian, in which only the θ dependent (μ dependent) part is considered. The propagator equation (eq 8) can be solved by expanding Z(μ;s) in the eigenfunctions ψα of the right-hand side operator in eq 8: N

1

Z(μ ; s ) =

∫−1

dμ′

∑ ψα(μ)ψα(μ′)e−ϵ s α

α=0

(9)

in which ϵα is the corresponding eigenvalue of ψα. The eigenfunctions ψα are in turn expanded in Legendre polynomials, which are the eigenfunctions of the angular Laplacian ∇⊥2.18 To obtain Z(μ;s), the eigenfunction expansion (eq 9) and the Legendre polynomial expansion of ψα are truncated at finite order N. Details for solving eq 8 are given in the Appendix. The mean-field model is combined with atomistic simulations to determine α for a given semiflexible polymer. Atomistic simulations of semiflexible oligomers of length S in the isotropic phase are used to determine the bending modulus κ. We apply a tension f to all oligomer backbones in the simulated isotropic phase to induce quadrupolar ordering q. Using our single-chain model, we self-consistently compute the nematic order parameter q(α) of a stretched semiflexible chain with the same stiffness κ, chain length S, and tension f. The value of α(T) at various temperatures is determined by fitting the calculated q(α) to the observed q in the atomistic simulations. MD Simulations of Chains under Tension. To estimate the temperature-dependent nematic coupling constant α(T) for P3HT, we combine our SCFT calculations with atomistic simulations of regioregular 3-hexylthiophene (3HT) 20-mers at different temperatures. The force field parameters for our atomistic simulations were developed by Huang and coworkers.19,20 Bonded interactions and partial charges on the atoms are determined using density functional theory (DFT) calculations. Lennard-Jones (LJ) parameters from the OPLSAA potential are used to describe the non-Coulomb nonbonded interactions between atoms. Huang et al. verified the atomistic model of P3HT by comparing the densities of 3HT monomers and crystal 3HT oligomers at 300 K with experiments.19 We also test the atomistic model by simulating a single crystal of P3HT at 300 K. In this simulation of the crystal, P3HT chains are bonded to themselves through the periodic boundary to represent chains of infinite length. The simulated crystal structure of P3HT, including the π−π stacking distance (3.64 Å) and the tilt angle of backbone thiophenes (22°), agree with X-ray and electron diffraction results.21,22 For our simulations of the isotropic phase, the initial configuration of 3HT 20-mers is built according to the crystal structure of P3HT measured by X-ray diffraction.21 Sixty-four chains are constructed in all-trans conformations with a π−π stacking distance of 3.83 Å and a separation distance of 16.8 Å parallel to the thiophene rings. Sulfur atoms on π−π stacked thiophene rings point in opposite directions (Figure 1). We obtain the isotropic phase of P3HT by melting the crystal at 700 K and 1 bar. The system is equilibrated when the density ρ fluctuates around its equilibrium value, and the order parameter q decays to zero. Starting configurations of isotropic P3HT at 650, 600, and 500 K are obtained by equilibrating at constant pressure and temperature (NPT) after quenching from 700 K. These temperatures are significantly above the

Figure 1. Initial configuration of 64 chains of 3HT 20-mers. Side chains and hydrogens removed for clarity.

experimental crystal melting temperature for 20-mers, observed to be about 450 K.23 Simulations of 20-mers in the isotropic phase at these elevated temperatures can be reliably equilibrated, as discussed further below. Using the equilibrated NPT samples as initial configurations, simulations at constant volume and temperature (NVT) are performed to determine the bending modulus κ for P3HT at different temperatures. NVT simulations are used for consistency with our simulations of chains under tension, for which the NVT ensemble is used to prevent deformation of the simulation box. To induce weak quadrupolar ordering, we apply an equal and opposite force in the x-direction to the ends of each P3HT chain, to induce a tension f of 2.2 kJ/nm along the chains. This value of f corresponds to about 2kT per end-to-end distance of a free chain. Starting as an isotropic melt, each simulation is equilibrated over a total time longer than 20τ, where τ is the autocorrelation time of order parameter q. Equilibrium is achieved in our simulations when q fluctuates around a steady value. We collect data over a time greater than 20τ for each equilibrated simulation. From the simulations of stretched chains at different temperatures, we can determine the temperature-dependent nematic coupling constant α(T), which allows us to predict the location of the IN transition TIN for P3HT. Multichain Free Energy. Our single-chain model provides a mean-field treatment of a multichain system which exhibits an isotropic-to-nematic (IN) transition. To determine the location of the transition, we employ the variational theorem to compute an upper bound to the multichain free energy24 F = ⟨H ⟩ − TS0

(10)

in which ⟨H⟩ is the Hamiltonian of the multichain system. S0 is the model ensemble entropy, defined by the single chain free energy F0 = ⟨ H0⟩ − TS0. The single-chain free energy F0 can be obtained using our single-chain model partition function Z0 = e−βF0 (eq 7). At the critical nematic coupling constant αc, the variational free energy F(q) develops a second minimum F(qc) at critical order parameter qc (qc < 1). Because F(qc) equals the isotropic free energy F(0), the multichain system has a first-order IN transition. To obtain the variational free energy F(q), we first compute the model ensemble entropy S0 as a function of quadrupolar order parameter q. We consider a single semiflexible chain aligned in an external quadrupolar field A. The aligning potential of the single chain now takes the form V(μ) = C

DOI: 10.1021/acs.macromol.5b00013 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules −AP2(μ). The order parameter q is an increasing function of the external field A, which can be computed numerically as a function of A. Operationally, we can compute both q and the single-chain free energy F0 as a function of A, thereby obtaining F0(q) as a function of q. The entropy S0(q) (defined with respect to a free chain with A = 0) is then obtained as −TS0(q) = F0(q) − F0(0) + Aq

(11)

For a multichain system, we write the variational free energy F(q) with respect to the isotropic phase using S0(q): 1 F(q) = − αq2 − TS0(q) 2

(12)

The 1/2 in the nematic coupling term avoids double counting in the multichain model. The variational free energy can be used to determine the critical nematic coupling constant αc at which a given semiflexible polymer undergoes an IN transition. Knowing the temperature-dependent nematic coupling constant α(T) for a given semiflexible polymer, we can predict the IN transition temperature TIN.

■

RESULTS Atomistic Simulations of P3HT. Our initial configurations are oriented single crystals of 64 chains of P3HT, with 20 monomers each (Figure 1). Equilibrating from the initial crystal phase, we observe that P3HT is an isotropic melt above 500 K. We detect nematic order in the simulations using the quadrupolar order parameter q (eq 3). In our equilibrated simulations without applied tension, q (eq 3) fluctuates about zero (Figure 2a), consistent with an isotropic phase. No remnant of the initial crystal orientation survives. From NVT simulations of the unoriented isotropic phase, we determine the persistence length of P3HT. We use the tangent−tangent correlation function ⟨t0·ts⟩ (Figure 3) to estimate the characteristic length of the exponential decay (⟨t0· ts⟩ = e−s/Np), the persistence length Np, and so determine κ (κ = kTNp). As shown in Figure 3, ⟨t0·ts⟩ is essentially temperature independent, with only slight variations in ⟨t0·ts⟩ at large s (which may reflect poor statistics). Thus, we conclude that κ is temperature independent. This is not unexpected: for molten P3HT, the dihedral angles can rotate nearly freely in the isotropic phase. The chain stiffness primarily arises from the small backbone deflection angles, which are insensitive to temperature.25 Even though P3HT is semiflexible, the conjugation along the P3HT backbone, which gives rise to delocalized electronic states observable with UV−vis spectroscopy, is not so easily disrupted. The conjugation is only efficiently broken for dihedral angles near 90° which are relatively rare given the P3HT dihedral potential. We determine the temperature-independent value Np = 7.1 (hence κ = 7.1kT) by averaging the exponential decay lengths of ⟨t0·ts⟩ for the first ten repeat units over different temperatures. The freely rotating chain model yields Np = 7. Fitting the entire range of ⟨t0·ts⟩ gives a slightly different value Np = 8.25 Our estimated Np agrees with the experimentally measured persistence length Lp of 3.0 ± 0.2 nm (corresponding to an Np of 7.5 ± 0.5) for P3HT in a marginal solvent.26 To directly observe the IN transition in simulations, isotropic P3HT needs to be slowly cooled and the nematic observed to form, or the P3HT crystal slowly heated to form a nematic which is then heated further until the isotropic phase forms.

Figure 2. Quadrupolar ordering of molten P3HT: (a) order parameters q for free chain simulations, (b) order parameters q for stretched chain simulations, and (c) autocorrelation functions of order parameters Cqq(t) for stretched chain simulations. Temperatures are 600 K (red), 650 K (green), and 700 K (blue).

Figure 3. Tangent−tangent correlation functions: 600 K (red), 650 K (green), and 700 K (blue).

From Figure 2c, we can estimate the equilibration time for simulations at a given temperature, in terms of the correlation time τ for order parameter fluctuations. At higher temperatures, τ is conveniently short5 ns at 700 K, up to 14 ns at 600 K. But at lower temperatures, the correlation time becomes quite long (over 200 ns at 500 K). Thus, it would be rather challenging to accurately locate the IN transition by directly observing the appearance or disappearance of nematic order. Instead of observing the IN transition directly in simulations, we predict the nematic phase and the coupling constant α of D

DOI: 10.1021/acs.macromol.5b00013 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 4. Snapshots of (a) a free chain simulation of P3HT and (b) a stretched chain simulation of P3HT at 600 K with tension f = 2.2 kJ/nm. Side chains and hydrogens removed for clarity.

P3HT from the tension-induced ordering in the isotropic phase. By applying uniaxial tension, we can induce weak quadrupolar ordering in P3HT in the isotropic phase. Figure 3 shows the snapshots of our MD simulations at 600 K for unstretched chains (a) and stretched chains (b), from which we see that P3HT backbones tend to align along the stretching axis (x-axis). The order parameter q becomes nonzero for the stretched P3HT chains (Figure 2b). Determining α. We analytically compute the quadrupolar order parameter q for a given α for P3HT (Figure 5a) using SCFT. To efficiently compute q(α), the eigenfunction and Legendre polynomial expansions are truncated at fourth order.

The results obtained with this truncation are indistinguishable from q(α) retaining higher order terms. The SCFT results for q(α) are fitted to the MD simulation results to determine α for P3HT at different temperatures (600, 650, and 700 K). We use the three high-temperature simulations because equilibrating at lower temperatures is increasingly expensive, indicated by the increasing autocorrelation time τ for the order parameter q (Figure 2c). The results for α at the three temperatures are fit to the form A + B/T, used previously by Olsen and co-workers, analogous to the form of the Flory− Huggins interaction parameter χ.13 The temperature-independent term A captures the entropic contributions, including the steric exclusions, where B denotes the enthalpic interactions. In this way, we obtain the temperature-dependent α(T) for P3HT as α(T) = −(1.63 ± 0.30) + (1550 ± 190) /T, in which we include the standard errors of the fitting parameters. For sufficiently long chains, SCFT calculations can be simplified using ground-state dominance (GSD). In the long chain limit, only the lowest eigenfunction (ground state) contributes to propagators and partition functions (see eqs 7 and 9). By comparing the order parameter qGSD(α) obtained using GSD to the full SCFT calculation of q(α) at 700 K (Figure 5a), we see that GSD is not a good approximation for oligomeric P3HT. The two results do not agree for small α, indicating the 3HT 20-mer is too short to be described well by GSD. For polymers with small α, the relative error Δq/q from using GSD is less than 10% only when the chain length S is much greater than Np (Figure 5b). The distribution of chain tangent vectors Π(μ) gives a more detailed picture of chain orientation than the order parameter q, which is simply the average of P2(μ) over the distribution Π(μ). (Here μ = cos θ, where θ is the angle between the tangent and the axis of stretching.) We compare the SCFT predictions for Π(μ) to the MD simulation results to further validate our method (Figure 6). The SCFT results agree almost perfectly with the tangent distributions Π(μ) from our MD simulations. The asymmetry of Π(μ) about the origin μ = 0 is a result of the applied tension, while the upward curvature is characteristics of induced nematic order (μ near ±1 is favored over μ = 0). With no applied tension, Π(μ) is symmetric about μ = 0, flat in the isotropic phase, and spontaneously concave up in the nematic phase. Locating the IN Transition. To predict the location of the IN transition, we first calculate the single-chain model entropy

Figure 5. Order parameters computed using SCFT. (a) Order parameter versus nematic coupling constant α: 600 K (red), 650 K (green), and 700 K (blue). The purple curve presents results at 700 K assuming ground state dominance (GSD). (b) Relative error Δq/q for using GSD versus S/Np. S is the chain length, and Np is the persistence length. q is obtained using fourth-order eigenfunction expansions (eq 9) and Δq = qGSD − q. E

DOI: 10.1021/acs.macromol.5b00013 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 6. Tangent distribution functions for stretched 3HT 20-mers at 700 K (blue), 650 K (green), and 600 K (red); symbols are simulation results, and curves are SCFT calculations.

S0(q) for chains with different stiffness κ as a function of an applied quadrupolar aligning field A (eq 11). Because κ characterizes the orientational correlations between the backbone segments along a chain, the entropic cost of anisotropic ordering is lower for chains with stiffer backbones (Figure 7). Figure 8. Free energy vs order parameter for 3HT 20-mers. (a) With various nematic coupling constants: red (α = 1.50 kT), green (α = 1.51 kT), and blue (α = 1.52 kT). (b) At different temperatures: 494 K (red), 493.3 K (green), and 493 K (blue).

Figure 7. Model entropy for semiflexible 20-mers with different stiffness κ.

Using the model entropy S0(q), we can compute the variational free energy F(q) of a semiflexible chain as a function of quadrupolar order q (eq 12). The variational free energy F(q) decreases with increasing nematic coupling α. At a critical value αc, F(q) develops a second minimum at finite order parameter qc, corresponding to the first-order IN transition (Figure 8a). We construct the IN phase boundary using the critical nematic coupling constant αc(κ) for polymers with various stiffness and chain lengths (Figure 9a). The critical nematic coupling constant decreases with increasing chain length or the chain stiffness. By projecting the 3D phase boundary onto the κ−α plane, we obtain the two-dimensional IN phase boundary for different chain lengths. The phase boundaries for different chain lengths are shown in Figure 9b, in which a crossover from rod-like to semiflexible chain behavior can be seen at a chain length of about one persistence length Np(βκ). The value of αc is estimated as 4.54kT for chain length of unity, in agreement with Maier− Saupe theory.27 The 2D phase boundaries are also presented in Figure 9c for chains with different values of S/Np (chain length in units of Np). For chains that are many persistence lengths long (S/Np ≫ 1), the critical nematic coupling constant αc scales as 1/κ. To understand why αcκ is a constant for semiflexible chains much longer than the persistence length Np, we make a simple argument as follows. The mean-field free energy of a multichain

Figure 9. IN phase boundary for semiflexible polymers. (a) Phase boundary for various values of chain length S and stiffness κ. (b) Phase boundary in κ−α plane for various chain lengths S. (c) Phase boundary in κ−α plane for various values of S/Np. Inset: phase boundaries for chains with S/Np ≫ 1.

system can be obtained from the single-chain partition function, which is the integral over chain configurations of the singleF

DOI: 10.1021/acs.macromol.5b00013 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

calorimetry (DSC).28 Hence, our predictions are consistent with DSC experiments that only observe a crystal melting peak for P3HT.14

chain Hamiltonian without applied tension. By expressing the arclength in units of Np(βκ), the single-chain Hamiltonian is written as βH =

∫0

S / Np

■

dz (1/2(dt /dz)2 − β 2ακ t(z) ·Q·t(z))

DISCUSSION In this paper, we report a hybrid method, combining selfconsistent field theory (SCFT) with atomistic simulations of semiflexible chains under tension, to estimate the nematic coupling constant α for semiflexible polymers. We employ a mean-field model to describe the uniaxial alignment of a semiflexible chain under tension. Using the mean-field model, we self-consistently compute the quadrupolar order parameter q(α). The value of α is determined by fitting the predicted q(α) to the simulation results for semiflexible chains under tension in the isotropic phase. In this way, we can obtain the nematic coupling without directly observing the IN transition in simulations or for polymers for which the IN transition is precluded by crystallization. Direct observation of the IN transition in atomistic simulations is challenging because equilibration times at lower temperatures become increasingly long. Thus, the isotropic phase must be cooled very slowly to identify the IN transition accurately, resulting in very expensive simulations. We have applied our method to determine the nematic coupling constant for P3HT. MD simulations of 20-mer chains under tension in the isotropic phase were carried out at 700, 650, and 600 K. At these temperatures, simulations can be equilibrated in less than 300 ns. Values of α at lower temperatures are calculated by extrapolating our results, fitted to the form A + B/T. With the nematic coupling in hand, we can predict the IN transition temperature for P3HT as a function of chain length, using self-consistent field theory. The predicted nematic phase behavior of P3HT agrees with experiments. Ho et al. observed a nematic phase under polarized optical microscopy for P3HT with length S between 43 and 73.14 For such chain lengths, the nematic temperature window is estimated to be about 18 K wide (Figure 10), so that the nematic phase should be accessible experimentally. For longer chains, the nematic phase is stable only in an increasingly narrow temperature range and thus is more difficult to access. Our method is robust to inaccuracies of atomistic force fields, even though the nematic coupling constant α is determined by fitting to atomistic simulation results. We expect the location of the IN transition to be affected primarily by the shape and stiffness of the chains. MD accurately describes local packing, but error in the potential can lead to incorrect overall densities. We have verified that our simulations give comparable quadrupolar alignment (with q decreasing by less than 5%) for chains under tension at 700 K, when the simulation density is artificially reduced by 5% (Figure 11a). Correspondingly, the value of α is reduced by 7% for the system with the reduced density. Likewise, error in the atomistic detailed potential can in principle lead to incorrect values for the chain stiffness. Nevertheless, we have shown that for conjugated polymers the stiffness arises primarily from the small deflection angle, not the details of the dihedral potential.25 To demonstrate this, we set the dihedral potentials to zero at 700 K, corresponding to freely rotating P3HT in the gas phase with Np = 7. In the melt, however, the interchain interactions of the chains with the flat dihedral potentials favor planar conformations along the

(13)

Hence, the free energy F is only a function of ακ, not α or κ independently. For long semiflexible chains, S/Np approaches infinity, and the free energy F develops a second minimum corresponding to the IN transition at a single critical value of α κ, no longer dependent on S/Np. We predict the IN transition temperature TIN for P3HT with various chain lengths S using the temperature-dependent nematic coupling constant α(T). From the variational free energy (eq 12), we determine the critical nematic coupling constant αc for P3HT of different lengths. The free energy of 3HT 20-mers is plotted in Figure 8b as an example. Because Np is temperature independent (Figure 3), TIN is determined by equating α(T) to αc. We obtain the IN phase boundary for P3HT with length S as TIN(S) = 535 K(1 − 1.64/S) by fitting the predicted TIN to the form A(1 + B/S). The uncertainty of the IN phase boundary, as a result of the uncertainties in α(T), is about 15 K. The experimental crystal melting temperature Tm(S) is estimated as 533 K(1 − 3.25/S) by fitting the experimentally observed Tm for P3HT of different molecular weights.14,23 The phase transition temperature TIN is above the crystal melting temperature (Figure 10), indicating the nematic phase should be accessible for P3HT at elevated temperatures.

Figure 10. P3HT phase diagram, with crystal (red), nematic (green), and isotropic phase (blue). The IN phase boundary (dashed blue curve) is TIN(S) =535 K(1 − 1.64/S), obtained by fitting our numerical results (blue circles). The crystal melting temperature (dashed red curve) is Tm(S) = 533 K(1 − 3.25/S), obtained by fitting reported Tm (red circles,14 red triangles.23).

To further compare our predictions with experiments, we predict the IN transition enthalpy HIN for P3HT. We expect a small HIN for semiflexible polymers. At the IN transition, the change in enthalpy HIN equals the loss in orientational entropy TΔS, which is only of order kT per Kuhn segment for semiflexible chains. Because the Kuhn segments are large objects, HIN per volume is small. We predict the IN transition enthalpy HIN using the model entropy S0 and the critical order parameter qc in our mean-field theory. The critical order parameter qc is about 0.12, and HIN is about 0.3 J/g for P3HT. The predicted HIN is much smaller than the enthalpy of fusion for P3HT (99 J/g), measured using differential scanning G

DOI: 10.1021/acs.macromol.5b00013 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

TIN(S) = 535 K(1 − 1.64/S). Using the estimated TIN, we predict that P3HT is nematic above the reported crystal melting temperatures.

■

APPENDIX By expanding the propagator Z(μ;s) in eigenfunctions (eq 9), we rewrite the diffusion equation (eq 8) as ϵαψα(μ) = −D∇⊥2 ψ (μ) + V (μ)ψ (μ)

(14)

The eigenfunctions ψα(μ) (eq 9) is further expanded in Legendre polynomials Pn(μ): N

ψα(μ) =

∑ cα ,n n=0

2n + 1 Pn(μ) 2

(15)

in which cα,n are the expansion coefficients for ψα(μ). Both eigenfunction expansion of Z(μ;s) and Legendre polynomial expansion of ψα(μ) are truncated at the same order N. The Legendre polynomials are orthogonal with respect to integration over all μ: ∫ 1−1dμ Pn(μ)Pm(μ) = 2δmn/(2n + 1). By taking the inner product of ((2m + 1)/2)1/2Pm(μ) with eq 14 and integrating both sides of the equation over μ, we write the diffusion equation as

Figure 11. Dependence of quadrupolar ordering on force fields at 700 K: (a) quadrupolar order parameters q; (b) tangent−tangent correlations functions. Our regular simulation with density ρ0 (hindered), a less dense system with density ρ = 0.95ρ0 (green), and freely rotating chains (blue).

ϵαcα , m = Dm(m + 1)cα , m +

∑ ⟨m|V (μ)|n⟩cα ,n n

(16)

Here the “|n⟩” represents ((2n + 1/2)) Pn(μ). The notation ⟨m|n⟩ denotes the inner product of |m⟩ and |n⟩, which by orthonormality is δmn. Likewise, ⟨m|n⟩ denotes the angular average “matrix element” of V(μ) between |m⟩ and |n⟩. For convenience, we write the second term on the right-hand side of eq 16 as Vmn. Using the recurrence relation of the Legendre polynomials Pn,18 we simplify Vmn (eqs 4 and 16) as 1/2

backbones, resulting in a persistence length Np of 8.3 (Figure 11b). For the artificially stiffened chains, the induced quadrupolar ordering q is about 8% higher (Figure 11a). The resulting nematic coupling constant α at 700 K is about 23% higher than the value of α for the hindered-rotating P3HT. Because the relative error in Np of the stiffened P3HT is rather large (17%), we expect our method is robust to milder errors in the dihedral potentials. Our method can be also extended to polymers other than the semiflexible linear homopolymers. The extension to alternating copolymers is straightforward, and low levels of long-chain branching would be only a perturbation. For alternating copolymers, an effective persistence length can be determined. Together with the size of the repeat unit, the effective persistence length can be used in our method for predicting α. For semiflexible polymers with low levels of long-chain branching, the joints have little effect on the overall aligning behaviors so that our method is still applicable. In summary, we compute the variational free energy F(q) of a multichain system using our mean-field model. We consider a single chain immersed in an external quadrupolar aligning field A, to compute the entropy S0 as a function of the order parameter q (eq 11). By adding a nematic coupling term to the entropy S0(q), we obtain F(q) with respect to free chains (eq 12), from which the IN transition can be located. We demonstrate our method by estimating the nematic coupling constant α(T) and the IN transition temperature TIN for P3HT. Combining the SCFT calculations with the atomistic simulations of stretched 3HT 20-mers at different temperatures, we determine α as −(1.63 ± 0.30) +(1550 ± 190)/T for P3HT. The IN phase boundary is obtained for semiflexible chains using the variational free energy. The critical nematic coupling constant αc is lower for longer or stiffer chains. Because the bending stiffness κ is temperature independent, we directly estimate TIN for P3HT with various lengths S, given by

⎛ 3(n + 1)(n + 2) δm , n + 2 Vmn = −αq⎜ ⎝ 2(2n + 1)(2n + 3) n(n + 1) 3n(n − 1) δm , n + (2n − 1)(2n + 3)) 2(2n + 1)(2n − 1) ⎞ ⎛ n+1 ⎞ n × δm , n − 2⎟ − f ⎜ δm , n + 1 + δm , n − 1⎟ ⎝ 2n + 1 ⎠ 2n + 1 ⎠ (17)

+

The right-hand side operator of eq 16 becomes an N × N matrix: Mmn = Dm(m + 1)δmn + Vmn, for which the eigenvalues are ϵα and the corresponding eigenvectors are {cα,m}. The eigenvalues ϵα and eigenvectors {cα,m} are used to construct the propagator Z(μ;s) (eqs 15 and 9) so that the tangent distribution Π(μ) (eq 6) and quadrupolar order parameter q (eq 3) can be determined self-consistently for a given α.

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail [email protected] (E.D.G.). *E-mail [email protected] (S.T.M.). Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Onsager, L. Ann. N. Y. Acad. Sci. 1949, 51, 627−659. (2) Khokhlov, A. R.; Semenov, A. N. Physica A 1981, 108, 546−556. (3) Khokhlov, A. R.; Semenov, A. N. Physica A 1982, 112, 605−614. (4) Grell, M.; Bradley, D. D. C.; Inbasekaran, M.; Woo, E. P. Adv. Mater. 1997, 9, 798−802.

H

DOI: 10.1021/acs.macromol.5b00013 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules (5) Wang, X.; Engel, J.; Liu, C. J. Micromech. Microeng. 2003, 13, 628. (6) Donald, A. M.; Windle, A. H.; Hanna, S. Liquid Crystalline Polymers; Cambridge University Press: New York, 2006. (7) Woltman, S. J.; Jay, G. D.; Crawford, G. P. Nat. Mater. 2007, 6, 929−938. (8) Redecker, M.; Bradley, D.; Inbasekaran, M.; Woo, E. P. Appl. Phys. Lett. 1999, 74, 1400−1402. (9) O’Neill, M.; Kelly, S. M. Adv. Mater. 2003, 15, 1135−1146. (10) Kinder, L.; Kanicki, J.; Petroff, P. Synth. Met. 2004, 146, 181− 185. (11) McCulloch, I.; Heeney, M.; Bailey, C.; Genevicius, K.; MacDonald, I.; Shkunov, M.; Sparrowe, D.; Tierney, S.; Wagner, R.; Zhang, W.; Chabinyc, M. L.; Kline, R. J.; McGehee, M. D.; Toney, M. F. Nat. Mater. 2006, 5, 328−333. (12) Morse, D. C.; Fredrickson, G. H. Phys. Rev. Lett. 1994, 73, 3235−3238. (13) Olsen, B. D.; Jang, S.-Y.; Lüning, J. M.; Segalman, R. A. Macromolecules 2006, 39, 4469−4479. (14) Ho, V.; Boudouris, B. W.; Segalman, R. A. Macromolecules 2010, 43, 7895−7899. (15) Park, M. S.; Aiyar, A.; Park, J. O.; Reichmanis, E.; Srinivasarao, M. J. Am. Chem. Soc. 2011, 133, 7244−7247. (16) Jiang, Y.; Chen, J. Phys. Rev. E 2013, 88, 042603. (17) Jiang, Y.; Chen, J. Phys. Rev. Lett. 2013, 110, 138305. (18) Arfken, G. B.; Weber, H. J.; Harris, F. E. Mathematical Methods for Physicists; Academic Press: New York, 2013. (19) Huang, D. M.; Faller, R.; Do, K.; Moulé, A. J. J. Chem. Theory Comput. 2010, 6, 526−537. (20) Schwarz, K. N.; Kee, T. W.; Huang, D. M. Nanoscale 2013, 5, 2017−2027. (21) Prosa, T. J.; Winokur, M. J.; Moulton, J.; Smith, P. Macromolecules 1992, 25, 4364−4372. (22) Kayunkid, N.; Uttiya, S.; Brinkmann, M. Macromolecules 2010, 43, 4961−4967. (23) Wu, Z.; Petzold, A.; Henze, T.; Thurn-Albrecht, T.; Lohwasser, R. H.; Sommer, M.; Thelakkat, M. Macromolecules 2010, 43, 4646− 4653. (24) Feynman, R. Statistical Mechanics: A Set of Lectures (The Advanced Book Program); Perseus Publishing: New York, 1998. (25) Zhang, W.; Gomez, E. D.; Milner, S. T. Macromolecules 2014, 47, 6453−6461. (26) McCulloch, B.; Ho, V.; Hoarfrost, M.; Stanley, C.; Do, C.; Heller, W. T.; Segalman, R. A. Macromolecules 2013, 46, 1899−1907. (27) Humphries, R. L.; James, P. G.; Luckhurst, G. R. J. Chem. Soc., Faraday Trans. 2 1972, 68, 1031−1044. (28) Malik, S.; Nandi, A. K. J. Polym. Sci., Part B: Polym. Phys. 2002, 40, 2073−2085.

I

DOI: 10.1021/acs.macromol.5b00013 Macromolecules XXXX, XXX, XXX−XXX