Shear Thinning from Bond Orientation in Model Unentangled Bottlebrush Polymer Melts

The rheology of molecular brushes remains challenging to control due to the multiple length scales and relaxation processes involved and the lack of direct observation of molecular conformation during flow. We use molecular dynamics simulations to determine the shear thinning of unentangled bottlebrush polymers with varying architecture, from linear chains to combs, to densely grafted bottlebrushes, to star-like and star polymers. We find shear thinning exponents in line with theoretical and experimental results and characterize the shape and orientation of bottlebrushes in steady-state flow. Many shape parameters derived from the gyration tensor show molecular alignment with the flow for all systems. Yet, the orientation of individual bonds is what most strongly correlates with the architecture-dependent shear-thinning exponents. In densely grafted bottlebrushes, the packing of side chains prevents alignment with the flow, causing a reduction in shear thinning. The molecular insight from our simulations is useful to tune the architecture of bottlebrushes to control their rheology.


Introduction
Bottlebrushes are a type of branched polymer composed of a linear backbone and multiple side chains grafted to it. Their characteristic architecture and the dense packing of side chains lead to unique properties that have been studied over the last few decades [1][2][3][4][5][6][7][8][9] . The shape and the size of bottlebrushes are determined by three molecular parameters: the length of the backbone, the length of the side chains, and the spacing between these side chains (grafting density). At equilibrium, a bottlebrush is more rigid and extended than a linear chain with the same backbone length due to the side chains' steric repulsion, leading to higher persistence length, stiffness, and increased packing efficiency 10 . The overall shape of a bottlebrush polymer is highly dependent on the interplay of architectural parameters.
Increasing grafting density results in more compact and cylindrical molecular shapes. For a fixed grafting density, increasing side-chain length stretches the backbone chain and the bottlebrush becomes more rigid 5 . This rich behaviour leads to highly tunable mechanical and thermal properties of bottlebrushes compared to linear polymer chains 11 and to softer and more flexible networks 12 . The control over the architectural composition makes bottlebrushes suitable candidates for multiple applications in biomedical sciences 13 , soft robotics 14 , pressure-sensitive adhesives 15 , capacitive pressure sensors 16 , and ultrasoft electronics 17 .
While the configuration of bottlebrushes at equilibrium is somewhat understood 11,18 , it remains challenging to predict their shape rearrangement under mechanical deformations and flow. This prediction is necessary to determine their rheological behaviour and to use them in technological applications such as 3D printing 19 . It is well-known, for example, that polymer melts exhibit shear thinning 20 due to the alignment of the polymer chains with the flow direction 21 . In branched polymers, the shear thinning is governed by multiple length scales and relaxation times 22,23 , which gives rise to unique rheology compared to linear chains 24-28 that is not trivial to predict from the molecular architecture. This is complicated by the challenges in the controlled synthesis and characterisation of branched polymers 29,30 . Several studies have been conducted to investigate the rheological properties of melts or solutions of bottlebrushes (synthetic or biological) such as their viscoelasticity 31,32 , shearinduced crystallinity 33,34 , and shear-thinning behaviour in aggrecan solutions 35  The primary goal of this work is to study the effect of bottlebrush architecture on the viscosity of the melt for unentangled chains under shear flow. We use non-equilibrium molecular dynamics (NEMD) simulations to calculate the flow curve in steady shear for different bottlebrush architectures at a fixed molecular weight and compare them with the limit cases of linear chains, comb polymers, and stars. We demonstrate that all polymer melts exhibit shear thinning, and we also observe a transition from star-like to bottlebrushlike structures. However, in the investigated frequency range, we find that the thinning exponent is determined by the orientation with the flow of individual bonds rather than chain orientation. Densely grafted bottlebrushes shear-thin the least because the packing of the side chains prevents their strong alignment with the flow direction, and this effect decreases with decreasing grafting density from the bottlebrush to the comb regime. Many other shape parameters that are derived from the principal components of the gyration tensor also show molecular alignment with the flow but do not directly correlate with the shear thinning behaviour as a function of molecular architecture. Our simulations provide important insight into the bond and chain configuration changes under shear flow that are crucial to predict and control the rheology of bottlebrush polymer melts in different architectural regimes.

Methods
In this work, we use coarse-grained molecular dynamics (CGMD) simulations under nonequilibrium conditions.

Model
We employ the bead-spring model with reduced Lennard-Jones (LJ) units where the energy (ε), mass (m), and diameter (σ) of a single bead are set to 1 leading to a unity time scale as τ = mσ 2 /ε 44 . The bonded interactions in the bead-spring model are given by the finite extensible nonlinear elastic (FENE) potential to represent the bonds between connected monomers, where k = 30ε/σ 2 is the spring constant, R 0 = 1.5σ is the maximum allowed bond length, and r is the distance between two bonded beads. Non-bonded interactions between the beads are given by the truncated LJ potential to represent the Van der Walls interactions up to a cutoff distance r c = 2.5σ, where σ * = 2 1/6 σ is the length scale at which the LJ potential attains its minimum value with depth ε. Additionally, we compare bottlebrushes to a system of linear chains with length L = 100 and a system of star polymers with f = 9 arms and M = 11 beads per arm. Note that these two additional architectures also satisfy the constraint we set on the molecular weight M w = 100 with 2, 000 chains. The chain parameters are summarised in Table 1.

Equilibration protocol
To prepare the polymer melt, we first randomly pack the desired number of polymer chains in the simulation box followed by an energy minimisation step. Then we equilibrate the system in the NPT ensemble (constant number of particles, pressure, and temperature) using a Nose-Hoover thermostat at temperature < T >= 1.0 and pressure < P >= 0.1.
We choose this pressure value since it is close to atmospheric pressure 47 . We run an initial equilibration for 5 × 10 7 simulation steps and additional 5 × 10 6 production steps with an integration step of 0.005τ ensuring that the initial end-to-end vector correlations of polymer chains are lost (see Fig. S1). We consider the longer subchains of bottlebrushes (either the backbone or the side chain) in the correlation calculations. For star polymers, we consider the chain segment from the core to a terminal arm bead. The correlations are averaged over all subchains and for all chains in the melt. We also have a final average over 3 independent simulation runs each of which starts with a different initial configuration of polymer chains at different random initial velocities.

Deformation protocol
After the creation of the polymer melt, we apply a simple shear with a fixed shear rate on the xy−plane where the directions of shear, gradient, and vorticity arex,ŷ, andẑ respectively.
The SLLOD equations of motion are integrated with Lees-Edwards boundary conditions at constant volume. We measure the average stress values σ xy of all beads in the system which is a function of the total strain γ =γt. We then calculate the viscosity η = σ xy /γ of each system after deforming them up to a total strain of γ = 100. The overall shearing procedure is the same as in Ref. 48 . We explore a range of shear rates fromγ = 10 −3 toγ = 0.2 similar to Ref. 49 observing unstable viscosity data for shear rates lower than 10 −3 as they also report in their work. We observe an unexpected drop in the pressure of the system when we try to simulate at lower shear rates. Also, at shear rates larger than 0.2, the bonds in the system become unrealistically large and cause the simulations to fail. To ensure the stability of our simulations at high shear rates, we reduce the integration time step to 0.001τ and check that the temperature fluctuates stably around the initially set value.

Definition of architectural parameters
We computed the following architectural parameters, most of which are a function of gyration tensor eigenvalues. The correlation between these parameters is shown in Fig. 5. We plot them as a function of the shear rate in Fig. S3 along with the corresponding eigenvalues in • Orientation resistance tan 2ϕ = <Gxy> <Gxx>−<Gyy>

Shear thinning of polymer melts
Here we show the shear thinning behaviour of polymer melts with different molecular architectures. We calculate the viscosity as where σ xy is the shear stress at steady flow given by the average per-monomer stress, andγ is the shear rate. The viscosity flow curves of polymer melts with different architectures are reported in Fig. 2 and show a clear shear thinning trend. We fit the data with a power-law model η = Kγ n−1 where n is the shear thinning exponent and K is the flow consistency index. Shear thinning exponents for all systems are reported in Table 2. The value of n ranges between 0 and 1 and smaller values of n indicate stronger shear thinning. The prediction of this exponent for unentangled linear chains is n = 0.5 which is known from the theory of polymer physics and confirmed by experiments 50,51 . We can confirm that the linear chains in our system are not entangled by the measurement of the scaling exponent n. Since the bottlebrushes in our systems are unentangled, we also choose unentangled linear chains to make a better comparison between different architectures. Note that for entangled systems, the initial drop in the viscosity during the start-up shear is attributed to the chain disentanglement 52 . Hence, we exclude this mechanism by choosing unentangled linear chains, as later proven by their shear-thinning behaviour (see Table 2 ). We would expect a kink in the viscosity curve aroundγ ∼ 10 −4 based on the relaxation times τ R of the longest subchain of our polymers (see Fig. S1). For shear rates lower than τ −1 R , we would reach the Newtonian plateau. Although we try to reach this plateau, the viscosity data fluctuate largely forγ < 10 −3 , and the simulations become unstable. Nevertheless, based on the relaxation times τ R in Fig. S1 and Rouse model arguments for unentangled polymers, we might expect a higher viscosity at lower shear rates for linear chains and other elongated architectures, given that the total molecular weight is fixed 50 . A crossover point is then expected at intermediate shear rates aroundγ ∼ 1/τ R , which is roughly what we observe in  respectively. Therefore, the increasing spacing between grafted side chains enhances the shear thinning. We discuss the molecular origin of these trends in the following sections.

Bond orientation parameter
Orientational order is the driving mechanism in the shear thinning of complex fluids 21 , but its definition for branched polymers is not trivial. In fact, one would expect the shear thinning to be correlated to some metric derived from the gyration tensor, but we report later in detail that this is not the case for our system. Instead, we find a strong correlation between shear thinning and the orientation of the bonds. To quantify it, we measure the bond orientation parameter < P 2 > which is defined as where P 2 is the second order Legendre polynomial, θ is the angle between the bond vector and the unit vectorx. The average < . > is taken over all bonds in a single chain and over all chains in the system. We also calculate P 2 for backbones and side chains separately (see Fig. S2). < P 2 > takes values between −0.5 and 1 where < P 2 >= −0.5 indicates a perpendicular alignment to the flow direction, and < P 2 >= 1 a perfect parallel alignment.
Under equilibrium conditions, < P 2 > becomes 0 indicating a random orientation without any directional preference. We also measure and report < P 2 > at equilibrium in Table 3 confirming that the orientation under equilibrium conditions is 0 as expected.
We show the values of < P 2 > at all shear rates in Fig. 3 with solid markers. < P 2 > converges to 0 at low shear rates, recovering the equilibrium value, whereas it increases with increasing shear rate as a result of an increased bond orientation. We expect P 2 to reach a plateau at the infinite shear rate limit imposed by both topological constraints and the chain stiffness 53 . Although the perfect alignment with flow indicates a value of 1 for P 2 , this limit is never reached because of the entropic contributions on the order of ∼ k B T . The interactions between different chains prevent them from moving on a straight line on the x-axis and hence aligning perfectly. Plateau values of P 2 have been previously reported up to 0.5 for squalane, a comb-like unentangled polymer, at different thermodynamic conditions 54 and up to 0.8 for rigid colloidal rods 53,55 . Note that reaching this plateau experimentally is not possible due to the high shear rates required, and even our simulations become numerically unstable beyond the data points presented here. Nevertheless, we can fit our data with the generalised logistic function (solid lines) that is given by the expression P 2 (γ) = a+ −a 1+(γ c ) 0.82 . The fit parameter a describes the behaviour at the infinite shear rate regime when the bond orientation reaches a plateau, and c is a characteristic shear rate at which bonds transition from their equilibrium configuration to an alignment with the flow. 1/c is roughly related to the inverse of the bond relaxation time, on the order of 10 2 τ , as we discuss in Fig. S1(c).
The exponent ofγ c determines how sharp the transition is from 0 to a. The value 0.82 is the best fit to data and is fixed for all systems. The fact that a single exponent is sufficient for all systems indicates that it does not strongly depend on the architecture, though we don't have a particular insight into this specific parameter.
In the range of shear rates studied, the progressive bond orientation with the flow governs the shear-thinning behaviour of all architectures: the higher the bond orientation, the stronger the shear thinning, see the inset of Fig. 3. We observe that bottlebrush orientation increases with increasing m for the same side-chain length (N SC = 4). The loosely-grafted bottlebrushes with m = 4 orient better than more densely-grafted bottlebrushes with m = 2 and m = 1. This trend is the same as their high shear rate viscosity values. The increase in

Shape Parameters
We expect that the orientation of the molecules with the flow plays an important role in determining the system rheology, and here we calculate different shape parameters by computing the radius of gyration tensor that is commonly used to determine the shape, size, and structure of polymer chains 56 . The components of the gyration tensor are given by where N is the total number of beads in a single chain, r ik is the Cartesian components (x, y, z) of the k th monomer and r (CM ) i is the i th component of the centre of mass (CM) of the chain. We diagonalise the gyration tensor to obtain its eigenvalues λ 2 max > λ 2 mid > λ 2 min in descending order. The trace of this tensor gives the radius of gyration R g for a polymer chain as R 2 g = λ 2 max + λ 2 mid + λ 2 min . We show the radius of gyration R 2 g and the diagonal components G ii in Fig. 4 which are normalised by their respective R 2 g0 values at equilibrium to identify the changes under shear with respect to equilibrium conditions. The equilibrium values of R 2 g0 are given in Table 3 along with the average bond orientations of different architectures. Note that the curves in Fig. 4 all start from 1 at equilibrium.
The increasing G xx and decreasing G yy (Fig. 4) indicate that all architectures extend in the flow direction and compress in the gradient direction leading to the orientation of the whole molecule in thex-direction. Note that linear chains have the highest increase in R 2 g and G xx , while the lowest increase is observed for star polymers. Moreover, we infer a compression along the vorticity direction from the decrease in G zz as expected from shear-thinning polymer melts 57 . We agree with Ref. 57 that increasing branching makes G zz insensitive to shear. This compression is more prominent for linear chains as they have more space to compress in the z-direction than other architectures. The star-like (5,19,1) and the sphere-like (10,9,1) bottlebrushes also show a decrease in G zz at low shear rates. For the remaining systems, G zz either stays constant or slightly increases with the increasing shear rate. This indicates that the chains resist the orientation along the flow direction and, as a result, expand in the vorticity direction. Note that G xx does not correlate with the viscosity directly. Some architectures orient better in the flow direction, yet they do not show a stronger shear thinning. What we infer from the analysis of the conformational properties is that the molecules can orient with the flow as a whole; however, shape parameters based on the gyration tensor are not the key predictors of the architecture dependence on shear thinning.

Correlations of architectural parameters with viscosity
One could think that it is not G ii , but some geometrical shape combination of them, e.g. asphericity, which correlates with the shear thinning. We test various architectural parameters which are a function of gyration tensor components with different weights and combinations.
Indeed, all shape parameters investigated behave as expected under flow. For instance, the asphericity increases and this increase is more prominent for stars and star-like systems than bottlebrush-like systems. However, no other shape parameter correlates as well with shear thinning as P 2 , across different architectures. We quantify this by the Spearman coefficient given in Fig. 5, and the full data sets of these parameters as a function of shear rate for all architectures are reported in Fig. S3-S4. We consider these parameters as a function of shear rate, which is an independent variable, and compute the correlation between two data sets by sorting them from low to high shear rate. The correlations between these parameters and the viscosity are reported in Table 4 in descending order of the absolute value. Although G yy correlates with viscosity as much as P 2 does, it cannot explain the ranking of viscosities for different architectures at high shear rates. On the other hand, P 2 is able to describe this trend for all architectures. A natural parameter for this study would be G xx as we apply the shear in the x-direction; however, it does not correlate with viscosity. Hence, the reason for high correlations in G yy is a secondary effect as the chains have to compress in the gradient direction when the shear rate is increased. There is no strong correlation between the viscosity and all the other quantities reported in Table 4. Figure 5: Spearman correlation coefficient between different parameters. We report the correlations between shape parameters and the viscosity that are given in Table 4 sorted by decreasing absolute value. We also see the correlations between different shape parameters. The red block in the middle of the map corresponds to high correlations of shape parameters expressed in terms of gyration tensor eigenvalues. Hence, we expect them to be somewhat correlated. Smaller values of the relative shape anisotropy (κ) determine how isotropic a configuration is. Similarly for asphericity (b) and acylindricity (c). Detailed definitions of these parameters are given in the Methods.
Our findings align with recent experimental results and provide a precise picture of the molecular rearrangements giving rise to the observed rheological behaviour of our melts. Like Zografos et al. 38 and Lopez et al. 27

Conclusions
We have studied the shear thinning of bottlebrush polymer melts under shear flow using non-equilibrium molecular dynamics simulations. We explore a wide range of architectures at fixed molecular weight, transitioning from linear chains to combs, to densely grafted bottlebrushes, to star-like and star polymers. We report a range of shear thinning exponents from a Rouse regime for linear and star-like systems to a Zimm regime for cylindrical bottlebrushes, observing a star-like to bottlebrush-like transition in line with recent experiments 38 . We observe directly the alignment of all molecules with the flow direction at an increasing shear rate, which is more pronounced for linear chains and for bottlebrush-like polymers than for stars due to their anisotropic shape. We quantify molecular alignment in the form of several parameters extracted from the gyration tensor such as the principal components along and perpendicular to the flow, and shape parameters such as asphericity or prolate- ness. Yet, in the range of shear rates considered, the alignment of molecular bonds explains the architecture effect on shear thinning, rather than the molecular alignment. Linear chains show the strongest shear thinning because all bonds can easily align with the flow, while this is hindered for the side chains of cylinder-like bottlebrushes, which then shear thin the least.

Code availability
In an effort to promote Open Science practices, the LAMMPS code used to perform our simulations is available on the GitHub page of our group https://github.com/giuntoli-group/

Characteristic Relaxation Time Scales
We estimate the relaxation times of different architectures by computing the end-to-end auto-correlation function in Fig. S1. For bottlebrushes, estimating this time is complicated since there are different relaxation times related to different parts of the chain. Hence, we consider the longer sub-chain in determining this relaxation time. We compute ⃗ r ee of all sub-chains that satisfy max(N BB , N SC ). After obtaining the auto-correlation function, the relaxation time is estimated as the time at which this function decays to 0.2. Figure S1: (a) Auto-correlation function of the longest chains in a molecule. For long bottlebrushes, ⃗ r ee refers to the end-to-end position vector of the backbone or the side chain depending on the number of beads on each sub-chain. For stars, it is the position vector from the core bead to the terminal bead in one arm. For linear chains, it is the end-toend position vector of the whole chain. The end result is obtained by averaging ⃗ r ee for all molecules in the system. (b) Auto-correlation function of the bond vectors ⃗ r b . We compute the correlations by considering all bonds in a given system. (c) The relation between the bond relaxation time τ b and the fit parameter c. τ b is the time at which the bond autocorrelation function reaches the value 0.2. Although there is no clear correlation between the two, our fit parameter is of the same order of magnitude as the bond relaxation time.

Bond orientation of backbones and side chains
We compute the P 2 parameter separately for bottlebrush backbones and side chains in Fig.   S2, excluding the linear chain and the star polymer that don't have both types of bonds. We obseve that both the backbone bonds and side chains bonds orient more easily as we start to decrease the grafting density of side chains without changing their length N SC . Figure S2: Bond orientation parameter P 2 of backbones (left) and side chains (right) are calculated separately. Backbone and side chain bonds orient more easily as the number of grafted side chains is decreased at fixed number of side chain beads N SC .

Architectural Parameters
We compute several architectural parameters based on the gyration tensor in Fig. S3. These parameters are a function of the gyration tensor eigenvalues. Figure S3: (a) Asphericity increases with increasing shear rate. The chains are more isotropic at low shear rates and they are deformed differently depending on the architecture.(b) Acylindricity. Bottlebushes are more cylindrical under equilibrium conditions and their shape is deformed under shear, whereas linear chains are more globular at equilibrium. They become more cylindrical as the chains are extended under high shear deformation (c) Prolateness. Larger values indicate more prolate, and smaller ones indicate more oblate shapes. (d) Relative shape anisotropy. All architectures are more isotropic at equilibrium which then is reduced with increasing shear rate. (e) Orientation resistance parameter Eigenvalues of the gyration tensor (see Fig. S4) help us identify the molecular size, shape, and symmetry. They also give insights on the molecular conformation. Figure S4: (a) Largest eigenvalue is directly proportional to the size of the polymer. It also defines the principal symmetry axis, i.e. the corresponding eigenvector, of the molecule. λ 2 max increases with shear rate for all architectures indicating an extension along the molecule's symmetry axis. Significantly larger λ 2 max compared to other eigenvalues indicates an elongated shape. (b) Middle eigenvalue (c) Smallest eigenvalue