Molecular shear heating and vortex dynamics in thermostatted two-dimensional Yukawa liquids
MMolecular shear heating and vortex dynamics in thermostatted two-dimensionalYukawa liquids
Akanksha Gupta ∗ and Rajaraman Ganesh † Institute for Plasma Research, HBNI, Bhat, Gandhinagar - 382428, India
Ashwin Joy
Department of Physics, Indian Institute of Technology Madras, Chennai - 600036, India
It is well known that two-dimensional macroscale shear flows are susceptible to instabilities lead-ing to macroscale vortical structures. The linear and nonlinear fate of such a macroscale flow ina strongly coupled medium is a fundamental problem. A popular example of a strongly coupledmedium is a dusty plasma, often modelled as a Yukawa liquid. Recently, laboratory experimentsand MD studies of shear flows in strongly coupled Yukawa liquids, indicated occurrence of strongmolecular shear heating, which is found to reduce the coupling strength exponentially leading to de-struction of macroscale vorticity. To understand the vortex dynamics of strongly coupled molecularfluids undergoing macroscale shear flows and molecular shear heating, MD simulation has been per-formed, which allows the macroscopic vortex dynamics to evolve while at the same time, “removes”the microscopically generated heat without using the velocity degrees of freedom. We demonstratethat by using a configurational thermostat in a novel way, the microscale heat generated by shearflow can be thermostatted out efficiently without compromising the large scale vortex dynamics. Inpresent work, using MD simulations, a comparative study of shear flow evolution in Yukawa liquidsin presence and absence of molecular or microscopic heating is presented for a prototype shear flownamely, Kolmogorov flow.
I. INTRODUCTION
Micron-sized dust grains get highly charged when im-mersed in a conventional plasma because of continuousbombardment of electrons [1]. The average charge on asingle dust grain is typically ∼ e − e , where e isthe absolute electronic charge. There are many examplesof plasma in nature wherein, large sized grains interactwith an ambient plasma and play an important role. Forexample, comets, planetary rings, white dwarf, earth’satmosphere and in laboratory conditions such as plasmaprocessing reactors, plasma torch and fusion devices [2].These charged grains or dust particles interact via ashielded Coulomb interaction or Yukawa potential asthe ambient plasma shields the grain charge. Suchplasmas can be characterized by two non-dimensionalparameters κ = a/λ D (where a is inter-grain-spacingand λ D = λ i λ e / (cid:112) λ i + λ e is Debye length of backgroundplasma, λ i , λ e are Debye length of electron and ion re-spectively) and coupling parameter Γ = Q d / (4 πε ak B T d )wherein Q d and T d are charge and temperature of grain.High charge on grain enhances the grain-grain potentialenergy and hence the coupling parameter Γ. Dependingupon κ and Γ values, such strongly coupled grains canexhibit solid-like (for Γ > Γ m ) [3–5], liquid-like (for1 < Γ < Γ m ) [6] and gas-like (for Γ << m is the liquid to solid phase transition point ∗ Institute for Plasma Research, India † [email protected] In laboratory experiments, grain dynamics can bevisualized (by unaided eye) and tracked by opticalcameras [3]. Using MD simulation, where the interactionbetween dust grain is modelled by shielded Coulomb orYukawa potential, various properties of grain medium,such as phase transition, instability, transport, graincrystallization, determination of transport coefficientssuch as shear and bulk viscosities [7], Maxwell relaxationtime [8], heat conduction, wave dispersion [9], selfdiffusion [10], fluid instability such as shear driven,Kelvin-Helmholtz instability [11] has been addressed.In recent past, strong peaks of temperature have beenobserved at the velocity shear location in macroscopicshear flows of grain generated by external laser-drive inYukawa liquid experiments [12, 13] and also in MD sim-ulations [14, 15]. It was seen that shear flows in Yukawaliquids lead to generation of heat at the microscopicscales resulting into the formation of a heat front [15].Using Kolmogorov flow [16–20] as a initial shear flow inYukawa liquids, it was also demonstrated that molecularshear heating destroys the vortex dynamics and reducethe coupling strength exponentially [14].In the present work, we investigate whether or not it ispossible to address macroscale vortex dynamics usingMD simulation and at the same time maintain the grainbed at the desired temperature. In the following, weconsider a prototype shear flow namely Kolmogorovflow, which has been studied in the context of laminar toturbulent transition [14] wherein it was shown that theaverage coupling strength decreases exponentially withtime due to molecular shear heating. We propose anddemonstrate here that using MD simulation and ther-mostat based on configurational space degree of freedom a r X i v : . [ c ond - m a t . s o f t ] M a y [21–24] (also called profile unbiased thermostat or PUT),it is possible to “remove” heat yet study macroscalevortex dynamics. Below we show the correspondencebetween PUT and a velocity scaling or profile basedthermostat (PBT) for the case of Kolmogorov shearflows.In the present work, we have introduced the method ofthermostatting namely, configurational thermostat. Inthe past, Rugh [21] and Butler[22] presented a method ofcalculating the temperature of a Hamiltonian dynamicalsystem. This method of calculating temperature onlydepends upon the configurational information of thesystem, hence named configurational temperature.Influenced by the concept of configurational temperaturea new method of thermostatting namely configurationalthermostat, has been introduced by Delhommelle andEvans [25, 26], Patra and Bhattacharya [27] and, Bragaand Travis [23, 24]. Due to its relative simplicityin implementation, we have chosen the Braga-Travisversion of PUT. This amounts to invoking appropriateLagrange multipliers that will efficiently couple thegrain bed to a configurational thermostat. Using thisPUT we demonstrate that the average coupling strengthcan be controlled without compromising the effects ofstrong correlations on the macroscopic shear flow andvortex dynamics. A detailed comparison of the evolutionand dynamics of Kolmogorov flow in the presence andabsence of molecular shear heating, its effect on lineargrowth rate, non-linear saturation and transition fromlaminar to turbulence flow will be presented.This paper is presented in the following fashion. In Sec.II and III, we describe configurational temperature andconfigurational thermostat used in our work. In Sec. IV,molecular shear heating effects on shear flow in presenceand absence of shear heating is described. In the last sec-tion Sec.V, we conclude and discuss about future work. II. CONFIGURATIONAL TEMPERATURE
Thermodynamic temperature is a estimate of theaverage or random kinetic energy of the particles in asystem. According to kinetic theory of gases, the kinetictemperature can be expressed as: k B T kinetic = 1 N d N (cid:88) i =1 m i v i k B is Boltzmann constant and m i , v i are massand instantaneous velocity of “ i th ” particle respectively, N and d are number of particles and dimensionality ofthe system. Apart from this definition of temperature,assuming ergodicity, one can also define temperature bya purely dynamical time averaging of a function whichis related to the curvature of energy surface. For exam- ple, H. H. Rugh [21], in 1997, presented a dynamical ap-proach for measuring the temperature of a Hamiltoniandynamical system. Rugh argues, using statistical ther-modynamics for any closed Hamiltonian system, that thekinetic and configurational temperature would be asymp-totically identical for a closed system. His definition ofconfigurational temperature is:1 k B T config = ∇ · ∇ H (cid:107) ∇ H (cid:107) (2)where H is the Hamiltonian of classical dynamical sys-tem. Butler [22] and others [28, 29] generalised Rugh’sidea for any function B ( Λ ) of phase space, such that k B T = (cid:68) ∇ H ( Λ ) · B ( Λ ) ∇ · B ( Λ ) (cid:69) (3)where Λ = { q , q ......q N , p , p ......p N } is the phasespace vector and ( q j , p j ) are 6 N generalized coordi-nates for conjugate positions and momenta respectively.The Hamiltonian of the system H ( Λ ) = K ( p j ) + U ( q j ), where K ( p j ) and U ( q j ) are kinetic and po-tential energy of the dynamical system respectively.In Eq. 3, B ( Λ ) can be any continuous and differen-tiable vector in phase space. For example, if onechooses B ( Λ ) as B ( Λ ) = B ( p , p , ..., p N , , , ..., N/ k B T = (cid:80) Ni =1 p i / (2 m i ). Similarly when B ( Λ ) = B (0 , , ..., , q , q , ..., q N ) = − ∇ U (say), Eq.3 gives con-figurational temperature in potential form k B T config = (cid:68) ( ∇ U ) ∇ U (cid:69) (4)In a closed system, these kinetic and configurational tem-peratures are expected to be asymptotically identical. Kinetic and configurational temperatures of strongly coupledYukawa liquid
Due to slow dynamics of grain medium, we consider am-bient plasma properties to be invariant and model onlygrain dynamics, which is strongly coupled. We first calcu-late the configurational temperature of the grain mediumwhere, grains interact via screening Coulomb or Yukawapotential : U ( r i ) = Q d π(cid:15) N (cid:88) j (cid:54) = i e − r ij /λ D r ij (5)where r ij = | r i − r j | is the inter particle distance of i th and j th particle. The N -body problem is then nu-merically integrated using our parallelized MD code [30].We have calculated configurational temperature by us-ing Eq. 4 where, interaction potential U ( r i ) is given byYukawa potential as in Eq. 5. Time, distance and en-ergy are normalized to inverse of dust plasma frequency √ ω − pd = ω − , mean inter-gain spacing a and averageCoulomb energy of dust particle Q d / (4 πε a ) respectively.Therefore, all physical quantities appearing hereafter inpresent paper are non-dimensional. In our simulations,for a given dust particles density ¯ n , the size of the systemis decided by the total number of particles. For ¯ n = 1 /π , N = 2500 and κ = 0 .
5, we first bring our system to thedesired Γ using a Gaussian thermostat [31]. After this welet the system evolve under micro-canonical conditionsand notice excellent agreement between temperatures ob-tained from both kinetic and configurational degrees offreedom (see Fig. 1). The absence of any noticeable driftin temperature even without a thermostat indicates goodnumerical stability of our time integration. C oup li ng pa r a m e t e r Γ ( t ) FIG. 1. color online: Kinetic (blue) and configurational(red) Γ( t ) extracted under micro-canonical conditions for theYukawa liquid previously equilibrated at Γ = 50 , , κ = 0 . III. CONFIGURATIONAL THERMOSTAT
In our previous work on shear flows [11, 14, 32], westudied the spatio-temporal evolution of instabilitiesas an initial value problem. As no attempt was madeto control the temperature of the liquid during thesimulation, the flow evolved under adiabatic conditionsand the shear heat generated due to viscosity remainedin the system. This led to a gradual increase in overalltemperature eventually resulting in short lifetimes oflarge scale vortex structures. Thus to address macroscalevortex dynamics with longer lifetimes it is highly de-sirable to “remove” this excess heat from the shearlayer without altering the physics of the problem. Inconventional MD simulations, thermostats are generally used to maintain the temperature of the system at adesired value in a canonical ensemble. For eg. in a typ-ical Gaussian thermostat [31], a Lagrangian multiplieris invoked for instantaneous velocities and equationsof motion are augmented (in the Nose-Hoover sense)with a velocity dependent non-holonomic constraint.While the trajectory of the system generated so forthstrictly conforms to the iso-kinetic ensemble it can beshown that the observed thermodynamic behavior of thesystem corresponds very well to that of the canonicalensemble in thermodynamic limit. As can be expected,such velocity scaling based thermostats can work only atlow shear rates and are immediately rendered useless athigh shear rates where secondary flows usually develop.It is the purpose of this paper to show that a PUT canbe efficiently applied to control the temperature of shearflows and below we present the details of our protocol.Once a macroscale flow is superimposed onto the ther-malized grain bed, the instantaneous particles velocitiescontain information regarding the “thermal” and the“flow (or average)” parts. Especially at high Reynoldsnumber regime where secondary flows usually develop,it becomes impossible to control temperature usingsuch PBTs which rely only on velocity scaling. Thuscontrolling the “thermal” component of velocity andletting mean component evolve is impossible usingthermostats which use augmented velocity equation asin a Gaussian thermostat.
Is it then possible to “ther-mostat” a system with N particles without modifyingthe instantaneous velocities of particles ? The answeris yes. As discussed earlier, a novel method of ther-mostatting, namely configurational thermostat has beenproposed amongst others, by Braga and Travis [23, 24],which controls the temperature by using augmentedequations of motions for the instantaneous particlepositions without disturbing the instantaneous velocityof particles. In this PUT, instead of a non-holonomicconstraint, a holonomic constraint is augmented to theequation of motion. The temperature of the system isthen calculated by using configurational definition [Eq. 4]that agrees well with the kinetic temperature calculatedusing Eq. 1. To better understand the subtle ways inwhich the configurational thermostat differs from its ki-netic counterpart, we describe both these schemes below.The equation of motion corresponding to kinetic Nose-Hoover thermostat are: ˙ r i = p i m i (6)˙ p i = − ∂U∂ r i − η v i (7)˙ η = 1 Q η (cid:18) N (cid:88) i =1 m i v i − k B T (cid:19) (8)where m i , r i , p i and T are mass, position, momentumof “ i ”-th particle and desired temperature respectively.Lagrange multiplier η is a dynamical variable and Eq. 6-Eq. 8 are the new augmented equations of motion. Q η is damping constant or effective mass [33, 34]. In thesame way, for “ i ”-th particle the augmented equationsof motion corresponding to configurational temperaturebased Nose-Hoover thermostat as defined by Braga andTravis [23, 24] are : ˙ r i = p i m i − µ ∂U∂ r i (9)˙ p i = − ∂U∂ r i (10)˙ µ = 1 Q µ (cid:18) N (cid:88) i =1 (cid:16) ∂U∂ r i (cid:17) − k B T N (cid:88) i =1 ∂ U∂r i (cid:19) (11)In above equations, the Lagrange multiplier µ is adynamical variable. In configurational thermostat, Q µ is an empirical parameter which behaves as the effectivemass associated with thermostat. The value of Q µ decides the strength of coupling between system andheat-bath. We find that the value of damping constantor effective mass Q µ is sensitive to the desired couplingparameter Γ values. For a desirable Γ, we find thatfor low values of damping constant, the system arrivesat the desired Γ faster and vice versa. In the limit Q µ → ∞ the configurational thermostat is de-coupledand formulation becomes micro-canonical ensemble.For the range of Γ values studied here, we find that Q µ = 2 × allows steady state with in time t = 100and hence is used throughout.Using a PUT, it is depicted in Fig.[2] that kinetic andconfigurational coupling parameter Γ follow the samebehaviour in canonical and in micro-canonical run aswell.In the following section, the evolution of shear flownamely, Kolmogorov flow, in Yukawa liquids in the pres-ence and absence of microscopic or molecular shear heatis presented. IV. KOLMOGOROV FLOW AS AN INITIALVALUE PROBLEM IN YUKAWA LIQUID WITHAND WITHOUT MOLECULAR HEATGENERATION
To study the shear flow evolution and vortex dynamicsfrom the microscopic dynamics in presence and absenceof heat generation phenomena, we have studied Kol-mogorov flow [16] as an initial value problem. Thissimply implies that a flow profile U is loaded only at t = 0 and no attempt is made to control the mean flow at
50 100 150 200 250 300 350 400 450 500050100150200250 Time (t) C oup li ng P a r a m e t e r [ Γ ( t ) ] PUT "OFF"PUT "ON"
FIG. 2. color online: Kinetic (blue) and configurational (red)Γ vs time. Parameters used: Q µ = 2 × , N = 2500 and κ = 0 . later times. At time t >
0, we use a PUT to maintain thedesired temperature. The loaded shear profile has theform U ( x, y )= U cos(2 πn x/L x )(1 + δ cos(2 πmy/L y ))ˆ y ,where the magnitude of equilibrium velocity U = 1,spatial period number n = 3 , magnitude of per-turbation δ = 0 .
01 , perturbed mode number m = 2[see Fig.3]. Coupling strength at time t = 0 isΓ = Γ( t = 0) = 50, for which calculated thermalvelocity is v th = (cid:112) / Γ = 0 .
2, which is much smallerthan the equilibrium velocity ( U = 1). Unlike earliersection Sec. III, here we have considered large number ofparticles N d = 62500 to study large-scale hydrodynamicphenomena using MD simulation. Due to the large sizeof the simulation box L x = L y = L = 443 .
12, we donot consider Ewald sums [35]. The non-dimensionalscreening parameter κ is 0 .
5. It is estimated that thesound speed of the system for Γ = 50 and κ = 0 . U . Hence the shear flow isconsidered to be “subsonic” in nature.To see the effect of molecular shear heating, we com-pare our present study with the earlier work [14]. Pre-viously, it has been demonstrated that after superpo-sition of the Kolmogorov flow, the coupling parameterwas found to weaken under adiabatic conditions due tomolecular shear heating. In the present paper, we showthat this average coupling parameter when coupled tothe configurational thermostat is approximately constantand close to desired Γ. In Fig.4, we have plotted ourearlier results along with the current results of averagecoupling parameter. In Fig. 4( a ) the system is thermallyequilibrated up to time t = 300 using a conventionalGaussian thermostat. In ( b ) we show the microcanonicalrun in the time interval 300 < t < t = 600 wesuperpose the shear flow just once and then observe thesystem both with PUT “ON” ( c ) and with PUT “OFF”( d ). As clearly seen from ( d ), the heat generated dueto shear flow remains within the system and weakens thecoupling strength Γ. We find this decay to be exponentialin time (see fit). This is in stark contrast to the regime( c ) where this excess heat has been “removed” from thesystem through heat-bath, which is turn facilitates Γ toremain constant. [We repeated the study by replacingthe Gaussian thermostat in a with a PUT and foundidentical results]. We also see from Fig.5 that both theconfigurational and kinetic coupling parameter are closeto the initial value Γ = 50. FIG. 3. color online: Kolmogorov velocity profile. Arrowsshow the direction of local flow.
Macroscopic quantities from microscopic information(process of “fluidization”):
A mesh-grid of size 55 ×
55 is superimposed on theparticles of the system to calculate the macroscopic or“fluid” variables. Average local fluid velocities along x and y directions are calculated as ¯ U x = (1 /N b ) (cid:80) N b i =1 v ix ,¯ U y = (1 /N b ) (cid:80) N b i =1 v iy , where v ix and v iy are individualinstantaneous particle velocities along x and y directionand N b is the total number of particles present in an indi-vidual bin. Each bin contains approximately 20 particles[ N b = N d / (55 × (cid:39)
20 with N d = 62500]. From aver-age local velocities, we calculate the average local vortic-ity ¯ ω ( x G , y G ) = ∇ × ¯ U and the average local temperature¯ T ( x G , y G ) = (2 / (cid:80) N b i =1 (cid:0) ( v ix − ¯ U x ) + ( v iy − ¯ U y ) (cid:1) /N b C oup li ng P a r a m e t e r [ Γ ( t ) ] (a) (b) (c)(d) FIG. 4. color online: Coupling parameter vs time. Thesystem is evolved (a) coupled to Gaussian thermostat (b) un-der micro-canonical conditions (c) with PUT “ON” and flowsuperposed at t = 600 (d) with PUT “OFF”. We provide afit (red) to show that weakening of Γ is indeed exponential intime. Parameters used: Γ = 50. C oup li ng pa r a m e t e r [ Γ ( t ) ] Kinetic Γ (t)Configurational Γ (t)No shear flow Shear flow FIG. 5. Coupling parameter vs time plot with PUT always“ON”. Shear flow is superimposed at time t = 1000. at the Eulerian grid location ( x G , y G ). In Fig.[6] timeevolution of “fluidized” vorticity after superposition ofshear flow both in the presence (Left panel) and absence(Right panel) of PUT. We clearly see that shear heatingdestroy vortex structures thus resulting in their shorterlifetimes compared to the case when PUT is “ON” wherethe lifetimes of these vortices become significantly longer.In Fig.7, we show our data for ¯ T ( x G , y G ) extracted atdifferent times with shear flow imposed. In the ab-sence of PUT (shown with symbols), the temperaturefirst increases and then eventually saturates at a par-ticular value at long times. This is in contrast to thecase when PUT is present where the temperature firstincreases (discussed later) and then saturates at the ini-tial temperature at long times. The value of initial tem-perature taken was 0.02 (Γ = 50). In Fig.[8] we haveplotted the y -component of velocity at various times inthe absence (Left panel) and presence (Right panel) ofPUT. We find that in both the cases the velocity profileshave not changed much between the linear and non-linearregimes of shear flow evolution. They remain qualita-tively similar and differ only quantitatively. The per-turbed x component of kinetic energy δE kx is obtainedfrom the expression below and we use it to calculate thegrowth rate shown in the Fig.[9]. (cid:12)(cid:12)(cid:12)(cid:12) δE kx ( t ) E kx (0) (cid:12)(cid:12)(cid:12)(cid:12) = (cid:82) (cid:82) [ v x ( t ) − v x (0)] dxdy (cid:82) (cid:82) v x (0) dxdy (12)We find that the calculated growth rate γ d in the absenceand presence of PUT are very close with the differencebeing only marginal ( < R = U l ¯ n/η , where l , η are theshearing length and initial shear viscosity of the flow re-spectively. Here, the initial value of shear viscosity η iscalculated using the Green-Kubo formalism [8, 37] beforeKolmogorov flow superimposed. It is depicted from figurethat for given value of Γ and κ , flow is neutrally stablebelow R < R c , where R c is critical value of Reynoldsnumber and for R > R c flow becomes unstable and even-tually turbulent. Such laminar to turbulent transition inour system might be a trans-critical bifurcation [38]. In-terestingly, we find that higher values of coupling param-eter Γ decreases the critical value of Reynolds number R c .Also the critical value of Reynolds number R c is foundto be independent of heat generation. It is evident fromFig.[10] that the growth-rate of perturbed mode ( m = 2)is not affected by molecular or microscopic heating, whichshows that the suppression of heat generation does notmodify the shear flow dynamics in early phase of simula-tion. V. CONCLUSIONS
In the present paper, taking Kolmogorov flow as aninitial condition, the role of molecular shear heating in strongly coupled Yukawa liquids has been investi-gated using configurational thermostat and results arecompared with our earlier work using micro-canonicalensemble [14]. In both the cases, we observe that thelaminar to turbulent transition of Kolmogorov shear flowcrucially depends upon the critical value of Reynoldsnumber R c in the presence and absence of heat gener-ation. Parametric study of growth-rate of perturbedmode over the range of Reynolds number shows theneutrally stable and unstable nature of Yukawa fluidsundergoing Kolmogorov flow for R < R c and R > R c case respectively. It is important to note that thecritical value of Reynolds number in the presence andabsence of heat generation is nearly same. Molecularshear heat is found to decrease the coupling strengthexponentially in time and hence destroys the secondarycoherent vortices. In this work, using the method ofconfigurational thermostat, it has been demonstratedthat the average or global temperature of the system canbe maintained at a desirable value in spite of molecularshear heating. This in-turn is found to help sustainthe secondary coherent vortices dynamics for a longertime. In the non-linear states obtained with PUT “ON”,spatially non-uniform profiles of temperature is observedin the regions of strong velocity shear. However, av-erage temperature of the system is controlled by theconfigurational thermostat. For example when, PUT is“ON”, it is observed that the peaks of local temperatureprofile at the shear flow location, are much lesser inmagnitude and global average temperature of the systemis maintained as compared to the case with PUT “OFF”.While it has been unambiguously demonstrated herethat configurational thermostats [Eq. 9-Eq. 11] areindeed effective in controlling the average or globaltemperature, it is desirable to be able to maintain thelocal and global temperatures at a same values. Workis under way to construct a novel set of equations ofmotion using an improved configurational thermostat,which will be reported in a future communication.We propose that the calculation of temperature usingthe definition of configurational temperature may bean important and useful alternative for those labora-tory experiments of strongly coupled dusty plasma inwhich instantaneous positions of particles are more ac-curately measured rather than instantaneous momenta.We strongly believe that this exercise should be under-taken using experimentally obtained instantaneous posi-tions { ¯ x i ( t ) } i =1 ,N and compare the same with conven-tional kinetic temperature. ACKNOWLEDGMENTS
Authors (A.G and R.G) acknowledge valuable discus-sions on configurational thermostat with Harish Charanand Rupak Mukherjee.
FIG. 6. color online: “Fluid” vorticity ( ω = ∇ × U ) contour plots. Color bars show the magnitude of local vorticity.Parameters used: perturbation mode m = 2, equilibrium spatial period number n = 3, Γ = 50, κ = 0 .
5, initial Reynoldsnumber R = 235 .
149 and shear velocity U =1. Left Panel: The micro scale heating quickly destroy the vorticity structureswhen PUT is ”OFF”. Right Panel: When PUT is ”ON”, vortex structures sustain for longer time.[1] N. N. Rao, P. K. Shukla, and M. Y. Yu, Planet. Space.Sci , 543 (1990)[2] U. de Angelis, Physics of Plasmas , 012514 (2006)[3] H. Thomas, G. E. Morfill, V. Demmel, J. Goree, B. Feuer-bacher, and D. M¨ohlmann, Phys. Rev. Lett. , 652 (Aug1994)[4] J. H. Chu and L. I, Phys. Rev. Lett. , 4009 (Jun 1994)[5] G. E. Morfill and A. V. Ivlev, Rev. Mod. Phys. , 1353(Oct 2009)[6] V. E. Fortov, V. I. Molotkov, A. P. Nefedov, and O. F.Petrov, Physics of Plasmas , 1759 (1999)[7] V. Nosenkol and J. Goree, Phys. Rev. Lett. (2004)[8] A. Joy and A. Sen, Phys. Rev. Lett. , 055002 (Feb2015)[9] H. Ohta and S. Hamaguchi, Phys. Rev. Lett. , 6026(2000)[10] H. Ohta and S. Hamaguchi, Phys. Plasmas (2000)[11] A. Joy and R. Ganesh, Phys. Rev. Lett. , 215003(May 2010)[12] W.-T. Juan, M.-H. Chen, and L. I, Phys. Rev. E ,016402 (Jun 2001)[13] Y. Feng, J. Goree, and B. Liu, Phys. Rev. Lett. (2012)[14] A. Gupta, R. Ganesh, and A. Joy, Physics of Plasmas , 103706 (2015)[15] A. Joy and R. Ganesh, Physics of Plasmas , 083704 (2011)[16] L. D. Meshalkin and Y. G. Sinai, J. Appl. Math. Mech , 1700 (1961)[17] A. M. Obukhov, Russian Mathematical Surveys , 113(1983)[18] L. D. Landau and E. M. Lifshitz, Fluid Mechanics (Perg-amon, Oxford)(1984)[19] D. H. Kelley and N. T. Ouellette, American J. of Physics , 267 (2011)[20] I. Bena, M. M. Mansour, and F. Baras, Phys. Rev. E ,5503 (1999)[21] H. H. Rugh, Physical Review Letters , 772 (feb 1997)[22] B. D. Butler, G. Ayton, O. G. Jepps, and D. J. Evans,Journal of Chemical Physics , 6519 (1998)[23] C. Braga and K. P. Travis, The Journal of ChemicalPhysics (2005)[24] K. P. Travis and C. Braga, The Journal of ChemicalPhysics , 014111 (2008)[25] J. Delhommelle and D. J. Evans, Molecular Physics ,1825 (2001)[26] J. Delhommelle and D. J. Evans, The Journal of Chemi-cal Physics , 6016 (2002)[27] P. K. Patra and B. Bhattacharya, Journal of ChemicalPhysics (2014)[28] Y. Han and D. G. Grier, The Journal of Chemical Physics , 064907 (2005) −250 −200 −150 −100 −50 0 50 100 150 200 2500.010.020.030.040.050.060.070.080.090.1 x G T e m p e r a t u r e t=0t=250t=500t=1000t=0t=250t=500t=1000 FIG. 7. y - averaged temperature profiles ( ¯ T ( x G , y G , t )) at various times after the shear flow is superimposed. Parameters usedare Γ = 50, U = 1 and κ = 0 .
5. Symbols show the temperature profile with PUT “OFF” and solid lines show temperaturevariation with PUT “ON”.[29] A. C. Bra´nka and S. Pieprzyk, Computational methodsin sci. and tech. , 119 (2011)[30] A. Joy and R. Ganesh, Phys. Rev. E , 056408 (Nov2009)[31] D. J. Evans and O. Morriss, Computer Physics Reports , 297 (1984), ISSN 0167-7977[32] A. J. and R. Ganesh, Phys. Rev. Lett. ,135001 (Mar 2011), http://link.aps.org/doi/10.1103/PhysRevLett.106.135001 [33] D. J. Evans and B. L. Holian, The Journal of ChemicalPhysics (1985)[34] G. J. Martyna, M. L. Klein, and M. Tuckerman, The Journal of Chemical Physics (1992)[35] G. Salin and J.-M. Caillol, Phys. Rev. Lett. , 065002(Jan 2002)[36] S. A. Khrapak and H. M. Thomas, Phys. Rev. E (2015)[37] J. P. Hansen and I. McDonald, Theory of Simple Liquids:With Applications to Soft Matter, 4th ed. (Academic,Oxford, 2013)[38] P. G. Drazin, Intoduction to Hydrodynamic Stabil-ity (Cambridge text in applied mathematics UniversityPress, Cambridge, England)(2002)[39] P. Hartmann, G. J. Kalman, Z. Donk´o, and K. Kutasi,Phys. Rev. E , 026409 (Aug 2005) x G -300 -200 -100 0 100 200 300 y - c o m p o n e n t o f v e l o c i t y ( v y ) -1-0.8-0.6-0.4-0.200.20.40.60.81 t=0t=500t=750t=1000 −300 −200 −100 0 100 200 300−1−0.8−0.6−0.4−0.200.20.40.60.81 x G y - c o m p o n e n t o f v e l o c i t y ( v y ) t=0t=500t=750t=1000 FIG. 8. Temporal evolution of y averaged velocity v y ( x G ) profile for Γ = 50, equilibrium velocity magnitude U = 1, screeningparameter κ = 0 . t here are shown after superimposition of shear flow. Time (t) l og δ E k Perturbed K.E [PUT "OFF"]Perturbed K.E [PUT "ON"]
FIG. 9. Perturbed kinetic energy in linear-log scale withand without PUT. Calculated growth rates from simulationare 5 . × − and 6 . × − for PUT “OFF” and “ON” caserespectively. −3 Reynolds number G r o w t h - R a t e Γ =50, PUT "OFF " Γ =50, PUT "ON" Γ =130, PUT "OFF" Γ =130, PUT "ON" R =235.149 FIG. 10. Growth-Rate vs initial Reynolds number (cid:16) R = U l ¯ nη (cid:17) plot showing trans-critical kind of bifurcation forscreening parameter κ = 0 .
5. The non-linear vortex struc-tures in Fig. 6 have been obtained for R = 235 ..