Theoretical prediction of an isotropic to nematic phase transition in bottlebrush homopolymer melts

Bottlebrush polymers—macromolecules with polymer side chains densely grafted to a polymeric backbone—can be prepared at extraordinarily high grafting densities using grafting-through polymerization of macromonomers with ring-opening metathesis polymerization (ROMP) chemistry (see, e.g., Refs. 1–4 and references therein). Since the advent of grafting-through synthetic approaches, which have greatly expanded the access to such materials, there has been a surge of interest in the properties of bottlebrush polymers. Bottlebrush polymers have been created with a wide range of composition and architecture variations, including those with one type of arm, bottlebrush random copolymers with random placement of two or more types of arms, or block or tapered bottlebrush copolymers. As the number of side chains and/or the length of side chains increases, the viscoelastic properties of a melt of such polymers change dramatically due to dilution of the backbone entanglements, which can lead to a material with a very low modulus, high strain at break, and extraordinary elasticity.^5–8 The bottlebrush architecture is remarkable in the number of handles it provides for manipulating material properties with or without changing the chemical composition.

There is an extensive body of research, both experimental and theoretical, studying the conformations of bottlebrush polymers in dilute solutions—effectively a single chain. Several studies have focused on understanding how the persistence length of bottlebrush polymer backbones varies with grafting density and molecular weight of the side-chains(SCs) and have suggested a transition from coil-like to rod-like backbone configurations^6,9–11 with a possible lyotropic behavior.⁹ Grafting-induced rigidification of the backbone is also implicated in bottlebrush block polymers, which predominantly self-assemble into lamellar phases with unusual scaling of the domain-spacing with the molecular weight.^4,7,12,13 Theoretically, Birshtein et al.¹⁴ and Fredrickson⁹ predicted that both the diameter D and persistence length, λ, of a densely grafted bottlebrush polymer in a dilute solution should increase with grafting density and side-chain length. Birshtein predicted that the ratio, λ/D, is approximately independent of side-chain length, while Fredrickson suggested a near linear increase, opening the door to lyotropic nematic behavior in solutions at higher concentrations. However, subsequent computer simulations of single chains by Saariaho et al.¹⁵ and experiments on dilute solutions by Rathgeber¹⁶ support the Birshtein picture that λ/D is approximately independent of side-chain length for densely grafted bottlebrushes and cannot be raised to a level (>10) expected to drive a nematic order at higher concentration in conventional semiflexible polymers.

On the experimental side, Schmidt¹⁷ as early as 1995 found evidence of the lyotropic nematic order in solutions of poly(methylacrylate)-graft-polystyrene in toluene at approximately 30% by weight, suggesting that further conformational stiffening of bottlebrushes (or contraction in diameter) arises from intermolecular excluded volume interactions as bottlebrushes are concentrated into the semidilute regime and beyond. More recent studies by Storm et al.¹⁸ observed the nematic order in a lyotropic bottlebrush system of DNA-based semiflexible chains and speculated that the aspect ratio λ/D will increase only at very high grafting densities of very long side chains. However, theory and simulation have addressed this ratio only at infinite dilution, whereas the forced interpenetration of bottlebrushes at a finite concentration will influence both the numerator and denominator. Given multiple experimental observations of the nematic order, we anticipate that interpenetration leads to an increase in λ/D over the infinite dilution value, possibly driven by a decrease in the effective diameter D. Clearly, the presence of lyotropic liquid crystalline behavior in bottlebrush solutions and the factors and precise criteria governing its emergence remain unresolved in the literature.

The case of bottlebrush polymer melts has received considerable attention by theorists and experimentalists in terms of chain conformations and rheology,^{5,6,8,16,19–23} but there are no reported attempts to locate and map out an isotropic-nematic phase transition. One recent report by Borodinov et al.²⁴ provided helium ion microscopy evidence of nematic order in an annealed, thin film bottlebrush polymer melt, supported by limited coarse-grained molecular dynamics (MD) simulations. The isotropic-nematic boundary was not accessed using either technique.

In this paper, we use field-theoretic simulations (FTSs) to address the effects of architecture variation on liquid crystalline order in bulk bottlebrush homopolymer melts. Reference 6 introduced four architecture classification regimes that were correlated with the chain conformations and the viscoelastic behavior of a melt: loosely grafted combs, densely grafted combs, loosely grafted bottlebrushes, and densely grafted bottlebrushes. Loosely grafted combs, characterized by grafting densities z < z₁, where z₁ = 1/N_SC, have grafts separated on the backbone by more than the graft/side-chain (SC) contour length (N_SC), while densely grafted combs satisfy z₁ < z < z^*, where $z*=NSC−1/2$⁠, having grafts separated by less than the SC contour length, but more than the SC diameter of gyration. Loosely grafted bottlebrushes correspond to z > z^* and have grafts spaced by less than the SC diameter of gyration, while densely grafted bottlebrushes have grafts on every backbone monomer. We stress that in this work, we do not study the densely grafted bottlebrush regime. Here, we address the transition from loosely grafted combs to densely grafted combs to loosely grafted bottle brushes upon increasing grafting density or SC length at a fixed backbone length (Fig. 1). We also study the transition from star-like conformations to comb and bottlebrush conformations as the backbone length is increased at a fixed grafting density for various side-chain lengths. As we already observed nematic order in the loosely grafted bottlebrush regime, we do not examine the case of densely grafted bottlebrushes.

Both atomistic and coarse-grained particle simulations of bottlebrush melts are extremely challenging due to the large molecular weight, the high density, and the long relaxation times involved. In contrast, field-theoretic simulations²⁵ (FTSs) are well suited to bottlebrush melts as their computational efficiency relative to particle based methods improves with both density and molecular weight. In this study, we use FTSs to show for the first time that variation of side-chain grafting density, side-chain length, and backbone length spanning the regimes illustrated in Fig. 1 can produce backbone stiffening in loosely grafted bottlebrushes sufficient to induce an isotropic-nematic phase transition. It is notable that in the chain model employed, both backbone and side-chain segments are inherently flexible with no resistance to bending.

This paper is organized as follows: Section II describes the coarse-grained particle model that we use to represent bottlebrush polymer melts and its equivalent representation as a statistical field theory. Section II A defines two nematic order parameters, constructed from gradients of the auxiliary field conjugate to monomer density and from the distribution of backbone end-to-end vectors. Section III A describes the mean-field self-consistent field theory (SCFT) approximation to our model, while Sec. III B outlines the procedure and parameters utilized for unapproximated field-theoretic simulations based on complex Langevin sampling and presents our results. Section IV provides a brief summary of our conclusions.

We consider a melt of n_BB bottlebrush homopolymers with backbone length N_BB, each containing n_SC side chains of length N_SC. The side chains are uniformly distributed along the backbone chain from 0.1N_BB to 0.9N_BB, creating n_SC + 1 subsegments in the backbone with the distance between grafts d = 0.8N_BB/(n_SC − 1). We assume here that all strands of the polymer, including the backbone, are intrinsically flexible and adopt continuous Gaussian chain statistics.

In practice, we describe the ensemble of bottlebrush polymers as a collection of n = n_BB(n_SC + 1) separate continuous Gaussian strands with branch-point tethering constraints imposed in the partition function. A microscopic operator for the density of segments is

where r_i(s) is the space curve of strand i, parameterized as a continuous linear filament with contour variable s, and the s integral runs [0, N_BB] for backbone strands and [0, N_SC] for side-chain strands. In an nVT ensemble, the overall segment density in the melt is ρ₀ = n_BB(N_BB + n_SCN_SC)/V, where V is the system volume, and the average volume of a segment is $v0=ρ0−1$⁠.

We begin with the following canonical partition function for a coarse-grained model of a polymer melt of bottlebrush homopolymers:

where the functional integrals $∫⋆Dri$ run over the possible space curves r_i(s) for the ith polymer strand, subject to an architecture-specific tethering constraint,

Here, $rjBB(s)$ are the space curves for the subset of strands that form bottlebrush polymer backbones, τ_k = 0.1N_BB + (k − 1)d are the grafting contour positions for side chains $k∈1,nSC$⁠, and $rkSC,j(s)$ are the space curves of side chains of polymer j. The total Gaussian stretching energy for all strands is

where b is the length of a statistical segment, which is here assumed to be uniform throughout each polymer. A Helfand compressibility penalty for the melt is imposed by

where $Γ(r)=2πa2−3/2e−r2/2a2$ is a normalized Gaussian smearing filter of width a that serves to regularize the field theory,^26–28 ⋆ is a spatial convolution, and ζ sets the compressibility of the melt.

The purely quadratic form of Eq. (5) permits a direct decoupling of interactions through a Hubbard-Stratonovich transformation given by

where $∫Dw$ is a functional integral over real-valued configurations of an auxiliary field w(r). By applying Eq. (6) with $J(r)=Γ⋆ρ^r−ρ0$ and $A−1(r,r′)=ζρ0δ(r−r′)$⁠, we find

where w(r) is an auxiliary chemical potential field that is conjugate to $Γ⋆ρ^$ and D_w is the Gaussian functional integral denominator from Eq. (6) that is a multiplicative contribution to the partition function.

Substituting Eq. (7), the canonical partition function becomes

where we have isolated the terms that depend on the polymer coordinates in the second row. By expressing spatial coordinates in units of reference length $l=b/6$⁠, the partition function becomes

where

is the Hamiltonian functional of the field theory. The normalized partition function of a single bottlebrush polymer in field w(r) is

where $Z0=∫Drexp−14∑lnSC+1∫0Nldsdrlds2$ is the ideal-chain partition function of a single bottlebrush polymer. The prefactor $Z0$ in Eq. (9) is related to Z₀ by the expression $Z0=Z0nBB/(DwnBB!)$⁠.

An efficient numerical approach for computing the single chain partition function Q of a bottlebrush polymer was reported in the supplementary information of Ref. 4. We reproduce the key expressions here. The approach relies on the definition of a chain propagator, q(r, s; [w]), which represents the statistical weight for a chain of length s to have its end at position r. This object is strictly analogous to the probability density p₀(r, s) that the chain contour s is at position r in the field free case.²⁹ Using propagators, the single partition function Q has an important factorization property that the statistical weight of a linear chain of length N can be built by joining the propagator for the end distribution of a chain with length s and the propagator of a “complementary” chain of length N − s. For our continuous Gaussian chain bottlebrush model, the forward propagators of the side chains satisfy the modified diffusion equation

with initial condition q_SC(r, 0) = 1 and s ranges from 0 to N_SC. The forward propagator of the backbone is obtained by dividing the backbone into n_SC + 1 sections with contour variable ranges τ_j ≤ s < τ_j+1, where τ_j = (0.1N_BB + (j − 1)d) is the branching point for the side chain j ∈ [1, n_SC], and τ₀ = 0 and $τnSC+1=NBB$ are the two ends of the backbone strand (see Fig. 2). Each section satisfies the diffusion equation

with initial conditions

The backward propagator of the backbone is also obtained by dividing the strand into n_SC + 1 segments, where $q¯BB(j)$ is the backward propagator for the segment from $s=NBB−τnSC+1−j$ to $s=NBB−τnSC−j$⁠. Each segment satisfies the diffusion equation of the forward backbone propagator [Eq. (13)] with the following initial conditions:

The sum of the backward propagators of the side chains can be solved by accumulating the initial conditions of all side-chain strands from the backbone branch points

and solving the side-chain diffusion equation [Eq. (12)] with this initial condition. We note that in this computational scheme, both the forward and reverse propagators for side chains are solved in a single pass regardless of n_SC. The computational effort is nearly ideal, scaling linearly with the backbone length at a fixed contour resolution, and is approximately independent of the number of side chains.

Finally, we can construct the single chain partition function and its first functional derivative,

which leads to an expression for $ρr$⁠, the monomer density operator that enters the thermodynamic forces in complex Langevin sampling of the $wr$ fields.²⁵ 

In field-theoretic simulations, the details of individual polymer chain conformations are eliminated in order to sample equilibrium configurations efficiently. However, since the auxiliary field w is a chemical-potential field conjugate to the microscopic segment density, one should expect any emergent orientation texture of the polymer strands to show a signature in w(r). Although methods have been reported for quantifying the chain orientation directly from the propagators (e.g., Ref. 30), here, we apply a simpler nematic order parameter evaluated directly from the w(r) field. Given a field configuration $wr$⁠, we define a local orientation unit vector as $nr=∇wr∇wr$⁠. Then, a rank-two local nematic tensor is

A global nematic order parameter can be obtained by averaging $S̃(r)$ over both space and field configurations so that $Sij=l3V∫dr ⟨S̃ij(r)⟩$⁠.³¹ The maximum eigenvalue of S_ij defines a nematic order parameter and the corresponding eigenvector is the orientation. However, in a large simulation cell, the fluid may have a nematic texture with multiple domains, each with a randomly oriented director, and the process of spatially averaging the nematic tensor over all such domains can eliminate the signature of the nematic order. In order to detect ordering, we therefore average the maximum eigenvalue instead of the tensor itself; we divide the simulation domain into cubic subcells of side c and volume c³, compute the maximum eigenvalue of S_ij for each subcell, and average the maximum absolute eigenvalue over all subcells contained in the simulation cell.

Note that the auxiliary field w is coupled in an unbiased way to both the backbone and side chains, which might be expected to make it difficult to detect the backbone nematic order. To directly validate our interpretation of $S̃$⁠, we also examine the orientation of the polymer backbone distribution. We study the distribution of the end segment of a polymer chain that has its first segment constrained to the center of the simulation cell, i.e., a chain propagator is computed with a delta-function initial condition for the chain start. Using this chain-end propagator distribution, we can compute the average length of the end-to-end vector, the asphericity of the distribution, and a propagator nematic tensor defined as

We define the asphericity of the distribution as

where λ₁, λ₂, and λ₃ are the eigenvalues of $S̃′$⁠. These quantities can individually signal the onset of backbone stiffening and the presence of the nematic order.

With the molecular model and field-theoretic partition function defined, we next discuss two methods of analyzing the field theory: self-consistent field theory (SCFT), which invokes a mean-field approximation, and approximation-free complex-Langevin sampling of the full partition function.²⁵ The latter approach stochastically generates importance-sampled w field configurations consistent with the statistical weight specified in Eq. (9) while tolerating the sign problem inherent to complex-valued field theories.

The self-consistent mean-field approximation assumes that the partition function of Eq. (9) is dominated by a single field configuration, w^⋆(r), for which the Hamiltonian functional is stationary, $δH[w]δw(r)w=w*=0$⁠. For our model,

and the saddle-point field configuration is $iw⋆(r)=ζρ0−1Γ⋆ρ$$(r)−1$⁠. Restricting to the case of periodic boundary conditions, the only physical saddle-point solution is homogeneous with $ρr=ρ0$ and w^⋆(r) = 0. Therefore, the mean-field approximation does not support the emergence of the nematic order, and we argue that field fluctuations are crucial to probe this phenomenon.

Our complex Langevin simulations use a cubic cell of side 64l with periodic boundary conditions. The w field is resolved on a 32 × 32 × 32 collocation mesh. The interactions between segments are specified by ζ = 1 and a = 2l. Unless noted otherwise, we use a backbone length of N_BB = 100 and the overall density of segments is fixed at ρ₀ = 0.7l⁻³ throughout. The so-called invariant degree of polymerization is $N¯≡Ntotb6/v02=63(0.7)2(NBB+nSCNSC)$⁠, where N_tot is the total number of segments in one bottlebrush polymer, and $v0=ρ0−1=l3/0.7$⁠. As noted previously, our bottlebrush architecture has grafted side chains between 0.1N_BB and 0.9N_BB with a graft spacing of $d=0.8NBBnSC−1$ and a grafting density of z = d⁻¹. The complex Langevin equations are propagated in time using the Petersen-Öttinger predictor-corrector algorithm³² with a time-step Δt in the range 0.001–0.005 chosen to provide a good balance between accuracy, stability, and efficiency. The propagator diffusion equations are solved pseudospectrally with operator splitting^33,34 using Δs = 0.01.

We begin with a fixed backbone length N_BB = 100 and vary the number of side chains from n_SC = 5 to 40, with steps of 5 (corresponding to a grafting density range of z = 0.05–0.4875), and the length of all side chains from N_SC = 5 to 50, with steps of 5.

This range of number and length of arms covers the transition from combs to loose bottlebrushes, as defined in Ref. 6. In the terminology used in Refs. 8 and 9, this range covers the transition between the low- and high-coverage limit, i.e., the limit when the diameter of gyration of the side chains becomes greater than the average backbone contour length between grafts, which implies that the chains start to interact, the side-chains stretch out normal to the polymer backbone to avoid contact with other chains, and the configuration of the bottlebrush changes from that of a coil toward a more extended rod.

Figure 3 shows the nematic order parameter derived from $S̃ij$ as a function of the grafting density and length of side chains. Our nematic order parameter is positive everywhere, but abruptly increases in value for sufficiently high grafting densities of long arms, signaling the onset of nematic order. Figure 4 shows the nematic order orientation as it is captured by the eigenvector corresponding to the largest eigenvalue of $S̃$ in a portion of the simulation box for a system in the isotropic and in the nematic phase, respectively. Indeed, for a system that our parameters characterize as nematic, we observe a multidomain nematic texture possessing a local orientational order within each domain. As further characterization of the shape and orientation of the bottlebrushes, we examine the chain-end propagator distribution for a bottlebrush initiated at the center of the simulation cell.

Figure 5 shows the probability amplitude q_BB(r, s = N_BB) for z = 0.05, N_SC = 5 and for z = 0.4875, N_SC = 50, which from Fig. 3 are expected to be in the isotropic and nematic phase, respectively. As expected, the propagator distribution becomes strongly anisotropic for systems that show the nematic order. In Fig. 6, we analyze the nematic order parameter from the chain-end distribution, $S̃′$⁠, and notice a decrease in the second and third eigenvalue and an increase in the largest eigenvalue. This is captured by the asphericity of the eigenvalues, which increases with increasing grafting density and length of side chains. In Fig. 7, we also demonstrate anisotropy from one more metric extracted from q_BB(r, N_BB): the average end-to-end vector of the backbone. Our results are in accordance with the theoretical predictions given in Ref. 6 for the end-to-end distance for loosely and densely grafted combs and loosely grafted bottlebrushes.

To place our predicted I-N boundary in context, in Fig. 8, we overlay the boundaries extracted from Figs. 3 and 6, as localized by the line of maximum variation of the order parameter, onto the grafting regimes defined in Ref. 6. More precisely, let $pzi,NSC$ denote the nematic order parameter at grafting density z_i and side-chain length N_SC. We define the relative variation of p as $(pzi,NSC−pzi−1,NSC)/pzi−1,NSC$⁠. The maximum variation of p at a fixed side-chain length, N_SC, occurs at the point of maximum relative variation over all grafting densities.

Our simulations range from loosely grafted combs to loosely grafted bottlebrushes (see Fig. 1). The isotropic-nematic phase boundaries mapped in Figs. 3 and 6 agree semiquantitatively with each other and lie fully within the loosely grafted bottlebrush regime. Our results therefore validate our interpretation of Eq. (20).

For this study, we imposed a fixed grafting density z = 0.4875 and varied the backbone length from N_BB = 38 to 158, in steps of 20. In each case, the sequence of side-arm grafting positions begins and ends 10 segments before the end of the backbone strand. We also varied the length of side chains from N_SC = 5 to 50, in steps of 5.

Figure 9 shows the nematic order parameter based on $S̃$ as a function of the backbone length and the length of side chains. Our results show that at a fixed grafting density, the location of the I-N phase transition depends on both the side-chain length and overall backbone length. In particular, for this range of side-chain lengths, we see an increase in our nematic order parameter occurring only when N_BB ≥ 58. We note that at a backbone length N_BB = 38 and a side-chain length N_SC > 35, the side-chain length is greater than the backbone length, N_SC > N_BB, suggesting a star-like conformation.⁴ Thus, we see an increase in our nematic order parameter occurring soon after entering the transition from star to bottlebrush conformations. We further observe that this nematic order parameter increase plateaus with increasing backbone length, as expected if the persistence length λ and the effective diameter D are independent of N_BB, but depend only on density, grafting density, and side-chain length.

Nematic order has been observed previously experimentally and supported computationally in Ref. 24 for dense bottlebrushes of a backbone length N_BB = 80 and a side-chain length N_SC = 10 with a grafting density of one side chain per backbone monomer. Such a system falls outside of the regimes considered here. Nonetheless, our study provides support for this observation and importantly shows that nematic ordering can occur even for loosely grafted bottlebrushes.

The structure and rheology of melts of bottlebrush macromolecules evidently depend sensitively on the side-chain grafting density, and the degree of polymerization of the side chains and of the backbone. Fully fluctuating field-theoretic simulations of bottlebrush homopolymer melts were carried out at a fixed backbone molecular weight, while varying the grafting density and the side-chain length, and at a fixed grafting density, while varying the backbone and the side-chain lengths. Our results, based on two different metrics for the orientational order, suggest that an isotropic to nematic transition can occur in large regions of the three parameter space, even for systems of fully flexible chains. We find that the isotropic to nematic transition is contained within the regime of “loosely grafted bottlebrushes,” which is readily accessible by modern synthetic polymer techniques. Our predictions should serve to guide the experimental exploration of liquid crystalline ordering in single component bottlebrush polymer melts.

E.P. was supported by NSF (Grant No. DMS-1913180). K.T.D. and G.H.F. were supported by NSF (Grant No. DMR-1822215). All simulations were performed using computational facilities of the Center for Scientific Computing at the CNSI and MRL: an NSF MRSEC (Grant No. DMR-1720256) and Grant No. NSF CNS-1725797.