text:
# Excluded volume effects on the structure of a linear polymer under shear flow

## I Introduction

Excluded volume (EV) and hydrodynamics interactions (HI) have been recognized as fundamental physical phenomena to explain the peculiarities of single chain behavior in dilute polymer solutions. For a solution at rest, their effects on the single chain structure and dynamics have been carefully investigated over the years both theoretically and experimentally. When a dilute solution is subjected to a flow field, the situation is much more complex. In addition to finite extensibility (FE) effects to be considered at high strain, the role of EV and HI effects is still largely debated. This situation certainly arises from the difficulties to solve the kinetic equations for the dynamics of the simple coarse-grained models when EV and HI are considered. Some progresses have been made concerning the HI effects in absence of EV interactions by the use of a self consistent preaveraging approximations , a non equilibrium version of the more famous equilibrium Zimm theory. Renormalization Group Theory techniques have also been applied to solve suitable dynamic equations for the chain . The main focus of those investigations was the rheological response of the system although some results on the chain structure where also provided. An attempt to introduce EV effects in the Rouse model has also been performed by perturbation theory expansion.

On the experimental side the situation is similar in that measurements of single chain properties in flow are considerably more difficult than at rest. However, several scattering experiments have been performed, mainly Light Scattering (LS), with the aim to investigate the global properties of the chain such as the orientation and the deformation of the coil by the flow . The effect of EV interaction over global chain properties has been discussed in a series of recent LS studies  where some systematic trends were reported. For dilute solutions in shear flow at a fixed reduced shear rate, a systematic weakening of the alignment of the main axis of the chain in the flow direction has been observed for increasing solvent quality. At a given reduced shear rate, the measured expansion ratio indicates that good solvent chains are less deformed than theta chains . Inspection of the flow effects on the internal structure of the coil can only be obtained by Small Angle Neutron Scattering (SANS) and indeed some data have been already published. However experimental limitations are here even more severe than for LS, and a quantitative study for various solvent qualities is still missing.

New experimental insight in the single chain behavior under flows and stretching forces has been provided recently by fluorescence techniques on biological molecules like DNA . Two papers on the single chain dynamics under shear flow have been published recently by S. Chu and coworkers  and by D. Wirtz and coworkers  in both cases employing DNA chains. In the former paper,the mean and time fluctuations of the fractional extension of DNA under shear flow in the flow direction was reported for reduced shear rates larger than one and culminating at $`\beta =76`$. In the second paper, the authors have shown that the dynamics of a DNA chain under weak shear flow ($`\beta 0.1`$) presents an unexpected wide range of relaxation times related to the particular initial configuration prior to the flow inset.

An interesting alternative to the experimental tools in the context of polymer physics are the numerical techniques such as Monte Carlo (MC), Molecular Dynamics (MD) and Brownian Dynamics (BD) simulations. They can be used to solve simple models for the system of interest and therefore test or suggest theoretical developments or stimulate new experiments. On the other hand simulations can help interpreting experimental data in a more complete and coherent fashion.

Many numerical studies of the behavior of a single chain under shear flow have already appeared by various techniques such as BD , MD  and even non-equilibrium MC . To mention a recent simulation work , BD studies for a chain without HI and EV found a reasonable agreement with the infinite chain Rouse model predictions but came to the disturbing conclusion that, while the relative elongation of DNA chain in the flow direction is correctly explained by the chain model, the same model fails to explain the global expansion ratio of a polystyrene chain which is experimentally much lower than predicted.

To our knowledge, the role of the excluded volume effect on the structure of the chain under SF has not been systematically studied by simulations. In the present paper we try to elucidate this effect by comparing the analytical results of the Rouse model (in the continuum limit) with BD data for a suitable microscopic model of flexible chain with excluded volume interactions and unavoidable FE features.

The structure of a single chain in shear flow is often discussed by exploiting an analogy with the well known structure of a single chain subjected to a pair of equal and opposite stretching forces applied to its ends. In this case the ”macroscopic” behavior of long chains, as measured by the elasticity law, is characterized by a linear hookean regime for weak forces followed by a non linear power law regime for strong forces and ultimately by specific finite extensibility effects. The crossover between the first two regimes is linked to the appearance of the so called “tensile” blobs in the chain structure.

A very clear signature of the blobs can be detected in the chain structure factor as a crossover from ideal behavior at small q’s to excluded volume statistics at large q’s. This has been shown long time ago by Neutron Scattering experiments for “thermal” blobs in the semi-diluted solutions at rest , and recently by MC simulation for “tensile” blobs in the stretched chain problem. In the latter case, the crossover is detected for q orthogonal to the external force and can be ascribed to the fact that scattering at small q’s in these directions comes from pairs of monomers which are far apart in the elongation direction and therefore do not interact directly through the EV potential.

For the chain in shear flow Onuki has formulated a blob model as a simple extension of the Pincus-deGennes theory for the stretched chain case. The blob size is related to the bare shear rate and the chain at high shear rates is imagined as a string of shear blobs aligned most of the time in the flow direction. Such a structure which corresponds to an aligned fluctuating rod, should yield a scattering signal for $`q`$ in the flow direction (flow and elongation directions provide essentially the same scattering signal at high enough shear rate) similar to the stretched chain signal for $`q`$ aligned in the stretching direction . For the stretched chain such behavior is described quantitatively by the Rouse model with a suitable choice of the chain elongation and the longitudinal and transverse fluctuations . In the present paper we test these ideas on the form factor of a single EV chain in flow for shear rates which are sufficiently high to reach the plausible shear blob regime but also sufficiently small to remain in the scaling regime, i.e. avoiding finite chain effects. Our strategy is to analyze to which extend the low q behavior resembles the single chain structure factor of ideal chains in the presence of shear flow, which can be derived exactly, and whether the high q regime reproduces the structure of unperturbed EV chains. In order to observe such phenomena taking place on different length scales, we need sufficiently long chains and therefore we have performed BD simulations of EV chains in absence and presence of shear flow.

The paper is organized as follows. In section II we set up the essential phenomenology concerning a chain in shear flow and define the relevant structural characteristics of the chain. Section III describes the procedure to solve the Rouse model in flow and provides the results for the specific SF case. This is not a totally new analytical development as we are aware of at least two previous works on it . The novelty here is however the calculation of the structure factor in the whole range of $`q`$ vectors and the use of such behavior to pinpoint the signature of shear blobs in the EV chain structure factor. In the next section IV we describe the model for EV chains and the BD technique used. In section V we discuss the results of BD simulations. Finally section VI is devoted to our conclusion and perspectives.

## II Phenomenological framework

When a dilute polymer solution is subjected to a simple shear flow, characterized by the velocity field $`𝐮=\dot{\gamma }y\widehat{𝐱}`$ where$`\dot{\gamma }`$ is the shear rate, the chains in the solution are oriented and deformed according to the flow. Those phenomena can be measured through the anisotropy arising in any tensorial quantity associated to the chain, such as for instance the gyration tensor $`𝐆`$ and the order parameter tensor $`𝐎`$ defined as

$`𝐆`$ $`=`$ $`{\displaystyle \frac{1}{2N^2}}{\displaystyle \underset{i,j=1}{\overset{N}{}}}<(𝐑_i𝐑_j)(𝐑_i𝐑_j)>`$ (1)
$`𝐎`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N1}{}}}{\displaystyle \frac{𝐮_i𝐮_i}{|𝐮|^2}}{\displaystyle \frac{\mathrm{𝟏}}{3}}`$ (2)
where $`𝐑_i`$ are the coordinates of the i-th monomer in the chain, $`𝐮_i=𝐑_{i+1}𝐑_i`$ and $`\mathrm{𝟏}`$ is the unit tensor. At equilibrium the system is isotropic and any tensorial quantity, say $`𝐀`$ reduces to a scalar quantity. Under shear flow, symmetry imposes $`A_{xz}=A_{zx}=A_{yz}=A_{zy}=0`$. Therefore there are at most four independent components of any tensor to be monitored, namely the three diagonal elements and the off-diagonal $`A_{xy}=A_{yx}`$. The orientational angle $`\chi _A`$, defined through the relation

$$cot(2\chi _A)=\frac{A_{xx}A_{yy}}{2A_{xy}}$$
(3)
measure the rotation around $`\widehat{𝐳}`$ of the principal axes (I,II,III) of the tensor $`𝐀`$ with respect to the flow axes (x,y,z). In shear flow, $`A_{xy}`$ starts linearly with $`\dot{\gamma }`$ while the first contribution to $`A_{xx}A_{yy}`$ is of order $`\dot{\gamma }^2`$. Therefore the linear (Newtonian) regime is characterized by $`\chi _A=\pi /4`$. Outside the linear regime $`\chi _A`$ decreases to zero for increasing shear rate.

The deformation ratios, defined by

$$\delta A_{\alpha \alpha }=\frac{A_{\alpha \alpha }(\beta )}{A_{\alpha \alpha }(0)}1$$
(4)
measures the amount of deformation of the chain in flow, either in the flow reference frame ($`\alpha =x,y,z`$), or in the molecular reference frame ($`\alpha =I,II,III`$). When defined in the flow reference frame, the first contribution to those quantities is of order $`\dot{\gamma }^2`$, so that the Newtonian regime is characterized by the absence of deformation. Outside the linear regime, the chain is elongated in the flow direction ($`\delta A_{xx}>0`$) and compressed in the in-plane gradient direction ($`\delta A_{yy}<0`$), while the neutral (out of plane) direction is only slightly decreased ($`\delta A_{zz}0`$).

## III Structure factor of the Rouse chain under shear flow

To compute the structure factor of an ideal gaussian chain (Rouse chain) we adopt the continuous chain model introduced by Edwards in which the only energy term is a quadratic potential acting between nearest neighbors along the chain. For such a model the equilibrium probability distribution in configurational space is

$$\mathrm{\Psi }[𝐑]exp\left(\frac{\chi }{2}_0^N𝑑n\left(\frac{𝐑_n}{n}\right)^2\right)$$
(5)
where $`\chi =3k_BT/d^2`$, $`𝐑_n`$ represents the positions of the chain monomers, $`0nN`$ is the continuous monomer index and $`d`$ is the single bond variance of the model. The steady state distribution in presence of a generic homogeneous solvent flow field can be obtained in terms of the normal modes of the model $`(𝐗_p,(p=0,\mathrm{}))`$. For homogeneous flows to which we limit our discussion, the normal modes are defined by the same transformation as at equilibrium. In those new coordinates the distribution is factorized over the normal modes

$$\widehat{\mathrm{\Psi }}(\{𝐗\},t)=\underset{p=0}{\overset{\mathrm{}}{}}\widehat{\psi }_p(𝐗_p,t)$$
(6)
In SF, the steady state distribution of the p-th mode is

$$\widehat{\psi }_p(𝐗_p)=\left(\frac{\chi _1}{2\pi k_BT}\right)^{3/2}\frac{p^3}{\sqrt{1+\dot{\gamma }^2\left(\frac{\tau _1}{2p^2}\right)^2}}e^{\frac{\chi _1}{2k_BT}p^2𝐗_p\beta _p^1𝐗_p}$$
(7)
where the matrix $`\beta _p`$ is defined as

$$\beta _p=\mathrm{𝟏}+\frac{\tau _1}{2p^2}(𝐤+𝐤^T)+2\left(\frac{\tau _1}{2p^2}\right)^2𝐤𝐤^T$$
(8)
and $`𝐤=𝐮`$ is the velocity gradient tensor. The derivation of eqs. (7) and (8) is given in appendix A.

The chain structure factor is defined as

$$S(𝐪,N,\dot{\gamma })=\frac{1}{N}_0^N_0^N𝑑n𝑑m<e^{i𝐪𝐑_{nm}}>_{\dot{\gamma }}$$
(9)
where $`𝐑_{nm}=𝐑_n𝐑_m`$ and $`<\mathrm{}>_{\dot{\gamma }}`$ indicates a statistical averages over the non equilibrium distribution eq.(6). In terms of normal modes we have

$$𝐑_{nm}=2\underset{p=1}{\overset{\mathrm{}}{}}\left(cos\frac{p\pi n}{N}cos\frac{p\pi m}{N}\right)𝐗_p$$
(10)
and using the model solution eq.(7) we obtain, after a gaussian integration,

$$<e^{i𝐪𝐑_{nm}}>_{\dot{\gamma }}=exp\left\{\frac{2k_BT}{\chi _1}𝐪\left[\underset{p=1}{\overset{\mathrm{}}{}}\frac{\beta _p}{p^2}\left(cos\frac{p\pi n}{N}cos\frac{p\pi m}{N}\right)^2\right]𝐪\right\}$$
(11)
Substituting the expression of $`\beta _p`$ from eq.(8) we get

$$<e^{i𝐪𝐑_{nm}}>_{\dot{\gamma }}=exp\left\{\frac{2k_BT}{\chi _1}\left[S_0(n,m)q^2+\frac{\tau _1}{2}S_1(n,m)𝐪(𝐤+𝐤^T)𝐪+\frac{\tau _1^2}{2}S_2(n,m)𝐪(𝐤𝐤^T)𝐪\right]\right\}$$
(12)
with

$$S_k(n,m)=\underset{p=1}{\overset{\mathrm{}}{}}\frac{1}{p^{2+2k}}\left[cos\frac{p\pi n}{N}cos\frac{p\pi m}{N}\right]^2k=0,1,2$$
(13)
In eq.(12) the first term in the exponent on the r.h.s. is the equilibrium contribution leading to the well known Debye function. The other two terms arise from the presence of the flow. The three series above can be summed analytically as explained in appendix B and the result is

$`S_0(a,b)`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{2}}|ab|`$ (14)
$`S_1(a,b)`$ $`=`$ $`{\displaystyle \frac{\pi ^4}{4}}(ab)^2\left[(a+b)\left(1{\displaystyle \frac{a+b}{2}}\right){\displaystyle \frac{|ab|}{3}}\right]`$ (15)
$`S_2(a,b)`$ $`=`$ $`{\displaystyle \frac{\pi ^6}{240}}(ab)^2\left\{|ab|^3+5(a+b)\left[(a+b)(2+a^2+b^2)3a^22ab3b^2\right]\right\}`$ (16)
where we have defined reduced indices $`a=n/N,b=m/N`$ and $`0(a,b)1`$.

Let us now introduce the reduced shear rate $`\beta `$ as the product of the bare shear rate and the global relaxation time of the chain at equilibrium. To compare the theory with the experiments we define $`\beta `$ as

$$\beta =\frac{M\eta _s[\eta ]\dot{\gamma }}{N_Ak_BT}$$
(17)
where $`M`$ is the molecular weight of the polymer, $`\eta _s`$ the solvent shear viscosity, $`[\eta ]`$ the intrinsic viscosity of the solution and $`N_A`$ the Avogadro number. In the Rouse model the intrinsic viscosity is $`[\eta ]=N_A/M\eta _s(N^2d^2\xi /36)`$ so that

$$\beta =\frac{\dot{\gamma }}{k_BT}\frac{N^2d^2\xi }{36}$$
(18)
which gives $`\tau _1\dot{\gamma }=12\beta /\pi `$. The chain structure factor in terms of the reduced shear rate is finally obtained as

$$\frac{S(𝐪,N,\dot{\gamma })}{N}=_0^1𝑑a_0^1𝑑be^{\frac{Nd^2}{6}\mathrm{𝐪𝐪}^T:𝐟(a,b,\beta )}$$
(19)
where $`𝐟(a,b,\beta )`$ is the following tensor

$$𝐟(a,b,\beta )=|ab|\mathrm{𝟏}+\frac{12\beta }{\pi ^4}\left[S_1(a,b)(\mathrm{\Lambda }+\mathrm{\Lambda }^T)+\frac{12\beta }{\pi ^2}S_2(a,b)(\mathrm{\Lambda }\mathrm{\Lambda }^T)\right]$$
(20)
and $`\mathrm{\Lambda }=\dot{\gamma }^1𝐤`$. Expressions (19) with (20) are our main interest here. The double integral over monomers cannot be performed analytically and we resort to numerical method. Before presenting specific results we want to make some further considerations. Firstly we note that the tensor $`𝐟`$ entering in the structure factor depends on $`N`$ and $`\dot{\gamma }`$ through $`\beta `$ only. This implies the following universality for the structure factor

$$\frac{S(𝐪,N,\dot{\gamma })}{N}=F(R_g𝐪,\beta )$$
(21)
where $`F`$ is a universal scalar function depending on the geometry of the applied flow and $`R_g=\sqrt{Nd^2/6}`$ is the equilibrium radius of gyration for the gaussian chain. Secondly we want to obtain the well known expression for the gyration tensor of the ideal chain under shear flow from which we can derive the orientation and deformation of the chain to be used below. The small $`q`$ expansion to order $`q^2`$ of the structure factor is directly related to the gyration tensor by

$$\frac{S(𝐪,N,\dot{\gamma })}{N}=1\mathrm{𝐪𝐪}^T:𝐆+O(q^4)$$
(22)
When applied to eq.(19) it provides

$$𝐆(N,\beta )=\frac{Nd^2}{6}_0^1𝑑a_0^1𝑑b𝐟(a,b,\beta )=\frac{Nd^2}{6}\left[\frac{1}{3}\mathrm{𝟏}+\frac{2\beta }{15}(\mathrm{\Lambda }+\mathrm{\Lambda }^T)+\frac{16\beta ^2}{105}(\mathrm{\Lambda }\mathrm{\Lambda }^T)\right]$$
(23)
Note that at fixed $`\beta `$ the gyration tensor has the same N dependence as at equilibrium (the scaling exponent is $`\nu =0.5`$). The orientation of the chain, as given by eq.(3), when defined through the gyration tensor is

$$cotg(2\chi _g)=\frac{G_{xx}G_{yy}}{2G_{xy}}=\frac{60\beta }{105}=\frac{\beta }{1.75}$$
(24)
implying that the orientational resistance $`m_g=1.75`$ and $`\beta `$-independent. This is the well known classical result.

Other used quantities in characterizing chains under SF are the deformation ratios as defined in section II. From eq.(23) we get $`\delta G_{xx}=16/35\beta ^2=0.457\beta ^2;\delta G_{yy}=\delta G_{zz}=0`$ as is well known. Note that the above values of the deformation ratios in the laboratory reference frame imply that in the molecular reference frame both elongation and compression directions present deformations as follows

$`\delta G_I`$ $`=`$ $`{\displaystyle \frac{2}{35}}\left[4\beta ^2+7\beta \sqrt{1+\left({\displaystyle \frac{4}{7}}\beta \right)^2}\right]`$ (25)
$`\delta G_{II}`$ $`=`$ $`{\displaystyle \frac{2}{35}}\left[4\beta ^27\beta \sqrt{1+\left({\displaystyle \frac{4}{7}}\beta \right)^2}\right]`$ (26)
with the following limiting behaviours

$`\delta G_I`$ $``$ $`{\displaystyle \frac{16}{35}}\beta ^2\beta >>{\displaystyle \frac{7}{4}}`$ (27)
$`\delta G_{II}`$ $``$ $`{\displaystyle \frac{49}{140}}=0.35\beta >>{\displaystyle \frac{7}{4}}`$ (28)
$`\delta G_I`$ $`=`$ $`\delta G_{II}{\displaystyle \frac{14}{35}}\beta \beta <<{\displaystyle \frac{7}{4}}`$ (29)
Now we go back to our original task, that is the calculation of the chain structure factor under shear flow. As already stated, the double integral over monomers in eq.(19) with eq.(20) cannot be performed analytically because of the exponential character of the integrand (to be compared with the double integral of the exponent which can be solved and provides the gyration tensor as in eq.(23)). We applied standard Romberg’s method  to numerically evaluate those integrals in a wide range of $`q`$’s and for given $`\beta `$. In figure 1 we report the chain structure factor in its universal form ($`S(𝐪)(qR_g)^2/3N`$ vs $`qR_g/\sqrt{3}`$), for various values of $`\beta `$ and for $`\widehat{𝐪}`$ along the elongation and the compression directions in the flow plane (molecular reference frame). We note as in the elongation direction the crossover to the high $`q`$’s equilibrium power law behavior strongly depends on the shear rate. The Kratky plot in the compression direction presents a $`\beta `$ dependent overshoot due to the compression of the corresponding gyration tensor component. This effects however saturates very soon (the curve for $`\beta =3.2`$ is already very close to the behavior for $`\beta =10`$ which is indistinguishable from the one at $`\beta =100`$) and the crossover to the high $`q`$’s power law behavior is in practice $`\beta `$ independent.

In comparison with the stretched ideal chain  we note the absence of oscillations in the elongation direction of the sheared chain signal. Such oscillations are characteristics of an aligned rod with uniform density of scattering centers. Our present results show definitively that, even at very high shear rates, the sheared chain cannot be thought in terms of such simplified model at variance with Onuki scheme .

Another important characteristics of the structure factor of the stretched chain, related to the presence of tensile blobs, is the Lorentzian shape of the scattering function at intermediate $`q`$’s as discussed by Benoit et al.. In this case one can easily show that, for large $`q`$, the scattering function for $`\widehat{𝐪}`$ along the stretching direction follows the equilibrium Ornstein-Zernike behavior $`S(q)=1/(1+q^2R_g^2/2)`$ if plotted in terms of the effective wave vector $`\stackrel{~}{q}^2=q^2+4/\xi _T^2`$, where $`\xi _T`$ is the tensile blob size. This led Pincus to formulate his scaling for stretched EV chains just changing the scaling exponent from the ideal to the EV value. Under shear flow the chain structure factor at intermediate $`q`$’s in the elongation direction has not a Lorentzian shape, but rather presents a behavior $`S(x)/N=1/(1+x^2/2+a(\beta )x^{\alpha (\beta )})`$ for $`x=qR_g/\sqrt{(}3)1`$, in analogy with the scattering signal from near critical fluids under shear. The fitting function $`a(\beta )`$ goes linearly from $`0`$ to $`73`$ in the interval $`0\beta 100`$, while $`\alpha (\beta )`$ goes almost linearly in $`log(\beta )`$ from $`0.62`$ to $`0.88`$ in the interval $`3.2\beta 100`$. Whether this form for $`S(q_I)`$ has some theoretical foundation remains unclear to us.

## IV Excluded volume chain model and Brownian Dynamics algorithm

The ideal chain model of previous section is only the starting point for a more realistic description of polymeric chains. When the aim is the study of dilute solutions at least two additional effects need to be taken into account namely the excluded volume and the hydrodynamic interactions between monomers. In this section we describe the model we have used to represent the excluded volume effect. In our previous work on stretched chains  we have used a model of linear chains with rigid bonds and no angle hindrance. The excluded volume interaction where modeled by hard sphere interactions. The diameter $`\sigma `$ of each monomer was chosen as $`\sigma =0.65`$ in units of the bond length. This model exhibits good solvent scaling and is very suitable for Monte Carlo calculations. Here we deal with non equilibrium systems for which the stationary distribution under shear is unknown and therefore we cannot apply the Monte Carlo sampling procedure. The only way to perform the calculation is by a dynamical technique such as Brownian Dynamics(BD). For this technique, the model of the previous work is not very appropriate because of the strong discontinuity in the interactions. We thus replaced the hard sphere by a Lennard-Jones potential truncated at the minimum and shifted in such a way to have a continuous potential

$`v_{EV}(r)`$ $`=`$ $`v_{LJ}(r)v_{LJ}(r_c)rr_c=2^{1/6}\sigma `$ (30)
$`=`$ $`0r>r_c`$ (31)
where $`v_{LJ}(r)=4ϵ[(\sigma /r)^{12}(\sigma /r)^6]`$. In our calculation $`ϵ`$ was the energy unit and $`\sigma `$ was chosen to be $`0.65`$ as in the previous work. Moreover the imposition of the bond constraints in BD requires an iterative technique such as SHAKE that for long chains is quite demanding in terms of computer time. We therefore prefer to use flexible rather then rigid bonds and represent the bond interaction by an harmonic potential with force constant $`\chi `$ and minimum position $`d`$: $`v_{bond}(r)=\chi (|𝐫|d)^2/2`$. Lennard-Jones interactions between nearest neighbors along the chain were not considered.

The equation of motion of the i-th monomer is

$$\xi \dot{𝐑}_i=𝐅_i+𝐟_i+\xi 𝐤𝐑_i$$
(32)
where $`𝐅_i`$ is the total force on monomer $`i`$ deriving from the potential energy, and the properties of the noise $`𝐟_i`$ are

$`<𝐟_i(t)>`$ $`=`$ $`0`$ (33)
$`<𝐟_n(t)𝐟_j(t^{})>`$ $`=`$ $`2k_BT\delta (tt^{})\delta _{ij}\mathrm{𝟏}`$ (34)
The above equations can be integrated numerically by the simple finite differences scheme

$$𝐑_i(t+h)=[\mathrm{𝟏}+h𝐤]𝐑_i(t)+\frac{h}{\xi }𝐅_i(t)+𝐰_i(t)+O(h^{3/2})$$
(35)
where the moments of the white noise $`𝐰_i`$ are

$`<𝐰_i(t)>`$ $`=`$ $`0`$ (36)
$`<𝐰_i(t)𝐰_j(t^{})>`$ $`=`$ $`2hD_0\delta (tt^{})\delta _{ij}\mathrm{𝟏}`$ (37)
and $`D_0=k_BT/\xi `$ is the single monomer diffusion coefficient. Higher order schemes has been proposed in the literature which allow to take a larger time step but require more operations per time step. A further advantage of the simple first order scheme over higher order ones, is the possibility to use a uniformly rather than normally distributed noise, provided that the zero and second moment are correctly chosen. Indeed, it is easy to show that the error introduced by such a substitution is of order $`h^{3/2}`$ or higher, that is beyond the precision of the scheme itself. Use of uniformly distributed noise rather than a gaussian noise saves about $`50\%`$ of the computer time. In absence of a systematic comparison on the overall efficiency of integration schemes we have preferred to use the simple one above.

In all our calculations we have fixed $`k_BT=1`$ and $`\xi =1`$ which lead to $`D_0=1`$. The bond force constant has been chosen to be $`\chi =40`$ and in order to obtain a unitary average bond length we had to fix $`d=0.945`$.

The presence of stiff springs, ensuring the rigidity of the bonds, and strongly repulsive core interactions, ensuring the self avoidance and avoiding the self crossing of the chain, impose a quite small time step for the stability of the numerical scheme. On the other hand the long relaxation time of polymers and its rapid increase with the number of monomers ($`\tau _NN^{2.2}`$ for excluded volume chains without HI) would require a time step as large as possible. This is the main limitation of any dynamical scheme when dealing with polymers. In the present work we have chosen a time step $`h=0.00025`$ in reduced units. We have tested in some cases that halving the time step does not change the results for the chain structure. Even with such small time step, numerical instabilities occasionally show up. These occur when two particles end a time step at a reciprocal distance considerably smaller then $`\sigma `$. The following time step is then driven by very large repulsive forces which stretch the bonds related to the pair of particle at a distance considerably larger then the minimum of the bond potential and so forth. To cure this occasional pathology we checked at the end of any time step the minimum distance between any pair of particles and if it is found to be smaller then $`0.43\sigma `$ we reject the time step and we proceed the time integration for a time interval $`h`$ in many smaller time steps. When the time interval $`h`$ is reached the original time step is restored. With the chosen value for the original time step we have found that such pathological events happen on average any $`2x10^5`$ time steps and that dividing the time step by ten is always enough to escape from the instability. This is a craft made version of more sophisticated adaptive time step integration schemes which are however more demanding in terms of computer time.

To further minimize the computer time per step we have used a Verlet neighbors list  which is very effective with such short range interactions.

## V Simulation results

### A Equilibrium scaling of the model

Before proceeding to study the chain model under shear flow we need to characterize its equilibrium behavior. In particular we need to know at which length the chains start following the known static and dynamic scalings and we need to determine the prefactor in the scaling laws in order to fix the values of the bare shear rate for a given chain length to work at fixed $`\beta `$. In figure 2 we report the radius of gyration $`R_g`$ and the end-to-end distance $`R`$ versus the number of links in the chain $`(N1)`$. The observed scalings are $`R_g=0.40448(N1)^{0.6}`$, and $`R=1.03(N1)^{0.6}`$ if we exclude $`N=9,20`$. Equilibrium structure factors for the various chain lengths are shown in figure 3a. In figure 3b we report the universal plot $`(qR_g)^2S(q)/N`$ versus $`(qR_g)`$. We also add the Debye curve for the ideal chain. We observe that, for $`N30`$, the EV chains follow indeed a universal $`q^{1/3}`$ power law as predicted on the basis of the static scaling and dimensional arguments .

By construction, the chain centre of mass diffusion coefficient is $`D_{cm}=D_0/NN^\nu ^{}`$ where $`\nu ^{}`$ is the dynamic exponent, which provides a characteristic time of the chain $`\tau _N=R_g^2/(6D_{cm})=(0.40448)^2(N1)^{1.2}N/6D_0N^{2\nu +1}`$. Here $`\nu =3/5`$ is the Flory scaling exponent for EV chains, and $`\nu ^{}=1`$ as HI are absent. In order to check the consistency between global and internal dynamics we extracted the characteristic time of the latter by inverting the chain dynamic structure factor $`I(q,t)=S(q,t)/S(q)`$ as we have previously done for MD chains . Indeed scaling and dimensional arguments predict that $`I(t,q,N)`$ be a universal function of $`(tq^x)`$ with $`x=2+\nu ^{}/\nu =11/3`$ in the present case. Figure 4 shows the behavior of $`(tq^{11/3})`$ vs $`q`$ for two values of $`I`$ and for various chain lengths. We observe that for $`N=20`$ the chain does not exhibit the correct scaling of internal times while from $`N=30`$ on, both static and dynamics are in the scaling laws regime. The data we show are results of quite long runs as reported in table I.

### B Chain structure at fixed $`\beta `$

In this subsection we describe the results for the chain structure under shear flow for various chain lengths at the same value of the reduced shear rate $`\beta `$ that we defined as $`\beta =\dot{\gamma }R_g^2/(6D_{cm})`$. We have limited our study to two values of $`\beta `$, namely $`\beta =3.2`$ as in our previous MD study and in several experiments, and $`\beta =10`$. In tables II and III, we report various technical details of our simulations and the results for the orientation and the deformation. We studied chains up to $`N=300`$ which require very long runs (up to $`310^9`$ time steps for the EV case, see table III).

As in our previous MD study of short chains (up to $`N=50`$), we analyze the results in the molecular reference frame on account of the observation that the orientational angle depends on $`N`$ and $`\dot{\gamma }`$ through $`\beta `$ only. This is also suggested by the Rouse model results of section III and by its extension when HI at the level of equilibrium preaveraging are considered . Our new results for EV chains at $`\beta =3.2`$ confirm this expectation within error bars, although short chains seems to be slightly less oriented as shown in figure 5. On the same figure we report the extinction angle related to the birefringence which also is found to be the same within error bars for all $`N`$ at fixed $`\beta `$.

Before discussing the results for our chain model with EV interactions it is instructive to consider the same coarse grained model with the EV interactions switched off ($`ϵ=0`$ in the LJ potential). Comparison of the chain structure factor of such a model with the continuous ideal model result eq.(19) (not shown here) validates the theory itself and the numerical procedure applied to compute the double integral. Moreover data for the orientation and the deformation of this model can be useful in understanding qualitatively the role of the finite extensibility effects due to the stiff harmonic potential which keeps the bond length nearly constant. Data for the orientation and deformations of ideal chains with $`N=100`$ and $`N=300`$ at $`\beta =3.2`$ and $`\beta =10`$ are reported in table II and compared to the Rouse model predictions. At $`\beta =3.2`$, the orientation for all chain lengths is in agreement with the long chain limit while the deformations for $`N=100`$ deviate considerably from the ones for $`N=300`$ (except in the direction II) which are instead in agreement with the Rouse prediction. Interestingly, the main eigenvalue $`G_I`$ of the gyration tensor follows an apparent power law $`N^{1.5}`$ qualitatively similar to the behavior observed for MD data on shorter chains with EV and HI. At $`\beta =10`$ both chains exhibit finite extensibility effects even on the orientation.

Now we turn to the EV results. In table III we collect the data for the orientation and the deformation of our chains. We also report data for the longest chains studied in ref.  by molecular dynamics. We note at $`\beta =3.2`$, that the 300-bead chain appears to be slightly more oriented and deformed than the shorter chains, although the noise is quite large. At $`\beta =10`$ instead, the 300-bead chain is considerably more oriented and deformed than the 100-bead chain, showing an important influence of the finite extensibility. Comparison with previous MD data for N=50 shows a general agreement for all quantities.

From these results and those of table II, it appears reasonable to assume that the N=300 results at $`\beta =3.2`$ are representative of the global properties of long chains, namely that the remaining small systematic finite chain effects are masked by statistical errors of a few percents (due to finite statistics). Comparing EV chains to $`\theta `$ chains at the same reduced shear rate, we observe for EV chains a less marked alignment with respect to the flow (larger $`\chi `$ value) and a smaller global deformation, in qualitative agreement with recent LS results .

An important question concerns the universal character of the structure factor under shear flow. It is well known that at equilibrium the quantity $`S(q,N)/N`$ is a universal function of $`qR_g`$ only, for any coarse-grained model and for long enough chains. In the case of the chain stretched at both ends, we have shown that different universalities in terms of $`qR_g`$ hold separately in the longitudinal and the transverse directions. Alternatively one could investigate self-similarities in terms of $`q\sqrt{G_\alpha }`$ where $`G_\alpha `$ are the diagonal components of the gyration tensor. For the chain stretched by its ends, it is well known that fluctuations of the end-to-end distance (and therefore $`G_\alpha `$) at fixed reduced force scale with $`N`$ as the equilibrium radius of gyration and therefore the above two alternatives are equivalent. In our previous study of linear chains under shear flow we observed apparent scaling exponents $`\nu _\alpha `$ for the eigenvalues of the gyration tensor with $`(N1)`$ which were different from the equilibrium one in the flow plane directions while in the neutral (out of plane) direction the equilibrium scaling was preserved within error bars. In the present study on longer chains we find that at $`\beta =3.2`$ the equilibrium scaling $`G_\alpha N^{2\nu }`$ is preserved in the neutral direction ($`\alpha =III`$) and in the in-plane compression direction ($`\alpha =II`$) (see figure 6). In the elongation direction an apparent scaling exponent higher than the equilibrium value is still observed ($`\nu _I0.7`$) although the noise is considerably higher than in the other two principal directions. The results for the chain without EV interactions mentioned above demonstrate that such apparent scaling is related to finite extensibility effects rather than to EV interactions. For this reason we investigate universality in terms of $`qR_g`$.

In figure 7a,b,c we plot, at $`\beta =3.2`$ and for $`30N300`$, $`(qR_g/\sqrt{3})^{5/3}S(q)/N`$ versus $`qR_g/\sqrt{3}`$ for $`\alpha =III,II,I`$ respectively and we compare with the MD data for N=50 at the same $`\beta `$ from ref.. For the latter we have used the observed equilibrium exponent $`\nu _{MD}=0.57`$ rather then the Flory classical value $`\nu _F=0.6`$. Note that in that model the solvent was explicitly considered so that deviations from the mean field value are not contradictory. We observe that, with the exception of $`N=30`$ in the neutral direction, universality is indeed observed within error bars in all directions. A general feature of those curves is the presence of a quite slow crossover back to the equilibrium behavior at large $`q`$’s which would be represented by an horizontal line. That $`N=300`$ is still too short to cover all length scales of the problem is clearly shown by the fact that even for this length the equilibrium behavior is not completely recovered before the chain finite extensibility shows up, a sign that finite chain length effects are still present. While in the neutral direction the effect of the flow is tiny, in the in-plane directions the flow strongly changes the scattering function. We note that the noise level is considerable higher in the elongation direction than in the other directions. Moreover statistical convergence of the results at low $`\beta `$ appears to be considerably slower than at high $`\beta `$ and at equilibrium (see also table III). This is compatible with the possible appearance of new characteristic relaxation times in weak flow as indicated by recent experiments .

Finally we note that the structure factor obtained in our previous MD simulation on a different chain model embedded in an explicit atomic solvent under shear flow , follows very closely in all directions the universal behavior predicted by the present BD simulations. As the former study includes automatically the hydrodynamics interactions  while the present calculation do neglect them from the start, it appears that HI play a secondary role on the chain structure under shear flow, once one works at a fixed value of the suitably defined reduced shear rate. Because of the limited length of the largest chain studied by molecular dynamics (N=50), we could not detect on the chain structure factors the slow trend back to the equilibrium behavior at large q’s, which explains our previous analysis in terms of anisotropic scaling laws.

In order to test the “shear” blob hypothesis we report in figures 8,9,10, $`(qR_g)^2S(q)/3N`$ versus $`qR_g/\sqrt{3}`$ for $`N=100`$ and $`N=300`$ at $`\beta =3.2`$ and for $`N=100`$ at equilibrium, and we compare with the ideal chain behavior. In all directions Rouse behavior describes the excluded volume chain up to $`qR_g/\sqrt{3}3`$. Similarly to the stretched chain case , we clearly see in the compression and the neutral directions a crossover from ideal statistics at low $`q`$’s to EV statistics at high $`q`$’s. The inset of EV statistics is at $`qR_g/\sqrt{3}9`$ in the two directions. Note that for $`N=100`$ the EV regime is completely missing and even with $`N=300`$ it is very narrow in the compression direction. In the elongation direction the noise level is quite high in the Rouse (small $`q`$) regime. As stated above, in this direction flow effects extends to higher $`q`$’s than in the other directions, above $`qR_g/\sqrt{3}=9`$. On the other hand, the finite extensibility takes place around $`qR_g/\sqrt{3}=20`$ for $`N=300`$ and therefore the EV regime cannot be observed in this direction for such chain length.

Onuki’s phenomenological model  is based on the assumption that below a characteristic length scale $`\xi `$ the effects of the flow are negligible. The length scale separation is dictated by the typical relaxation time of chain density fluctuations at that scale. At a given shear rate $`\dot{\gamma }`$, the crossover occurs at a scale $`\xi =An_c^\nu `$ such that the longest relaxation time of the blob is $`\tau _c=\xi ^2/(6D_c)=\dot{\gamma }^1`$ where $`D_c=D_0/n_c^\nu ^{}`$ is the diffusion of the blob as a whole and $`\nu ^{}`$ is the dynamic scaling exponent. Working at fixed reduced shear rate $`\beta =\tau _N\dot{\gamma }`$, the number of monomers per blob results in $`n_c=N/\beta ^{1/(2\nu +\nu ^{})}`$ so that the number of blobs per chain is $`N_b=N/n_c=\beta ^{1/(2\nu +\nu ^{})}`$ and $`\xi =AN^\nu /\beta ^{\nu /(2\nu +\nu ^{})}=R_g/\beta ^{\nu /(2\nu +\nu ^{})}`$. For our model ($`A=0.40448,\nu =3/5,\nu ^{}=1`$) at $`\beta =3.2`$ we have $`N_b=1.7,n_c=177`$ and $`\xi =9.02`$ for $`N=300`$. The crossover from ideal to EV statistics in figure 8 is observed around $`qR_g/\sqrt{3}=4.0`$ which for $`N=300`$ provides $`q_c=0.565/\xi 2\pi /\xi `$. Therefore the shear blob hypothesis is compatible with our results although a more firm test on the basis of longer chains studied at various values of $`\beta `$ is still missing. To finally prove that our chains are too short for this purpose we show in figure 11 the Kratky plots in the three principal directions for $`N=300`$ and $`\beta =3.2`$ and $`\beta =10.0`$. It is evident that at $`\beta =10.0`$ the finite extensibility effect interferes with the EV regime giving rise to anomalous behaviours.

## VI Conclusions

In order to study the structure of a single chain subjected to a homogeneous and steady shear flow and to elucidate the effects of excluded volume interactions we have computed the scattering function of the continuous Rouse chain model in flow and we have performed BD simulations of a chain model of Fraenkel springs with fully developed EV interactions.

Concerning global properties, our simulations confirm that at a given reduced shear rate, EV chains are less deformed and less oriented along the flow lines than chains at the $`\mathrm{\Theta }`$ point .

We have shown that, even in the ideal case, the chain structure in shear flow cannot be explained by the known behavior of a chain stretched at its ends. In particular we have shown that the structure factor of the Rouse model under SF for the scattering vector lying along the in plane elongation direction is very different from the one of the Rouse chain stretched at its ends with $`q`$ along the stretching direction. It does not present the oscillations typical of permanent extension, but it is a monotonous curve qualitatively similar to the equilibrium Debye curve (see figure 10). In the elongation direction, the onset of the well known $`q^2`$ behavior signaling the static scaling of ideal chains is shifted to higher $`q`$’s than at equilibrium and is strongly dependent on the shear rate. In the compression direction it is at the same location as at equilibrium and it is almost insensible to shear rate.

Adding the EV interactions, our results suggest that the universality of $`S(𝐪,N,\dot{\gamma })`$ in terms of $`x=𝐪R_g`$ and $`\beta `$ only ($`R_g`$ is the equilibrium radius of gyration and $`\beta `$ is computed on the basis of the equilibrium characteristic relaxation time), exhibited by the Rouse model at any flow intensity, is preserved by the EV interactions provided the chains are long enough. Studying chains between 30 and 300 monomers we indeed observe that the eigenvalues of the gyration tensor in the compression (II) and in the neutral (III) directions follow scaling laws with N with the equilibrium Flory exponent. In the elongation direction a larger apparent scaling exponent is instead measured (around 0.7), in agreement with previous MD data which led us to propose an anisotropic scaling picture for intermediate reduced shear rates. From the results for the same chain model without EV interactions we infer that this anomalous scaling arises from finite extensibility effects rather than EV interactions so as the equilibrium scaling should be recovered even in the elongation direction for longer chains. Unfortunately the dynamical character of the BD algorithm and the unfavorable scaling with N of the chain relaxation time does not allow us to study chains longer than $`N=300`$.

The general agreement found between our present data (where HI are absent from the model) and the previous non-equilibrium MD data for short chains (where HI are automatically included) strongly suggests that HI effects play a minor role on the chain structure under shear flow, provided we compare results at the same reduced shear rate, which incorporates implicitly the HI effect through the longest relaxation time of the chain at equilibrium. The EV and the finite extensibility effects are by far the most relevant effects to be taken into account.

For our longest chain at $`\beta =3.2`$, the structure factor in the neutral (III) and in the compression (II) directions exhibit a clear crossover from the ideal behavior at small $`q`$’s, well predicted by the Rouse model, to the isotropic equilibrium EV behavior at high $`q`$’s, characterized by a power law $`S(q)q^{1/\nu }`$ with $`\nu =\nu _{eq}=3/5`$. This behavior is usually related to the presence of the so called “blobs”, in the present case“shear blobs”. Blobs were detected in the chain structure factor either experimentally for semi-dilute solutions  or by Monte Carlo simulations for a single EV chain stretched at its ends . The presence of blobs in a chain under shear flow has been predicted theoretically  but had never been observed so far, neither by SANS nor by simulation. The results of our BD simulations are compatible with the existence of blobs (see section V B) but, because of finite chain size effects, a quantitative test of this hypothesis would require chain lengths larger than N=300. It is interesting to note that the largest chain considered here has a size which corresponds to a polystyrene chain of the order of 300 Kuhn segments of about $`18\AA `$, which yields a molecular mass of 300.000 Daltons .

The present analysis should be tested by scattering experiments on dilute solutions for chains in good solvent subjected to steady shear flow. It would be interesting to study at low $`qR_g`$, either by Light or by Neutron scattering, whether the Rouse model indeed describes the global orientation and deformation of chains for large deformations, i.e. at $`\beta >>1`$. Remaining discrepancies could be interpreted as a direct measure of FE and HI influence on these structural properties far from the equilibrium structure. Ultimately, the most spectacular experimental result we seek for remains a SANS study of a suitable “polymer-good solvent” pair which would allow us to follow the form factor of a polymer of a few million Daltons under shear flow at large $`\beta `$ over a broad $`q`$ regime. This would allow a direct observation of the crossover from ideal chain to EV chain statistics across a crossover diffusion vector $`q_c`$ characterizing the shear blob dimension in a particular direction.

Concerning future simulation work, a more quantitative test of HI effects on the chain structure is in order and could be performed by BD on short chains.

## VII Acknowledgments

We thank M. Baus, B. Dünweg, P. Lindner and R.Winkler for useful discussions and suggestions.

## A Ideal chain under steady shear flow

In this appendix we derive the steady state distribution function of the continuous gaussian chain model under an homogenous shear flow, eqs. (6), (7) and (8).

The dynamics of the model under a generic flow field, specified by the velocity gradient tensor $`𝐤`$, is defined by the Langevin equation

$$\xi \frac{𝐑_n(t)}{t}=\chi \frac{^2𝐑_n(t)}{n^2}+𝐟_n(t)+\xi 𝐤(t)𝐑_n(t)$$
(A1)
with the boundary conditions

$$\left[\frac{𝐑_n(t)}{n}\right]_{n=0}=\left[\frac{𝐑_n(t)}{n}\right]_{n=N}=0$$
(A2)
$`𝐟_n(t)`$ is the gaussian random force acting on the chain and is defined by

$`<𝐟_n(t)>`$ $`=`$ $`0`$ (A3)
$`<𝐟_n(t)𝐟_m(t^{})>`$ $`=`$ $`2k_BT\delta (tt^{})\delta (nm)\mathrm{𝟏}`$ (A4)
where $`<\mathrm{}>`$ indicates an average over the distribution of the random noise and $`\mathrm{𝟏}`$ is the unit tensor. As it is well known the above Langevin dynamics corresponds to the following diffusion equation (Smoluchowsky) for the distribution function in the chain configurational space

$$\frac{\mathrm{\Psi }}{t}=_0^N𝑑n_{𝐑_n}\left[\dot{𝐑}_n\mathrm{\Psi }\right]$$
(A5)
with the dynamics of monomers given by

$$\dot{𝐑}_n=\frac{\chi }{\xi }\frac{^2𝐑_n}{n^2}\frac{k_BT}{\xi }_{𝐑_n}[ln\mathrm{\Psi }]+𝐤𝐑_n$$
(A6)
Let us introduce now the normal modes defined as

$$𝐗_p(t)=\frac{1}{N}_0^N𝑑ncos\left(\frac{p\pi n}{N}\right)𝐑_n(t)p=0,1,2,\mathrm{}.$$
(A7)
in terms of which the chain positions are given as

$$𝐑_n(t)=𝐗_0(t)+2\underset{p=1}{\overset{\mathrm{}}{}}cos\left(\frac{p\pi n}{N}\right)𝐗_p(t)0nN$$
(A8)
The completeness relations for the normal modes basis are

$`\delta _{pq}={\displaystyle \frac{2}{N}}{\displaystyle _0^N}𝑑ncos\left({\displaystyle \frac{p\pi n}{N}}\right)cos\left({\displaystyle \frac{q\pi n}{N}}\right)`$ (A9)
$`N\delta (nm)=1+2{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}cos\left({\displaystyle \frac{p\pi n}{N}}\right)cos\left({\displaystyle \frac{p\pi m}{N}}\right)`$ (A10)
Note that, for homogeneous flows which we consider here, the normal modes of the dynamics are the same as at equilibrium.

The Langevin equations for the normal modes are

$$\xi _p_t𝐗_p(t)=\chi _p𝐗_p(t)+\stackrel{~}{𝐟}_p(t)+\xi _p𝐤𝐗_p(t)$$
(A11)
with the following definition of the symbols

$`\xi _0`$ $`=`$ $`N\xi ,\xi _p=2N\xi p=1,2,3\mathrm{}`$ (A12)
$`\chi _p`$ $`=`$ $`{\displaystyle \frac{6\pi ^2k_BT}{Nd^2}}p^2={\displaystyle \frac{2\pi ^2\chi }{N}}p^2=\chi _1p^2`$ (A13)
$`<\stackrel{~}{𝐟}_p(t)>`$ $`=`$ $`0,<\stackrel{~}{𝐟}_p(t)\stackrel{~}{𝐟}_q(t^{})>=2\xi _pk_BT\delta _{pq}\delta (tt^{})\mathrm{𝟏}`$ (A14)
The distribution function of normal modes is factorized

$$\widehat{\mathrm{\Psi }}(\{𝐗\},t)=\underset{p=0}{\overset{\mathrm{}}{}}\widehat{\psi }_p(𝐗_p,t)$$
(A15)
and one can easily derive the following diffusion equation for the $`pth`$ mode ($`p>0`$)

$$_t\widehat{\psi }_p(𝐗_p,t)=_{𝐗_p}\left[(𝐤𝐗_p\frac{\chi _p}{\xi _p}𝐗_p)\widehat{\psi }_p\frac{k_BT}{2\xi }_{𝐗_p}\widehat{\psi }_p\right]$$
(A16)
The center of mass motion is instead described by the diffusion equation for the zero mode

$$_t\widehat{\psi }_0(𝐗_0,t)=_{𝐗_0}\left[𝐤𝐗_0\widehat{\psi }_0\frac{k_BT}{2\xi }_{𝐗_0}\widehat{\psi }_0\right]$$
(A17)
The solution of eqs.(A16) is a multivariate guassian

$$\widehat{\psi }_p(𝐗_p,t)=\left(\frac{1}{(2\pi )^3det[\alpha _p(t)]}\right)^{1/2}e^{\frac{1}{2}𝐗_p\alpha _p^1(t)𝐗_p}$$
(A18)
where the time dependent covariance matrix $`\alpha _p(t)`$ is given by (see eq.(7.161) in ref.)

$$\alpha _p(t)=_{\mathrm{}}^t𝑑t^{}\frac{2k_BT}{\xi _p}𝐄(t,t^{})𝐄^T(t,t^{})e^{\frac{2p^2}{\tau _1}(tt^{})}$$
(A19)
Here $`\tau _1=\xi _1/\chi _1`$ is the relaxation time of the first normal mode and the tensor $`𝐄(t,t^{})`$ is the deformation tensor of the external field defined by

$$𝐄(t,t^{})=\mathrm{𝟏}+_t^{}^t𝑑t^{\prime \prime }𝐤(t^{\prime \prime })$$
(A20)
For the steady shear flow to which we are interested in this work we have

$$𝐤=\left(\begin{array}{ccc}0& \dot{\gamma }& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)$$
(A21)
$$𝐄(t,t^{})=\left(\begin{array}{ccc}1& \dot{\gamma }(tt^{})& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)$$
(A22)
so that the covariance matrices $`\alpha _p`$ results time independent and takes the form

$$\alpha _p=\frac{k_BT}{\chi _p}\left[\mathrm{𝟏}+\frac{\tau _1}{2p^2}(𝐤+𝐤^T)+2\left(\frac{\tau _1}{2p^2}\right)^2𝐤𝐤^T\right]=\frac{k_BT}{\chi _1}\left(\frac{\beta _p}{p^2}\right)$$
(A23)
where we have defined the matrix $`\beta _p`$ through the second equivalence (see eq.(8)). Noting that

$$det(\alpha _p)=\left(\frac{k_BT}{\chi _1p^2}\right)^3\left[1+\dot{\gamma }^2\left(\frac{\tau _1}{2p^2}\right)^2\right]$$
(A24)
the distribution function of the $`pth`$ normal mode under steady shear flow results in eq.(7). It is easy to check that eq.(7) reduces to the standard equilibrium distribution for $`\dot{\gamma }=0`$.

## B Derivation of the sums

In this appendix we outline the method to compute the sums in eqs. (13)

$$S_k(a,b)=\underset{p=1}{\overset{\mathrm{}}{}}\frac{1}{p^{2+2k}}\left[cos(p\pi a)cos(p\pi b)\right]^2k=0,1,2$$
(B1)
where $`0(a,b)1`$. It is evident that $`S_k(a,b)=S_k(b,a)0`$ and the equality holds for $`a=b`$ only. Developing the square inside the sum in eq.(B1) and using standard relations among trigonometric functions we obtain

$$S_k(a,b)=\sigma _k(0)+\frac{1}{2}\left[\sigma _k(a)+\sigma _k(b)\right]\sigma _k\left(\frac{a+b}{2}\right)\sigma _k\left(\frac{ab}{2}\right)$$
(B2)
with

$$\sigma _k(x)=\underset{p=1}{\overset{\mathrm{}}{}}\frac{cos(2\pi px)}{p^{2+2k}}$$
(B3)
so that the problem reduces to compute $`\sigma _k(x)`$ for $`k=0,1,2`$. Note that $`\sigma _k(x)`$ is an even and periodic function of $`x`$ with unitary period. Eq.(B3) expresses the Fourier series a function $`f_k(x)`$ with these properties, after removing the $`p=0`$ term. Therefore we have to seek for a function $`f_k(x)`$ that satisfies the following conditions

$`2p^{2+2k}{\displaystyle _0^1}𝑑xf_k(x)cos(2\pi px)`$ $`=`$ $`1p1`$ (B4)
$`{\displaystyle _0^1}𝑑xf_k(x)sin(2\pi px)`$ $`=`$ $`0`$ (B5)
Applying repeatedly integration by parts we easily obtain the following identities

$`{\displaystyle _0^1}𝑑xf_k(x)cos(2\pi px)`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi p)^2}}{\displaystyle _0^1}𝑑xf_k^{}(x)_xcos(2\pi px)`$ (B6)
$`=`$ $`{\displaystyle \frac{1}{(2\pi p)^4}}{\displaystyle _0^1}𝑑xf_k^{\prime \prime \prime }(x)_xcos(2\pi px)`$ (B7)
$`=`$ $`{\displaystyle \frac{1}{(2\pi p)^6}}{\displaystyle _0^1}𝑑xf_k^𝚟(x)_xcos(2\pi px)`$ (B8)
having imposed the boundary conditions ($`f_k(0)=f_k(1)`$) for the function and for its derivatives up to order five.

Let us now consider the three $`k`$ values separately. For $`k=0`$ we have

$`1`$ $`=`$ $`{\displaystyle \frac{1}{2\pi ^2}}{\displaystyle _0^1}𝑑xf_0^{}(x)_xcos(2\pi px)`$ (B9)
$`0`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi p)^2}}{\displaystyle _0^1}𝑑xf_0^{}(x)_xsin(2\pi px)p1`$ (B10)
$`f_0(0)`$ $`=`$ $`f_0(1)`$ (B11)
where we have used the first of the equalities in eq.(B6) and a similar one for the sinus function. These relations can be satisfied by choosing $`f_0(x)=c_0+c_1x+c_2x^2`$ with $`c_0=0,c_1=\pi ^2`$ and $`c_2=\pi ^2`$ which provides $`f_0(x)=\pi ^2(x^2x)`$. The function $`\sigma _0(x)`$ is then

$$\sigma _0(x)=f_0(x)_0^1𝑑xf_0(x)=\pi ^2(x^2x\frac{1}{6})$$
(B12)
and, after eq.(B2)

$$S_0(a,b)=\frac{\pi ^2}{2}|ab|$$
(B13)
which is the wanted result for $`k=0`$.

For $`k=1`$ we must solve the problem

$`1`$ $`=`$ $`{\displaystyle \frac{1}{8\pi ^4}}{\displaystyle _0^1}𝑑xf_1^{\prime \prime \prime }(x)_xcos(2\pi px)`$ (B14)
$`0`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi p)^4}}{\displaystyle _0^1}𝑑xf_1^{\prime \prime \prime }(x)_xsin(2\pi px)p1`$ (B15)
$`f_1(0)`$ $`=`$ $`f_1(1)`$ (B16)
where the second equalities in eq.(B6) has been used. These relations can be satisfied choosing $`f_1^{\prime \prime \prime }(x)=c_3+c_4x`$ with $`c_4=8\pi ^4`$ and $`c_3=4\pi ^4`$. Integrating three times over $`x`$ and imposing the boundary conditions on $`f_1^{\prime \prime },f_1^{}`$ and $`f_1`$ we obtain $`f_1(x)=\pi ^4(2x^3x^2x^4)/3`$ which provides

$$\sigma _1(x)=f_1(x)_0^1𝑑xf_1(x)=\frac{\pi ^4}{3}(\frac{1}{30}x^2+2x^3x^4)$$
(B17)
Using eq.(B2) we finally obtain

$$S_1(a,b)=\frac{\pi ^4}{4}(ab)^2\left[(a+b)\left(1\frac{a+b}{2}\right)\frac{|ab|}{3}\right]$$
(B18)
In the case $`k=2`$ the relation to be satisfied are

$`1`$ $`=`$ $`{\displaystyle \frac{2}{(2\pi )^6}}{\displaystyle _0^1}𝑑xf_1^𝚟(x)_xcos(2\pi px)`$ (B19)
$`0`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi p)^6}}{\displaystyle _0^1}𝑑xf_1^𝚟(x)_xsin(2\pi px)p1`$ (B20)
$`f_2(0)`$ $`=`$ $`f_2(1)`$ (B21)
A procedure similar to the previous cases provides $`f_2(x)=\pi ^6x^2(2x^46x^3+5x^21)/45`$ and therefore

$$\sigma _2(x)=\frac{\pi ^6}{45}\left[\frac{5}{7}+x^2(2x^46x^3+5x^21)\right]$$
(B22)
and

$$S_2(a,b)=\frac{\pi ^6}{240}(ab)^2\left\{|ab|^3+5(a+b)\left[(a+b)(2+a^2+b^2)3a^22ab3b^2\right]\right\}$$
(B23)
Corresponding author:
Carlo Pierleoni
Physics Department, University of L’Aquila
Via Vetoio, Coppito 67010-L’Aquila (Italy)
phone: +39+0862+433056
fax: +39+0862+433033
email: carlo.pierleoni@aquila.infn.it