Shear Banding in Simulations of Polymer Melts

Introduction

Polymeric liquids are strongly entangled when the polymer molecular weight $M$ greatly exceeds a critical value $M_c$ [1, ?]. The dynamics of entangled polymers are well-described by tube models such as the Doi-Edwards (DE) model [1, 3, 4], in which polymers diffuse within a ‘tube’ of $Z$ distinct entanglements, and reset their conformation on the scale of the terminal ‘reptation’ time $\tau_d \sim3Z^2$.

The DE model predicts a non-monotonic steady state relation between shear stress $\sigma_{xy}$ and shear rate $\dot{\gamma}$ with a stress maximum when $\dot{\gamma} \simeq1/\tau_d$. This can lead to instabilities such as spurt from a pipe [5, 6] and shear banding [7, 8, 9, 10, 11, 12], whose simplest signature is a steady state constant ‘stress plateau’ in the shear stress as a function of shear rate beginning around $\dot{\gamma} \simeq1/\tau_d$.

Most experiments indicate a slight increase in stress rather than a plateau. [13] To eliminate this inconsistency in DE theory, Marrucci and Ianniruberto [14, 15, 16, 17] incorporated polymer stretch and its relaxation due to nonequilibrium convected constraint release (CCR), which speeds up the terminal time. In CCR, polymers stretched by flow retract more quickly, which accelerates equilibrium (thermal) constraint release, leads to a more disordered tube contour, and increases the shear stress. CCR is controlled by a parameter $\beta$: for large enough $\beta$ the increased stress eliminates non-monotonicity and thus eliminates shear banding [3, 18].

Some highly entangled polymer solutions demonstrate phenomena consistent with non-monotonic constitutive curves, including discontinuous shear rate jumps in stress-driven flow [19] and banded velocity profiles in startup [20], large-amplitude oscillatory [21], and steady-state [22, 23] shear flows. These results are debated due the confounding factors of wall slip [24, 25] and edge fracture [26, 27], with some studies with similar materials and conditions not observing banding [28, 29, 30]. This may be because the formation of shear bands is history-dependent [12, 31, 32, 33, 34, 35]. Wang [36] attributed flow inhomogeneities to disentanglement in the higher-shear-rate band, while other studies [13, 37, 38, 39, 40] have coupled shear banding to concentration fluctuations.

We are not aware of incontrovertible experimental evidence for shear banding in polymer melts. Observations such as the spurt effect [5] and long stress plateaus [41] suggest banding-like behavior, and phenomena such as transient banding [42] and other flow inhomogeneities [43, 44] are again complicated by wall slip [24, 25] and edge fracture [26, 27], which makes steady-state observations difficult at high $Z$ [45]. By contrast, simulations show clear evidence for steady-state banding in polymer melts. Cao and Lichtkman [46], Mohagheghi and Khomami [47, 48], and Ruan et al. [49] observed shear banding in bead-spring polymer simulations of steady-state shear. Most simulations of polymers under shear use the SLLOD algorithm [50], which imposes a uniform shear flow and does not allow shear banding [46].

To understand these phenomena, we consider flow-induced disentanglement, which is encapsulated in CCR. The number of entanglement “kinks” $Z_k$ measured by primitive path analysis [51, 52] has been shown in simulations [49, 53, 54, 55, 56] to decrease in strong flows with $\dot{\gamma} \gtrsim1/\tau_R$, where $\tau_R \sim Z^2$ is the Rouse time, agreeing with CCR predictions. To model this behavior, Ianniruberto and Marrucci [15] applied the CCR mechanism to the Doi-Edwards model and derived equations of motion coupling the reduced number of entanglements $\nu(t) \equiv Z_k(t)/Z_{k,\mathrm{eq}}$ to the polymer orientation tensor, reproducing the steady state reduction of entanglements $\nu(\dot{\gamma})$ found in simulations by Baig, Mavrantzas, and Kroger [53]. More recently, D’Olata and Olmsted [57] (DO) derived a thermodynamically-consistent version of this model coupling CCR to the Rolie-Poly model, and showed, by comparing with simulations, that the CCR parameter $\beta$ is larger for stiffer polymers [56].

In this paper we compare predictions of the DO model for shear banding with observations of shear banding in Kremer-Grest bead-spring simulations [58]. The DO model predicts banding for sufficiently-entangled polymers $Z > Z_c(\beta,\alpha)$, where the critical entanglement number $Z_c$ is larger for larger CCR parameter $\beta$ and smaller relaxation anisotropy $\alpha$ (defined below). The region of banding found in the simulations agrees surprisingly well with the model, given that the DO (and Rolie-Poly) model is a simplified single-mode version of the more precise microscopic GLaMM model[3]. We speculate that physical polymers whose dynamics are controlled by dihedral rotations have larger $\beta$, which may explain why shear banding is not seen in polymer melts.

Theory

DO Model

The DO model couples the evolution of the local tube conformation tensor $\mathbf{A}(t)$, roughly following the Rolie-Poly (RP) model[59], to the reduced number of primitive-path kinks $\nu(t) = Z_k(t)/Z_{k,\mathrm{eq}}$. These quantities, and the local stress tensor $\boldsymbol{\sigma}$, evolve according to the following equations:

Here, $\overset{\triangledown}{\mathbf{A}}$ is the upper convected derivative

$\widetilde{\tau}_d(\lambda,\beta)$ is a non-equilibrium reptation time given by

$\zeta_Z\equiv Z_{k,\mathrm{eq}}/Z$ is the ratio of the number of kinks to the standard entanglement number $Z$; $\lambda\equiv\sqrt{\operatorname{tr}\mathbf{A}/3}$ is the relative stretch of the polymer; and $f(\lambda)$ is the Cohen approximation[60] to the Finitely-Extensible Nonlinear Elastic (FENE) spring force

The maximum relative stretch $\lambda_{\max}\equiv\sqrt{N_{\mathrm{eK}}}$, where $N_{\mathrm{eK}}$ is the number of Kuhn steps per entanglement. This term diverges as $\lambda\to\lambda_{\max}$, but in the limit $\lambda_{\max}\to\infty$, $f(\lambda)=1$, and the behavior of the chain becomes Gaussian.

Previous studies indicate that $Z_{k,\mathrm{eq}}\approx2Z$[53], where $Z$ is derived from the linear modulus $G_0$[1]; the DO model assumes an exact relation $Z_{k,\mathrm{eq}}=2Z-1$[57]. The relation between the Rouse time $\tau_R$ and the (equilibrium) reptation time $\tau_d$ is given by

where the Likhtman relation

incorporates equilibrium constraint release and contour length fluctuations[61].

The CCR parameter $\beta$ is a positive parameter of order 1, which determines the strength of convective constraint release. The relaxation anisotropy parameter $\alpha$ controls a Giesekus-like term in the model[62], which allows for nonzero second normal stress differences not permitted by the RP model. Its value should lie between 0 (for isotropic relaxation) and 1 (for maximum anisotropy), and it determines the ratio of the normal stress differences in the low-shear-rate limit: $\alpha=-2\lim_{\dot{\gamma}\to0}(N_2/N_1)$[57], where $N_1$ and $N_2$ are the first and second normal stress differences. Existing data on polymer solutions indicate $\alpha\approx0.6$[63, 65], while data on melts suggest $\alpha>0.5$[66, 71] but is not conclusive. We discuss this issue in more detail in Section IV.

To account for fast fluctuations at scales smaller than an entanglement, we add a Newtonian stress $2\eta_p\mathbf{D}$ to the stress term from the original DO model, where $\eta_p$ is a background viscosity and $\mathbf{D}=\frac{1}{2}[\nabla\mathbf{v}+(\nabla\mathbf{v})^{\mathsf{T}}]$ is the fluid’s symmetric strain rate tensor. Following Adams, Fielding, and Olmsted[18], we parameterize this background viscosity by $\varepsilon=\eta/\eta_p$, where $\eta_p=G_0\tau_d$ is the low-shear-rate viscosity of the entangled liquid.

For melts, the background viscosity is $\eta_e\equiv c\zeta N_{\mathrm{eK}}b_K^2/36$, which is the (Rouse) viscosity for shorter chains of $N_{\mathrm{eK}}$ Kuhn steps of length $b_K$ with Kuhn segment number density $c$ and monomer friction constant $\zeta$[1]. The plateau modulus of the full chains is $G_0=ck_BT/N_{\mathrm{eK}}$, so the viscosity ratio becomes

For dilute solutions, $\eta$ is the solvent viscosity, so $\varepsilon$ no longer has this exact form, but still scales as $\eta\sim\eta_sZ^{-3}$.

The DO model has five dimensionless parameters: $Z$, $\beta$, $\alpha$, $\lambda_{\max}$, and $\varepsilon$. For melts, $\varepsilon$ is determined by $Z$, while $Z$, $\lambda_{\max}$, and $\alpha$ can be determined from linear or weakly-nonlinear rheometry along with the modulus and time scales. This leaves $\beta$ as, strictly, the only free parameter.

Predictions of Banding

To determine when the model predicts banding, we numerically solve for $\nu$ and $\mathbf{A}$ in steady state homogeneous simple shear flow to calculate the shear stress $\sigma_{xy}$ as a function of the shear rate $\dot{\gamma}$ (Figure 1). We express this shear rate in terms of two different Weissenberg numbers in this paper: the reptation Weissenberg number $\mathrm{Wi}_d = \dot{\gamma}\tau_d$ based on the terminal relaxation time, and the Rouse Weissenberg number $\mathrm{Wi}_R = \dot{\gamma}\tau_R$ which specifically governs disentanglement. These constitutive curves are monotonic for some parameter combinations, but become non-monotonic for $1/\tau_{d,\mathrm{eq}} \lesssim\dot{\gamma} \lesssim1/\tau_R$ for others, which coincides with a reduction in entanglements $\nu(\mathrm{Wi}_d) < 1$. When $Z_k = \nu Z_{k,\mathrm{eq}}$ drops below a value of about 3 (horizontal dashed lines in Figure 1), the concept of an entanglement tube becomes poorly-defined and we expect the model to fail.

We use these constitutive curves to construct a “phase diagram” for monotonicity, and thus for banding (Figure 2). The monotonicity boundary depends weakly on $\lambda_{\max}$ and $\epsilon$, even for extreme variations $\lambda_{\max} \in[2,\infty)$ (solid vs. dashed lines) and $\epsilon\in(0,\epsilon_{\mathrm{melt}}]$ (solid vs. dotted) (see Appendix C for a discussion).

Increasing $Z$ promotes banding[1, 18, 20] by increasing the separation between the Rouse and reptation times, broadening the range of shear rates in which the stress can decrease (compare the green curves in Figure 1ab, for which all other parameters are identical). Strengthening CCR (increasing $\beta$) inhibits banding, as proposed by Ianniruberto and Marrucci [72]. Figure 1a demonstrates the mechanism: stronger CCR increases disentanglement (smaller $\nu$) when $\dot{\gamma} \sim1/\tau_R$, which decreases polymer alignment and increases $\sigma_{xy}$. This value of $\dot{\gamma}$ is typically near the local stress minimum, so a sufficient increase of $\sigma_{xy}$ there can eliminate the nonmonotonicity.

Significantly, a larger relaxation anisotropy $\alpha$ promotes banding. As evident in Figure 1b, increasing $\alpha$ decreases the shear stress for a range of shear rates above $1/\tau_d$, which promotes nonmonotonicity. Mechanically, increasing $\alpha$ accelerates the relaxation of $\mathbf{A}$ along directions in which it is elongated, so the polymer component of $\sigma_{xy}$ decays more quickly as $\dot{\gamma}$ increases.

Model and Methodology

To test these predictions, we performed non-equilibrium molecular dynamics simulations of semiflexible Kremer-Grest (KG) bead-spring polymers[58] using the large scale atomic molecular massively parallel simulator (LAMMPS) software[73]. The model consists of beads of mass $m$ and diameter $\sigma$ connected with FENE bonds, which interact with a Lennard-Jones potential with interaction strength $\epsilon_{\mathrm{LJ}}$ cut off at $r_c = 2^{1/6}\sigma$. To vary the entanglement length, a 3-body cosine angle interaction $k_\theta\epsilon_{\mathrm{LJ}}[1+\cos(\theta)]$, with $k_\theta= 0$, $1.5$, and $2.0$ was included. Chains lengths $N_m = 200$ were studied, depending on the bending stiffness. The number density is $0.85\sigma^{-3}$ for all systems.

Most of the systems are in a cubic box with periodic boundary conditions, with equilibrium configurations prepared following the procedure described Auhl et al. [74]. To determine $\beta$, some systems were then replicated in the $x$ direction four times and equilibrated using the double bridging algorithm[74] before applying shear. See Appendix A for a list of all systems studied.

The nonequilibrium MD simulations were performed using two different procedures. The first uses the SLLOD equations of motion[50, 75, 76, 77] with a Nose-Hoover thermostat[78, 79] and a temperature damping constant of $\tau\equiv\sqrt{m\sigma^2/\epsilon_{\mathrm{LJ}}}$. This algorithm imposes a linear velocity profile in the shear plane and thus prevents any shear banding. In the second procedure, the simulation cell is deformed at a constant rate in the shear plane using the fix deform command in LAMMPS, and the beads are coupled to a pairwise DPD thermostat[80, 81, 82] to control the temperature. This thermostat acts only between neighboring beads with a distance $r_c$ and damping coefficient $\gamma= m/\tau$. As it does not impose a uniform velocity profile, this procedure allows shear banding. Both procedures give the same shear stress under all conditions, even when the system forms shear bands using the second method while the first method suppresses the shear bands.

Due to computational limitations arising from the need to measure time scales ranging from bond relaxation times to that of the entire polymer chain, the strain rates $\dot{\gamma}$ currently accessible in this study are $10^{-6}$ to $10^{-2}/\tau$. The time step for all shear simulations was $0.008$–$0.01\tau$, and all simulations were performed at constant volume. We compute the entanglement Rouse time $\tau_e$, number of Kuhn steps per entanglement $N_e^K$, and ratio of Kuhn step to packing length $b_K/p$ for these systems using formulas from Everaers et al. [83], then compute $\tau_d$ and $\tau_R$ using (7) and the relation $\tau_R = Z^2\tau_e$; these quantities are given in Table 1 and Table 2. We determine $Z$ from equilibrium simulations by $Z = \langle L_{\mathrm{pp}}\rangle^2/\langle R^2\rangle$, where $\langle R^2\rangle$ is the mean squared end-to-end length and $\langle L_{\mathrm{pp}}\rangle$ is the mean primitive path length measured by the Z1+ code[52].

Model	Ref.	$k_\theta$	$N_{\mathrm{mon}}$	$Z$	$\beta$	Bands	$b_{\mathrm{K}}/p$	$N_{\mathrm{eK}}$
F-KG	[56], †	0	$500^{*}$	6	$0.03 \pm0.01$	N	2.8	39
F-KG	†	0	1000	13	$0.03 \pm0.01$	N	2.8	39
F-KG	†	0	1500	19	$0.03 \pm0.01$	Y	2.8	39
F-KG	†	0	2000	26	$0.03 \pm0.01$	Y	2.8	39
SF-KG	†	1.5	200	6	$0.17 \pm0.08$	N	6.8	8.8
SF-KG	†	1.5	300	9	$0.17 \pm0.08$	N	6.8	8.8
SF-KG	†	1.5	$400^{*}$	12	$0.17 \pm0.08$	Y	6.8	8.8
SF-KG	†	2	200	8	$0.19 \pm0.06$	N	9.9	5.5
SF-KG	[49], †	2	$300^{*}$	12	$0.19 \pm0.06$	Y	9.9	5.5
SF-KG	†	2	400	16	$0.19 \pm0.06$	Y	9.9	5.5
MK	[48]	2	$200^{*}$	7	$0.28 \pm0.02$	N	11.2	7.9
MK	[48]	2	250	9	$0.28 \pm0.02$	Y	11.2	7.7
MK	[48]	2	400	14	$0.28 \pm0.02$	Y	11	8
UA-PE	[85]	–	400	5	$1.25 \pm0.11$	–	12	6
UA-PE	[86]	–	700	8.7	$1.11 \pm0.08$	–	12	6.2
UA-PE	[87]	–	1000	12.9	$0.98 \pm0.16$	–	10	6.4

$k_\theta$	$N_{\mathrm{mon}}$	$N_{\mathrm{chains}}$	$\frac{L_{y,z}}{\sigma}\times10^{-1}$	$\frac{L_x}{L_{y,z}}$	$\frac{\tau_e}{\tau}\times10^{-3}$	$\frac{\tau_R}{\tau}\times10^{-5}$	$\frac{\tau_d}{\tau}\times10^{-6}$	Type
0	500	500	6.7	1	10.9	4.5	1.9	Neither
0	500	16000	13.3	4	10.9	4.5	1.9	SLLOD
0	1000	2000	13.3	1	10.9	18.1	24.5	Both
0	1500	1000	12.1	1	10.9	40.8	101.7	DPD
0	2000	1000	13.3	1	10.9	72.5	272.3	DPD
1.5	200	250	3.9	1	1.7	0.6	0.2	SLLOD
1.5	200	2000	7.8	1	1.7	0.6	0.2	DPD
1.5	300	3200	10.4	1	1.7	1.4	1.1	DPD
1.5	400	400	5.7	1	1.7	2.6	3.2	SLLOD
1.5	400	1600	9.1	1	1.7	2.6	3.2	DPD
1.5	400	3200	11.5	1	1.7	2.6	3.2	DPD
1.5	400	12800	11.5	4	1.7	2.6	3.2	SLLOD
1.5	800	400	7.2	1	1.7	10.3	35.6	SLLOD
2	200	2000	7.8	1	1.2	0.8	0.5	DPD
2	300	400	5.2	1	1.2	1.7	2.2	Both
2	300	3200	10.4	1	1.2	1.7	2.2	Both
2	300	10800	15.6	1	1.2	1.7	2.2	DPD
2	300	12800	10.8	4	1.2	1.7	2.2	SLLOD
2	400	3200	11.5	1	1.2	3.1	6.0	DPD
2	400	12800	11.5	4	1.2	3.1	6.0	SLLOD

To measure disentanglement $\nu= Z_k/Z_{k,\mathrm{eq}}$ as a function of shear rate, we analyzed the simulations using the Z1+ code[52] to compute $Z_k$. The Rouse Weissenberg number $\mathrm{Wi}_R = \dot{\gamma}\tau_R$ at which $Z_k$ begins to drop below its equilibrium value in the DO model depends strongly on $\beta$, so by fitting its prediction of $\nu$ vs $\mathrm{Wi}_R$ to our simulations, as shown in Figure 3, we are able to determine appropriate values of $\beta$. The $\nu$ vs $\mathrm{Wi}_R$ curves we calculate assuming non-shear banded flow, so we perform this fitting on systems using the SLLOD shearing procedure, however measuring disentanglement on DPD simulations which do not shear band gives very similar results. In this fitting process, we set $\lambda_{\max} = \sqrt{N_e^K}$ in the DO model and arbitrarily choose $\alpha= 0.5$, as the latter has little effect on $\nu$.

This technique of determining $\beta$ was introduced by Dolata et al. [56], using the earlier Z1 code[84]. However, the simulation sizes used by Dolata et al. were small enough that at high shear rates, the flow-aligned polymer chains would wrap around the periodic boundary conditions and interact with periodic copies of themselves. This artificially decreased the value of $Z_k$ measured by Z1 and Z1+, which by default do not count kinks which result from a chain interacting with itself; the hollow symbols in Figure 3 are data from Dolata et al. and show the magnitude of this increased disentanglement when compared with the new, larger systems. Z1+ is less susceptible to this issue than Z1, and replicating the system in the shear direction before running Z1+ mostly mitigates the issue, but to be certain we only perform our $\beta$ fitting on simulations which are longer in the shear direction than the maximally-extended chain length.

As a result of this discrepancy, the values of $\beta$ we find for our simulations are somewhat lower than reported in Dolata et al. [56] for analogous systems, but the main trends reported still hold: $\beta$ is greater for chains with greater stiffness $k_\theta$, but does not substantially depend on chain length for systems of the same stiffness (compare $k_\theta= 2.0$ with chains of length 300 and 400 beads in Figure 3, which agree within error bars). Consequently, we assume $\beta$ to be constant for any given $k_\theta$, including for chain lengths whose disentanglement we have not measured.

To supplement our own simulations, we use the same method to fit $\beta$ in simulations by Mohagheghi and Khomami[47, 48] (MK), who investigated banding in a slightly different bead-spring model with three different lengths of chain, and also measured disentanglement using the Z1 code. Their simulations were all performed

using a DPD procedure, and two of the three lengths exhibited banding, so we used their reported $Z_k$ vs $\mathrm{Wi}_{\mathrm{R}}$ data for the non-banding system to fit $\beta$.

Table 1 collects the fit values of $\beta$ for the KG systems studied and the three systems studied by Mohagheghi and Khomami [47], along with the estimated values of $Z$ and whether they exhibit banding. The normal stress ratio $-\alpha/2$ cannot be determined with confidence, but is independent of $Z$ and $k_\theta$, as discussed in Section IV.

Figure 4 shows how the inferred value of $\beta$ depends on chain stiffness, quantified by the ratio of the Kuhn length $b_{\mathrm{K}}$ to the packing length $p$. Dolata et al. [56] suggested that $\beta$ should be a simple function of $b_{\mathrm{K}}/p$, which they approximated by a power law $\beta\sim(b_{\mathrm{K}}/p)^{1.9}$. We find that a slightly less steep power law $\beta\propto(b_{\mathrm{K}}/p)^{1.6}$ agrees with our KG simulations and the results from Mohagheghi and Khomami. However, Dolata et al. also computed $\beta$ and $b_{\mathrm{K}}/p$ for a series of United-Atom Polyethylene (UA-PE) simulations performed by Nafar Sefiddashti, Edwards, and Khomami [85, 86, 87, 88], which do not agree with the computed power law. These simulations are small enough that they likely suffer from the same self-entanglement issue as Dolata et al., inflating the measured value of $\beta$, but this effect may not be enough to account for the discrepancy.

Observing Banding

We observe steady-state shear banding in many of the shear simulations performed using the second shear procedure, depending on the length and stiffness of the chains and the rate at which they are sheared. Figure 5 shows velocity profiles for a system with 3200 chains of 400 beads with $k_\theta= 2.0$, corresponding to $Z = 16$, at six different shear rates. All were started from an equilibrium configuration and the shear velocities are averaged over 50 evenly-spaced bins in the gradient direction after the system reached steady-state. At $\dot{\gamma}\tau= 10^{-3}$, the velocity profile is completely uniform, while at $10^{-4}$, it develops several small bands, and at $10^{-4.5}$ and $10^{-5}$ it clearly separates into one fast and one slow band. At $10^{-5.5}$, the velocity gradient is small enough that local fluctuations become prominent, but there are still distinct bands at higher and lower shear rates, while at $10^{-6}$ only uncorrelated fluctuations remain, marking the end of shear banding. These simulations lack the overall stress gradients present in experimental Taylor-Couette and cone-plate geometries, so configurations with multiple bands are more stable than in experiment.

Other lengths and stiffnesses of chain follow the same pattern: flow is uniform at high shear rates, separates into bands at lower shear rates, and becomes uniform again at lower shear rates, for large $Z$. This is in line with the behavior predicted by a non-monotonic constitutive curve, in which banding occurs at shear rates near the non-monotonic region of the constitutive curve. Yet others do not band at any shear rate, as would be predicted by a monotonic constitutive curve. Figure 6 shows which simulations exhibit banding, overlaid on the phase-space boundaries from Figure 2. For each stiffness $k_\theta$ (and corresponding value of $\beta$), chains with $Z$ below some particular value $Z_c(k_\theta)$ never exhibit banding (hollow symbols), whereas those with $Z > Z_c(k_\theta)$ band at some shear rates (solid symbols). In all cases, $Z_c(k_\theta)$ is close to that predicted by the DO model for corresponding values of $\beta$ and plausible values of $\alpha$, despite the DO model being only a single-mode version of the more accurate microscopic GLaMM model[3]. However, our uncertainty in $\alpha$ limits the specificity of the DO model’s predictions, and the most likely range of $\alpha$ indicated by literature data is higher than would best agree with $Z_c(k_\theta)$. Furthermore, the simulations show a slight trend for $Z_c(k_\theta)$ to decrease for larger $\beta$, opposite the trend predicted by the model. We can think of no mechanism which would produce such a trend, given the role of CCR in preventing banding.

fig:7 shows the predicted reptation Weissenberg number $\mathrm{Wi}_d$ for banding, as a function of $Z$, compared to the simulation data for $k_\theta= 2.0$ with the fit value of $\beta= 0.19$; equivalent plots for other values of $k_\theta$ are in Appendix B. The top row of points in this figure represent the system whose velocity profiles are shown in Figure 5, with uniform flow at the highest and lowest shear rates and banding in between. Banding should occur at shear rates roughly in the shaded region of the plot, between the dash-dotted lines whose position depends on $\alpha$. For any non-monotonic constitutive curve, such as that shown in fig:7, uniform flow is strictly unstable in the range where $\sigma_{xy}$ decreases as a function of $\dot{\gamma}$, and should separate into bands in a process analogous to spinodal decomposition in fluid-fluid phase separation. These bands must exhibit the same selected shear stress $\sigma_{xy}$ to be stable, just as phase separated fluids must have the same pressure. This selected stress could be determined by augmenting the DO model with gradient terms[11, 89, 90] and solving for the stress at which the interface between shear bands is stationary, though we have not done so.

The selected stress plateau must be between the local maximum (“top jumping”) and local minimum (“bottom jumping”) on the constitutive curve, indicated by horizontal dotted and dashed lines in fig:7. We approximate the true selected stress as the average of these two values (the dash-dotted line). Each of these potential plateau stresses intersects the stable parts of the constitutive curve at two points, indicated by vertical lines in fig:7 and plotted in fig:7 as a function of $Z$ for two different values of $\alpha$. These two points are the predicted shear rates in the two bands, and in an infinitely-large system banding should occur exclusively at shear rates between them.

The predicted reptation Weissenberg number $\mathrm{Wi}_d$ that demarcate the banding region are consistently lower than the banding shear rates observed in the simulations, regardless of $\alpha$. We have identified several factors that could contribute to this discrepancy: (1) high-shear-rate breakdown of tube models from instability and disentanglement, (2) history dependence of banding behavior, and (3) finite-size effects in our simulations.

The DO model and other tube models are expected to break down in the high-shear-rate limit, as the loss of the majority of the entanglements shown in Figure 1 eliminates the constraints that define the tube, and as faster dynamical modes are activated that are not captured in these models. Furthermore, Lewy et al.[91] have shown that the RP model can develop 2D instabilities at these shear rates, which if present in the DO model would in-

caption: Range of reptation Weissenberg numbers for which the DO model predicts shear banding for a melt of chain lengths $Z$, with $\alpha= 0.5$ or $0.25$, $\lambda_{\max} \to\infty$, and $\beta(k_\theta=2) = 0.21$ (fit in Table I). Symbols represent SF-KG simulations with $k_\theta= 2.0$ and whether they band in all, some, or no runs. The blue line is the range of $\mathrm{Wi}_{\mathrm{d}}$ where banding was measured by Ruan et al. [49], for the same KG parameters. Line styles indicate method of calculating banding range, as illustrated in top panel and explained in text. Crosshatched region indicates where banding is expected in an infinite system but may be unobservable or absent due to finite-size effects.

:::

validate our assumption of locally-uniform shear.

Unbanded homogeneous flow has also been shown to be metastable in some conditions where banding is the most stable state[11, 33], so it could require much longer simulations or different initial conditions to induce banding. The presence or magnitude of banding can depend on the shear history of the system. For instance, the inset of Figure 5 shows two systems which have reached steady state at the same final shear rate, but have different shear histories: the solid curve started directly from equilibrium, but the dashed curve was sheared at a rate with more pronounced bands first, resulting in more pronounced bands. We also see the reverse effect, where initially shearing at a much higher (non-banding) shear rate can prevent banding at a rate where it otherwise occurs, consistent with metastability. For example, as shown in fig:13a, starting from a non-banded high shear rate of $\dot{\gamma}\tau= 10^{-3}$ and reducing to $\dot{\gamma}\tau= 10^{-5}$ (which forms bands when starting from rest), no shear band is observed, at least over the time scales accessible. The reverse process, in which we start at a lower shear rate and increase to a higher rate, generally gives the same result as starting from rest in terms of presence/absence of banding (as shown in fig:13b of the appendix). Continuously ramping the shear rate from rest to $\dot{\gamma}\tau= 10^{-5}$ gives a similar long-time velocity profile to the abrupt startup. Similarly, the shear band disappears when increasing from $\dot{\gamma}\tau= 10^{-5}$, (which has a shear band) to $\dot{\gamma}\tau= 10^{-3}$ (which does not).

Finite size effects can also prevent or obscure the formation of bands, particularly at low shear rates. When computing the banding range predicted by the model from its constitutive curve, we assume that the fast and slow bands will have fixed shear rates $\dot{\gamma}_{\mathrm{fast}}$, $\dot{\gamma}_{\mathrm{slow}}$ at the selected stress plateau (on the constitutive curve), and that the bands will obey the lever rule[12]: they will resize themselves so that the overall shear rate $\dot{\gamma}_{\mathrm{net}} = \phi_{\mathrm{fast}}\dot{\gamma}_{\mathrm{fast}} + \phi_{\mathrm{slow}}\dot{\gamma}_{\mathrm{slow}}$ is equal to the applied shear rate, where $\phi_{\mathrm{fast}} = \ell_{\mathrm{fast}}/L_y$ and $\phi_{\mathrm{slow}} = \ell_{\mathrm{slow}}/L_y = 1 - \phi_{\mathrm{fast}}$ are the thicknesses of the two bands as a fraction of the overall sample thickness $L_y$. Expressed in terms of Weissenberg numbers, we thus have

This assumption breaks down in systems with a finite thickness $L_y$ when the width of one of the bands becomes smaller than molecular sizes. To estimate when this effect becomes relevant for simulations of various sizes, we assume that the high shear rate branch cannot be observed in practice when $\ell_{\mathrm{fast}}$ drops below $\sim3$ monomers in size or $\phi_{\mathrm{fast}} < 0.03$ (in our simulations, the box sizes (Table II) are all of order $100\sigma$). The crosshatched regions in fig:7 (and Figure 9 in Appendix B) correspond to shear rates at which $\phi_{\mathrm{fast}}$ is between 0 and 0.03. In principle, the region where $\phi_{\mathrm{slow}} \equiv1 - \phi_{\mathrm{fast}}$ is less than 0.03 should also be shaded, but due to the logarithmic scale of the $\mathrm{Wi}_{\mathrm{d}}$ axis, this region is negligibly small.

Normal Stress Ratio

In non-Newtonian fluids, the deviation of the stress from the isotropic pressure is frequently quantified in terms of normal stress differences

for stresses $\sigma_{xx}, \sigma_{yy}, \sigma_{zz}$ in the shear, gradient, and vorticity directions respectively. Most complex fluids reduce to Newtonian fluids in the linear regime at very low shear rates ($\mathrm{Wi}_{\mathrm{d}} \ll1$), so both normal stress differences scale quadratically with the shear rate $\dot{\gamma}$ as $\dot{\gamma} \to0$[92]. Consequently, the ratio $-N_2/N_1$ (which is positive for polymer fluids because $N_2$ is typically negative) is predicted to approach a constant as $\dot{\gamma} \to0$. The value of this constant is not clear; in the DE model, the limiting value is $1/7$ or $2/7$ depending on whether or not the Independent Alignment Approximation holds[1]; in the IM model, the limiting value is $1/4$[93]; and in the DO model, it is $\alpha/2$[57]. Identifying this limiting value from experiments or simulations would fix the value for $\alpha$ in our model, improving the precision of our predictions.

Outside the linear regime most models predict that $N_2/N_1$ decreases with increasing shear rate $\mathrm{Wi}_d > 1$, in the shear thinning regime. Due to the quadratic dependence of both normal stress differences on the shear rate, the forces which must be measured to determine the normal stresses are orders of magnitude lower in the linear regime compared to the thinning regime, which makes it difficult to measure both regimes with the same apparatus. Various groups [63, 64, 65] have measured the normal stress ratio in well-entangled polymer solutions in the linear regime using flush-mounted pressure sensors and cone-partition-plate rheometers, and have found a nearly constant ratio of $-N_2/N_1 \approx0.3$ (corresponding to $\alpha\approx0.6$). Magda and Baek [65] also measured the normal stresses at shear rates above $1/\tau_d$ and found that $-N_2/N_1$ decreased as predicted, although instabilities such as edge fracture prevented their highest-concentration samples from reaching steady state in this regime.

Conversely, groups measuring the normal stress ratio in melts have mostly focused on the shear-thinning regime, despite the fact that polymer melts are even more susceptible to edge fracture and similar instabilities than solutions. Schweizer and co-workers [66, 67, 68, 69] and Vlassopoulos and co-workers [70, 71] measured the normal stress difference in several polystyrene melts, using cone-partitioned-plate rheometers to delay the onset of edge fracture. Data from these groups is shown in Figure 8; due to the large

uncertainty in the values of $N_2$ at low shear rates, the low-shear-rate limit cannot be determined with confidence, but the data suggest that $\alpha_{\mathrm{melt}} \gtrsim\alpha_{\mathrm{solution}}$. To better compare data between different sources, we recalculated the reptation time for each of the melts in question using the entanglement Rouse time $\tau_e$ and the entanglement molecular weight $M_e$ derived by Li, Mavromanolakis, and Vlassopoulos [71], as well as using their time-temperature superposition curve to estimate the temperature dependence of $\tau_e$, but applying the Likhtman relation ((8)) in the derivation rather than the simpler expression $\tau_d = 3Z\tau_R$ used in their calculations.

Figure 8 shows all of the normal stress data from polymer melts and simulations of polymer melts that we could find in the literature. More information on these simulations is Table III. Most of the data, both experimental and simulation, is for large Weissenberg numbers, and none of the data show the plateau expected at low shear rates where both $N_1, N_2 \sim\dot{\gamma}^2$.

Source	Polymer	Length	$M_w/M_n$	$T$	$Z$
[71]	PS	283 kDa	1.01	160°C, 170°C	16.6
[71]	PS	143 kDa	1.11	130°C, 160°C	8.4
[66]	PS	158 kDa	2.85	190°C	9.3
[67]	PS	200 kDa	1.06	175°C	11.8
[68]	PS	206 kDa	1.06	160°C, 170°C, 180°C	12.1
[53]	UA-PE	400 CH$_2$	1	450 K	6
[86]	UA-PE	700 CH$_2$	1	450 K	8.6

Measurements of the normal stress ratio in simulations are similarly focused on the shear-thinning regime. Low-shear-rate simulations require much more simulation time to reach steady state than high-shear-rate ones, and measurements of $N_1$ and $N_2$ are much more susceptible to fluctuation noise due to their smaller magnitude, so the linear regime is impractical to reach with current methods and hardware. Outside of these restrictions the normal stress differences are much more straightforward to measure in simulation than in experiment, and edge fracture is eliminated due to periodic boundary conditions. Baig, Mavrantzas, and Kroger [53] and Nafar Sefiddashti, Edwards, and Khomami [86] measured the normal stress ratio in their simulations of UA-PE (which has dihedral angles and rigid bond lengths and angles), which we reproduce in Figure 8, along with measurements from our own KG simulations (which has only springs and no prescribed bond or dihedral angles).

The resulting $-N_2/N_1$ vs $\mathrm{Wi}_d$ curves collapse reasonably well between experiment and the two simulation types, in line with the DO model’s prediction that the normal stress ratio primarily only depends on $\mathrm{Wi}_d$ and $\alpha$. This suggests that the value of $\alpha$ may be universal between different polymer melts, but the data at low $\mathrm{Wi}_d$ is insufficient to determine its value. The decay in $-N_2/N_1$ occurs over several orders of magnitude of $\mathrm{Wi}_d$, whereas the DO model predicts a more rapid drop off (solid lines); more detailed models may match this decay better.

The UA-PE simulations diverge somewhat from the KG data at high shear rates, which we speculate is due to sub-entanglement-scale differences between the two polymer representations (UA-PE vs KG); the divergence begins when $\mathrm{Wi}_d \sim10^3$, at which point $\dot{\gamma} \sim\tau_e^{-1}$.

Conclusions

The DO model makes specific predictions for the conditions under which entangled polymer melts and solutions should exhibit shear banding. These predictions agree qualitatively and semi-quantitatively with molecular dynamics simulations of bead-spring polymers. Shear banding is controlled by $Z$, the normal stress ratio $-\alpha/2$, and the CCR parameter $\beta$. The normal stress ratio is notoriously difficult to determine for polymer melts, although its value appears to be universal. Regardless of the value of $\alpha$, for the small values of $\beta\simeq0.03$–$0.3$ that best fit the simulations the DO model predicts that shear banding should occur for chains with $Z \gtrsim10$, which agrees with simulations. The overestimate by the model of the range of shear rates that support banding suggests that the model is inadequate at high shear rates for significant degree of disentanglement, as noted by Wang [94]. The prediction that the entanglement cutoff for banding $Z_c$ should increase with increasing $\beta$ and thus increasing $k_\theta$ does not appear to hold for the range of stiffnesses in our simulations. Given the role that convective constraint release plays in avoiding nonmonotonic constitutive curves, it is unexpected that more rapid disentanglement should be correlated with a higher tendency to shear band rather than a lower one. There may be some confounding effect that contributes to this discrepancy, such as the qualitative difference in the bending potentials of flexible and semiflexible KG chains. However, there may also be a flaw in the methodology of measuring $\beta$: the model’s description of high-shear-rate entanglement dynamics may be inadequate, or the identification of kinks with entanglements may break down under flow or alignment. A similar polymer-specific variation of constraint release was found by Matsumiya and Watanabe [95], who measured equilibrium constraint release experimentally. In contrast, they found that polymers with stronger equilibrium constraint release tended to be less stiff, which if applicable to CCR would be consistent with the trend in Figure 6. However, it is also possible that the detailed microscopic mechanisms of equilibrium constraint release and convected constraint release (termed ‘entanglement stripping’ by Hawke et al. [96]) differ significantly. Analysis[56] of data from UA-PE models[85, 86, 87, 88] for polyethylene shows a substantially larger $\beta\simeq1$. This estimate of $\beta$ is likely inflated due to finite system size leading to self-entanglement, but these effects may not be sufficient to account for the discrepancy from bead-spring simulations. United-atom polyethylene has a single dihedral angle for every interior atomistic monomer on the chain, compared to the three degrees of freedom for each (coarse-grained) ‘monomer’ in the KG bead-spring model. This structural difference may account for the increased strength of convective constraint release above that predicted by the chain stiffness, and the richer microscopic detail of fully-atomistic polymers with bending and stretching energies could support similar large values for $\beta$. For these values, the predicted crossover value for banding $Z_c(\beta,\alpha)$ is much larger and more sensitive to $\alpha$, perhaps as high as $Z_c = 20$–$300$ depending on $\alpha$, which may contribute to the lack of banding observed in experimental melts of comparable stiffness. If the trend in $\beta$ vs. $b_k/p$ shown in Figure 4 holds for experiments, and if the predicted trend that higher $\beta$ leads to higher $Z_c$ is accurate for larger differences in $\beta$, we expect that flexible polymers such as PDMs (for which $b_k/p = 2.8$[83]) would be more likely to exhibit banding.

Acknowledgments

We thank Joe Peterson, Ben Dolata, Thomas O’Connor, Marco Galvani Cunha, and Shi-Qing Wang for discussions; and Martin Kröger for generous help with and advice about the Z1+ code. PDO and LLN thank the National Science Foundation (DMREF-2118769), and PDO thanks Georgetown University and the Ives Foundation for financial support. This work was performed, in part, at the Center for Integrated Nanotechnologies, an Office of Science user facility operated by the U.S. Department of Energy (DOE) Office of Science. Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. DOE National Nuclear Security Administration under contract no. DENA-0003525. The views expressed in this article do not necessarily represent the views of the U.S. DOE or the United States government.

Author Declarations

The authors have no conflicts to disclose.

Data Availability

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

Simulation Details

Table 2 lists the specific systems sizes of all the simulations we have performed, along with the three relevant timescales $\tau_e$, $\tau_R$, $\tau_d$ of each system and which of the shear procedures has been used for that system. For systems with non-cubic simulation boxes, $L_{y,z}$ is the box length in the gradient and vorticity directions and $L_x$ is the length in the shear direction, so $L_x/L_{y,z}$ is the aspect ratio of the box.

Additional Banding Range Plots

Fig. fig:7 of the main text compares the range of banding shear rates between theory and simulation for chains with bending stiffness $k_\theta= 2.0$; Figure 9 shows the equivalent data for $k_\theta= 0$ and $k_\theta= 1.5$. The simulations with $k_\theta= 1.5$ agree with the theory better than $k_\theta= 2.0$, but those with $k_\theta= 0$ agree considerably worse, in line with their greater discrepancy in $Z_c$.

Effects of $\varepsilon$ and $\lambda_{\max}$

Figure 10 illustrates the effects of varying $\lambda_{\max}$ and $\varepsilon$ on the constitutive curves of the DO model. Neither parameter has much effect on monotonicity, but they do affect the constitutive curve at relatively high shear rates, which are typically outside the non-monotonic regime.

The maximum stretching ratio $\lambda_{\max}$ has a weak effect because the departure from linear elastic behavior only becomes relevant to steady-state banding once the tube has stretched considerably, which only happens for shear rates exceeding $1/\tau_R$. For smaller values of $Z$, $1/\tau_R$ is closer to the start of the plateau regime, so the slight increase in stress due to smaller $\lambda_{\max}$ is more likely to eliminate a decreasing regime. The value of $\lambda_{\max}$ does have a substantial impact on tube disentanglement at high shear rates, so it is always relevant for fitting $\beta$ from disentanglement data in simulations.

The viscosity ratio $\varepsilon$ becomes relevant only when the

associated stress (black dotted line in Figure 10) becomes significant relative to the tube stress. While this can occur in the plateau regime and affect monotonicity, it only does so for relatively large $\varepsilon$ and $\varepsilon$ as $\varepsilon$ scales as $Z^{-3}$, is only relevant for small $Z$. Figure 11 illustrates this: the boundaries between monotonic and non-monotonic constitutive curves in $\varepsilon$-$\beta$ space become independent of $\varepsilon$ for $\varepsilon\lesssim10^{-4}$. The values of $\varepsilon_{\mathrm{melt}}$, illustrated by the symbols along the boundaries, are well below that point for all $Z \gtrsim10$.

Figure 11 shows several alternate slices through the monotonicity phase space of the model, illustrating the effect of $\varepsilon$ more clearly. fig:11a is analogous to Figure 3 of Adams, Fielding, and Olmsted [18], with a logarithmic scale for $\varepsilon$. fig:11b illustrates the effects of varying $\alpha$ in this perspective, but also allows for direct comparison to several versions of the Rolie-Poly model, discussed below.

Rolie-Poly Comparison

The DO model is directly based on the Rolie-Poly model, with $Z$ and $\beta$ corresponding directly to equivalently-named parameters in that model. It incorporates several additional features, such as a dynamical equation of motion for disentanglement, finite extensibility, anisotropic relaxation, and conformation-dependent relaxation; if we eliminate some of these features by setting $\alpha= 0$ and $\lambda_{\max} \to\infty$, it should be quite similar to the RP model. Figure 12 shows how phase boundaries in Figure 2 compare to those predicted by the Rolie-Poly model (dotted), along with two variations of the RP model.

The first variation (dashed) is the ‘Likhtman’ Rolie-Poly model, in which the relationship $\tau_d = 3Z\tau_R$ in the original model has been replaced with the more detailed Likhtman relation $\tau_d = 3Zf_L(Z)\tau_R$ derived by Likhtman and McLeish [61], which is used in the DO model. The difference between the two becomes negligible in the high-$Z$ limit, but applying the correction makes the two models identical when $\beta= 0$.

The second variation is due to Chen et al. [97], who investigated both models in the high-$Z$, nonstretching limit.

By expanding both models in terms of $1/Z$, they found that the two models become identical in this limit when $\alpha= 0$, except for a rescaling in the value of $\beta$:

(This is a corrected version of eq:50 of [97]; per private communication, the term in parentheses was incorrectly given as $1-\frac{1}{3}\beta_{\mathrm{DO}}^{2}$). By replacing $\beta$ in the RP model with this rescaled value, in addition to applying the Likhtman correction factor, we get a Rescaled Likhtman Rolie-Poly model, the dash-dotted line in Figure 12, whose monotonicity behavior is close to the DO model with $\alpha= 0$.

Banding profiles

Figure 13 presents steady state profiles for a system of 3200 chains of 400 beads with $k_{\theta} = 2.0$ for different shear histories

Normal Stress Ratio Data

Table III lists some properties of the polymers in the previous experiments and simulations measuring the normal stress ratio which are plotted in Figure 8.