Blends of Semiflexible Polymers: Interplay of Nematic Order and Phase Separation

Mixtures of semiflexible polymers with a mismatch in either their persistence lengths or their contour lengths are studied by Density Functional Theory and Molecular Dynamics simulation. Considering lyotropic solutions under good solvent conditions, the mole fraction and pressure is systematically varied for several cases of bending stiffness κ (the normalized persistence length) and chain length N. For binary mixtures with different chain length (i.e., NA=16, NB=32 or 64) but the same stiffness, isotropic-nematic phase coexistence is studied. For mixtures with the same chain length (N=32) and large stiffness disparity (κB/κA=4.9 to 8), both isotropic-nematic and nematic-nematic unmixing occur. It is found that the phase diagrams may exhibit a triple point or a nematic-nematic critical point, and that coexisting phases differ appreciably in their monomer densities. The properties of the two types of chains (nematic order parameters, chain radii, etc.) in the various phases are studied in detail, and predictions on the (anisotropic) critical behavior near the critical point of nematic-nematic unmixing are made.


Introduction
Semiflexible polymers are macromolecules with linear chemical architecture in the simplest case, which display considerable bending rigidity along the chain backbone [1][2][3][4]. In lyotropic solutions of semiflexible polymers, the competition between orientational and translational entropy can induce a transition from an isotropic (i) to a nematic (n) phase. Depending on the chemical nature of the monomeric units, the persistence length p of such macromolecules can vary over several orders of magnitude: For synthetic polymers, p typically lies in the range from 1 nm to 30 nm [5], while the distance between repeat units along the backbone, b , is typically on the order of several angstroms. Much larger values for these lengths may occur for biologically relevant polymers, e.g., p ≈ 50 nm for double stranded DNA [6] and p ≈ 17 µm for filamentous actin [7]. Note that p does not only depend on the specific polymer chemistry [8] but may also depend on external factors, e.g., the nature of the solvent [6], and the density of monomeric units in a nematically ordered solution or melt [9].
The ability to control the materials properties of polymeric systems is highly desirable for many practical applications. It is common to do this by mixing two chemically different polymers, e.g., blending poly(phenylene ether) resins with polystyrene results in materials with high heat resistance and strong mechanical stability [10]. The statistical thermodynamics of such polymer blends is usually studied via the Flory-Huggins (FH) solution theory [11][12][13][14][15][16], which describes the polymer mixture using a lattice model. Here, the entropy of mixing favors homogeneous systems even for very long chains, whereas unmixing is driven by the chemical incompatibility of the polymers, typically quantified through the FH χ parameter. One central approximation of the FH solution theory is the mean-field description of the energy of mixing, which completely disregards the (local) chain connectivity. As a consequence, this contribution is identical for regular solutions, polymer solutions, and blends [15]. Further, chain stiffness does not enter at all in the standard FH solution theory, so that it cannot distinguish between fully flexible and stiff polymers. Therefore, alternative approaches have been developed to describe rod-coil mixtures, where the flexible coil-like polymer is often modeled as an effective soft spherical particle [17][18][19][20]. Although such coarse-grained descriptions can provide important insights on the (qualitative) phase behavior, they completely lose information on the scale of monomeric units, and are also unsuitable for describing blends of semiflexible polymers, which can neither be described as strictly hard rods nor as soft spheres.
In this work, we consider lyotropic solutions containing two types (A, B) of semiflexible polymers with finite persistence lengths B p > A p b . This problem has been previously considered by a few analytical approaches, yielding interesting explicit results only for special limiting cases. Semenov and Subbotin [21] generalized the Onsagertype [22] treatment of single-component lyotropic dilute polymers solution [23][24][25][26]  b , they predicted interesting (qualitative) phase diagrams containing both mixed and coexisting nematic phases [21]. Experimentally, the coexistence of two nematic phases n 1 and n 2 was observed in mixed virus suspensions [27], which can be modeled theoretically by mixtures of hard rods differing in diameter [28,29]. In recent MD simulations, Zhou et al. studied the behavior of mixtures of semiflexible (ring) polymers in spherical [30] and ellipsoidal [31] containers, finding a confinement-induced phase separation of the two species. In thermotropic systems, n 1 -n 2 coexistence was also predicted for solutions of two kinds of rigid rods with suitable enthalpic interactions [32], and experimental observations of n 1 -n 2 unmixing were also reported for mixtures of sidechain liquid crystalline polymers with small molecule liquid crystals [33], but all such temperature-driven phase transitions are out of consideration here. Thus, we also do not discuss in detail the theory of Liu and Fredrickson [34], which was based on a Landau expansion in terms of both the nematic order parameter and deviations of the local volume fraction of monomeric units from its average. Since Landau expansions require that order parameters are small, this approach is not suitable deeply in the nematic phase. In addition, it was assumed that the transitions were driven by standard enthalpic interactions throughout, i.e., i-i unmixing due to a standard FH χ parameter, and nematic ordering due to a Maier-Saupe-like term [35].
In the present work, we focus on systems where transitions are solely driven by entropic interactions, focusing on cases where contour and persistence lengths are of the same order. Preliminary results of our work were already presented as a letter [36]; in a complementary study [37], we have described related results for shorter chains, with emphasis on a comparison between polymers with finite stiffness and strictly hard rods, and the classification of various contributions to the Gibbs excess free energy within the Density Functional Theory (DFT) framework. In Section 2, we briefly characterize the employed models and methods. Following our previous work on one-component solutions of semiflexible polymers [9,[38][39][40], we use both DFT calculations and Molecular Dynamics (MD) simulations for closely related models. Section 3 gives an overview of the possible phase diagrams, emphasizing their description in different statistical ensembles, as predicted by DFT. In Section 4, we analyze in more detail selected systems that were studied by both DFT and MD. For one example, we show how a stiffness mismatch leads to an effective χ parameter of entropic origin, by applying an approach pioneered by Fredrickson et al. [41] and Kozuch et al. [42]. We also discuss the anisotropic character of critical fluctuations for the n-n critical point. Finally, Section 5 summarizes our results and briefly mentions pertinent experiments [43,44].

Models and Recorded Observables
In the context of DFT [9,[36][37][38][39][40], it is convenient to describe the polymer chains as sequences of N tangent hard spheres of diameter σ. Bending stiffness is introduced via a bending potential U bend (θ ijk ), which depends on the angle θ ijk between subsequent bond vectors a i = r i − r j and a j = r k − r j , with monomer positions r µ , µ = 1, . . . , N. This potential is chosen as where we anticipated that κ ≡ ε bend /(k B T) 1, hence small angles θ ijk dominate. The dimensionless parameter κ controls the persistence length p of the chains, which is traditionally defined from the exponential decay of bond orientational correlations, a i · a i+s ∝ exp(−s b / p ) [15,45], with bond length b = σ here. However, this exponential decay with s is doubtful in various contexts [8]; in particular, it does not apply for large s when nematic order is present. But even then, the case s = 1 can be used to define an effective persistence length eff p : For simplicity, the superscript "eff" in Equation (2) will be omitted in the following. The nematic order parameter S is defined as the largest eigenvalue of the traceless tensor Q αβ (α, β = x, y, z) with unit vector u i along a i . The nematic director is defined from the eigenvector belonging to the largest eigenvalue of Q αβ . In a system with a single kind of polymer, the average · · · in Equation (2) is taken over all bond angles of all chains, and over all unit vectors in Equation (3). Correspondingly, in a two-component system with two kinds of polymers A and B (we use σ A = σ B = σ and A b = B b = σ, but have κ A < κ B ), only A chains are included in Equations (2) and (3) for computing A p and S A , respectively, and likewise for the B chains.
In the MD simulations, all monomeric units have identical masses m A = m B = m. Excluded volume interactions between beads are described by the Weeks-Chandler-Andersen (WCA) potential [46] where r is the center-to-center distance between a pair of monomers, and ε sets the energy scale for the repulsion. Neighboring beads along a chain are bonded through the finitely extensible nonlinear elastic (FENE) potential [47] with spring constant k = 30 ε/σ 2 and maximum bond extension r 0 = 1.5 σ to prevent unphysical chain crossing. All parameters in Equations (4) and (5) are chosen the same for both chain types, resulting in a bond length b ≈ 0.97 ± 0.03 σ. MD simulations were performed either in the N VT ensemble (N being the total number of monomeric units) with a Langevin thermostat [47], or in the N PT ensemble [48,49] using the HOOMD-blue software package (v. 2.9.4) [50]. A time step of ∆t = 0.005 τ MD was chosen to numerically integrate the equations of motion, with intrinsic MD time unit τ MD = √ mσ 2 /ε. To achieve thermal equilibrium and reach the desired statistical accuracy, 10 8 to 10 9 integration steps were performed. Initial configurations were created by placing the A and B chains in rod-like configurations, perfectly stretched out in the z-direction (all bond angles being zero). The systems contained typically 10 5 to 10 6 monomeric units, depending on the desired average density ρ and composition X B . Periodic boundary conditions were employed in all Cartesian directions, and the size of the simulation box was chosen distinctly larger than the polymer contour lengths L to avoid significant finite size effects. Figure 1 shows examples of snapshot pictures of the simulated systems. In the MD simulations, mean square end-to-end distances R 2 and gyration radii R 2 g are also obtained for both chain types. In nematic phases, components parallel ( R 2 , R 2 g ) and perpendicular ( R 2 ⊥ , R 2 g ⊥ ) to the director must then be distinguished. We chose the coordinate system such that the nematic director lies parallel to the z-axis, so that R 2 ≡ R 2 z and R 2 ⊥ ≡ ( R 2 x + R 2 y )/2 (and analogous for R 2 g ). We also consider the single-chain structure factor [15,45] For small |q|, we have the expansions [51] Again, indices A and B need to be introduced in the mixture for the two types of polymers appropriately. Finally, we mention that both the tangent hard sphere chain and the bead spring model with bending potential U bend (see Equation (1)), can be considered as discretized versions of the Kratky-Porod worm-like chain model [52,53]. The latter model is recovered in the limit N → ∞, b → 0, and σ → 0, while keeping the contour length L = (N − 1) b fixed.

DFT Methods
For semiflexible polymers, DFT expresses the free energy of the system as a functional of the orientational distribution function f (ω). Here, ω is a shorthand for the two polar angles ϑ and φ describing the orientation of a molecule, which is defined by the unit vector belonging to the smallest eigenvalue of the moment of inertia tensor [54]. The theory extends Onsager's theory for lyotropic solutions of hard rods [22] and describes the nonuniform molecular density as a product of the average density ρ mol = N p /V and f (ω), with number of polymers in the system N p and system volume V. The Helmholtz free energy per molecule is then split into an ideal contribution and an excess contribution where V excl (ω, ω ) is the effective excluded volume between two semiflexible polymers. One can simplify the description further noting that V excl depends only on the relative angle γ between the two chains. While V excl (γ) can be computed analytically for two isolated hard rods [22], the conformational degrees of freedom of semiflexible polymers necessitate a numerical approach. Following Fynewever and Yethiraj [54], we performed additional two-chain Monte Carlo (MC) simulations to determine V excl (γ) for the AA, AB, and BB chain pairs (see Refs. [37,54] for more details). Note that computing V excl (γ) from two interacting chains at infinite dilution and using it in the approximation based on the second virial coefficient, Equation (10), is self-consistent only in dilute solutions. To extend the approach to larger densities ρ mol , the prefactor of V excl (γ) is enhanced by Parsons-Lee [55,56] rescaling. Clearly, this procedure is purely heuristic, but extensive comparisons with MD simulations have shown that the accuracy of these approximations is reasonable at least for single-component solutions [9,38,39]. Finally, note that the nematic order parameter S p derived in this coarse-grained representation refers to the chain ordering as a whole, and not to the order parameter S associated with individual bond vectors.

Nematic-Nematic Phase Separation in Binary Mixtures of Semiflexible Polymers Described by the Random Phase Approximation
In our preliminary communication [36], we have presented evidence that a homogeneously mixed binary nematic phase can phase separate into A-and B-rich nematic phases at a critical point (see Figure 1). Thus, one expects that associated critical fluctuations should occur in the homogeneous phase near this critical point. In this section, we present a phenomenological (mean-field type) treatment of these critical phenomena within the random phase approximation (RPA) [12,15,57,58], to guide the analysis of the DFT and MD numerical work.
The RPA considers the collective structure factor G coll (q), which experimentally is accessible by scattering under wavevector q. In the case of an incompressible mixture, RPA relates G coll (q) to the single-chain structure factors G A (q) and G B (q) as Here, χ eff is the appropriate effective FH interaction parameter, and φ A and φ B are the volume fractions of A and B chains, respectively, with φ A + φ B = 1. In a symmetric mixture (N A = N B = N) that is nematically ordered, Equations (6)- (8) and (11) yield where the quantity G −1 coll (0) describes the volume fraction fluctuations, The correlation lengths ξ , ξ ⊥ of these volume fraction fluctuations parallel and perpendicular to the director of the nematic ordering field are Equation (13) implies that the critical composition is φ A = φ B = 1/2 and criticality occurs for χ eff = χ c = 2/N [11,12,15]. Therefore, one finds the critical amplitudesξ ,ξ ⊥ then being related to the averages of the corresponding mean square gyration radii componentŝ For a well-ordered nematic system, a chain occupies a region that can be approximately described as a cylinder of height H = N b and radius R = 1/(π b ρ) 1/2 . This estimate comes from the condition that the volume of the cylinder πR 2 H contains only the N monomeric units of the chain, so the density inside of the cylinder is just the average density [38,39]. Hence, for dense systems, ξ is of order L, while ξ ⊥ is of order b . However, for very stiff chains with p L b , we know that the i-n transition occurs already for densities ρ of order 1/N, and then also the density, where n 1 -n 2 unmixing sets in, is of the same order. Then, we predictξ ⊥ ∝ R ∝ b √ N, whileξ always is of order b N. Thus, RPA predicts critical fluctuations of a very anisotropic character, and this consideration motivates some of the analyses that will be discussed in the next section. Figure 2 shows phase diagrams in the space of the intensive thermodynamic variables pressure, P, and chemical potential difference, ∆µ, for N A = N B = N = 32. When we fix κ B and vary κ A , there is a particular "multicritical" value κ M A (estimated here as κ M A = 20.5 for κ B = 128 and N = 32), where the topology of the phase diagram changes: For κ A < κ M A , the (first-order) transition between two distinct nematic phases ends at a triple point, where (first-order) i-n phase separations into A-or B-rich nematic phases set in. However, for κ A > κ M A , the first-order n 1 -n 2 phase separation ends at a critical point instead. It is interesting that, for κ A < κ M A , the two first-order lines for the i-n transitions meet in the (P, ∆µ) plane at the triple point under some angle, while, for κ A = κ M A , all transition lines seem to have a common tangent at the triple point. Note also that the value κ M A where the phase diagram topology changes is not universal but depends on both κ B and N [36]; qualitatively similar results were found in related work for shorter chains (N = 16) [37], but there the critical point occurred at distinctly larger pressures, where the unmixing of nematic phases could be potentially preempted by the appearance of smectic and/or crystalline phases. The polymers with N = 32 studied here are just large enough so that one can expect n 1 -n 2 critical points at monomer densities that are low enough so that simulation studies, and perhaps also experiments, become feasible.  These special features have their counterparts in the phase diagrams, where one chooses X B rather than ∆µ as a variable (Figure 3). At the multicritical point κ A = κ M A , the phase boundary of the i-n coexistence region has a horizontal tangent on the nematic side and touches the n 1 -n 2 two-phase region at its critical point. For κ A > κ M A , however, there is a single lens-shaped i-n coexistence region, which extends from the pure A system (X B = 0) to the pure B system (X B = 1). In this representation, the n 1 -n 2 coexistence curve has an approximate parabolic shape near the critical point, which occurs at the minimum of this curve in the (P, X B ) plane. Such a parabolic shape reflects a mean-field critical exponent, as expected for DFT. For κ A > κ M A , this coexistence region no longer interferes with i-n coexistence. With increasing P, the n 1 -n 2 coexistence will end when other phases (smectic liquid crystalline or crystalline solid phases) come into play, but such phases can not be captured by the present DFT treatment.

DFT Predictions of Possible Phase Diagrams
In mean-field theories of critical unmixing of binary mixtures, the response function [59] not only diverges at the critical point (X c B , P c ), but also along the whole "spinodal curve" X sp B (P), which touches the coexistence curve at the critical point. In the phenomenological FH theory discussed in Section 2.3, where χ rather than P was used as a control parameter, the spinodal curve is simply given when we require G −1 coll (0) = 0 in Equation (13), since C ∝ G coll (0) [59,60]. Although the concept of a spinodal curve is doubtful outside of meanfield theories [61], we have included its location in Figure 3c nevertheless. Generally of interest, however, is the so-called "Widom line", describing the locus of maxima of C(X B , P); a short piece of this line is also shown in Figure 3c. Figure 3d shows the molecular nematic order parameters S p A and S p B of the two types of chains plotted versus the pressure P, for a mole fraction X B = 0.5. A clear advantage of DFT is that it is straightforward to discuss the phase behavior of the considered systems in various ensembles of statistical mechanics, which is particularly useful for making contact with experiments: The (osmotic) pressure P of a polymer solution is usually not readily accessible, and one rather uses the polymer concentrations as variables.
In our implicit solvent model, these concentrations translate into the monomer number densities ρ A and ρ B , respectively. Alternatively, we may take the total density ρ = ρ A + ρ B and the mole fraction X B = ρ B /ρ as variables to draw the phase diagram (Figure 4a,b).
While coexisting phases in equilibrium must be at the same temperature and pressure, their densities can differ, as clearly shown in Figure 4a,b. Note also that the critical point no longer coincides with the minimum of the two-phase coexistence curve of n 1 -n 2 phase separation in the (X B , ρ) phase diagram. While the tie lines, connecting the coexisting phases, are almost parallel to each other for n 1 -n 2 coexistence, this is not true for i-n coexistence: In this case, the tie lines must get perpendicular to the X B -axis when X B → 0 and X B → 1, but are much flatter in between. To avoid a too confusing picture, we have not shown any tie lines in Figure 4c, where the phase diagram is shown in the (ρ A , ρ B ) plane; in that representation, three phases must coexist for any state within the triangle formed by the three dotted lines enclosing the triple region.
In Figure 5a, we present the response function C (Equation (18)) in the homogeneously mixed nematic phase, plotted versus X B for various pressures smaller than the critical pressure P c . The coordinates of the maxima of these curves yield the location of the "Widom line" in the (X B , P) plane, as shown in Figure 3c. The log-log plot of the inverse height of this maximum, C −1 max , reveals the expected critical variation, C −1 max ∝ P c − P (Figure 5b).    Figure 3a in this representation correspond to a three-phase triangle, enclosed by the three dotted straight lines, with triangles marking its corners. The region at small densities near the origin up to the full curves contains the isotropic phase; for large ρ B and small ρ A , there is a homogeneous nematic phase rich in B chains; for large ρ A and small ρ B , there is a homogeneous nematic phase rich in A chains. In between the two dashed green curves, n 1 -n 2 two phase coexistence occurs. The white region underneath the three-phase triangle is the isotropic-B-rich nematic two-phase region, and above it the isotropic-A-rich nematic two-phase region.

Selected MD Results for the Phase Behavior of Mixtures of Semiflexible Polymers
We start this section with a rather simple special case, i.e., two types of chains with the same stiffness, κ A = κ B = 32, but different chain lengths, N A = 16 and N B = 32. We have studied this system in cubic boxes, containing 12544 A and 6272 B chains (X B = 0.5). We systematically varied the pressure P and measured the average number density of monomeric units. Figure 6a shows the resulting equation of state, compared with the corresponding isotherms of the pure systems (X B = 0 and X B = 1). In the mixed systems, a sudden increase of the density from ρ ≈ 0.342 σ −3 to ρ ≈ 0.366 σ −3 takes place at P ≈ 0.215 − 0.220 k B Tσ −3 , where the nematic order parameters S A and S B (Figure 6b) and end-to-end distances (Figure 6c) indicate the transition, as well.  While pure systems in the N PT ensemble should exhibit a sharp first-order i-n transition, where the density ρ and the nematic order parameter S have a well-defined jump, this is not the case for mixtures: There, the phase diagram must have in general the shape of a lens, similar to the i-n two-phase coexistence region of Figure 3c. We would have a unique transition pressure P trans only for a truly intensive thermodynamic variable, such as ∆µ, as in Figure 2; however, since X B is formed from densities ρ A and ρ B of extensive variables, we rather have a two-phase coexistence region again for fixed X B , with P i (X B ) < P < P n (X B ). The fact that we seem to observe a unique, well-defined transition pressure in Figure 6 instead needs to be interpreted by the hypothesis that P i (X B ) and P n (X B ) are so close together that we cannot resolve the difference.
DFT respects these general rules from thermodynamics, but it suffers from other problems: The formulation which we used in this work does not give any information on chain linear dimensions, since it is just based on the distribution f (ω) of coarse-grained chains (see Section 2.2). As a consequence, also the nematic order parameter S p computed from DFT (Figure 3d) refers to the average orientation of the whole chain, and not to the orientation of individual bonds, as defined in Equation (3) and shown in Figure 6. It has been demonstrated [62] that there is a significant difference between the nematic order parameter defined from bond orientations and its counterpart from the orientation of whole chains. Thus, we do not have any DFT counterparts to the simulation data of Figure 6b,c; with respect to Figure 6a, even away from the i-n transition, the agreement is only qualitative but not quantitative.
An interesting question concerns the comparison of data for the nematic order parameters and chain radii in the nematic phase for the mixed system (X B = 0.5) with their pure counterparts: We find that S B (X B = 0.5) < S B (X B = 1), i.e., the admixture of the shorter chains, which exhibit somewhat less nematic order, weakens the order of the longer chains slightly. This is also evident from the fact that S A (X B = 0.5) < S B (X B = 0.5) throughout the nematic phase. On the other hand, we also have S A (X B = 0.5) > S A (X B = 0) for those pressures where nematic order occurs. Those chains which order better (here, the longer chains) act like a nematic ordering field on the chains which have less tendency to order (in the pure phase). We shall see later that the same phenomenon occurs when we examine mixtures of chains that have identical N but differ in stiffness.
Throughout the isotropic phase, we find here R 2 g,A = 18.15 σ 2 and R 2 g,B = 66.4 σ 2 , irrespective of P and X B , and likewise for the mean square end-to-end distances, R 2 A ≈ 182 σ 2 and R 2 B ≈ 669 σ 2 . Since the persistence length is about twice as large as the contour length for the A chains, only a rather modest increase of R 2 A can be observed at the i-n transition (Figure 6c). In contrast, the longer B chains exhibit a marked jump of their mean square end-to-end distances at the i-n transition. Further note that R 2 B for P = 1 k B Tσ −3 has already reached 96% of the contour length, at which this quantity must saturate for large enough pressures. Now, we turn to some of the cases that were already studied by DFT in Section 3, where chains have the same chain length N A = N B = 32, but differ in stiffness, e.g., κ A = 24 and κ B = 128. From DFT, we would expect a transition from the isotropic phase into the i-n two-phase coexistence region with increasing pressure, then a nematic phase to which both types of chains contribute, followed by n 1 -n 2 unmixing at still higher pressures (Figure 3c). Figure 7a shows the density versus pressure curve of the equimolar mixture (X B = 0.5) from MD simulations, compared to the corresponding pure systems. For the latter, the i-n transitions show up as little kinks at P trans ≈ 0.24 k B Tσ −3 and 0.07 k B Tσ −3 for the A and B chains, respectively, with corresponding monomer number densities being ρ i ≈ 0.37 σ −3 and ρ n ≈ 0.38σ −3 for species A, while ρ i ≈ ρ n ≈ 0.21 σ −3 for species B. Unlike DFT, we cannot resolve the density jump at the transition with meaningful accuracy in the MD simulations. Nevertheless, the corresponding DFT results for the transition pressures, i.e., P trans = 0.145 k B Tσ −3 and 0.035 k B Tσ −3 for the pure A and B systems, respectively, are clearly far outside of the MD error bars. Thus, there is almost a discrepancy by a factor of two in the pressure scale. However, the underestimation by DFT with respect to the densities is much smaller (Figure 4). When we examine the ρ versus P curve for the mixture, the i-n transition cannot be recognized from the data at all (Figure 7a). This behavior can be expected when the two-phase coexistence region is not very narrow (as it was in Figure 6), but rather wide, as found by DFT for the present case (Figure 3c). Then, the density variation in this two-phase region is described by the lever rule, where X B = X i B (P) and X B = X n B (P) describe the curves limiting the two-phase region at the isotropic and nematic side, respectively. At these curves, ρ(X B , P) has at most a (small) discontinuity of slope, which clearly could not be seen in the simulation. In order to obtain a simulation counterpart (Figure 7b) to the i-n miscibility gap predicted from DFT in Figure 3c, we had to carry out MD simulations for various choices of X B , and analyze the compositions observed in suitably placed x × L y × L z subsystems with x L x = 3L y ; "suitably placed" means that simulation snapshots were inspected to check for states where two-phase coexistence occurs with interfaces perpendicular to the x-axis (see Figure 1b). The subsystems were then centered far from these interfaces, and compositions (as well as other observables) were recorded separately for the isotropic and nematic subsystems. To avoid systematic errors, we have checked that the diffusion of the domain walls in the x-direction was small enough during the time interval in which subsystem averages were taken. However, rather large statistical fluctuations of all observables are inevitable in such simulations of two-phase coexistence. Further note that we did not succeed to prepare reasonably stable states of two-phase coexistence for pressures P ≤ 0.085 k B Tσ −3 and P ≥ 0.175 k B Tσ −3 , but rather metastable homogeneous nematic or isotropic states took over in those pressure regions. Hence, the connections of the observed phase boundaries in Figure 7b to the pure system phase transitions at X B = 0 and X B = 1 are tentative interpolations only.
In principle, one could improve the results by running considerably larger systems and/or longer simulation times, which would, however, require a prohibitively large investment of additional computational resources. Thus, while MD simulations yield in principle the exact statistical mechanics of the considered model system, one has to be aware of its practical limitations. Nevertheless, our MD simulations confirm the DFT prediction that a rather wide i-n miscibility gap occurs for intermediate values of X B (Figure 3c), unlike the system shown in Figure 6. Figure 7b suggests that, for X B = 0.5, the two-phase region extends from about P ≈ 0.08 k B Tσ −3 to about 0.14 k B Tσ −3 , compatible with Figure 7a.
Since DFT predicted for the case κ A = 24, κ B = 128 a region of n 1 -n 2 unmixing for P > P c ≈ 0.2 k B Tσ −3 (Figure 3c), or ρ > ρ c ≈ 0.43 σ −3 (Figure 4b), we extended our MD studies to significantly larger pressures (and corresponding densities) than shown in Figure 7, in order to search for MD evidence of n 1 -n 2 phase separation in this model. To this end, we performed additional simulations in elongated boxes of size 128 σ × 64 σ × 64 σ, with hard walls placed at x = −64 σ and x = +64 σ to stabilize nematic order. Starting configurations were generated by orienting all chains along the z-axis, and placing the A and B chains in the right (x > 0) and left half (x < 0) of the simulation box, respectively. However, we found that the initially separated A-and B-rich phases completely decay by interdiffusion. Even for ρ = 0.65 σ −3 , a state deep in the nematic phase for both pure A and pure B systems, the system still develops towards a homogeneous mixture in the region far from the confining walls. At still larger densities, the system is too slowly relaxing to clearly establish equilibrium, and smectic order starts to take over in the pure B phase [40].
Since the occurrence of n 1 -n 2 unmixing can also be seen already in the homogeneous phase by a growth of the response function C(X B , P) as P → P c (Equation (18)), cf. Figure 5a, we attempted to estimate this quantity in the density regime where we still can equilibrate well with MD, namely 0.42 σ −3 ≤ ρ ≤ 0.57 σ −3 . To mimic a grandcanonical ensemble in our MD simulations where the total numbers of A and B chains are strictly fixed, we study fluctuations of X B , as described by C(X B , P), in small subsystems of the total system. In any such subsystem, the volume fraction is not conserved: Because the remainder of the system acts like a reservoir on the subsystem, one essentially realizes a grandcanonical ensemble for the subsystem. The feasibility of such an approach has been demonstrated previously for both the Lennard-Jones fluid at constant density [63] and the Ising/lattice gas model [64,65].
For this subsystem analysis, it is most convenient to work in the N VT ensemble and choose a cubic simulation box with L x = L y = L z . For an Ising system, it would be adequate to choose subsystems of the same shape, i.e., a cubic volume × × with L z . However, in the present case of nematically ordered systems, we must be aware of an extremely strong anisotropy: Since the A and B chains are stretched out over a length of the order of L ≈ (N − 1) b along the z-axis, all volume fraction fluctuations in the z-direction are strongly correlated; on a mean-field level, this effect was already described via the distinction of the correlation lengths ξ and ξ ⊥ (see Equations (14)- (17)). To test whether there is a tendency in favor of phase separation in the lateral directions, x and y, we found it advantageous to study quasi-two dimensional subsystems perpendicular to the z-axis.
The resulting fluctuations C(X B , P) = G coll (X B , 0) are plotted in Figure 8a versus 1/ , since one expects a (1/ )-finite size correction due to the subsystem boundary [63][64][65]. The data are compatible with this expectation but also show that there is only a very modest increase of G coll (X B = 0.5, 0) with increasing density (Figure 8a). This is also true for other choices of X B (Figure 8b), confirming our conclusion that the blend κ A = 24, κ B = 128 does not undergo critical n 1 -n 2 unmixing at densities where one still has the nematic phase. Next, we studied systems with slightly more flexible A chains (κ A = 20, κ B = 128), which, according to our DFT calculations (cf. Figure 3a), no longer have a mixed nematic phase extending over the full range of investigated pressures (0.27 k B Tσ −3 < P < 0.75 k B Tσ −3 ) and densities (0.4 σ −3 < ρ < 0.6 σ −3 ). Indeed, now the fluctuation near X B = 0.5 increases distinctly stronger with increasing density (Figure 8c) than for the case κ A = 24. However, it turns out to be extremely difficult to reach a decent precision of these data; hence, we could still not locate the critical point of n 1 -n 2 unmixing of the system with this method. One problem comes from the fact that critical slowing down is known to lead to a systematic bias (underestimation) of G coll (q = 0) when the runs are not long enough [66]. This problem is particularly severe for critical binary mixtures, where a finite size scaling of the relaxation time for interdiffusion τ int ∝ 4 is predicted [67]. However, we found it possible to study n 1 -n 2 unmixing for κ A = 20, κ B = 128, using the N VT ensemble and elongated boxes with L y = 32 σ, L z = 64 σ, N = 196,608, and a choice of L x such that the desired density was realized (e.g., L x = 192 σ for ρ = 0.5 σ −3 ). The time step here was chosen as 0.0025 τ MD , and runs of 2.5 × 10 9 MD steps were carried out for compositions X B = 1/3, 1/2 and 2/3. Starting again the simulations with chains oriented along the z-axis and initially separated in pure A and B domains (according to the chosen X B ), equilibrium was reached for κ A ≤ 20 for the densities of interest. We then determined the phase diagrams of Figures 9 and 10 by analyzing the resulting density profiles ρ A (x) and ρ B (x). From Figure 9, it is obvious that the topology of the phase diagram found in MD disagrees with its DFT counterpart for κ A = 20 (Figure 3a), but rather resembles the DFT phase diagram for κ A = 24 (Figure 3c). Likewise, the MD phase diagram for κ A = 16 ( Figure 10) is qualitatively similar to the DFT phase diagram found for κ A = 20 (Figure 3a). Thus, the massive discrepancies between DFT and MD phase diagrams for choices of κ A = 24 and κ A = 20 should not be interpreted as a complete breakdown of the DFT method: Rather, the latter is only somewhat inaccurate with respect to the prediction of the value κ M A , where the phase diagram topology changes. For κ A < κ M A , the phase diagram exhibits a triple point, while, for κ A > κ M A , two separate coexistence regions occur, i.e., an i-n region at lower pressures and a n 1 -n 2 coexistence region (ending in a critical point) at higher pressures. DFT predicts κ M A = 20.5 for κ B = 128, while MD rather suggests κ M A = 18 ± 1. Apart from this quantitative mismatch in κ M A , the general sequence of phase diagram changes with increasing κ A predicted by MD and DFT is the same. Further, these findings are also compatible with the predictions obtained by Semenov and Subbotin [21] for the limiting case of very large contour lengths. Hence, we conclude that the general features of the phase behavior predicted here are rather robust.
An important aspect of phase coexistence in the studied systems is that only the pressure P in the coexisting phases is identical, while the density ρ is not, as evidenced from the nonzero slopes of the tie lines drawn in Figures 9b and 10b. This density inhomogeneity could also be seen very well in the density profiles ρ(x) in the direction perpendicular to the interfaces: Due to the periodic boundary conditions, there must be two interfaces, and the slope of the tie lines tells us that the density of the phase n 2 is larger than both the density of the phase n 1 and the isotropic phase. Due to this density variation at the i-n 2 interfaces, there could occur some excess density ("interfacial adsorption") associated with these interfaces, which might affect our analysis of the MD results (as a finite size effect of order 1/L x ). However, a more detailed analysis of our results has shown that this effect (and other finite size effects) is negligible. Rather, we could show that our results on phase coexistence are fully compatible with consequences of the Gibbs phase rule. In fact, when the two A-rich (a) and B-rich (b) phases coexist in our simulation volume, the volume consists of two parts, V = V a + V b , since, per definition, interfaces have zero volume. The particle numbers N A and N B of the two kinds of monomers can be likewise split into particle numbers in the two phases with X B = N B /(N A + N B ). The phase boundaries in Figures 9 and 10 were based on the four partial densities ρ Aa = N Aa /V a , ρ Ab = N Ab /V b , ρ Ba = N Ba /V a , ρ Bb = N Bb /V b , namely ρ a = ρ Aa + ρ Ba , ρ b = ρ Ab + ρ Bb , and X Ba = ρ Ba /ρ a = N Ba /N a , as well as Note that the volume fraction V b /V must not be confused with the mole fraction X B = N B /N . Due to Equations (20) and (21) and the relation for V, the four partial densities are not all independent of each other. We have found it convenient to formulate this dependence via two equations for the volume fraction of phase b, V b /V = r 1 = r 2 , with ratios r 1 and r 2 defined as with partial densities ρ B = N B /V and ρ A = N A /V. Tables A1 and A2 quote the partial densities for X B = 0.5, from which Figures 9 and 10 were constructed, as well as the estimates r 1 and r 2 of the corresponding volume fraction of the B-rich phase. Results for X B = 1/3 and 2/3 were compatible with these results and are therefore not shown. We now discuss the nematic bond order parameters of the two types of chains, S A and S B , respectively, extracted as usual [35,39] as the largest eigenvalue of the tensor Q αβ (Equation (3)). In the pure A and B systems, a discontinuous transition occurs at P ≈ 0.38 k B Tσ −3 and P ≈ 0.07 k B Tσ −3 , respectively, where the nematic order parameter jumps from zero to about S A ≈ 0.52 and S B ≈ 0.76 [9,39], respectively (Figure 11a). In contrast, the mixed systems exhibit a more gradual behavior, as shown in Figure 11a. As expected from the phase diagram (Figure 10a), S B (X B = 0.5, P) starts to rise steeply as soon as the i-n 2 two-phase coexistence region has entered around P ≈ 0.10 k B Tσ −3 . Unlike the pure B system, the increase of S B (X B = 0.5, P) is continuous, since the mole fraction of the nematic phase grows continuously from zero to one as the two-phase coexistence region is crossed. In addition, S A (X B = 0.5, P) starts to become nonzero together with S B (X B = 0.5, P), due to A chains that are dissolved in the nematic B-rich phase (at a very small mole fractions, cf. Figure 10a), and which align due to the nematic ordering field exerted by surrounding majority of B chains. At P ≈ 0.3 k B Tσ −3 , the rise of S A (X B = 0.5, P) is already rather steep: We interpret this phenomenon as a "capillary nematization" effect of the nematic B-rich domains on the remaining A-rich domains, which, in a truly macroscopic system, would still be in the isotropic phase up to a pressure of about P ≈ 0.35 k B Tσ −3 , according to the phase diagram shown in Figure 10a. This finite size effect will be discussed in more detail below.
In bulk (macroscopic) systems, the gradual increase of the order parameters can be interpreted in terms of the lever rule Here, the two boundaries of the i-n two-phase coexistence region were denoted as X i B (P) and X n B (P). If we could cross this region at fixed P by varying X B , the variation of S A (or S B , respectively) with X B would be simply linear. However, this clearly is not true when we vary P at fixed X B (different from X B = 0 or 1, of course): Then, S A and S B are nontrivial functions, reflecting both the variation of the phase boundaries X i B (P) and X n B (P), as well as of S A,B (X n B (P)). The behavior described in Figure 11a and Equation (22) would be observed experimentally by methods that average over all domains in the system, e.g., scattering, optical birefringence, etc. In the simulations, it is possible to resolve the order parameters separately in the two kinds of domains: Then, S B (P) measured in the nematic B-rich domain rises much faster with P (Figure 11b), similar to the result of the pure B phase. Likewise, S A (P) measured in the isotropic A-rich domains stays zero up to P ≈ 0.35 k B Tσ −3 , where also the A-rich phase becomes nematic (Figure 11b). An analogous interpretation explains the behavior of the end-to-end distance measured separately in the coexisting domains ( Figure 11c). Equation (23) applies only for pressures P < P t because, for larger pressures, we enter a different two-phase coexistence region between two different nematic phases. The generalization to this case is  Figure 9), and, for about 0.2 k B Tσ −3 < P < 0.7 k B Tσ −3 , we have a uniformly mixed nematic phase, while, for about P > 0.7 k B Tσ −3 n 1 -n 2 , unmixing has occurred.
In the simulation, the order parameter of the less stiff chains increases indeed steeply around P ≈ 0.35 k B Tσ −3 : According to theory, there should be a clear jump of S A (X B , P) at P t , since S A (X i B (P < P t )) ≡ 0 in the isotropic phase, while S A (X n 1 B (P ≥ P t )) > 0 in the phase n 1 . There is also a discontinuous increase from the composition of the isotropic phase (X i B (P t )) to the composition X n 1 B (P t ) of the A-rich nematic phase. The resulting singular behavior of S A (X B , P) and S B (X B , P) predicted by Equations (23) and (24) is somewhat smeared out in the simulations, presumably due to finite size effects. In addition, the polymer end-to-end distances (Figure 11c) do not show related singularities either. Thus, it clearly would be desirable to study the system using much larger N and much better statistics, which was, however, computationally infeasible for us. One reason for unexpectedly large finite size effects is recognized when we examine the order parameters S i A , S i B , S n 1 A , S n 1 B , S n 2 A , and S n 2 B that belong to the various phases: We find that, in the case where we have i-n 2 two-phase coexistence in the simulation box, the order parameter S i B grows gradually with P and is of order 0.5 already when P = P t is reached. We interpret this finding as a "capillary nematization" effect of the B chains, which have S n 2 B ≈ 0.9 for P near P t already in the n 2 phase, in the isotropic phase adjacent to the i-n 2 interfaces. As P → P t , a kind of nematic wetting layer grows in the isotropic phase at the i-n 2 interface. If the linear dimensions L x → ∞, this effect will become negligibly small, but it is not for the L x values accessible to us. Another interesting feature is that, in the mixed nematic phase, the order parameter of the less stiff phase (S A ) always exceeds the corresponding value of the pure A phase at the same pressure, and likewise for the stiffer component (S B ) it is slightly smaller than for its pure counterpart. Thus, the more ordered stiffer chains in the direct neighborhood of less stiff chains enhance the order of the latter; conversely, the less stiff chains somewhat perturb the order of the stiffer ones, but the local nematic order of the mixture is very nonuniform since S A (X B = 0.5) is distinctly smaller than S B (X B = 0.5). A similar effect has already been noted for mixtures of chains which have the same stiffness but differ in chain length (where the longer chains are more ordered), as in Figure 6b,d. Figure 11d shows a counterpart to Figure 11b for the case κ A = 20, where a uniformly mixed nematic phase occurs. While we expect a discontinuous transition from S A = 0 in the isotropic phase to a nonzero value in the uniformly mixed nematic phase at the transition from the i-n two-phase region, a continuous splitting of both S A and S B into their different values in the coexisting nematic phase is expected at (X c B , P c ). Indeed, the data are compatible with this prediction. We see that S Ba in the A-rich nematic phase grows gradually from zero to distinctly nonzero values, when the transition to the nematic phase is approached: Again, we interpret this precursor effect as a finite size effect, due to "capillary nematization" at the i-n interfaces of the domain of the isotropic phase which has a finite extent in x-direction (according to Table A1, the volume fraction taken by the isotropic phase is about 0.39 for ρ = 0.35 σ −3 (P ≈ 0.172 k B Tσ −3 )).
Finally, we make contact between the n 1 -n 2 critical point found here (Figures 9 and 11d) and the phenomenological theory of Section 2.3: Can one understand the χ parameter postulated there from the stiffness asymmetry? To answer this question, we apply the approach of Kozuch et al. [42], who had considered only isotropic mixtures of rather flexible long polymers with slightly different stiffnesses.
Noting that Equation (13) also follows from the FH free energy with volume fractions φ A + φ B = 1 [11,12,15], one can see that the excess term χ eff φ A φ B relative to the entropy of mixing terms of the ideal mixture can only result from the difference in bending energies U A bend (θ ijk ) and U B bend (θ ijk ), since both types of chains are identical in all other respects for N A = N B . Since the average energy E = ∂ βF /∂β (with β ≡ 1/(k B T)), one can apply thermodynamic integration methods, based on the appropriate difference in bending energy. Denoting the total free energy of the mixture with X A = X B = 0.5 as F AB (κ A , κ B ), the sought after excess free energy is Equation (26) leads to since the number of bonds in the B chains is only one half of the total number of bonds per chain in the mixed system. Hence, Here, U B bend is the bending energy per monomer of a B chain, . . . AB denotes the thermal average in a 1:1 AB mixture, while . . . B refers to an average in a pure system of B chains. The desired free energy difference then is We have performed this thermodynamic integration for a system at ρ = 0.42 σ −3 in a cubic box of size L x = L y = L z = 64 σ with N = 111,232 monomeric units. We started the numerical integration in a state κ A = κ B = 128, and then lowered the stiffness of the A component to smaller values (Figure 12). One sees that the integrand is extremely small if both chains have similar stiffness, but it rises steeply for large stiffness mismatches. The resulting effective χ parameter is shown in the inset of Figure 12. Criticality (according to mean-field theory) would be reached for κ B /κ A ≈ 6, i.e., κ c A ≈ 21. However, while the pressure P ≈ 0.29 k B Tσ −3 for ρ = 0.42 σ −3 and κ A = 16 clearly falls inside the i-n 2 coexistence region, for κ A = 20, the critical pressure P c ≈ 0.7 k B Tσ −3 corresponds to ρ c = 0.58 σ −3 ; hence, this system at ρ = 0.42 σ −3 is in the mixed nematic region. Since χ eff can be expected to increase with increasing density, we cannot estimate χ eff for the data shown in Figure 10, but the order of magnitude of χ eff expected on the basis of Figure 12

Conclusions
In this paper, we have considered binary mixtures of semiflexible polymers in a common good solvent. We have studied the case where both constituents have the same stiffness but differ in chain length and the case where both constituents have strictly the same chain length but have a large stiffness disparity. We have considered here only very simple models appropriate for lyotropic solutions, where nematic order is purely entropydriven, and applied two complementary computational approaches: (i) Density Functional Theory, based on a tangent hard sphere model of chains where stiffness is controlled by a bond angle potential (Equation (1)), and (ii) Molecular Dynamics simulation, based on standard bead-spring models augmented with the same bond angle potential. We considered binary mixtures of chains with N A = 16 and N B = 32 monomeric units (of diameter σ and bond length between the units b ≈ σ) with persistence length p / b = 32, as well as mixtures with N A = 16 and N B = 64 with p / b = 128, both at mole fraction X B = 0.5. It was found that the i-n coexistence regions are similarly narrow as for the corresponding pure systems. In the mixed nematic phase, the nematic order parameter S A of the shorter chains is always less than S B but larger than S A for the corresponding pure system.
The second case focused on chains of the same length, N A = N B = 32, where the stiffer component is almost rod-like ( B p / b = 128), while the persistence length of the less stiff chain was much smaller (16 ≤ A p / b ≤ 26). These systems also have a strong tendency for nematic order, but unlike the former cases exhibit relatively wide i-n two-phase regions. Again, the chains in nematic phases adapt to their environment; the less stiff chains are rather strongly stretched when they are dissolved in a stiffer matrix, whereas the inverse effect is much less pronounced. The stiffness mismatch between the two constituents has been shown to cause an unmixing tendency, which can be described by an effective Flory-Huggins χ parameter, although there are no attractive forces between monomeric units of different chains present whatsoever. We have computed this effective χ parameter for one density as function of κ B /κ A (for κ B = 128) through thermodynamic integration, and have shown that mean-field theory based on this χ parameter would yield roughly the correct ratio κ B /κ A for which this phase separation is actually observed. Using the random phase approximation, one predicts that critical correlations are very anisotropic in dilute and semi-dilute solutions: The correlation length parallel to the nematic director scales proportional to chain length N, while it only is proportional to √ N in the perpendicular direction. A thorough numerical investigation of this new type of critical behavior is very challenging, however, and must be left to the future. There has been great recent interest in the critical behavior of anisotropic systems in the Ising universality class [68,69], since the principal directions of order parameter correlations beyond mean-field theory may deviate from the expectations based on mean-field theory, and one can show that then the principle of "two-scale factor universality" [70] is violated.
This critical n 1 -n 2 unmixing can be seen for an intermediate range of κ B /κ A only: When κ B /κ A is not large enough, the critical point would be at unphysically large monomer densities, where the nematic phase no longer exists at all (and smectic or crystal phases have taken over). On the other hand, when the ratio increases, a multi-critical point emerges (see Figure 3b), where the critical point touches the i-n coexistence region. For still larger κ B /κ A , it would be metastable, and the equilibrium phase diagram then exhibits a triple point (Figures 2a,b, 3a and 10a). For κ B = 128, DFT has predicted the multi-critical point where the phase diagram topology changes to occur for κ M A = 20.5, while MD, rather, suggests κ M A = 18 ± 1. We do not attribute this discrepancy to the minor differences between the chain models used by DFT and MD but, rather, hold responsible the approximations used by DFT, where only the orientation ω of a rod-like effective chain enters as local freedom, local density is not a variable, and the effective pairwise interaction V excl (ω, ω ) used in Equation (10) is approximate. Nevertheless, the guidance provided by DFT to yield an overall picture of phase behavior and lead the interpretation of MD results is extremely valuable.
Unfortunately, we are not aware of real systems to which the present model calculations could be compared directly, although there have been studies of mixtures of two different semiflexible polymers in a common solvent. For example, Russo and Cao [44] obtained the phase diagram of poly(-γ-benzyl-α, L-glutamate) (PBLG) and Nylon 6 in m-cresol as a solvent. PBLG in various solvents has broadly been used as a model material of a liquid-crystal forming polymer in lyotropic solutions [1,2]. However, for the system of Reference [44], the case N A = N B = N studied here has been out of focus, and, furthermore, the relevance of hydrogen-bond interactions was stressed, which makes the problem much more complicated.
Funding: We thank the German Research Foundation (DFG) for support under Project Numbers BI 314/24-2, NI 1487/4-2, and NI 1487/2-2. S.A.E. thanks the Alexander von Humboldt Foundation for support. The authors gratefully acknowledge the computing time granted on the supercomputer Mogon (hpc.uni-mainz.de). A.M. thanks the COST Action No CA17139 supported by COST (European Cooperation in Science and Technology) and its Bulgarian partner FNI/MON under KOST-11.

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

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Table A1. Partial densities and ratios r 1 and r 2 for a system with N A = N B = 32, κ A = 20, and κ B = 128.