Nonlinear evolution of the magnetized Kelvin-Helmholtz instability: from fluid to kinetic modeling
Abstract
The nonlinear evolution of collisionless plasmas is typically a multi-scale process where the energy is injected at large, fluid scales and dissipated at small, kinetic scales. Accurately modelling the global evolution requires to take into account the main micro-scale physical processes of interest. This is why comparison of different plasma models is today an imperative task aiming at understanding cross-scale processes in plasmas. We report here the first comparative study of the evolution of a magnetized shear flow, through a variety of different plasma models by using magnetohydrodynamic, Hall-MHD, two-fluid, hybrid kinetic and full kinetic codes. Kinetic relaxation effects are discussed to emphasize the need for kinetic equilibriums to study the dynamics of collisionless plasmas in non trivial configurations. Discrepancies between models are studied both in the linear and in the nonlinear regime of the magnetized Kelvin-Helmholtz instability, to highlight the effects of small scale processes on the nonlinear evolution of collisionless plasmas. We illustrate how the evolution of a magnetized shear flow depends on the relative orientation of the fluid vorticity with respect to the magnetic field direction during the linear evolution when kinetic effects are taken into account. Even if we found that small scale processes differ between the different models, we show that the feedback from small, kinetic scales to large, fluid scales is negligable in the nonlinear regime. This study show that the kinetic modeling validates the use of a fluid approach at large scales, which encourages the development and use of fluid codes to study the nonlinear evolution of magnetized fluid flows, even in the colisionless regime.
pacs:
I Introduction
In typical laboratory, space and astrophysical conditions, the nonlinear dynamics of magnetized plasmas is driven by the energy injected at large, fluid scales. The energy then cascades self-consistently towards smaller and smaller scales, until kinetic effects come into play. From a theoretical/modeling point of view, space plasmas represent a laboratory of excellence to study the physics of fundamental plasma processes, because of the wealth of in-situ diagnostics of improving quality accumulating in the form of electromagnetic profiles and particle distribution functions. The plasma turbulent state, routinely observed by satellites in the solar wind, is an archetype of plasma multi-scale behavior Bruno and Carbone (2005).
Many plasma processes naturally lead to a multi-scale dynamics, as for example at the interface between two different plasma regions, where large scale, fluid instabilities develop self-consistently and act as an energy source. This is the case for instance of the solar wind-magnetosphere interface, which plays a key role in the context of space weather modeling and forecasting. The connection between the solar wind and magnetosphere is mediated through the magnetosheath and magnetopause boundaries, that strongly depend on the solar wind properties and their variability. At the transition region between the solar wind flowing plasma and the Magnetosphere plasma at rest at low latitude, nearby the equatorial plane, the velocity shear between the two plasmas is an efficient source for the development of the Kelvin-Helmholtz instability (hereafter KHI).
Satellite measurements have supplied clear evidence of rolled-up vortices at the flank of the Earth’s Magnetopause Fairfield et al. (2000); Hasegawa et al. (2004) where the KHI has been invoked to provide a mechanism by which the solar wind enters the Earth’s magnetosphere, in particular to account for the increase of the plasma transport. Indeed, during northward magnetic field periods, when magnetic reconnection is considered as inefficient, at low latitude where the solar wind and magnetospheric fields nearly parallel, a relevant mixing between the solar wind and the magnetospheric plasma is observed, even larger than during southward configurations. KHI driven by a velocity shear can grow at low latitude since the magnetic field, nearly perpendicular to the plane where the instability develops, does not inhibit its development. This provides an efficient mechanism for the formation of a mixing layer and for the entry of the solar plasma into the magnetosphere, explaining the efficient transport during northward solar wind periods. Other space observations have provided support for the development of the KHI in the environment of Saturn Masters et al. (2010) and Mercury Sundberg et al. (2011, 2012).
After saturation, the KH rolled-up vortices drive the formation of gradients at the ion inertial length and/or the ion Larmor radius up to electronic scales, and act as a source of secondary instabilities. In these conditions, previous studies have demonstrated that the nonlinear evolution of KH vortices enables the occurrence of magnetic reconnection driven by large-scale vortex motions, of interest for space plasmasFaganello, Califano, and Pegoraro (2008a, b); Faganello et al. (2010); Henri et al. (2012), as well as for astrophysical plasmas jetsKeppens et al. (1999); Baty, Keppens, and Comte (2003). The vortex formation process indeed drags the magnetic field component parallel to the solar wind direction into the flow. As a result, the magnetic field is more and more stretched inside the vortices until it reconnects, redistributing the initial kinetic energy into accelerated particles and heating. Moreover, the density jump between the magnetosheath and magnetospheric plasmas drives fluid-like secondary instabilities such as the Rayleigh-Taylor instability Matsumoto and Hoshino (2004); Faganello, Califano, and Pegoraro (2008c); Tenerani et al. (2011). On top of that, the downstream increase of the magnetosheath velocity leads to super-magnetosonic regimes for which the KH vortices act as obstacles to the plasma flow, generating quasi perpendicular shocks Palermo et al. (2011a, b); Henri et al. (2012), thus modifying the transport properties of the plasma of the solar wind-magnetosphere. It is thus crucial to establish the role of these different secondary instabilities on the dynamics of the system, since they strongly influence the increase of the width of the mixing layer and its internal dynamics that are the most important factors for the evolution at the flank of the Earth’s Magnetosphere.
To summarize, the nonlinear evolution of the large scale fluid vortices is a fundamental plasma physics process driving the development of secondary instabilities. The energy then self-consistently cascades towards smaller and smaller scales, where the dynamics become kinetic, playing a significant role in the transport properties of the global system.
Even if numerical studies of the nonlinear evolution of magnetized shear flows have been carried mainly by the means of fluid models (ideal/resistive MHD, Hall MHD, two-fluid), the increase of computational power has recently enabled to adress the problem of the kinetic modelling of shear flows in collisionless plasmas through hybrid PIC Cowee, Winske, and Gary (2009, 2010) and full PIC simulations Wilber and Winglee (1995); Nakamura, Hasegawa, and Shinohara (2010); Nakamura et al. (2011). Low resolution simulations of the KHI have also been computed as a test problem to benchmark Vlasov codes Umeda et al. (2010). One of the main difficulty in the kinetic modeling of shear flows however remains the choice of the initial conditions. Indeed few kinetic equilibriums are known for shear flow configurations Cai, Storey, and Neubert (1990).
To progress beyond the current state of the art in the nonlinear study of collisionless plasma evolution, the primary kinetic effects and physical processes at play should be understood through the complementary use of different models, from fluid to kinetic. At the end, only the kinetic modeling can validate or not the choice of a fluid approach. This is why a multi-model study is necessary to shade light on the fluid modeling of the nonlinear evolution of colisionless plasmas.
We decided to focus in this work on the comparison of the evolution of a magnetized shear flow through the development of the KHI, using different plasma codes/models. However, we underline that the numerical modeling of plasma dynamics is a fundamental problem of major interest in present plasma physics research. Therefore, this study is of broad interest and is not limited to the KHI itself and related nonlinear dynamics.
The paper is organized as follows. The different models and codes used in this study are described in section II and the configuration of the system under study is described in section III. The results are presented in section IV: first, the kinetic relaxation to the initial fluid sheared flow equilibrium is discussed in section IV.1, then the linear growth rates of the magnetized Kelvin-Helmholtz instability are compared in section IV.2, the nonlinear phase of the magnetized Kelvin-Helmholtz is then studied for the different models in section IV.3. We finally discuss and conclude this work in section V.
Ii Models and codes
The primary goal of this comparative study is to use different plasma models to solve a given physical problem, namely the evolution of a shear flow in a magnetized plasma. Multiscale properties of magnetized plasmas arise when nonlinear processes, hardly modeled by analytical studies, are at play. This is why numerical tools are needed to efficiently integrate in time the equations describing the plasma dynamics. On the other hand, benchmarking different numerical algorithms of the same model is out of the scope of this paper.
ii.0.1 Hierarchy of plasma models
Model | Closure | Closure | Length | |
on ions | on electrons | scales | ||
(i) | MHD | Isothermal/Adiabatic | ||
(ii) | Hall-MHD | Adiabatic | , | |
(iii) | Two-fluid | Adiabatic | Adiabatic | ,, |
(iv) | Hybrid | no closure | Isothermal | ,, |
(v) | Full kinetic | no closure | no closure | ,,,, |
The multiscale intrinsic nature of collisionless magnetized plasmas has lead to the development of a variety of plasma descriptions starting from a N-body description, through a mean field kinetic Vlasov description, to fluid descriptions (multi-fluid and magnetohydrodynamic). The underlying idea is that unresolved phenomena at small scales can be averaged out, leading to models much easier to handle than a N-body description.
In this work, we have restricted our study to the following models: (i) magnetohydrodynamic (MHD), (ii) Hall-MHD, (iii) two-fluid, (iv) hybrid model with kinetic ions and massless electrons fluid, (v) full kinetic description of ions and electrons.
The full kinetic model (full PIC code) contains all the relevant plasma length scales: the ion and electron inertial lengths and , respectively, as well as the ion and electron Larmor radius and , respectively. The hybrid model (hybrid PIC code) contains the ion length scales and , while the two-fluid model contains the ion and electron inertial lengths and . Ideal and resistive MHD are both transparent to the plasma lengths scales. This is summarized in table 1, third column.
Concerning the closure of the models (summarized in table 1), the fluid codes (MHD, Hall-MHD and two-fluid) use either an isothermal or an adiabatic closure or respectively, the hybrid code uses an isothermal closure for electrons, while no closure is needed for the species described by a kinetic model: this is case for ions for the hybrid PIC code and for electrons and ions for the full PIC code.
We stress that even if the kinetic models (full PIC in section IV.2 and hybrid PIC in section IV.3) should be considered the ”models of reference”, the initial equilibrium setup is chosen fluid-like, because of the difficulty of finding a Vlasov equilibrium, so that it does not describe a kinetic equilibrium for these two models. This is why the inter-model comparison should be done carefully, taking into account the fact that the force equilibrium can readjust during the initial phase.
The numerical algorithms used to solved the different plasma models are described in the next sections. (i) Two different codes are used to solve the MHD equations, they are described in sections II.0.2 and II.0.3, (ii) and (iii) the code solving the Hall-MHD and the two-fluid models is described in section II.0.4, (iv) the hybrid and (v) full kinetic models are solved using two different PIC algorithms described in sections II.0.5 and II.0.6 respectively. These different codes have all been extensively tested and previously used to produce peer-reviewed scientific publications, so their validation is assumed and out of the scope of this paper.
ii.0.2 MHD Stagger code
One of the MHD codes used in this work, the ‘Stagger Code’, is a grid based, resistive and compressible MHD code. The code incorporates an adaptive hyper-resistivity and -viscosity scheme for enhanced control of dissipation introduced in fastmode waves, advective motion and shocks — see for example Baumann, Galsgaard, and Nordlund for an overview of the dissipation scheme and further code features. The MHD variables are defined on staggered grids, and the code conserves mass, momentum and to machine precision. Interpolation of variables between different staggered grids is handled by using 5th order interpolation. Spatial derivatives are computed using 6th order accurate differential operators. The time integration of the MHD equations is performed using an explicit 3rd order low storage Runge-Kutta procedure Williamson (1980). For the MHD simulations conducted here, the resistive MHD equations solved in the Stagger Code are:
(1) | |||||
(2) | |||||
(3) | |||||
(4) | |||||
(5) | |||||
(6) | |||||
(7) | |||||
(8) |
Here, is the mass density, the bulk velocity, the pressure, the current density, the magnetic field, the thermal energy. is the viscosity, is a thermal conductivity, and is the resistivity. is the velocity shear tensor, the viscous stress tensor. is a diffusive heat flux, needed for numerical stability. From Eqs. (14) & (15) in Ref. Baumann, Galsgaard, and Nordlund using parameters , and , the resistivity is for both the adiabatic and isothermal cases with the setup described below. This initial resistivity will adjust during the instability growth and deviate from an approximate constant. This hyper-resistivity/-viscosity scheme will consequently yield larger effective values of the resistive and viscous terms, locally on the grid, when and where numerical critical structures appear. For the simulations with isothermal conditions, the energy equation is not active in the code and (with ). The ideal gas law with is assumed for the adiabatic case.
ii.0.3 Mpi-Amrvac
The second code used in this work is the MHD module of the MPI-AMRVAC software Keppens et al. (2012). The code uses a finite-volume discretization, combining explicit multi-step timestepping schemes with a variety of shock-capturing spatial discretizations. The code is designed to solve equations of generic form
(9) |
and in this work, we employ two variants of the physical equations, which are (1) the isothermal, ideal MHD equations where conservative variables include density, momentum and magnetic field; and (2) the full set of ideal MHD equations where in addition a total energy density equal to is evolved. For the isothermal case, the pressure appearing in the momentum equation is at any time set from in dimensionless units, such that it maintains the same equal uniform temperature conditions as initially in the full MHD case. The ratio of specific heats is set to for the latter run. Although the code allows inclusion of physical sources S, such as viscosity or resistivity, in this work, we omit any explicitly added physical source that would introduce deviations from the pure conservation form. This means that any reconnection of the magnetic field is entirely due to the nonlinearities in the shock-capturing discretizations. To make the comparison possible with all other fixed grid simulations, we use MPI-AMRVAC in a domain decomposition approach, using only 1 grid of overall size , but employing 3072 grid blocks of size to exploit the MPI parallelism. The actual scheme used is an overall third-order Total Variation Diminishing (TVD) Lax-Friedrichs scheme, using a three-step Runge-Kutta variant combined with the third order C̆ada limiter Čada and Torrilhon (2009). The control is handled through a diffusive approach, introducing a (non-physical) source term diffusing monopole errors, which is handled in a split fashion. For the explicit time-advance, a courant parameter of 0.9 is used throughout.
ii.0.4 Two-fluid code
The Two-fluid codeFaganello, Califano, and Pegoraro (2009) is based on a two-fluid, ion-electron plasma approach including electron inertia effects in a fluid framework. The dimensionless equations of the model are obtained by using ion characteristic quantities: the ion mass , the Alfvén velocity and the ion inertial scale length (in dimensionless units the electron inertial length reads ). The density and motion equations read:
(10) |
(11) |
where is the fluid velocity, is the total pressure, with the thermal pressure, and the density (we assumed quasi-neutrality). We adopt adiabatic closures for ions and electrons:
(12) |
where . In the following we take . The dimensionless sound velocity is defined as . The electric field is calculated by means of a generalized Ohm’s law including electron inertia effectsValentini et al. (2007):
The two-fluid code can also be run as a Hall-MHD code by imposing . Finally the magnetic field is calculated by solving the Faraday equation (Eq. 2) and the current is given by (we neglect the displacement current). These equations are integrated numerically in a 2D () slab geometry , using fully 3D fields (the so called 2.5D geometry). This code is based on a standard third-order Adams-Bashforth method for temporal discretization. It uses fast Fourier transform routines for spatial derivatives along the periodic -direction and sixth-order compact finite difference scheme with spectral like resolution for spatial derivative along the inhomogeneous -direction (Lele, 1992). Numerical stability is achieved by means of filters, a spectral filter along the periodic -direction and a sixth-order spectral-like filtering scheme along the inhomogeneous -direction (Lele, 1992).
ii.0.5 Hybrid PIC code
The hybrid code is based on a current advance method and a cyclic leapfrog algorithm Matthews (1994). In the model, ions are treated by using a particle in cell scheme while electrons are represented by a massless, isothermal, charge neutralizing fluid. The code self-consistently solves equations of motion for ions
(13) | |||||
(14) |
together with Faraday’s law for magnetic field (Eq. 2) and a generalized Ohm’s law for electric field in the form
E |
here , are positions and velocities of particles of species , and are total ionic charge and current densities, B is the magnetic field, is vacuum permeability and is a resistivity parameter. The pressure of electrons is obtained as where the electron temperature is set as initial condition and remains constant during the simulation (the electron fluid is isothermal). The electron density is computed from the total ionic charge density using the assumption of quasineutrality , here stands for elementary charge.
ii.0.6 Full PIC code
iPIC3D is a fully kinetic, fully electromagnetic Particle-in-Cell code (Markidis, Lapenta, and Rizwan-uddin, 2010). It implements the moment implicit method (Brackbill and Forslund, 1982). In this code, the second order formulation of Maxwell’s equations for the electric field is discretized implicitly in time:
(15a) | |||
Once the electric field is calculated, the magnetic field is advanced in time solving the discretized Faraday’s law: | |||
(15b) |
The Maxwell’s equations are differenced in space on a uniform cartesian grid and the simple box scheme is used for the spatial differentiation of spatial operators in the field equations (Eqs.(15a)-(15b)). The linear system arising from Eqs.(15a)-(15b) is solved using the Generalized Minimal RESidual (GMRes) solver. The equations of motion of particles are differenced in time using the implicit midpoint integration rule:
(16a) | |||||
(16b) |
, and are the electric and magnetic field, calculated at the midpoint of the orbit and is the average of and . The iPIC3D code solves the particle equation of motion (Eq.(16)) by an iterative method based on a fixed number of predictor-corrector iterations.
Iii Description of the simulations setup
The different models and corresponding codes previously described are used to integrate the linear and nonlinear evolution of a magnetized shear flow plasma unstable to the KHI (section IV.2 and IV.3 respectively).
The simulation setup is identical for the different models. We consider a 2D physical space, with 3D vector fields and an initial MHD equilibrium (note: not a Vlasov equilibrium). In order to avoid additional/spurious effects due to different implementations of boundary conditions, a double shear layer is considered in order to impose periodic boundary conditions. The initial flow configuration is shown in Fig. 1.
In the following, all quantities are normalized to ion quantities: the ion gyro-frequency, , the ion inertial length, , and the Alfvén velocity . The size of the numerical box is given by and . The number of grid points in the XY-plane of the simulation is and . We use 1024 particles per cell in the hybrid PIC code and 200 particles per cell in the full PIC code.
The initial velocity field contains a periodic double shear layer where the velocity varies from to . The shear layers are located in and . The velocity profile reads:
where the maximum velocity field strength is and the shear length scale is given by (in units), which implies in terms of the electron inertial length for the two-fluid and full-PIC simulations, using a mass ratio . The initial current is taken to be zero, . The initial magnetic field is , where and , so that .
Note that the (mostly out-of-plane) magnetic field is constant, always pointing in the same direction. On the contrary, the direction of convective electric field varies from one layer to the other: it is directed towards (resp. away from) the shear layer at (resp ). Note also that the vorticity is parallel to the out-of-plane magnetic field (, where ) at , while it is anti-parallel () at .
The initial electron and ion pressures are isotropic, with , corresponding to a total thermal pressure . The initial density is constant in the simulation box, , so that the electron and proton temperatures are equal, . Quasineutrality is imposed at the beginning of the simulation for the full PIC code, while it is assumed in the other models. Finally, in the two-fluid and full kinetic models, the proton-to-electron mass ratio is .
An initial incompressible perturbation on the velocity field is imposed in the fluid models (not needed in the hybrid and full PIC codes) as follows:
with
where
and is such that and are random phases.
When studying the linear phase of the KHI (section IV.2), random phases are considered that may vary from one simulation to the other, since this does not modify the evaluation of the linear growth rate. On the contrary, identical random phases are considered for all runs, but different from one layer to the other, when studying the nonlinear phase of the instability (section IV.3), since the initial phases influences the nonlinear evolution of the full system.
The plasma beta is , so that the ion inertial length and the ion gyroradius are equal. The setup has been chosen so that the initial shear length is close to the ion inertial length and the ion gyroradius, in order to study the transition between the fluid and kinetic regimes. The grid length has been chosen so that the fastest growing mode is m=2, this way two rolled-up vortices form in the simulation box.
Iv Comparative study
The kinetic relaxation oberved with kinetic models is first studied in section IV.1. We then discuss the linear part of the KHI in section IV.2 and its nonlinear part in section IV.3.
iv.1 Initial kinetic relaxation
In both kinetic simulations (hybrid and full PIC), the initial profiles are (asymmetrically) modified at the velocity shear location on the ion gyration time scale, before the KHI develops. This is due to the kinetic relaxation in response to the initial fluid (not Vlasov) equilibrium Cai, Storey, and Neubert (1990). This phenomenon is particularly evident in the case and becomes stronger as the plasma beta increases. This effect, illustrated by the average mean (fluid) velocity shear computed for times (Fig. 2, top panels), may explain part of the differences between the fluid and kinetic linear growth rates, since the fastly modified initial shear profile will respond differently to the imposed fluid profile.
Furthermore, Fig. 2 also shows that the shear layer is modified more significantly in the full kinetic simulation than in the hybrid simulation. To quantify this feature, the layer with broadens from the initial (fluid) thickness to and in the hybrid and the full (2) kinetic simulations respectively. While the layer with broadens to and in hybrid and full kinetic simulations respectively. Note that both simulations are computed for the same ion beta . Here is computed as the distance from the layer center ( or ) at which the velocity field reaches the value .
This effect on the mean velocity goes together with a modification of the plasma density nearby the center of the velocity shear layer, due to finite Larmor radius effects. A ‘bump’ for (resp. a ‘dip’ for ) in the plasma density is generated in both the hybrid and the full PIC simulations, as shown in Fig. 2, bottom panels, where the blue and red lines respectively show the hybrid and the full PIC simulations and the black dashed lines show the fluid (MHD and two-fluid) profile.
Since the initial condition is not a Vlasov equilibrium, a disturbance on the plasma density is excited and travels at constant speed through the simulation domain along the -direction, as shown by two kinetic simulations (hybrid PIC and full PIC) in Fig 3. This localized travelling perturbation, of magnetosonic nature, excited at a given shear layer eventually reaches the other layer and may modify the dynamics of the system. This artificial effect is further increased by the double periodic boundary conditions that do not allow to get rid of these density perturbations. We do not observe in our hybrid and full simulations a clear relaxation of this disturbancy.
Still considering the growing phase of the instability, we observe an initial deformation of the distribution functions in both the full and hybrid PIC simulations, which is interpreted in terms of the kinetic relaxation from an initial setup that is an MHD equilibrium but not a kinetic equilibrium. At the velocity shear layers location, the initially gyrotropic plasma becomes agyrotropic after a few gyro-periods (so that temperatures in and directions are different), as illustrated in Fig. 4 where we show the result of the full PIC simulation. Note in the right panel the generation of a temperature anisotropy: corresponding to the bottom layer . The agyrotropy leads to the modification of the full pressure tensorCerri, S.S. et al. (2013) which is properly included only in the kinetic models. It is unclear whether this behavior affects significantly the evolution of the instability. Note also that the ion tracer show a different thickness in the two different shear layer (central panel), due to the cumulative effect of the (fluid) velocity shear and (particle) ion gyromotion Nakamura, Hasegawa, and Shinohara (2010).
iv.2 Linear growth rate
The different codes first run with the previously described initial conditions until the end of the linear phase of the instability.
In Fig. 5 we show the growth rates for the first three unstable modes, computed using the different models, with the growth rate of the KHI (over-plotted with a black continuous line) calculated from the linear MHD theory with an adiabatic closure, using the linear MHD LEDAFLOW code Nijboer et al. (1997). In the figure the cross and square symbols refer to the relative orientation of the vorticity with respect to the magnetic field. The fastest growing mode in this setup is the mode , which corresponds to a wavenumber . The different models show a good agreement at large wavelength (mode m=1) and differ at smaller wavelengths (modes m=2 and m=3). This is a first signature of the expected disagreement between different plasma descriptions, when kinetic and ion inertial scales are reached.
Hereafter, we will refer to the MPI-AMRVAC code (section II.0.3) to model ideal MHD, and to the Stagger code (section II.0.2) to model resistive MHD. The growth rates obtained using the two MHD codes with an adiabatic closure are the same and identical to the growth rate resulting from the linear MHD (adiabatic) theory. In a compressible fluid, part of the available energy also goes to compression. This is why the growth rate is slightly smaller when the closure is isothermal, since the compressibility is larger.
The MHD, Hall MHD and two-fluid models show very similar results during the linear stage of the instability, before the nonlinear dynamic of the vortices create gradients at scales of the order of the ion inertial length. The growth rates obtained from the Hall-MHD and two-fluid simulations are however slightly larger than those obtained from the MHD simulations. The Hall effect modifies the effective magnetic tension since the magnetic field lines are now frozen-in in the electron fluid and no more in the ion fluid. Thus, the stabilizing effect of the magnetic tension is reduced in the presence of the Hall term, explaining why the KHI growth rate is slightly larger when using the Hall-MHD or the two-fluid model.
Squares (resp. crosses) are used in Fig. 5 to represent the growth rate computed at the shear layer where (resp. ). The relative orientation of the velocity field vorticity compared to the background magnetic field does not modify the KH growth rate in the MHD, Hall-MHD and two-fluid regimes, as shown in Fig. 5. Indeed, the fluid models used in this study do not contain intrinsic information on the particle gyration. On the contrary, models that include ion kinetic effects (Hybrid and full PIC codes represented by yellow and red lines in Fig. 5) show that the relative orientation of and modifies the linear dynamics of the KHI. Note that different fluid models that take into account Finite Larmor Radius (FLR) corrections in a fluid description Passot, Sulem, and Hunana (2012) are expected to exhibit the same asymmetry in the KHI growth rate, observed here with hybrid and full kinetic models, as shown in previous studies with FLR-MHD modelsHuba (1996). In particular, we observe that both hybrid and full PIC calculations show the same qualitative feature: the growth rate for is smaller than that for , as observed also in previous PIC simulations Wilber and Winglee (1995); Nakamura, Hasegawa, and Shinohara (2010). This result is in contrast with previous analytic calculations including FLR effects in an extended MHD description Nagano (1978, 1979), where the growth rate was found to be larger for than for . It is particularly striking that all kinetic simulations exhibit the contrary. This is most probably due to the lack of initial kinetic equilibrium that implies a fast readjustment of the kinetic plasma which modifies the effective setup before the instability starts (see section. IV.1). The discrepancies observed in the growth rate computed from the hybrid and full kinetic descriptions are most probably due to the treatment of electrons as a massless isotropic isotermal fluid in the first case and kinetically in the second case. In this setup, electrons and ions have the same initial temperature so that the total pressure is initially equally divided between them. In this case, the electron compressibility is also expected to play a key role. The high level of noise in PIC codes, which implies a greater uncertainty on the calculation of the KHI growth rates, could also account partly for the growth rate discrepancies using the kinetic models.
Note that the ion inertia should also affect the growth rates depending on the relative orientation of and Fujimoto and Terasawa (1991). However, at least in the range of parameters used in this study, the Hall term do not give any measurable difference in the growth rate between the and cases. On the other hand, the hybrid and full kinetic models show different growth rates for different orientation of the vorticity and the magnetic field. Therefore, since in these models both the Hall and the FLR effects are included, we conclude that the difference observed in the growth rates for the two orientation is dominated by the FLR effects.
Interestingly, including ion inertial effects (spatial scale ) in the model increases the growth rates, as seen when comparing (adiabatic) MHD growth rates to (adiabatic) Hall-MHD and two-fluid growth rates. On the contrary, including ion kinetic effects by means of a kinetic code (in particular the ion gyroradius ) decreases the growth rates whatever the relative orientation of and .
The observed difference between fluid and kinetic growth rates could be a consequence of different effects: (i) the absence of initial kinetic equilibrium, which could induce significant modifications in the shear layer, (ii) a difference in the plasma compressibility between the fluid and kinetic models, (iii) the influence of finite ion Larmor radius effects.
iv.3 Nonlinear evolution
In this section, we consider the following models: isothermal ideal MHD, isothermal resistive MHD, adiabatic ideal MHD, adiabatic resistive MHD, two-fluid and hybrid models. The full kinetic (PIC) model is not used in this part because of the huge computational time required to model electrons kinetically. A full kinetic simulation of the KHI is a computational challenge that will be tackled in future works. To properly compare the nonlinear evolution of the KHI obtained from different models, we consider identical initial velocity perturbations in the different runs by imposing the same first random phases , , however different from one layer to the other. The initial amplitude of the perturbations is now set to .
Comparative nonlinear evolution at large scales
The nonlinear evolution of the KHI is shown in Figs. 6 and 7 (layers centered at and , characterized by and respectively), through the evolution of the density and of the magnetic field lines. For sake of clarity, both layers are plotted separately on different pictures, Figs. 6 and 7 respectively; however we recall that the simulations are periodic in both directions. Each row corresponds to a single model, while time evolves from left to right.
At the beginning of the nonlinear phase (, first column in Figs. 6 and 7), two fully rolled-up vortices develop along each shear layer. The differences between the two shear layers are mainly due to the different phases used in the initial perturbations. At this early nonlinear stage, before smale scale processes develop, no appreciable discrepancies arise between the different models for a given shear layer.
The classical evolution of the hydrodynamic KHI, in the absence of magnetic field, is an inverse cascade characterized by the merging of the different vortices generated from the fastest linear growing mode, up to the formation of a single large vortex in the numerical boxBaty, Keppens, and Comte (2003). In the magnetized case reported here, the onset of vortex pairing is shown at (second column in Fig. 6), on the first shear layer, , with . The merging from two rolled-up vortices to a single one then develops at successive times and (third and fourth columns). On the contrary, in the case (Fig. 7), no vortex merging is observed. Instead, at (second column), the vortices remain aligned with the initial shear layer direction at . At later times, the onset of vortex pairing is observed at (third column) but small scale processes disrupt the vortices before the pairing is completed, leading to the formation of a mixing layer instead of a large single vortex.
To summarize the global non linear dynamics observed in the two different layers, we observe that in the first layer the final state is given by a large scale, coherent single vortex, while in the second layer, the final state is given by a turbulent layer. Both cases are well captured at large scale by all the considered models, from ideal MHD to hybrid kinetic.
In Figs. 6 and 7, the -direction is represented in order to give a better feeling on the matching (or mismatchings) between the successive models. We recall that the numerical box is periodic. Note that the different models perfectly connect one another at the early stage of the KH nonlinear evolution (first column) and connect well at large scales during the nonlinear evolution, showing the agreement between all models at large scales. Note also that some discrepancies arise at small scale in the late stage of the KH nonlinear evolution (last column in particular), as can be seen when looking at the ‘boundaries’ between two successive models.
Comparative nonlinear evolution at small scales: the essential role of magnetic reconnection
At scales much smaller than the hydrodynamic scale , the magnetic field is stretched and compressed by the plasma motion, forming current sheets between and inside the vortices, at scales of the order or smaller than the ion inertial length (Figs. 8 and 9). In the ideal and resistive MHD simulations, strong current sheets are formed but magnetic reconnection does not start at early stage. On the contrary, the same current sheets modeled by the two-fluid and the hybrid models become tearing unstable and develop magnetic reconnection. Such a difference at small scales is due to the Hall effect, included in the two-fluid and the hybrid models, but absent from the ideal and resistive MHD models. The Hall effect allows for magnetic reconnection to develop at a faster rate in the two-fluid and the hybrid models, while the MHD models rely on the explicit (resistive MHD code) or the numerical (ideal MHD code) resistivity for reconnection to occur Birn et al. (2001).
In the nonlinear phase of the magnetized KHI, two regions are typically subject to reconnection: (i) the compressed region that separates two merging vortices and (ii) the sheared magnetic regions inside the vortices themselves. The Hall effect enables the tearing instability to develop faster in both regions: (i) between the two pairing vortices, as seen by the X-points and magnetic islands formed at at (first column) in Fig. 8 and at at (second column) in Fig. 9 for the two-fluid and hybrid models (bottom two panels) in contrast to the surviving magnetic shear in the two MHD models (top two panels) at the same times; (ii) inside the single vortices, as seen by the chains of magnetic islands (e.g. at , at (second column) in Fig. 8 and at at (second column) in Fig. 9. Note that vortex-induced reconnection still occurs in resistive MHD in both regions but at later times, , at and at in Fig. 9. In the hybrid simulation, the reconnection process seems to be triggered earlier with respect to the two-fluid simulations; as seen by the early formation of chains of magnetic islands inside the vortices in Fig. 8 (bottom left panel, ). This may be due to the intrinsic higher level of noise in particle simulations that seeds the tearing instability at a higher initial level, possibly through non-modal transient amplification, and makes it saturate much earlier.
V Discussions and conclusions
In this paper we have reported the first comparison of the magnetized Kelvin-Helmholtz instability, during its linear and nonlinear evolution. We have used several different plasma fluid and kinetic models: an isothermal/adiabatic ideal/resistive MHD, Hall-MHD, a two-fluid, a hybrid kinetic and a full kinetic model.
In the linear stage of the KHI, the MHD models do not care about the relative orientation of the vorticity with respect to the magnetic field direction, while the fluid simulations that include ion inertia (Hall-MHD, two-fluid) are insensitive to it; on the other hand,
kinetic simulations (Hybrid and Full PIC) clearly exhibit different growth rates depending on this relative orientation, showing that FLR effects dominate ion inertia, within the parameter range used in this study.
With the parameters used in this benchmark, the fastest growing mode m=2 builds up two rolled-up vortices in the nonlinear stage of the instability, that eventually merge into a single vortex or are disrupted before merging because of secondary magnetic reconnection.
At the end of the simulations, even if differences are present between the two sides of the layer, a common feature arises: although the KH vortices are large-scale MHD structures at the beginning of the nonlinear phase, their nonlinear evolution eventually leads to the formation of smaller length scales, both inside the large merged MHD vortex (left side) or in the mixing layer (right side). Even if small scale processes are strongly dependent on the choice of the model, interestingly, the global large-scale picture is well captured by all the considered models.
Note however that the degree of plasma mixing and the development of turbulence may strongly vary from one model to another at later times due to the nature of the small scale process at play.
The final stage of the instability is significantly different between the two shear layers, and all models capture the same final state, at least using a coarse graining view. Since the differences between the evolution of the two layers are also captured by the MHD model, we conjecture that such differences are not driven by kinetic or inertial effects, in particular not by the relative orientation of the vorticity with respect to the magnetic field. Previous works in the context of the MHD modeling of the KHI have shown that the nonlinear interactions differ dramatically, as influenced through the imposed phase differences between linear modes Baty, Keppens, and Comte (2003). In the setup described here, the initial phases of the first perturbed modes are fixed among the various models but are different for both shear layers. The origin of the two different final states in the two shear layer could rely upon the fact that the linear perturbations on both sides differ, as influenced by the initial phases, and can hence nonlinearly evolve differently.
The difference arising for a given model between the two different shear layers (see the differences between
Figs. 6 and 7 considering a same row) is much more important than
the difference arising between different models on the same shear layer (see the last columns of
Figs. 6 and 7). This indicates that the large
scale, fluid structure is much more influenced by the choice of the initial phases than by the small scale processes at play.
Moreover, the fluid models are shown to capture the large scale dynamics as well as the kinetic model, at least in the regime of parameters used in this study.
This indicates that the feedback of small (inertial and kinetic) scales to the large scale dynamics does not appear as a dominant process even in the nonlinear phase of the magnetized shear flow.
This study emphasizes the importance of the fluid behavior in the nonlinear evolution of magnetized shear flows, which determines the large scale structures and the saturation of the vortices. This result should strongly encourage the development of fluid codes to study the nonlinear dynamics of magnetized plasmas at large, fluid scales.
Note however that the complementary use of kinetic models is necessary to carry out the validation of the fluid approach at large scales.
We must also stress the important consequences of the different closures used in the different codes.
The MHD and two-fluid codes both use a standard adiabatic closure, using a polytropic law with polytropic index (on electrons also in the two-fluid model). The electrons are treated as an isothermal fluid in the hybrid PIC code.
On the contrary, no explicit evolution law is imposed on the ions (hybrid PIC, full PIC) and electrons (full PIC),
so that the behavior of effective ion (and electron for the full PIC code) compressibility remains a priori unknown
and may change in space and time. We recall that the compressibility plays a key role in the development of the KHI.
In particular, it modifies the growth rate of the instability. Finding a clever closure on the fluid codes,
in accordance with the compressibility found in the kinetic simulations, will enable to properly compare the growth rate,
in a first step, and the nonlinear dynamics of the KHI, in a second step.
In order to accurately describe the cross-scale, nonlinear evolution of collisionless plasmas, the coupling between different plasma models is showing recently an increasing interest from the plasma community Lapenta, Giovanni et al. (2013).
In this context, comparisons of numerical simulations from different models, as the one described in this paper,
represent a necessary step before coupling codes from different plasma models.
In order to properly couple different codes, it is a necessity to make sure that the different physical ingredients, introduced by the different models, describe the common features of interest at the large scales.
We have shown in this work that the large, fluid scales are little disturb by the small, kinetic scales, since fluid and kinetic simulations capture the same behavior at large scales. Such a result is likely to strongly encourage the development of multi-scale code coupling for colisionless plasmas.
We stress here that the numerical modeling of plasma dynamics is a fundamental problem of major interest in present plasma physics research.
Therefore, this study is of broad interest and is not limited to the KHI itself and related nonlinear dynamics.
We thus underline two important aspects of the fluid and kinetic modelling of the nonlinear evolution of magnetized plasma.
From a fluid point of view, the question of the fluid closure that plays an important role in the dynamics is shown not to be trivial in a magnetized plasma.
From a kinetic point of view, we have illustrated the necessity to find a correct initial setup that takes into account the kinetic effects at play, such as finite Larmor radius effects. These fundamental problems need to be addressed in future works.
Acknowledgment The research leading to these results has received funding from the European Commission’s Seventh Framework Programme (FP7/2007 2013) under the grant agreement SWIFF (project n 263340, www.swiff.eu). This work was supported by the Italian Super-computing Center CINECA under the ISCRA initiative. This work was supported by the HPC-EUROPA2 project (project number: 228398) with the support of the European Commission - Capacities Area - Research Infrastructures. UNIPI acknowledge the access to the HPC resources of CINECA made available within the Distributed European Computing Initiative by the PRACE-2IP, receiving funding from the European Community’s Seventh Framework Programme (FP7/2007-2013) under grant agreement n RI-283493 and project number 2012071282. Work at ASI and IAP, ASCR was also supported by projects RVO: 67985815 and RVO: 68378289. This work was supported by NASA award NNX11A1164G. UCPH acknowledges computer resources provided by the Danish Center for Scientific Computing (DCSC), now part of the Danish e-Infrastructure Collaboration (DeIC).
References
- Bruno and Carbone (2005) R. Bruno and V. Carbone, Living Reviews in Solar Physics 10 (2013), 2.
- Fairfield et al. (2000) D. H. Fairfield, A. Otto, T. Mukai, S. Kokubun, R. P. Lepping, J. T. Steinberg, A. J. Lazarus, and T. Yamamoto, Journal of Geophysical Research (Space Physics) 105, 21159 (2000).
- Hasegawa et al. (2004) H. Hasegawa, M. Fujimoto, T.-D. Phan, H. Rème, A. Balogh, M. W. Dunlop, C. Hashimoto, and R. TanDokoro, Nature 430, 755 (2004).
- Masters et al. (2010) A. Masters, N. Achilleos, M. G. Kivelson, N. Sergis, M. K. Dougherty, M. F. Thomsen, C. S. Arridge, S. M. Krimigis, H. J. McAndrews, S. J. Kanani, N. Krupp, and A. J. Coates, Journal of Geophysical Research (Space Physics) 115, 7225 (2010).
- Sundberg et al. (2011) T. Sundberg, S. A. Boardsen, J. A. Slavin, L. G. Blomberg, J. A. Cumnock, S. C. Solomon, B. J. Anderson, and H. Korth, Planetary and Space Sciences 59, 2051 (2011).
- Sundberg et al. (2012) T. Sundberg, S. A. Boardsen, J. A. Slavin, B. J. Anderson, H. Korth, T. H. Zurbuchen, J. M. Raines, and S. C. Solomon, Journal of Geophysical Research (Space Physics) 117, A04216 (2012).
- Faganello, Califano, and Pegoraro (2008a) M. Faganello, F. Califano, and F. Pegoraro, Physical Review Letters 101, 105001 (2008a).
- Faganello, Califano, and Pegoraro (2008b) M. Faganello, F. Califano, and F. Pegoraro, Physical Review Letters 101, 175003 (2008b).
- Faganello et al. (2010) M. Faganello, F. Pegoraro, F. Califano, and L. Marradi, Physics of Plasmas 17, 062102 (2010).
- Henri et al. (2012) P. Henri, F. Califano, M. Faganello, and F. Pegoraro, Physics of Plasmas 19, 072908 (2012).
- Keppens et al. (1999) R. Keppens, G. Tóth, R. H. J. Westermann, and J. P. Goedbloed, Journal of Plasma Physics 61, 1 (1999).
- Baty, Keppens, and Comte (2003) H. Baty, R. Keppens, and P. Comte, Physics of Plasmas 10, 4661 (2003).
- Matsumoto and Hoshino (2004) Y. Matsumoto and M. Hoshino, J. Geophys. Res. 31, 2807 (2004).
- Faganello, Califano, and Pegoraro (2008c) M. Faganello, F. Califano, and F. Pegoraro, Physical Review Letters 100, 015001 (2008c).
- Tenerani et al. (2011) A. Tenerani, M. Faganello, F. Califano, and F. Pegoraro, Plasma Physics and Controlled Fusion 53, 015003 (2011).
- Palermo et al. (2011a) F. Palermo, M. Faganello, F. Califano, F. Pegoraro, and O. Le Contel, Journal of Geophysical Research (Space Physics) 116, 4223 (2011a).
- Palermo et al. (2011b) F. Palermo, M. Faganello, F. Califano, and F. Pegoraro, Annales Geophysicae 29, 1169 (2011b).
- Cowee, Winske, and Gary (2009) M. M. Cowee, D. Winske, and S. P. Gary, Journal of Geophysical Research (Space Physics) 114, 10209 (2009).
- Cowee, Winske, and Gary (2010) M. M. Cowee, D. Winske, and S. P. Gary, Journal of Geophysical Research (Space Physics) 115, 6214 (2010).
- Wilber and Winglee (1995) M. Wilber and R. M. Winglee, Journal of Geophysical Research (Space Physics) 100, 1883 (1995).
- Nakamura, Hasegawa, and Shinohara (2010) T. K. M. Nakamura, H. Hasegawa, and I. Shinohara, Physics of Plasmas 17, 042119 (2010).
- Nakamura et al. (2011) T. K. M. Nakamura, H. Hasegawa, I. Shinohara, and M. Fujimoto, Journal of Geophysical Research (Space Physics) 116, 3227 (2011).
- Umeda et al. (2010) T. Umeda, J.-I. Miwa, Y. Matsumoto, T. K. M. Nakamura, K. Togano, K. Fukazawa, and I. Shinohara, Physics of Plasmas 17, 052311 (2010).
- Cai, Storey, and Neubert (1990) D. Cai, L. R. O. Storey, and T. Neubert, Physics of Fluids B 2, 75 (1990).
- (25) G. Baumann, K. Galsgaard, and Å. Nordlund, Solar Physics 284, 467 (2013).
- Williamson (1980) J. H. Williamson, Journal of Computational Physics 35, 48 (1980).
- Keppens et al. (2012) R. Keppens, Z. Meliani, A. J. van Marle, P. Delmont, A. Vlasis, and B. van der Holst, J. Comp. Phys. 231, 718 (2012).
- Čada and Torrilhon (2009) M. Čada and M. Torrilhon, Journal of Computational Physics 228, 4118 (2009).
- Faganello, Califano, and Pegoraro (2009) M. Faganello, F. Califano, and F. Pegoraro, New Journal of Physics 11, 063008 (2009).
- Valentini et al. (2007) F. Valentini, P. Trávníček, F. Califano, P. Hellinger, and A. Mangeney, J. Comp. Phys. 225, 753 (2007).
- Lele (1992) S. K. Lele, J. Comp. Phys. 103, 16 (1992).
- Matthews (1994) A. P. Matthews, J. Comput. Phys. 112, 102 (1994).
- Markidis, Lapenta, and Rizwan-uddin (2010) S. Markidis, G. Lapenta, and Rizwan-uddin, Mathematics and Computers and Simulation 80, 1509 (2010).
- Brackbill and Forslund (1982) J. Brackbill and D. Forslund, Journal of Computational Physics 46, 271 (1982).
- Cerri, S.S. et al. (2013) S. S. Cerri, P. Henri, F. Califano, D. Del Sarto, M. Faganello, and F Pegoraro, “Extended fluid models: pressure tensor effects and equilibria”, Phys. Plasmas (to be published).
- Nijboer et al. (1997) R. J. Nijboer, B. v. D. Holst, S. Poedts, and J. P. Goedbloed, Computer Physics Communications 106, 39 (1997).
- Passot, Sulem, and Hunana (2012) T. Passot, P. L. Sulem, and P. Hunana, Physics of Plasmas 19, 082113 (2012).
- Huba (1996) J. D. Huba, Geophysical Research Letters 23, 2907 (1996).
- Nagano (1978) H. Nagano, Journal of Plasma Physics 20, 149 (1978).
- Nagano (1979) H. Nagano, Planetary Space Science 27, 881 (1979).
- Fujimoto and Terasawa (1991) M. Fujimoto and T. Terasawa, Journal of Geophysical Research (Space Physics) 96, 15725 (1991).
- Birn et al. (2001) J. Birn, J. F. Drake, M. A. Shay, B. N. Rogers, R. E. Denton, M. Hesse, M. Kuznetsova, Z. W. Ma, A. Bhattacharjee, A. Otto, and P. L. Pritchett, Journal of Geophysical Research (Space Physics) 106, 3715 (2001).
- Lapenta, Giovanni et al. (2013) Lapenta, Giovanni, Pierrard, Viviane, Keppens, Rony, Markidis, Stefano, Poedts, Stefaan, Sebek, Ondrej, Trávnícek, Pavel M., Henri, Pierre, Califano, Francesco, Pegoraro, Francesco, Faganello, Matteo, Olshevsky, Vyacheslav, Restante, Anna Lisa, Nordlund, Åke, Frederiksen, Jacob Trier, Mackay, Duncan H., Parnell, Clare E., Bemporad, Alessandro, Susino, Roberto, and Borremans, Kris, J. Space Weather Space Clim. 3, A05 (2013).