Electron-scale shear instabilities: magnetic field generation and particle acceleration in astrophysical jets

Strong shear flow regions found in astrophysical jets are shown to be important dissipation regions, where the shear flow kinetic energy is converted into electric and magnetic field energy via shear instabilities. The emergence of these self-consistent fields make shear flows significant sites for radiation emission and particle acceleration. We focus on electron-scale instabilities, namely the collisionless, unmagnetized Kelvin-Helmholtz instability (KHI) and a large-scale dc magnetic field generation mechanism on the electron scales. We show that these processes are important candidates to generate magnetic fields in the presence of strong velocity shears, which may naturally originate in energetic matter outburst of active galactic nuclei and gamma-ray bursters. We show that the KHI is robust to density jumps between shearing flows, thus operating in various scenarios with different density contrasts. Multidimensional particle-in-cell (PIC) simulations of the KHI, performed with OSIRIS, reveal the emergence of a strong and large-scale dc magnetic field component, which is not captured by the standard linear fluid theory. This dc component arises from kinetic effects associated with the thermal expansion of electrons of one flow into the other across the shear layer, whilst ions remain unperturbed due to their inertia. The electron expansion forms dc current sheets, which induce a dc magnetic field. Our results indicate that most of the electromagnetic energy developed in the KHI is stored in the dc component, reaching values of equipartition on the order of $10^{-3}$ in the electron time-scale, and persists longer than the proton time-scale. Particle scattering/acceleration in the self generated fields of these shear flow instabilities is also analyzed.


Abstract.
Strong shear flow regions found in astrophysical jets are shown to be important dissipation regions, where the shear flow kinetic energy is converted into electric and magnetic field energy via shear instabilities. The emergence of these self-consistent fields make shear flows significant sites for radiation emission and particle acceleration. We focus on electron-scale instabilities, namely the collisionless, unmagnetized Kelvin-Helmholtz instability (KHI) and a large-scale dc magnetic field generation mechanism on the electron scales. We show that these processes are important candidates to generate magnetic fields in the presence of strong velocity shears, which may naturally originate in energetic matter outburst of active galactic nuclei and gamma-ray bursters. We show that the KHI is robust to density jumps between shearing flows, thus operating in various scenarios with different density contrasts. Multidimensional particle-in-cell (PIC) simulations of the KHI, performed with OSIRIS, reveal the emergence of a strong and large-scale dc magnetic field component, which isnot captured by the standard linear fluid theory. This dc component arises from kinetic effects associated with the thermal expansion of electrons of one flow into the other across the shear layer, whilst ions remain unperturbed due to their inertia. The electron expansion forms dc current sheets, which induce a dc magnetic field. Our results indicate that most of the electromagnetic energy developed in the KHI is stored in the dc component, reaching values of equipartition on the order of 10 −3 in the electron time-scale, and persists longer than the proton time-scale. Particle scattering/acceleration in the self generated fields of these shear flow instabilities is also analyzed.

Introduction
Relativistic jets are found in a wide range of extreme astrophysical scenarios like active galactic nuclei (AGN) and gamma-ray bursts (GRBs) [Bridle (1984), Mirabel (1999)]. The energetic outflows of plasma associated with astrophysical jets represent massive sources of free-energy for collisionless plasma instabilities to operate. The onset of plasma instabilities play a central role in dissipating the jet's kinetic energy into electric and magnetic turbulence [Gruzinov & Waxman (1999), Medvedev (1999)] resulting in particle acceleration to ultra-high energies and nonthermal radiation emission. A deep understanding of these processes and their interplay is challenging, requiring full kinetic simulations to address their highly nonlinear nature. First principle modeling of these processes are, however, computationally intensive due to the wide range of temporal and spatial scales involved. Therefore, full kinetic simulations demand massive computational resources and advanced numerical and visualization techniques.
Much attention has been devoted to relativistic shocks, which are thought to be a strong mechanism for particle acceleration. Such shocks arise from the collision and bulk interpenetration of different velocity plasma shells, due to either intermittencies or inhomogeneities of the ejecta. The Weibel [Weibel (1959)] and the purely transverse two stream instabilities ] act as the dissipation mechanism in these scenarios, and are critical for shock formation. A vast number of fully kinetic simulations have focused on shock formation settings, where long-lived equipartition magnetic field generation via the Weibel instability has been observed , , Frederiksen (2004), Nishikawa (2005)]. A Fermi-like particle acceleration process has also been identified in simulations of long-term evolution of collisionless shocks [Spitkovsky (2008), Martins et al.(2009)]. These previous works have only considered shearless flows.
However, in addition to bulk plasma collision sites, the transition layers of shear flows have also been probed [Gruzinov (2008)] and shown to constitute important dissipation regions [Alves and Grismayer (2012), Grismayer and Alves (2013), Liang (2013)]. Increasing evidence has pointed to a general stratified organization of the structure of jets in AGN and GRBs [Granot (2003), Rieger (2004)], where different internal shear layers can occur; rotating inner cores vs. axially moving outer shells, or fast inner cores vs. slower outer shells. Moreover, external shear layers, resulting from the interaction of the jet with the interstellar medium, may also be considered. In these scenarios, collisionless shear instabilities such as the Kelvin-Helmholtz instability (KHI) [D'Angelo (1965), Gruzinov (2008), Zhang et al. (2009)] play a role in the dissipation of the jet kinetic energy into electric and magnetic turbulence [Zhang et al. (2009), Alves and Grismayer (2012), Liang (2013)]. In fact, the combined effect of shear flow with collisionless shock formation has not yet been addressed, and may also lead to interesting novel phenomenology since density inhomogeneities generated by shear instabilities can also constitute important scattering sites for particle acceleration. Recent fully kinetic simulations of shear flow settings have probed the self-consistent evolution of the electron-scale KHI, demonstrating that the operation of kinetic effects are responsible for the generation of large-scale, equipartition magnetic fields [Alves and Grismayer (2012), Grismayer and Alves (2013)]. Nonthermal particle acceleration has also been investigated in hybrid electron-positron-ion shear flows , Liang (2013)], with different pair/ion ratio compositions, showing spectral features similar to those found in GRBs.
In laboratory experiments, scenarios where the unmagnetized KHI can be triggered are now being examined both in the collisional [Harding (2009), Hurricane (2012] and in the collisionless regimes [Kuramitsu (2012)] (the latter is explored in this paper). In this work we focus on electron-scale processes triggered by velocity shears, namely the unmagnetized KHI and a dc magnetic field generation mechanism. In Section 2, we develop the linear theory for the cold unmagnetized KHI, and analyze the impact of density contrast between sharp shearing flows. We find the onset of the KHI is robust to density contrasts, allowing for a strong development in various density contrast regimes (inner shears with low density contrasts, and outer shears with high density contrasts). We then extend the analysis to finite shear gradients, where we find that KHI growth rate decreases with increasing shear gradient length. Particle-in-cell (PIC) simulations are performed to verify the theoretical predictions. At late times, PIC simulations reveal the formation of a large-scale, dc magnetic field extending along the entire shear surface between flows, which is not predicted by the linear KHI theory. This dc magnetic field is the dominant feature of the magnetic field structure of the instability at late times. In Section 3, we find that the dc magnetic field results from kinetic effects associated with electron mixing between shearing flows, which is driven by the nonlinear development of the cold KHI. The dc magnetic field generation is discussed and an analytical model is developed that captures the main features of the dc magnetic field evolution and saturation. In Section 4, we analyze the dynamics of the electrons in the self-generated fields. The electrons are scattered in the self consistent electric and magnetic fields generated by the KHI, and are accelerated to high energies. We discuss the particle energy spectra resulting from the development of these shear instabilities, and we investigate the mechanism underlying the acceleration of energetic particles using advanced particle tracking diagnostics.

The cold, unmagnetized, electron-scale KHI
The KHI is a well known instability that is driven by velocity shear. This instability was first derived for neutral shearing fluids within the hydrodynamic framework, where the flows interact via pressure gradients [Chandrasekhar (1961), Drazin & Reid (1981)]. The KHI in charged fluids (plasmas) has also been studied within the MHD framework [D'Angelo (1965), Thomas & Winske (1991)], where the shearing flows also interact via electric and magnetic fields, in addition to pressure gradients. In both these frameworks, the length (and time) scales involved are much larger than the kinetic scales associated with the particles that make up the neutral or charged fluid. In this work, we study the KHI undergone by the electron fluid component of the plasma; the ion dynamics, due to their large inertia, are neglected and are assumed to be unperturbed during the development of the electron-scale KHI. The physics underlying the development of the KHI at the electron-scale is different from the more usual hydrodynamics and MHD forms of the instability, and leads to interesting features that are not observed in more macroscopic frameworks. It is important to note that the KHI can occur at various scales (from electron-kinetic to MHD scales), and that the cross-scale connection and interplay of these instabilities remains to be understood. In this Section, we present the linear two-fluid theory of the electron-scale KHI for an initially cold and unmagnetized plasma shear flow. We generalize for arbitrary velocity and density profiles, and derive analytical solutions for step-like velocity shear and density profiles. The theoretical results are then compared and verified with PIC simulations.

Linear two-fluid theory
In the case of an initially unmagnetized plasma in equilibrium, there is no need for a pressure term to balance out the magnetic force. The equilibrium of the system is then naturally obtained by taking the cold limit of the plasma. In order to describe the linear regime of the electron KHI, we employ the relativistic fluid theory of plasmas. The equations that constitute this theoretical framework are ∂ρ ∂t The equations are written in SI units. Eq. (1) and Eq. (2)  n  (x) Figure 1. Theoretical setting for a 2D shear flow, with arbitrary velocity and density profiles v 0 (x) and n 0 (x), respectively and v are the linear momentum and velocity vectors, where γ = (1 − v 2 /c 2 ) −1/2 is the relativistic Lorentz factor; c is the speed of light, m e and e are, respectively, the electron mass and electron charge, and 0 is the electric permittivity of vacuum. We assume a two-dimensional (2D) cold relativistic shear flow with initial velocity and density profiles described by, respectively ( Figure 1). Due to the 2D assumption, the system lies in the xy plane and sustains electric and magnetic fields of the form: Since we are first interested in the linear evolution of the system, we linearize all physical quantities: The subscripts 0 and 1 denote zeroth and first-order quantities, respectively. External electric and magnetic fields are absent and therefore the zeroth-order quantities of these fields are zero. Since the structures produced by the instability emerge along the y direction, we look for solutions of the form: The ions are assumed to be infinitely massive and thus free streaming, and consider only perturbations in the electron dynamics. The linearized equation of continuity for the electron fluid reads ∂ ∂t Substituting the solution form of Eq. (8) into n 1 and v 1 , we arrive at: The linearized equation of motion of the electrons is given by The zeroth and first-order momentum are, respectively, (12) into Eq. (11) and solving for v 1 , we arrive at Combining Eq. (13) and Eq. (14) with Eq. (10) we compute the perturbed current density J 1 = −e (n 0 v 1 + n 1 v 0 ), We now couple these current densities to Maxwell's equations in order to close our system of equations. The linearized form of Eq. (3) and Eq. (4) are written as: These two equations, Eq. 17 and Eq. 18 are combined by taking the curl of Eq. (17) and substituting in Eq. (18), Splitting Eq. (19) into its components and inserting the candidate plane-wave solutions of the form of Eq. (8), we obtain Next, we insert the current densities, Eq. (15) and Eq. (16) into the equations Eq. (20) and Eq. (21) which, after some algebra, leads to the following equation describing the linear electromagnetic eigenmodes of the system: where the functions A, B, and C are: Here, ω p = n 0 e 2 /γ 0 0 m e denotes the relativistic electron plasma frequency (which is a function of x as it depends on the plasma density profile n 0 (x)). For general density and velocity fields, Eq. (22) may only be solved numerically. However, analytical solutions may be obtained for special settings where Eq. (22) is simplified. We now derive an analytical solution of Eq. (22) for such a setting. 2.1.1.
Step velocity shear and density profiles We consider the following step-function velocity shear profile, and step-function density profile, The values v 0 and n ± are constants. This setting translates into two counter propagating flows with different densities which shear at the plane x = 0 ( Figure 2) and generalizes the standard configuration of equal density flows. Inserting these profiles into Eq. (22), we note that the functions A and C are step-like functions, and that the function B is proportional to δ(x) since it contains the derivative of the density profile (embedded in the plasma frequency, ω p ). We begin by integrating Eq. (22) for x > 0 and x < 0 separately and later join the two solutions at the discontinuity plane x = 0. Applying the well known dielectric boundary conditions to our system, we deduce that E y , being the component of the electric field tangential to the dielectric interface, must be continuous, i.e., E y1 (0 + ) = E y1 (0 − ) ≡ E y1 (0). Thus, for x = 0, the functions A and C are constants and B = 0, leading evanescent wave solutions: where k ⊥ = k 2 + ω 2 p+ /c 2 − ω 2 /c 2 and ω p± is the electron plasma frequency of the n ± plasma. The dispersion relation is finally deduced from the derivative-jump of the electric field at the discontinuity plane. To obtain the derivative jump-condition we perform the standard procedure of integrating Eq. (22) over the interval − < x < , and then take the limit → 0. The first term of Eq. (22) is trivially integrated and the integration of the third term yields 0 in the limit → 0. The second term, however, is the product between a Heaviside step-function and a Dirac delta function δ(x), and is to be evaluated as follows: The derivative jump-condition is thus given by, Finally, manipulating Eq. (28) we obtain the following dispersion relation where β 0 = v 0 /c, ω = γ 0 ω/ω p+ and k = γ 0 kv 0 /ω p+ (where ω p± = n ± e 2 /γ 0 0 m e ) are respectively the normalized frequency and wave number in the dispersion relation. Although this model contains ideal velocity and density profiles, it is a useful tool to provide insights into the behaviour of the KHI in the presence of density contrasts between shearing flows. The density contrast is embedded in the dispersion relation through the density ratio, n + /n − . In the density symmetric limit, n + /n − = 1, Eq. (29) reduces to a biquadratic equation in ω , and we recover the analytical solution presented in [Gruzinov (2008)]: Eq. (30) gives the growth rate of the unstable modes, and is plotted in which corresponds to the KHI dispersion relation obtained in the ideal hydrodynamic model for a symmetric shear flow in the absence of surface tension and gravity [Chandrasekhar (1961)]. The two fluids plasma model dispersion relation differs here from the classical hydrodynamics results by introducing a cut-off at k = 1. There is, therefore, a maximum value of the curve that corresponds to the growth rate (Γ max ) of the fastest growing mode (k max ). These quantities satisfy ∂ k Γ = 0 and are given by: and The real part of ω vanishes over the range of unstable modes meaning that the unstable modes are purely growing waves which is consistent with the symmetry of the system. Note that these electron-scale unstables modes occur when the plasma is considered to be cold, i.e., v th v 0 , whereas compressible MHD or Hydro modes in an initially unmagnetized plasma are only unstable for v 0 < √ 2c s = √ 2v th m e /m i , which correspond to very slow (or very hot) flows [Miura and Pritchett 1982]. Therefore, shear flow instabilities in initially unmangetized conditions with fast drift velocities (relative to the temperature) can only develop on the electron-scale.
In the case of a density jump (n + /n − > 1), Eq. (29) has to be solved numerically. Figure 3 illustrates the effect of the density asymmetry on the growth rate of the unstable modes for multiple values of n + /n − . The values of the density ratio are changed assuming n + fixed so that the normalizing frequency, ω p+ /γ 0 , and wave number, ω p+ /(v 0 γ 0 ), which determine the axes scales of Figure 3, remain constant. We also consider that n + corresponds to the denser flow. Thus, larger density ratios are achieved by lowering the value of n − . The qualitative evolution of the unstable modes is independent of the value of n + /n − , indicating that the general features of the instability are maintained. When n + /n − > 1 the frequency ω acquires a real part over the range of unstable modes leading to propagation ( Figure 4). The drifting character of the unstable modes results from an unbalanced interaction when each flow has different densities. The dispersion relation Eq. (29) in the small k limit reduces to where r = n − /n + . This asymptotic result can be verified in Figure 3 and Figure  4. However this result do not coincide with the dispersion relation obtained in the ideal hydrodynamics model when there is a density jump, Γ hydro = 2k √ r/(1 + r). As we noticed before the two growth rates only coincide in the small k limit when r = 1. In the regime n + /n − > 1, the unstable oscillations develop differently in each flow due to their different densities. The growing oscillations are more strongly manifested in the lower density flow (n − ) and will thus drift in the direction of the n − bulk flow. On the other hand, in the density symmetric regime, the surface interaction between flows is balanced; the unstable modes develop equally in each flow, leading to the development of purely growing waves, as previously discussed. The typical growth rate of the KHI (Γ max /ω p+ ), as was observed in Figure 3, slows down as n + /n − increases. This is because the shear surface current sheets decrease as n − is lowered. In the limit n − → 0 (n + /n − → ∞), we obtain a free streaming plasma in vacuum where the development of the KHI is inhibited, Γ max /ω p+ → 0, as expected. The scaling relations of the KHI with n + /n − are shown in Figure 5. In the similar density regime, n + /n − ≈ 1, the growth rate scales as Γ max /ω p+ ∝ (n + /n − ) −1/4 for both relativistic and non-relativistic shears. In the high density contrast regime, n + /n − 1, the growth rate scales as Γ max /ω p+ ∝ (n + /n − ) −1/3 for non-relativistic shears, and Γ max /ω p+ ∝ (n + /n − ) −1/2 for highly-relativistic shears. Note also that in the n + → 0 limit (ω + p → 0), a scenario where the n − plasma 1 10 100 1000 Figure 5. Scalings of the growth rate of the KHI (Γ max /ω p+ ) with the density ratio between shearing flows. The blue and red curves characterize non-relativistic (γ ≈ 1) and highly-relativistic (γ 1) settings, respectively.
streams in vacuum, the KHI shuts down. At large n + /n − regimes, the KHI dominates over other common plasma instabilities in unmagnetized scenario such as the Weibel and Two-Stream instabilities. The growth rates of the Weibel [Silva et al.(2002)] and Two-Stream instabilities [O'Neil (1971)] scale as Γ Weibel /ω p+ ∝ (n + /n − ) −1/2 and Γ 2−stream /ω p+ ∝ (n + /n − ) −1/3 , respectively. Both growth rates decay more rapidly with n + /n − than the growth rate of the KHI. The physics of large density-contrast settings will thus be mainly determined by the evolution of the KHI. Hence, in realistic astrophysical settings with high density contrasts, where various plasma instabilities are triggered simultaneously, magnetic field generation can also be attributed to the development of the KHI.

Effect of mobile ions
The influence of mobile ions on the theory previously shown is easily incorporated. The ion fluid obeys the same equations as the electron fluid and only the charge and the mass of the ions are the physical parameters that could impact the dispersion relation. For the sake of clarity, we will restrict ourselves to initial equal density plasma. Following the exact same derivation as aforesaid, we obtain a differential equation with the same form of Eq.(22) where the plasma frequency needs to be renormalized, ω 2 p → ω 2 p (1 + m e /m i ). In the case of heavy ions, m i m e the effect can be considered negligible. On the other hand, for an electron-positron plasma, the transverse wavenumber k ⊥ (see Eq.(26)) is rescaled with the plasma frequency (that is multipled by √ 2) and so is the wave number k max associated to the maximum growth rate that peaks at Γ max = 1/2.

Comparisons with PIC simulations
Numerical simulations were performed with OSIRIS [Fonseca et al.(2003), ], a fully relativistic, electromagnetic, and massively parallel PIC code. We have simulated 2D systems of shearing slabs of cold (v 0 v th , where v th = 10 −3 c is the thermal velocity) unmagnetized electron-proton plasmas with a realistic mass ratio m p /m e = 1836 (m p is the proton mass), and evolve it until the electromagnetic energy saturates on the electron time scale. We explored a subrelativistic shear flow scenario with v 0 = 0.2c. The setup of the numerical simulations is prepared as follows. The shear flow initial condition is set by a velocity field with v 0 pointing in the positive x 1 direction, in the upper and lower quarters of the simulation box, and a symmetric velocity field with −v 0 pointing in the negative x 1 direction, in the middle-half of the box. Note that the coordinates (x 1 , x 2 , x 3 ) used in the PIC simulations correspond to the cartesian coordinates (y, x, −z) of the theory presented in the previous sections. Initially, the systems are charge and current neutral, and the shearing flows have equal densities. The simulation box dimensions are 10 × 10 (c/ω p ) 2 , where ω p = (ne 2 / 0 m e ) 1/2 is the plasma frequency, and we use 20 cells per electron skin depth (c/ω p ) in the longitudinal direction and 200 cells per electron skin depth (c/ω p ) in the transverse direction. Periodic boundary conditions are imposed in every direction and we use 36 particles per cell. In order to ensure result convergence, higher numerical resolutions and more particles per cell were tested.
2.2.1. Equal density shear flows We begin by analyzing a subrelativistic shear scenario (v 0 = 0.2 c) where the counter streaming flows have equal densities. The evolution of the electron density of the system is depicted in Figure 6, where the signature roll-up dynamics at the end of the linear phase of the KHI is observed. The protons of the system remain unperturbed (free-streaming) at these time scales due to their inertia. The wavelength of the growing perturbations in the electron density measure 2 c/ω p , which corresponds to the wavelength of the fastest growing mode given by Eq. (33). The magnetic field structure excited by the instability is shown in Figure 7. The first inset of    The physical picture underlying the growth and evolution of the dc mode of the magnetic field will be discussed later in Section 3.

Different density shear flows
The density contrast effects predicted by the theoretical two fluid model have also been verified with numerical simulations. Figure 9 shows the development of the electron density structures for a density contrast setting with n + /n − = 10. The KHI modulations that eventually turn into vortices are strongly manifested in the lower density plasma cloud (represented by the blue flow in Figure 4). The typical length of these modulations is larger than those of the density symmetric case, as predicted by the theoretical model, measuring λ 3.3c/ω p+ . This value agrees with the theoretical wavelength of the fastest growing mode, λ max = 3.1c/ω p+ . The self-generated magnetic field structure is represented in figure 10, where the asymmetry in the evanescent behaviour of the surface mode in the different density regions can be observed. The growth rate of the instability is lowered with respect to the equal-density case and is in good agreement with the linear theory ( figure 8 (b)).

Finite velocity shear gradient
The analytical treatment of the effect of a finite velocity shear gradient (smooth velocity shear profile) on the development of the electron-scale KHI is not trivial. The details of the underlying mathematics and numerics can be found in the Appendix.  We analyse the development of the electron-scale KHI for smooth velocity shear profile given by Electron-scale shear instabilities: magnetic field generation and particle acceleration in astrophysical jets16 analytical solution for the growth rate of the instability can be found for small k ⊥ L (see Appendix), which reads where Γ 0 max /ω p = 1/8 is the growth rate obtained for a step velocity profile. The growth rate of the step velocity profile is recovered for L = 0, and decreases for increasing k ⊥ L. The wavenumber corresponding to the maximum growth rate follows a similar trend by slightly decreasing when the parameter k ⊥ L increases. For arbitrarily large k ⊥ L, an exact numerical solution for the dispersion relation of the electron-scale KHI can be found, and we discuss a numerical scheme in the Appendix.
The above analytical and numerical results have been verified with PIC simulations. The setup of the simulations is identical to those previously described, only replacing the discontinuous velocity profile by the smooth function v 0 (x) = V 0 tanh(x/L). This profile is also used in the numerical algorithm to solve the dispersion relation in the finite gradient shear scenario. The measurement of the maximum growth rate in the simulation is done by following in time the peak of the Fourier spectrum of one of the field structures during the linear phase of the instability. Figure 11

dc magnetic field generation in unmagnetized shear flows
For the sake of completeness, we review in this Section the main results of the dc magnetic field generation mechanism in unmagnetised shear flows which are outlined in [Grismayer and Alves (2013)]. We then present a more detailed analysis of the equipartition fields and the dependence of the equipartition number on the dimensionality of the model. In the previous Section, numerical simulations showed the growth of a dc (k = 0) magnetic field mode (Figure 7 (c) and Figure 10 (c)), which is not predicted by the linear fluid theory (Figure 3), Γ (k = 0)) nor has it been previously identified in MHD simulations and only kinetic simulations [Alves and Grismayer (2012), Grismayer and Alves (2013), ] have been able to capture this mode. The growth of the dc magnetic field mode results from a current imbalance due to electron mixing across the shear interface, while the ion flows remain almost unperturbed due to their inertia. The orientation of the dc magnetic field peak is determined by the proton current structure. The mixing arises due to the deformation of the electron interface between the two flows which, in the linearized fluid calculations, is not accounted for and, in zeroth order, remains fixed. Alternatively, we find that the physics describing the formation of a dc mode can be modeled in a 1D reduced theory where an initial temperature drives the mixing effect.

Warm shear flow
We discuss here the temperature effect in a shear flow scenario, and its role in the generation of a dc magnetic field mode along the shear. For the sake of simplicity, and without loss of generality, we assume a simple sharp velocity shear transition between two plasmas with equal temperatures. We consider that the temperature is sufficiently high such that the electron thermal expansion time scale is much faster than the electron expansion induced by the onset of the fluid KHI in an equivalent cold scenario. The theoretical setting of the system is illustrated in Figure 12-a. We consider only the electron thermal velocity, neglecting the thermal velocity of the protons due to their inertia. Since we are interested in describing the dc phenomena, all derivatives along the x direction vanish, reducing the system to the 1D problem displayed in Fig. 12-b. We, therefore, consider the purely one-dimensional case where the particles can move along x, as in Figure 12-b. Initially all the fields are zero and we assume a warm initial plasma with a tangential shear flow identical to the one described in Sec.2.1 with an initial temperature such as v th v 0 . This setting is not in Vlasov equilibrium and it is clear that the thermal expansion of the electrons across the shear surface (ions are assumed to be cold and free streaming) leads to an imbalance of the current neutrality around the shear surface, forming a dc magnetic field in z direction. The initial corresponding electron distribution function reads The situation can be seen as two thermal plasmas with shearing counter propagating fluid velocities. The thermal expansion of the electrons (the ions, due to their inertia, are assumed free streaming) across the shear will transport an electron current on the order of en 0 v 0 , with a characteristic width of v thx t. This should lead, at early times, to the formation of a field B z around the shear of width v thx t and magnitude of µ 0 en 0 v 0 v thx t. This is the underlying physical picture of the dc magnetic field growth. Due to the dimensionality of the problem, it is clear that E z , B x , B y remain zero. The reduced set of equation is where The formal solution of the Vlasov equation Eq. (40) is where x 0 , v x0 and v y0 denote the position and velocities of an electron at t = 0 and f 0 = dv z0 F 0 . At early times, if we assume that the induced fields are sufficiently small that we can neglect the change of momentum of the electrons, the distribution can be solved along the free streaming orbits, i.e., For the sake of simplicity, we divide the initial electron distribution in two parts, , corresponding to the two initially separated flows. In the approximation of free streaming orbits, the electron currents read 2πv th represents the Maxwellian velocity distribution, we obtain for the electron currents The total current is obtained by adding the unperturbed proton currents, this yields The magnetic field is then simply given by the Maxwell-Ampere equation, Eq. (38), where the displacement current is neglected We verify then that the thermal expansion, transporting currents across the shear, induces a magnetic field that grows linearly with time. Its typical width, on the order √ 2v thx t and its peak B z (x = 0) = 2/πµ 0 en 0 v 0 v thx t, is in agreement with our previous estimates. The associated magnetic energy growing in the system is given by with du |u|erfc(|u|) − e −u 2 / √ π 2 0.156. As was pointed out, this derivation is only valid as long as the orbits of the electrons do not diverge much from the free streaming orbits, i.e., as long as the electric and magnetic fields that develop in the expansion process do not affect the free motion of the particles. In fact, there are two phenomena that affect the growth of the magnetic field. First, the electrons will eventually feel the induced magnetic field which tends to push more electrons across the shear via the v 0 × B z force. This will increase the rate of current transported and, thus, will increase the growth rate of the magnetic field. One should note that, at first, only the magnetic field acts on the electrons since the electric field E x remains zero in our model: the initial temperature is uniform in space and therefore, there are as many electrons crossing from the left than from the right. The charge neutrality is then conserved in the system. A crude estimate of the time at which our model breaks can be made considering that only the slow electrons, initially around the shear, will experience a strong velocity change due to the peaked shape of the ] 10 -2 10 -3 10 -2 10 -1 a4 b4 Figure 13. Evolution of the electron phasespace (insets a1-3 and b1-3) and dc magnetic field peak (insets a4 and b4); log-log is used to display linear dc peak evolution in inset a4, and log-linear is used to display exponential evolution in inset b4. Left: 1D warm shear flow with v 0 = 0.2c and v th = 0.016c. Right: 2D cold shear flow for the same v 0 . The blue (red) color represents the electrons with a negative (positive) drift velocity v 0 . The self-consistent dc magnetic field is represented by the solid curve, whereas the dashed curve represents the magnetic field given by the theoretical model. magnetic field. Therefore, the model will break down approximatively when an electron initially at rest (around x = 0) acquires a velocity change on the order of v thx , which corresponds to a strong distortion of the Maxwellian distribution function around the shear. This can be written for the velocity change as It follows that the model is valid until ω p t ∼ (2π) 1/4 c/v 0 . Second, the growth of the magnetic field induces an electric field E y through the Maxwell-Faraday equation Eq. (37). We can estimate through the characteristic time t and length of the problem v thx t, the magnitude of the field, E y ∼ v thx B z . Inserting the value of E y into the Maxwell-Ampere equation, we find the displacement current term ∂E y /∂t leads to a (v thx /c) 2 correction to the dc magnetic field peak at early times. However, the displacement current tends to increase the electron current on either sides of the shear interface, eventually building the dc magnetic field side wings observed at later times ( Fig. 13 a3-b3).
In order to verify our analytical calculations and to further investigate the phase at which the electrons deviate from their free streaming orbits, we have simulated a shear flow between warm electron-proton plasma slabs with a realistic mass ratio m p /m e = 1836 until the dc magnetic field structure saturates on the electron time scale. The Debye length is resolved in the 1D simulations (∆x = λ D ) and we have used 1000 particles per cell. These PIC results have been presented in [Grismayer and Alves (2013)] but are reproduced here for sake of completeness. Fig. 13 (a1-a3) shows the time evolution of the xp x phase space and the magnetic field for v th = 0.016c and v 0 = 0.2c. At earlier times ω pe t = 9 ( Fig. 13 (a1)) an excellent agreement between the model and the simulation is observed. According to Eq. 51, the model breaks down for times larger than ω pe t10. The deviation from the Maxwell equilibrium of the distribution function in the shear region is shown at ω pe t = 17. The model underestimates the magnitude of the magnetic field and one can clearly observe the distortion of the distribution function in the field region. As the magnetic field grows, the Larmor radius (r L ) of the electrons crossing the shear interface decreases. When the minimum r L,min (associated to the peak of B DC ) becomes smaller than the characteristic width of the magnetic field l DC , the bulk of the electrons becomes trapped by the magnetic field structure. This is illustrated in Fig. 13 (a3) at ω pe t = 55. The magnetic trapping prevents the electron bulk expansion across the shear (that drives the growth of the magnetic field), saturating the magnetic field. An estimate of the saturation can be obtained by equating r L,min ∼ l DC . From Eq. (49), it is possible to write the magnetic field as B DC (x, t) = 4πen 0 β 0 w(x, t), where w(0, t) should be interpreted as the characteristic width of the field. With l DC ∼ w(0, t), r L,min = mv 0 γ 0 /eB DC (0, t), we find that l DC ∼ c √ γ 0 /ω pe giving the saturation level of the magnetic field as Electron-scale shear instabilities: magnetic field generation and particle acceleration in astrophysical jets23 This scaling has been verified for 1D simulations (see Fig. 3 in [Grismayer and Alves (2013)]) for which the best fit function matching the simulation is eB sat DC /m e cω pe = 1.9β 0 √ γ 0 .

Cold shear flow
In the absence of an initial temperature, an alternative mechanism is needed to drive the electron mixing across the shear surface that, in turn, generates the dc field. This mechanism is the cold fluid KHI that has been thoroughly discussed in Sec. 2. In fact, in the warm shear flow scenario, both the cold fluid KHI and the electron thermal expansion can contribute to the generation of the dc field. This occurs when the typical length of the dc field due to the thermal expansion (l DC ) after a few e-foldings of the cold fluid KHI (T KHI−growth = N e−foldings /Γ max , where N e−foldings is on the order of 10) is on the order of the relativistic electron skin depth, i.e., v th T KHI−growth ∼ √ γ 0 c/ω pe .
Therefore, the cold fluid KHI dominates the electron mixing in the limit For a two dimensional cold plasma undergoing the KHI, the electron distribution function can be written as where v xf l , v yf l correspond to the velocity field solutions of the fluid theory. In this case, the self-generated KHI fields play the role of an effective temperature that transports the electrons across the shear surface, while the protons remain unperturbed, inducing a dc component in the current density, and hence in the fields. We then have to solve the evolution of the distribution function and show that the current density J y , averaged over a wavelength λ = 2π/k , has a non zero dc part. We follow the same approach as before and calculate the average distribution function defined as: To obtain analytical results we will assume that the linearly perturbed fluid quantities are purely monochromatic, which is equivalent to assume that after a few e-foldings, the mode corresponding to k = k max dominates with a growth rate of Γ = Γ max . We then write v yf l v 0 (x) and v xf l =v xf l sin(k y)e −k ⊥ |x|+Γt , wherev xf l , the amplitude of the velocity perturbations at t = 0, is associated to the small thermal fluctuations (small enough to ensure that the thermal expansion is negligible over T KHI−growth ). Inserting v xf l , v yf l into Eq. (55), we obtain where We observe that the development of 2D cold KHI reveals close similarities with the 1D hot model previously described. In the 2D KHI, averaging the distribution in the direction of the flow shows that the perturbation gives rise to a spread in v x that may be interpreted as an effective temperature. The spread in v x decays exponentially away from the shear and grows exponentially with time. The mean velocity is zero and the effective temperature associated to this distribution function is defined as One can then expect a similar physical picture as in the hot shear scenario and, as a result, the emergence of dc components in the fields which are induced by the development of the unstable KHI perturbations. The evolution of the phase space in Fig. 13 illustrates the similarity between the warm 1D (insets a1-3) and cold 2D (insets b1-3) scenarios. The challenge in this scenario is to determine how such a distribution function expands across the shear surface due to the complexity of the orbits in the fields structure (multidimensional fields with discontinuities at x = 0). One can however overcome this difficulty by solving the expansion along approximate orbits. This procedure, although not self-consistent, gives rich qualitative and quantitative insight regarding the features of the current that develops around the shear interface. In the xp x phase space, the electrons describe outward-spiraling growing orbits, since they are drifting across the standing growing perturbation. In the region where the electron mixing occurs, we assume electron orbits given by x ∼ x 0 + (v x0 /Γ)e Γt and v x ∼ v x0 e Γt where x 0 and v x0 are the position and velocity of a particle at the time t 0 when the instability begins.
where xΓ ∈ [−v 0 max , v 0 max ] and v 0 max (t) = v max (x = 0, t) that represents the maximum velocity of a particle that was originally in the vicinity of the shear. The limits of the integral in Eq. (59) represent the deformation of the boundary between the two flows on a characteristic distance of v 0 max /Γ as the instability develops. In the fluid theory, the boundary remains fixed, precluding the development of the dc mode. We then find the total current density by summing the proton contribution and integrating to obtain the induced dc magnetic field: with ζ = Γx/v 0 max . The peak of the dc magnetic field is located at x = 0 where the expression above reduces to B DC (0, t) = 8eβ 0 n 0 v 0 max (t)/(πΓ) and thus grows at the same rate as the KHI fields. One can verify in Fig. 13 (b1-b2) that Eq. (61) shows reasonable agreement with the 2D simulations. This derivation neglects the dc Lorentz force on the electron trajectories, which makes this model valid as long as the induced dc fields remains small compared to the fluid fields associated to the mode k max . The peak of the B DC field is proportional to v 0 max (t) that represents the maximum value of the fluid velocity v xf l , which obeys v xf l = −ie(E xf l + β 0 B zf l )/m e γ 0 (ω − k v 0 ). From the linear fluid theory, one can compute the ratio B zf l /E xf l from which we deduce that for subrelativistic shears (γ 0 ∼ 1), |v xf l | ∼ ec 2/7|B zf l |/m e ω pe v 0 implying B DC ∼ (8/ √ 7π)B zf l and that for ultra-relativistic shears (γ 0 1), |v xf l | ∼ e|B zf l |/m e c √ 2ω pe γ 3/2 0 yielding B DC ∼ (4/π)B zf l . We conclude that the induced dc magnetic field is always on the same order as the fluid fields and thus its consequences to KHI development cannot be neglected. As the dc field evolves, electrons start to get trapped and we expect a level of saturation similar to the saturating level obtained in the 1D model. This has been verified by the simulations. The comparisons between the saturation level of the 1D, 2D and 3D simulations are shown in Fig. 3 in [Grismayer and Alves (2013)], also verifying the β 0 √ γ 0 scaling. One can also compute the equipartition number related to the magnetic field at saturation. Using Eq. (52) one finds, which is similar to the equipartition number found for the Weibel instability [Medvedev (1999)]. Our derivation for the saturation of the dc magnetic field allows also to recover the empirical estimate of [Alves and Grismayer (2012)] already given for such a shear scenario. When a smooth shear is considered, the electron KH still develops as we have shown in Sec. 2. We verified that the initial electron transport across the shear, due the development of the instability, is the mechanism triggering the magnetic field generation, therefore validating the physics captured by our model. At saturation the dc magnetic field has a typical width on the order of the initial shear gradient length.
Keeping the same arguments we have used to derive the approximate value of the dc field at saturation, i.e., r L,min ∼ L implies withL = Lω pe /c √ γ 0 . Such a scaling has been verified for Lω pe /c 1 and the comparison between the crude estimate and the simulations is presented in Fig. 14. Interestingly, the dc magnetic field remains stable beyond the electron time scale and persists up to 100s ω −1 pi (see Fig. 1 in [Grismayer and Alves (2013)]). Eventually the protons will drift away from the shear surface due to the magnetic pressure, broadening the dc magnetic field structure and lowering its magnitude.

Particle acceleration
We investigate in this Section the acceleration of particles due to the development of electron-scale shear instabilities. Particle acceleration in shear flows has been previously investigated by many authors, mainly related to astrophysical scenarios [Berezhko & Krymskii (1981), Jopikii & Morfill (1990), Webb (1989), Ostrowski (1990), Rieger & Mannheim (2002), Rieger & Duffy (2005)]. The shear acceleration mechanism ([Rieger & Duffy (2005)]) is based on the idea that energetic particles may gain energy by systematically scattering off of moving small-scale magnetic field irregularities. These irregularities are thought to be embedded in a collisionless shear flow such that their velocities correspond to the local flow velocity. In the presence of velocity shear, the momentum of a particle travelling across the shear changes and the acceleration process essentially draws on the kinetic energy of the background flow. In the shear flows we discussed in the previous sections, fluid and kinetic effects lead to the emergence of organized electric and magnetic fields that are maintained in the shear region up to ion time scale. The electrons flowing in the shear region experience strong acceleration in these fields and also emit strong radiation while gyrating in the dc magnetic field. This acceleration process therefore differs from the shear acceleration mechanism of [Rieger & Duffy (2005)].
In order to investigate the acceleration of electrons in the shear due to the selfgenerated fields, we performed a 2D simulation of a relativistic cold shear flow, γ 0 = 3, v th /c = 10 −3 . The simulation domain dimensions are 250 × 2000(c/ω p ) 2 , resolved with 10 cells per electron skin depth (∆x 1 = ∆x 2 = 0.1c/ω p ) and 36 particles per cell per species are used. The shear flow initial condition is set by a velocity with +p 0 e 1 for x 2 > 0 and a symmetric flow with −p 0 e 1 for x 2 < 0. We impose periodic and absorbing boundary conditions in the x 1 and x 2 directions, respectively. The transverse direction x 2 has been extended up to 2000c/ω p with absorbing boundary condition, in order to avoid particle recirculation over the shear region which tends to produce unphysical additional acceleration. The dimension of the simulation box allows us then to follow the evolution of the system until ω p t ∼ 1000, time at which some particles approach the boundaries of the box.
The growth rate and fast growing mode in the early development of the relativistic electron-scale KHI was found to agree with the theoretical predictions of Sec. 2. As explained in Sec. 3, the nonlinear development of the instability drives the mixing between electron populations at the shear interface and generates dc components in the fields. At full saturation of the instability, the persistent electric and magnetic field components are on the order of √ γ 0 m e ω p c/e. During the early stage of the instability the oscillating fields are responsible for acceleration and deceleration of the electrons. This results in a slight temperature increase of the plasma and the electron distribution function widens around the mean energy γ 0 . Once the instability saturates, the electrons can experience strong acceleration due to the long-lived electric field structures in the shear region. Furthermore, the dc magnetic field, that only remains intense in the shear region, provides one of the mechanisms to curve the electron trajectories and hence takes part in the thermalization and isotropization of the electron distribution function. The energy distribution function of the electrons is plotted in Fig. 15 at at ω p t = 0 and ω p t = 1000. The final distribution can be separated in three distinct parts. The low energy part, 1 < γ < 5, corresponds to a thermal Juttner distribution, f (γ) ∼ γ 2 e −µγ with µ ∼ 1/γ 0 . The medium energy part exhibits a power law γ −5 up to γ ∼ 25. Finally after the elbow of the distribution one notices a hot tail that extends up to γ = 80.
To gain deeper insight into the acceleration process of the most energetic electrons in the shear region, we followed various electron trajectories whose final energy lied in the second and third part of the distribution shown in Figure 15. Figure 16 shows segments (on every inset the electron is tracked during ω p t ∼ 100) of two electron trajectories. The two electrons are initially in the vicinity of the shear, as shown in Fig. 16-a. The black and white background represents the total electric field magnitude. The growing middle structure represents the dc part of the field whereas the modulated patches on both sides originate from the saturated unstable modes. One can clearly see that the most of the acceleration occurs when the electrons cross the electric field patches. The magnitude of the electric field patches is mainly due to the transverse component E 2 , whereas the acceleration takes places in the x 1 direction. After being accelerated, the particle can cross the shear to be finally reaccelerated on the patches of the other side or definitively leave the shear region. The time evolution of the energy of the two tracked electrons is shown in Fig. 17. As explained previously, the energy of the electrons does no increase significantly until ω p t ∼ 300 which corresponds to the saturation of the KHI. Soon thereafter both electrons experience a strong energy kick, ∆γ ∼ 30 at around ω p t ∼ 50. At ω p t 430, the electrons have acquired their maximum energy that then remains constant, indicating they have left the field dominated shear region. The electron denoted by e 1 crosses the shear without really being affected by the transverse component of the electric field while the electron e 2 experiences a strong deceleration before eventually leaving the shear region.
One can understand the acceleration mechanism if one sees the process as relativistic electrons surfing on the electric field patches, flowing with the plasma that, can be assumed constant in magnitude. First it is necessary to quantify the magnitude of the three electromagnetic components. Using Eq. (15), Eq. (17) and Eq. (20), we obtain the following ratio in the limit γ 0 1: |E 1 /E 2 | ∼ 1/ √ 2γ 0 and |E 2 /B 3 | ∼ 1. If we assume that those ratio still hold during the non-linear phase, then the orbits of the electrons, approximatively streaming in the x 1 direction with the Lorentz factor of γ 0 , are mainly governed by E x2 and B x3 . The equations of motion read dp In the case E 2 = B 3 , the solution for the trajectory of the electron is given by [Landau & Lifshitz (1975)] with α = mc 2 γ − cp 1 and = mc 2 . One can clearly see that no matter what is the initial condition, for long enough time, p 2 /p 1 ∼ 2α/cp 2 which means that the momentum increases most rapidly in the direction perpendicular to E 2 and B 3 , i.e., x 1 . This confirms the observations from the simulations where the acceleration is mostly directed towards the x 1 axis where the component of the electric field is transverse. The maximum energy gain an electron undergoes while encountering an electric field patch is given by We treat the electric field as constant and we assume that its magnitude is on the order of its saturation value, E sat ∼ √ γ 0 m e ω p c/e. The maximum length of an electric field patch is typically 1/k max = 8/3γ 3/2 0 c/ω p . Since the patches are drifting with the initial flow, one must take into account the difference in speed between the electron being accelerated and the flow: ∆v = v e − v 0 . Assuming γ e γ 0 1 (γ e is the Lorentz factor associated to the electron speed v e ), we obtain For the parameters of our simulations, the cut-off of the spectra is typically on the order of γ 4 0 and in agreement with the energy spectra cut-off obtained in [Liang et al.(2013)] .The energy kick experienced by the electrons around ω p t = 400 in Fig.17 is about ∆E ∼ 40 for both electrons, which is smaller than ∆E max since they are encountering here smaller patches . As was previously illustrated with electrons e 1 and e 2 , the energy kick does not determine the final energy of an electron before it leaves the shear region. In fact, depending on its transverse momentum, the electron can either experience another acceleration or a strong deceleration. This erratic feature allows to understand the wide range of energies among which the electrons are distributed in Fig.  15.

Conclusions
Shear instabilities in plasmas are usually studied within the framework of magnetohydrodynamics where the plasma is considered as a magnetized fluid and where the typical time scale is governed by ion motion. We have shown in this work that electron scale physics leads to a variety of new effects when one considers an initially unmagnetized cold shearing collisionless plasma.
The collective electron dynamics can be described in first approximation by using a two-fluid model that allows to take into account electron inertial physics, not captured in MHD models. In this fluid framework, we have presented the derivation of the equations for the linear development of the longitudinal KHI, assuming arbitrary velocity and density profiles. This framework was applied to a special case, where the velocity and density profiles were simple step-functions, allowing analytical solutions to the equations. The model provided new insights into the effect of density-contrasts between shearing flows, namely that the development of the KHI is robust to density jumps, making it ubiquitous in astrophysical settings. We also observed that the unstable modes begin to drift when the density symmetry is broken. In large density-contrast regimes, the KHI dominates over other common astrophysical plasma instabilities such as the Weibel and Two-Stream instabilities. The case of a finite shear profile has also been investigated in detail. A smooth velocity profile (non step-like) induces a phase mixing of the eigenmodes of the system that results into a damping term. The combined effect of the instability due to the shear with the damping term suggests that the maximum growth rate is a decreasing function of the shear gradient length. All of these results obtained in the limit of linearized fluid equations have been accurately verified by 2D PIC simulations.
PIC simulation results have also demonstrated the emergence of a large-scale dc magnetic field after the onset of the electron-scale KHI . However this dc field is not captured by the two-fluid KHI theory nor MHD model. We have shown that the emergence of the dc magnetic field is intrinsically associated with electron-ion shear flows. The dc magnetic field is induced through the formation of dc current sheets driven by the expansion of electrons in the shear region due to either a thermal expansion or the development of the cold fluid electron-scale KHI perturbations. We have presented an analytical description of the formation of the dc field in agreement with 1D, 2D, and 3D PIC simulations. The dc magnetic-field saturation on the electron time scale is independent of the type of the expansion and persists beyond ion time scales.
Finally, we addressed the particle acceleration physics due to scattering in the selfgenerated fields of electron-scale instabilities triggered in unmagnetized shear flows. An understanding of how electrons are accelerated is essential if we are to fully interpret observations since it is presumably the radiation from energetic electrons that is most often observed from astrophysical sources. To address this issue, we have tracked the most energetic electrons in our simulations to identify the acceleration mechanism. It was found that the kinetic energy, initially stored in the drift, was mainly redistributed thermally. The bulk of the energy distribution energy of the electrons has a typical temperature comparable with the Lorentz factor of the flow. Nevertheless, the energy distribution also displays a non-thermal tail. The most energetic electrons, that make up the power-law tail of the distribution, are accelerated while surfing close to the speed of light on electric and magnetic field patches, self consistently developed by the electron-scale KHI, which are carried by the flow. This results in an efficient acceleration mechanism where electrons can reach energies on the order of γ 4 0 , where γ 0 is the initial Lorentz factor of the flow. In this appendix we present the detailed mathematics and numerics underlying the study of the effect of a finite shear gradient on the development of the electron-scale KHI, i.e., non-step like transition regions between the two flows. The equation that is required to solve in order to calculate the effect of a finite shear gradient on the dispersion relation is Eq. (22). For the sake of simplicity, we consider the case of a constant plasma density profile (n + = n − ). For a uniform density profile, Eq. (22) reads where the functions A, B and C are now given by where k 2 ⊥ = k 2 + ω 2 p /c 2 − ω 2 /c 2 , ω = ω r + iγ and v 0 is a function of x. The boundary conditions are such that the field E y1 should vanish at |x| → ∞. It is useful to obtain an integral formulation of the latter equation in order to relate the conditions of instability of the system to the function v 0 . To achieve this, let us first multiply Eq. (69) by the conjugate of E and integrate over x to obtain that can also be written as by separating the real and the imaginary parts of Eq. (71). If we look for an instability condition, one assumes that γ > 0 and it is then clear that (|∂ x E| 2 + k 2 ⊥ |E| 2 ) > 0, which implies that in Eq.(72) −γ 2 + (ω r − kv 0 ) 2 should be positive somewhere in order for the integral to vanish. To satisfy the second condition, i.e, Eq.(73), we see that the function ω r − kv 0 cannot be strictly positive or negative. For instance, ω r − kv 0 can be an odd function. If v 0 is an odd function, then ω r = 0. We will focus on the case of a symmetrical shear flow profile (ω r = 0). Returning to the first condition, Eq.(72), we see that k 2 v 2 0 > γ 2 somewhere and if |v 0 (x)| ≤ V 0 (V 0 the maximum absolute value of the flow velocity) and therefore the unstable modes should lie underneath the curves defined by γ = ±kV 0 . This is typically verified when one considers a tangential velocity shear, v 0 (x) = V 0 sgn(x), as in Section 2.1.1 where we obtained that Γ ≤ kV 0 .
If we return to Eq. (69) and we expand the first term in Eq. (69), the differential equation becomes This differential equation can be encountered in various physical plasma scenarios. It has been extensively studied in the case electrostatic oscillation in a non uniform cold plasma by [Barston (1964), Sedlacek (1971] and in the case of excitation of magnetic surface modes in MHD configurations by [Chen & Hasegawa (1974), Zhu & Kivelson (1988]. In these previous works, the frequency spectrum and the eigenfunctions of the differential equation have been calculated and three main regimes were identified depending on the function . When is an everywhere non-constant function, the spectrum is purely continuous and consists of those values of ω that satisfy the equation (x) = 0. In this case, the differential equation is always singular and the eigenfunctions can be constructed from well-known theorems on solutions of such equations near regular singularities. The well-behaved solution of the equation is obtained by an integral superposition over the whole spectrum of these eigenfunctions. On the other hand, as we already saw in Section 2.1.1, a plasma with a discontinuous velocity profile and therefore a discontinuous dielectric constant , exhibits different behavior. As a result of the analysis of [ Barston (1964)], the introduction of a jump discontinuity in (x) is necessary for the existence of well-behaved modes (so called surface modes) and the existence of a dispersion relation. The question that [Sedlacek (1971)] addresses is how to unite these two antagonistic pictures, which amounts to understand what happens to the surface wave when the velocity profile is smooth (or when the dielectric constant changes continuously from a region to another) and whether some collectives modes remain.
If the velocity profile is odd and varies continuously from −V 0 to V 0 , there is a local point in x at which the shear flow satisfies the resonant condition kv 0 (x) = ω r = 0. At this point the spatial resonance in the smooth profile modifies the growth rate of the collective modes. One can interpret the modification of the growth rate by the introduction of damping and spatial dispersion to the originally surface wave of a discontinuous plasma. Physically, at the resonant point ,the eigenmodes face phase mixing which results in a wave damping.
The theoretical framework to solve Eq. (74) was given by [Sedlacek (1971)]. Formally, one needs to construct the Green's function G(x, x , ω) of the differential equation where the mode solution of the electric field is given by where E 0 y1 (x , ω) denotes the initial perturbation and F the integration path that requires time causality and that is deformed around the singularities. The time asymptotic solution of the field comes from the contribution due to the singularities of the Green's function. More precisely, the contributions near the isolated singularities correspond to the collectives modes while the contribution of the integration along the branch cuts leads to the continuous spectrum [Sedlacek (1971)]. The general theory to obtain the Green's function of the differential Eq. (74) relies on a theorem demonstrated by [Friedman (1956)]. The Green's function can be expressed in terms of two linearly independent solutions ψ 1 and ψ 2 of the homogeneous Eq. (74), with ψ 1 satisfying the boundary condition at x = −∞ and ψ 2 at x = +∞. The Green's function is written as G(x, x , ω) = J −1 [ψ 1 (x, ω)ψ 2 (x , ω)H(x −x)+ψ 2 (x, ω)ψ 1 (x , ω)H(x −x)], (77) where H is the Heavyside function and J the conjunct of the two solutions defined by which is a function independent of the variable x (see [Friedman (1956)] for the details). It is clear that the singularities of the Green's function G are given by the values of ω for which J = 0. To obtain the new dispersion relation, one needs to therefore obtain an analytical expression for J and then find the zeros of this latter function. One problem that arises is that is usually unlikely that ψ 1 and ψ 2 can be expressed in terms of standard functions for a given dielectric constant (x), which poses difficulties in obtaining further analytical results. However, when the function (x) has a linear profile given by the eigenmode of this equation has been given already by [Sedlacek (1971)]. For k ⊥ L 1, where k ⊥ = k ⊥ (L = 0) as a first approximation, the dispersion relation is given approximately by 1 which reduces to the prior dispersion relation Eq. (29) when k ⊥ L → 0. From this crude description of the evolution of the dielectric constant in the transition region, one gets that the maximum growth rate decreases linearly with the small parameter k ⊥ L where Γ 0 max /ω p = 1/8 stands for the growth rate obtained for a discontinuous profile and where the second term in the r.h.s of Eq. (81), √ 3πk ⊥ L/8, can be interpreted as the damping factor arising from the finite shear gradient. The wavenumber corresponding to the maximum growth rate follows a similar trend by slightly decreasing when the parameter k ⊥ L increases.
Although the basic physical effects of the smooth transition on the growth rate have been discussed above for small values of k ⊥ L, one aims to calculate with greater accuracy the evolution of the maximum growth rate when the characteristic length of the velocity profile is varied. Since it is not possible to obtain analytical expression for J, a numerical solution of Eq. (74) is required. The underlying numerical algorithm can be outlined as follows. For a given wavenumber k, a guess value for ω is chosen. By following the evolution of ω as the characteristic shear gradient length is continuously varied, one can choose the analytical solution for the sharp shear case (equation (32)) as the initial guess value for ω; the initial guess value for higher shear gradient length should then be based on the previously calculated value for a smaller shear gradient length. For a given velocity profile, we numerically solve the differential equation Eq. (74) where the dielectric constant is evaluated for k and the guess value ω. The numerical solution of the electric field is then injected into the integral equation Eq.(71), where we look now for the value of ω such that the integral of Eq.(71) vanishes. We then iterate the process until the value of ω converges.