Nonequilibrium dynamics of coupled oscillators under the shear-velocity boundary condition
Hidetsugu Sakaguchi



Interdisciplinary Graduate School of Engineering Sciences, Kyushu University,  Kasuga, Fukuoka 816-8580, Japan
§ ABSTRACT 
Deterministic and stochastic coupled oscillators with inertia are studied on the rectangular lattice under the shear-velocity boundary condition. Our coupled oscillator model exhibits various nontrivial phenomena and there are various relationships with wide research areas such as the coupled limit-cycle oscillators, the dislocation theory, a block-spring model of earthquakes, and the nonequilibrium molecular dynamics. We show numerically several unique nonequilibrium properties of the coupled oscillators.  We find that the spatial profiles of the average value and variance of the velocity become non-uniform when the dissipation rate is large. The probability distribution of the velocity sometimes deviates from the Gaussian distribution. The time evolution of kinetic energy becomes intermittent when the shear rate is small and the temperature is small but not zero. The intermittent jumps of the kinetic energy cause a long tail in the velocity distribution.
Nonequilibrium dynamics of coupled oscillators under the shear-velocity boundary condition
    Hidetsugu Sakaguchi
    September 9, 2024
==========================================================================================



§ INTRODUCTION

Coupled phase oscillators called the Kuramoto model have been intensively studied by many authors as a simple solvable model of collective synchronization <cit.>.  There are many generalized models for the Kuramoto model. One model includes the inertia term or the second derivative of the phase variable <cit.>. The coupled phase oscillators on square or cubic lattices was called oscillator lattices and the phase transition via the collective synchronization was studied in the finite-dimensional systems <cit.>.
 
In the previous paper, we demonstrated that the vortex motion in the oscillator lattices is closely related to the dislocation motion in solids <cit.>. The dislocation motion is important to understand the mechanical properties of the solid under the shear stress <cit.>. 
The correspondence between the coupled phase oscillators and lattice dynamics is as follows. The phase variable in the oscillator lattice with inertia is interpreted as the one-dimensional displacement in the z direction in a three-dimensional crystal, and the displacement is assumed to be uniform in the z direction. That is,  the displacement z_i,j,k at the (i,j,k) site takes the same value for all k's, and is expressed as z_i,j. The vortex in the oscillator lattice corresponds to the screw dislocation in the crystal. The shear stress on the boundaries corresponds to the external force at the boundaries of the oscillator lattices. We found that the vortex begins to move if the shear stress is beyond a critical force, which corresponds to the Peierls force in the dislocation theory. The spontaneous generation of vortices and the complex motion with nonzero frequencies correspond to the slip motion or the plastic flow in solids. We further studied the oscillator lattices under external noises, and found that the vortex motion occurs under the critical force owing to the fluctuation effect, which corresponds to the dislocation motion at a finite temperature <cit.>. 

The coupled oscillators under the shear boundary conditions are also related to the mechanics of earthquakes and faulting <cit.>. The block-spring models of earthquakes have been studied by many authors <cit.>. In the block-spring model, each block has a mass and is coupled with the neighboring blocks through a linear spring. The shear force is applied from the top plate moving with a small velocity. A large slip motion in the block-spring model corresponds to an earthquake. 
 

The deterministic chaos and fluctuation in the oscillator lattice under the shear stress is related to the nonequilibrium statistical mechanics. The nonequilibrium statistical mechanics has been numerically studied with the molecular dynamics by many authors <cit.>. Various nonequilibrium properties such as the shear viscosity <cit.> and heat conductivity <cit.> have been investigated by numerical simulations of Newton's equation of motion of many particles. The relationship between the reversible equations and irreversible behavior was discussed by several authors. <cit.>. 
In the molecular dynamics of fluids, interacting pairs of molecules change with time. In our oscillator lattice models, each oscillator interacts only with the nearest neighbors, but a large displacement occurs by the phase slip between the neighboring oscillators. In this paper, we will discuss several nonequilibrium properties and velocity distributions in the coupled oscillators under the shear velocity boundary conditions. 



§ DETERMINISTIC COUPLED OSCILLATORS UNDER THE SHEAR-VELOCITY BOUNDARY CONDITION

In this section, we consider deterministic coupled oscillators with inertia on the rectangular lattice:  

    d^2 z_i,j/dt^2=K_x∑_i^'=i-1,i+1(z_i^',j-z_i,j)+K_y∑_j^'=j-1,j+1sin (z_i,j^'-z_i,j)-ddz_i,j/dt,

where z is the displacement in the z direction  at the lattice point of (i,j),  d is the coefficient of resistance force in proportion to the velocity, and (i^',j^')'s are the nearest neighbor sites of the (i,j) site on the rectangular lattice of L_x× L_y. The linear and sinusoidal couplings are assumed respectively in the x- and y-directions. The coupling strength K_x and K_y are set to be 1 in this paper. 
Periodic boundary conditions are imposed at i=1 and i=L_x. In the y direction, the boundary conditions corresponding to the shear velocity are imposed, that is, z_i,j=0 at j=0 and z_i,j=v_0 t at j=L_y+1. 
The numerical simulation was performed with the Runge-Kutta method with timestep Δ t=0.005. For d=0 and v_0=0, the coupled oscillator system has the Hamiltonian:

    H=∑_i,j1/2 (dz_i,j/dt )^2+1/2∑_i,j∑_i^'=i-1,i+1(z_i^',j-z_i,j)^2+∑_i,j∑_j^'=j-1,j+1{1-cos(z_i,j^'-z_i,j)}.


For L_x=L_y=1, a single oscillator obeys the equation

    dz/dt   =    v,
    dv/dt   =   sin(v_0t-z)-sin z-dv.

Equation (3) is a conservative (dissipative) system at d=0 (d 0) since the phase space volume decays as e^-dt. Figure 1(a) shows the Poincaré map (stroboscopic mapping) in the phase space of ( mod(z,2π),v) at t=2π n/v_0 (n: integer) for v_0=0.1 and d=0. The initial condition is z(0)=0 and v(0)=0.05001. A chaotic dynamics is observed owing to the periodic forcing term sin(v_0t-z). Figure 1(b) shows the velocity distribution of v. The average velocity  is 0.05. The velocity is confined between -2.78 and 2.88. 




For L_x=300 and L_y=1, the one-dimensional coupled oscillators obey 

    dz_i/dt   =    v_i,
    dv_i/dt   =   K{sin(v_0t-z_i)-sin z_i}-dv_i+K(z_i+1-2z_i+z_i-1).

Figure 2(a) shows the time evolution of the average value of the kinetic energy: E=(1/L_x)∑_i=1^L_x(1/2)v_i^2 at v_0=0.01, 0.1, and 1 for d=0. The kinetic energy increases with time. The time evolution can be approximated as a power law E∝ t^α with α≃ 0.27, 0.35, and 0.41 respectively for v_0=0.01, 0.1, and 1. The increase of kinetic energy or heating occurs under the shear-velocity boundary conditions. The heating is not observed at L_x=1 as shown in Fig. 1. We have checked that the heating occurs for L_x≥ 4 at v_0=0.1. The critical size for the heating depends on v_0, that is, the critical size L_xc was evaluated as 3, 8, and 13 respectively for v_0=1, 0.05, and 0.03. The critical size decreases as v_0 is larger. 
Although the heating occurs due to the nonequilibrium boundary conditions, the mechanism is not clear. Figure 2(b) shows the probability distributions of v at v_0=0.1 near t=50000 and 100000, The probability distributions are well approximated at the Gaussian distribution, and the variance increases with time. The Gaussian distribution of the velocity corresponds to the Maxwell distribution of the ideal gas in the statistical mechanics. 

Figure 3(a) shows the time evolution of E at v_0=0.1 and d=0.001 and 0.05. The increase of the kinetic energy stops owing to the dissipation term -dv_i  with d=0.001. 
The kinetic energy is fluctuating at d=0.001 and 0.05. The average value of E and the fluctuation amplitude is the same order at d=0.05.  
Figure 3(b) shows the relationship between d and the average kinetic energy ⟨ E⟩ at v_0=0.1. The dashed line is ⟨ E⟩∝ t^-α with α=0.35. The numerical result of the time evolution of E as E∝ t^α at d=0 suggests that E satisfies an approximate equation dE/dt=β E^1-1/α at d=0.   
By adding the dissipation term of -2dE, the time evolution of the kinetic energy E is assumed to be

    dE/dt=β E^1-1/α-2dE.

The stationary value of E is expressed with 

    E∝ d^-α.
 
Figure 3(b) shows that the power law ⟨ E⟩∝ d^-α is satisfied for d<0.001. ⟨ E⟩ deviates from the power law for d>0.001. 
Figure 3(c) shows the probability distributions P(v) at d=0.001 and 0.05. 
The dashed lines are P(v)=1/√(2πσ^2)exp{-(v-v_0/2)^2/(2σ^2)} with σ^2=2.35 and P(v)=1/(2γ)exp(-|v-v_0/2|/γ) with γ=0.39. (The dashed line  P(v)=1/√(2πσ^2)exp{-(v-v_0/2)^2/(2σ^2)} is hardly seen because of the overlapping with the numerically obtained P(v).) For small d, P(v) is well approximated at the Gaussian distribution, however, P(v) deviates from the Gaussian distribution and closer to the exponential distribution. Figure 3(d) shows the numerically obtained kurtosis q=⟨(v-v_0)^4⟩/(⟨ (v-v_0)^2⟩)^2-3 as a function of d. The kurtosis is 0 for the Gaussian distribution and q is a quantity expressing the deviation from the Gaussian distribution. Figure 3(d) shows that the deviation from the Gaussian distribution appears for d>0.002. 
The deviation of P(v) from the Gaussian distribution at d=0.05 might be due to the large fluctuation of the kinetic energy compared to the average value, although the mechanism is not well understood. 




Next, we show numerical results for a rectangular system of L_x× L_y=300× 50.
Figure 4(a) shows the time evolutions of the variance of v_j for i=1,2,⋯,L_x at j=5,30, and 45 for d=0 and v_0=0.5. The variance increases with time and do not depend on j. The dashed line denotes the power law of t^α with 0.35. The variance increases faster than a power law, but it is confirmed that the heating occurs also in this rectangular system. Figure 4(b) shows the profile of the average velocity ⟨ v_j⟩ at d=0 and v_0=0.5. The average velocity changes linearly in j. The linear shear flow is characteristic in the Newton fluid with the linear viscosity. 

Figure 5(a) shows the profiles of the average velocity ⟨ v_j⟩ at d=0.0005, 0.002, 0.005, and 0.01 for v_0=0.5. The average velocity changes almost linearly in j at d=0.005 and the difference ⟨ v_j+1⟩-⟨ v_j⟩ or the shear rate between the jth and (j+1)th layers is nearly constant.  However, the shear rate increases with j at d=0.002. The shear rate takes a small value for j<33 and takes a large vale for j>33 at d=0.005. The average velocity changes rapidly from 0 to v_0=0.5 near j=L_y at d=0.01. That is, the shear rate is almost zero for j<46 and takes a very large value for j>46 at d=0.01.  The localization of shear rate  near the boundary j=L_y is a phenomenon similar to the one observed in the coupled noisy oscillators without inertia found in the previous paper <cit.>. We showed that the shear localization occurs when v_0 is relatively large, and the vortex density is  high in the boundary region where the shear rate of the average velocity takes a large value and the vortex density takes a lower constant value in the bulk region where the average velocity is almost zero. The phenomena of shear or strain  localization are observed in various complex fluids and plastic materials such as the plug flow in the Bingham fluid <cit.>, the shear banding <cit.>, and the formation of shear zones in rocks <cit.>, although the mechanism might be different. The reason of the shear localization in our model is not well understood yet. 
Figure 5(b) shows the profiles of the variance of the velocity v_j at d=0.0005, 0.002, 0.005, and 0.01. The variance of the velocity corresponds to the temperature. The temperature is almost uniform for d=0.0005, but the temperature is higher near j=L_y for d=0.01. This is the excess heating due to the shear of the velocity profile near j=L_y. Even if the average velocity is zero, the variance of the velocity takes a nonzero value and changes slowly for j<46 at d=0.01, that is, the oscillators are not stationary. 

Figure 6(a) shows the time evolution of the kinetic energy per each oscillator E=(1/L_xL_y)∑ (1/2)v_i,j^2  at d=0.01. Intermittently, E takes a large value, which is related to the excess heating. Figure 6(b) shows the probability distribution of the velocity v_j at j=20 and j=48 for d=0.01. The velocity distribution at j=20 is well approximated at the Gaussian distribution exp{-v^2/(2σ^2)}/√(2πσ^2) with σ^2=0.22. Figure 6(c) is the velocity distribution P(v) in the semi-logarithmic scale at j=48 for d=0.01. The dashed line is 1/√(2πσ^2)exp{-(v-0.35)^2/(2σ^2)} with σ^2=0.785. The velocity distribution at j=48 is not mirror symmetric around the peak value and deviated from the Gaussian distribution probably due to the shear velocity. 


§ COUPLED NOISY OSCILLATORS UNDER THE SHEAR-VELOCITY BOUNDARY CONDITION

The coupled oscillator model at a finite temperature is assumed to be  

    d^2 z_i,j/dt^2=K_x∑_i^'=i+1,i-1(z_i^',j-z_i,j)+K_y∑_j^'=j+1,j-1sin (z_i,j^'-z_i,j)-ddz_i,j/dt+ξ_i,j(t),

where ξ_i,j(t) is the Gaussian white noise satisfying ⟨ξ_i,j(t)ξ_i^',j^'(t^')⟩=2dTδ_i,i^'δ_j,j^'δ(t-t^'). Here, T is a quantity corresponding to the temperature. The coupling constants are assumed to be K_x=K_y=1. If the periodic boundary conditions are imposed for all the four boundaries: i=1, i=L_x, j=1, and j=L_y, the thermal equilibrium distribution is realized. The probability distribution of the velocity v_i,j is the Gaussian distribution

    P∝ e^-H/T∝ e^-∑_i,jv_i,j^2/(2T),

where H is the total energy of Eq. (2). 
The shear-velocity boundary conditions induce a nonequilibrium state.
Figure 7(a) shows the average velocity ⟨ v_j⟩ as a function of j at T=2, d=0.01 and v_0=0.5. By the effect of finite temperature, the velocity profile is more delocalized than the one shown in Fig. 5(a). The dashed blue line is an approximate line 0.5exp{-0.185(50-j)} for the velocity profile. Figure 7(b) shows the variance profile of the velocity v_j. 
The variance is almost T=2 for j<35, which is almost equal to the strength of the noise, however, the variance is slightly larger than 2 near the boundary j=51. Figure 7(c) is the probability distributions at j=20 and 48. The probability distribution can be well approximated at the Gaussian distribution in this finite-temperature system. When T is large, the stochastic motion is dominant as shown in Fig. 7. 

The length scale of the shear localization is evaluated by a quantity λ=∑_j=1^L_y(L_y+1-j)⟨ v_j⟩ /∑_j=1^L_y⟨ v_j⟩.
Figure 8(a) shows λ for several d's at T=2 and v_0=0.5, and Fig. 8(b) shows λ for several v_0's at T=0.5 and d=0.01. The length scale λ decreases with d and v_0.






When v_0 and T are small, a more regular motion appears. Figure 9(a) shows the profile of the average velocity ⟨ v_j⟩ at d=0.01, T=0 and v_0=0.01. Figure 9(b) shows the time evolutions of z_i,j between t=390000 and 430000 at (i,j)=(n/2,1), (n/2,2), (n/2,3), and (n/2,4) for d=0.01, T=0 and v_0=0.01. z_i,j increases stepwise except for j=1. That is, the phase slip of 2nπ occurs almost simultaneously.  
The period of the stepwise increase is around 2500, and the phase slip of 4π occurs at j=2 and the phase slip of 8π occurs for j≥ 3. The average velocity is zero at j=1 and takes the same value for j≥ 3, and 
the velocity profile has a jump from 0 to 0.01 at j=2 as shown in Fig. 9(a).  Figure 9(c) shows the time evolution of the kinetic energy E for T=0, v_0=0.01, and d=0.01. The kinetic energy increases periodically due to the periodic phase slips, whose period is around 2500.  

Figure 10(a) shows the profile of the average velocity ⟨ v_j⟩ at d=0.01, T=0.05 and v_0=0.01. The average velocity changes continuously with j owing to the finite temperature effect and the profile might be approximated at the linear one. 
Figure 10(b) shows the time evolution of Δ z=√(∑_j∑_i(z_i,j-∑_i z_i,j/L_x)^2/(L_xL_y)) at d=0.01, T=0.05 and v_0=0.01. Here, ∑_i(z_i,j-∑_i z_i,j/L_x)^2/L_x is the variance ⟨ (Δ z_i,j)^2⟩ of the displacement z_i,j within the jth layer, and Δ z is the root mean square of the variance for j=1,2,⋯, L_y: Δ z=√(∑_j ⟨ (Δ z_i,j)^2⟩/L_y). Δ z is fluctuating by the noise effect, but sometimes jumps to a fairly large value and then decays stepwise. Figure 10(c) is the time evolution of the kinetic energy E. The kinetic energy E jumps when Δ z jumps. When Δ z and E jump, phase slips occur between some neighboring layers, which will be shown in Fig. 11(a). Figure 10(d) shows the probability distribution of v_j at j=20 in the semi-logarithmic scale. The dashed line is the Gaussian approximation 1/√(2πσ^2)exp{-(v-0.00345)^2/(2σ^2)} with σ^2=0.05025. The Gaussian approximation is good near v=0, however, the deviation from the Gaussian distribution is observed for large |v|. The acceleration of v_j at the jumps of E might be related with the velocity distribution with the long tail. 

Nontrivial spatio-temporal dynamics is observed when Δ z and E jumps.  Figure 11(a) shows the profiles of z_i,j as a function of j for i=25 at t=12000+25n (n=0,1,⋯,30), which corresponds to the second peak region of Δ z shown in Fig. 10(b). For example, the lowest line represents a profile of the displacement z_i,j at i=25 and t=12000, and the second lowest one represents z_i,j at i=25 and t=12025. There is a discontinuity at j=25 at the profiles of z_i,j for t≤ 12375, which was generated by the phase slip event near t=5900 which corresponds to the first peak in Fig. 10(b). Figure 11(a) shows that a phase slip occurs and the discontinuity of z is created at j=9 near t=12400. And then phase slips occur at j=34, and 38 near t=12800. The jump of Δ z near t=12000 shown in Fig. 10(b) is caused by these successive phase slips. 
The spatial variation Δ z of the displacement z_i,j decreases stepwise after the jump as shown in Fig. 10(b). The stepwise decrease of Δ z is understood through the spatio-temporal dynamics of z_i,j as follows.  
Figure 11(b) shows the time evolution of z_i,j at six points (i,j)=(1,36), (51,36), (101,36), (151,36), (201,36), and (300,36) in the same layer of j=36 from t=13000 to 15000. It is observed that z_i,j's tend to be aligned  with time. The stepwise merging of z_i,j makes the spatial variation of z_i,j decrease in time, which corresponds to the stepwise decrease of Δ z shown in Fig. 10(b). Figure 11(c) shows the time evolution of z_i,j at the same layer j=36. The lowest line is the profile z_i,j as a function of i at j=36 and t=14500, and the second lowest one shows the profile of z_i,j at j=36 and t=14525. There are vortex and antivortex solutions in Eq. (1) because of the sinusoidal coupling in the y direction. The kink and antikink structures seen in Fig. 11(c) represent the vortex or antivortex.  
Figure 11(c) expresses that  z_i,j tends to be uniform through the collisions of vortex and antivortex, which corresponds to the stepwise merging  of z_i,j shown in Fig. 11(b) and the stepwise decrease of Δ z shown in Fig. 10(b). 

The intermittent time evolution of the coupled oscillators at v_0=0.01 and T=0.05 as shown in Figs. 10 and 11  can be qualitatively explained as follows. The small perturbations δ z_i,j in Eq. (1) obeys      

    d^2 δ z_i,j/dt^2=K_x∑_i^'=i+1,i-1(δ z_i^',j-δ z_i,j)+K_y∑_j^'=j+1,j-1cos (z_i,j^'-z_i,j)(δ z_i,j^'-δ z_i,j)-ddδ z_i,j/dt.

If the difference z_i,j+1-z_i,j between the neighboring layers increases with time and is larger than π/2, cos(z_i,j+1-z_i,j) becomes negative and the perturbations grow exponentially owing to the linear instability, which induces the phase slips. If the phase slips occur, the kinetic energy jumps and the spatial variation Δ z in the x direction grows. After the phase slip, the spatial fluctuations tend to decay through the collisions of vortex and antivortex. After Δ z and E decay to the thermal fluctuation level at the finite temperature T=0.05, the differences z_i,j+1-z_i,j between the neighboring layers increase further again, and phase slips occur at points (i,j) where cos(z_i,j+1-z_i,j)<0, which are generally different from the previous phase-slip points. This process repeats many times. These intermittent time evolutions are observed at small values of v_0 and T. If the temperature or the shear rate v_0/L_y is large, the spatial irregularity owing to the temperature or deterministic chaos increases and the intermittent behavior of E becomes unclear. 


The three parameters v_0, T, and d determine the nonequilibrium dynamics in the coupled noisy oscillators Eq. (5). Although the whole parameter range is not yet investigated and the nonequilibrium dynamics is not theoretically understood, we investigate numerically the excess kinetic energy characterizing the nonequilibrium state by changing some control parameters. The kinetic energy jumps intermittently from the thermal equilibrium level T/2 as shown in Fig. 10(c). The excess kinetic energy is the temporal average of the difference of the kinetic energy from T/2: ⟨ E⟩ -T/2. Figure 12(a) shows ⟨ E⟩ -T/2 for several T's at d=0.01 and v_0=0.01. The excess kinetic energy weakly decreases with T for T>0.1. Figure 12(b) shows  ⟨ E⟩-T/2 for several d's in the double-logarithmic scale at v_0=0.01 and T=0.01. The dashed line is 0.00005/d. The excess kinetic energy decreases as 1/d for small d. 



§ SUMMARY

We have studied coupled oscillators with inertia under the shear velocity boundary conditions on the rectangular lattice. We have found that the kinetic energy increases with time when d=0 or the dissipation is absent. We have calculated the probability distribution of the velocity. When d=0 or sufficiently small, the velocity distribution is well approximated at the Gaussian distribution or the Maxwell distribution. We have observed the deviation from the Gaussian distribution when d is large. The deviation from the Gaussian distribution is characteristic of the nonequilibrium state far from equilibrium. We have found that the velocity profile changes from the linear one to the localized one near the boundary as d is larger. In the coupled oscillators under a finite temperature, the localization becomes weaker, and the probability distribution of the velocity becomes closer to the Gaussian distribution. When v_0 is small and T is small but not zero, the coupled oscillators are a weakly nonequilibrium stochastic system. We have observed intermittent jumps of the kinetic energy. The jumps correspond to the phase slips in the y direction, and an annealing process of spatial fluctuations of the displacement z_i,j in the x direction occurs through the vortex-antivortex collisions. 

Our coupled oscillator model under the shear velocity boundary condition is an interesting system, since the model exhibits various nonlinear-nonequilibrium phenomena and is related to various research areas such as the coupled limit-cycle oscillators, the dislocation theory in solids, a block-spring model of earthquakes, and the nonequilibrium molecular dynamics. We have shown several numerical results in this paper, however, the theoretical understanding is not sufficient, which is left to future study. 
99
Kuramoto Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence (Springer, New York, 1984). 
Bonilla J. A. Acebrón, L. L. Bonilla, C. J. Perez Vincente, F. Ritort, and R. Spigler, Rev. Mod. Phys. 77, 137 (2005). 
Strogatz S. Strogatz, Physica D 143, 1 (2000). 
Tanaka H. Tanaka, A. J. Lichtenberg, and S. Oishi, Phys. Rev. Lett. 78, 2104 (1997).
Sakaguchi H. Sakaguchi and T. Matsuo, J. Phys. Soc. Jpn. 81, 074005 (2012).
Sakaguchi2 H. Sakaguchi, S. Shinomoto, and Y. Kuramoto, Prog. Theor. Phys. 77, 1005 (1987).
Hong H. Hong, H. Chate, H. Park, and L. H. Tang, Phys. Rev. Lett. 99, 184101 (2007).
Sakaguchi3 H. Sakaguchi, Phys. Rev. E 105. 054211 (2022).
Read W. T. Read, Dislocations in crystals (McGraw Hill, 1953).
Cottrell A. H. Cottrell Dislocations and plastic flow in crystals (Clarendon Press, 1953).
Sakaguchi4 H. Sakaguchi, Phys. Rev. E 106, 054154 (2022).
Scholz C. H. Scholz, The Mechanics of Earthquakes and Faulting (Cambridge University Press, Cambridge, 2002).
Langer J. M. Carlson and J. S. Langer, Phys. Rev. Lett. 62, 2632 (1989).
Sakaguchi5 H. Sakaguchi and S. Kadowaki, J. Phys. Soc. Jpn. 86, 074001 (2017).
Evans D. J. Evans and G. P. Morriss, Statistical Mechanics of Nonequilibrium Liquids (Academic Press, San Diego, 1990).
Ashurst W. T. Ashurst and W. G. Hoover, Phys. Rev. A 11, 658 (1975). 
Liem S. Y. Liem, D. Brown, and J. H. R. Clarke, Phys. Rev A 45, 3706 (1992).
Lepri S. Lepri, R. Livi, and A. Politi, Phys. Rev. Lett. 78, 1896 (1997).
Holian B. L. Holian, W. G. Hoover, and H. A. Posch, Phys. Rev. Lett. 59, 10 (1987).
Bingham E. C. Bingham, Fluidity and Plasticity (MacGraw-Hill, New York, 1922).
Salmon J. B. Salmon, A. Colin, S. Manneville, and F. Molino, Phys, Rev. Lett. 90, 228303 (2003). 
Rutter E. H. Rutter, Techtonophysics 303, 147 (1999).
 

