Bottlebrush polymers are a class of highly branched macromolecules that show promise for applications such as self-assembled photonic materials and tunable elastomers. However, computational studies of bottlebrush polymer solutions and melts remain challenging due to the high computational cost involved in explicitly accounting for the presence of side chains. Here, we consider a coarse-grained molecular model of bottlebrush polymers where the side chains are modeled implicitly, with the aim of expediting simulations by accessing longer length and time scales. The key ingredients of this model are the size of a coarse-grained segment and a suitably coarse-grained interaction potential between the non-bonded segments. Prior studies have not focused on developing explicit forms of such potentials, instead, relying on scaling arguments to model non-bonded interactions. Here, we show how to systematically calculate an interaction potential between the coarse-grained segments of bottlebrush from finer grained explicit side chain models using Monte Carlo and Brownian dynamics and then incorporate it into an implicit side chain model. We compare the predictions from our coarse-grained implicit side chain model with those obtained from models with explicit side chains in terms of the potential of mean force, the osmotic second virial coefficient, and the interpenetration function, highlighting the range of applicability and limitations of the coarse-grained representation. Although presented in the context of homopolymer bottlebrushes in athermal solvents, our proposed method can be extended to other solvent conditions as well as to different monomer chemistries. We expect that our implicit side chain model will prove useful for accelerating large-scale simulations of bottlebrush solutions and assembly.

## I. INTRODUCTION

Bottlebrush polymers are a class of branched macromolecules, consisting of numerous short chains densely grafted onto a linear polymer. They are rapidly emerging as attractive building blocks for soft materials that can be widely used in many applications.^{1–9} Recent applications of this class of macromolecules include solvent-free soft elastomers,^{10–16} surface modification for improved lubrication and antifouling,^{17,18} photonic materials,^{19–23} pressure sensitive adhesives,^{24} self-healing polymeric networks,^{25} and micellization for targeted drug delivery.^{6,26} The wide adoption of the bottlebrush architecture stems from two distinct factors: (i) the presence of densely grafted side chains adds a non-trivial molecular thickness that reduces intermolecular entanglements and helps in rapid self-assembly^{27,28} and (ii) many choices are available for the molecular design, e.g., degree of polymerization of the backbone and side chains,^{29,30} grafting density,^{31} and monomer chemistry,^{32} which can be independently varied to tune material properties at the macroscopic level. Understanding the role of architectural parameters is crucial to improved designing of bottlebrush materials. In conjunction with experimental characterization and analytical theory, computational studies can provide key insights into the underlying physics^{33} and help us correlate bulk material properties with the molecular architecture.

Bottlebrush polymer molecules are most commonly modeled using a coarse-grained bead–spring representation,^{34–38} although several recent studies have used different models as well, such as density-functional theory,^{39} self-consistent field theory,^{40,41} and polymer reference interaction site models.^{35,42} In these bead–spring models, each bead of the backbone and the side chains of bottlebrushes both typically resolve down to the level of a Kuhn step. However, resolving the side chains results in very large degrees of freedom at the molecular level significantly drives up the computational cost and makes it difficult to access the large length scale properties characteristic of multi-chain bottlebrush materials. For instance, a bottlebrush with 500 Kuhn steps on the backbone and 30 Kuhn steps/side chain with a grafting density of unity will require 15 500 beads for a single molecule. Clearly, this approach becomes restrictive for non-dilute bottlebrush solutions and melts, where the total number of molecules is large. This motivates us to seek a molecular model in which the side chains can be coarse-grained away in an appropriate fashion—we will call such models *implicit side chain* (ISC) models, in contrast to the more common *explicit side chain* (ESC) models. A typical simulation snapshot and schematic of the ESC model and a conceptual picture of the ISC model are shown in Fig. 1.

In our previous works,^{43–45} we develop an ISC model for bottlebrushes based on a wormlike cylinder. This models a bottlebrush molecule using four descriptors: (i) a contour length *L*, (ii) a Kuhn length *λ*^{−1}, (iii) a diameter *d*, and (iv) an excluded volume parameter *B*.^{43} These four values parameterize a bottlebrush viewed as a wormlike cylinder with a uniform diameter *d* rather than the individual architectural features, such as the length of the backbone or the side chains.^{43} The parameters are determined from equilibrium averaged conformations obtained from molecular simulations using an ESC model and systematically fitted to analytical expressions for a wormlike cylinder.^{43,46} The diameter *d* is approximately four times the gyration radius of the side chains, and the contour length *L* is found to be shorter than the actual contour length of the backbone due to local undulations within a tube populated by the side chain monomers.^{43} For computational purposes, this model is discretized as a sequence of touching beads with an angle-dependent potential between successive bonded beads and a hard sphere like repulsive interaction between non-bonded beads.^{43} This ISC representation was successful in capturing equilibrium and near-equilibrium properties, such as end-to-end distance, radius of gyration, hydrodynamic radius, and intrinsic viscosity, and accessing much longer bottlebrush lengths than is possible with an ESC model.^{43} Furthermore, this same model is also able to accurately capture the stretching response of a bottlebrush under low force.^{44}

Rooted in these efforts, we applied this ISC model to study the self-assembly behavior of bottlebrush block copolymers in solutions.^{45} We derived an intermolecular interaction potential from scaling arguments, with two parameters to be tuned to match experimental observations.^{45} The intermolecular interaction is a soft potential, in contrast to commonly used hard-sphere potentials such as Weeks–Chandler–Anderson (WCA) potential.^{47} Despite the fact that the simulation results based on the scaling-based potential are consistent with experimental measurements,^{45} we would prefer to use the ESC model itself to determine the interaction between the coarse-grained beads of bottlebrushes using the ESC model. Meanwhile, with the derived potential in mind, we can directly compare the scaling relationships with simulations using the ESC model. It is beneficial for future studies to establish a systematic coarse-graining workflow so that similar ISC models for bottlebrush under different thermodynamic conditions (solvent qualities, monomer chemistries, etc.) can be generated.

In this paper, we investigate the interaction between two bottlebrushes and establish a protocol for incorporating these interactions into a coarse-grained, implicit side-chain representation of bottlebrush polymers. First, we use the ESC representation to calculate the potential of mean force (PMF), osmotic second virial coefficient, and interpenetration function between two bottlebrushes. We use a corresponding ISC representation that incorporates this molecular information (e.g., using the PMF to calculate interactions between coarse-grained beads, bending potential, etc.) and show that the two models give near-quantitatively matched PMF profiles and second virial coefficients. We note that we only consider the athermal condition in this work; however, our methods are not specific to this thermodynamic condition, and thus, the general methodology for understanding PMF between bottlebrush coarse-grained beads that we present here can be extended to different thermodynamic conditions. This can be useful for generating a coarse-grained potential of bottlebrushes, which is applicable in large-scale simulations of bottlebrush systems.

## II. SIMULATION MODELS FOR BOTTLEBRUSH POLYMER MOLECULES

We consider two classes of models: (i) Explicit Side Chain (ESC) model, where each side chain is itself a bead–spring chain, and (ii) Implicit Side Chain (ISC) model, where the side chain degrees of freedom are coarse-grained out. Our goal is to calculate the PMF between bottlebrush polymer molecules using the ESC model and subsequently determine the interaction between the coarse-grained beads of the ISC model which generates the same intermolecular PMF, the latter being calculated based on an ISC representation alone.

### A. Explicit side chain model

Bottlebrush polymers molecules with explicit side chains are modeled by a set of beads connected by finitely extensible nonlinear elastic (FENE) springs. The backbone consists of *N*_{bb} beads where each backbone bead carries a side chain comprising *N*_{sc} beads [see Fig. 1(a)]. The motion of each bead *i* with position **r**_{i} is governed by a potential *U* for the system, which is given by

The potential consists of two contributions: the bonding potential *U*_{s} is between beads that are connected, and the excluded volume potential *U*_{ev} captures the pairwise repulsion between all beads. All beads interact via the purely repulsive WCA potential,^{47} which depends on the distance *r*_{ij} = |**r**_{i} − **r**_{j}| between those beads,

where *σ* = 2*a* is the length parameter, *ɛ* = *k*_{B}*T* is the energy parameter (*k*_{B} is the Boltzmann constant and *T* is the absolute temperature), and *r*_{cut} = 2^{1/6}*σ*. In the case where bead *i* and bead *j* are connected, the form of the spring potential representing bonded interaction is

where *k*_{s} = 30*ɛ*/*σ*^{2} and *r*_{max} = 1.5*σ*. This model has been used in our previous works^{46,48} to study bottlebrush polymers with poly(norbornene) backbone and poly(lactic acid) side chains, where we demonstrated excellent agreement with experimental observations upon suitable parameterization for the chemistry-specific details.

### B. Implicit side chain model

The ISC model is a coarse-grained description of the ESC model that we discussed in Sec. II A. A detailed development of this model was presented in our earlier work,^{43} and we adopt that overall formalism with several modifications to account for PMF data from the ESC model in the pairwise potential. Briefly, the ISC model treats a bottlebrush polymer molecule as a wormlike cylinder, characterized by a contour length *L*, Kuhn length *λ*^{−1}, diameter *d*, and an excluded volume parameter *B*.^{43} These parameters are obtained by a systematic fitting procedure using conformational data determined from molecular simulations using the ESC model.^{43} The wormlike cylinder is discretized into a sequence of beads of diameter *d*, connected by rigid bonds and arranged along the centerline of the molecule to generate a representation suitable for molecular simulations [see Fig. 1(b)]. The number of beads is proportional to the contour length of the backbone *L*_{bb} = *l*_{b}(*N*_{bb} − 1), where *l*_{b} = 1.94*a* is the equilibrium bond length in the ESC model.^{43} By placing coarse-grained beads subsequently along a straight backbone, we have *N*_{CGB} = *L*_{bb}/*d*. However, the shrinking of the backbone due to local fluctuations further causes the contour length of the bottlebrush to be shorter than that of the backbone. This gives *L* = *m*_{L}(*N*_{bb} − 1), where *m*_{L} is a constant that accounts for this contour length disparity (also given in previous work^{43}). The bond length is thus *l*_{CGB} = *L*/*N*_{CGB} = *m*_{L}*d*/*l*_{b}, and the fixed length of the bonds ensures that the contour length of an ISC chain does not fluctuate. By doing so, we obtain an integer number of coarse-grained beads to represent the corresponding ESC chain, while the subsequent beads are slightly overlapping. Thus, the total energy of the system using the ISC model is described by

To account for the molecular stiffness embodied in the Kuhn length, we use a bending potential of the form

where *k*_{θ} = (2*λd*)^{−1} is the bending constant and $\theta ij=cos\u22121u\u0302j\u22c5u\u0302i$ is the complementary bond angle between the unit vectors $u\u0302i$ and $u\u0302j$ along consecutive backbone bonds. We still need the pairwise interaction *U*_{CGB} between the coarse-grained beads to fully specify the ISC model. This will be extracted from the PMF data determined from the ESC model, as will be described in Sec. V.

To emphasize the key difference between the above described ISC model and that introduced in previous work,^{43,45} our choice to represent the molecule with overlapping beads enables the incorporation of the PMF into the ISC model interaction potential (see Fig. 2 and discussion in Sec. II C). This is a consequence of how we set up the PMF calculation for the ESC model and later convert it to the interaction potential in the ISC model, noting that the shrinking of the bottlebrush contour is due to backbone fluctuations in the ESC model. We expect that this will help us better model bottlebrush polymers, especially in the case of inter-molecular interactions where the local shape of the molecule under the ISC representation is crucial.

### C. Overview of workflow

We outline the workflow of this investigation in Fig. 2 to describe the interrelation between the various calculations we perform to both parameterize and demonstrate consistency between the ISC from the ESC models. We first use the ESC model to calculate PMF, shown in the left panel of Fig. 2 [including (a) and (b)]. This involves both calculating the molecular PMF in a traditional manner, where the molecule conformation is unconstrained [see Fig. 2(b)], and calculating the segmental PMF where the two bottlebrush molecules are aligned parallel to each other with straight backbones [see Fig. 2(a)]. In the latter case, the length of backbones are chosen carefully such that the ESC bottlebrush can be mapped to an integer number of coarse-grained beads, which is listed for various architectures in Table I. The segmental PMF can be normalized by this integer number to obtain the pairwise interaction between coarse-grained segments of the ISC model. Using this pairwise interaction as input, we then finalize the ISC model where the bond length is properly shortened, as displayed in Fig. 2(c). This enables us to sample the conformations purely based on the ISC model. We perform the molecular PMF calculation for the ISC model [see Fig. 2(d)] and relate it back to the initial point by comparing the results (PMF profiles, second virial coefficients, and interpenetration functions) with those from the ESC model. This comparison eventually serves as a metric to inform us the suitability of the ISC model that we build.

## III. CALCULATION METHODS

The PMF between two polymer molecules *U*(*r*) in dilute solution is defined as^{49}

where $U12r$ is the total *inter*molecular potential between two molecules with their centers of mass held at a distance *r* apart [see Fig. 3(a)] and the angle brackets denote an average over all possible orientations and all possible pairs of molecules. Computational studies using bead–spring models^{50–54} typically evaluate $U12r$ by summing over the intermolecular interactions between the beads belonging to each molecule, i.e., without any intramolecular contributions; the operation indicated by the angle brackets is performed by averaging over multiple independent pairs of Monte Carlo (MC)-generated conformations and relative orientations. This method of direct evaluation via MC (which we call the **single-chain method**) is particularly suitable for athermal chains in dilute solutions, so we will use this technique for the molecular PMF between bottlebrushes as well.^{55} For the segmental PMF calculation [see Fig. 3(b)], we adopt an alternative approach as a supplement—we perform dynamical simulations with two bottlebrushes present in the box and relax the side chain conformation simultaneously (which we call the **two-chain method**). In the long distance regime, the two methods provide identical PMF results (see Fig. S1 in the supplementary material); however, it can be difficult to sufficiently sample at short distances due to the strongly repulsive WCA potential and large overlap of the molecules. The benefit of two-chain method is that since the two molecules are relaxing at the same time, all these “heavily overlapping events” will be naturally avoided. This method also facilitates the extension of this workflow to more complex environments, such as polymer chains in the presence of ions or other components, as discussed by Withers *et al.*^{56}

In the molecular PMF calculation, we generated a set of equilibrium conformations using a combination of Brownian Dynamics^{57} (BD) without hydrodynamic interaction augmented by global MC moves. The BD position update occurs following the discretized stochastic differential equation,

where **R** is the column vector of length 3*N* (*N* being the total number of beads) comprising all bead positions, *U* = *U*_{s} + *U*_{ev} is the total interaction potential, *δt* is the time step, and **x** is a column vector of length 3*N* containing random numbers drawn from a normal distribution that has a mean of 0 and a standard deviation of 1, which represents displacement due to Brownian motion. The right-hand side of the above equation was evaluated using the positions at time *t* and to obtain the bead positions at time *t* + *δt* on the left-hand side. In Eq. (7), the distance is non-dimensionalized with *a*, time with bead diffusion time *τ* = *ζa*^{2}/*k*_{B}*T* (*ζ* is the bead friction coefficient), and energy with *k*_{B}*T*. Two different MC moves were used when necessary (see the supplementary material of Ref. 46): (i) backbone pivot and (ii) side chain double bridging, both subject to acceptance via the Metropolis criterion.^{58} The BD time step was 10^{−4}*τ*, with a pivot move every 50 BD steps and a side chain double bridging cycle every 11 BD steps. For every (*N*_{bb}, *N*_{sc}) pair, six trajectories were generated, each trajectory running for $>108$ BD time steps. This procedure generated $>6\xd7104$ independent conformations on average for each bottlebrush molecule. An example of bottlebrush conformation using this model is shown in Fig. 1(a). The scheme we use to compute molecular PMF in Eq. (6) is as follows: We choose a pair of conformations from the set of equilibrium conformations and label them as molecule 1 and molecule 2. The center of mass of molecule 1 is then held fixed at the origin and that of molecule 2 is placed with its center at a distance *r* along the *x*-axis. Pairwise overlap between the beads of molecules 1 and 2 is checked using a hierarchy of axis-aligned bounding boxes.^{59–61} For overlapping beads, we calculate the contribution of that pair toward $U12r$. For each value of *r*, the relative orientation of molecule 2 is changed four different times by rigid body rotation to average out the angular dependence. However, this rotation is not necessary if the number of available pairs is large. After one sweep over all values of *r* for a given pair, the same procedure is performed over as many pairs as available. Note that in the molecular PMF calculation, no simulation box or periodic boundary conditions are used, as seen in Fig. 3(a). There is no dynamical evolution of the molecules and the calculation is purely based on the equilibrium conformations that are priorly generated. The above described strategy is also applied to both the ESC model and the ISC model for PMF calculation.

In the segmental PMF calculation between two aligned bottlebrushes, no MC moves are used, since the BD update for side chain beads is sufficient for generating independent conformations. The other aspects (e.g., BD time step and length of the simulation) are all kept the same as those used in the molecular PMF calculations. Two bottlebrush molecules were initialized inside a simulation box such that they cross the periodic boundaries along the *y*-direction, as seen in Fig. 3(b). In this case, we are essentially modeling two infinitely long and aligned bottlebrushes, calculating the interaction force between segments with a length dictated by the simulation box size in the *y*-direction. The backbone beads are fixed in a fully extended conformation, with the two backbones separated by a distance *r* in the *x*–*z* plane. Those two molecules are either using (i) two copies of configurations from the same configuration set or (ii) two independent molecules that are placed in the box simultaneously; these correspond to PMF calculation methods 1 and 2, as previously described. We vary *N*_{sc} from 2 to 20 and choose *N*_{bb} to be such that the corresponding number of coarse-grained beads *N*_{CGB} is close to an integer number. The chosen values of *N*_{bb} as well as the corresponding *N*_{CGB} can be found in Table I. The box size along *y* is set to be *L*_{y} = *N*_{bb}*σ*, which ensures that the backbone bond crossing the periodic boundaries along this direction has the same length as all other bonds, satisfying the periodic boundary condition. The box size along *x* and *z* directions is taken to be large enough to prevent interaction between the periodic images in those directions. Typically, the box extends by more than ten times the side chain radius of gyration beyond the two molecules along *x* and similarly along *z* directions. The average force between two bottlebrush molecules at each distance is obtained after a sufficiently long simulation, and the segmental PMF is calculated via numerically integrating the force using trapezoidal rule. The discretization of the distance is small enough to ensure a good numerical approximation from this integration. Beyond the longest distance that is calculated, the force will be linear extrapolated to obtain a residual triangular area, which contributes a small amount of energy to the overall PMF profile.

We focus exclusively on a segmental PMF for aligned backbones to maximize the points of interaction between bottlebrush molecules to benefit the sampling of the PMF. It is also possible to consider unaligned bottlebrush backbones, for example, with the two backbones positioned orthogonally. For computational simplicity, we do not perform this calculation; in part, we justify this by noting that the difference between a parallel vs orthogonal orientation is still represented even in the ISC due to the connectivity between beads. We also argue based on the results in Sec. V that this is ultimately unnecessary, as pair potentials calculated from aligned bottlebrush PMFs give rise to excellent matching between ISC and ESC models.

Once molecular PMF *U*(*r*) is obtained, the osmotic second virial coefficient *A*_{2} was obtained using the formula^{49}

where *N*_{A} is the Avogadro number, ⟨*S*^{2}⟩ is the mean square radius of gyration, M is the molecular weight, and $\rho =r/\u27e8S2\u27e91/2$. For the bottlebrush molecules studied here, *M* = *Nm*_{0}, where *m*_{0} is the molecular weight of a segment and $N=Nsc+1Nbb$ is the total number of beads. Note that ⟨*S*^{2}⟩ can be directly determined from the set of equilibrium conformations. The interpenetration function *ψ* was obtained from *A*_{2} using the formula^{49}

## IV. POTENTIAL OF MEAN FORCE FROM THE ESC MODEL

We first focus on the results for the PMF between unconstrained bottlebrush polymers. The benchmark results of linear polymers are shown in Fig. S2 of the supplementary material. We highlight two example series of data in Figs. 4 and 5. First, we fix the backbone length as *N*_{bb} = 320 and vary the side chain length from *N*_{sc} = 0–32, as shown in Fig. 4. Here, the case with *N*_{sc} = 0 serves as a linear polymer baseline for comparison. As the side chain length increases, the PMF profile exhibits stronger repulsion between bottlebrush molecules. The upward shift in the PMF is mainly induced by the increased steric repulsion between side chains from different bottlebrushes, which is a consequence of increasing side chain length. However, the extent of the PMF is well-described at all side-chain lengths *N*_{sc} by the radius of gyration $\u27e8S2\u27e91/2$, suggested by the converging behavior of the curves at $r/\u27e8S2\u27e91/2\u22482.0$.

Second, we fix the side chain length as *N*_{sc} = 14 and vary the backbone length from *N*_{bb} = 20–640, shown in Fig. 5(a). In this case, the bottlebrush conformation is closer to a spherical object at the low backbone length limit, while at the long backbone length limit, it acts more like a linear polymer. Therefore, when normalized by $\u27e8S2\u27e91/2$, the length scale is quite different among different *N*_{bb}. With *N*_{bb} = 20 and *N*_{bb} = 40, at the distance $r/\u27e8S2\u27e91/2=0.8$, the two bottlebrush molecules overlap significantly because the side-chains must interpenetrate due to their star-like structure. This gives rise to a marked increase in the PMF due to the lack of acceptable configurations. To contrast, for *N*_{bb} > 80, the PMF profile resembles those in Fig. 4 and goes to a finite value of *U*(*r*) because the coiled structure of the bottlebrush at these lengths allows them to overlap without significant side-chain interpenetration. To evaluate this change in molecular shape in a quantitative manner, we use the volume fraction *N*/*V* within pervaded volume $V\u223c\u27e8S2\u27e93/2$ as a metric relevant for bottlebrush–bottlebrush overlap. As can be seen in Fig. 5(b), *N*/*V* generally decreases as the backbone length increases, indicating a transition of the bottlebrush conformation from a spherical object to a linear-like chain. This is also concomitant with the change in the PMF profile, where the shape of the PMF profile resembles that of the linear polymer in the long-chain limit, while having a hard sphere-like shape in the short-chain limit.

Figure 6 shows the second virial coefficient *A*_{2} for linear and bottlebrush polymers as a function of backbone and side chain lengths. For linear chains, scaling theory predicts^{62} *A*_{2} ∼ *N*^{3ν−2}, which for athermal chains (*ν* = 0.588) reduces to *A*_{2} ∼ *N*^{−0.236}. Fitting to our data for *A*_{2} for long chains, we obtain a scaling exponent of −0.29. This value is slightly lower, indicating that our simulations are still not in the long-chain limit, even for *N*_{bb} = 1280. To compare, earlier work by Dautenhahn and Hall^{52} reported an exponent of −0.272 and Yethiraj and Hall^{50} reported −0.25 after extrapolating their data to *N* → *∞*. Our prediction is in good agreement to these earlier results. Compared to linear chains, the second virial coefficient of bottlebrushes are systematically lower, with *A*_{2} decreasing with an increase in *N*_{sc} with *N*_{bb} held fixed. We attribute this to the increased size of the bottlebrush chain as the side-chain interactions induce an increasing molecular rigidity, leading to less compact structures. This is also consistent with the general trend of decreasing *A*_{2} with an increase in *N*_{bb} in all curves.

Figure 7 shows the interpenetration function *ψ* as a function of overall and backbone molecular weight for both linear chains and bottlebrush polymers. This quantity, as suggested by the name, is a good measurement of how much two molecules are interpenetrating each other at a similar normalized distance (with respect to radius of gyration). For sufficiently long molecules, *ψ* of both linear chains and bottlebrush polymers can be seen to approach an asymptote. Renormalization group theory^{63} predicts that for a self-avoiding chain, *ψ* → 0.219 as *N* → *∞*. We show this as a dashed line in Fig. 7, noting that our results are indeed consistent with this asymptoting value. In the case of bottlebrushes, there is a non-monotonic trend in *ψ* with a peak at short backbone lengths—this reflects the change in the overall shape of a bottlebrush polymer molecule from a spherical star to a cylindrical rodlike shape before transitioning to a chain-like conformation as the backbone length is increased at a fixed side chain length. This trend is related to our observations connecting the intrachain monomer volume fraction to the molecular shape, in that the interpenetration is a way to quantify the effect of short-range interactions on bottlebrush–bottlebrush interactions. The results of segmental PMF between aligned bottlebrushes are shown in Fig. 8, with both the original data from BD simulations (two-chain method) and a series of fitting based on scaling arguments. For all the side chain lengths we examine here, two-chain method provides a rather smooth curve throughout a large range of length scale. We show in Fig. S1 of the supplementary material that the single-chain method is quantitatively matching but becomes unreliable at the low distance regime due to the inability to adequately sample highly interpenetrating conformations. To better understand the trend in these segmental PMF profiles, we apply the scaling analysis that we developed in our previous work to fit the data to an analytical form.^{45} Here, we adopt the following expression:

where the exponent $m=13ln\u206128ln(2/d)$. The constant *C* includes all molecular information that is not a function of *N*_{sc}, as well as the proportionality constant inherent to the scaling analysis that was used to develop this expression. We choose *d* to be the point at which the PMF goes to $\u223c0.1kBT$ (0.2*k*_{B}*T* for *N*_{sc} = 20) to reflect that the scaling argument does not account for long-range fluctuations in interactions due to side-chains that can “stretch” beyond *d*. We thus expect long interaction tails, which we observe as the difference between the scaling theory and simulation data at large *N*_{sc}. We choose *C* = 3.2 as an overall fitting parameter of order unity that leads to excellent matching over all side-chain values *N*_{sc} (Fig. 8). For short side chains (*N*_{sc} = 2 and *N*_{sc} = 4), there is a small degree of deviation that we attribute to the inherent limitations of the scaling picture in describing short chains.

The results we obtain here suggest that the scaling-based theory we derived in our previous work^{45} is already a reasonable approximation of inter-bottlebrush interactions, and in particular, the simple form of this potential is enabling for large-scale simulations that mainly focus on materials properties at a large length scale, such as block bottlebrush self-assembly. This provides a connection between this current study of segmental PMFs from the ESC model and experimentally observable properties; the results in this paper are consistent with pairwise interactions that are used to predict scattering and phase behavior in block bottlebrush polymers.^{45} Now, we have gone further and can show how to incorporate ESC-determined PMFs directly into the ISC model without needing to assume a theoretical form to make predictions of bottlebrush solution properties.

## V. POTENTIAL OF MEAN FORCE FROM THE ISC MODEL

We can use the PMF as an effective pairwise potential in a coarse-grained ISC representation of a bottlebrush so that it is informed by the ESC calculations. We fit a power-law curve to the segmental PMF data from Sec. IV and extrapolate this to extend over all values of *r*. This PMF can then be rescaled and directly incorporated as a pair potential into the ISC model, for which standard MC is performed to generate equilibrium bottlebrush conformations.

The PMF from both the ISC and ESC model can be compared by establishing the correspondence between the backbone length and the number of coarse-grained beads. The PMFs for both ISC and ESC models are plotted in Fig. 9 for *N*_{sc} = 14 using the *N*_{CGB} values listed in Table I to find the values of *N*_{bb} associated with the ISC representations. We include analogous plots of *N*_{sc} (Fig. S3) in the supplementary material. In general, for all *N*_{sc} we considered, we observe near-quantitative matching between the ISC and ESC models. Differences are most noticeable when the backbone is relatively short (i.e., small *N*_{CGB}) in the ISC model. This is attributed to two reasons: (i) end effects, which were not part of the ESC segmental PMF parameterization, become significant in this limit, and (ii) the value of $\u27e8S2\u27e91/2$ that we use as a normalization constant for the *x*-axis is least accurate in this same limit. This latter point is rationalized by the use of only the center of the coarse-grained beads as a proxy for the mass distribution in these chains, despite the non-negligible distribution of mass throughout the coarse-grained bead. While these deviations are noticeable, they mostly become negligible at longer bottlebrush lengths where there is a clear consistency between these two approaches. Indeed, the main motivation to develop the ISC model is to simplify the calculation of high molecular weight bottlebrushes, which is where the agreement between ISC and ESC PMFs is the most quantitative.

We can further consider a comparison of the second virial coefficients and interpenetration functions between ISC and ESC models. As shown in Fig. 10, *A*_{2} from the ISC model is systematically smaller than that of the ESC model. This is mainly attributed to the error in the calculation of the molecular size, which will be enhanced when $\u27e8S2\u27e93/2$ enters the calculation. This disparity is largely normalized out in the PMF profiles such that the integral in Eq. (8) is essentially identical between the two models, and the disparities primarily arise in the prefactor. Thus, the systematic offset in *A*_{2} indicates that one of the biggest disadvantages of the ISC model is that the size of the molecule is underestimated, but this would be a general consequence for all other similar coarse-grained models. In addition, the ISC calculation of *A*_{2} is not monotonically decreasing at the short-chain limit, which we attribute to the aforementioned inaccuracies at the short-chain limit. This is even more notable when calculating the interpenetration functions, plotted in Fig. 11, where *ψ* shows a significantly larger peak at small values of *N*_{bb}. Notably, however, these significant differences become quite small at the long-chain limit where this approach works well.

## VI. CONCLUSION

In this paper, we studied the PMF between bottlebrushes using both an ISC and the ESC model. The molecular PMF profile between bottlebrushes shows a similar trend as the profile of a linear polymer when the bottlebrush is long but shows a qualitative difference when the bottlebrush is short. This trend in the PMF profile is, in part, related to the shape change of a bottlebrush molecule, where it behaves like a spherical object when the backbone is short but resembles a linear chain when the backbone is long. We calculated the segmental PMF between aligned bottlebrushes and found that this PMF can be described by the previously derived scaling form in a precise manner. Finally, by mapping the segmental PMF between aligned bottlebrushes to a potential between coarse-grained segments, we incorporated ESC information directly into an ISC bottlebrush model. We used this ISC model to perform a molecular PMF calculation and found near-quantitative agreement between the ISC and the ESC models. We note that the ISC model is most inaccurate when bottlebrushes have short backbones; however, this disadvantage occurs in the portion of the parameter space where the advantage of the ISC model—the reduction of computational cost vs the ESC model—is not as significant. There are a number of ways in which this model could be improved to address this limitation for short bottlebrushes, for example, explicitly parameterizing the chain ends interactions in the ESC model, although this comes at additional computational cost. For longer bottlebrushes, where the ESC model is particularly costly, the incorporation of these PMF calculations is a useful way to inform much more efficient ISC models that are useful for simulating multi-chain bottlebrush materials. In the future, this methodology holds promise as a way to include other useful interactions (charges, crosslinks, *χ*-driven interactions, etc.) informed by chemistry considerations.

## SUPPLEMENTARY MATERIAL

See the supplementary material for (i) a comparison of the single-chain and two-chain methods for calculating the PMF, (ii) PMF calculations for polymers with different lengths, and (iii) comparisons of the ISC and ESC PMFs for a wide variety of side chain and backbone lengths.

## ACKNOWLEDGMENTS

This work was supported by the National Science Foundation under DMREF Award No. DMR-1727605 and the ACS Petroleum Research Fund under Award No. 61500-ND7.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

T.P. and S.D. contributed equally to this work.

## DATA AVAILABILITY

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