Molecular dynamics study of shear-induced long-range correlations in simple fluids

We investigate long-range correlations (LRCs) induced by shear flow using the molecular dynamics (MD) simulation. We observe the LRCs by comparing the MD results with the linearized fluctuating hydrodynamics (LFH). We find that the MD result has large finite-size effects, and it prevents the occurrence of LRCs in small systems. We examine the finite-size effects using a sufficiently large system consisting of more than ten million particles, and verify the existence of shear-induced LRCs without ambiguity. Furthermore, we show that MD result is quantitatively consistent with the LFH solution for the large system. As we reduce the system size $L$ or increase the shear rate $\dot{\gamma}$, the hydrodynamic description gradually breaks down in the long-wavelength region. We define a characteristic wavenumber $k^{\rm vio}$ associated with the breakdown and find the nontrivial scaling relations $k^{\rm vio} \propto L^{-\omega}$ and $k^{\rm vio} \propto \dot{\gamma}$, where $\omega$ is an exponent depending on $\dot{\gamma}$. These relations enable us to estimate the finite-size effects in a larger-size simulation from a smaller system.


I. INTRODUCTION
For equilibrium systems with short-range interactions, long-range correlations (LRCs) appear in certain situations, such as for a critical point and for an ordered phase with spontaneous symmetry breaking. Near the critical point, the correlation length diverges and the correlation function exhibits a long-range nature [1]. In the ordered phase with spontaneous breaking of the continuous symmetry, the so-called Nambu-Goldstone mode [2][3][4] appears and leads to the LRC.
The mechanism and nature of nonequilibrium LRCs have been established using phenomenological models such as fluctuating hydrodynamics [24][25][26][27][28] and stochastic lattice gases [29][30][31]. These coarse-grained models enable the nonequilibrium fluctuations to be examined in terms of the violation of a detailed-balance condition, conservation law, and anisotropy. However, they do not produce the mechanism of nonequilibrium LRCs from molecular-scale dynamics. There are few theoretical attempts to study nonequilibrium LRCs from the underlying Hamiltonian dynamics. Then, how LRCs arise from the molecuar-scale dynamics remains poorly understood.
In the hydrodynamic description, nonequilibrium fluctuations consist of two terms: where A(r) and B(r) are density fields (e.g., density fluctuations δρ(r) or velocity fluctuations δv a (r)), and c 1 and c 2 are appropriate constants. The first term proportional to the delta function implies that A(r) and B(r) are uncorrelated on the hydrodynamic scale. Then, the second term represents the LRCs, which are generally absent in equilibrium fluids.
The fluctuating hydrodynamics provides a phenomenological model for describing the fluctuations at the hydrodynamic scale. This model is widely used to study shearinduced LRCs. One characteristic behavior of shearinduced LRCs is the crossover between two power-law decays [21,26]. For example, the spatial correlation of the density fluctuation δρ(r)δρ(r ′ ) decays according to |r − r ′ | −1 for short-distance scales and crosses over to the stronger decay |r −r ′ | −11/3 for long-distance scales. Similarly, the spatial correlation of the velocity fluctuations v(r) · v(r ′ ) crosses over from |r − r ′ | −1 to |r − r ′ | −5/3 . Another important prediction from the fluctuating hydrodynamics is the existence of shear-induced corrections to the pressure P and shear viscosity η. These corrections arise from the nonlinear coupling of the LRCs. Kawasaki and Gunton initially found these corrections by using the projection operator method and the mode-coupling theory [52]. They were subsequently derived from the fluctuating hydrodynamics [15,21,53]. The shear-induced correction to the pressure P depends on the Reynolds number Re, and is given in the two limits as P − P eq ∝ Lγ 2 for Re ≪ 1, where P , L, andγ are the pressure, system size, and shear rate, and P eq is the pressure in the limitγ → 0. For the shear viscosity, the corresponding behavior is given by where η is the viscosity and η eq is the viscosity in the limitγ → 0. After these results had been obtained by the fluctuating hydrodynamics or kinetic theory, numerous MD simulations attempted to verify them. The basic idea was to probe the shear-induced LRCs by observing Eq. (2) or (3). The results remain controversial. Earlier simulation results [32][33][34] were interpreted in favor of the nonanalytical shear-rate dependence. In particular, Evans and coworkers [35][36][37][38] calculated the shear viscosity at the Lennard-Jones triple point and observed Eq. (3). However, more sophisticated simulations [41][42][43][44][45] support the assertion that the shear-induced correction behaves asγ 2 , not asγ 3/2 . Furthermore, the MD simulations of Sadus and coworkers [46,47] found that the exponent of pressure varies continuously between 1.2 and 2.0 depending on the temperature and density. More recently, Ortiz de Zárate et al. [15] reported that the shearinduced correction has two different origins, from shortand long-range scales. The long-range correction comes from the nonlinear coupling of the LRCs, which is calculated by the fluctuating hydrodynamics. The short-range correction is a molecular-scale effect and is independent of the LRCs. Ortiz de Zárate et al. estimated the magnitude of the short-range correction using kinetic theory and demonstrated that it yields non-negligible contributions. Their argument suggests the possibility that previous MD simulations captured the short-range correction. Thus, we find it difficult to extract the shear-induced LRCs from the shear-rate dependence of pressure P and shear viscosity η.
Another direction for probing the existence of LRCs is through direct observations, such as Eq. (1). Two groups studied the LRCs along this direction: Otsuki and Hayakawa [48,49] and Varghese et al. [44,45]. Otsuki and Hayakawa initially succeeded in observing the power-law decay of density and velocity fluctuations in a granular particle system, and found that the exponent α is close to the value predicted by the fluctuating hydrodynamics [48]. Their simulation size was insufficient for quantitatively examination of large-distance correlations beyond 10σ, where σ is the diameter of the particles. Subsequently, Varghese et al. performed a mesoscale simulation based on the multiparticle collision dynamics [45]. They successfully observed the shear-induced LRCs, and reported the behavior that is quantitatively consistent with the fluctuating hydrodynamics. However, the multiparticle collision dynamics is not based on microscopic interactions and cannot describe the molecular-scale behavior.
In this paper, we directly observe the LRCs by comparing the MD results with the linearized fluctuating hydrodynamics (LFH). We find that the MD result has large finite-size effects, and it prevents the occurrence of LRCs in small systems. We examine the finite-size effect using a sufficiently large system consisting of more than ten million particles, and show the existence of shear-induced LRCs without ambiguity.
Furthermore, we verify that our MD result is quantitatively consistent with the LFH solution for the large system. However, as we reduce the system size or increase the shear rate, the MD result gradually deviates from the LFH solution in the long-wavelength region. As a quantitative description of how the deviation increases, we define the characteristic wavenumber k vio such that the prediction from the fluctuating hydrodynamics is valid for k > k vio . We find that k vio has a nontrivial scaling dependence on the system size and shear rate.
The remainder of this paper is organized as follows. In Sec. II, we briefly review the analysis results based on the fluctuating hydrodynamics. In Sec. III, we explain the setup of the MD simulations. The main part of this paper is Sec. IV, where the MD result is presented and compared with the LFH solution. Section V gives our concluding remarks and discussions.

II. HYDRODYNAMIC DESCRIPTION OF SHEAR-INDUCED LONG-RANGE CORRELATIONS
The fluctuating hydrodynamics provides a powerful analytical tool for describing the nonequilibrium LRCs. Here, we briefly review the established results regarding shear-induced LRCs.

A. Model
We consider an isothermal fluid with a uniform temperature T defined in a three-dimensional region The isothermal fluid is described by two fluctuating fields, namely the density ρ(r, t) and the velocity v(r, t). The time evolution of ρ(r, t) and v(r, t) is given by [54] ∂ρ ∂t where Π ij (r, t) is the momentum flux tensor, written as Here, η 0 is the bare shear viscosity, ζ 0 is the bare bulk viscosity, p(r, t) is the pressure, and s ij (r, t) is the Gaussian random noise tensor satisfying We study the nonequilibrium steady state characterized by the average density field and velocity field. i.e., A schematic illustration of the steady state is presented in Fig. 1. In analyzing the fluctuating hydrodynamics, we focus on the bulk region and neglect boundary effects. Therefore, we do not have to specify the boundary condition. This is crucially different from the setup adopted in the MD simulations. Boundary effects are inevitable in the MD simulations because the nonequilibrium steady state is maintained using the Lees-Edwards boundary condition as explained in the next section. We will revisit this difference in Sec. IV, where we compare the MD result with the LFH solution.

B. Spatial correlation of momentum field
In the fluctuating hydrodynamics, the spatial correlation of the momentum is defined as where δg i (r, t) is given by In the MD simulations, we define a counterpart of this correlation function in terms of phase-space variables. To avoid confusion, we introduce the superscript FH to denote the fluctuating hydrodynamics. The existence of nonequilibrium LRCs is identified by the power-law decay of the correlation function.
The steady state under shear flow has translational symmetry [55]. This is expressed in terms of the correlation function as C FH ij (r, r ′ ) = C FH ij (r +a, r ′ +a), where a is an arbitrary constant vector. It is useful to introduce the Fourier transform of the correlation function We restrict our interest to the two correlation functions C FH xx (k) and C FH zz (k) at k y = k z = 0, and denote these as C FH xx (k x ) and C FH zz (k x ). Note that linear approximations of C FH yy (k x ) are not affected by the shear flow [26,45]. Therefore, we do not discuss C FH yy (k x ) in this paper. From Eqs. (4) and (5), we derive the integral expressions for C FH xx (k x ) and C FH zz (k x ) under the linear approximations: with We call Eqs. (13) and (14) the LFH solution. C FH xx (k x ) and C FH zz (k x ) correspond to the longitudinal and transverse momentum fluctuations, respectively. Therefore, Eq. (14) for C FH zz (k x ) does not contain ζ 0 , and is the same as that for an incompressible fluid. In contrast, C FH xx (k x ) is strongly affected by the compressibility of the fluid. These expressions were initially derived in Ref. [26]. Appendix A provides a brief sketch of the derivation; for further details, see Ref. [48].
From the LFH solution of Eqs. (13) and (14), we can see the existence of the shear-induced LRCs. First, aṡ γ → 0, Eqs. (13) and (14) reduce to This means that the correlation in the real space is given by the delta function. The correlation length is interpreted to be of the molecular scale. Forγ > 0, Eqs. (13) and (14) have nonequilibrium corrections, which lead to the LRCs. The asymptotic expression of Eq. (14) in the long-wavelength region is calculated as for k x ≫ k cross x , and for k x ≪ k cross x . Here, k cross x determines the crossover scale between Eqs. (17) and (18), and is given by These expressions imply that an additional correlation proportional to k −4 x appears at short-distance scales and crosses over to k −4/3 x at large-distance scales. Such power-law behavior in the Fourier space corresponds to an algebraic decay in the real space. We can repeat the same discussion for C FH xx (k x ) and derive the LRC [26]. We use the LFH solution of Eqs. (13) and (14) to probe the existence of shear-induced LRCs in the MD simulation. Additionally, we quantitatively examine the validity of the LFH solution. Note that expressions such as Eqs. (13) and (14) provide a starting point for explaining various phenomena coming from shear-induced LRCs. For example, Lutsko and Dufty [53] derived a nonequilibrium correction to the shear viscosity in the form of Eq. (3). Similarly, Wada and Sasa [21] and Ortiz de Zárate et al. [15] derived the shear-rate dependence of pressure for incompressible fluids. Therefore, it is important to establish the LFH solution quantitatively from the molecular-scale dynamics.

A. Model
We consider an N particle system that is con- The dynamics is given by where (r i , p i ) is the position and momentum of the ith particle, m is the mass, U (r) is the interparticle interaction, and f th i is the force acting on the ith particle from a thermostat. We use the Weeks-Chandler-Andersen (WCA) potential as the interparticle interaction, which is the Lennard-Jones potential with the cutoff-length r LJ c = 2 1/6 σ, i.e., where θ(r) is the Heaviside step function and σ is the diameter of the particle.
To maintain a constant temperature under the shear flow, we use the dissipative particle dynamics (DPD) thermostat [56], which is given by Here,r ij is a unit vector in the direction and . γ and T represent the friction and the temperature of thermostat, respectively. Because the DPD thermostat satisfies the fluctuation-dissipation relation, our model relaxes to equilibrium when no external forces are imposed. The cutoff length r DPD c is set to 2.0σ. Note that the DPD thermostat obeys Newton's third law, which ensures momentum conservation [57]. This is why we apply the DPD thermostat. Indeed, one of the origins of nonequilibrium LRCs is the conservation law [5,6].
The shear flow is realized using the Lees-Edwards boundary condition [58,59] along the z-axis. Along the x-and y-axes, we impose standard periodic boundary conditions. Thus, the velocity profile in the steady state is realized asṽ Note that the Lees-Edwards boundary condition violates the momentum conservation law along the x-direction.
The total amount of momentum along the x-direction depends on the number of atoms that leave the lower and upper boundaries. This is because, when one atom leaves the lower (upper) boundary z = −L z /2 (z = L z /2) with velocity u, the corresponding atom is introduced from the upper (lower) boundary with velocity u ±γL zêx . However, except at the boundaries, the local conservation law still holds. Moreover, in the steady state, because the net mass transfer via the lower or upper boundary is balanced, the violation is sufficiently small and the time-averaged total momentum must be zero. Thus, we expect that the effect of the violating the conservation law through the Lees-Edwards boundary condition will be sufficiently small [60].

B. Parameters
In the numerical simulations, all quantities are measured by the Lennard-Jones units (m, σ, ǫ). In particular, the time is measured by τ unit = mσ 2 /ǫ. All the MD simulations are performed by LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [61,62]. The time integration is calculated by the velocity Verlet algorithm. The timestep is set to 0.0025, 0.00375, or 0.005 depending on the shear rate and the system size. We fix the temperature and friction of the thermostat to T = 1.0 and γ = 1.0, respectively. The density ρ 0 ≡ N/L x L y L z is fixed to 0.78.
The transport coefficients take almost the same value in all simulations. In particular, we use η 0 = 1.74 and ζ 0 = 14.04 to compare the MD result with the LFH solution, which are calculated from the Green-Kubo formula (see Appendix B for details).

C. Observation method
All observations are performed in the nonequilibrium steady state, which is prepared by different methods depending on the system size L z . For L z ≤ 512, we start from the initial state in which the particles are randomly located with zero overlaps. We then perform the relaxation run for about 10 times the relaxation time. The relaxation time is estimated from the relaxation of the velocity profile (see Appendix C for details). For L x = 1024, L y = 32, and L z = 512 (N = 13 086 228), the relaxation time is about 3000, and the relaxation run with a timestep 0.0025 takes 100 hours using 16 nodes of the ISSP supercomputer (AMD EPYC 7702, 64 cores × 2 per node).
For L z > 512, we adopt a locally relaxed state as the initial state and perform the relaxation run for about three times the relaxation time. For example, for L x = 1024, L y = 32, and L z = 1024, the initial state is prepared by combining two different relaxation states for L x = 1024σ, L y = 32σ, and L z = 512σ.
After the relaxation run, we observe the correlation function of the momentum fluctuation: where · represents the time average in the steady state and the ensemble average over different noise realizations.
Here, δĝ i (k) (i = x, y, z) is the Fourier transform of the momentum density field with the mean flow subtracted, which is expressed as with We rewrite δg(r) in terms of the microscopic density field ρ(r) = N i=1 mδ(r − r i ) and momentum fieldg(r) = N i=1 p i δ(r − r i ) as δg(r) =g(r) −ρ(r) ṽ(r) .
By comparing Eq. (29) to Eq. (11), we find that C MD ij (k;γ) is the microscopic counterpart of C FH ij (k;γ). We also introduce the relative deviation ∆ ij (k x ) to verify the validity of the LFH solution: In a region where the relative deviation is large, the LFH solution cannot be applied to describe the MD result. For a quantitative discussion, we introduce the criterion ∆ ij (k x ) < 0.1 for the applicability of the LFH solution.
We then define the characteristic wavenumber k vio x as the largest wave number satisfying ∆ ij (k x ) > 0.1. For the wavenumber region k x > k vio x , the descriptions given by the fluctuating hydrodynamics are quantitatively valid.

IV. MAIN RESULTS
A. Nonequilibrium LRC Figure 2 presents the results forγ = 0.02, L x = 1024, L y = 512, and L z = 1024 (N = 26 172 456), which is the largest system size that we examined. The blue lines in the left-and middle-hand panels represent the equilibrium value from Eq. (16), and the deviations from this value give the shear-induced correction. The black lines show the LFH solutions of Eqs. (13) and (14). The MD result clearly exhibits shear-induced corrections, and is in quantitative agreement with the LFH solution except in the long-wavelength region.
We can identify the nonequilibrium LRC from the power-law behavior of Eqs. (17) and (18), as explained in Sec. II B. The right-hand panel of Fig. 2 shows double-log plots of C MD zz (k x ) and C FH zz (k x ) (red and black lines, respectively). From the LFH solution, the crossover scale between the k −4 x and k −4/3 x behavior is k cross x = 0.086. The MD result is quantitatively consistent with the LFH solution around the crossover region. Thus, we conclude that the MD result exhibits the nonequilibrium LRC.
In the long-wavelength region, there is a qualitative difference between C MD xx (k x ) and C FH xx (k x ). Specifically, as shown in the left-hand panel of Fig. 2, the MD result monotonically increases from the equilibrium value as k x → 0, whereas the LFH solution monotonically decreases from the equilibrium value. We expect that the boundary effect has a strong influence on the longwavelength behavior, and thus causes this difference. The boundary effect is neglected in the LFH solution, as explained in Sec. II A. We now study how the MD result is affected by changes in the shear rate and system size.

B. System-size dependence of nonequilibrium LRC
We first examine the MD result for the various system sizes. Figure 3 shows that the MD result does not depend on L y . In contrast, we can observe strong L zdependence in Fig. 4, where the MD result approaches the LFH solution as L z increases. In Figs. 4-(a) and -(d), we present the L z -dependence of C MD xx (k x ) and C MD zz (k x ). Figure 4-(d) shows that the nonequilibrium LRC of C MD zz (k x ) gradually grows from the equilibrium value as L z increases. In contrast, Fig. 4-(a) shows that the nonequilibrium correction of C MD xx (k x ) is positive for the small system sizes of L z = 32, 64, and 128. This behavior is inconsistent with the LFH solution; the correction of C FH xx (k x ), which is the second term of Eq. (13), is always negative. As L z increases, the correction dips into the negative region and the positive correction region becomes smaller. We can infer that the positive correlations for smaller L z mainly come from finite-size and boundary effects.
We now examine how the MD result approaches the LFH solution as L z increases. The deviations ∆ xx (k x ) and ∆ zz (k x ) are plotted in Figs. 4-(b) and -(e). Moreover, the characteristic wavenumber k vio for C xx (k x ), and for C zz (k x ). These are depicted by the blue lines in Figs. 4-(c) and -(f). Note that Eqs. (31) and (32) are the quantitative relations and enable us to estimate the finite-size effects. Furthermore, we consider the dependence of the scaling relation on the shear rateγ. We plot k vio x as a function of L z for severalγ in Fig. 5. The figure shows that the scaling form k vio x ∝ L −ω z holds, regardless of the value oḟ γ. For C zz (k x ), ω takes values of 0.34-0.45 depending onγ. In contrast, for C xx (k x ), ω is close to 0.27, and is largely insensitive toγ.

C. Shear-rate dependence of nonequilibrium LRC
We now examine the shear-rate dependence of the LRCs for a fixed system size. In Fig. 6, we plot the correlation C zz (k x ) and the deviation ∆ zz (k x ) for various values ofγ from 0.005 to 0.06. We observe that the deviation increases monotonically asγ increases from 0.01 to 0.06 in Fig. 6-(g). However, the result forγ = 0.005 does not exhibit this tendency. We can infer that the LRC does not fully develop whenγ = 0.005 because L z is too small.
The inset of Fig. 6-(g) shows k vio x as a function ofγ. Clearly, k vio x is linearly dependent onγ fromγ = 0.01 tȯ γ = 0.06. By fitting this with the functional form Aγ +B, we obtain the quantitative relation Similar behavior can be observed for C xx (k x ),

D. Kinetic temperature and pressure
Finally, we study the kinetic temperature T ke , which is defined as and the pressure P z along the z-direction, which is defined as Here, f ij,z is the z-component of the intermolecular force between particles i and j. These quantities are plotted in Fig. 7. Previous MD simulations [32-43, 46, 47, 50, 51] have explored theγ-dependence of the kinetic temperature and pressure to probe the LRC, as explained in the Introduction. Following these studies, we fit the simulation data to the form Aγ B + C and obtain P z = 4.8985γ 1.9831 + 6.0638. (38) Both exponents are close to 2. We now consider whether the long-or short-range contributions dominate our result, as suggested in Ref. [15]. To this end, we decompose the corrections into the contributions from the short-and the long-range scales: Our result suggests the dominance of the short-range scale as follows. The correction proportional toγ 2 in the fluctuating hydrodynamics is linearly dependent on the system size L: which comes from the k −4 x -behavior of Eq. (17). However, our result is almost independent of L z , as shown in the inset of Fig. 7, although our result catches the k −4 xtail. Thus, our result supports the assertion that theγ 2dependence comes from the short-range scale instead of the nonequilibrium LRC. However, recall that the LFH solution is not valid in the long-wavelength region, as shown in Fig. 4. Therefore, further theoretical studies on the short-range corrections are required to form a final conclusion.

V. CONCLUDING REMARKS AND DISCUSSION
Let us compare our result to those reporeted by Otsuki and Hayakawa [48] and Varghese et al. [45]. These previous studies directly observed the shear-induced LRC in particle-based simulations. Table I presents a comparison of their setups with that of our simulations. The system size in our study is about 10 times larger than that in the previous studies. As a result, we can systematically study the finite-size effect on shear-induced LRCs. We showed that the LFH solution is quantitatively consistent with the MD result when L z is sufficiently large. Conversely, for smaller L z , the MD result deviates from the LFH solution in the long-wavelength region. Such deviations were also observed in the previous studies [45,48]. Varghese et al. proposed that these deviations originated from the density-dependence of viscosity. However, our simulations have clarified that the derivations are caused by an insufficient system size.
Furthermore, we examined how the deviations depend on the system size and the shear rate. As a quantitative examination, we introduced the characteristic wavenumber k vio x associated with the breakdown of the hydrodynamic description. k vio x determines the applicable wavenumber region of the LFH solution as k x > k vio x . We then found two scaling relations, k vio x ∝ L −ω z at fixeḋ γ and k vio x ∝γ, for a fixed L z . The interesting point is that the finite-size effect is nonnegligible in a large region. For example, k vio x is obtained from Fig. 2 as k vio x ≃ 0.0237 for L z = 1024 andγ = 0.02. In the real space, the corresponding x vio ≡ 2π/k vio x is about 265. Therefore, if we consider the system with L x = L y = L z = 1024, the LFH solution breaks down in about three-quarters of the region, 0.26L x < x < L x , where large finite-size effects exist. Note that the magnitude of the finite-size effects is related to the value of the exponent ω. The scaling relation k vio x ∝ L −ω z can be rewritten as x vio ∝ L ω z . By noting that the breakdown of the LFH solution occurs in the region L ω z < x < L x , we find that a smaller ω yields finite-size effects in a larger region. Actually, ω = 0.27-0.45 as our model is quite small.
The question to be asked is the origin of such large finite-size effects. We can infer that they come from the Lees-Edwards boundary condition and the nonlinearity of the fluctuating hydrodynamics. Future work should analyze these effects in the fluctuating hydrodynamics. As a related problem, it is interesting how the hydrodynamic description predicts the exponent ω.
Finally, we remark on the utility of the quantitative relations for k vio x , such as Eqs. (31)- (34). They enable us to estimate the finite-size effects in larger-size simulations from smaller-size simulations. For example, we can use the estimation to observe the k −4/3 x -tail of C zz (k x ). We could not observe this tail as shown in the right-hand panel of Fig. 2 because the hydrodynamic description breaks down before C zz (k x ) exhibits the k Appendix A: Brief sketch of derivation of Eqs. (13) and (14) To derive the integral expressions in Eqs. (13) and (14) for C xx (k x ) and C zz (k x ), we use two approximations. Here, we briefly sketch their derivation while focusing on these approximations.
The first approximation is to neglect the nonlinear fluctuations. ρ(r, t), p(r, t), and v(r, t) are expanded around the zero-order solution as ρ(r, t) = ρ 0 + δρ(r, t), where c T is the isothermal speed of sound. By substituting Eq. (A1) into Eqs. (4) and (5) and neglecting the higher-order terms of δρ and δv, we have Under this approximation, the momentum correlations are connected with the velocity correlations as The second approximation is used when we decompose the longitudinal and transverse waves. Here, it is convenient to introduce the oblique coordinate used by Lutsko and Dufty [26]: The vectors {ê (a) (k)} a=1,2,3 are a set of orthogonal unit vectors given byê (1) where k = |k|, k ⊥ = k 2 − k 2 z , andk ⊥ = k 2 − k 2 z /k. The time evolution ofξ(k, t) is immediately obtained by substituting the inverse transformation of Eq. (A4) into Eq. (A2) as where α = 1, 2, 3 and the matrices B, C, and D(k) are given by To decompose the longitudinal and transverse waves in Eq. (A8), we need to solve the eigenvalue problem: whereζ (i) (k) and λ i (k), respectively, are the ith eigenvector and eigenvalue (i = 1, 2, 3, 4). At zero shear rate, the eigenvalue problem in Eq. (A15) reduces to the diagonalization of the matrix −ikB +k 2 C, which can be solved exactly. However, for a finite shear rate, the longitudinal and transverse waves obtained at zero shear rate are strongly coupled and, as a result, the eigenvalue problem in Eq. (A15) is difficult to solve exactly. Therefore, we use the perturbation expansion with respect to the wave vector k: and calculate the solution to O(k 2 ). This approximation is the second one that is used to obtain Eqs. (13) and (14). The calculation of the perturbation expansion is lengthy but straightforward, and so the detailed steps are omitted in this paper. The solution of Eq. (A8) is written using the eigenvector as and a (i) (k, t) is given as the solution of the following equation We can solve Eq. (A19) without any approximations to obtain the integral expression for a (i) (k, t) a (i) (k, t) = a (i) (k(−t), 0) exp − where k(−t) = (k x , k y , k z +γtk x ).

Appendix B: Measurement of viscosity
There are two viscosities, η 0 and ζ 0 , in the fluctuating hydrodynamic equation. We use the Green-Kubo formula to measure them in the MD simulations [63,64]. The Green-Kubo formula provides the microscopic expression of the transport coefficient, and is a useful tool for computing the transport coefficient in the MD simulations.
The Green-Kubo formula for the viscosity is given by [65]: whereP αβ (t) is the microscopic expression of total stress tensor and δP (t) = 1 3 αP αα (t)− 1 3 αP αα (t) . The summation in the expression for η 0 is taken over the off-diagonal elements.