Dynamical structure of entangled polymers simulated under shear flow

TL;DR

This paper uses computer simulations to analyze how entangled polymers' microscopic structures respond dynamically to shear flow, revealing directional differences in relaxation and diffusion behaviors.

Contribution

It introduces an extended blob model simulation to study the dynamical structure of entangled polymers under shear flow, providing new insights beyond rheological descriptions.

Findings

Relaxation speeds up along the flow direction.

Relaxation slows down along the shear gradient.

Complex diffusion behavior in vorticity direction.

Abstract

The non-linear response of entangled polymers to shear flow is complicated. Its current understanding is framed mainly as a rheological description in terms of the complex viscosity. However, the full picture requires an assessment of the dynamical structure of individual polymer chains which give rise to the macroscopic observables. Here we shed new light on this problem, using a computer simulation based on a blob model, extended to describe shear flow in polymer melts and semi-dilute solutions. We examine the diffusion and the intermediate scattering spectra during a steady shear flow. The relaxation dynamics are found to speed up along the flow direction, but slow down along the shear gradient direction. The third axis, vorticity, shows a slowdown at the short scale of a tube, but reaches a net speedup at the large scale of the chain radius of gyration.

Figures

Equations

$$\zeta \frac{\partial {\mathbf{R}}_c(t,s)}{\partial t}$$

$$+\frac{Nv}{{\lambda }^3}\sum _{c^{\prime }=1}^C{\int }_0^{1}d{s}^{\prime }F[{\mathbf{R}}_c(t,s)-{\mathbf{R}}_{c^{\prime }}(t,s^{\prime })]$$

$$+\zeta \dot{\gamma }\left({\mathbf{R}}_c(t,s)\cdot {\mathbf{e}}_y\right){\mathbf{e}}_x$$

$$+\sqrt{2k_{B}T\zeta }\text{lower}0.4pt\text{vrule}height=6.94444pt,depth=-5.55559pt{\mathbf{W}}_c(t,s)$$

$$\tau =\frac{6\pi {\eta }_s{b}^3}{k_{B}T},$$

$$F(\mathbf{r})=\frac{k_{B}T\mathbf{r}}{{\lambda }^2}{e}^{-{\mathbf{r}}^2/2{\lambda }^2}.$$

$$\lambda \propto {\rho }^{-3/4}\text{and}N\propto {\rho }^{5/4}{M}_w$$

$$V=L^3=v_0\left(\frac{4\pi }{3}\right)NC{\lambda }^3$$

$$g_{\alpha \beta }(t)=⟨[\mathbf{R}(t)-\mathbf{R}(0){]}_{\alpha }[\mathbf{R}(t)-\mathbf{R}(0){]}_{\beta }⟩$$

$$\mathbf{R}(t)\to \mathbf{R}(t)-{\mathbf{e}}_x\sum _{t^{\prime }/\Delta t=0}^{t/\Delta t}\left(\dot{\gamma }\Delta t\right){\mathbf{e}}_y\cdot {\mathbf{R}}_{\text{cm}}(t^{\prime })$$

$${\tau }_e=\frac{a^4\zeta }{k_{B}T{b}^2},{\tau }_d=\frac{\zeta N^3{b}^4}{{\pi }^2{k}_{B}T{a}^2}$$

$$F(\mathbf{q},t)=\frac{1}{J}⟨\sum _{j,j^{\prime }=1}^J\exp [i\mathbf{q}\cdot ({\mathbf{R}}_j(t)-{\mathbf{R}}_{j^{\prime }}(0))]⟩$$

$$E(t)=F_{\text{Wi}}(\mathbf{q},t)/F_{\text{Wi}}(\mathbf{q},0)-F_0(\mathbf{q},t)/F_0(\mathbf{q},0)$$

$$\mathbf{R}\to \mathbf{R}+L[(\gamma n_y+n_x){\mathbf{e}}_x+n_y{\mathbf{e}}_y+n_z{\mathbf{e}}_z],$$

$$\gamma (t)={\int }_{-\infty }^t\dot{\gamma }(t^{\prime })d{t}^{\prime }$$

$$\gamma \to \gamma -\text{round}(\gamma )$$

$$\mathbf{r}=x{\mathbf{e}}_x+y{\mathbf{e}}_y+z{\mathbf{e}}_z.$$

$${\hat{\mathbf{E}}}_1$$

$${\hat{\mathbf{E}}}_2$$

$${\hat{\mathbf{E}}}_3$$

$$g_{ij}={\hat{\mathbf{E}}}_i\cdot {\hat{\mathbf{E}}}_j=\left(\begin{array}{ccc}1& \gamma & 0\\ \gamma & 1+{\gamma }^2& 0\\ 0& 0& 1\end{array}\right)$$

$$\hat{\mathbf{r}}={\hat{x}}^1{\hat{\mathbf{E}}}_1+{\hat{x}}^2{\hat{\mathbf{E}}}_2+{\hat{x}}^3{\hat{\mathbf{E}}}_3$$

$${\hat{x}}^1$$

$${\hat{x}}^2$$

$${\hat{x}}^3$$

$$x$$

$$y$$

$$z$$

$$U(\mathbf{r})=\exp \left[-\frac{{\mathbf{r}}^2}{2{\lambda }^2}\right]=\exp \left[-\frac{\left(x^2+y^2+z^2\right)}{2{\lambda }^2}\right].$$

$$\hat{U}(\hat{\mathbf{r}})=\exp \left[-\frac{{\left({\hat{x}}^1+\gamma {\hat{x}}^2\right)}^2+{\left({\hat{x}}^2\right)}^2+{\left({\hat{x}}^3\right)}^2}{2{\lambda }^2}\right]$$

$$\hat{\mathbf{f}}(\hat{\mathbf{r}})=-\hat{\nabla }\hat{U}(\mathbf{r})=-\left[{\hat{\mathbf{E}}}^1\frac{\partial \hat{U}}{\partial {\hat{x}}^1}+{\hat{\mathbf{E}}}^2\frac{\partial \hat{U}}{\partial {\hat{x}}^2}+{\hat{\mathbf{E}}}^3\frac{\partial \hat{U}}{\partial {\hat{x}}^3}\right]$$

$${\hat{\mathbf{E}}}^1$$

$${\hat{\mathbf{E}}}^2$$

$${\hat{\mathbf{E}}}^3$$

$$g^{ij}={\hat{\mathbf{E}}}^i\cdot {\hat{\mathbf{E}}}^j=\left(\begin{array}{ccc}1+{\gamma }^2& -\gamma & 0\\ -\gamma & 1& 0\\ 0& 0& 1\end{array}\right)$$

$${\hat{f}}_1$$

Full text

††preprint: A18.04.0117

I Introduction

Entangled polymers are a network of long, interpenetrating chains, as illustrated in Fig. 1. The motion is driven by random thermal fluctuations, while the excluded volume repulsion prevents the chains from crossing each other. The resulting dynamics display a rich variety of phenomena, distinct from most other materials. Despite its macroscopically complicated behaviour, the physics of the polymer network is mostly governed by one parameter: the number of entanglements $Z=N/N_{e}$, which is proportional to the chain length $N$. The entanglement length $N_{e}$ determines the onset of multi-chain effects, and its value may depend on the chemical composition Cho, Ahn, and Lee (2004), temperature Cavaille et al. (1987), and concentration Kasaai, Charlet, and Arul (2000); Wingstrand et al. (2015). However, these parameters typically do not influence the macroscopic properties, but only shift the time and the length scales. This universality suggests that most of polymer physics might be described in terms of one fundamental model, perhaps containing only one parameter, $Z$.

Theoretical descriptions usually simplify the full multi-chain network with just one chain immersed in the averaged confining field of all the other chains. Currently the main reference is the GLaMM theory Graham et al. (2003), rooted in the earlier Doi and Edwards theory Doi and Edwards (1988), itself based on de Gennes work de Gennes (1971). It is widely applied in rheology Han (2007) where a strong shear flow results in a drop of viscosity (shear thinning), and the emergence of normal stress differences (the Weissenberg effect). The theory also corroborates structural measurements, in particular small angle neutron scattering (SANS) and neutron reflectometry (NR). These experiments have established that the polymer shape elongates in the shear flow direction $x$, while shrinking perpendicular to it Chennevière et al. (2016); Korolkovas et al. (2017); Yearley et al. (2010); Bent et al. (2003); Graham et al. (2006); Clarke et al. (2010); Muller, Pesce, and Picot (1993); Blanchard et al. (2005); Boué (1987); Kirkensgaard et al. (2016), along both shear gradient $y$ and vorticity $z$. However, there is some disagreement with the theory on the exact relaxation pathways as evidenced by a recent kinetic SANS experiment Wang et al. (2017).

Rheology and SANS provide a fundamentally static molecular picture: the stress tensor and the average distance between labeled monomers. The time evolution of these quantities can also be measured in kinetic experiments Wang et al. (2017). The next level of understanding is the dynamical structure, which tells us how fast the various segments of the polymer are fluctuating and how correlated are they. While numerous dynamical experiments have been published Gooßen et al. (2015); Schleger et al. (1998); Wischnewski et al. (2002); Richter et al. (1993a, b); Pyckhout-Hintzen, Wischnewski, and Richter (2016), they are mostly limited to equilibrium fluctuations. By contrast, practical and industrial use of entangled polymers usually involves flow, with examples in food thickeners Rao (2014), paint formulations Chhabra and Richardson (2011), plastics injection Bartczak (2005) and extrusion Bartczak (2010), fibre spinning Shenoy et al. (2005), and enhanced oil recovery from shale rock Wever, Picchioni, and Broekhuis (2011). Non-equilibrium entanglement dynamics also play a biological role in areas like the cytoskeleton Tsang et al. (2017) during cell division, and the cartilage June et al. (2009) during joint movement.

The most complicated, but potentially very rewarding experiments aim to measure polymer dynamics while in the non-equilibrium state. Such information quantifies the molecular motion directly, rather than having to infer it from the change of static structure, as in kinetic SANS. However, due to practical limitations, dynamical experiments are rare, although not impossible. One example is video microscopy of entangled DNA, which can record the shape fluctuations of an individual stained chain as it undergoes flow Teixeira et al. (2007). A very broad range of chain conformations were detected, raising doubts about the validity of the mean-field assumption used in tube theories. Another option is superposition rheology (SR), where a small oscillatory force is applied on top of a steady shear flow, probing the mechanical relaxation as it is altered by shear. Using solutions of entangled polyisobutene (PIB), the authors have observed a speedup of the relaxation spectrum, and explained it with a macroscopic Wagner model Vermant et al. (1998). In the case of entangled polyisoprene (PI) solutions, a slight speedup along the shear gradient axis was reported using dielectric spectroscopy (DS) Watanabe, Matsumiya, and Inoue (2003), which probes the motion of the whole molecule. The measured effect is much smaller than the tube theory prediction. By contrast, a clear slowdown has been observed with nuclear magnetic resonance (NMR) spectroscopy Böhme and Scheler (2011) probing internal motions of an entangled polypropylene (PP) solution. Recently, neutron spin echo (NSE) was attempted under shear flow, reporting no change along the vorticity axis within the accessible shear rates using an entangled polydimethylsiloxane (PDMS) melt Kawecki et al. (2017).

Perhaps due to the scarcity of reliable experiments, theoretical work has largely neglected the question of non-equilibrium dynamics. Multi-chain computer simulations is a useful addition, since they do not involve any mean field assumptions and can verify, refute, or complement the effective single chain theories. However, just like in experiments, most published simulations have focused on the average steady-state Aoyagi and Doi (2000); Kröger and Hess (2000); Padding and Briels (2003); Baig, Mavrantzas, and Kröger (2010); Dorgan, Rorrer, and Maupin (2012); Halverson et al. (2012); Xu, Chen, and An (2014) and occasionally the transient response Kröger, Luap, and Muller (1997); Cao and Likhtman (2015); Xu et al. (2018). More recent simulations have started to probe the non-equilibrium relaxation dynamics, for both short Kim et al. (2011); Huang, Gompper, and Winkler (2012) and long entangled chains Nafar Sefiddashti, Edwards, and Khomami (2015); Mohagheghi and Khomami (2016); Sefiddashti, Edwards, and Khomami (2017). Some disagreements with tube theories have been pointed out, questioning their validity under strong shear, especially as far as dynamics is concerned.

Clearly, more experiments are needed to resolve the controversies. In this paper we provide the multi-chain simulation data for comparison with anticipated future experiments. We show the specific 3D pathway of how the diffusion morphs from isotropic at short time scales, to anisotropic at the longest time scales. An intriguing behaviour is found along the vorticity axis, where the change of dynamics is non-monotonic: it slows down on the short scale, but speeds up on the long scale, helping to reconcile experimental evidence. Our results will encourage new experiments and contribute to further developments in tube theory.

II Simulation method

The simulation is based on a previously published algorithm Korolkovas, Gutfreund, and Barrat (2016) for equilibrium systems, but now extended for bulk shear flow by imposing a Couette velocity profile, described in detail in Appendix A. The code in MATLAB and CUDA is available for download from a Zenodo repository Korolkovas (2017). We numerically solve the motion of $c=1,\,2,\,\ldots,\,C$ chains in a box, each described by a continuous curve $\mathbf{R}_{c}(t,s)$, with variables $t$ for time and $0<s<1$ for monomer index. There are $N$ degrees of freedom per chain, governed by the first order equation of motion:

[TABLE]

It describes the motion of a continuous Rouse chain (Eq. (4.9) in Ref. Doi and Edwards (1988)), with an added excluded volume field that we have implemented in Ref. Korolkovas, Gutfreund, and Barrat (2016), as well as a shear field that is described in Appendix A of this article. Here $\zeta=6\pi\eta_{s}bN$ is the friction coefficient of the chain center of mass, defining the unit of time:

[TABLE]

where $b$ is the spring length. The thermal noise is modeled by a Wiener process $\braket{\mathbf{W}_{c}(t,s)\mathbf{W}_{c^{\prime}}(t^{\prime},s^{\prime})}=\mathbb{I}\,\delta_{cc^{\prime}}\,\delta(t-t^{\prime})\,\delta(s-s^{\prime})$. The repulsive excluded volume force between the chains is described by the blob model, which is well approximated Bolhuis et al. (2001) by a Gaussian function:

[TABLE]

For a semi-dilute polymer solution de Gennes (1979), the blob radius $\lambda$ and the number of blobs per chain $N$ scale as

[TABLE]

where $\rho$ is the density and $M_{w}$ is the molecular weight. The solvent is implicitly contained in the blob particle and its effect is taken into account via the viscosity $\eta_{s}$ (Eq. (2)). This model can also be extended for polymer melts, $\rho\rightarrow\rho_{0}$, although in that case the mapping of Eq. (4) does not apply and the absolute values of $\lambda$ and $N$ depend on the chemical species and temperature.

Entire studies have been devoted to relate the coarse interaction parameters with the atomistic details of a particular chemical species Sun and Faller (2005). However, within the numerical pre-factor, the physical properties on the large scale are found to be universal across vastly different polymer models, ranging from lattice-based Dorgan, Rorrer, and Maupin (2012), to hard bead-and-spring Kremer and Grest (1990), to our soft blobs Korolkovas, Gutfreund, and Barrat (2016). We do not attempt a one-to-one correspondence with any particular experiment, and for a generic model we pick the natural choice of the excluded volume $v=\lambda^{3}=b^{3}$. With these settings, the excluded volume and the spring forces have the same strength, enabling a maximum time step. Other simulations often use a bending potential Cao and Likhtman (2012), which has the advantage of decreasing the entanglement length $N_{e}$. In our case, we could obtain a similar effect by increasing the chain stiffness $\lambda/b>1$, but that would also require increasing $v$ to prevent the contour from shrinking, and this reduces the time step since $\Delta t\propto 1/v$, thus offsetting the benefit of a smaller $N_{e}$.

While the chain crossings are not strictly forbidden, their probability can be suppressed sufficiently so that highly entangled dynamics emerge, as explained in our equilibrium publication Korolkovas, Gutfreund, and Barrat (2016). First, the continuous equation of motion (1) is sampled with $J$ discrete points, and choosing $J\gg N$ guarantees that the gaps between those points remain much smaller than the blob radius. In other words, there are more interaction centers ($J$) than degrees of freedom ($N$), enabling the time step $\Delta t$ to be on the large scale of the $\lambda$-sized blobs $N$, rather than on the much shorter scale of the $J$-points. A special case of $J=N$ would correspond to standard bead and spring simulations. Second, the instantaneous random force must be much weaker than the excluded volume force. Note that we are working with a first order equation of motion, which is overdamped and does not conserve momentum, unlike the second order equation of motion. In this case, we can decrease the random force magnitude by 1/\mathchoice{{\hbox{\displaystyle\sqrt{M,}}\lower 0.4pt\hbox{\vrule height=6.83331pt,depth=-5.46667pt}}}{{\hbox{\textstyle\sqrt{M,}}\lower 0.4pt\hbox{\vrule height=6.83331pt,depth=-5.46667pt}}}{{\hbox{\scriptstyle\sqrt{M,}}\lower 0.4pt\hbox{\vrule height=4.78333pt,depth=-3.82668pt}}}{{\hbox{\scriptscriptstyle\sqrt{M,}}\lower 0.4pt\hbox{\vrule height=3.41666pt,depth=-2.73334pt}}}, provided that we only update its direction once per $M$ steps, effectively cutting off high frequency fluctuations, while maintaining the mean square displacement fixed in the long run. If the cutoff time $M\Delta t$ remains much smaller than the entanglement time $\tau_{e}$ (Eq. (8)), the physics on the long time scale are not affected, and a value of $M=120$ was found to be a suitable compromise. A similar argument is made in the modified Kremer-Grest models, where a bending potential is introduced, justified by the fact that it does not influence the physics, as long as the chain persistence length remains much shorter than the entanglement length Bačová et al. (2013).

When shear is applied, the probability of chain crossings obviously increases. However, the increase remains negligible for shear rates below the inverse Rouse time: $\dot{\gamma}\tau N^{2}\lesssim 1$, where the main effect is the reorientation of the entanglement network. We do not apply a higher shear, since the chain backbone begins to stretch at that point. We use chains of length $N=512$, shorter than the longest one in our equilibrium study. The appeal of working with a coarse model is that high frequency fluctuations are not interesting on the large scale Ilg, Öttinger, and Kröger (2009), which automatically precludes an excessively strong shear. Rheological experiments seldom go beyond $\text{Wi}=10$, and there is plenty of physics to be learned even in a modestly non-Newtonian regime. To the contrary, the ultra high shear rates of order ${10}^{8}\text{\,}\mathrm{s}^{-1}$, routinely considered in atomistic polymer simulations Kröger and Hess (2000), are never realized experimentally, and may only occur during a cataclysmic event like an asteroid impact Ramesh (2008). Moreover, strong shear is known to destroy polymers by rupturing covalent bonds Nghe, Tabeling, and Ajdari (2010).

In general, computer simulations may suffer from a number of technical biases: the chains are too stiff, the box is too small, the system is not properly equilibrated, and there may be chain crossings. These biases all result in faster apparent dynamics, while trying to minimize them makes the program run slower, limiting the accessible length and time scales. We have strived to reach a reasonable compromise between these limitations. Dedicated studies are necessary to fully assess chain crossings Padding and Briels (2001), going beyond the scope of the present article. In the meanwhile, we caution that our results may indeed have a slight bias, more so at the highest shear rates and the longest time scales.

III Results

The simulation was used to explore the changes in relaxation dynamics of polymers when subject to a steady shear flow. We have simulated $C=64$ chains of length $N=512$ in a cubic box of volume

[TABLE]

with a numerical pre-factor $v_{0}=2$ to allow some freedom of movement for the blobs. The microscopic parameters are set to $v=\lambda^{3}=b^{3}$. This melt of blobs is moderately entangled, with about $Z=N/N_{e}=8.5$ entanglements per chain at equilibrium, as determined topologically using the Z1 code Karayiannis and Kröger (2009). Seven simulation runs where performed with all parameters kept fixed except for the shear rate $\text{Wi}=\dot{\gamma}\tau_{d}=0,\,0.3,\,1,\,3,\,10,\,32,\,\text{and}\,100$. The computations were carried out on an Nvidia Quadro P5000 GPU card for about 4 days per run, generating trajectories of $20\tau_{d}$ relaxation times each. To avoid any start-up effect in the subsequent data analysis, we have discarded the first $\tau_{d}$ steps of each trajectory, so the data reported in this article is for a steady-state shear.

At high shear the system response becomes non-linear, and its topological network is disrupted as evidenced by a decreasing $Z$ number. This disentanglement trend was examined with primitive path analysis in Ref. Baig, Mavrantzas, and Kröger (2010). In the present study we aim to characterize the non-linear steady state from a dynamical standpoint. That is, we look at fluctuations around the deformed structure. Among the many quantities that could be considered, we choose diffusion, because it is clearly visible from the simulation perspective, and it can be measured experimentally with neutron spin echo (NSE) Kawecki et al. (2016). The most topologically confined part of the molecule is the central monomer $\mathbf{R}(t,s=0.5)\equiv\mathbf{R}(t)$, and its diffusion under shear shows the strongest deviation from equilibrium. It is quantified by the mean square displacement (MSD), defined as

[TABLE]

plotted in Fig. 2. In the shear geometry there are four distinct components: $\alpha\beta=(xx,\,yy,\,zz,\,xy)$, respectively referring to shear, gradient, vorticity, and the off-diagonal term in the shear-gradient plane. Concerning the diffusion along the $x$-axis, the raw trajectory is overwhelmed by a huge ballistic component $g_{xx}\propto t^{2}$ of the macroscopic shear flow, so we subtract it before the data analysis:

[TABLE]

The big sum is the total displacement imparted on the center of mass (cm) of every chain, over the course of time from [math] to $t$. In other words, we transform the trajectory to the oblique coordinate system where the effect of advection is removed (see Eq. (18)), retaining only the molecular diffusion.

In equilibrium, the MSD of long entangled chains is known to scale as $t^{1/4}$ (reptation), whereas unentangled Rouse chains follow $t^{1/2}$, while a random walk has $t^{1}$. To highlight the entanglement character, it is customary to divide the MSD by the Rouse law $t^{1/2}$, so that the topologically confined motion has a negative slope, as shown in Fig. 2. The negative portion of the curve is demarcated by two peaks, corresponding to the entanglement and the disentanglement times:

[TABLE]

The above formulas are derived for highly entangled polymers at equilibrium. The idea is that the monomer starts to feel the confining topological network only after having diffused a square distance $g_{\alpha\beta}(\tau_{e})=a^{2}$, where $a$ is the tube radius. It then slowly escapes the initial network, and after traveling a square distance equal to the chain size, $g_{\alpha\beta}(\tau_{d})=Nb^{2}$, it is fully relaxed in a new conformation.

While there is no consensus on how to extend Eqs. (8) to the non-equilibrium regime, it is clear from Fig. 2 that the two scales $\tau_{e}$ and $\tau_{d}$ continue to exist and with increasing Wi are gradually moving away from their equilibrium positions. We use these two dynamical peaks to quantify the extent of topological confinement in terms of time $(\tau_{d}/\tau_{e})_{\alpha\beta}$ and MSD $g_{\alpha\beta}(\tau_{d})/g_{\alpha\beta}(\tau_{e})$ ratios, shown in Fig. 3.

The largest effect of shear is along the flow direction $x$, where $g_{xx}(t)$ increases with shear at every timescale $t$, to the point where the entanglement character (the negative slope) disappears completely between $\text{Wi}=3$ and $10$. An opposite, but weaker effect is found perpendicular to the flow. The shear gradient axis, $y$, shows a slowdown of dynamics at all timescales $t$. The ratio $(\tau_{d}/\tau_{e})_{yy}$ (see Fig. 3) reveals that the monomer remains topologically confined for up to five times longer along this axis. On the other hand, if we look at the ratio $g_{yy}(\tau_{d})/g_{yy}(\tau_{e})$, which quantifies the distance traveled to escape the confinement, it appears to be insensitive to shear. Next, let us inspect the cross-correlated displacement $g_{xy}(t)$ shown on the semilog plot in Fig. 2, normalized by $t$ to reveal its characteristic times $\tau_{e}$ and $\tau_{d}$. Curiously, the confinement time span $(\tau_{d}/\tau_{e})_{xy}$ increases similarly to what is seen in the $g_{yy}$ component, while the distance span $g_{xy}(\tau_{d})/g_{xy}(\tau_{e})$ decreases rapidly just like in the $g_{xx}$ component.

The vorticity $z$ is perpendicular to the shear plane and shows the weakest effect. Here the change of dynamics is not monotonic, slowing down on a short time scale, but then speeding up to an extent where the terminal diffusion coefficient becomes slightly higher than in equilibrium. Noteworthy is that both the time $(\tau_{d}/\tau_{e})_{zz}$ and the distance $g_{zz}(\tau_{d})/g_{zz}(\tau_{e})$ ratios remain roughly constant.

Another way to visualize chain dynamics is the coherent intermediate scattering function:

[TABLE]

It can be measured with NSE and quantifies how likely are two monomers initially at a distance $1/q$ apart to stay correlated after a time $t$. A recent simulation Hsu and Kremer (2017) has investigated this function in detail for strongly entangled polymers in equilibrium. Here we focus on the change that may be observed under shear flow, quantified by the scattering excess

[TABLE]

which is the difference between the intermediate scattering function at rest and under shear, shown in Fig. 4. The excess $E(t)$ can range between $(-1,1)$, where positive (negative) values indicate that there is a slowdown (speedup) of the dynamics. To reduce random noise, the raw data has been smoothed with a moving average filter of width $120000\Delta t=0.12\tau_{d}$. We focus on the two axes perpendicular to shear, the $q_{y}$ and $q_{z}$, as they can be accessed by experiments like NSE and dielectric spectroscopy. A significant dynamical change begins at around $\text{Wi}=1$, where both $y$ and $z$ components display a slowdown starting at low $qR\lesssim 1$ and long $t\gtrsim\tau_{d}$, the magnitude of the effect being $50\text{\,}\mathrm{\char 37\relax}$ greater for the $y$-axis. (The small speedup visible at short $t$ is probably a simulation artifact due to numerical anisotropy imposed by the oblique grid, as explained in Appendix A). At a strong shear of $\text{Wi}=32$, the slowdown not only grows in magnitude, but also shifts to shorter time $t$, even four orders of magnitude shorter than the shear rate itself, demonstrating that a macroscopic shear flow can alter microscopic motion at the scale of barely one entanglement. While the $y$ component is slowed down at all times and distances, the situation is different for the $z$-axis. Its slowdown trend is reversed at around the disentanglement time $\tau_{d}$, where the low $q_{z}$-vectors start speeding up with respect to equilibrium.

IV Discussion

The largest structural change occurs in the flow direction $x$. To reveal the dynamics of those new structures, one must subtract the macroscopic shear flow, which is easier in simulation than experiment. Currently the only technique able to measure the hidden dynamics is the parallel superposition rheology. Although its data analysis is far from perfect Unidad and Ianniruberto (2014), experiments on entangled PI solutions Cleyn and Mewis (1981); Vermant et al. (1998) and more recently wormlike micelles Kim et al. (2013), both report a strong speedup of stress relaxation modulus, in qualitative agreement with our simulation results.

The amplitude of structural change is smaller along the shear gradient $y$ and even smaller along the vorticity $z$. Nevertheless, those changes may be easier to detect experimentally, since they are not overwhelmed by the bulk flow. Dielectric spectroscopy can access the autocorrelation of the end-to-end vector $\Phi(t)=\braket{\mathbf{R}_{\text{ete}}(t)\cdot\mathbf{R}_{\text{ete}}(0)}$ of the polymer chain. Data on entangled PI solutions in the $y$-axis shows a slight speedup of relaxation dynamics, although the change is much smaller than predicted by tube theory Watanabe, Matsumiya, and Inoue (2003). Even though we have not analyzed the simulated end-to-end vector separately, the experimental trend of overall speedup disagrees with our simulation, since we found a clear slowdown for all $q_{y}$ vectors.

Another tool to probe polymer dynamics under shear is NMR spectroscopy. It can access intramolecular dynamics, corresponding to an entanglement strand, or high $q$ region. A slowdown of dynamics in the shear gradient direction $y$ was reported Böhme and Scheler (2011), this time in accord with our simulation data.

Recently, neutron spin echo experiments under in situ shear became possible Kawecki et al. (2016), measuring the intermediate scattering function $F(q_{z},t)$ along the $z$ (vorticity) axis. A PDMS melt sheared at $\text{Wi}=0.4$ showed no change in scattering Kawecki et al. (2017). Although higher shear rates could not be reached at this stage, the null result agrees with our simulation, which shows only minuscule change below $\text{Wi}=1$.

Labelled polymers can be filmed under a microscope to visualize their motion in real space. Observations of flowing entangled DNA Teixeira et al. (2007) and semiflexible actin filaments Harasim et al. (2013) have been described in terms of tumbling dynamics, which are abrupt fluctuations between strongly elongated and collapsed states. This can be understood with a double potential well, where the chain spends most of the time confined by the flowing network, but can occasionally flip around by $180\text{\,}\mathrm{\SIUnitSymbolDegree}$, where it finds another elongated state of an equal probability. The jump must be rapid since being misaligned with the flow is energetically unfavorable. We believe that this tumbling behaviour is reflected in our simulated scattering excess along the $z$ axis (Fig. 4). The elongated chain state corresponds to the slowdown seen at short $t$ and high $q$ (local confinement), while the occasional rapid jumps are seen as faster dynamics at long $t$ and low $q$ (global speedup). In agreement with our simulation, a clear speedup along the $z$ axis is also found in superposition rheology applied perpendicular to shear flow Cleyn and Mewis (1981); Vermant et al. (1998); Kim et al. (2013), which is sensitive to molecular motion at long time and distance scales.

To conclude, the simulation results are in qualitative agreement with most of the dynamics experiments conducted under shear so far. We are looking forward to the emerging technique of rheo-NSE spectroscopy to verify the local slowdown for the $z$-axis as predicted by the simulation.

V Acknowledgements

The authors thank Jean-Louis Barrat for his valuable advice. Gunnar K. Palsson, Anton Devishvili, and Franz A. Adlmann have contributed to the development of the simulation software. The Carl-Tryggers stiftelse as well as the Swedish research council are acknowledged for financial support, grant CTS 16:519. The Quadro P5000 used for this research was donated by the NVIDIA Corporation.

VI Appendix A. Oblique coordinate system

To simulate bulk shear flow, we need to impose special boundary conditions, invented by Lees and Edwards, as illustrated in Fig. 5. A further complication is that we need to use the Fourier transform to speed up the excluded volume force calculation, and may also need it to evaluate the polymer structure factor for comparison with SANS. The technique was first developed in Ref. Kobayashi and Yamamoto (2011), and here we adapt their approach to our situation.

In a cubic box of length $L$, every particle $\mathbf{R}$ is replicated at points

[TABLE]

where the $n_{x}$, $n_{y}$, and $n_{z}$ are integers $-\infty,\,\ldots,\,-1,\,0,\,+1\,\ldots,\,+\infty$. The box strain

[TABLE]

is incremented at every time step, and the oblique-periodic images can shift many times over their own length during a long run. For numerical stability reasons, the strain is always folded to

[TABLE]

restricting it to the symmetric range $-0.5<\gamma<0.5$.

A particle position in the rectangular coordinate system $(\mathbf{e}_{x},\mathbf{e}_{y},\mathbf{e}_{z})$ can be decomposed into its Cartesian components $(x,y,z)$:

[TABLE]

In terms of the rectangular unit vectors, the oblique coordinate system of strain $\gamma$ is defined by the covariant basis:

[TABLE]

In other words, the metric tensor of the oblique system is given by

[TABLE]

As illustrated in Fig. 6, the position of a particle in the oblique system, denoted by $\hat{\mathbf{r}}$, is expressed as

[TABLE]

where its coordinate components in the oblique system can be found by solving the equation $\hat{\mathbf{r}}=\mathbf{r}$:

[TABLE]

The simulated polymer network is periodic in these coordinates. The inverse is

[TABLE]

The interaction potential is a scalar quantity, in our case a Gaussian:

[TABLE]

A scalar remains invariant regardless of a coordinate system, so to express it in the oblique coordinates, we simply plug in the values from Eq. (19):

[TABLE]

The force, in the oblique system, is given by the gradient:

[TABLE]

which is expanded in terms of the contravariant basis vectors:

[TABLE]

These are defined to be conjugate to the covariant basis vectors: $\hat{\mathbf{E}}_{i}\cdot\hat{\mathbf{E}}^{j}=\delta_{i}^{j}$, and in particular $\hat{\mathbf{E}}_{2}\cdot\hat{\mathbf{E}}^{1}=(\gamma\mathbf{e}_{x}+\mathbf{e}_{y})\cdot(\mathbf{e}_{x}-\gamma\mathbf{e}_{y})=0$. Morevover, the contravariant metric tensor

[TABLE]

is clearly the inverse of the covariant one: $g_{ij}g^{ij}=1$. Let us be very clear on what are the specific components of the force, in terms of the oblique contravariant basis:

[TABLE]

Just like the position vector, the force vector can be expanded in either the rectangular or the oblique systems:

[TABLE]

which leads to its rectangular components

[TABLE]

Having presented all the necessary formulas, we now outline the algorithm required to evaluate the excluded volume force in the oblique periodic system.

We start off with a set of all particle positions $(x,y,z)$ given in the usual rectangular coordinates. A fixed strain value $\gamma$ is given as well. 2. 2.

The positions are transformed into the oblique system: $(x,y,z)\rightarrow(x-\gamma y,y,z)$, and the periodic boundary conditions are applied: $\hat{\mathbf{r}}=\hat{\mathbf{r}}-L\,\texttt{floor}\,(\hat{\mathbf{r}}/L)$. We are now folded inside the central oblique box, which is periodic along its covariant unit vectors. 3. 3.

A density histogram $\hat{\rho}(\hat{\mathbf{r}})$ is constructed, which holds the number of particles in every oblique pixel $\hat{\mathbf{r}}=L\left(k_{1}\hat{\mathbf{E}}^{1}+k_{2}\hat{\mathbf{E}}^{2}+k_{3}\hat{\mathbf{E}}^{3}\right)$ where the indices run from $k=1,\,2,\,\ldots,K=\texttt{round}(L/\Delta x)$, and $\Delta x$ is the pixel width. 4. 4.

The density is Fourier-transformed into the oblique reciprocal space:

[TABLE]

Notice that the components of the reciprocal lattice vector $\hat{\mathbf{q}}=(\hat{q}_{1},\hat{q}_{2},\hat{q}_{3})$ come out in the contravariant basis defined in Eq. (23), because this basis has the unique property that

[TABLE] 5. 5.

Steps 3-4 are repeated for the force kernel. First, the force array $\mathbf{\hat{f}}(\mathbf{\hat{r}})$ is declared in the real oblique space, and then it is Fourier-transformed to the reciprocal oblique space to give $\mathbf{\hat{f}}(\mathbf{\hat{q}})$. 6. 6.

The force on every oblique pixel $\mathbf{\hat{r}}$ is found by the convolution theorem:

[TABLE] 7. 7.

Lastly, the oblique force is transformed to the rectangular system (the inverse of step 2, but for the contravariant basis): $(\hat{F}_{1},\hat{F}_{2},\hat{F}_{3})\rightarrow(\hat{F}_{1},\hat{F}_{2}-\gamma\hat{F}_{1},\hat{F}_{3})=(F_{x},F_{y},F_{z})$. These force components are then used in the equation of motion, Eq. (1), solved in the rectangular system as usual.

The strain position is updated at every time step, $\gamma\rightarrow\gamma+\dot{\gamma}\Delta t$, and wrapped periodically: $\gamma\rightarrow\gamma-\texttt{round}(\gamma)$. In addition, it must be noted that the use of the discrete convolution introduces some degree of numerical anisotropy due to the finite size of the pixels. In equilibrium, the pixels are rectangular, and the numerical anisotropy is symmetric around every axis, which does not result in a net bias. When the pixels are oblique, however, the anisotropy tilts and produces a net numerical stress, which is unphysical and is not related to the shear flow. In fact, the numerical stress would persist even at equilibrium, if we would choose to simulate a system in a fixed oblique periodic system of $\gamma\neq 0$. To mitigate this parasitic effect, we keep track of the total accumulated strain

[TABLE]

and for every new timestep we choose the image of $\gamma$ which would bring the accumulated strain $\gamma_{\text{sum}}$ as close to zero as possible. For instance, in the example of Fig. 5, the strain is $\gamma=0.232$, but one could equally well describe the system using a strain $\gamma^{*}=\gamma-1=-0.768$. Therefore, our code propagates two timesteps using the value $\gamma$, and the third step would use the conjugate $\gamma^{*}$, because $|3\gamma-1|<|3\gamma|$.