Bottlebrushes are an emerging class of polymers, characterized by a high density of side chains grafted to a linear backbone that offer promise in creating materials with unusual combinations of mechanical, chemical, and optoelectronic properties. Understanding the role of molecular architecture in the organization and assembly of bottlebrushes is of fundamental importance in polymer physics, but also enabling in applications. Here, we apply field-theoretic simulations to study the effect of grafting density, backbone length, and side-chain (SC) length on the structure and thermodynamics of bottlebrush homopolymer melts. Our results provide evidence for a phase transition from an isotropic to a nematic state with increasing grafting density and side-chain length. Variation in the backbone length is also observed to influence the location of the transition, primarily for short polymers just above the star to bottlebrush transition.

## I. INTRODUCTION

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.^{9} 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.*^{14} and Fredrickson^{9} 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.*^{15} and experiments on dilute solutions by Rathgeber^{16} 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^{17} 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.*^{18} 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.*^{24} 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*_{1}, where *z*_{1} = 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*_{1} < *z* < *z*^{*}, where $z*=NSC\u22121/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^{25} (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.

## II. MOLECULAR MODEL

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.1*N*_{BB} to 0.9*N*_{BB}, creating *n*_{SC} + 1 subsegments in the backbone with the distance between grafts *d* = 0.8*N*_{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 *ρ*_{0} = *n*_{BB}(*N*_{BB} + *n*_{SC}*N*_{SC})/*V*, where *V* is the system volume, and the average volume of a segment is $v0=\rho 0\u22121$.

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

where the functional integrals $\u222b\u22c6Dri$ run over the possible space curves **r**_{i}(*s*) for the *i*th 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.1*N*_{BB} + (*k* − 1)*d* are the grafting contour positions for side chains $k\u22081,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 $\Gamma (r)=2\pi a2\u22123/2e\u2212r2/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 $\u222bDw$ is a functional integral over real-valued configurations of an auxiliary field *w*(**r**). By applying Eq. (6) with $J(r)=\Gamma \u22c6\rho ^r\u2212\rho 0$ and $A\u22121(r,r\u2032)=\zeta \rho 0\delta (r\u2212r\u2032)$, we find

where *w*(**r**) is an auxiliary chemical potential field that is conjugate to $\Gamma \u22c6\rho ^$ 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=\u222bDrexp\u221214\u2211lnSC+1\u222b0Nlds\u2009drlds2$ is the ideal-chain partition function of a single bottlebrush polymer. The prefactor $Z0$ in Eq. (9) is related to *Z*_{0} 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*_{0}(**r**, *s*) that the chain contour *s* is at position **r** in the field free case.^{29} 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.1*N*_{BB} + (*j* − 1)*d*) is the branching point for the side chain *j* ∈ [1, *n*_{SC}], and *τ*_{0} = 0 and $\tau 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\xafBB(\u2009j)$ is the backward propagator for the segment from $s=NBB\u2212\tau nSC+1\u2212j$ to $s=NBB\u2212\tau nSC\u2212j$. 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 $\rho r$, the monomer density operator that enters the thermodynamic forces in complex Langevin sampling of the $wr$ fields.^{25}

### A. Nematic order parameter

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=\u2207wr\u2207wr$. Then, a rank-two local nematic tensor is

A global nematic order parameter can be obtained by averaging $S\u0303(r)$ over both space and field configurations so that $Sij=l3V\u222bdr\u2009\u27e8S\u0303ij(r)\u27e9$.^{31} 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*^{3}, 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\u0303$, 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 *λ*_{1}, *λ*_{2}, and *λ*_{3} are the eigenvalues of $S\u0303\u2032$. These quantities can individually signal the onset of backbone stiffening and the presence of the nematic order.

## III. SAMPLING THE FIELD THEORY

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.^{25} 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.

### A. The mean-field approximation

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, $\delta H[w]\delta w(r)w=w*=0$. For our model,

and the saddle-point field configuration is $iw\u22c6(r)=\zeta \rho 0\u22121\Gamma \u22c6\rho $$(r)\u22121$. Restricting to the case of periodic boundary conditions, the only physical saddle-point solution is homogeneous with $\rho r=\rho 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.

### B. Beyond mean-field theory: Complex Langevin sampling

Our complex Langevin simulations use a cubic cell of side 64*l* 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* = 2*l*. Unless noted otherwise, we use a backbone length of *N*_{BB} = 100 and the overall density of segments is fixed at *ρ*_{0} = 0.7*l*^{−3} throughout. The so-called invariant degree of polymerization is $N\xaf\u2261Ntotb6/v02=63(0.7)2(NBB+nSCNSC)$, where *N*_{tot} is the total number of segments in one bottlebrush polymer, and $v0=\rho 0\u22121=l3/0.7$. As noted previously, our bottlebrush architecture has grafted side chains between 0.1*N*_{BB} and 0.9*N*_{BB} with a graft spacing of $d=0.8NBBnSC\u22121$ and a grafting density of *z* = *d*^{−1}. The complex Langevin equations are propagated in time using the Petersen-Öttinger predictor-corrector algorithm^{32} 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.

#### 1. Results: Variation of bottlebrush grafting density

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\u0303ij$ 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\u0303$ 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\u0303\u2032$, 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\u2212pzi\u22121,NSC)/pzi\u22121,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).

#### 2. Results: Fixed grafting density

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\u0303$ 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.^{4} 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.

## IV. CONCLUSIONS

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.

## ACKNOWLEDGMENTS

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.