Kelvin-Helmholtz instability at proton scales with an exact kinetic equilibrium
A. Settino, F. Malara, O. Pezzi, M.Onofri, D. Perrone, F. Valentini
aa r X i v : . [ phy s i c s . s p ace - ph ] M a y Draft version May 20, 2020
Typeset using L A TEX twocolumn style in AASTeX62
Proton Kinetics in Kelvin-Helmholtz instability
A. Settino, F. Malara, O. Pezzi,
2, 3
M.Onofri, D. Perrone, and F. Valentini Dipartimento di Fisica, Universit`a della Calabria, 87036 Rende (CS), Italy Gran Sasso Science Institute, I-67100 LAquila, Italy INFN/Laboratori Nazionali del Gran Sasso, I-67100 Assergi (AQ), Italy TAE Technologies Inc., PO Box 7010, Rancho Santa Margarita, CA 92688, USA ASI – Italian Space Agency, via del Politecnico snc, 00133 Rome, Italy
Submitted to ApJABSTRACTThe Kelvin-Helmholtz instability is a ubiquitous physical process in ordinary fluids and plasmas,frequently observed also in space environments. In this paper, kinetic effects at proton scales in thenonlinear and turbulent stage of the Kelvin-Helmholtz instability have been studied in magnetizedcollisionless plasmas by means of Hybrid Vlasov-Maxwell simulations. The main goal of this work is topoint out the back reaction on particles triggered by the evolution of such instability, as energy reacheskinetic scales along the turbulent cascade. Interestingly, turbulence is inhibited when Kelvin-Helmholtzinstability develops over an initial state which is not an exact equilibrium state. On the other hand,when an initial equilibrium condition is considered, energy can be efficiently transferred towards shortscales, reaches the typical proton wavelengths and drives the dynamics of particles. As a consequenceof the interaction of particles with the turbulent fluctuating fields, the proton velocity distributiondeviates significantly from the local thermodynamic equilibrium, the degree of deviation increasingwith the level of turbulence in the system and being located near regions of strong magnetic stresses.These numerical results support recent space observations from the Magnetospheric MultiScale missionof ion kinetic effects driven by the turbulent dynamics at the Earth’s magnetosheath (Perri et al. , 2020,JPlPh, 86, 905860108) and by the Kelvin-Helmholtz instability in the Earth’s magnetosphere (Sorriso-Valvo et al. , 2019, PhRvL, 122, 035102). INTRODUCTIONThe Kelvin-Helmholtz instability (KHI) is a phe-nomenon that can develop in both fluids and plasmas, inconfigurations where velocity shears are present. Dur-ing KHI, perturbations are generated in form of a chainof vortices, located along the shear layer, which growin time starting from infinitesimal fluctuations. In thecase of a magnetized plasma, the magnetic field has astabilizing effect with respect to KHI. Typically, a con-figuration is unstable when the jump in the bulk velocityacross the shear layer is larger than a threshold, whichis of the order of the component of the Alfv´en veloc-ity parallel to the bulk velocity (Chandrasekhar 1961).When unstable modes reach a sufficiently large ampli-tude, they start interacting among them, fragmenting
Corresponding author: Adriana [email protected] and generating structures at increasingly smaller scales.Moreover, vortices tend to merge forming larger coher-ent structures, moving part of the fluctuating energy tolarger scales. These phenomena lead to a final turbu-lent state where part of the kinetic energy associatedwith the velocity shear is dissipated. Therefore, KHIrepresents a way for a fluid or a plasma to give rise toa turbulent scenario and to convert large-scale motionenergy into heat.KHI has been considered in many natural systems,such as in terrestrial, heliospheric and astrophysi-cal contexts. For instance, (i) KHI has been ob-served at planetary magnetospheres (Kivelson & Chen1995; Seon et al. 1995; Fairfield et al. 2000, 2003;Hasegawa et al. 2004, 2006; Nykyri et al. 2006); (ii)it has been invoked to explain the penetration of solarwind into cometary ionospheres (Ershkovich & Mendis1983); (iii) it has been considered in turbulence mod-els at the interface between fast and slow solar wind
Settino et al. streams (Roberts et al. 1991, 1992); and (iv) it hasbeen observed in the solar corona at the surface ofcoronal mass ejections (Foullon et al. 2011). Moreover,the role of KHI has been also studied in the gener-ation of astrophysical jets in relativistic magnetizedplasmas (Hamlin & Newman 2013) or at the interfacebetween the accretion disk and the magnetosphere of aslowly rotating magnetized star (Lovelace et al 2010), aswell as in black holes and neutron stars (Li & Narayan2004). KHI is thought to be responsible for the plasmatransport across the Earth’s magnetopause, during pe-riods of both northward and southward orientation ofthe interplanetary magnetic field (Foullon et al. 2008;Kavosi & Raeder 2015).The unprecedented high-resolution observations con-ducted by the NASA Magnetospheric MultiScale (MMS)mission, launched in March 2015, have allowed to in-spect KHI onset at kinetic scales (Stawarz et al. 2016;Hwang et al. 2020).
In-situ measurements, supportedalso by numerical simulations, suggest that magneticreconnection induced by KHI breaks the frozen-incondition, thus favoring the solar-wind plasma entryinto the Earth’s magnetosphere (Nakamura et al. 2017;Eriksson et al. 2016; Sisti et al. 2019). Moreover, pri-mary and secondary KHI have been associated to thegeneration and shaping of flux ropes (Hwang et al. 2020;Zhong et al. 2018). The interconnection between turbu-lence development and KHI at the non-linear stagehas been recently studied by comparing MMS ob-servations with both magnetohydrodynamics (MHD)(Hasegawa et al. 2020; Nakamura et al. 2020) and hy-brid kinetic simulations (Franci et al. 2019). Finally,KHI is supposed to be dawn-dusk asymmetric owingto the different vorticity at the two flanks. This mid-latitude asymmetry has been investigated by means ofsimultaneous in-situ observations of THEMIS and MMSsatellites (Lu et al. 2019).KHI in magnetized plasmas has been widely studied invarious configurations. Several theoretical studies havebeen carried out within the MHD framework. The lin-ear stage of the instability, when unstable modes growexponentially in time, has been investigated for differ-ent spatial profiles of the bulk velocity u and density,and different orientations of the magnetic field B withrespect to u (see e.g. Axford 1960; Walker 1981; Miura1982; Contin et al. 2003). Moreover, the interplay withtearing instability has been also considered for inhomo-geneous magnetic field profiles (Wesson 1990). Disper-sive or kinetic effects come into play when the shearlayer thickness is of the order of ion length scales (ioninertial length and/or Larmor radius). These phenom-ena affect the growth rate of unstable modes, which in this regime depends on the relative orientation betweenmagnetic field and vorticity (Nagano 1979; Huba 1996;Cerri et al. 2013).The nonlinear evolution of KHI has been numericallystudied in a large number of investigations, using bothfluids (MHD, Hall-MHD and two-fluid) and kinetic ap-proaches. In fluid simulations, it has been shown thatviscosity generates momentum transfer between flows onthe two sides of the shear layer (Miura 1982). Moreover,in the case of perpendicular magnetic field, if the sim-ulation box is larger than the vortex length, an inversecascade takes place where KHI-generated vortices mergeforming structures at larger scales (i.e., vortex pairing)(Miura 1997, 1999a). This effect has been proposedas a way to follow the time evolution of KHI in non-periodic configurations (Mills et al. 2000; Wright et al.2000), such as at the Earth’s magnetopause (Miura1999b). In the fully nonlinear regime, secondary in-stabilities can develop, such as Rayleigh-Taylor, sec-ondary KH, or kink-like instabilities, which can com-pete with the pairing process leading to the disruption ofvortices (Matsumoto & Hoshino 2004; Nakamura et al.2004; Faganello et al. 2008). Furthermore, in config-urations where the in-plane magnetic field componentchanges sign across the shear layer, magnetic reconnec-tion can couple with KHI, thus creating a magneticconnection between the two sides of the shear layer,with consequences on the transport properties. How-ever, even when the in-plane magnetic component keepsthe same sign, reconnection takes place during the non-linear stage, leading to the formation of complex mag-netic topologies (for a detailed discussion see the reviewby Faganello & Califano (2017) and references therein).These phenomena take part to the more general prob-lem of reconnection in small-scale structures generatedby turbulence (Servidio et al. 2011a,b, 2012).In cases when KHI develops in collisionless plasmas atscales of the order of ion scales, such as in the Earth’smagnetosphere, kinetic simulations appear to be moresuitable than fluid approaches (Pritchett & Coroniti1984; Matsumoto & Hoshino 2006; Cowee et al. 2009;Matsumoto & Seki 2010; Nakamura 2010, 2011, 2013;Henri et al. 2013; Karimabadi et al. 2013). The in-terplay of KHI with other kind of instabilities, suchas lower-hybrid drift instability, has been very re-cently considered by Dargent et al. (2019). Further-more, kinetic effects can be important during thenonlinear stage of the instability, when vortices mixand a turbulent state develops. Indeed, kinetic sim-ulations of turbulence at ion scales have highlightedthe formation of small-scale structures in the phys-ical space closely related to the generation of out- roton kinetics in Kelvin-Helmholtz instability
Settino et al. to a larger extent, the nonlinear development of the in-stability, giving origin to a more developed turbulenceand larger values for the current density. These resultsare relevant in the perspective of correctly evaluatingthe spectral energy transfer and dissipation generatedby the KHI.The plan of the paper is the following: in Section 2 wedescribe the initial setup of the simulations with a focuson the equations of the model, the characteristics of theDFs and the perturbations introduced. An insight intothe derivation of the EE solution is also provided. InSection 3 we discuss simulations results, directly com-paring the EE and SM data. Finally, we give the con-clusions in Section 4. SIMULATION SETUP AND INITIALCONDITIONSTo perform the numerical analysis of the KHI, retain-ing kinetic effects at proton scales, we employed theHVM numerical code (Valentini et al. 2007). The HVMalgorithm solves numerically the Vlasov equation forthe proton DF, self-consistently coupled to the Maxwellequations for electromagnetic fields, while electrons aretreated as a massless fluid. We considered two shared-flow initial conditions, with two different initial protonDFs: the exact shared-flow HVM equilibrium distribu-tion, f ( EE ) , derived in Malara et al. (2018) and a SM dis-tribution, f ( SM ) . Results obtained starting from thesetwo initial conditions will be discussed and compared inSection 3.The HVM equations are numerically solved in a 2.5D-3V phase-space domain, that is, fully three-dimensionalin velocity space while, in physical space, all vectorshave three components depending only on two variables( x, y ). Quasi-neutrality condition is assumed and thedisplacement current is neglected in the Amp`ere equa-tion, in such a way to discard light waves.In dimensionless units, HVM equations are: ∂f∂t + v · ∇ f + ( E + v × B ) · ∇ v f = 0 (1) E = − u × B + 1 n j × B − n ∇ P e (2) ∂ B ∂t = −∇ × E ; ∇ × B = j (3)being f = f ( x, y, v x , v y , v z ) the proton DF, ∇ =( ∂ x , ∂ y ), ∇ v = ( ∂ v x , ∂ v y , ∂ v z ), E and B respectivelythe electric and magnetic fields, n and u respectivelythe proton density and bulk velocity, computed as thefirst two velocity moments of f , and j the total cur-rent density. For the electron pressure, we assume anisothermal equation of state, P e = nT e , in the case of SM initial condition, where n e = n p = n for the quasi-neutrality assumption and T e = ¯ T , being ¯ T the protontemperature far from the shear. On the other hand,for the EE initial condition, as extensively discussedin Malara et al. (2018), we need to relax the electronclosure in order to maintain the equilibrium, by treat-ing the electron pressure, P e , as a further independentquantity determined by the following equation: (cid:20) ∂∂t + ( u e · ∇ ) (cid:21) (cid:18) P e n γ e (cid:19) = 0 (4)where γ e = 5 / u e = u − j /n is the electron bulk velocity.In Eqs. (1)-(4), time is scaled by the inverse protoncyclotron frequency, Ω cp , velocities by the Alfv´en speed, v A = B / p π ¯ nm p (where B is the background mag-netic field, ¯ n the proton density away from the shearregions and m p the proton mass), lengths by the protonskin depth, d p = v A / Ω cp , the magnetic field by B , theelectric field by v A B /c , the density by ¯ n , and the elec-tron pressure by ¯ nm p v A . The development that followswill be expressed in terms of the above dimensionlessquantities.Spatial domain D = [0 , L x ] × [0 , L y ] ( L x = L y = L = 100) is discretized on a uniformly spaced grid with N x = N y = 256 grid points; periodic boundary con-ditions have been implemented in the spatial domain.Velocity-space domain is discretized on a uniform gridwith N v j = 71 ( j = x, y, z ) grid points in each direction.Vanishing boundary conditions have been implemented: f ( | v j | > v max ) = 0, being v max = 7 v th and v th = ( ¯ T ) / the proton thermal speed; the proton plasma beta is β = 2 v th /v A = 2.The unperturbed configuration is characterized by: (i)a sheared bulk velocity field u = u ( x ) e y , that is directedalong y and varies in the x direction; (ii) a perpendicu-lar uniform magnetic field B = B e z , with B = 1; and(iii) an electric field E = E ( x ) e x , whose profile E ( x ) isrelated to the bulk velocity u ( x ). In the above expres-sions e x , e y and e z are the unit vectors in the directionsof the three Cartesian axes.In the case of SM configuration the bulk velocity hasthe form u = U y ( x ) e y , where the function U y ( x ) de-scribes the double shear profile: U y ( x ) = U (cid:20) tanh (cid:18) x − x ∆ x (cid:19) − tanh (cid:18) x − x ∆ x (cid:19) − (cid:21) (5)Here, x = L x / x = 3 L x / x = 2 . U = 2 v A is thevelocity jump. We point out that the velocity shear hasbeen replicated along the x direction to satisfy periodicboundary conditions. roton kinetics in Kelvin-Helmholtz instability y component of the proton bulkvelocity, u y , as a function of x for the SM distributionfunction is reported as red dots in Fig. 1, where thepresence of the two shear layers is clearly visible. Figure 1.
The y component of the proton bulk velocity asa function of x for the distribution f ( EE ) (black curve) and f ( SM ) (red dots). Exact solution
Beside the SM DF, we considered the stationary so-lution found by Malara et al. (2018) for the system ofHVM equations. In the following, we briefly revisit thederivation and properties of such solution, while moredetails can be found in Malara et al. (2018). An exactstationary solution of the Vlasov equation, Eq. (1), canbe written as a function of constants of single particlemotion. Therefore, we consider the motion of a protonin the above electric and magnetic fields. The followingrelation between the x -position and the y -component ofthe particle velocity is found (in dimensionless units): v y ( t ) = W − x (6)where W is a constant determined by initial conditions.The particle motion in the x -direction corresponds tothat of a nonlinear oscillator, whose effective potentialenergy has the form: U eff ( x ) = Φ( x ; x ) + 12 ( x − W ) + 12 v y (7)where Φ( x ; x ) = − R xx E ( x ′ ) dx ′ is the electrostaticpotential which vanishes at a given position x , and v y = W − x . In Eq. (7) energies are normalizedto m p v A . Assuming that the bulk velocity profile is uni-form away from the shear layers, which corresponds toa uniform electric field, from Eq. (7) it follows thatthe particle motion along x is periodic within a poten-tial well. Therefore, we can define the guiding cen-ter position x c and velocity v yc as the average x po-sition and average v y velocity, respectively: x c = h x i t , v yc = h v y i t = W − x c . In particular, the point x whereΦ( x ; x ) is null is chosen as x c . The reduced energy isdefined as: E ( x, v x , v y , v z ) = 12 (cid:0) v x + v y + v z (cid:1) + Φ( x ; x c ) − v yc (8)The total energy (kinetic + potential) and v yc are bothconstants of motion. Therefore, E is another constantof motion, equal to the total energy minus the kineticenergy associated with the drift motion. We define adistribution function f ( EE ) ( x, v x , v y , v z ) = C exp (cid:20) − E ( x, v x , v y , v z ) v th (cid:21) (9)with C and v th constants. Since f ( EE ) is a combinationof constants of motions, it is an exact stationary solutionof the Vlasov equation. It can be shown that far from theshear layer f ( EE ) reduces to a shifted Maxwellian cen-tered around the drift velocity ( − E/B ) e y . The den-sity n associated to f ( EE ) is spatially uniform every-where except in the regions corresponding to the veloc-ity shears (see, for details, Malara et al. 2018).In the general case, the explicit form of f ( EE ) is nu-merically calculated on the grid in the 4D phase space { x, v x , v y , v z } . For each grid point the particle trajec-tory is integrated until it closes in the v x v y plane, cal-culating the corresponding values for the constants ofmotion: the guiding center position x c = h x i t and ve-locity v yc = h v y i t ; the kinetic energy; and the potentialΦ( x ; x c ). Those values are used to calculate f ( EE ) at thegiven grid point. Results show that the bulk velocity isdirected along y , i.e. u = u ( x ) e y , and hence the term − u × B in Eq. (2) is directed along x . Choosing a formfor the electric field, Eq. (2) is exploited to determinethe electron pressure P e . In particular, we adopted theexpression E ( x ) = − B U y ( x ), where U y ( x ) is the bulkvelocity associated with the SM distribution function inEq. (5).In Fig. 1 the corresponding profile of u y ( x ) from theexact equilibrium solution is plotted (black curve). Itcan be seen that the bulk velocity profiles correspond-ing to the exact solution and to the shifted Maxwellianare very close to each other. Nevertheless, we will showthat the time evolution of the KHI is different in the two Settino et al. cases. Larger differences are found in the density profile,which is uniform in the case of the shifted Maxwellian,while in the case of the exact solution it has a maxi-mum and a minimum, respectively localized at the twoshears. Moreover, f ( EE ) exhibits a clear temperatureanisotropy in regions close to the shears, being elongatedin a direction transverse to the background magneticfield, while reduces to a shifted Maxwellian far from theshears (Malara et al. 2018).2.2. Initial perturbation At t = 0, we perturbed the initial configurationthrough a broadband spectrum of bulk velocity fluctu-ations. Such perturbations have only y spatial depen-dence and are generated in the form of random noise.We excited the first 32 modes in the spectrum withrandom phases. For both EE and SM simulations, wesummed to the unperturbed function ( f ( EE ) or f ( SM ) ,respectively) the perturbation, shaped as a Maxwellianfunction shifted in the v x and v y directions; that is f ( EE ) = f ( EE ) + f or f ( SM ) = f ( SM ) + f , where f is defined (in scaled units) as follows: f ( y, v ) = n ( πβ ) / exp ( − [ v x − u x ( y )] β + − [ v y − u y ( y )] β − v z β ) ; (10)here, n = 0 .
01 is the amplitude of the perturbation, and u x = P i =1 cos ( k y,i y + ψ i ), u y = P i =1 sin ( k y,i y + φ i ),where k y,i = i π/L and ψ i and φ i random phases. Forboth the initial perturbed EE and SM distributions, theproton density n and the bulk velocity u can be writtenas: n = n + n ; u = n u + n u n + n (11)where n and u are the density and the bulk velocityof f , respectively. For small amplitude perturbations,the above equation for the bulk velocity can be Taylorexpanded in series of n /n leading to: u = u + n n u (12)Since n is constant and n is uniform (except in theshear regions for f ( EE ) ), the perturbed part of the bulkvelocity is largely shaped by u , for both EE and SMcases. NUMERICAL RESULTSDuring the linear phase of the instability, the exponen-tial growth of the energy E k y of the velocity Fourier com-ponents perturbed at t = 0 is observed. In Fig. 2, we report the time evolution of E k y = h| u k y ( x, t ) | i x , where u k y ( x, t ) is obtained by Fourier transforming u ( x, y, t )along the y direction and h· · · i x indicates average over x ∈ [0 , L/ Figure 2.
Time evolution of the first 5 Fourier componentsof the spectral kinetic energy E k y for the EE (top panel) andSM (bottom panel) simulations. Growth rates γ estimated during the early exponen-tial phase are plotted in Fig. 3 for both EE (black dots)and SM (red triangles) simulations, as functions of k y .Left and right panels report the growth rates of the firsteight Fourier components of the energy E k y averagedover x ∈ [0 , L/
2) and x ∈ [ L/ , L ), respectively. As it roton kinetics in Kelvin-Helmholtz instability Figure 3.
Growth rate γ of the first 8 Fourier components of E k y , as a function of k y for the EE (black dots) and SM (redtriangles) simulations. In the left (right) panel E k y has been averaged over x ∈ [0 , L/
2) ( x ∈ [ L/ , L )). can be appreciated from the two panels in Fig. 3, thedevelopment of the instability is not symmetric on thetwo shears: in particular, the fastest growing mode isnot the same at the two shears, as first and second mostunstable Fourier components are switched from left toright panel. Such asymmetry can be reasonably due todifferences in the sign of ω · B (the proton vorticitybeing ω = ∇ × u ) at the two shears (positive in cor-respondence of the left shear and negative at the rightone) (Henri et al. 2013).In the time evolution of the system, the initial expo-nential growth is followed by nonlinear saturation andlater by a transition to turbulence. This can be appreci-ated in Fig. 4, which shows the contour plot of | j | for theEE simulation at three different times. The left panelof this figure corresponds to the time of the late linearphase of the instability and displays the formation ofvortical structures in the shear regions (here the asym-metry between left and right shear is remarkable); inthe middle panel, corresponding to the nonlinear satura-tion phase, vortices in both shears start merging and fi-nally collapse in two distinct large-scale structures (rightpanel), in which short-scale filaments, whose size is fewproton skin depth, are generated.In order to quantify the level of turbulence in thesystem, in Fig. 5 we report the time evolution of themean squared current density h| j | i ( h· · · i meaning av-erage over the whole spatial domain D ) for the EE (blackcurve) and SM (red curve) simulations. Here, significantdifferences are recovered between EE and SM cases: infact, generation of turbulence seems to be inhibited inthe case of the SM initial condition, for which the sat-uration value of h| j | i is about one order of magnitude lower than in the EE case. To better point out this ef-fect, in Fig. 6 we show the omni-directional spectra ofmagnetic (top panel) and kinetic (bottom panel) energy,evaluated at the time in correspondence of the verticalblack (red)-dashed line in Fig. 5 for the EE (SM) sim-ulation. These spectra (Kolmogorov expectation k − / is indicated by a blue dashed line as a reference) clearlyshow a larger energy content (about an order of magni-tude in the inertial range) for the EE case (black lines)as compared to the SM case (red lines). Moreover, al-though in both EE and SM cases the spectral energyis peaked at low wavenumbers, a Kolmogorov-like spec-trum is observed for about a wavenumber decade.The inhibition of turbulence generation occurring inthe simulation with the SM initial condition may be dueto the fact that f ( SM ) is not an exact equilibrium DFfor the HVM equations in presence of a velocity shear.Indeed, as discussed in previous works (see, for exam-ple, Cerri et al. (2013) and Malara et al. (2018)), thisfeature naturally induces large-scale oscillations in theproton density, bulk speed, and also higher order mo-ments; these oscillations may lock the energy at lowwave numbers, preventing it from efficiently cascadingtowards smaller scales.As the HVM code retains kinetic effects on protons,the question arises whether the development of turbu-lence across d p produces deformations of the protonDF. In the following, we seek for local deviations fromMaxwellianity and for the generation of sharp gradientsin the proton velocity distribution. We recall that, at t = 0, f ( EE ) departs from a Maxwellian in the shearregions, while the perturbed SM initial condition, be-ing setup as a sum of two distinct Maxwellians, is not Settino et al.
Figure 4.
Contour plot of | j | at three different times in the EE simulation. Left panel corresponds at time t = 80, in thelinear phase of the evolution of the instability, middle panel corresponds to t = 180, when vortical structures start merging and,finally, right panel is at the end of the simulation, where vortices have collapsed in two large-scale structures and thin currentfilaments have been generated. Figure 5.
Time evolution of the mean squared current den-sity h| j | i for the EE (black curve) and SM (red curve) simu-lation; vertical black (red) dashed line indicates the time atwhich the maximum level of turbulence is reached in the EE(SM) simulation. a Maxwellian. Then, we investigate: (i) if distortionsfrom the Maxwellian shape increase as turbulence de-velops, and (ii) if these distortions remain confined inthe shear regions. In order to quantify deviations froma Maxwellian, we employ the non-Maxwellianity indi-cator introduced in Greco et al. (2012) and defined as: ǫ ( x, y, t ) = 1 n sZ ( f − g ) d v (13)where g is the Maxwellian DF associated with f , i.e.,which has the same velocity moments (density, bulk ve-locity and temperature) as f . In the left panel of Fig. 7 we show the time evolu-tion of h ǫ i ( h· · · i meaning average over the whole spa-tial domain D ), for both EE and SM simulations. At t = 0 h ǫ i starts from a non-zero value for both EEand SM simulations, as anticipated above. Both quan-tities then grow in time, indicating efficient generationof non-Maxwellian features during the EE simulation,while saturating after the initial growth in the case ofthe SM simulation. However, the saturation level of h ǫ i is larger for the EE case with respect to the SM one,this suggesting that the generation of non-Maxwellianfeatures in the DF is much more efficient in the formercase. In the right panel of the same figure, we presentthe scatter plot of h ǫ i versus h| j | i , showing that the in-crease of the non-Maxwellianity indicator appears to bewell correlated in time with the increase of the level ofturbulence in the system (time Pearson correlation coef-ficient is C t ≃ .
99) for the EE simulation (black dots).On the other hand, for the SM case (red dots) the cor-relation between the two quantities is not as high as inthe previous case ( C t ≃ . ǫ and | j | at a fixed instant of time. In Fig. 8, we reportthe contour plot of ǫ (left panel) and | j | (middle panel)at the time of the maximum level of turbulence in thesystem (vertical black-dashed line in Fig. 5). The spa-tial features of the two quantities are very similar, withpeaks of ǫ concentrated inside the vortical structures of | j | ; the horizontal cuts (right panel) of ǫ (red) and | j | roton kinetics in Kelvin-Helmholtz instability Figure 6.
Omnidirectional magnetic (top panel) and kinetic(bottom panel) energy spectra for the EE (black curve) andSM (red curve) simulation, taken at the time correspondingto the peak of h| j | i (see Fig. 5). Kolmogorov expectation k − / is indicated in both panels by a blue-dashed line as areference. (black), taken along the horizontal white-dashed linesin the left and middle panels of this figure, confirm thatthe peaks in the non-Maxwellianity indicator thicken in-side the vortical current structures.In Fig. 9 we show how the proton DF looks like at thetime of the maximum level of turbulence in the systemand at the spatial point where ǫ is maximum (black dotin the left panel of Fig. 8). Here, the 2D contour plots of f are shown in the ( v y , v z ) plane for v x = 0 (left panel),in the ( v x , v z ) plane for v y = 0 (middle panel) and in the( v x , v y ) plane for v z = 0 (right panel). Left and mid-dle panels of this figure display generation of significant temperature anisotropy, while the right panel shows pe-culiar deformations, with the generation of modulationsand sharp velocity gradients, driven by the interactionof particles with the field fluctuations.A deep analysis of both anisotropy and agyrotropy ofthe proton DF allows to provide quantitative informa-tion on the features observed in the contour plots of Fig.9. In the left panel of Fig. 10, we plot the spatial varia-tion of the anisotropy index A ( x, y ) = 1 − T ⊥ /T k , where T ⊥ and T k are the temperatures in the direction trans-verse and parallel to the local magnetic field, respec-tively. This plot is calculated at the time correspond-ing to the maximum level of turbulence in the system.Negative (positive) values of A correspond to T ⊥ > T k ( T ⊥ < T k ). In the right panel of the same figure, we showthe agyrotropy parameter √ Q , linked to the off-diagonalterms of the pressure tensor , where Q is defined as: Q = P xy + P xz + P yz P ⊥ + 2 P ⊥ P k ; (14)here, P ij are the components of the pressure tensor inthe reference frame in which one of the axes is along thelocal magnetic field (see Swisdak 2016, for more details).The agyrotropy parameter ranges from 0 to 1, where √ Q = 0 and √ Q = 1 correspond to fully gyrotropicconfigurations and maximum agyrotropy, respectively.It can be easily noticed that iso-contours of A and √ Q exhibit a pattern similar to those visible in the contourplots of ǫ and | j | (left and middle panels in Fig. 8), re-spectively. Indeed, A reaches its highest values in thecenter of each vortex (same as ǫ ), while √ Q achievesits highest value at the edges of the vortices (similar to | j | ). Moreover, in correspondence of the maximum valueof ǫ , √ Q displays highly non gyrotropic features of theproton DF and A > T ⊥ > T || (Malara et al. 2018).Finally, it is interesting to look at the shape of thevelocity DF in correspondence of three different valuesof A , corresponding to A > A = 0 and A <
0. In Fig.11, we report the three-dimensional velocity iso-surfaceof the proton DF at spatial point ( x, y ) = (53 . , .
4) for
A > x, y ) = (53 . , .
04) for A = 0 (mid-dle panel) and at the coordinate ( x, y ) = (53 . , . A < A = 0 (middle panel), where complex structures in the0 Settino et al.
Figure 7.
Left: time evolution of the non-Maxwellian indicator averaged over the whole spatial domain h ǫ i for the EE (blackcurve) and the SM (red curve) simulation. Right: scatter plot of h ǫ i as a function of h| j | i for the EE (black dots) and SM (reddots) simulation. Figure 8.
Contour plot of ǫ (left panel) and of | j | (middle panel) at the time of the maximum level of turbulence in the EEsimulation. In the right panel, cuts of | j | (black curve) and of ǫ (red curve), taken along the horizontal white-dashed paths inleft and middle panels. velocity space are visible, in spite of the temperatureisotropy. SUMMARY AND CONCLUSIONSIn this paper we have studied the nonlinear and turbu-lent stage of the KHI and the related kinetic effects pro-duced on particles, by means of Hybrid Vlasov-Maxwellsimulations at proton scales. In particular, we have con-sidered an unperturbed configuration where the back-ground magnetic field is perpendicular to the shear flow.This configuration is Kelvin-Helmholtz unstable regard-less of the value of velocity jump across the shear layer,at least in the MHD case.In kinetic descriptions of KHI, the unperturbed con-figuration has been often represented by means of a SMDF, though this is not a stationary solution. We haveshown that, when an exact equilibrium solution is cho-sen to initialize the system, relevant effects on the dy-namics of KHI appear. To highlight this point, we have compared KHI simulations with two different onsets forthe DF, namely (i) the EE solution and (ii) a SM DF,which is not a stationary solution. Due to spurious fluc-tuations, the non-stationary solution tends to inhibitturbulence that develops during the nonlinear phase ofKHI. The enhancement of turbulent activity in the EEsimulation can be observed in the spectra, where themagnetic and kinetic energy in the inertial range of theEE simulation are roughly one order of magnitude largerthan in the SM case. Moreover, considering the meansquare current density, which is mainly determined bysmall scales, we found that h| j | i reaches a much higherlevel in the EE simulation when compared with the SMcase. This is a further indication of an enhanced turbu-lent activity in the EE case.Differences between the two cases have also been foundduring the linear stage of KHI. Growth rates of unstablemodes have different values in the EE and SM cases, alsoaccording to the relative vorticity-magnetic field orienta- roton kinetics in Kelvin-Helmholtz instability Figure 9.
Two-dimensional contour plots of the proton DF for the EE simulation at the time of the maximum level of turbulence(black dashed line in Fig. 5 ) and at the spatial point of the maximum of ǫ (black dot in Fig. 8). ( v y , v z ) plane is reported inthe left panel, ( v x , v z ) plane in the middle panel and ( v x , v y ) plane in the right panel. Figure 10.
Two-dimensional contour plots of the temperature anisotropy index A (left) and agyrotropy parameter √ Q (right)evaluated at the time of maximum activity for the EE simulation. vvv x yz B vvv xyz B v vv xyz B Figure 11.
Iso-surface plots of three proton DFs for the EE simulation at the same time of Fig. 9. From left to right the DFshave been evaluated at the spatial points corresponding to
A > A = 0, A <