Understanding changes to microstructural dynamics under nonlinear deformations is critical for designing flow processes of entangled polymeric fluids, motivating the development of experimental methods to probe strain- and rate- dependent modifications to relaxation mechanisms. Although orthogonal superposition rheometry (OSR) holds promise as such a probe, the ability to interpret the superposition moduli accessible by OSR in the context of entangled polymer dynamics remains an open question. To fill this gap, we report model OSR predictions using detailed microstructural models for both monodisperse and polydisperse entangled polymers, i.e., the Rolie-Poly and the Rolie-Double-Poly models, respectively, which account for reptation, chain retraction, and convective constraint release. By combining numerical calculations with a perturbation analysis, we demonstrate that for polymers that can be described by a single-mode model, the OSR superposition moduli at different shear rates and frequencies can generally be collapsed onto a single master curve, with rate-dependent shift factors that depend on the nonlinear rate-dependent modification of polymer conformation and relaxation rates without changing the dominant relaxation mechanisms. We systematically study how the OSR moduli are sensitive to the shape and dispersity of the molecular weight distribution. We discuss the generality of our results for a broad class of constitutive models and suggest an analogy to Laun’s rule to relate OSR moduli to the first normal stress difference. Our results provide a foundation to guide the design and interpretation of future experiments and demonstrate that orthogonal superposition rheometry often probes features in nonlinear dynamics more directly than conventional rheometry techniques.

## I. INTRODUCTION

Flow processing of soft materials and complex fluids usually involves nonlinear deformations, which can nontrivially modify both the associated microstructural configuration and dynamics. This is particularly important for entangled polymeric fluids, where entanglements produce nonlinear relaxation mechanisms that are significantly modified when polymers stretch and orient in flow, leading to various flow instabilities and complications in polymer processing. Various recent flow protocols have been proposed to highlight the nonlinear contributions to the relaxation: one is large amplitude oscillatory shear (LAOS) [1] and another is superposition rheology, which was originally proposed by Tanner and Simmons [2,3]. Although in LAOS the frequency dependence offers the possibility to probe different time scales, the method relies on subjecting the sample to a very complex kinematic history [4]. Alternatively, superposition rheology still enables frequency dependent measurements while retaining simpler kinematics. As a result, the superposition of a small-strain oscillatory motion onto a steady or transient shear flow can provide a clearer insight into the effects of flow on the mechanisms underlying the nonlinear response of rheologically complex fluids [5]. The oscillatory motion can be imposed either parallel ($\gamma \u02d9\u2225$, parallel superposition, PSR) or perpendicular ($\gamma \u02d9\u22a5$, orthogonal superposition, OSR) to the direction of motion for the steady shear flow. Of these, orthogonal superposition offers the advantage that the two flow fields are not coupled and, as we show in this work, the corresponding moduli can be more directly related to the microstructural configuration under flow. OSR provides measurement of the frequency- and shear-rate-dependent nonlinear viscoelastic superposition moduli (storage and loss moduli, $G\u22a5\u2032$ and $G\u22a5\u2033$), which are analogous to their linear viscoelastic (LVE) counterparts that characterize equilibrium dynamics. So, in principle, OSR is a very promising experimental method for studying the coupling of microstructural dynamics and nonlinear rheological response in polymeric liquids.

In practice, however, it has been difficult to achieve OSR measurements experimentally due to instrumental challenges and measurement sensitivity. The earlier devices for OSR by Simmons [3], Mewis and Schoukens [6], and Zeegers *et al.* [7] have a complex and delicate mechanical design and only cover a rather narrow viscosity range. Vermant *et al.* [8] described an orthogonal superposition technique based on a simple modification of the force rebalance transducer in a commercial rheometer and introduced the design of a double wall Couette cell with open bottom for the inner cylinder to avoid pumping flow. These new instrument developments finally enabled the implementation of sensitive orthogonal superposition on commercial rheometers. A recent experimental and simulation study thoroughly investigated the calibration procedures and corrections needed for the nonidealized flow field in the OSR setup on commercial rheometers [9]. Orthogonal superposition rheometry (OSR) has been used to study structural changes during flow in complex fluids, including polymeric fluids [10,11], starlike and wormlike micelles [5,12], and both in Brownian dynamics simulations [13] and experiments for colloidal suspensions [4,14–17].

Because OSR has only recently become experimentally available on commercial rheometers, relatively little is known about how to interpret the nonlinear viscoelastic results in the context of entangled polymer dynamics. Specifically, there are very few theoretical and computational studies that would provide a fundamental basis for interpreting orthogonal superposition measurements. To date, most of the experimental studies have interpreted the data by analogy to the interpretation of linear viscoelastic measurements. However, though proposed analogies between linear viscoelasticity and orthogonal viscoelasticity have been used (e.g., using empirical rules to estimate chain length [18–20]), these analogies have not been theoretically verified. This provides strong motivation for a deeper theoretical and computational study to investigate what is fundamentally being measured in the OSR experiments.

To better understand superposition experiments, the resulting moduli should be described by suitable rheological constitutive models. Lacking a generally valid nonlinear rheological constitutive equation, previous studies have used various phenomenological models. Yamamoto proposed a strain rate-dependent relaxation spectrum, which demonstrates that the parallel moduli reflect the coupling between shear and superimposed flow, while orthogonal superposition moduli do not [21]. Wong and Isayev predicted the orthogonal moduli with the Leonov model [22]; Kwon and Leonov later corrected the predictions [11]. The corrected predictions agree qualitatively well with the early experimental work by Simmons on polymer solutions [10]. Kim *et al.* used the Giesekus model to predict parallel and orthogonal superposition moduli and compared them with experiments on wormlike micelles [5].

For studying entangled polymers in superposition flows, only two molecularly based models have been used. The first involved a prediction of parallel superposition moduli by Unidad and Ianniruberto [23] using the differential constitutive equation accounting for convective constraint release (CCR) proposed by Marrucci and Ianniruberto: the double-convection-reptation model with chain stretch (DCR-CS model) [24]. However, the interpretation of parallel superposition results is complicated by the coupling of the steady shear flow and oscillatory flow. The second study by Mead [25] uses the monodisperse and polydisperse MLD models to predict PSR and OSR moduli [26–28]. A potential weakness of the MLD model is that the stretch and orientation dynamics are treated using separate dynamical equations, which can lead to problems even in monodisperse rheology predictions [29]. As pointed out by Boudara *et al.* [30], the decoupling approximation is even more problematic for polydisperse blends because couplings between constraint release and chain retraction must be readmitted in an *ad hoc* fashion and do not arise naturally from the model. Furthermore, the studies by both Mead as well as Unidad and Ianniruberto neglected chain stretch. So far, no studies have systematically studied the effect of polydispersity on OSR moduli using models that treat stretch and orientation in the same equation and probed a sufficiently wide range of shear rates and frequencies so that significant shifts in the moduli can be observed with increasing shear rate.

The current investigation addresses these limitations of previous studies by using molecularly based models that treat orientation and stretch in the same equation for predicting OSR moduli. Specifically, we use the Rolie-Poly model [31] for monodisperse polymers and the Rolie-Double-Poly model [30] for polydisperse polymers. These models account for the detailed nonlinear relaxation processes of entangled polymers including reptation, chain retraction, and convective constraint release. Furthermore, by combining a perturbation analysis and numerical calculations, we can obtain OSR moduli across a wide range of frequencies and shear rates and systematically vary the degree of chain polydispersity. After the detailed modeling study, we will discuss our results in a broader context to help design and interpret future OSR experiments. We find that the relationships governing OSR moduli for the Rolie-Poly model also work for a broader class of constitutive models. We suggest an analogy to Laun’s rule [32] to relate OSR moduli to the first normal stress difference. Additionally, we derive expressions to extract model parameters from experiments without the need for detailed model fitting. By making an analogy to the linear relaxation spectrum, we show the utility of OSR experiments to obtain material structural information that would be difficult to obtain from conventional rheological measurements. In this way, studying orthogonal superposition computationally not only provides an important comparison to future experiments but also provides better sensitivity for testing models for nonlinear polymer processing.

## II. THEORY

### A. Storage and loss moduli in linear viscoelastic measurement

Storage and loss moduli have been used ubiquitously to describe viscoelastic materials, such as entangled polymers. The storage modulus $G\u2032$ is proportional to the energy storage and, therefore, reflects the elasticity in the system. The loss modulus $G\u2033$ is proportional to the dissipation or loss of energy, reflecting the viscous nature of the system [33]. The most common approach to obtain information about the deformation-induced structural rearrangement within a complex fluid is to measure the two moduli using small amplitude oscillatory shear (SAOS). In this regime, the material response is independent of the strain amplitude, and this is known as the linear viscoelastic behavior [34].

The storage and loss moduli at equilibrium can be used to estimate the average contour length of entangled polymers and wormlike micelles, as demonstrated in previous works [18]. Specifically, based on the polymer network theory, the number of entanglements of a chain is given by the ratio of the plateau modulus $GN0$ and the minimum of the loss modulus $Gmin\u2033$. The entanglement length, $le$, is evaluated in the context of a particular chain model (e.g., a flexible chain or a semiflexible chain) [19,20]. Thus, one can estimate the average chain length, $Lc$, as follows:

### B. Rolie-Poly model for monodisperse polymers

Before discussing the model for polydisperse polymers, we first take a moment to review the classic Rolie-Poly model of Likhtman and Graham [31], which has its origin in the full GLaMM model. The GLaMM model was originally proposed by Graham *et al.* [35] as a full chain molecular theory for entangled monodisperse linear polymer chains under fast deformation. The GLaMM model accounts for an accumulation of stress via affine deformation and relaxation of stress by reptation, chain retraction, and convective constraint release (CCR). Reptation was included in the original tube model introduced by Doi and Edwards [19] and describes the curvilinear motion of a polymer chain in a tube formed by entanglements with surrounding chains. When the deformation becomes much faster than the inverse of the Rouse time, the chain can become stretched, and as a result, the length of the chain and the occupied tube exceed their equilibrium configuration. When the strain stops, the chain retracts along the deformed tube until it regains its equilibrium contour. Convective constraint release (CCR) is the release of entanglement constraints due to chain retraction from affine deformation in nonlinear flows. The effect of CCR diminishes for shear rates larger than the inverse of the Rouse time.

Although the GLaMM model is successful in predicting the rheology of fast flows, it is computationally prohibitive for nonviscometric flow calculations and for OSR calculations that span several orders of magnitude of frequency and shear rate. Hence, as a simplified version of the GLaMM model, Likhtman and Graham derived a one-mode approximation of the GLaMM model to produce a differential constitutive model for entangled monodisperse polymer chains: the Rolie-Poly model (for Rouse linear entangled polymers) [31]. The time evolution of the contribution to the stress tensor due to the polymer, $\sigma $, is calculated via the conformation tensor for an entanglement segment, $A$, as

and

Here, $GN0$ is the plateau modulus, $\kappa $ is the velocity gradient tensor, and the stretch ratio of the polymer chains, i.e., the ratio of the current to equilibrium-with-no-flow chain contour lengths, is defined as

The reptation term (with $\tau d$) is a single time relaxation toward equilibrium. It is also possible to include the effect of contour length fluctuation (CLF) with the reptation term by modifying the reptation time [30,36]. However, as shown by Boudara *et al.* [30], the inclusion or absence of CLF only affects the quantitative agreement of the model with experiment without affecting the qualitative trends. To keep the physical picture simple, we will not include CLF in our calculations. The CCR term is also a relaxation toward equilibrium but with the rate dependent on the amount of stretch $\lambda $ and the Rouse time $\tau R$. Retraction is relaxation to zero stress, also with the rate dependent on $\lambda $ and $\tau R$. $\beta CCR$ specifies the CCR contribution and is analogous to $c\nu $ in the full GLaMM model. $c\nu $ determines the number of retraction events necessary to result in one tube hop of a tube diameter with a value of $c\nu \u22641$. [35] To fit the GLaMM model to experiments, the authors used $c\nu =0.1$, which approximately corresponds to $\beta CCR=1$. $\rho $ is an additional fitting parameter. $\beta CCR=1$, $\rho =\u22120.5$ are the optimal values to fit the Rolie-Poly model prediction to the full GLaMM theory, thus these are generally used in the Rolie-Poly model and also in this work.

The Rolie-Poly model incorporates changes in the conformation due to orientation and stretch in the same equation. Thus, the model avoids anomalous shear thickening behavior at shear rates corresponding to chain stretch predicted by models that decouple stretching and orientation [29]. Furthermore, since Rolie-Poly model is a single mode approximation of the full GLaMM model, it has sound microstructural basis while still making the calculation of OSR moduli computationally feasible.

### C. Rolie-Double-Poly (RDP) model for polydisperse polymers

To develop a constitutive model for polydisperse entangled polymer, Boudara and co-workers used ideas from both the Rolie-Poly model [31] and the double reptation approximation [37] to formulate the RDP model [30]. The “RDP model” is used to refer the “Rolie-Double-Poly” model, where “Double” signifies both “double reptation” and “double poly” as in “POLYdisperse POLYmers.” The RDP model can be formally derived as a greatly simplified approximation of the detailed molecular theory for bidisperse polymer melts of Read and co-workers [38].

In the RDP model, the total polymer contribution to the stress is the sum of the contributions coming from each species $i$, weighted by their volume fraction, $\varphi i$,

with

The conformation tensor $Ai$ accounts for the stresses that come from the interaction of the species $i$ with itself and the other $n\u22121$ species,

where $Aij$ is the stress conformation tensor on the $i$-chains coming from their entanglements with the $j$-chains,

In the original Rolie-Poly model, $\tau d$ specifies the terminal stress relaxation time of the polymer. However, in the RDP model, the effects of reptation are explicitly separated out from thermal constraint release (CR) given by the parameter $\beta th$. Due to the interaction of different chains, reptation of a smaller chain causes the constraint release of a longer chain. Thus, CR needs to appear in a separate term. To ensure that the RDP model converges to the Rolie-Poly model in the monodisperse limit, the reptation term and the CR term are divided by two. This choice is physically sound because the reptation term in the Rolie-Poly model is the effective reptation term that actually includes both reptation and CR.

## III. METHODS

### A. Numerical approach

#### 1. Monodisperse calculation

To evaluate the orthogonal superposition moduli, it is necessary to calculate the transient stress response under a combination of axial steady shear flow and orthogonal oscillation. The velocity gradient tensor is

where $\gamma \u02d9$ is the steady state shear rate and $\gamma \u02d9\u22a5(t)=\gamma 0,\u22a5\omega cos\omega t$ is the oscillatory superposition term. In this work, we kept $\gamma 0,\u22a5=0.05$ to remain in the linear viscoelastic region, as determined in the strain sweep (see supplementary material Fig. S1 [39]). The six chain conformation tensor components of Eq. (3) need to be calculated numerically by transiently solving the coupled differential equations. First, the equations are nondimensionlized using the Rouse time $\tau R$ and the plateau modulus $GN0$. Specifically, the Deborah number $De$ is defined as the dimensionless frequency ($De=\omega \u22c5\tau R=\omega ~$), the Weissenberg number is the dimensionless shear rate ($Wi=\gamma \u02d9\u22c5\tau R$), the dimensionless time is $t\u2032=t/\tau R$, and nondimensionalized stress tensor is $\sigma ~=\sigma /GN0$. Furthermore, the ratio of the reptation time and the Rouse time is $\tau d/\tau R=\theta =3Z$, where $Z$ is the number of entanglement points per chain. In this work for the monodisperse calculations, we take $Z=20$ for a well entangled polymer system.

The six transient equations for the conformation tensor components of Eq. (3) were solved using the MATLAB built-in solvers “ode45” and “ode15s.” “ode15s” is applied for the low frequency range to take advantage of the significant acceleration of the computation using adaptive mesh refinement. However, for the high frequency range, the adaptive mesh refinement caused large fluctuations in the moduli for adjacent frequency values. Consequently, for the high frequency range the very robust solver “ode45” is applied, which allows a user-defined time discretization. For this work, a logarithmic time spacing enables a sufficient sampling of the oscillations, across the entire frequency range simulated ($De=10\u22124$ to $100$).

The shear stress in the orthogonal direction $(\u22a5)$ is used for the determination of $G\u22a5\u2032$ and $G\u22a5\u2033$ using

where $\sigma ~0,\u22a5$ and $\gamma 0,\u22a5$ are the amplitudes of the dimensionless shear stress and the applied strain, whereas $\delta $ defines their phase shift. Furthermore, the term $\eta s\omega ~$ denotes the solvent contribution, which is purely viscous and, therefore, only appears in $G\u22a5\u2033$, with the scaled solvent viscosity $\eta s=0.001$ in this work.

To determine the values of $\delta $ and $\sigma ~0,\u22a5$ from the simulation, we first need to identify the start of the steady state in the main shear flow. The shear stress $\sigma ~xy$ is plotted as a function of time [Fig. 1(a)]. Within the plateau, the stress continues oscillating with a very small amplitude due to the superposed oscillation. We specify the plateau criterion as the smaller value of $10\u22123$ and $10\u22122\omega ~$ to account for the dependence of the amplitude of the stress oscillation on the frequency. The start of the plateau is defined as the first point in the transient shear stress with a moving standard deviation that is smaller than the specified plateau criterion. Once the time for the start of the plateau for $\sigma ~xy$ is specified, the same time is used to identify the start of steady state for the stress in the orthogonal direction $(\sigma ~yz)$, which is plotted in Fig. 1(b) along with the applied strain rate $(\gamma \u02d9\u22a5)$. The gray section of the figure illustrates the transient start-up before the harmonic oscillation. After neglecting the start-up phase, $\sigma ~yz$ and $\gamma \u02d9\u22a5$ are fitted by sine functions. Since the strain rate is a cosine function, it needs to be shifted by $\pi /2$ for the fitting of a sine function. Equations (12) and (13) define the two functions, including their fitting parameters $c1$, and $c2$,

In the equations, $\sigma ~0,\u22a5$ and $\gamma 0,\u22a5$ specify the amplitudes of the orthogonal stress and the applied strain in the orthogonal direction, respectively. Furthermore, the difference of the resulting coefficients defines the phase shift between the two wave functions: $\delta =c1\u2212c2$. With the phase shift and the amplitude values, the OSR moduli $G\u22a5\u2032$ and $G\u22a5\u2033$ are calculated for each combination of shear flow strength $Wi$ and superposed oscillation frequency $De$ according to Eqs. (10) and (11).

Finally, before using the numerical method above to make predictions of $G\u22a5\u2032$ and $G\u22a5\u2033$, we verified that $\gamma 0,\u22a5=0.05$ is small enough to be in the linear viscoelastic region by conducting a strain sweep at a fixed oscillation frequency [see supplementary material Fig. S1(a) [39]]. Furthermore, we validated the numerical method by comparing steady state predictions with predictions reported with the original Rolie-Poly model formulation [31] and comparing orthogonal superposition predictions using a Giesekus model with analytical results from the work by Kim and co-workers [5] (see supplementary material Fig. S2 [39]).

#### 2. Polydisperse calculation

The $G\u22a5\u2032$ and $G\u22a5\u2033$ calculations for polydisperse polymers follow the same workflow as the monodisperse calculation. The only difference is that several different conformation tensors $Aij$ need to be calculated before the summation, given by Eqs. (5) and (7), is computed. For the nondimensionalization, the Rouse time ($\tau R$) of the weight-average chain length ($Z\xaf=20$) is used to define the dimensionless shear rate $Wi$, the dimensionless frequency $\omega ~$, and the dimensionless time $t\u2032$. Furthermore, the continuous length distribution must be discretized into different species, each with its own chain length and volume fraction. However, the computational expense significantly increases with the number of species by a factor of $N2$. For this study, a maximum $N=32$ is used. To systematically study the effect of polydispersity on the OSR moduli, we assumed a lognormal distribution for the volume fraction distribution with varying degrees of polydispersity while keeping the weight-average molecular weight ($Mw$) the same. The lognormal distribution is frequently seen in polymers used for rheological studies [40]. The continuous distribution is defined as

with $s2=ln\u2061(Ip)$ and $\mu =ln\u2061(Z\xaf/Ip)$. $\varphi (Z)\delta Z$ gives the volume fraction of chains with entanglement number falling in the range $Z\xb1\delta Z/2$ for asymptotically small $\delta Z$. We also define $Z\xaf=\u222b\varphi ZdZ$ as the mean number of entanglements per chain and $Ip=Mw/Mn=Z\xaf\u222b\varphi (Z)/ZdZ$ as the polydispersity index. We then partition the continuous distribution into $N$ sections, where the $i=1,2,3,\u2026N$ section is bracketed by chains with entanglement numbers $Z^i\u22121$ and $Z^i$. In the approach used in this work and first developed by Peterson [41,42], the partitions are drawn at lines of equal volume fraction. In our view, this is the simplest discretization scheme that guarantees convergence to the correct cumulative distribution when $N\u2192\u221e$. Written in terms of the cumulative distribution function $\Phi 0(Z)$, the distribution scheme is formally given by

For the lognormal distribution, $\Phi 0(Z)$ can be inverted to obtain the set of $Z^i$,

Now defining $Zi$ as the weight average molecular weight of chains in the range between $Z^i\u22121$ and $Z^i$, we approximate the continuous distribution by a discrete distribution of $N$ components with entanglement numbers $Z1,Z2,Z3,\u2026,ZN$ and equal volume fractions $\varphi 1=\varphi 2=\varphi 3=\cdots =1/N$. This strategy ensures that the mean entanglement number $Z\xaf=\u2211i\varphi iZi$ is not changed by discretization. The end result for the discretization scheme is

### B. Perturbation analysis

For orthogonal superposition, the oscillatory flow is assumed to be in the linear viscoelastic regime. Thus, the strain amplitude $\gamma 0,\u22a5\u226a1$ and $\gamma 0,\u22a5\omega ~\u226aWi$. $Wi$ is the shear rate of the mean flow nondimensionalized by the Rouse time $\tau R$. In this case, we can use a perturbation analysis with a small parameter $\u03f5$ being the ratio of the magnitudes of the strain amplitudes of the orthogonal and main flow: $\u03f5=\gamma 0,\u22a5|\gamma \u02d9|$. We shall assume a regular asymptotic solution structure and retain only the first order perturbation to the steady state solution. The superscript “(0)” denotes the constant, order-one steady state solution and the superscript “(1)” denotes the order-$\u03f5$ perturbation to the steady state value. With these notations, the stress tensor components and the chain stretch can be expressed as an expansion about the steady state levels: $\sigma ij(t)=\sigma ij(0)+\u03f5\sigma ij(1)(t)$ and $\lambda (t)=\lambda (0)+\u03f5\lambda (1)(t)$. Furthermore, as in the numerical calculation, we nondimensionalize the stress by the plateau modulus $GN0$: $\sigma ~ij=\sigma ij/GN0$, which also equals the conformation tensor components $Aij$ according to Eq. (2). Thus, we have $Aij(t)=Aij(0)+\u03f5Aij(1)(t)$. We will first write the mathematical framework for the perturbation analysis using the monodisperse Rolie-Poly model. A similar procedure is used for the RDP model. For all six equations of the stress components in Eq. (3), we collect the $O(1)$ terms $(0)$ and the $O(\u03f5)$ terms $(1)$. For the purpose of calculating the OSR moduli, we only need to solve the $O(\u03f5)$ equation for $Ayz(1)$,

We can write this equation completely in terms of $\lambda (0)$ by obtaining the steady state solution for $Ayy(0)$,

The equation for $Ayz(1)$ can then be rewritten as

with

Equation (21) can be solved analytically for the OSR storage and loss moduli,

Here, only the polymeric contribution is included. To include the solvent contribution, we can add $\eta s\omega ~$ as in Eq. (11). In Eqs. (23) and (24), $Ayy(0)(\gamma \u02d9)$ and $B(\gamma \u02d9)$ are evaluated numerically from steady state Rolie-Poly equations without superposition. The values are then used in the equations to calculate $G\u22a5\u2032$ and $G\u22a5\u2033$ analytically. The perturbation analysis gives a speedup of more than 100 times while offering the same accuracy compared to the numerical calculation of OSR moduli (see supplementary material Fig. S3 [39] for a comparison of the two results).

A similar procedure is used for the perturbation analysis using the RDP model for polydisperse polymers. In this case, $Aijyy(0)(\gamma \u02d9)$ and $Bij(\gamma \u02d9)$ of each entanglement pair need to be determined from the steady state calculation. The OSR moduli are then sums of contributions from all entanglement pairs,

At this point, we reflect on the way in which we derived these relationships and their potential generality. Results for the orthogonal moduli in terms of the $yy$-component of the polymer conformation tensor come about from the convection of the normal stress through the upper-convected time derivative to the $yz$-component of the stress, which is the direction of the oscillation. Regardless of the constitutive model that is used, this is the procedure and the mathematical propagation that gives us a relationship similar to Eqs. (25) and (26). In fact, the linear proportionality between the plateau value of $G\u22a5\u2032$ and $Aijyy(0)(\gamma \u02d9)$ not only holds for Rolie-Poly and RDP models but also for any model with the upper-convected time derivative of the conformation tensor $A$ and stress relaxation terms that are linear in the conformation tensor and any linear superposition of these models. The results are valid, regardless of how the chain conformation is quantified, i.e., whether the end-to-end vector or the orientation tensor is used. A similar analysis on the DCR-CS model by Marrucci and Ianniruberto [24] (see supplementary material Sec. V [39]) and the MLD model by Mead [25] produces the same relationship.

## IV. RESULTS

### A. Monodisperse case

We first present numerical calculations to obtain the orthogonal superposition moduli using the Rolie-Poly model. Figure 2 shows the computational results of a frequency $De$ sweep for the oscillatory shear at different shear rates of the steady flow, i.e., Weissenberg numbers, $Wi=\gamma \u02d9\u22c5\tau R$, ranging from 1 to 100. Initial studies of OSR on entangled polymers were limited to relatively low $Wi$, but recent advances in instrumentation have extended the accessible range to cover the range of $Wi$ we study. These results are for a monodisperse polymer with an entanglement number of $Z=20$, which is the same as what accompanied the original Rolie-Poly model formulation [31]. Consequently, the ratio of the two relaxation time scales $\theta =\tau d/\tau R=60$. We find two distinctive trends in the moduli behavior with increasing $Wi$. First of all, a horizontal shift is observed, which can be seen best in the rightwards shifting of the $G\u22a5\u2032$–$G\u22a5\u2033$ crossover point. Second, $G\u22a5\u2032$ decreases monotonically with increasing main flow shear rates $Wi$. This effect can be seen best by the plateau values of the $G\u22a5\u2032$ curves. Comparing $G\u22a5\u2032$ and $G\u22a5\u2033$ at $Wi=0.001$ to the linear viscoelastic moduli $G\u2032$ and $G\u2033$ in the absence of flow [Fig. S1(b) in supplementary materials [39]] verifies that results at low flow rates approximate the result at equilibrium.

Furthermore, evaluating the contribution of each stress relaxation mechanism (reptation, retraction and CCR) in the $yz$-component of the stress, $\sigma ~yz$, provides a better understanding of the causes of the shifts in the OSR moduli. Figure 3 shows the fractional contribution of relaxed orthogonal stress terms separated by relaxation mechanism, at a frequency of $De=0.001$. The shape of the plots is independent of the choice of the frequency. As $Wi$ increases, contributions of the stretch-based phenomena (retraction and CCR) increase, with CCR developing a maximum and then decreasing. In contrast, reptation starts at 100% for very low shear rates and decreases with $Wi$ to a value close to 0%. This is because the stress relaxed by reptation approaches a constant when the inverse of the disentanglement time $\tau d$ is reached. Consequently, the reptation contribution in the Rolie-Poly model remains constant, while retraction and CCR start dominating relaxation for increasing $Wi$. In the first part of their increase, the same slope is observed. However, as $Wi$ continues to increase the chains can retract only partially from affine stretch, and convective constraint release is less effective and the relative magnitude of the CCR contribution decreases. At the same frequencies, the retraction term develops a shoulder. Finally, at sufficiently large $Wi$, retraction dominates the relaxed stress with a value close to 100%.

By examining the numerical results alone, it is difficult to pinpoint what physical phenomena directly cause the horizontal and vertical shifts. If these were purely linear viscoelastic measurements under different temperatures and concentrations, researchers would typically perform a time-temperature superposition or time-concentration superposition and force all of these curves to collapse. Here, we show that we can achieve such a collapse with shift factors. Furthermore, these shift factors can be predicted based on the perturbation analysis results in Eqs. (23) and (24). By assuming the superposed oscillation is a small perturbation to the main shear flow ($\gamma 0,\u22a5\u226a1$ and $\gamma 0,\u22a5De\u226aWi$), we obtained analytical expressions for $G\u22a5\u2032$ and $G\u22a5\u2033$ in Eqs. (23) and (24). Using the two parameters $Ayy(0)(\gamma \u02d9)$ and $B(\gamma \u02d9)$ identified in the analytical expressions, we scaled the numerical results from Fig. 2 after subtracting out the solvent contribution ($\eta s\omega ~$). This procedure completely collapses all the $Wi$-dependent moduli onto a single set of curves. Specifically, we divided the OSR moduli by $Ayy(0)(\gamma \u02d9)$ and divided the dimensionless frequency by $B(\gamma \u02d9)$ [Fig. 4(a)]. The collapse of the curves suggests that the steady shear flow only speeds up the relaxation and decreases the effective oscillatory driving force, without changing the relaxation mechanism.

We can further quantify the magnitude of the vertical shift by plotting the plateau value of dimensionless $G\u22a5\u2032$ as a function of $Wi$ [Fig. 4(b)]. In fact, in the plateau, we have $lim\omega ~\u226b1G\u22a5\u2032GN0=Ayy(0)(\gamma \u02d9)$. With increasing $Wi$, $G\u22a5\u2032$ monotonically decreases. Around $Wi=1$, a shoulder shape appears in the plateau [also in $Ayy(0)(\gamma \u02d9)$], similar to the shoulder shape in the retraction contribution in Fig. 3. This shoulder again arises due to the competition of retraction and convective constraint release. Since $Wi$ is defined in terms of the Rouse time $\tau R$, $Wi=1$ represents the start of significant chain stretch, which results in increased contribution of the stress relaxation by chain retraction. We can also quantify the horizontal shift by plotting the frequency at the crossover of $G\u22a5\u2032$ and $G\u22a5\u2033$ as a function of $Wi$ [Fig. 4(c)]. With increasing $Wi$, the crossover frequency monotonically increases. In linear viscoelastic measurements, the inverse of the crossover frequency gives a characteristic relaxation time of the system: $\tau =1/\omega c$. So, Fig. 4(c) indicates the effective relaxation time of the polymers is decreasing with increasing $Wi$. By examining Eqs. (23) and (24), we find the crossover frequency equals $B$. To get a better understanding of what $Ayy(0)(\gamma \u02d9)$ and $B(\gamma \u02d9)$ represent physically, we can return the $G\u22a5\u2032$ and $G\u22a5\u2033$ results in Eqs. (23) and (24) to dimensional form,

with

The effective relaxation time $\tau eff$ represents the speedup of the stress relaxation in the flow. Its asymptotic value changes from $\tau eff=\tau d$ at low $Wi$ before the onset of chain stretch to $\tau eff=\tau R/2$ for very large chain stretch. So, $Ayy(0)(\gamma \u02d9)$ represents the nonlinear shear-dependent chain conformation in the $yy$-direction, which is orthogonal to the main flow direction, and $B(\gamma \u02d9)$ represents the nonlinear shear-dependent relaxation rate.

As discussed in Sec. II A, for conventional linear viscoelastic measurements, the plateau modulus $GN0$ is used to estimate the entanglement length while the ratio of the plateau modulus and $Gmin\u2033$ is used to estimate the number of entanglements. Thus, LVE $G\u2032$ and $G\u2033$ moduli are used to estimate average chain length, and previous investigators have assumed that the corresponding OSR moduli can be interpreted in the same way, thus providing a way to explore the effect of a steady shear flow on chain scission for wormlike micelles and other “living” polymers [5,43]. However, we see that the decrease in $G\u22a5\u2032$ with increasing $Wi$ in Fig. 4(b) and the decrease in the ratio of the plateau and $Gmin\u2033$ in Fig. 2 occur without any modification to the properties of the chain; for example, those that might be extracted from a scaling analysis of the linear viscoelasticity (i.e., the number of entanglements is fixed at 20 per chain). The plateau value and the crossover frequency change due to chain stretch and orientation, as illustrated in the contribution of stress relaxation in Fig. 3. Therefore, changes in the plateau modulus and relaxation time cannot be interpreted simply as due to changes in the scaling variables (e.g., entanglement length, entanglement density, etc.). This is consistent with a recent Brownian dynamics simulation for star polymers [44], in which Metri and Briels showed that the shear relaxation modulus from orthogonal superposition is not the same as the stress relaxation modulus from step strain (as is the case for linear viscoelastic experiment at equilibrium), even for weak shear flow in the OSR experiment.

Although the data interpretation from the OSR measurement is complicated even for a monodisperse polymer, the perturbation analysis nevertheless offers physical insights to pinpoint the material properties that cause the change in the OSR moduli. Before moving to the polydisperse calculations, we first demonstrate the generality of the perturbation analysis by applying it to a different model, the DCR-CS (double-convection-reptation model with chain stretch) by Marrucci and Ianniruberto [24]. As shown in supplementary material Sec. V [39], the resulting expressions for the OSR moduli are the same as those from the Rolie-Poly model [Eqs. (23) and (24)]. The differences are in the expressions for $Ayy(0)(\gamma \u02d9)$ and $B(\gamma \u02d9)$. Even though these expressions are different, the similarity in the functional dependence of $G\u22a5\u2032$ and $G\u22a5\u2033$ on $Ayy(0)(\gamma \u02d9)$ and $B(\gamma \u02d9)$ indicates that the dependence of the plateau shift on the $yy$-component of the chain conformation and the horizontal shift on the effective relaxation time are model-independent and represent the physical phenomena that determine the OSR moduli for linear polymers of fixed chain length.

### B. Polydisperse case: A lognormal distribution

To systematically study the effect of polydispersity on the OSR moduli, we used a discretized lognormal molecular weight distribution, as discussed in Sec. III A 2. Specifically, we used 32 bins to represent the continuous MW distribution with equal volume fraction in each bin. The polydispersity index was varied from 1.01, which is nearly monodisperse, to 5, which has a long tail. The weight-average molecular weight (weight-average number of entanglements per chain, $Z\xaf$) is fixed at 20. Figure 5(a) plots all of the molecular weight distributions used in this work. First, we compare the numerical calculations at two representative cases of the degree of polydispersity: $PDI=1.01$ and $PDI=2$ in Figs. 5(b) and 5(c), respectively. For $PDI=1.01$, the results are nearly identical to the monodisperse results in Fig. 2. For $PDI=2$, the effects of polydispersity are more evident. Specifically, the shoulder in the storage modulus becomes broader at all $Wi$ and the plateau value decreases less. The gradual evolution of the polydisperse moduli reflects the broad relaxation spectrum in the highly polydisperse system (i.e., many relaxation times). As we will show below, the weakened softening of the plateau for the polydisperse case compared to the monodisperse case arises from the nonlinear dependence of $Aijyy(0)$ on the chain length $Zi$. The presence of the short chains causes the plateau to decrease less than the monodisperse case.

To better quantify the horizontal and vertical shifts in the OSR moduli, we can plot the plateau value and the location of the apparent crossover frequency as a function of $Wi$ for different $PDI$. However, having to account for $N2$ number of chain pairs for all combinations of $De$ and $Wi$ significantly increases the computational time for the polydisperse case compared to the monodisperse case. Thus, once again, we turn to a perturbation analysis to obtain analytical expressions for the OSR moduli, which makes prediction of OSR moduli at finely resolved $Wi$ and $De$ possible. The procedure for the perturbation analysis and the main results are described in Sec. III B. We verified that the perturbation analysis produced identical results as the full numerical calculations. So, subsequently, we only report results from the perturbation analysis.

Figure 6 quantifies changes in the OSR moduli and crossover frequency as a function of $Wi$ for varying degrees of polydispersity: $PDI=1.01,1.05,1.2,2,5$ and for the monodisperse result (shown as the solid black curve). The crossover frequency $G\u22a5\u2032$ and $G\u22a5\u2033$ increases monotonically with increasing $Wi$ for all cases except for $PDI=5$ [Fig. 6(a)]. For $PDI=5$ and $Wi>10$, the decrease of the crossover frequency at large $Wi$ is a result of significant stretching of the longest chains in the system. Such extreme stretching is unphysical for finitely extensible chains and can be removed with the inclusion of a Warner spring [51] (see supplementary material Fig. S7(a) [39]). Figure S7 compares the linear spring calculation with the Warner spring calculation. Subfigures (a)–(d) are organized similarly to those in Fig. 6 here. Predictions using the Warner spring only start to differ significantly from predictions for the linear spring for $Wi>10$ and $PDI\u22652$. On the one hand, across most of the parameter space (spanning different shear rates and $PDI$) that is accessible experimentally, the predictions for the linear and Warner spring are very similar. On the other hand, the Warner spring adds extra complications associated with an additional parameter $\lambda max$, the maximum chain stretch. Thus, we will only discuss results for the linear spring in the main paper. At low dispersity, the moduli value at the crossover decreases monotonically with $Wi$ [Fig. 6(b)]. However, for $PDI=2$ and $5$, the moduli at the crossover first increase then decrease with $Wi$. The $G\u22a5\u2032$ plateau decreases monotonically with increasing $Wi$ for all cases [Fig. 6(c)]. With increasing $PDI$, the curves deviate more and more from the monodisperse curve. Additionally, the shoulder shape in the monodisperse result that arises from the competition between retraction and CCR becomes smoothed out with increasing polydispersity. Interestingly, the monodisperse result for the moduli at the crossover has the same shape as the monodisperse plateau value, with the former being half of the magnitude. This effect is better illustrated in Fig. 6(d), in which the ratio of the crossover moduli and the plateau moduli is plotted as a function of $Wi$ for the monodisperse case and all polydisperse cases. The monodisperse case is a flat line at 0.5 for all $Wi$ which indicates that for a monodisperse polymer, the moduli value at the crossover is exactly half of the plateau value for all $Wi$. Examining the perturbation result for the monodisperse case [Eqs. (23) and (24)] illustrates this point more clearly. At the crossover, $G\u22a5\u2032=G\u22a5\u2033$ and $B(\gamma \u02d9)=\omega ~$. Thus, the moduli value at the crossover equals $Ayy(0)(\gamma \u02d9)/2$. However, the $Wi$-independent, constant moduli ratio is not preserved in the polydisperse cases. With increasing polydispersity, the ratio decreases and becomes more nonmonotonic [Fig. 6(d)].

The reason for all the changes in the polydisperse results compared to the monodisperse result, as illustrated in Figs. 6(a)–6(d), arises from the fact that the OSR moduli for the polydisperse case are summations of the contributions from all entanglement pairs and the OSR moduli are affected by the interactions between chains of different length [Eqs. (25) and (26)]. In particular, instead of being solely determined by the $Ayy(0)(\gamma \u02d9)$ component of a single chain length as in the monodisperse case, the polydisperse plateau modulus is determined by $\u2211i=1n\u2211j=1n\varphi i\varphi jAijyy(0)(\gamma \u02d9)=Ayy(0)(\gamma \u02d9)$, which represents the dimensionless weight-average chain conformation in the $yy$-direction of the entire population of chains. However, a weight-average of $Bij(\gamma \u02d9)$ for all entanglement pairs does not yield the overall effective relaxation rate (i.e., the crossover frequency) for the polydisperse case. This is because $Bij(\gamma \u02d9)$ appears in the denominator of the individual moduli terms that are summed together. Thus, for polydisperse systems, although there is a qualitative correlation of the data with the degree of polydispersity, as shown in Fig. 6 especially in (d) for the ratio of the moduli, no simple scaling exists to explain the shifts in the moduli, as in the monodisperse case.

One may wonder to what extent the changes in the OSR moduli in the polydisperse case can be understood by simple superposition of results for the discrete chains as opposed to needing the full double reptation physics included in the RDP model. Here, we make a comparison of the RDP results for $PDI=2$ shown in Fig. 6 and the results of the simple superposition of the monodisperse results for each discrete species in the polydisperse calculation. Thus, for the linear superposition case, we have $G\u22a5\u2032GN0=\u2211i=1n\varphi iAiyy(0)\omega ~2Bi2+\omega ~2$ and $G\u22a5\u2033GN0=\u2211i=1n\varphi iAiyy(0)Bi\omega ~Bi2+\omega ~2$, where $Aiyy(0)(\gamma \u02d9)$ and $Bi(\gamma \u02d9)$ are calculated from the perturbation analysis for the monodisperse case for chains of length $Zi$ and the values for $Zi$ are obtained from the discretization of the lognormal distribution. The subfigures (a)–(d) in Fig. 7 are organized similarly to the subfigures in Fig. 6. Although the two sets of predictions share similar qualitative trends, the simple superposition cannot quantitatively capture the changes in the crossover frequency and crossover moduli at low $Wi$. Thus, the full double reptation physics is needed to quantitatively describe the OSR moduli. We also note that the amount of discrepancy between the simple superposition results and the RDP results increases with increasing $PDI$. Nevertheless, since the plateau value at large $Wi$ ($Wi=100$) is nearly identical for the simple superposition and the RDP model, the simple superposition results can be used to explain the weakened softening of the plateau value for the polydisperse case compared to the monodisperse case, as observed in Figs. 5(c) and 6(c). As shown in supplementary material Fig. S7 [39], the $yy$-component of the conformation tensor $Aiyy(0)$ depends nonlinearly on the chain length $Zi$. Thus, for a discretized distribution of the lognormal distribution with a weight-average the same as the chain length in the monodisperse case, the weight-average of $Aiyy(0)$ is not the same as $Ayy(0)$ for the monodisperse case.

In conclusion, the presence of chains of different length results in the smoother transition to the plateau. The presence of short chains and the nonlinear dependence of $Aijyy(0)$ on the chain length results in the weakened softening of the plateau compared to the monodisperse case. Although a simple superposition of the monodisperse results can capture the qualitative trends observed in the polydisperse OSR moduli, to obtain quantitative predictions, the double reptation physics is needed.

## V. DISCUSSION

Given the results just presented, we will now provide perspective on how this work could potentially help aid the design and interpretation of orthogonal superposition experiments.

### A. An analogy to Laun’s rule

By using a combination of numerical calculations and perturbation analysis, we showed for both the monodisperse case and the polydisperse case, the plateau of the orthogonal modulus is completely determined by the $yy$-component of the polymer conformation in the steady imposed shear, namely,

where $Ayy(0)(\gamma \u02d9)=\u2211i=1n\u2211j=1n\varphi i\varphi jAijyy(0)(\gamma \u02d9)$ is the weight-averaged result for the polydisperse case. As we discussed in the theory Sec. III B, Eq. (29) is generally valid for any model with the upper-convected derivative of the conformation tensor $A$ and relaxation terms that are linear in the conformation tensor or any linear superposition of these models. Given this result and its generality, we can use orthogonal superposition measurements to obtain estimates of $Ayy(0)(\gamma \u02d9)$. Experimentally, this is accomplished by shifting $G\u22a5\u2032$ and $G\u22a5\u2033$ at different shear rates to achieve a complete overlap, in the same manner performed routinely for linear viscoelastic measurements for the purpose of time-temperature or time-concentration superposition, one can obtain shift factors $\alpha (\gamma \u02d9)=G\u22a5\u2032(\omega ,\gamma \u02d9)G\u22a5\u2032(\omega ,\gamma \u02d9=0)=G\u22a5\u2033(\omega ,\gamma \u02d9)G\u22a5\u2033(\omega ,\gamma \u02d9=0)$. We can use the vertical shift factors to obtain $Ayy(0)(\gamma \u02d9)$:

It is worth noting here that in time-temperature superposition, the conventional nomenclature for the vertical shift factor is $b$. Since we already used $B$ to represent the dimensionless crossover frequency, which is related to the horizontal shift, we are using $\alpha $ to represent the vertical shift factor.

Since a rheological constitutive equation for an incompressible fluid only specifies the stress to within an arbitrary isotropic component, we cannot use $Ayy(0)(\gamma \u02d9)$ to directly determine the normal stress component in the $yy$-direction, $\sigma yy$. However, we can attempt to make predictions of the first normal stress difference, $N1$. In particular, for the Rolie-Poly model, we can solve for the steady state value of $Axx(0)(\gamma \u02d9)$ using the $O(1)$ equation in the perturbation analysis (see Sec. VIII in supplementary material for the detailed derivation [39]). Then, using results for $Ayy(0)(\gamma \u02d9)$ and $B(\gamma \u02d9)$, we obtain

If we recall that $Wi=\gamma \u02d9\u22c5\tau R$ and $B(\gamma \u02d9)$ represents the dimensionless crossover frequency, $B(\gamma \u02d9)=\omega ~c,\gamma \u02d9=\omega c,\gamma \u02d9\u22c5\tau R$, we get

Experimentally, this means we can use the plateau value of $G\u22a5\u2032(\gamma \u02d9)$, i.e., $GN0Ayy(0)(\gamma \u02d9)$, the applied shear rate, $\gamma \u02d9$, and the shear-rate dependent crossover frequency, $\omega c,\gamma \u02d9$, to completely determine the first normal stress difference $N1(\gamma \u02d9)$. In the limit where a single-mode model accurately describes the experimental data, one can use this relation to determine $N1(\gamma \u02d9)$. We can further rearrange Eq. (32) to get

Previously, several relations have been proposed to connect the first normal stress difference or the $yy$-component of the stress to more easily measurable quantities in shear flow. Specifically, Holroyd *et al.* [45] derived analytical solutions of the Rolie-Poly model and found an expression to relate the recoverable strain to the $yy$-component of the stress, $\gamma rec=\sigma xy0/\sigma yy0$. This theoretical result has been verified in molecular dynamics simulations of entangled polymers.[46] In LAOS measurements, Lee *et al.* found $N1(t)=2G(t)\gamma rec2(t)$, where $G(t)$ is the recoverable strain-dependent elastic modulus and $\gamma rec$ is the recoverable strain [47]. Tanner and Williams [48] used the phenomenological K-BKZ model [49,50] to derive a relationship between the OSR moduli and the first normal stress difference $lim\omega \u21920[G\u22a5\u2032(\omega ,\gamma \u02d9)\omega 2]=N1(\gamma \u02d9)2\gamma \u02d92$. This relation only works in the low-frequency limit and is specific to the K-BKZ model. The well-known Laun’s rule [32] relates the first normal stress difference to the linear viscoelastic moduli and is predicted to work for any general viscoelastic model,

Our result is an alternative set of relations that may prove more useful depending on where the sensitivity of the measurement lies. Previous rules provide relationships that are generally valid at asymptotically small frequencies and shear rates. We can now make a similar comparison at asymptotically large frequencies, which is usually the region with better instrumental sensitivity. Notably, the normal stress component is obtained as a by-product of measuring the plateau of the orthogonal modulus without having to build extra instrumentation to measure the normal stress.

### B. Determination of constitutive parameters

From the perturbation analysis, we found the horizontal shift of the OSR moduli is determined by the effective relaxation rate $B(\gamma \u02d9)$. In the context of the Rolie-Poly model (i.e., for monodisperse polymers), since both $Ayy(0)(\gamma \u02d9)$ and $B(\gamma \u02d9)$ only depend on the ratio of relaxation time $\theta $, the CCR parameter $\beta CCR$, and chain stretch due to the steady shear flow $\lambda (0)$, OSR experiments can be used to estimate chain stretch and the CCR parameter $\beta CCR$. Specifically, one can first use linear viscoelastic measurements to estimate the Rouse time $\tau R$ and the reptation time $\tau d$. Then, one can subtract the solvent contribution from the orthogonal superposition measurement (in the values for $G\u22a5\u2033$). Next, $G\u22a5\u2032$ and $G\u22a5\u2033$ at different shear rates are shifted numerically to achieve a complete overlap with associated shift factors $\alpha (\gamma \u02d9)$ and $B(\gamma \u02d9)$. Finally, we can use the crossover frequencies to obtain the relaxation rate $B(\gamma \u02d9)$. For the simplest case of $\beta CCR=1$ as we assumed in the previous calculations, we can obtain chain stretch $\lambda (\gamma \u02d9)$ from the crossover frequency

In order to obtain both $\lambda (\gamma \u02d9)$ and $\beta CCR$, we need to use expressions for $Ayy(0)(\gamma \u02d9)$ and $B(\gamma \u02d9)$ that include $\beta CCR$,

By fitting the shear rate dependent experimental values of $Ayy(0)(\gamma \u02d9)$ and $B(\gamma \u02d9)$ to Eqs. (35) and (36), we can obtain $\lambda (\gamma \u02d9)$ and $\beta CCR$. Furthermore, using Eq. (30) and converting $B(\gamma \u02d9)$ to dimensional form [$B(\gamma \u02d9)=\omega c,\gamma \u02d9\tau R$], we obtain expressions for $\lambda (\gamma \u02d9)$ and $\beta CCR$ in terms of the vertical shifting factor $\alpha (\gamma \u02d9)$ and crossover frequency measured in experiments,

and

Typically, the CCR parameter is assumed to be a constant and independent of shear rate. In the case that model predictions are used to generate the OSR moduli, as in the present analysis, Eq. (38) can be used to test the self-consistency of $\beta CCR$. More interestingly, in the case of experimental data where the apparent value(s) of the nonlinear constitutive model parameter(s) is unknown and in the limit where a single mode model accurately describes the experimental data, one can use this shifting procedure in order to extract and fit model parameters of the desired constitutive equation. In order to successfully estimate model parameters from OSR experiments, we, therefore, give the following recommendations: (1) The shear rates need to span a large enough range to cover both the case of no chain stretch and the case of significant chain stretch ($\gamma \u02d9\tau R>1$). (2) The frequency range needs to span from values below the crossover frequency at equilibrium to frequency values corresponding to the plateau of $G\u22a5\u2032$ in flow.

### C. Physical interpretation of OSR moduli

We now discuss the implications of our analysis on the more general issue of the physical interpretation of OSR moduli measured in experiments in the absence of a particular constitutive model. In light of Fig. 6, the plateau modulus for both the monodisperse and polydisperse cases decreases with increasing $Wi$ and the moduli shifts horizontally without any modification to the properties of the chain; e.g., those that might be extracted from a scaling analysis of the linear viscoelasticity. For example, in previous orthogonal superposition experiments on wormlike micelles by Vermant and co-workers [5] and Rothstein and co-workers [43], the authors assumed by analogy with linear viscoelasticity $\xi m(\gamma \u02d9)=(kBT/G\u22a5\u2032(\gamma \u02d9)|\omega \u2192\u221e)1/3$ to estimate the apparent, shear-rate dependent mesh size, $\xi m$, from the plateau modulus of orthogonal superposition moduli, $G\u22a5\u2032(\gamma \u02d9)|\omega \u2192\u221e$. Vermant and co-workers linked the systematic decrease of the plateau modulus to shear-induced changes in the mesh size of the micellar network. Rothstein and co-workers went further to suggest the change in the mesh size is related to the breakage of micelles in flow. However, our analysis of OSR for the Rolie-Poly model clearly shows, even in the context of regular polymers, the plateau modulus decreases with increasing shear flow due to stretch-based relaxation mechanisms (CCR and retraction), in the absence of any chain breakage. A comparison between the experiments on WLMs by Rothstein and co-workers and our RDP model calculation with a lognormal distribution of $PDI=2$ shows a similar change in the plateau modulus $(\u223c20%)$ over the comparable range of $Wi$, suggesting that changes in entanglement dynamics are important and must be considered in this range. Additionally, the effect of polydispersity adds another complication for any analysis that tries to extract chain length in flow from the plateau of OSR moduli. As such, one cannot simply use empirical relationships for estimating length that are developed for equilibrium linear viscoelastic measurements, such as those described in Sec. II A, to interpret OSR results. This knowledge from simulations is very useful for interpreting OSR experiments.

In addition to the applicability of linear viscoelastic scaling analysis, one may wonder about the potential for using OSR moduli to infer or understand polydispersity effects given the present analysis. We identified new parameterization of the OSR moduli to qualitatively relate to the degree of polydispersity. Specifically, the ratio of the moduli at the crossover to the plateau modulus as a function of shear rate becomes increasingly nonmonotonic and deviates from the constant value of 0.5 with increasing degree of polydispersity [Fig. 6(d)]. Although Fig. 6 focuses on the result for the lognormal distribution, the qualitative trend also holds for other continuous, broad distributions, except for a bi-modal distribution with sharp peaks. In particular, for $PDI=1.01$, the shape of the molecular weight distribution is nearly identical for lognormal and Gaussian; similarly, for $PDI=2$, the shape of the molecular weight distribution is very similar for the lognormal and exponential. Although in our analysis we primarily focused on the lognormal distribution, by changing the degree of polydispersity, we are also changing the shape of the molecular weight distribution.

Conventionally, viscoelastic measurements of polydisperse systems have been used to extract information about the distribution of relaxation times, i.e., the discrete linear response relaxation spectrum. In Figs. 7 and S9 [39], we show that the nonlinear decomposition of the orthogonal modulus into its relative contributions $Aijyy(\gamma \u02d9)$ is sensitive to the spectrum of the nonlinear relaxation times. In this way, these measurements could be used as a nonlinear analog to linear viscoelastic measurements in terms of obtaining the spectra of nonlinear relaxation rates. The idea of extracting a shear rate-dependent relaxation spectrum from OSR measurement was proposed by Yamamoto [21]. However, this approach has a drawback because one has to assume a functional form of the relaxation spectrum beforehand. For the Rolie-Poly and RDP models, we can analytically derive the results for the superposition moduli that naturally contain a sum of the individual chain pairs in the system [Eqs. (25) and (26)]. In this case, the relaxation rate of individual modes is proportional to $Aijyy(\gamma \u02d9)$. Thus, if the experimental system can be accurately described by these relationships, then one can fit experimental data to these two equations and just leave $Aijyy(\gamma \u02d9)$ and $Bij(\gamma \u02d9)$ as adjustable parameters for some number of modes. Then, what one measures is the discrete nonlinear relaxation spectrum and the resulting $Aijyy(\gamma \u02d9)$ and $Bij(\gamma \u02d9)$ in the nonlinear case are equivalent to $G$ and $\tau $ for the multimode Maxwell model in the linear viscoelastic limit. The advantage of this approach is that we do not need to presuppose anything about what $Aijyy(\gamma \u02d9)$ and $Bij(\gamma \u02d9)$ are, as long as a general multimode model with the upper-convected time derivative of the conformation tensor and relaxation terms linear in the conformation tensor describes the material of interest. For the specific case of the RDP model, we actually know what $Aijyy(\gamma \u02d9)$ and $Bij(\gamma \u02d9)$ are, whereas if we presuppose a form for the relaxation spectrum, we would not necessarily know what $Aijyy(\gamma \u02d9)$ and $Bij(\gamma \u02d9)$ are. In other words, our results reaffirm the approach by Yamamoto that one can define some effective nonlinear relaxation spectrum, but beyond this, it gives proof to the idea that the relaxation spectrum should be related to the underlying constitutive behavior of the fluid. Going a step further, for the case of the Rolie-Poly and the RDP models, we can obtain analytical expressions for the discrete nonlinear relaxation spectrum from the expressions for the OSR moduli [Eqs. (25) and (26)].

To summarize the previous discussions, our study emphasizes the critical importance of comparing OSR experiments with this type of detailed OSR predictions using rheological models to understand the physical meaning of changes in orthogonal superposition moduli as is now routinely done for other types of nonlinear rheological measurement protocols.

## VI. CONCLUSIONS

In this work, we combine numerical calculations and a perturbation analysis using detailed microstructural models (Rolie-Poly and Rolie-Double-Poly models) to study orthogonal superposition for monodisperse and polydisperse entangled linear polymers. A perturbation analysis allows identification of key parameters ($Ayy(0)(\gamma \u02d9)$ and $B(\gamma \u02d9)$) that determine the shifts of the OSR moduli. By quantifying the shifts in an OSR experiment and using the expressions for $Ayy(0)(\gamma \u02d9)$ and $B(\gamma \u02d9)$, we can estimate the average chain stretch $\lambda (\gamma \u02d9)$, the CCR parameter $\beta CCR$, and the first normal stress difference $N1$.

By systematically varying the degrees of polydispersity while keeping the weight-average molecular weight constant in a lognormal distribution of the molecular weight, we find qualitative signatures in the OSR moduli that correlate with the degree of polydispersity. Furthermore, the presence of short chains and the nonlinear dependence of $Aijyy(0)$ on the chain length results in the weakened softening of the plateau compared to the monodisperse case. Although a simple superposition of the monodisperse results can capture the qualitative trends observed in the polydisperse OSR moduli, to obtain quantitative predictions, the double reptation physics is needed. By making an analogy with the discrete relaxation spectrum obtained from LVE moduli, we showed that OSR moduli can be used to extract a nonlinear relaxation spectrum, which reveals a potential significant advantage of OSR measurement over conventional rheological measurements.

We also found that the contributions of different relaxation mechanisms and polydispersity to the shifts of OSR moduli complicate the analysis of estimating changes in chain length in flow. Thus, we cannot simply use conventional scaling analyses developed for LVE measurements to interpret OSR results.

The perturbation analysis and simulation scheme developed in this work are quite general and can be applied for other rheological models and material systems. In fact, we found the result for $yy$-component of the chain conformation determining the plateau value of the orthogonal moduli to be a general result for any model with the upper-convected time derivative of the conformation tensor and relaxation terms that are linear in the conformation tensor. Our result can be used as an analogy to Laun’s rule for relating OSR moduli to the first normal stress difference.

In summary, we found that orthogonal superposition gives very useful information about nonlinear material moduli under flow, which can provide better sensitivity for testing constitutive models for nonlinear polymer processing. Results in our work have important implications for the design and interpretation of future orthogonal superposition experiments.

## ACKNOWLEDGMENTS

This research was primarily supported by U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under No. DE-SC0020988. The authors made use of computational facilities purchased with funds from the National Science Foundation (No. CNS-1725797) and administered by the Center for Scientific Computing (CSC). The CSC is supported by the California NanoSystems Institute and the Materials Research Science and Engineering Center (MRSEC; NSF No. DMR 1720256) at UC Santa Barbara. The authors wish to thank Dr. Joseph Peterson, Dr. Patrick Corona, and Michael Burroughs for helpful discussions.

## REFERENCES

*Rheology: Principles, Measurements, and Applications*, Advances in Interfacial Engineering Series (VCH, New York, 1994).