Exact hybrid-Vlasov equilibria for sheared plasmas with in-plane and out-of-plane magnetic field
F. Malara O. Pezzi F. Valentini








Dipartimento di Fisica, Università della Calabria, 87036, Rende (CS), Italy 

 
The hybrid Vlasov-Maxwell system of equations is suitable to describe a magnetized plasma at scales of the order of or larger than proton 
kinetic scales. An exact stationary solution is presented by revisiting previous results with a uniform-density shear flow, directed either 
parallel or perpendicular to a uniform magnetic field, and by adapting the solution to the hybrid Vlasov-Maxwell model. A quantitative 
characterization of the equilibrium distribution function is provided by studying both analytically and numerically the temperature 
anisotropy and gyrotropy and the heat flux. In both cases, in the shear region, the velocity distribution significantly departs from local 
thermodynamical equilibrium. A comparison between the time behavior of the usual “fluid-like” equilibrium shifted Maxwellian and the exact 
stationary solutions is carried out by means of numerical simulations of the hybrid Vlasov-Maxwell equations. These hybrid equilibria can be 
employed as an unperturbed states for numerous problems which involve sheared flows, such as the wave propagation in inhomogeneous 
background and the onset of the Kelvin-Helmholtz instability.


Exact hybrid-Vlasov equilibria for sheared plasmas with in-plane and out-of-plane magnetic field
    F. Malara, O. Pezzi and F. Valentini
    December 30, 2023
================================================================================================





§ INTRODUCTION



Shearing flows in plasma are found in many natural contexts, like, for instance, the interface between the solar wind and planetary 
magnetospheres (e.g., <cit.>), the interaction region between fast and slow streams of the solar 
wind <cit.>, or astrophysical jets <cit.>. 

A magnetized plasma configuration with a shearing flow is stable if the jump Δ u of the bulk velocity u accross 
the shear is lower than a certain threshold, which is typically of the order of the background Alfvén velocity component parallel to 
u <cit.>. In this case, waves possibly propagating in the plasma are affected by the velocity shear 
through different effects. In particular, small scales can be progressively generated in the wave pattern in the direction of the bulk 
velocity gradient by a phase-mixing mechanism; this effect is similar to phase-mixing of Alfvén waves which propagate in a static 
background with an inhomogeneous Alfvén speed profile (e.g., <cit.>). The formation of structures at increasingly smaller scales can locally change the nature of waves, for instance converting an Alfvén wave into a Kinetic Alfvén wave 
<cit.>. Finally, the inhomonegeity associated with the shear couples propagating modes, generating energy 
transfers among different kind of waves <cit.>.

In the opposite case of sufficiently large Δ u (or small parallel magnetic field), a shearing flow in a plasma is unstable 
and undergoes the Kelvin-Helmholtz (KH) instability <cit.> which leads to the formation of vortices located at the shear layer. The dynamics 
can be more complex in configurations where one component of the magnetic field changes sign across the shear layer 
<cit.>, which generates an interplay between KH instability and magnetic reconnection. 
When the width of the shear layer becomes of the order of the ion inertial 
length, dispersive effects can modify the properties of the KH instability <cit.>. 
The presence of KH vortices has been revealed all along the flanks of the 
low-latitude Earth magnetopause <cit.> (see also <cit.> for a review of the KH 
instability in the magnetosphere context). Moreover, the development of a mixed KH-tearing mode instability has been considered to take 
place in cometary plasma tails <cit.>. 
A large Reynolds numbers, KH instability can become a source of turbulence through 
nonlinear interactions among vortices or secondary instabilities <cit.>. As an obervational example, we can mention a turbulent layer observed in the Earth's magnetopause <cit.>, 
where KH vortices forms in consequence of the solar wind motion relative to the magnetosphere. 

From a theoretical point of view, prior to the study of either wave evolution or KH instability in shearing-flow plasma, 
there is the problem of setting up a stationary configuration with a shear flow, where the above-mentioned phenomena can occur. 
This topic does not present any difficulty in 
the framework of the Magnetohydrodynamics (MHD), where a large variety of stationary flows can be envisaged, provided that a null total 
force acts on the magnetofluid at all positions. In contrast, the situation is totally different in the framework of a kinetic theory. In 
fact, building up a stationary particle distribution function (DF) representing a shear flow in the presence of a background magnetic 
field represents a non-trivial problem. Such an issue is particularly relevant in contexts where the width of the shear layer is of the 
order of typical kinetic scales, such as the ion inertial length or Larmor radius. Indeed, in these cases, kinetic aspects become relevant 
and can affect some properties of the stationary configuration, like profiles of the bulk flow, temperature, heat flux, etc., modifying 
them with respect to what can be deduced from fluid approaches. For instance, this is the typical situation which is encountered 
when one tries to model the shear flow at the Earth's magnetopause <cit.>.

The problem of deducing a stationary configuration with a shear flow for a magnetized plasma within a kinetic description has received a 
certain attention in the literature. In the case of uniform magnetic field, Cai et al. <cit.> have deduced a set of stationary ion 
distribution functions giving origin to 1D shear flows directed perpendicularly to the magnetic field. In particular, they have pointed out 
a different behaviour that is found according to the orientation of the flow vorticity with respect to the magnetic field. A similar topic 
has been treated by other authors: in the uniform B case, a form for a stationary shearing kinetic configuration has been deduced 
within the problem of calculating properties of ion-cyclotron modes <cit.>; the same problem has been considered 
in the non-uniform B case, deriving a set of specific profiles for bulk velocity and magnetic field <cit.>. 
More recently, a form of stationary distribution function has been considered in the case of a non-uniform parallel magnetic field and 
bulk velocity <cit.>. In the above approaches the fully kinetic problem has been considered, calculating the distribution 
functions of both ions and electrons. 

Despite of the existence of such results, in many studies of the KH instability a shifted local Maxwellian distribution has been 
considered as unperturbed plasma state (e.g., Ref. <cit.>). A shifted Maxwellian has the advantage to be more 
manageable because it allows to easily choose profiles of density, bulk velocity, temperature and magnetic field. On the other hand, at 
shear scales comparables with ion kinetic scales a shifted Maxwellian does not represent a stationary state: in fact, initializing the 
system with a shifted Maxwellian leads to the generation of oscillations <cit.>, whose amplitude becomes larger when decreasing the scale 
of the shear. Though this aspect could be considered as not relevant in the formation of the final turbulent state, it could somehow 
affect the linear stage of the instability. Moreover, a shifted Maxwellian does not describe the differences which arise in shears with a 
different vorticity orientation with respect to the magnetic field <cit.>, as happens in the dusk or dawn flanks of the 
magnetopause, where the vorticity to magnetic field orientation is reversed. 

More importantly, when studing wave propagation in shearing flows, the use of exact stationary states instead of a shifted Maxwellian is 
crucial. In fact, spurious oscillations associated with a non-exact stationary state would superpose to waves making difficult 
to single out effects due to waves.

When one is interested in describing phenomena at scales comparable with ion scales, a 
successful approach is represented by the so-called hybrid Vlasov-Maxwell (HVM) model, in which ions are 
kinetically described while electrons are treated as a massless fluid <cit.>. It represents a sort of compromise between the 
usual “MHD-like” coarse-grained fluid description and an exceedingly complex fully kinetic approach. The HVM model has been 
adopted for describing several phenomena occurring at scales where the kinetic ion physics starts to play a significant role into the 
plasma dynamics <cit.>. 
Within the HVM framework, a method to derive ion distribution functions representing approximated kinetic stationary solutions 
has been recently presented by Cerri et al. <cit.>. It is based on the “extended fluid model approach”, which considers finite 
Larmor radius effects in the determination of the ion pressure tensor, and it has been recently adopted to describe temperature 
anisotropies due to the shear velocity <cit.>. However, it should be noticed that the solutions proposed by Cerri et al. 
<cit.> are not exactly stationary, even in the framework of the HVM theory. In fact, oscillations are still present, even if with 
amplitudes definitely smaller than those recovered in the shifted Maxwellian case. 

In the present paper we derive exact stationary solutions for the set of the HVM equations describing a 
magnetized plasma with an arbitrary shearing flow u profile in two different configurations, namely, with a uniform magnetic 
field B parallel or perpendicular to u. The derivation of our solutions is inspired by full-kinetic 
solutions previously obtained in analogous configurations <cit.>, which have been adapted to 
the HVM model. In particular, an explicit analytical expression for the solution is found in the parallel B case, while 
in the case of perpendicular B the solution is calculated by a numerical procedure which integrates single-particle trajectories.
The interest of these solution is twofold: first, they are exactly stationary, thus can be safely used as unperturbed states either in 
wave propagation models and in instability studies; second, due to the properties of the HVM model, they realistically represents 
situations in which the width of shear layer is of the order of or larger than ion kinetic scales, avoiding the complexity of a 
fully kinetic treatment.

The plan of the paper is the following: in Section II we introduce the equations of the HVM model; in Sections III and IV we derive 
and discuss the stationary solution in the cases of parallel and perpendicular magnetic field, respectively; in Section V we present the results of numerical simulations where we analyze the behaviour of the solution in comparison with that of a shifted Maxwellian; finally, in Section VI we summarize the results. 
















§ EQUATIONS OF THE MODEL

We consider a quasi-neutral magnetized plasma composed of kinetic protons and a massless fluid of isothermal electrons (the current 
analysis easily extends to heavier ions). We are interested in describing shears occurring at scales larger or comparable with proton 
kinetic 
scales: i.e. l≳ d_p ∼ρ__Lp and τ≳Ω_cp^-1, being d_p=V__A/Ω_cp the proton skin depth;
ρ__Lp=v_th,p/Ω_cp the proton Larmor radius; Ω_cp=e B_0/m_p c the proton cyclotron frequency; 
V__A=B_0/√(4π n_0m_p) the Alfvén speed; v_th,p=(k_B T_0/m_p)^1/2 the proton thermal speed. Note that β_p=2 
v_th,p^2/V__A^2; m_p, e, n_0 and T_0 are respectively the proton mass, charge, density and temperature; while B_0, c and 
k_B are the magnetic field typical value, the light speed and the Boltzmann constant. 

At these scales, the proton dynamics is successfully modeled by the hybrid Vlasov-Maxwell (HVM) equations: 

    ∂ f/∂ t +  v·∇ f + e/m_p(  E + v/c× B) ·∂ f/∂ v =0 
          ∂ B/∂ t= -c ∇× E    ;     j=c/4π∇× B
           E=-1/c(  u× B)+ 1/en(  j× B/c) -
    1/en∇ p_e

where f=f( x, v,t) is the proton distribution function in the phase space (  x,  v). The electric E( 
x,t) and magnetic B( x,t) fields are respectively determined by the generalized Ohm's law and by the Faraday and Ampere laws, 
by neglecting the displacement current. In Eqs. (<ref>–<ref>), n( x,t) = ∫ d^3 v f( x, v,t) is the 
proton number density, u ( x,t) = ∫ d^3 v  v f( x, v,t)/n( x,t) is the proton bulk speed 
and j( x,t) is the current density. Electrons are a massless fluid, whose density is equal to that of ions n_e=n for the 
quasi-neutrality condition. The electron pressure p_e is determined by imposing a closure assumption for the electron dynamics (such as 
isothermal or adiabatic state equation). 

Our goal is to build up a stationary hybrid equilibrium state for a sheared flow in presence of a homogeneous background magnetic field. The 
shear is directed along y and depends on x, i.e. u_0=u_0(x) e_y, and spatial variations occur only along x. Two different 
cases are discussed: i) magnetic field B parallel case to u_0, i.e. B=B_0 e_y; ii) magnetic field 
B perpendicular to u_0, i.e. B=B_0 e_z. These two cases will be investigated separately in the 
following of the paper. 

In the parallel case, no electric field is needed to set the initial equilibrium and electrons are a massless isothermal fluid: T_e= 
const→ p_e = k_B n T_e. In the perpendicular case, an equilibrium electric field, 
directed along x, is needed to equilibrate the u_0× B term in Eq. (<ref>): E= E(x)  e_x. In this case we 
also need to relax the electrons closure, by treating the electron pressure p_e as a further independent quantity determined by the 
following equation:

    [ ∂/∂ t + ( u_e ·∇) ] ( p_e/n^γ_e) = 0

where γ_e is the electron adiabatic index. Last equation implies that the electron temperature is not homogeneous. Therefore:

    E(x) = - u_0 B_0/c - 1/n ed p_e/dx

The presence of such electric field introduces a charge separation, which - in principle - is not taken into account within the HVM model 
(n_e=n). However, the discrepancy from the quasi-neutrality condition is extremely small.

To conclude this section we remark that the first attempt for modeling a plasma configuration with a velocity shear at kinetic scales 
essentially extends the fluid, large scales, equilibrium at smaller scales (e.g., Ref. <cit.>). 
These equilibria, based on the local thermodynamical hypothesis <cit.>, assume the following form for the proton VDF 
(sheared Maxwellian):

    f_SM(x, v,t) =  n_0/(2π)^3/2 v_th,p^3exp[ - v_x^2 + ( v_y -u_0(x))^2 + v_z^2/2 
    v_th,p^2]

where v=(v_x, v_y, v_z), n_0 is the proton density and u_0=u_0(x)  e_y is the bulk velocity, being u_0(x) a given 
function.  It can be easily verified that f_SM is not a stationary solution of HVM 
equations in both cases (parallel and perpendicular) discussed above.

In the following of the paper, we will revisit and formalize, for the HVM equations, the derivation of kinetic stationary equilibria
for a sheared flow, by considering two different geometrical configurations. We anticipate that our approach is easier with respect to the 
ones adopted in previous works <cit.> and implementing our equilibrium in the HVM 
code is hence quite simple. This allows 
to perform hybrid kinetic simulations with a “proper” (and simple) hybrid kinetic equilibrium. We also remark that, the sheared Maxwellian 
DF is often adopted also for analyzing phenomena occurring at kinetic scales in sheared flows, such as KHI. This choice can be easily 
justified for the investigation of a particular class of phenomena, where nonlinearities play a crucial role in developing turbulence, such 
as KHI (e.g., <cit.>). One can correctly argue that, for these processes, the final difference that would occur starting with the 
sheared Maxwellian or with the “proper” equilibrium DF is minimal and does not affect the final dynamical state. However, for another 
class of processes (such as phase mixing of a linear wave in a velocity shear <cit.>), starting from a 
correct equilibrium is crucial for correctly investigating the phenomenon itself. In this direction, several works have been focused on the 
extension the MHD-like fluid equilibrium to the microphysics <cit.> and our work gives a 
further contribution in this direction.



§ STATIONARY SOLUTION FOR THE PARALLEL CASE: U_0  ||  B



Here we revisit the derivation of the stationary solution of Eqs. (<ref>–<ref>), in the case in which velocity shear and 
homogeneous magnetic field are parallel and the electric field is vanishing. The method we adopt is similar to the Harris approach for 
deriving kinetic equilibria corresponding to localized current sheets <cit.> and it is based on the determination of the proton DF 
as a function of the motion constants, derived from the Lagrangian description of single-particle dynamics. A similar derivation can be 
found in Ref. <cit.> in the framework of fully kinetic theory; here, we reconsider the same problem to adapt the solution 
to the HVM model. Moreover, we discuss in a deeper detail the properties of the derived stationary DF, illustrating the spatial
profiles of its moments (bulk velocity, temperatures, heat flux).

For the present geometry, the Lagrangian of the single-particle dynamics is:

    ℒ = m_p/2( v_x^2 + v_y^2 + v_z^2) - e/m_p v_zB_0 x

where A=A_z  e_z = B_0 x  e_z is the vector potential associated with the magnetic field. From the integrals of 
motion of ℒ, which are the generalized momenta P_y, P_z and the energy E, we can define three auxiliary constants:

    k_1 = P_y =m_p v_y 
       k_2 = - P_z c/e B_0 = x - v_z/Ω_cp
       k_3 = E = m_p/2( v_x^2 + v_y^2 + v_z^2)

To write the proton DF f_eq,||, we consider the following combinations of the above constants, having the dimension of a velocity, α_1 = {2/m_p[ E - P_y^2/2 m_p] }^1/2 and α_2 = P_y/m_p - U(x - 
v_z/Ω_cp), being U=U(k_2) the arbitrary shear, function of the motion integral k_2. 

Since α_1 and α_2 are motion integrals, each generic function F(α_1^2 + α_2^2) is a stationary 
solution of the Vlasov equation <cit.>. The function F is determined by imposing that f_eq,|| reduces to a Maxwellian with 
density n_0 and thermal speed v_th,p in absence of the shear (U=0), hence:

    f_eq,||(x, v) = n_0 /(2π)^3/2 v_th,p^3exp{-1/2v_th,p^2[v_x^2 + (v_y-U(x - 
    v_z/Ω_cp))^2 + v_z^2 ]}

This solution is a stationary equilibrium for the HVM set of equations Eqs. (<ref>–<ref>), for the case of a shear 
parallel to the background magnetic field. Indeed Eqs. (<ref>–<ref>) are also exactly satisfied, since the 0-th order moment 
of f_eq,|| gives a homogeneous density n_0 and the bulk speed u is along y (see the next subsection for further details). 
The DF f_eq,|| is a Maxwellian-like function shifted along v_y, in which however the amount of the shift depends on both 
the position x and the velocity v_z (through the argument of U, k_2=x -v_z/Ω_cp). For the same reason, the bulk speed u_y 
does not coincide with U in the general case. 



 §.§ Properties of the stationary distribution function f_eq,||

To characterize the physical properties of f_eq,||, we focus here on its shape in velocity space and on the evaluation of its moments, 
in particular density, bulk velocity, temperature and heat flux.

To display the shape of f_eq,|| in the velocity space, we choose the tanh-like shear profile of U, routinely adopted for 
investigating the KHI instability:

    U(k_2) = U_0tanh(k_2/Δ x)

where we choose U_0=2V__A, β_p=2v_th,p^2/V__A^2=4 and Δ x is the width of the shear function U. The width of the 
sheared bulk-velocity u(x) can be different from Δ x. By using the latter expression of U, we compute f_eq,|| in Eq. 
(<ref>), by discretizing the four dimensional phase space through N_x=512 grid points in the one-dimensional spatial domain 
(x∈[-L/2,L/2]) and N_v=141 grid points in each velocity direction (v_j∈[-v_max,v_max], being j=x,y,z and 
v_max=7v_th,p), while we chose L=50 d_p and Δ x = 2.5 d_p. Figure <ref> show the iso-surface plots of f_eq,|| in 
velocity space at several spatial positions across the shear: x/d_p=-2.34 (a), x/d_p=0.0 (b) and x/d_p=2.34 (c). The red tube in Fig. 
<ref> indicates the magnetic field direction. Far from the shear (not explicitly reported) the distribution function is a shifted 
Maxwellian, while across the shear it becomes significantly stressed, resembling potato-like shapes with non-null heat flux, temperature 
anisotropy and agyrotropy with respect to its principal axes.

We also calculate the moments of the DF f_eq,||. We already anticipated that the density associated with the DF given by Eq. (<ref>) 
is uniform, i.e. n(x) = n_0. The bulk velocity can be easily evaluated, starting from u ( x) = ∫ d^3 v  
v f_eq,||( x, v)/n_0. It can be easily shown that u_x and u_z are null, while u_y is:

    u_y(x)= 1/(2π)^1/2 v_th,p∫_-∞^∞ U( x - v/Ω_cp) 
    exp( -v^2/2 v_th,p^2)dv

The bulk velocity u_y(x) associated with f_eq,|| does not coincide with the function U(x), rather it is the result of the 
convolution between U and a Gaussian function. 

Despite the inverse process - i.e. the determination of U(x) and of f_eq,|| for a given profile of the bulk velocity u_y(x) 
- requires the inversion of the convolution in Eq. (<ref>), it is still possible to deduce some simple results. If Δ x is 
the characteristic spatial length of U(x), then U(x-v/Ω_cp) - considered as a function of v - varies over the scale Δ v = 
Ω_cpΔ x. The Gaussian factor inside the integral of Eq. (<ref>) represents a windowing function selecting a v 
interval of width ∼ v_th,p; we indicate such a windowing function by:

    W_G(v)=1/(2π)^1/2 v_th,pexp( -v^2/2 v_th,p^2)


a)
In the large scale limit Δ x ≫ρ_Lp, i.e. Δ v ≫ v_th,p, the profile of the windowing function is relatively 
unimportant and W_G(v) can be successfully approximated with the simpler square window W_S(v), centered in v=0, with width and 
amplitude √(2π)v_th,p and 1/√(2π) v_th,p, respectively. W_G(v) and W_S(v) have the same value at v=0 
and the same integral in the interval ( -∞ , +∞ ). Within the approximation W_G(v) ≃ W_S(v), we obtain:

    u_y (x) ≃U̅(x) = 1/√(2π)ρ__Lp∫_-√(π/2)ρ__Lp^√(π/2)ρ__Lp 
    U(x-ξ)dξ

which represents the running average of U(x) performed over the interval x-√(π/2)ρ__Lp,x+√(π/2)ρ__Lp. In the large scale limit, the bulk velocity u_y(x) is approximately given by the function U(x) smoothed over 
an interval of amplitude √(2π)ρ__Lp centered at the position x. This result can be easily understood by thinking that the 
protons gyromotion in the plane perpendicular to B mixes up the v_y velocities of the single protons, thus smoothing the profile 
of U(x) over a length scale which is of the order of the Larmor radius. 

b)
In the opposite small scale limit, , i.e. Δ x ≪ρ__Lp,  Eq. (<ref>) can be rewritten as follows:

    u_y(x)= 1/(2π)^1/2∫_-∞^∞ U(ρ__Lpφ) exp[ -1/2(φ 
    -x/ρ__Lp)^2 
    ]dφ

where φ=k_2/ρ__Lp and U(ρ__Lpφ) as a function of φ varies on a scale much smaller than unity. If 
U(ρ__Lpφ) describes a shear layer corresponding to a bulk velocity which varies in the range [-U_0,U_0], it can be 
approximated with the Heavyside function H(φ):

    U(ρ__Lpφ) ≃ U_0 H(φ)=
    
        U_0        if φ>0 
    
       -U_0       if φ<0

After some algebraic steps, Eq. (<ref>) reduces to

    u_y(x)=
    
        U_0        if  x≫ρ__Lp
    
       -U_0       if  x≪ -ρ__Lp    ,     for Δ x ≪ρ__Lp

where, for simplifying the integrals, we considered that, for x≫ρ__Lp (x≪ -ρ__Lp), the Gaussian is essentially 
located in the positive (negative) part of the φ axis, where H(φ)=1 (H(φ)=-1). Hence, in spite of the small scale 
of variation of the function U(x) (Δ x ≪ρ__Lp), the bulk velocity u_y(x) varies on a scale comparable with the 
proton Larmor radius ρ__Lp. It is not possible to construct, in the current configuration, shear layers with a width smaller than 
the proton Larmor radius. This is again due to the proton gyromotion which mixes up the parallel velocity of single particles on a 
transverse scale of the order of ρ__Lp.

To directly display the shape of the bulk velocity, we consider the shear function U(x) [Eq. (<ref>)] and we numerically compute 
u_y(x), for Δ x = 2.5 d_p  ≃ 1.77 ρ__Lp. Figure <ref> reports, the spatial profile of the function U(x) 
(black solid line) and the corresponding bulk velocity u_y (red dashed line). It is clear to note that significant differences between U 
and u_y are recovered. We also verified that in the large scale limit the windowing function does not play a significant role and hence, 
u_y(x) ≃ U(x), while in the small scale limit protons arrange themselves to produce a bulk velocity u_y(x) varying over a 
scale comparable with the proton Larmor radius, nevertheless the shear function U varies over scales much smaller than ρ__Lp 
(not reported here).

To further characterize the moments of f_eq,||, we consider the variance matrix, defined by:

    σ_ij(x) = 1/n(x)∫[ v_i - u_i(x) ] 
    [ v_j - u_j(x) ] f_eq,||(x, v) d^3  v   ;    i,j = x,y,z

which is related to the proton temperature by T_0 = m_p∑_j=1^3 σ_jj/3k_B. Since the magnetic field B is uniform and 
directed along y, the proton temperatures parallel and perpendicular to B, i.e. in the local B frame (LBF), are defined, 
respectively, by T_||=m_p/k_Bσ_yy and T_⊥=m_p/k_B (σ_xx+σ_zz)/2, so that T_0 = 
(T_||+2T_⊥)/3.

The analytical evaluation of anisotropy and agyrotropy at the center of a symmetric shear, presented in App. <ref>, indicates 
that the equilibrium DF is anisotropic and agyrotropic. By means of the numerical evaluations of the variance matrix elements, we 
can extend the analytical computation and consider not only the center of the shear. We numerically diagonalize the matrix σ, thus 
rotating the DF into the minimum variance frame (MVF). The eigenvalues of σ, corresponding to the temperatures in the MVF, are 
indicated by: λ^(3) < λ^(2) < λ^(1). 

The top panel of Fig. <ref> shows the temperature anisotropy ratio η and η^*, while the bottom panel indicates the 
agyrotropy ratio ζ and ζ^*, for the shear U(x) given by Eq. (<ref>) in the case Δ x=2.5d_p, both in the MVF 
(black solid) and in the LBF (red dashed). Temperature anisotropy and agyrotropy have been evaluated as follows: (a) temperature anisotropy 
in the MVF η=(λ_2+λ_3)/2λ_1; (b) temperature anisotropy in the LBF η^*=(σ_xx+σ_zz)/2σ_yy; 
(c) agyrotropy in the MVF ζ =λ_3/λ_2; (d) agyrotropy in the LBF ζ^* 
=min(σ_xx,σ_zz)/max(σ_xx,σ_zz). In the LBF, the DF is anisotropic close to the shear, while no temperature 
agyrotropies are recovered. On the other hand, in the MVF, the DF displays strong anisotropies as well as agyrotropies close to the velocity 
shear.

We finally characterize the DF f_eq,|| by computing the heat flux for unit of mass:

    q_j(x) = 1/2∫[v_j - u_j(x] ) [  v -  u(x) ]^2 f_eq,||(x, v) d^3  v   ;    
    j= x,y,z

where the shear function U(x) in Eq. (<ref>), with Δ x= 2.5 d_p, is adopted. Figure <ref> reports the three 
components of the heat flux q_x (black-solid line), q_y (red-dashed line) and q_z (blue-dotted line), as a function of x/d_p. The 
equilibrium DF f_eq,|| is such that a non-vanishing heat flux is recovered at x/d_p≃ 0 in the two components q_y and q_z, 
which tends to zero away from the shear.



§ STATIONARY SOLUTION FOR THE PERPENDICULAR CASE: U_0 ⊥  B



We revisit here the derivation of the stationary solution of Eqs. (<ref>–<ref>), in the case in which velocity shear and 
homogeneous magnetic field are perpendicular, while the electric field is E = E(x)  e_x. The method here adopted is based 
on the dynamics of a single proton in the electric E and magnetic B field. Particle trajectories have been already studied 
to build up a stationary solution in a fully kinetic (ion+electron) description in previous studies <cit.>. 
In particular, in the analytical description of single particle dynamics we follow the same method as Ganguli et al. <cit.>, but 
deriving further general results which are important for setting up our numerical description of particle dynamics. We will also derive a 
form for the proton distribution function which is different from that in Ref. <cit.> (except in the particular case of linearly 
growing electric field). In our case the derived DF is furthermore supplemented with a form for the electron pressure profile, which allows 
us to obtain an exact stationary state for the whole set of HVM equations.

Our derivation starts from considering the single-particle motion. The proton motion along the parallel z direction is decoupled from 
the motion in the transverse plane. Therefore, we focus on the particle motion in the xy plane, described by:

    dv_x/dt   =   Ω_cp v_y + e/m_p E(x) 
    dv_y/dt   =    -Ω_cp v_x 
    dx/dt   =    v_x

The particle motion depends on the electric field profile E(x), which indirectly determines the profile u(x) of the bulk velocity. We 
consider a situation where u(x) varies crossing one or more shear layers, but becomes essentially uniform far from the shear layers: 
E(x)= E_+∞ (E(x)=E_-∞) for large positive (negative) x. Thus, in the homogeneous region particles drift along y with a 
uniform drift velocity v_d,±∞ = c( E_±∞× B)/B^2.

From Eqs. (<ref>), it is easy to obtain:

    v_y    =    -Ω_cp x + W_0  
    d^2 x/dt^2   =    -Ω_cp^2 x + e/m_p E(x) + Ω_cp W_0

where W_0 is a constant determined by initial conditions and Eq. (<ref>) is a nonlinear oscillator equation for x(t). We integrate 
Eq. (<ref>) in the interval [ x_0,x], by considering that dv_x/dt = d/dx(v_x^2/2) and by rewriting 
v_0y=W_0-Ω_cp x_0 (being x_0 an arbitrary position corresponding to v_0y):

    1/2 m_p v_x^2 + e ϕ(x) + 1/2 m_p Ω_cp^2 (x-x_0)^2 -m_pΩ_cp v_0y (x-x_0) = e_0

where ϕ(x)= -∫_x_0^x E(x') dx' is the electrostatic potential, which vanishes at x=x_0, and e_0=m_p v_x0^2/2 is a constant, 
v_x0 being the value of v_x at x=x_0. 

Eq. (<ref>) expresses the energy conservation for a particle with mass m_p following a 1D motion in the effective potential energy 
U_ eff(x) = U_E(x) + U_Ω(x), where U_E(x)=eϕ(x) is the electrostatic potential energy and 

    U_Ω(x) = 1/2 m_p Ω_cp^2 (x-x_0)^2 -m_pΩ_cp v_0y (x-x_0) =
    1/2 m_p Ω_cp^2 ( x- W_0/Ω_cp)^2 + 1/2 m_p v_0y^2

The term U_Ω(x) hence corresponds to an attractive force towards the position W_0/Ω_cp.

For large |x-x_0| the term U_Ω(x) dominates in the determining the effective potential energy U_ eff(x). Indeed, since 
E(x) becomes constant for sufficiently large values of |x-x_0|, we have that U_E ≃ -E_+∞(x-x_0) [U_E ≃ 
-E_+∞(x-x_0)] for large positive (negative) values of x-x_0, while U_Ω(x) is quadratic with respect to (x-x_0). The 
motion of the particle along x is hence confined inside a potential well: x_m ≤ x(t) ≤ x_M, where U_ eff(x_m)=U_ 
eff(x_M)=e_0 and the particle moves back and forth in the interval [ x_m , x_M ], with vanishing v_x at x_m and x_M. 
In other words, x(t) and v_x(t) are periodic function with period τ and, from Eq. (<ref>), follows that v_y(t) is also 
periodic with period τ. Therefore, the particle follows a closed trajectory in the v_xv_y plane. Notice that y(t) is not 
necessarily a periodic function and the particle trajectory in the xy plane is, in general, an open curve. The details of the motion along 
x, like the period τ, depend both on the specific form of the electric field profile E(x) and on the particle initial conditions, 
which determine the constant quantity W_0. However, the periodicity of variables x(t), v_x(t) and v_y(t) holds for any form of 
E(x) and for any initial condition. This result is crucial for the setup of the numerical method we employed to calculate a stationary 
proton DF f_eq,⊥ for an arbitrary electric field profile.

Since the particle motion in the v_xv_y plane is periodic, the time average of the velocity over the period τ provides the drift 
velocity in the particle motion. Therefore, we define the guiding center velocity v_c as: 

    v_c = ⟨ v⟩_τ = 1/τ∫_0^τ v(t)  dt

The x-component of the guiding center velocity is trivially vanishing, hence: v_c=v_cy e_y = ⟨ v_y ⟩_τ  
e_y. We also define the guiding center position x_c as the position where the particle velocity component v_y is equal to the 
guiding center velocity: v_y=v_cy; hence, from Eq. (<ref>),  we find

    x_c = (W_0-v_cy)/Ω_cp

Note that: (i) equation (<ref>) implies that a single value is admitted for x_c; (ii) the guiding center position x_c 
represents also the time-averaged particle x-position: x_c=⟨ x ⟩_τ.

We point out that our particle guiding center definition is different from that used in previous studies. In fact, in Refs. <cit.> the 
guiding center position is defined as a point where v_y=u_y, implying that a given particle can have more than one guiding center (see 
Ref. <cit.> for a discussion). In contrast, in our approach a single guiding center is defined for each particle, regardless of 
the specific electric field profile E(x) and of the particle initial condition. In Ref. <cit.> the guiding center position is 
defined as the position where the effective potential energy is minimum. Also this definition differs from ours, except for particular 
profiles of the linearly growing electric field (see App. <ref>).

To build the stationary DF for the HVM Eqs. (<ref>)-(<ref>), we consider the particle total energy, which is a 
constant of motion:

    ℰ=K+U'_E = 1/2m_p ( v_x^2 + v_y^2 + v_z^2 ) + U'_E

where the electric potential energy U'_E(x) is re-defined such that U'_E(x_c)=0, i.e. it vanishes at the guiding center position x_c 
of the considered particle U'_E(x)=-e ∫_x_c^x E(x')  dx'. This choice can be justified by the following argument. Let us consider 
the particular case of a uniform electric field E(x)=E_0, corresponding to proton circular orbits in the v_xv_y plane, with uniform 
drift velocity v_cy=-cE_0/B. In such a case, the potential energy has the form U'_E(x)=-e E_0 (x-x_c). Since x_c=⟨ x 
⟩_τ, it follows that ⟨ U'_E ⟩_τ=0, i.e. it has the same value for all particles, regardless of their average 
position x_c. This is in accordance with the macroscopic invariance of fluid properties with x, which characterizes this particular 
case. In contrast, a potential energy U_E which vanishes at a fixed position x_0 (equal for all the particles) would give ⟨ U_E 
⟩_τ=eE_0(x_c-x_0), i.e., an average potential energy which systematically varies with the average position x_c of particles.

We define also the quantity ℰ_0=ℰ-m_p v_cy^2/2 which represents the part of the particle energy not due to the 
drifting motion. 
Of course, ℰ_0 is another motion constant. Since ℰ is constant, its value can be evaluated at x=x_c, where 
U'_E(x_c)=0 and 
where, by definition, v_y=v_cy. Thus:

    ℰ=1/2m_p {[v_x(x=x_c)]^2 + v_z^2 } + 1/2m_p v_cy^2

implying 

    ℰ_0=1/2m_p {[v_x(x=x_c)]^2 + v_z^2 }


To understand how to define the proton DF, we again consider the particular case of constant electric field and we require that, in this 
case, the DF is a shifted Maxwellian:

    f_SM( v) = C exp[ -v_x^2+(v_y-v_cy)^2+v_z^2/2 v_th,p^2]

where v_th,p is the thermal speed and C is a normalization constant. The uniform-E case above analyzed can be also reproduced in 
the local approximation description, whose details are reported in App. <ref>, by setting α_0=0 
(ω^2=Ω_cp^2), which implies that v_x^2+(v_y-v_cy)^2=[ v_x(x=x_c)]^2. Therefore, in this case, by comparing 
Eq. (<ref>) and Eq. (<ref>), we conclude that:

    f_SM( v) = C exp( -ℰ_0/m_p v_th,p^2)

The above considerations lead us to the following “ansatz": we hypothesize that a stationary proton DF representing a shearing flow for 
any electric field profile E(x) can have the following implicit form:

    f_eq,⊥(x, v)=C exp[ -ℰ_0(x,v_x,v_y,v_z)/m_p v_th,p^2]

where ℰ_0(x,v_x,v_y,v_z)=ℰ-m_p v_cy^2/2 is the single-proton energy not due to the drifting motion. As showed above, 
Eq. (<ref>) is a shifted Maxwellian for an uniform electric field. 

Of course, such a conjecture must be verified a posteriori. This can be done first in the particular case of the local approximation, 
reported in App. <ref>, where we also derive the explicit form of the equilibrium DF. Then, in the case of a 
generic electric field profile E(x) the same quantities will be calculated by employing a numerical technique. It is important to 
highlight that, since the quantity ℰ_0 is a constant of the particle motion, f_eq,⊥ is a stationary solution of the 
Vlasov equation [Eq. (<ref>)], provided that both the electric and magnetic fields are temporally constant <cit.>.

The numerical method employed to generate the equilibrium DF f_eq,⊥, for a generic electric field profile E(x), is described 
in the following. We assume that f_eq,⊥ has the form given by Eq. (<ref>), being 

    ℰ_0(x,v_x,v_y,v_z)=m_p/2( v_x^2+v_y^2+v_z^2 ) +U'_E(x,v_x,v_y)-m_p/2 v_cy^2(x,v_x,v_y)

Here, v_cy^2(x,v_x,v_y) indicates the guiding center velocity of a particle which is located at the position x, with velocity 
(v_x,v_y) at a given time t, while U'_E(x,v_x,v_y) is the electric potential of the same particle. Since ℰ_0 is a motion 
constant, the time t when ℰ_0 is evaluated can be arbitrarily chosen. Hence, in Eq. (<ref>), (x,v_x,v_y,v_z) can be 
interpreted as the position and velocity of a single particle at the initial time of its motion. The evaluation of the last two terms in 
Eq. (<ref>) requires however to integrate the single particle motion Eqs. (<ref>). This has been done by the following 
numerical procedure:

(i) We consider a 1D-3V phase space, composed by a spatial coordinate x∈[0,L], discretized with N_x grid points, and three 
velocity coordinates v_m ∈[ -v_max,v_max]   m=x,y,z, discretized with N_v grid points along each direction. We 
numerically integrate N_x× N_v^2 particle motion Eqs. (<ref>), using each point (x_i,v_x,j,v_y,k) of the subgrid as 
initial condition: x(t=0)=x_i, v_x(t=0)=v_x,j, v_y(t=0)=v_y,k (i,j,k are indexes which span along x, v_x and v_y, 
respectively). The v_z component, whose index is l; is neglected since the motion is trivial along z. The time integration of Eqs. 
(<ref>) has been carried out through a 3-order Adam-Bashforth scheme, being the time step Δ t chosen to maintain the CFL 
condition.

(ii) Since the exact trajectories in the v_x v_y plane are necessarily closed, each integration is stopped when the corresponding 
orbit in the v_x v_y plane is completed. The corresponding time T represents the orbit period. 

(iii) We calculate, for each orbit, the velocity and x-position of the guiding center: v_cy,ijk=⟨ v_y ⟩_T and 
x_c,ijk=⟨ x ⟩_T. The electric potential associated with the particle initial position is calculated as 
ϕ_ijk=-∫_x_c,ijk^x_i E(x)  dx, where the integral is numerically evaluated. Then, the DF at a given point of the phase 
space is evaluated as

    f_eq,⊥(x_i,v_x,j,v_y,k,v_z,l) = C exp[ -1/2 v_th,p^2( v_x,j^2 + v_y,k^2 + v_z,l^2 - v_yc,ijk^2 
    ) 
    -e ϕ_ijk/m_p v_th,p^2]

From this expression, the moments of the distribution function (density, temperatures, bulk velocity and heat flux components) are 
numerically evaluated. In particular, all the moments vary only in the x direction and the bulk velocity u(x) is directed in the 
y direction.

(iv) The resulting u_y(x) depends on the chosen profile of the electric field E(x). However, for an arbitrary electric field profile 
E(x), the bulk velocity does not coincides with the local E× B drift velocity. On the other hand, 
E and u appear in the generalized Ohm's law (<ref>), which must be consistently satisfied. This is obtained by 
choosing a profile for the electron pressure p_e(x) such that

    dp_e/dx=-qn(x) [ Bu_y(x)/c+E(x) ]

Then, in the general case the electron pressure p_e is not uniform, except in the case of a linear electric field profile (See App. 
<ref>). By adopting this closure for the electron pressure, it is easy to show that the considered configuration is a stationary 
solution of the whole set of HVM Eqs. (<ref>)-(<ref>) supplemented by Eq. (<ref>) for the electron pressure. 




 §.§ Properties of the stationary distribution function f_eq,⊥


To analyze the properties of the equilibrium DF f_eq,⊥, we focus on the following shape for the electric field: 

    E_x(x) = - E_0 tanh( x/Δ x)

representing a shear layers of amplitude E_0=1 (in scaled unit) and width Δ x=2.5 d_p. By using the latter expression of E_0,x, 
we compute f_eq,⊥ in Eq. (<ref>), by discretizing the four dimensional phase space through N_x=512 grid points in the 
one-dimensional spatial domain (x∈[-L/2,L/2]) and N_v=141 grid points in each velocity direction (v_j∈[-v_max,v_max], being 
j=x,y,z and v_max=7v_th,p), while we chose L=50 d_p and Δ x = 2.5 d_p. It is worth to note that the difference between 
the electric field from which we compute the equilibrium [Eq. (<ref>)] and the term  - u× B/c evaluated using the 
equilibrium DF mean speed is of the order of 10^-5. Despite this quantity is small at the initial time, one needs to take 
care of it by self-consistently solving Eq. (<ref>) to maintain the equilibrium.

Figure <ref> reports iso-contour of the proton DF f_eq,⊥ in the velocity space. Panels (a) to (c) refer to 
different positions across the shear: x/d_p=-2.34 (a), x/d_p=0.0 (b) and x/d_p=2.34 (c). The red tube in Fig. <ref> 
indicates the magnetic field direction. We note that, against the parallel case, here the equilibrium DF is less stressed and exhibit a 
bi-Maxwellian like structure, elongated in a direction transverse to the magnetic field direction. As for the parallel case, far from the 
shear, the DF f_eq,⊥ reduces to the shifted Maxwellian, while - across the shear - it exhibits a clear temperature anisotropy.


Figure <ref> reports the temperature anisotropy (top panel) and agyrotropy (bottom panel) ratios both in the MVF (black 
solid) and in the LBF (red dashed). Temperature anisotropy and agyrotropy have been evaluated as follows: (a) temperature anisotropy in the 
MVF η=(λ_2+λ_3)/2λ_1; (b) temperature anisotropy in the LBF η^*=(σ_xx+σ_yy)/2σ_zz; (c) 
agyrotropy in the MVF ζ =λ_3/λ_2; (d) agyrotropy in the LBF ζ^* 
=min(σ_xx,σ_yy)/max(σ_xx,σ_yy). Note that the definitions in the LBF are different from the parallel case, 
since the orientation of the magnetic field has been changed. If in the parallel case the equilibrium distribution function was 
characterized by stronger anisotropies in the MVF frame, here the situation is opposite. Indeed, in the MVF, the DF is strongly anisotropic 
at the shear, but it does not present significant non-gyrotropic features. On the other hand, in the LBF frame, the DF is significantly 
anisotropic as well agyrotropic.

Finally, we characterize the DF by evaluating the heat flux [Eq. (<ref>)]. Figure <ref> reports the three components of the 
heat flux q_x (black-solid line), q_y (red-dashed line) and q_z (blue-dotted line), as a function of x/d_p. Clearly, a non-vanishing 
heat flux is recovered at x/d_p ≃ 0 in the y component.



§ HYBRID VLASOV-MAXWELL SIMULATIONS OF THE EQUILIBRIUM


In this section we numerically test that distribution functions derived in the previous sections f_eq,|| and f_eq,⊥ are 
effectively stationary equilibria for the HVM set of equation, which are solved numerically in the so-called HVM code <cit.>, by also comparing these equilibria with the sheared Maxwellian case f_SM. 

The HVM code solves numerically the set of Eqs. (<ref>–<ref>) in dimensionless form through a Eulerian algorithm described 
in detail in Ref. <cit.>. In the parallel case, since T_e is homogeneous and constant, Eq. (<ref>) is trivial. 
Dimensionless HVM equations are obtained by scaling velocities by the Alfvén speed V__A, lengths by the proton skin depth d_p and 
time by the inverse proton cyclotron frequency Ω_cp^-1. Since the problem is intrinsically one-dimensional in physical space, we 
restrict our numerical runs to a phase space of reduced dimensionality (1D in physical space and 3D in velocity space). The code assumes 
periodic boundary conditions in the spatial coordinate x∈[0,L], while the DF f(x, v,t) is set equal to zero for |v|>v_max in 
each velocity direction and at each spatial position. 



 §.§ Parallel case

For this case, we discretized the four dimensional phase space through N_x=512 grid points in the one-dimensional spatial domain and 
N_v=141 grid points in each velocity direction, while v_max=7 v_th,p. 

We performed two simulations (S1 and S2) keeping fixed the background magnetic field B_0=B_0  e_y (B_0=1 in scaled units) and 
the proton plasma beta β_p=2v_th,p^2/V__A^2=4, but changing the initial proton distribution function. In both simulations the 
system dynamics is followed up to a time t_max=40Ω_cp^-1 and no perturbations are introduced.

We first consider the SMDF f_SM [Eq. (<ref>)] as initial condition for the simulation S1. In this case we set: 

    u_0(x)=U_0[tanh(x-L/4/Δ x)-tanh(x-3L/4/Δ x) -1 ]

with U_0=2V__A, Δ x=2.5 d_p and L=100 d_p.

Then, we performed a second simulation S2, using as initial condition the stationary DF f_eq,|| [Eq. (<ref>)], with

    U(x-v_z)=U_0[ tanh(x-L/4-v_z/Δ x)-tanh(x-3L/4-v_z/Δ x) - 1]

and compared the results of the two simulations. We point out that the expressions in Eqs. (<ref>) and (<ref>) 
describe a smooth jump in velocity at the position x=L/4 along the x-direction; this jump has been mirrored at x=3L/4, in 
order to satisfy the periodic boundary conditions.

As expected, the simulation S1 clearly indicates that the initial distribution SM DF is not an equilibrium and, as a consequence, its 
velocity moments display an oscillatory behaviour with a period equal to the ion gyroperiod. Left column of Fig. <ref> displays the 
contour plots of δ u_y,%(x,t) = (u_y(x,t)-u_y(x,0))/U_0×100 [panel (a)] and δ 
T_%(x,t)=(T(x,t)-T(x,0))/T^∞×100 [panel (c)] in the (x,t) plane, for the simulation S1, where the x range has been 
set to focus on the left half of the spatial box and where T^∞ = v_th,p^2= β_p/2 = 2 in dimensionless units. Significant 
oscillations (about 30%) are recovered in both u_y and T, localized around the shear.

On the other hand, in the simulations S2 with the DF f_eq,||, the system remains at equilibrium (no significant oscillations are 
visible). To better point out the differences between S1 and S2, we considered the temporal profiles of δ u_y,%(x,t) and δ 
T_%(x,t), evaluated at a fixed spatial position x=x^*. The point x^* corresponds to the spatial point where each quantity 
exhibits the largest amplitude oscillations in simulation S1 (vertical black-dashed lines in panels (a) and (c) of Fig. <ref>), with 
respect to the initial condition, i.e. δ u_y,%(x^*,t)=max_x{δ u_y,%(x,t)} and δ T_%(x^*,t)=max_x{δ 
T_%(x,t)}. These temporal profiles are reported in panels (b) and (d) of Figure <ref>, as black (f_SM) and red 
(f_eq,||) curves, respectively. Here, one can realize that no significant oscillations in the signals are recovered in the case of the 
DF f_eq,||, as compared to the case of the SMDF.

To further characterize the differences between the two cases S1 and S2, we also computed the L_2 norm of δ u_y,%(x,t) and 
δ T_%(x,t), defined as L_2(g(x,t)) = √(∫ (g-g_0)^2 dx / L), being g a generic function and g_0=g(x,0). 
Clearly L_2(g(x,t)) is a function of time t. In Figure <ref>, we show, in semi-logarithmic plot, the temporal evolution of 
L_2(δ u_y,%) (top row) and L_2(δ T_%) (bottom row). As it can be appreciated from the plots in this figure, significant 
oscillations with respect to the initial configuration are present in the case of the SMDF (black-solid curves), confirming that this 
distribution is not a HVM equilibrium. On the other hand, no oscillations are visible for the case of the DF f_eq,|| (red-solid curves), 
in which the small departure from the initial configuration (about 10^-4%), is presumably due to the numerical error in the calculation 
of the velocity moments of f_eq,||. 



 §.§ Perpendicular case


For this case, we discretized the four dimensional phase space through N_x=512 grid points in the one-dimensional spatial domain and 
N_v=141 grid points in each velocity direction, while v_max=7 v_th,p. Two simulations have been performed to compare the SM DF 
f_SM (S3) and the equilibrium DF f_eq,⊥ (S4), while the background magnetic field is B_0=B_0  e_z (B_0=1 in 
scaled units) and the ion plasma beta β_p=2v_th,p^2/V__A^2=4. In both simulations the system dynamics is followed up to a time 
t_max=40Ω_cp^-1 and no perturbations are introduced. The initial electric field considered for these simulations is the one 
given by:

    E(x) = E_0[ 1 - tanh( x-L/4/Δ x) + tanh( x-3L/4/Δ x)]

with E_0=1 (in scaled units), Δ x=2.5 d_p and L=100 d_p. Note that the shear has been mirrored to hold periodic boundary 
conditions.

As expected, the simulation S3 indicates that the initial distribution SM DF is not an equilibrium. However, with respect to the 
parallel case, its velocity moments does not display an oscillatory behaviour but some propagating structure is also recovered. Left column 
of Fig. <ref> displays the contour plots of δ u_y,%(x,t) = (u_y(x,t)-u_y(x,0))/U_0×100 [panel (a)] and δ 
T_%(x,t)=(T(x,t)-T(x,0))/T^∞×100 [panel (c)] in the (x,t) plane, for the simulation S3, where the x range has 
been set to focus on the left half of the spatial box and where T^∞ = v_th,p^2= β_p/2 = 2 in dimensionless units. 
Disturbances from the equilibrium, (about 1-2%) are recovered in both u_y and T, mainly localized around the shear but also showing 
a propagating structure.

On the other hand, in the simulations S4 with the equlibrium DF, the system remains at equilibrium (oscillations significantly smaller 
than in S3 are in fact recovered in S4). To better point out the differences between S3 and S4, we considered the temporal profiles of 
δ u_y,%(x,t) and δ T_%(x,t), evaluated at a fixed spatial position x=x^*. The point x^* corresponds to the spatial 
point where each quantity 
exhibits the largest amplitude oscillations in simulation S3 (vertical black-dashed lines in panels (a) and (c) of Fig. 
<ref>), 
with respect to the initial condition, i.e. δ u_y,%(x^*,t)=max_x{δ u_y,%(x,t)} and δ 
T_%(x^*,t)=max_x{δ T_%(x,t)}. These temporal profiles are reported in panels (b) and (d) of Figure <ref>, as 
black and red curves, respectively. Here, one can realize that much smaller oscillations in the signals are recovered in the case of the 
DF f_eq,⊥ (red), as compared to the case of the SMDF (black).

To further characterize the differences between the two cases S3 and S4, we also computed the L_2 norm of δ u_y,%(x,t) and 
δ T_%(x,t), defined as L_2(g(x,t)) = √(∫ (g-g_0)^2 dx / L), being g a generic function and g_0=g(x,0). 
Clearly L_2(g(x,t)) is a function of time t. In Figure <ref>, we show, in semi-logarithmic plot, the temporal evolution 
of L_2(δ u_y,%) (top row) and L_2(δ T_%) (bottom row). As it can be appreciated from the plots in this figure, 
significant departures from the initial configuration are present in the case of the SMDF (black-solid curves), confirming that 
this distribution is not a HVM equilibrium. On the other hand, in the case of the DF f_eq,⊥ (red-solid 
curves) much small departures from the initial configuration (about 10^-2%) are observed, presumably due to the numerical error in the 
calculation of the velocity moments of f_eq,⊥.



§ SUMMARY AND CONCLUSIONS


In this paper, we have derived exact solutions for the system of Hybrid Vlasov-Maxwell equations which represents a stationary 
shearing flow with a uniform magnetic field directed either parallel or perpendicular to the plasma bulk velocity. Plasmas supporting 
shearing flows are found in many situations and a kinetic description is necessary whenever the shear width is of the order of kinetic 
scales, like, for instance, in the case of the Earth's 
magnetopause <cit.>. 

The interest of building up stationary solutions can be related to the problem of describing the propagation and evolution of waves 
in a plasma with a stable shearing flow. The interaction 
between waves and the background inhomogeneity associated with the shearing flow moves the wave energy towards small scales, where
kinetic effects are more effective. Moreover, the presence of a shearing flow can generate wave coupling, with an energy transfer among 
different wave modes. It is clear that, in order to properly study wave propagation, it is necessary that the background structure 
remains stationary; otherwise, a time evolution intrinsic of the background state would superpose to waves, making difficult to 
single out the wave contribution in the overall time evolution.
Another possible application of exact shearing flow solutions can be found in the study of the Kelvin-Helmoholtz instability, which takes 
place in unstable shearing velocity configurations. In fact, though the turbulent stage following the instability saturation should be 
quite insensitive to the details of the initial state, only a stationary unperturbed configuration allows to  
properly describe the linear stage of the instability. 
Therefore, in both cases an exact stationary distribution function is preferable to the simpler shifted Maxwellian DF.

Stationary solutions in various configurations have been described in previous studies of the fully 
kinetic case, i.e., involving the full set of ion and electron Vlasov-Maxwell equations. However, the fully kinetic treatment is 
quite complex and such solutions have rarely been employed in numerical simulations like, for instance, in investigations of the KH instability. In this respect, the set of Hybrid Vlasov-Maxwell equations represents a good compromise, because it correctly describes 
a plasma at scales of the order of or larger than the ion scales, but avoiding the complexity of a fully kinetic treatment. 
In the framework of Vlasov-Maxwell equations, Cerri et al. <cit.> have developed a method to derive approximately stationary 
solutions. The solutions presented here situate in the same framework, but have the advantage to be exactly stationary. 

The starting point of our derivation are previous studies where stationary DFs are derived in the fully kinetic framework (Refs. 
<cit.>) which we revisited and adapted to the hybrid Vlasov-Maxwell description. In particular, we have examined 
the special cases in which the magnetic field is uniform and either parallel or perpendicular to the bulk velocity. In the former 
case the stationary solution have a simple analytical form which can be directly used in numerical simulations; in the latter case, 
the explicit construction of the distribution function is obtained through a simple numerical procedure which integrates particle orbits 
throughout the relevant phase space. In the case of parallel B an isothermal electron fluid have been assumed. In contrast, 
in the case of perpendicular B a nonuniform electron pressure p_e is necessary in order to satisfy the generalized Ohm's law. 
As a consequence, an adiabatic equation for the electron fluid has been added to the set of equations. This aspect represents a novelty 
for the hybrid Vlasov-Maxwell approach, in which an isothermal electron fluid has been routinely assumed. 

The main properties of these solutions have been examined, calculating the associated profiles of bulk velocity, temperatures and heat flux. 
In the shear region, the ion distribution functions are distorted with respect to shifted Maxwellians, with stronger distortions 
for more localized shears. In particular, marked anisotropy and agyrotropy in the ion temperature are generated, and none of the DF 
principal axes is aligned to B. Moreover, a non-vanishing heat flux is present, directed in the plane perpendicular 
to the inhomogeneity direction x. We found that the width of the velocity shear cannot be smaller than ion Larmor radius; this can be 
justified by considering the ion gyromotion which mixes up single particle velocities on the scale of the Larmor radius.
 
The HVM code <cit.> has been employed to verify to what extent the derived configurations remain stationary when used 
as initial conditions in numerical simulations. The time behaviour has been compared with that obtained in the case of a shifted Maxwellian 
distribution function. We found that in the case of our solutions the deviation from the initial condition remains much smaller 
(two orders of magnitude for perpedicular B and more than three orders of magnitude for parallel B) than in the case 
of the shifted Maxwellian. The small deviation from exact stationarity of the former case are probably due to numerical errors in 
the HVM code and, for perpendicular B, also in the procedure integrating particle orbits.

We are planning to use these results for studying the problem of Alfvén wave evolution in a shearing flow plasma. Moreover, we are 
currently working to extend them to the case of an obliquely-directed magnetic field, a configuration commonly observed in the 
Earth's magnetosphere.



§ ACKNOWLEDGEMENT

This work has been supported by the Agenzia Spaziale Italiana under the contract ASI-INAF 2015-039-R.O “Missione M4 di ESA: Partecipazione 
Italiana alla fase di assessment della missione THOR”.




§ ANISOTROPY EVALUATION AT THE CENTER OF THE SHEAR IN THE PARALLEL CASE



In this Appendix, we analytically calculate, in the case of the shear parallel to the background magnetic field, the variance matrix for 
the equilibrium DF f_eq,|| of Eq. <ref>. We focus on the case where U represents a shear layer across which the bulk velocity 
varies from a value -U_0 to U_0. In the spatial positions far from the velocity shear f_eq,|| reduces to a shifted Maxwellian, then 
the variance matrix becomes diagonal σ_ij^∞ = v_th,p^2 δ_ij and 
T_0^∞=T_||^∞=T_⊥^∞, where the upper index “∞" identifies values calculated far from the shear 
layer. We focus the center x=0 of a symmetric shear layer, i.e. U(k_2) is an odd function of k_2, where we expect to find the 
strongest departures from a Maxwellian. In this case, u_y(x=0)=0, while the variance matrix is:

    σ_ij(x=0)= 1/(2π )^3/2 v_th,p^3∫ v_i v_j exp{ -1/2 v_th,p^2[ v_x^2 + ( v_y + U( 
    v_z/Ω_cp) )^2 +v_z^2 ] } d^3 v

In the large scale limit Δ x ≫ρ__Lp (Δ v ≫ v_th,p) and by considering that the typical value for the velocity 
v_z is v_th,p, we can retain the first-order Taylor expansion term of U(v_z/Ω_cp): U(v_z/Ω_cp) ≃ω_0/Ω_cp v_z, being ω_0 ≡. dU/dk_2|_k_2=0

Within this approximation, all the integrals involved in Eq. (<ref>) can be easily calculated and the resulting form is:

    σ(x=0) = 
    [ [                   v_th,p^2                          0                          0;                          0 ( 1+ω_0^2/Ω_cp^2) v_th,p^2         -ω_0/Ω_cp v_th,p^2;                          0         -ω_0/Ω_cp v_th,p^2                   v_th,p^2 ]]


Diagonalizing σ(x=0) implies to rotate the DF into the minimum variance frame (MVF), The eigenvalues of σ(x=0), corresponding 
to the temperatures in the MVF, are λ^(3) < λ^(2) < λ^(1), whose explicit expressions are:

    λ^(3)   =    v_th,p^2 ( 1- ω_0/Ω_cp√(1 + ω_0^2/4Ω_cp^2) + 
    ω_0^2/2Ω_cp^2) 
    λ^(2)   =    v_th,p^2 
    λ^(1)   =    v_th,p^2 ( 1+ ω_0/Ω_cp√(1 + ω_0^2/4Ω_cp^2) + 
    ω_0^2/2Ω_cp^2)

The corresponding eigenvectors are given by:

    ξ^(3)   =   ( √(1 + ω_0^2/4Ω_cp^2)-ω_0/2Ω_cp)  e_y +  e_z 
    ξ^(2)   =    e_x 
    ξ^(1)   =    -( √(1 + ω_0^2/4Ω_cp^2)+ω_0/2Ω_cp)  e_y +  e_z

From the above expressions we can deduce the following informations: 

(i) At x=0, σ_xx = σ_zz, therefore the two temperatures in the directions orthogonal to B are equal, 
i.e. the DF is gyrotropic in the LBF at x=0.

(ii)
By comparing σ_ij^∞ and σ_ij(x=0) in Eq. (<ref>), we note that, at x=0, σ_xx and σ_zz 
keep the same value of σ_ij^∞. Therefore, at the center of the shear, the perpendicular proton temperature remains constant 
(T_⊥(x)=T_⊥^∞), while the parallel temperature increases (T_||(x=0) ≥ T_||^∞). As a consequence, 
at x=0 the perpendicular to parallel proton temperature ratio is 

    T_⊥(x=0)/T_||(x=0) = 1/1+ω_0^2/(2Ω_cp^2) < 1

i.e., the proton parallel temperature is larger than the perpendicular one.

(iii)
The eigenvectors ξ^(m) (m=1,2,3) give the directions of the principal axes of the DF in the velocity space. Far from the shear 
layer the DF is isotropic and the directions of principal axes are arbitrary. At x=0, the principal axis ξ^(2) is in the v_x 
direction, which corresponds to the direction of spatial variation of the bulk velocity u, while the other two principal axes 
ξ^(1) and ξ^(3) are in the v_yv_z plane. The angle γ between ξ^(1) and the v_y axis, which gives the 
directions 
corresponding to the maximum width of the DF, is:

    tanγ = -1/√(1+ω_0^2/4Ω_cp^2)+ω_0/2Ω_cp

In the limit of slowly varying bulk velocity (ω_0 ≪Ω_cp), we get γ≃ -45^∘. Finally, it is worth noting 
that, since the three eigenvalues are all different at x=0, the DF is not gyrotropic with respect to any of the three principal axes in 
the MVF. 

It is interesting to extend the results illustrated above through the numerical evaluation of the temperature in the MVF and in the LBF, 
for several values of the shear width Δ x, where the shear function U(x) is given by Eq. (<ref>). Table <ref> 
reports the values of λ^(i) (with i=1,2,3), T_||, T_⊥ and γ, at the center of the shear x=0. As 
expected, in the large scale limit Δ x ≫ρ__Lp, the results obtained numerically are in good accordance with the 
analytical predictions. By decreasing the width of the shear function Δ x, the analytical calculations tend to diverge from  
numerical evaluations. Note also that, as the shear width Δ x becomes smaller, stronger anisotropies in the LBF (i.e. bigger 
T_||) and in the MVF (i.e. larger ratios between the eigenvalues λ^(i)) are recovered at the center of shear.



§ STATIONARY SOLUTION IN THE LOCAL APPROXIMATION IN THE PERPENDICULAR CASE



In the present Appendix, we report the evaluation of the equilibrium DF, in the perpendicular case, within the so called ”local 
approximation“, i.e. in the simplified case of linearly growing electric field. The single particle motion is calculated by solving Eqs. 
(<ref>), in the particular case in which the electric field is linear E(x)=E_0+α_0(x-x_0), with E_0 and α_0 
constant. This profile for the electric field is not fully realistic, since |E(x)| grows without limit for increasing |x-x_0|. 
However, it can be considered as a local approximation of an electric profile E(x) around a given position x_0, being 
α_0=(dE/dx)(x_0). 

In this case, Eq. (<ref>) reads:

    d^2 x/dt^2 + ω^2 x = e/m_p( E_0 - α_0 x_0) + Ω_cp W_0

where ω^2 = Ω_cp^2 - e α_0/m_p. If Ω_cp^2 > eα_0/m_p, Eq. (<ref>) describes an harmonic 
oscillator of solution:

    x(t) = R_0 sin( ω t + φ) + 1/ω^2[ e/m_p( E_0-α_0 x_0 ) +Ω_cp 
    W_0 ]

with R_0 and φ the amplitude and the phase of the motion, respectively. The constant term in Eq. (<ref>) 
represents the time-averaged x position, i.e. the guiding center position x_c:

    x_c  = 1/ω^2[ e/m_p( E_0-α_0 x_0 ) +Ω_cp W_0 ]

and then

    x(t) = R_0 sin( ω t + φ) + x_c 
        v_x(t)= R_0 ωcos( ω t + φ) 
        v_y(t)= -Ω_cp x(t)+W_0 = -R_0 Ω_cpsin( ω t + φ) + v_cy

where v_cy  = ⟨ v_y(t) ⟩_t = -Ω_cp x_c + W_0. Using Eq. (<ref>) into the v_cy expression, one obtains:

    v_cy = -e/m_p Ω_cp^2( Ω_cp E_0 + α_0 v_0y)

where v_0y=-Ω_cp x_0+W_0 is the particle streamwise velocity component at the position x=x_0. Note also that 
v_0y=v_cy+Ω_cp(x_c-x_0). By using these last expressions, one obtains:

    v_cy=-e/m_p Ω_cp[ E_0+α_0 ( x_c - x_0 )] = -c E(x_c)/B

Hence, the guiding center moves along y with the local E× B drift velocity calculated at the guiding center position 
x_c. Then, in the local approximation case the particle orbit in the v_x v_y plane is an ellipse elongated along v_x (v_y) for 
α_0 <0 (0 < α_0 < m_p Ω_cp^2/e), while it reduces to a circle in the case of uniform electric field α_0=0. In 
the case α_0 ≥ m_p Ω_cp^2/e, x(t) increases linearly or exponentially in time, causing the breakdown of the local 
approximation, as already noticed in Ref. <cit.>. 











Note that, if α_0=0 (i.e. ω^2=Ω_cp^2), then v_x^2+(v_y-v_cy)^2=R_0^2 ω^2. Indicating by 
t_n=(nπ-ϕ)/ω (n integer) an instant of time when x(t_n)=x_c, then v_x(x=x_c)=v_x(t=t_n)=± R_0ω [Eqs. 
(<ref>)]. Hence, v_x^2+(v_y-v_cy)^2=[ v_x(x=x_c)]^2. The same argument show that v_x(x=x_c)=± R_0ω also 
when the electric field is non-uniform (α_0  0). Then, the energy ℰ_0 [Eq. (<ref>)] is:

    ℰ_0=m_p/2( R_0^2 ω^2 + v_z^2 )

In order to explicitly write the form for f_eq,⊥(x, v) [Eq (<ref>)], we manipulate Eq. (<ref>) by using Eqs. 
(<ref>) and Eq. (<ref>) and by expressing x_c in terms of the particle position and velocity through Eqs. 
(<ref>):

    ℰ_0 = m_p/2[ v_x^2+ω^2/( Ω_cp-cα_0/B )^2( v_y + c/B E(x) 
    )^2 + v_z^2 ]

This expression, which is the argument in the exponential of the DF f_eq,⊥, suggests that f_eq,⊥ is a shifted bi-Maxwellian 
with different temperatures T_u=T_y and T_⊥ u=T_xz parallel to and in the plane perpendicular to the bulk flow, respectively. 
The temperature ratio is T_u/T_⊥ u=(Ω_cp-cα_0/B)^2/ω^2; using the explicit expression for ω we find:

    T_u/T_⊥ u=1-1/Ω_cpc α_0/B=1-e α_0/m_p Ω_cp^2

Therefore, T_u > T_⊥ u (T_u < T_⊥ u) when E(x) decreases (increases) with increasing x. Note that, in order to 
have a positive temperature T_u, the spatial derivative of the electric field has an upper limit: α_0 < m_p Ω_cp^2/e. 
This condition is the same which avoids the breakdown of the local approximation in the single ion dynamics, as found in the previous 
section. Of course, for a uniform electric field (α_0=0), we have T_u = T_⊥ u.

By requiring that the uniform density n^LA (x)=n_0, it can be easily shown that:

    f_eq,⊥^LA(x, v) = n_0/(2π)^3/2 v_th,p^3( T_⊥ u/T_u)^1/2
    e^-1/2v_th,p^2{ v_x^2+
    T_⊥ u/T_u[ v_y + c/B E(x) ]^2 + v_z^2 }

which depends on two arbitrary constants n_0 and v_th,p, while the ratio T_u/T_⊥ u is given by the equation 
(<ref>).

The bulk velocity associated with the DF is u^LA=∫ v  f_eq,⊥^LA  d^3 v/n^LA. The bulk velocity 
components u^LA_x and u^LA_z are both vanishing, while the component u^LA_y= - c E(x)/B, indicating that the bulk velocity 
coincides with the local E× B drift velocity.

The considered DF is an exact stationary solution of the HVM equations Eqs. (<ref>–<ref>), with Eq. (<ref>) for the 
pressure closure. The DF f_eq,⊥^LA is a stationary solution of the Vlasov equation [Eq. (<ref>)], because it is a 
function only of the motion constant. The electric field profile is linear and its profile is correctly given by the u× 
B term in Eq. (<ref>) (then the stationarity holds for an electron pressure p_e constant and uniform). Since the electric field is 
irrotational, the magnetic field is stationary. Therefore, the considered configuration is an exact stationary solution of the whole set of 
the HVM equations.






99
fujimoto98
M. Fujimoto, T. Terasawa, T. Mukai, Y. Saito, T. Yamamoto, S. Kokubun, J. Geophys. Res., 103, A3, 4391 (1998).

hasegawa04
H. Hasegawa, M. Fujimoto, T. D. Phan, H. Reme, A. Balogh, M. W. Dunlop, C. Hashimoto, R. TanDokoro, R., Nature, 430, 755 (2004).

sckopke81
N. Sckopke, G. Paschmann, G. Haerendel, B. U. Sonnerup, S. J. Bame, T. G. Forbes, E. W. Hones, C. T. Russell, J. Geophys. Res., 86, 
A4 2099 (1981).

brunocarbone13
R. Bruno, V. Carbone, Liv. Rev. in Solar Phys. 2 (2013).

hamlin13
N. D. Hamlin, W. I. Newman, Phys. Rev. E, 87, 043101 (2013).

chandrasekhar61
S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, Dover (1961).

heyvaerts83
J. Heyvaerts, E. R. Priest, Astron. Astrophys., 117, 220 (1983).

pucci14
F. Pucci, M. Onofri, F. Malara, Astrophys. J., 796, 43 (2014).

vasconez15
C. L. Vásconez, F. Pucci, F. Valentini, S. Servidio, W. H. Matthaeus, F. Malara, Astrophys. J., 815, 7 (2015).

pucci16
F. Pucci, C. L. Vásconez, O. Pezzi, S. Servidio, F. Valentini, W. H. Matthaeus, F. Malara, J. Geophys. Res., 121 1024 (2016).

valentini17
F. Valentini, C. L. Vásconez, O. Pezzi, S. Servidio, F. Malara, F. Pucci, Astron. Astrophys., 599, A8, (2017).

hofmann75
I. Hofmann, Plasma Phys. 17, 143 (1975).

einaudi86
G. Einaudi, F. Rubini, Phys. Fluids 29, 2563 (1986).

chen90
X. L. Chen, P. J. Morrison, Phys. Fluids B 2, 495 (1990).

faganello10
M. Faganello, F. Pegoraro, F. Califano, L. Marradi, Phys. Plasmas 17, 062102 (2010).

chacon03
L. Chacón, D. A. Knoll, J. M. Finn, Phys. Lett. A, 308, 187 (2003).

chen93
S. H. Chen, M. G. Kivelson, J. T. Gosling, R. J. Walker, A. J. Lazarus, J. Geophys. Res., 98, 5727 (1993).

kivelson95
M. G. Kivelson, S. H. Chen, Geophys. Monogr. Ser., 90, edited by P. Song, P., B. U. Ö. Sonnerup and M. F. Thomsen, p. 257, AGU, 
Washington, D. C. (1995).

fairfield00
D. H. Fairfield, A. Otto, T. Mukai, S. Kokubun, R. P. Lepping, J. T. Steinberg, A. J. Lazarus, T. Yamamoto, J. Geophys. Res., 105, 21 
(2000).

faganello17
M. Faganello, F. Califano, J. Pkasma. Phys., 83, 535830601 (2017).

malara89
F. Malara, G. Einaudi, A. Mangeney, J. Geophys. Res., 94, 11813 (1989).

matsumoto04
Y. Matsumoto, M. Hoshino, Geophys. Res. Lett., 31, L02807 (2004)

nakamura04
T. K. Nakamura, D. Hayashi, M. Fujimoto, I. Shinohara, Phys. Rev. Lett., 92, 145001 (2004).

matsumoto06
Y. Matsumoto, M. Hoshino, J. Geophys. Res., 111, A05213 (2006).

matsumoto07
Y. Matsumoto, K. Seki, J. Geophys. Res., 112, A06223 (2007).

faganello08
M. Faganello, F. Califano, F. Pegoraro, Phys. Rev. Lett., 100, 015001 (2008).

henry13
P. Henri, S. S. Cerri, F. Califano, F. Pegoraro, C. Rossi, M. Faganello, O. Šebek, P. M. Trávníček, P. Hellinger, J.
T. Frederiksen, A. Nordlund, S. Markidis, R. Keppens, G. Lapenta, Phys. Plasmas, 20, 102118 (2013).

rossi15
C. Rossi, F. Califano, A. Retinó, L. Sorriso-Valvo, P. Henri, S. Servidio, F. Valentini, A. Chasapis and L. Rezeau, Phys. Plasmas 22, 122303 (2015).

dimare18
F. Di Mare, L. Sorriso-Valvo, A. Retinò, H. Hasegawa, F. Malara, Transition to turbulence in the Kelvin-Helmholtz instability at the terrestrial Magnetopause, submitted to J. Geophys. Res. (2018).

cai90
D. Cai, L. R. O. Storey, T. Neubert, Phys. Fluids B: Plasma Physics, 2, 75 (1990).

ganguli88
G. Ganguli, Y. C. Lee, P. J. Palmadesso, Phys. Fluids 31, 823 (1988).

nishikawa88
K.-I. Nishikawa, G. Ganguli, Y. C. Lee, P. J. Palmadesso, Phys. Fluids 31, 1568 (1988).

mahajan00
S. M. Mahajan, R. D. Hazeltine, Phys. Plasmas 7, 1287 (2000).

roytershteyn08
V. Roytershteyn, W. Daughton, Phys. Plasmas 15, 082901 (2008).

karimabadi13
H. Karimabadi, V. Roytershteyn, M. Wan, W. H. Matthaeus, W. Daughton, P. Wu, M. Shay, B. Loring, J. Borovsky, 
E. Leonardis, S. C. Chapman, T. K. M. Nakamura, Phys. Plasmas, 20, 012303 (2013).

cerri13
S. S. Cerri, P. Henri, F. Califano, D. Del Sarto, M. Faganello, F. Pegoraro, Phys. Plasmas 20, 112112 (2013).

valentini07 
F. Valentini, P. Travnicek, F. Califano, P. Hellinger, A. Mangeney, J. Comput. Phys. 225, 753 (2007).

servidio12 S. Servidio, F. Valentini, F. Califano, P. Veltri, Phys. Rev. Lett. 108, 045001 (2012)

matthaeus14 W. H. Matthaeus, S. Oughton, K. T. Osman, S. Servidio, M. Wan, S. P. Gary, M. A. Shay, F. Valentini, V. Roytershteyn, 
H. Karimabadi, S. C. Chapman, Astrophys. J. 790, 155 (2014)

franci15 L. Franci, A. Verdini, L. Matteini, S. Landi and P. Hellinger, Astrophys. J. Lett. 804, L39 (2015). 

servidio15 S. Servidio, F. Valentini, D. Perrone, A. Greco, F. Califano, W. H. Matthaeus and P. Veltri, J. Plasma Phys. 81, 
325810107 (2015)

valentini16 F. Valentini, D. Perrone, S. Stabile, O. Pezzi, S. Servidio, R. De Marco, F. Marcucci, R. Bruno, B. Lavraud,
J. De Keyser, G. Consolini, D. Brienza, L. Sorriso-Valvo, A. Retinó, A. Vaivads, M. Salatti and P. Veltri, New J. Phys. 18, 125001 
(2016)

cerri16 S.S. Cerri, F. Califano, F. Jenko, D. Told and F. Rincon, Astrophys. J. Lett. 822, L12 (2016).

pezzi17a O. Pezzi, T.N. Parashar, S. Servidio, F. Valentini, C.V. Vásconez, Y. Yang, F. Malara, W.H. Matthaeus & P. Veltri, 
The Astrophys. Journal 834, 166 (2017).

cerri14 S.S. Cerri, F. Pegoraro, F. Califano, D. Del Sarto and F. Jenko, Phys. Plasmas 21, 112109 (2014).

delsarto16 D. Del Sarto, F. Pegoraro and F. Califano, Phys. Rev. E 93, 053203 (2016).

delsarto18
D. Del Sarto and F. Pegoraro, MNRAS, 475, 181 (2018).

valentini11 F. Valentini, D. Perrone, and P. Veltri, Astrophys. J. 739, 54 (2011)

valentini14 F. Valentini, S. Servidio, D. Perrone, F. Califano, W. H. Matthaeus, and P. Veltri, Phys. Plasmas 21, 082307 
(2014)

harris62
E. G. Harris, Nuovo Cimento, 23, 117 (1962).

krall73
N. A. Krall, A. W. Trivelpiece, Principle of Plasma Physics, McGraw Hill, New York, (1973). 

cerri17 S.S. Cerri, S. Servidio and F. Califano, Astrophys. J. Lett. 846, L18 (2017).

pezzi17b 
O. Pezzi, T.N. Parashar, S. Servidio, F. Valentini, C.V. Vásconez, Y. Yang, F. Malara, W.H. Matthaeus & P. Veltri, J. Plasma Phys. 83, 905830105 (2017).

pezzi17c 
O. Pezzi, F. Malara, S. Servidio, F. Valentini, T.N. Parashar, W.H. Matthaeus and P. Veltri, Phys. Rev. E 96, 023201 (2017).



naitou80
S. Naitou, S. Tokuda, T. Kamiura, J. Comput. Phys, 38, 265 (1980).

































