Tethered Semiflexible Polymer under Large Amplitude Oscillatory Shear

The properties of a semiflexible polymer with fixed ends exposed to oscillatory shear flow are investigated by simulations. The two-dimensionally confined polymer is modeled as a linear bead-spring chain, and the interaction with the fluid is described by the Brownian multiparticle collision dynamics approach. For small shear rates, the tethering of the ends leads to a more-or-less linear oscillatory response. However, at high shear rates, we found a strongly nonlinear reaction, with a polymer (partially) wrapped around the fixation points. This leads to an overall shrinkage of the polymer. Dynamically, the location probability of the polymer center-of-mass position is largest on a spatial curve resembling a limaçon, although with an inhomogeneous distribution. We found shear-induced modifications of the normal-mode correlation functions, with a frequency doubling at high shear rates. Interestingly, an even-odd asymmetry for the Cartesian components of the correlation functions appears, with rather similar spectra for odd x- and even y-modes and vice versa. Overall, our simulations yielded an intriguing nonlinear behavior of tethered semiflexible polymers under oscillatory shear flow.


Introduction
The structural and rheological properties of polymers are strongly affected by external fields, such as shear or extensional flows. As such, this is well established, and the external fields provide means to control the behavior of polymer solutions, melts, and networks [1][2][3][4][5]. Similarly, it has been shown that the macroscopic rheological behavior of a polymer solution, e.g., shear-rate dependent viscosities, normal-stress differences, and shear thinning, are tightly linked to the microscopic polymer conformational and dynamical properties [5][6][7][8][9][10][11]. Hence, insight into the behavior of individual polymers is fundamental in the strive to unravel the macroscopic nonequilibrium polymer properties.
Direct observation of the nonequilibrium properties of single molecules, also termed "molecular rheology" [5], in experiments [12][13][14][15][16][17][18][19] and simulations [10,13, has provided valuable visual illustrations of polymer conformations and has helped to characterize their nonequilibrium properties in terms of their deformation, orientation, relaxation dynamics, and rheology, both for free and tethered molecules in shear and extensional flow. These studies, however, typically probe the linear viscoelastic properties, which are usually insufficient to fully characterize the nonlinear aspects of polymers under flow.
In order to probe nonlinear properties of complex fluids, the large-amplitude oscillatory shear (LAOS) method was developed [19,[45][46][47][48]. Here, in contrast to small-amplitude oscillatory shear, the stress response is typically no longer sinusoidal, but of rather complex shape. Such strong, time-dependent flows will affect the conformational properties of individual polymers to a yet unresolved extent. Studies on the dynamics of single polymers under large-amplitude oscillatory extensional flow yield qualitatively different stretching flow-rate curves (Lissajous curves) as a function of the extension rate and the oscillation frequency [19], and illustrate the complex interplay between time-dependent flows and polymer conformations. Here, further studies are desirable to resolve the time-dependent conformational properties of single polymers under large-amplitude oscillatory flows, aspects which, so far, have not been addressed by simulations.
In this article, we perform non-hydrodynamic, Brownian-type simulations of individual polymers exposed to large-amplitude oscillatory shear. The ends of the polymer are fixed and the polymer dynamics is constraint to the xy-plane of a Cartesian reference frame. We consider stiffer polymers only, with the persistence lengths L p /L = 1/2 and 2, respectively (L p is the persistence and L the polymer contour length). In any case, the polymers exhibit cyclic conformational modulations, specifically at higher Weissenberg numbers (Wi), which is the product of the applied shear rate and the longest polymer relaxation time. The tethering of the ends leads to a more-or-less linear oscillatory response at small Weissenberg numbers, where a polymer moves back and forth like grass swaying in the wind. With increasing Wi, polymers (partially) wrap around the fixation points and shrink. In general, the probability of the center-of-mass position is largest on a limaçon-type curve, as a consequence of the periodic excitation, however, with a non-uniform probability. Interestingly, the center-of-mass autocorrelation function normal to the line connecting the tethering points exhibits frequency doubling with respect to the imposed shear frequency, as a consequence of the non-crossability of the polymer and the tethering points. This reflects the symmetry breaking of the polymer dynamics during an oscillation cycle. Our findings illustrated the complex nonlinear interplay of polymer internal degrees of freedom and external periodic oscillations.
The rest of the paper is organized as follows. In Section 2, the polymer model and the coupling to the shear flow are described. Section 3 presents results for conformational and dynamical properties. Our findings are summarized in Section 4.

Model and Method
The two-dimensional linear polymer chain is composed of N beads of mass M, with its ends, r 1 , r N , tethered at r 1 = (−H/2, 0) T and r N = (H/2, 0) T , respectively (cf. Figure 1). The contour length L = (N − 1)r 0 , where r 0 is the bond length, is fixed and L > H. The interactions between the beads are defined in terms of the potential U = U bond + U bend + U ex , comprising bond, bending, and excluded-volume interactions. The bonds between consecutive beads are described by the harmonic potential: where r i is the position of bead i (i = 1, . . . , N) and κ h is the elastic constant. Bending restrictions are captured by the potential: with κ the bending rigidity and ϕ i the angle between two consecutive bond vectors. In the limit κ → ∞, the persistence length is given by L p = 2κr 0 /k B T, with T the temperature and k B Boltzmann's constant. Bead overlapping and bond crossings are prevented by the shifted and truncated Lennard-Jones potential: where r is the distance between two non-bonded beads, and Θ(r) is the Heaviside function; Θ(r) = 0 for r < 0 and Θ(r) = 1 for r > 0. The dynamics of the beads is described by Newton's equations of motion, which are integrated by the velocity-Verlet algorithm with time step ∆t p [49,50].
The polymer is coupled to a Brownian heat bath implemented via the Brownian (or random) multiparticle collision dynamics (B-MPC) approach [51][52][53]. Hence, no hydrodynamic effects are considered in the present work. MPC consist of streaming and collision steps, where collisions occur in regular time intervals of length ∆t [52,54]. During streaming, the dynamics of the beads is described by Newton's equations of motion. In the collision step, the velocities of the beads change in a stochastic manner. In B-MPC, the Brownian interaction of a bead with the surrounding fluid is implemented by a stochastic collision with a phantom particle, taking its momentum from the Maxwell-Boltzmann distribution of variance Mk B T and a mean, which is zero in absence of shear. In the presence of oscillatory shear in the xy-plan the mean momentum is Mv s i (t), with the shear velocity: at the time t;γ is the shear rate and ω the frequency (cf. Figure 1). The collision is implemented via the stochastic rotation dynamics variant of MPC [52,55,56]. Here, the relative velocity of a bead, with respect to the mean of the velocities of the bead and related phantom particle, is randomly rotated by angles ±α. We choose the following parameters for the simulations: α = 130 • , M = 5, ∆t = 0.1t u , with the time unit t u = Mr 2 0 /(5k B T), κ h r 2 0 /(k B T) = 4 × 10 3 , /(k B T) = 1, σ = r 0 , N = 101, hence, the polymer length is L = 100r 0 , and ∆t p = 10 −2 ∆t. The value of κ h guarantees that the polymer length is constant within 1% for all considered systems. Two bending rigidities are considered, corresponding to the persistence lengths L p /L = 0.5 and 2.
Simulations of free polymers yield the longest relaxations times τ r /t u = 1.9 × 10 6 and τ r /t u = 3.9 × 10 6 for the two stiffness values, determined from the end-to-end vector correlation function [11,40]. The strength of the shear flow is characterized by the Weissenberg number Wi = τ rγ , for which the values Wi = 10, 25, 50, and 100 are considered. The frequency ω is related to the Deborah number De = τ r ω, where we set De = 10.
The polymer, with H = L/5, is initialized with beads along a semi-elliptical contour with the major axis along the y-direction and minor axis along the x-axis. The polymers are equilibrated up to 5 × 10 6 t u > τ r , and data are collected up to the longest simulated time t L = 10 8 t u . Figure 1. Sketch of the tethered bead-spring polymer exposed to oscillatory linear shear flow.

Results
The oscillatory flow dragged the beads along and, at least at small shear rates, the polymer moved back and forth like grass swaying in the wind. This is illustrated in Figure 2. The stiffer polymer (L p /L = 2) closely maintained its shape at low shear rates (Figure 2b Figure 3 shows the time-dependence of the x-coordinate of the center-of-mass position for various Weissenberg numbers. For smaller Wi 25, entropic effects were strong, and x cm was only partially following the external flow. Perturbations were stronger when the polymer got trapped by the fixation points and some time was needed to disentangle it (cf. Supplementary Movies for L p /L = 0.5, 2 and Wi = 25). There were in-phase periods with a small phase shift, which were interrupted by time intervals with a center-of-mass motion decoupled from the flow. A stronger flow (Wi ≈ 50) enhanced the in-phase periodic motion, but the center-of-mass dynamics was phase shifted and seemed to be no longer harmonic. The modulations of the approximately periodic peaks became more pronounced for Wi = 100, and an original single peak split into two peaks, with the minimum of the second peak close to zero. Hence, x cm exhibited an approximate doubling of the frequency, an aspect which is more closely discussed in Section 3.2.2 in the context of normal modes. Overall, we found a highly nonlinear response of the polymer to the external excitation. This modified the polymer conformational and dynamical properties (cf. Supplementary Movies for L p /L = 0.5, 2 and Wi = 100).  Figure 4 depicts the probability distribution of the center-of-mass position for various Weissenberg numbers. The probability was high for positive y cm -values, specifically at lower Wi. This was related to the chosen initial conditions, with a polymer always in the half-plane y > 0. Evidently, the polymers were too stiff to restore isotropy normal to the shear-flow direction. The anisotropy was maintained at higher Wi. However, the most probable center-of-mass position shifted gradually to larger |x cm |. At large Weissenberg numbers, the probability increased in the vicinity of x cm = 0. As illustrated in Figure 2g,h, the polymers were wrapped around the tethering points by the flow. Overall, the probability of the center-of-mass position was highest on a limaçon-type curve, however, with a non-uniform distribution. A limaçon is defined as a curve formed by the path of a point fixed to a circle, when that circle rolls around the outside of a circle of equal radius. In our case, the shear flow (partially) rotated (oscillated) the semiflexible polymer of more-or-less circular shape, which looked like rolling, and the center-of-mass followed a limaçon.

Center-of-Mass Properties
The polar coordinate representation of a limaçon is r = a + b cos θ, and the parameter representation for our reference frame is [57]: with the off-set y 0 . The fits of Equations (5) and (6) to our simulation results for L p /L = 2 are displayed in Figure 4. For Wi = 10, the limaçon was hardly distinguishable from an ellipse. With increasing Weissenberg number, we approached a limaçon, and for Wi = 100, the limaçon turned into a cardioid, where a = b. The limaçon curves for Wi = 50, and especially Wi = 100, did not fully agree with the simulation data for y cm /r 0 10. The reason was the conformational freedom of the polymer and shear-induced shape changes, implying changes in the radii of the rolling circles underlaying the mathematical limaçon construction. These changes were most pronounced while the polymer explored regions of high shear rate.
The distribution functions for the center-of-mass Cartesian coordinates of the polymers are displayed in Figure 5. The distribution function for P(x cm ) (Figure 5a,b) is symmetric with respect to the center between the tethered ends. At small shear, distinct "off-center" peaks were present as a consequence of the projection onto the x-axis. With the increasing Weissenberg number, the peaks became more pronounced and shift closer to the end positions. This is also reflected in the variance x 2 cm , which increased with Wi. At high Weissenberg numbers, a peak appeared in the center, reflecting the polymer "wrapping" around the fixed ends. The asymmetry in the initial condition of y cm is also reflected in Figure 5c,d, with a pronounced peak at positive y cm . The probability of smaller y cm increased with increasing Wi, and for Wi 100, a peak appeared at y cm < 0, consistent with the high probability in the vicinity of x cm = 0 of Figure 4g,h. The snapshots in Figure 2g,h show stretching and alignment along the x-direction of the polymer-part between the fixed end point and the point where wrapping around the other fixed end appears. For such conformations, the major parts of the polymer were at y < 0. The wrapping combined with the reversal of the polymer advective dynamics resulted in a slow lateral dynamics, resulting in a high probability with y cm < 0. Yet, the average y cm was always positive, as shown in Figure 6, but decreased quickly with increasing Wi. A fit with the logarithmic Weissenberg-dependence y cm = −βr 0 ln Wi yielded β = 4.9 and 7.6 for L p /L = 1/2 and 2, respectively. Evidently, the drop is more pronounced for the stiffer polymer.
The stronger conformational changes with increasing strength of the shear flow imply a shrinkage of the polymer. The mean square radius of gyration shrinks by 10-15%, with respect to non-sheared conformation.

Normal Mode Expansion
We studied the internal polymer conformational and dynamical properties via the mode amplitudes A n (t) = (A nx (t), A ny (t)) T of the eigenfunction expansion of a polymer with fixed ends [58], withx the unit vector along the x-axis and the wave numbers q n = nπ/(N − 1) (n = 1, . . . , N − 1). The mode amplitudes are: in terms of the bead positions S ix = r ix − (2i − N − 1)H/2(N − 1) and S iy = r iy (i = 1, . . . , N). The mean, S , and mean square, S 2 , values of the components of S are displayed in Figure 7 for L p /L = 2. The behavior was qualitatively similar for L p /L = 1/2. The shape of S x was centrosymmetric with respect the polymer center i = 50. Any deviation was a consequence of statistical inaccuracy. The magnitude | S x | decreased significantly with the increasing Weissenberg number, and the mean was close to zero for Wi = 100. The component S y was always positive and the largest amplitude was naturally in the polymer center. With increasing Wi, the amplitudes shrunk, and S y was close to zero for all beads. The mean square values S 2 y (Figure 7, right) decreased with increasing Wi. However, S 2 x changed in a nonmonotonic manner, and the fluctuations were larger for Wi = 50 than for Wi = 100. This was a consequence of the symmetry-breaking shear flow with respect to the x-axis and the wrapping of the polymers around the tethering points.
on the mode number for the two stiffness values and the various Weissenberg numbers. For small Wi, we obtained a 1/n 4 dependence as characteristic for stiff polymers [59,60]. This even applied over a range of mode numbers 4 < n 40 at higher shear. However, the variances of the lower-mode amplitudes deviated from this dependence, with a weaker mode-number dependence for small n, and a rapid drop from mode n = 3 to n = 4. This reflected the large-scale conformational changes by the wrapping of the polymers.

Dynamical Properties
Since the large-scale properties of the polymer were modified most by the shear flow, we more closely considered the dynamics of the mode amplitudes for the modes n = 1 and 2 by the mode-autocorrelation function C n (t) = A n (t) · A n (0) . Results for various Weissenberg numbers and the persistence length L p /L = 2 are presented in Figure 9. The simulation data were analyzed by fitting the exponentially damped periodic function: where T 0 = 2π/ω is the period of the applied oscillation, γ characterizes the damping, and Ω accounts for variations of ω. The correlation C 1x (t) = A 1x (t)A 1x (0) , in Figure 9a, decreased at short times and increased again for t/T 0 > 1/2. Fitting yielded a factor Ω ≈ 1.2, i.e., C 1x (t) followed roughly the shear flow. At high Weissenberg numbers, an additional weak modulation appeared, which, however, did not change the primary frequency. The decay of the correlation function C 1y (t) = A 1y (t)A 1y (0) showed a strong Weissenberg number dependence. As depicted in Figure 9b, C 1y (t) depended only very weakly on time for Wi = 10. With the increasing Weissenberg number, the decay rate, γ, increased, and the correlation function showed damped oscillations, which were most pronounced for Wi = 100. Interestingly, the characteristic frequency was twice the externally applied frequency (Ω ≈ 1.8). The correlation function C 2x (t) exhibited a very similar time and Weissenberg number dependence as C 1y (t) (cf. Figure 9c), in terms of drop at the various Wi, as well as the characteristic frequency. The correlation C 2y (t) (cf. Figure 9d) was very similar to C 1x (t), where the fitted expression closely followed the simulation data. Here, we found Ω ≈ 1, i.e., the external frequency, but the decay was stronger with γ ≈ 1.4. Moreover, within both sets of correlation functions, the exponential decay depended only weakly on the Weissenberg number. It was easily recognized that there was an odd-even asymmetry between the xand y-correlation functions. This asymmetry was also obtained for higher mode numbers. Theoretical models of flexible and semiflexible polymers predict an exponential decay of the normal-mode correlation functions of the form e −t/τ n , where the τ n are the relaxation times [6,61]. In the presence of external fields, e.g., shear flow, the time dependence is modified and the normal-mode correlation functions do not necessarily decay exponentially anymore [11], but the exponential factor is still determined by the relaxation times. We did find an initial exponential decay, but no clear mode-number dependence. Figure 8 shows a deviation of the smaller mode numbers from the dependence τ n ∼ 1/n 4 , but the relaxation times of the modes n = 1 and 2 were still different. This is not reflected in Figure 9, where the decay of A 1x (t)A 1x (0) and A 2y (t)A 2y (0) differs by approximately a factor of 2.5, which is roughly the ratio of the mode numbers. Hence, the oscillatory flow field strongly affected the relaxation times, at least at higher Weissenberg numbers. This reflects and is consistent with a rather nonlinear response of the polymer with respect to the external periodic excitation.
The autocorrelation functions for the center-of-mass Cartesian coordinates displayed in Figure 10 help to understand the appearance of the frequency doubling. Evidently, both Cartesian components exhibited a periodic motion, with typically larger amplitudes for higher Weissenberg numbers. Clearly, the x-component shows the same frequency as the applied flow, with some modulations appearing for the highest value of Wi. In contrast, the y-component revealed frequency doubling for Wi 25. The Supplementary Movies illustrate the different dynamical features of x cm and y cm . In one period, the x-position of the center of mass, let us say, moved from the maximum positive value to the minimum value (half period) and, by the oscillation, back again to the maximum x-value. The oscillation was similar to a complete "rotation" (period), and no difference between oscillation and rotation was visible. The behavior of the amplitude of y cm was different. While in a cycle, x cm decreased, y cm increased first, reached its maximum for x cm = 0, and decreased then again. Instead of being able to complete a full cycle, the oscillation and the tethering constrained y cm to move back along a similar path as in the first half of the period, i.e., y cm (t + T 0 /2) increased, reached again a maximum and decreased then again to a value close to the initial value. Hence, y cm exhibited two maxima during a period, whereas x cm exhibited only one. As a consequence, the correlation function y cm (t)y cm (0) showed twice the frequency than the respective correlation function of the x-coordinate. The primary reason was the self-avoidance by the tethering points, which forced a reversal of the dynamics along the y-direction.

Summary and Conclusions
We have analyzed the nonequilibrium properties of semiflexible polymers confined in two dimensions with tethered ends exposed to oscillatory shear. The applied Brownian multiparticle collision dynamics algorithm neglected hydrodynamic interactions. It coupled the polymer dynamics in a stochastic (Brownian) manner to the local flow field and allowed for fluctuations, since the average flow field was imposed only.
For small shear rates, low Weissenberg numbers (Wi 10), the tethering of the ends led to a more-or-less linear oscillatory response, where the polymers moved back and forth like grass swaying in the wind. With increasing Wi, the polymers (partially) wrapped around the fixation points and a more complex, nonlinear response emerged. In fact, the wrapping significantly changed the polymer conformations, and the overall size, measured by the radius of gyration, shrunk. Dynamically, the probability of the polymer center-of-mass position was largest on a spatial curve resembling a limaçon, although the distribution was inhomogeneous. Of course, this was a consequence of tethering and periodic excitation of the polymer, moving it cyclically back and forth, which led to an approximate rolling motion of the center-of-mass, the origin of the limaçon.
At high Weissenberg numbers, we found shear-induced modifications of the mode-spectrum of the mode-amplitude correlation functions. Since we considered semiflexible polymers only, the mode spectrum exhibited the characteristic 1/n 4 dependence for low Wi. At high Wi, the lower modes, specifically the modes n 3, showed a weaker n dependence. This is reflected in a particular dynamical behavior, where a frequency doubling appeared for the normal-mode amplitudes along the xand y-direction. However, the doubling appeared alternating for even and odd modes-the correlation function A 1y (t)A 1y (0) shows twice the frequency of A 1x (t)A 1x (0) , whereas A 2x (t)A 2x (0) is similar to A 1y (t)A 1y (0) . The frequency doubling is also reflected in the autocorrelation function of the y-coordinate of the center-of-mass position. It was a consequence of the hindered and truncated "rotational" motion of the polymer by the tethering points.
In the current simulations, hydrodynamic interactions (HI) have been neglected. Such interactions are of major importance for flexible polymers, since they significantly change the relaxation time spectrum, but are less relevant for semiflexible polymers, where they provide approximately logarithmic corrections only [62]. Hence, we expected only minor differences between the presented results for free-draining polymers and simulations, accounting for HI as long as the polymer motion follows instantaneously the fluid flow.
Our studies revealed an intriguing behavior of tethered polymers under oscillatory flow, which affects their macroscopic rheological behavior. Specifically, the overall shrinkage of the polymer reduced the viscosity. The consequences of the frequency doubling on macroscopic properties need to be further analyzed.