Fully kinetic large scale simulations of the collisionless Magnetorotational instability
Giannandrea Inchingolo Thomas Grismayer Nuno F. Loureiro Ricardo A. Fonseca Luis O. Silva
astro-ph.HE physics.plasm-ph





GoLP/Instituto de Plasmas e Fusão Nuclear, Instituto Superior Técnico, Universidade de Lisboa, 1049-001 Lisboa, Portugal
Plasma Science and Fusion Center, Massachusetts Institute of Technology, Cambridge, MA 02139, USA


GoLP/Instituto de Plasmas e Fusão Nuclear, Instituto Superior Técnico, Universidade de Lisboa, 1049-001 Lisboa, Portugal


Plasma Science and Fusion Center, Massachusetts Institute of Technology, Cambridge, MA 02139, USA


GoLP/Instituto de Plasmas e Fusão Nuclear, Instituto Superior Técnico, Universidade de Lisboa, 1049-001 Lisboa, Portugal
DCTI/ISCTE Instituto Universitário de Lisboa, 1649-026 Lisboa, Portugal


GoLP/Instituto de Plasmas e Fusão Nuclear, Instituto Superior Técnico, Universidade de Lisboa, 1049-001 Lisboa, Portugal


We present two-dimensional particle-in-cell (PIC) simulations of the fully kinetic collisionless magnetorotational instability (MRI) in weakly magnetized (high β) pair plasma.
The central result of this numerical analysis is the emergence of a self-induced turbulent regime in the saturation state of the collisionless MRI, which can only be captured for large enough simulation domains.
One of the underlying mechanisms for the development of this turbulent state is the drift-kink instability (DKI) of the current sheets resulting from the nonlinear evolution of the channel modes.
The onset of the DKI can only be observed for simulation domain sizes exceeding several linear MRI wavelengths.
The DKI, together with ensuing magnetic reconnection, activate the turbulent motion of the plasma in the late stage of the nonlinear evolution of the MRI.
At steady state, the magnetic energy has an MHD-like spectrum with a slope of k^-5/3 for kρ<1 and k^-3 for sub-Larmor scale (kρ>1).
We also examine the role of the collisionless MRI and associated magnetic reconnection in the  development of pressure anisotropy.
We study the stability of the system due to this pressure anisotropy, observing the development of mirror instability during the early-stage of the MRI.
We further discuss the importance of magnetic reconnection for particle acceleration during the turbulence regime.
In particular, consistent with reconnection studies, we show that at late times the kinetic energy presents a characteristic slope of ϵ^-2 in the high-energy region.

 




§ INTRODUCTION
 
Accretion disks are astrophysical structures in which a gas or a plasma rotates around a massive central object, such as a black hole or a neutron star, under the effect of the gravitational force <cit.>. In particular, since the gravitational force decreases as the distance from the central object increases, the angular velocity of the plasma is, therefore, lower far from the central object. 
This property of accretion disks induces the development of the so-called magnetorotational instability (MRI) <cit.>, through the action of which an initial seed (weak) magnetic field is exponentially amplified on a time scale comparable with the typical rotational period of the disk. 
Theoretical arguments and numerical simulations suggest that the saturation amplitude of the MRI is such that there is an approximate equipartition between the kinetic and magnetic energies (i.e., the plasma β, initially very large, saturates at values around 1) <cit.>.

The current understanding of the MRI stems largely from MHD theory and simulation <cit.>. <cit.> discussed the formation of large-scale coherent structures — channel flows — in the early nonlinear regime of the MRI.
These structures have been shown to be unstable to parasitic instabilities:  
Kelvin-Helmholtz, tearing, kink and pinch have all been suggested as possible modes that may play a crucial role in the disruption of the channel flows and subsequent activation of a turbulent stage <cit.>.

In radiatively inefficient accretion flow models for accretion onto compact objects, the accretion proceeds via a hot, low-density plasma with the proton temperature larger than the electron temperature (see <cit.> and <cit.> for reviews).
In order to maintain such a two-temperature flow, the typical collision rate must be much smaller than the accretion rate.
This suggests that the standard MHD approach for the description of the dynamics of such accretion disks may be insufficient, and a kinetic description is required instead. 
Indeed, several theoretical studies of collisionless MRI (e.g. <cit.>) have shown the development of pressure anisotropies during the evolution of the MRI when kinetic effects are taken into account. 
Fundamentally, this is due to the fact that in typical accretion disks, the growth rate of the MRI is much smaller than the ion cyclotron frequency, and so the magnetic moment μ = mv_⊥^2/2B, where v_⊥ is the component of the velocity perpendicular to the magnetic field B, ought to be conserved. 
The amplification of the magnetic field produced by the MRI therefore leads to an increase in v_⊥. 
The absence of significant collisions implies that v_⊥ and v_∥ will thus become different, originating a pressure anisotropy.
Several studies have shown that this anisotropy activates various kinetic instabilities, including, e.g., the mirror and the firehose <cit.>.
It is conjectured that these instabilities may significantly affect the nonlinear development of the MRI, and critically impact the transport of momentum and energy in accretion disks <cit.>.

The generation of a pressure anisotropy is not exclusive to the ions; indeed, it is expected that electrons will also develop a non-unity ratio of v_⊥ and v_∥, and thus trigger their own pressure anisotropy instabilities <cit.>. Their effect on MRI development and saturation, and the interplay between ion and electron scale instabilities, is not currently understood. 

Addressing these questions requires first-principles, fully kinetic simulations — as does the detailed understanding of energy partition and dissipation.
Unfortunately,  the typical range of scales and frequencies of the collisionless MRI is such that global kinetic simulations of accretion disks are impossible with present day computational resources. 
However, the simulation of a local portion of the disk (the local shearing-box approximation) is just about feasible, and might be sufficient to gain insight into how the collisionless MRI behaves. 
In their pioneering PIC numerical studies, using a shearing co-rotating framework, Riquelme et al. showed that the MRI generates pressure anisotropies in both low <cit.> and high β <cit.> regimes, and argued that to be the reason for their observation of the mirror instability. 
Nonetheless, those simulations did not reach the saturation stage of the MRI within the limits of validity of the shearing co-rotating model used. 
Using a co-rotating framework with shearing periodic boundary conditions, Hoshino performed two-dimensional <cit.> and three-dimensional <cit.> PIC simulations of electron-positron plasma for weakly magnetized, non-relativistic, collisionless MRI with  β≫ 100, confirming the generation of the mirror instability in the linear regime of the MRI and emphasizing the role of magnetic reconnection in the saturation of the MRI for the high β regime.
Due to the large computational cost of these studies, Hoshino's simulations were limited to the analysis of relatively small simulation domains.
Despite yielding saturation of the MRI, these simulations did not reach a turbulent stage.

Observing the saturation stage of the MRI with realistic ion-electron mass ratios and sufficiently large simulation domains is presently an insurmountable challenge for PIC simulations. 
A compromise can be found by using hybrid-kinetic codes, which treat the ions kinetically, but retain a fluid description of the electrons. 
Using this approach, <cit.> demonstrated that the saturation of the kinetic MRI proceeds via a steady-state turbulent regime.

This work presents an investigation that is complementary to these studies.
We perform ab-initio two-dimensional PIC simulations of collisionless MRI in a pair plasma. 
The PIC description of both species intrinsically includes all kinetic effects; but our study misses potential effects that critically depend on dimensionality and scale separation between the two species.
Such a compromise is imposed by the rather stringent computational limitations that characterize  this problem.
We investigate the dependence of our results on the size of simulation domain and observe, for the first time in fully kinetic studies, the turbulent saturation of the MRI, provided that the simulation domain is sufficiently large compared to the wavelength of the linearly most unstable MRI mode. 
When that is the case, we witness the onset of the drift-kink instability (DKI) in the nonlinear regime. The combined effect of this instability with magnetic reconnection of the channel flows, appears to be the key ingredients to trigger a turbulent regime.
 
This paper is organized as follows.
In Section <ref> we describe the shearing co-rotating framework implemented in our numerical code OSIRIS <cit.>.
The nonlinear evolution of the MRI is analyzed in Section <ref> describing the generation of the drift-kink instability in our larger simulations and how this influences the generation of subsequent turbulent regime in the plasma.
Section <ref> focuses on the generation of pressure anisotropy and the consequent mirror instability.
In section <ref> we investigate the turbulent saturation regime; the detailed analysis of the energy distribution and particle acceleration in this stage is reported in Section <ref>. 
Section <ref> summarizes the results obtained in the paper and discusses  future work.



§ SHEARING CO-ROTATING FRAME
 
Our numerical calculations employ the shearing co-rotating frame developed by <cit.>. For completeness, we summarize here the key aspects of this model; the reader is referred to that reference for details.

In order to study the evolution of the collisionless accretion disk, we investigate the two-dimensional (poloidal) x-z plane, where x is the radial direction of the accretion disk and z the vertical direction, parallel to the rotation axis. 
The y direction, perpendicular to the simulation plane, represents the transverse (toroidal) direction of the accretion disk. 
At equilibrium, the accretion disk follows Keplerian orbits around the central mass, where the plasma at each radial position x_0 rotates with angular velocity Ω=Ω(x_0)ẑ∝1/x_0^2. 
In our simulation, we study the evolution of a local portion of the disk, centered around the equilibrium position x_0, and require that the simulation domain in the radial direction, L_x, be small compared to the equilibrium position L_x≪ x_0. 
In this approximation, the shearing velocity of the Keplerian disk, v_0 = Ω×x_0, can be linearized to v_0 = -3/2α xŷ, where α=(dΩ/dx)|_x_0.

In order to include a differential shearing velocity v_0, the standard approach, used both in MHD simulations (see for example <cit.>) and kinetic simulations <cit.>, consists in the implementation of shearing periodic boundary conditions along the radial direction. 
An alternative method has been proposed by <cit.>, whereby one performs a Galilean transformation of equations (<ref>-<ref>) implementing shearing coordinates.

To conduct our numerical simulations, we follow this last approach, modifying the PIC code OSIRIS <cit.> to include a local shearing co-rotating framework.
In this particular frame, the shearing term v_0 appears explicitly in the equations.
The implementation of the shearing co-rotating frame requires a series of approximations that we now discuss.
The co-rotating reference frame is non-inertial, so Maxwell's equations become <cit.>

    ∇·E    =     4πρ+2α·B/c-v_0/c·∇×B, 
    ∇·B    =     0, 
    ∂B/∂ t    =     -c∇×E, 
    ∂E/∂ t    =     c∇×B - 4 πJ + v_0/c×∂B/∂ t
           -∇×(v_0×(E-v_0/c×B)),

where α =αẑ is the angular frequency of the accretion disk.
In order to simplify these equations, we limit our analysis to the non-relativistic case in which we can neglect the last two terms in Equation (<ref>) and all terms proportional to v_0 in Equation (<ref>).
Strictly speaking, the non-relativistic approximation would also require us to neglect the displacement current. 
However, the details of the PIC numerical algorithm require us to keep it to update the electric field. There is no inconsistency between keeping the displacement current and neglecting the terms proportional to v_0 provided that the equilibrium position x_0 is chosen to be large enough.
With the non-relativistic approximation, the Maxwell's equations (<ref>-<ref>) simplify to the usual ones.

In the co-rotating frame, the motion of the plasma is affected by the Coriolis force. In the case of a Keplerian disk, this is given by the well-known expression

    d p/dt=q(E+v×B/c)-2α×p,

where p and v are the particle momentum and velocity and q is its charge. 
The expression for the Coriolis force is valid in the cold limit, where the fluid velocity |u| is small compared to the shear velocity, |u|≪ |v_0|. 

The equation of motion (<ref>) is valid in any non-relativistic co-rotating frame.
The non-relativistic restriction imposed above for the simplification of the Maxwell's (<ref>-<ref>) and momentum (<ref>) equations refers to the main bulk velocity of the plasma.
There is no restriction on the motion of single particles, which can in principle be relativistic.
The set of equations for the non-relativistic case remains valid also in the presence of relativistic particles, as long as the fluid motion of the plasma remains non-relativistic.
For further discussion on the neglected relativistic effects in a generic co-rotating frame, we refer the reader to <cit.>.

To move from the co-rotating frame just described – in the non-relativistic limit – to the final shearing co-rotating frame, we apply a Galilean transformation, following the approach described in the Appendix of <cit.>.
The shearing, co-rotating Maxwell's equations become

    ∇·B    =     0, 
    	
    ∇·E    =     4πρ, 
    
        ∂B/∂ t    =     -c∇×E - 3/2α B_xŷ, 
    
        ∂E/∂ t    =     c∇×B-4πJ -3/2α E_xŷ,

and the equation of motion transforms to

    d p/dt=q(E+v×B/c)-2α×p+3/2α p_xŷ,

where we neglected the terms proportional to ∂/∂ y in Eq. (<ref>-<ref>) since we restrict our study to two-dimensions (the poloidal (x-z) plane).
As we present our results below, we will attempt to discuss how the extension to a three-dimensional setup might affect them, or not.

To verify the validity of the results obtained in the two-dimensional shearing co-rotating frame, we performed benchmarks against the linear theory of collisionless MRI <cit.>, adapted for pair plasmas; this is reported in Appendix <ref>.



 §.§ Simulation setup

We start with a non-relativistic, isotropic, weakly magnetized pair plasma (e_+-e_-) with β=8π(p_+ + p_-)/B^2_0=100, where the pressure of each species is related to their respective thermal velocity v_th,±=(3k_BT_±/m)^1/2 by p_± = (1/2)m nv^2_th,±.
As discussed in the Introduction, the choice of a pair plasma enables us to simulate several orbital periods (2π/α) of the accretion disk and much larger simulation domains  that would not be possible with realistic mass ratios.

The external magnetic field is set to be vertical to the accretion disk, i.e., B_0= B_0ẑ.
Its initial value is set using the corresponding Alfvén speed v_A,0=B_0/√(4π m n) and fixed to v_A,0/c=1.43×10^-2. 
The orbital frequency α is expressed in terms of the initial cyclotron frequency Ω_0=eB_0/mc and fixed to α/Ω_0 = 1/11. 
The simulation domains L_x and L_z are normalized by λ_0=2π v_A/α (approximately the wavelength of the fastest growing MRI mode <cit.> and vary from 2 λ_0 to 16 λ_0 (λ_0  is related to the plasma skin depth d=c/ω_p by λ_0 = 2π (v_A/c)(ω_p/α) d).
Time is normalized to the orbital period P_0 = 2π/α. 
Other numerical parameters are summarized in Table <ref>. The spacial resolution is chosen such as to simultaneously resolve both the skin depth d and the Larmor radius ρ (with ρ=√(β)d and β>1) during the evolution of the simulation. 





§ EFFECT OF SIMULATION DOMAIN SIZE
 


In this section, we analyze the time evolution of the collisionless MRI, with particular focus on two key ingredients of its nonlinear evolution: magnetic reconnection and pressure anisotropy generation. 

Fig.<ref> shows the time evolution of the domain-averaged magnetic energy fluctuation <δ B^2> for different simulation domains. 
We highlight four different phases that represent, respectively, the transition between the linear and the nonlinear regime (labeled `1'), and three meaningful stages of the nonlinear evolution that will be further analyzed below.

For simulation domains larger than 4λ_0, we observe that our numerical simulations have converged. The noticeable differences in the nonlinear evolution between the small and the large simulation domains indicate that this is a key parameter in determining the dynamics, as will be further documented below.
The amplitude of the magnetic energy in the saturation regime of the instability decreases by almost one order of magnitude from small to large simulation domains until reaching convergence when the size of the simulation domain L≥8λ_0.


Figure <ref> shows the time evolution of the domain-averaged Alfvén velocity v_A for different simulation domains.
We observe that the evolution of v_A depends on the simulation domain size: for small domains (run A), the Alfvén speed reaches values of v_A∼3c, in agreement with <cit.>. The conclusions for the small box case are then the same as those of Riquelme et al.: for small domains, we do not observe a saturation of the MRI within the limit of validity of the shearing co-rotating framework. 
However, when the simulation domain size is increased, we observe that the exponential growth of the MRI saturates at values of the Alfvén velocity below the speed of light, within the validity range of our model; this is in agreement with our general observation that the size of simulation domain is critical to correctly capture the nonlinear dynamics of the collisionless MRI.
In addition, we checked a posteriori the magnitude of the neglected terms in the Maxwell’s equations <ref>-<ref> compared to the terms that we keep, and found that, on average, those ratios are less than 1% for the run D (16x16 box).

We will focus the analysis on our largest simulation, corresponding to a simulation domain of 16λ_0×16 λ_0,  run D in Table <ref>.


Figure <ref> shows the time evolution of the domain-averaged magnetic and kinetic energy and the β parameter for this run.
During the amplification of the magnetic field produced by the MRI, the β parameter decreases from its starting value of 100, reaching the equipartition value at time T_orbit = 3.1 P_0. Interestingly, however, the value of β does not remain around unity and increases up to 10 during the saturation regime of the instability.
[In figure <ref> we observe that the beta parameter is practically constant throughout the late stages of evolution of the instability. However, both the magnetic and kinetic energy slowly rise during this period. In our setup, there is a continuous source of kinetic energy injection (the shearing term in our equations). The constancy of beta indicates that the system has reached an equilibrium between the kinetic and magnetic energies. However, in a true steady-state, the energy injected should match the energy dissipated, and thus both the kinetic and the magnetic energy should be constant (on average); this is not what we observe. It is possible that to attain a real steady state we would have to run the simulations for much longer, this would require significant computing resources that we currently do no have. Another possibility is that this secular growth is caused by insufficient energy dissipation in our code. Another possibility is that this is a manifestations of residual numerical effects such as numerical heating intrinsic to the PIC algorithm.]

As we will describe later, the magnetic reconnection that is activated along the phase 3 is the mechanism responsible for the growth of β, progressively converting the magnetic field energy into kinetic energy. 
This behavior is observed also in our smaller simulation domains (not shown here), and was also manifest in the previous work of <cit.>. 

 
In Fig. <ref> we display snapshots of the average plasma density n=(n_++n_-)/2, the module of the in-plane magnetic field B_p = √(B_x^2+B_z^2), the transverse magnetic field B_y and the module of the in-plane current density j_p=√(j_x^2+j_z^2) for run D at different times.
T_orbit = 2 P_0 represents the transition between the linear regime and the nonlinear regime of the MRI in the sense that the growing perturbation of the magnetic field δ B starts to be on the order of the initial external magnetic field B_0. 
This regime is characterized by the formation of coherent structures — channel flows  <cit.> — which, like the linear mode that they evolve from (see Fig.<ref> in the Appendix), are on MHD scales. 
It is worth commenting in passing that this shows that channel flows are a robust MHD solution, still observed in the fully kinetic, collisionless regime that we explore here.

Channel flows are an exact solution of the nonlinear MHD equations <cit.>; as such, their amplitude would grow unbounded unless they are disrupted by parasitic instabilities.
We will now focus on how this occurs in run D.
The channel flows confine the plasma density between regions of positive and negative transverse magnetic field B_y, with the formation of current sheets in the plane of the simulation, as shown by the in-plane current density j_p in Fig. <ref>. 
From time T_orbit = 2 P_0 to T_orbit = 2.8 P_0, we observe that the thickness of the current sheets decreases with the growth of the transverse component B_y of the magnetic field. 
These current sheets subsequently develop a radial modulation, as illustrated in Fig. <ref>, right panel, which is a magnification of the red box identified in the plot of j_p at time T_orbit = 2.8 P_0 in Fig. <ref>.
For contrast, the left panel of Fig. <ref> shows a snapshot of j_p from our smallest simulation (run A) at the same time.
This comparison shows that the amplitude of this modulation is significantly more pronounced in the large simulation domain. 
 

For times T_orbit > 2 P_0, the plasma possesses well-defined current sheets along the radial direction that carry both J_x and J_y currents, with comparable magnitudes. 
These current sheets are surrounded by shearing magnetic fields B_x and B_y,  with B_x ∼ B_y. 
In principle, this configuration allows for the development of both tearing and drift-kink modes <cit.> in the x-z plane of the simulation, along the x direction.[The drift kink investigations of <cit.> consider a configuration where the only component of the magnetic field would correspond to our B_x. 
The channel flows whose stability we are discussing are threaded by both B_x and B_y, with B_x/B_y∼ 1, and thus it is not immediately obvious that those results on the DK instability are still valid here. 
However, note that both B_x and B_y are modulated in the z-direction, and we will find that the DK instability is located at values of z where B_y≈ 0, legitimizing our comparison with the aforementioned theories. ] 

The characteristic thickness of the current sheets before the onset of the modulation can be measured from the simulation to be δ∼ (0.2-0.3)λ_0 ∼ 20 c/ω_p —
see Fig. <ref>. 
From the simulation, we can also measure the average wavenumber of the modulation to be approximately k=2π/λ∼ 0.045 ω_p/c. 
These measurements are consistent with previous numerical studies of the drift-kink instability, where k_DKδ∼ 1 is expected <cit.>. 

Importantly, observe that λ_DK>λ_0, which partially accounts for the need to have large simulation domains.
In particular, we notice that the amplification of the DK modulation is larger for larger domains.
For small domains, the growth of the magnetic field proceeds unhindered until the current sheets become unstable to reconnection; in such cases, modulation of the current sheet due to the DKI is small or non-existent. 
This has important consequences for the subsequent nonlinear dynamics. 
In small domains, where DKI is mostly absent, the motion of the magnetic islands formed once the channel flows break is mostly confined in the radial (x) direction of the simulation. 
When the size of the simulation domain is increased, the different current sheet structures and the larger amplitude of the DK modes activate a non-uniform motion of the current sheets along the vertical (z), as well as radial, directions. Additionally, there is a much larger variety of island sizes produced. The combination of these different effects results in a transition to fully turbulent dynamics that is absent in smaller simulation domains.



§ PRESSURE ANISOTROPY DRIVEN INSTABILITY


Fig. <ref> shows the spatial distribution of the pressure anisotropy at time T_orbit = 2 P_0, represented by the parameter Δ = p_⊥/p_∥-1. 
The pressure anisotropy grows with the growth of the magnetic field induced by the MRI. 
In particular, the regions of maximum anisotropy are where the magnetic field forms filaments oblique to the direction of the external magnetic field, as shown in Fig. <ref> b) and c). 
These results are consistent with previous numerical studies of the mirror instability <cit.>.

Figure <ref> shows the distribution of the pressure anisotropy p_⊥/p_∥ as a function of β_∥ = 8π p_∥/B^2, where B is the total magnetic field, at different times. 
The solid line marks the stability threshold of the mirror instability, p_⊥/p_∥ -1/2 > √(1/4+1/2β_∥) <cit.>.
The dashed line denotes the stability threshold of the firehose instability,  p_⊥/p_∥ -1 < 2/β_∥ 2there was a mistake in the threshold, missing a factor 2<cit.>.
At time T_orbit = 2 P_0, the distribution of the pressure anisotropy is mostly above the threshold of the mirror instability. 
The combined action of further magnetic field amplification due to the MRI and of the mirror instability moves the pressure anisotropy to regions of lower β_∥ and within the stability margins, as obtained at T_orbit=2.8P_0. 
At this time, there is thus a balance between pressure anisotropy generation by the MRI, and its destruction by the mirror instability <cit.>.

Note that in a realistic disc, one expects α/Ω_0 ≪ 1/β ; the growth rate of the mirror instability would therefore always be much larger than that of the MRI, implying that the signature of the mirror instability should appear earlier in time than what we obtain in our simulations.
Due to numerical constraints, this condition is not verified with the initial parameters of our simulations. During the μ-conserving phase of the MRI, however, the growth of the magnetic field simultaneously reduces the β parameter and increases Ω_0, such that the above condition becomes verified. This delay of the effects of the mirror instability explains why its saturation only occurs at ∼2.8 orbits.

When the MRI starts to saturate at T_orbit = 3.1 P_0, the distribution of the pressure anisotropy is completely within the stability bounds. 
However, an interesting feature emerges at this stage:  a secondary peak of the distribution arises at p_⊥/p_∥∼ 1.
At time T_orbit = 5.4 P_0, this peak has grown, and we observe strong violation of the mirror stability boundary at β_∥∼3. 
As discussed in Sec. <ref>, at this time the motion of the plasma is turbulent and magnetic reconnection plays a crucial role in the dynamics and coalescence of plasma islands. 
We think reconnection is the origin of this violation of the mirror stability threshold, as we now explain.

In two-dimensional geometry, previous numerical studies of reconnection <cit.> show that the particles are accelerated mostly along the direction perpendicular to the reconnection plane (y in our configuration).
The typical duration of a reconnection event in our simulations is the Alfvén time τ_A=l/v_A, where l is the length of a given current sheet and v_A the local Alfvén velocity.
Along the direction perpendicular to the plane of reconnection, the dominant contribution to the acceleration of a particle is due to the reconnection  electric field. 
The amplitude of this electric field is estimated to be (in normalized units) E_y∼0.1v_A/cB <cit.>, where the factor 0.1 represents the characteristic relativistic reconnection rate <cit.>. 
With this assumption, we can estimate the velocity gain along  y to be Δ(γ v_y)/τ_A∼ e/m E_y, where γ is the Lorentz factor. 
This yields a final proper velocity of  u_y=γ v_y∼ 0.1v_A l/d, where d=c/ω_p is the skin depth.
In our simulations, the typical length of the current sheet is l≳100 c/ω_p and thus, on average, we expect v_y ≈ 10 v_A.
In the reconnection plane (i.e., the plane of the simulation), instead, reconnection-accelerated particles typically move at the proper Alfvén velocity v_x∼ v_A <cit.>. 
 This suggests a mechanism for the generation of velocity anisotropy, with v_y>v_x. 

With the above estimates, one can predict the growth time for the mirror instability resulting from this reconnection-generated pressure anisotropy to be <cit.> τ_M∼(1+(l/10d)^2)τ_A∼ 100τ_A. 
This is much longer than the typical reconnection event, explaining, we believe, why the mirror stability boundaries are violated at this stage of our simulations.
 
In this regard, our results are different from <cit.>, where, instead, the mirror instability threshold remains a solid boundary constraining the nonlinear dynamics. 
It is conceivable that the differences between our results and theirs stem from additional constraints imposed by the two-dimensional geometry that we use; but, alternatively, it is also possible that  we have uncovered an effect that critically depends on a kinetic treatment of electrons, which Kunz et al. (2016) do not do.
The extension of our work to fully three-dimensional geometries requires extraordinary computational resources and must thus be left for future work.



§ MRI TURBULENCE


Figure <ref> shows the magnetic energy spectrum for both the in-plane B_p and transverse B_y components at different times.
The energy spectrum is defined as <B_j^2> = ∫ dΩ_k(k/2π)^2|B_j|^2, where k=√(k_x^2+k_z^2) and Ω_k = tan^-1(k_x/k_z).
During the linear regime (subplot a) at time T_orbit = 2.0 P_0, the energy spectrum shows a peak at k_MRI∼0.04 ω_p/c, that corresponds to the maximum wavelength of the collisionless MRI in our system (see appendix <ref>). 
In the in-plane B_p energy spectrum, we can also observe a secondary peak at k_Mirror∼0.15 ω_p/c which is consistent with the maximum wavelength of the mirror instability <cit.>. This coexistence of the MRI and the mirror instability is clearly visible in Figure <ref>.

At T_orbit = 2.8 P_0, the energy spectrum shows a well-defined power law distribution for high k, with a -3 slope for both the in-plane and transverse component of the magnetic field. 
This slope occurs at kρ≪1 (recall from Figure <ref> that β≈ 1 at this time, so d≈ρ and thus also k d≪ 1).
It is not obvious why we observe this power-law behavior at this stage, since it occurs before the transition to fully-developed MRI turbulence. 
It could conceivably be the result of mirror instability-driven turbulence, except that that is predicted to yield a -5/3 slope at the fluid scales <cit.>, different from what we observe.
A tentative explanation is that this power-law behavior is due to the effect of magnetic reconnection at those scales, as suggested by recent analytical predictions <cit.> (with the caveat that these predictions were made for ion-electron plasmas, not pair plasmas). Visual evidence for magnetic reconnection occurring already at this time is discernible in the contour plots of Figure <ref>, second row. 

At time T_orbit = 3.1 P_0, the role of reconnection in the nonlinear dynamic of the system is more pronounced, with the disruption of the MRI channel flows and activation of large scale turbulence.
The transverse component of the magnetic energy spectrum maintains the slope of -3 at fluid scales.
The in-plane component of the magnetic field, instead, shows a transition from a slope consistent with the familiar -5/3 at large scales, to -3.

At time T_orbit = 5.4 P_0, the plasma is in a fully turbulent state, as shown in the last row of Fig. <ref>. The last subplot in Fig <ref> shows that the spectrum has a characteristic slope of k^-5/3 for scales bigger than the Larmor radius (kρ <1), and retains the -3 slope at sub-Larmor scales. Such a slope in the kinetic range is consistent with expectations from kinetic Alfvén (or perhaps whistler) wave turbulence <cit.>, but also with reconnection-mediated turbulence  <cit.>.



§ PARTICLE ACCELERATION


In Fig. <ref> we plot the evolution of the kinetic energy distribution at different times. 
At T_orbit = 2 P_0, corresponding to the late linear evolution of the MRI, the energy distribution is still thermal. 
At times T_orbit = 2.8 P_0 and T_orbit = 3 P_0, corresponding to the activation of reconnection, we observe the development of a hot tail in the energy distribution, consistent with the claim made in Section <ref> that reconnection is efficient at the generation of high-energy particles. 
At late times, during the steady state turbulence regime, we see that the kinetic energy distribution has evolved to exhibit a clear power law slope of -2 at high energies.
Recent numerical studies of particle acceleration via relativistic magnetic reconnection <cit.> show that the slope of the particle acceleration depends of the ratio between the the typical length of the current sheets l and the magnetization parameter σ. For a hot plasma, as is the case here, the magnetization parameter can be defined as σ = B^2/(4π n(γ m_ec^2+2.5T)), where T is the plasma temperature.<cit.>
Using for these parameters the values obtained in our simulations close to current sheets (l∼100-200 c/ω_p and σ∼3-5), the predicted slope can vary between -2 and -2.5 <cit.>, in good agreement with our simulations.

The role of reconnection in generating non-thermal tails in the nonlinear development of the MRI has been previously pointed out by <cit.> who, however, obtained a shallower slope of -1, possibly due to the smaller simulation domain employed there.



§ SUMMARY

In this work, we numerically analyzed the fully kinetic nonlinear evolution of two-dimensional MRI in a collisionless, high β pair plasma. 
Our main results are the following. 
(i) The amplitude of the channel flows generated in the early nonlinear regime is limited by the onset of the drift-kink instability (DKI), and subsequent magnetic reconnection. 
Both play a critical role in the transition to a regime of fully-developed turbulence. 
Importantly, we observe that simulations with smaller domain sizes yield insignificant DKI; as a result, such simulations fail to develop substantive turbulent dynamics.
One of our noteworthy conclusions, therefore, is that sufficiently large simulation domains (compared to the wavelength of the most unstable MRI mode) are required to reach saturation.[This conclusion does not explicitly contradict any previous numerical study that we are familiar with. The two-dimensional kinetic MHD simulations reported in <cit.> fail to saturate, but they are performed in simulation domains of the order of  ∼λ_0. For such a domain, our simulations also do not saturate.]
(ii) Reconnection leads to significant velocity anisotropy. It is more efficient at generating such anisotropy than the mirror instability is at destroying it. 
As a result, the nonlinear turbulent plasma state significantly violates the mirror stability boundary.  
(iii) During the initial phase of the nonlinear MRI, the magnetic energy spectrum presents a characteristic slope of -3 for scales larger that the Larmor radius.
We interpreted this result as the generation of a turbulent regime at these scales driven by magnetic reconnection.
Subsequently, the turbulence is activated also at larger scale and the magnetic energy spectrum presents a -5/3 slope for scales larger than the Larmor radius and a slope -3 for sub-Larmor scales.
(iv) A energetic particle spectrum is obtained, well described at high energies by a power-law with slope -2.
 
Two obvious limitations of our work are the reduced mass ratio employed (we consider a pair plasma) and the dimensionality (2D instead of 3D). Going forward, we aim to address both of these issues to assess the extent to which they affect the conclusions drawn here.
In particular, pressure-anisotropy generation by magnetic reconnection may be strongly impacted by the dimensionality of the setup, given that in 3D we expect the current sheets to be oriented obliquely to the x-z plane that we simulate here. 
In addition, the introduction of a significant mass ratio between the two plasma species may lead to interesting multiscale effects, such as the generation of pressure anisotropy-driven instabilities both at electrons and ions scales, that are absent from our simulations and cannot be captured via hybrid simulations.

At a qualitative level, the MRI dynamics evidenced by our simulations resembles that observed in MHD studies: the linear growth period is followed by a stage of channel flow formation; these channel flows disrupt due to parasitic instabilities and lead to a fully turbulent saturated state. Quantitatively, however, there are differences: for example, in the way that kinetic instabilities (i.e. mirror, drift kink) regulate certain stages of the evolution; the fact that our dominant parasitic mode seems to be magnetic reconnection instead of the usual Kelvin-Helmholtz of MHD, and how this is intimately linked to efficient particle acceleration; and, of course the details of the energy spectrum itself. Whether these differences matter in terms of transport and other macroscopic properties of the system remains to be understood.

This work was supported in part by the Fundação para a Ciência e a Tecnologia (FCT) under the grant PD/BD/105855/2014 and by the European Research Council under the Advanced Grant In-Pairs n.695088.
We would like to acknowledge the assistance of high performance computing resources (Tier-0) provided by PRACE on Marenostrum4 based in Spain.
Simulations were performed at the IST cluster (Lisbon, Portugal), and the Marenostrum4 supercomputer (Barcelona, Spain).
NFL was supported by NSF CAREER award no. 1654168.
GI acknowledges M. W. Kunz, M. Hoshino and C. Ruyer for valuable discussions.





§ LINEAR REGIME OF MRI

In this appendix, we derive a one-dimensional two-fluid model for the linear regime of the MRI in the shearing co-rotating framework for pure vertical initial magnetic field[This can be generalized to include a finite azimuthal field B_y <cit.>; in that case, the dispersion relation becomes sensitive to kinetic effects and would provide a complementary test of our algorithm. Here we adopt this setup because it is the one used for our simulations.], and use the analytical results as a benchmark of our numerical code. 

In the non-relativistic limit of the shearing co-rotating frame,  Faraday's and Ampére's equations are (see Eq. <ref> - <ref>):

    ∂B/∂ t    =     -c∇×E - 3/2α B_xŷ -3/2α ct∂E/∂ y×x̂, 
    ∂E/∂ t    =     c∇×B-4πJ -3/2α E_xŷ  -3/2α ct    ∂B/∂ y×x̂.


In the limit that all fluctuations have wavelengths much larger than the ion Larmor radius and frequencies much smaller than the ion cyclotron frequency, the momentum and the continuity equations become:

    m_jn_j∂v_j/∂ t = 2α m_jn_jv   _z,jŷ - 1/2α m_jn_jv_y,jẑ + q_jn_j(E + v_j/c×B), 
    ∂ n_j/∂ t + ∇·(n_jv_j) =0,

where the suffix j indicates both the particles species (electrons and ions). 
The current is computed from Ampére's law as:

    J = ∑_j q_jn_jv_j.


We consider an external magnetic field B = B_0ẑ and a quasi-neutral equilibrium (n_0,e≈ n_0,i≡ n_0). 
Linearizing Eq. (<ref> - <ref>) and seeking solutions of the form exp(γ t + ikz), 
 we obtain the one-dimensional linear dispersion relation of MRI in the limit of weak magnetic field (v_A→0), low rotational frequency (X ≫ 1), cold (β = 0) and pair plasma (R = 1):

    4ν^4+4(1+ζ^2)ν^2+ζ^4-6ζ^2=0,

where we have adopted the following normalization:

    ν=γ/α, X=ω_ci/α, R=m_i/m_e,  v_A=ω_ci/ω_pi,ζ=kv_A/α.

Eq. <ref> is the same as obtained in <cit.> in the limit of pair plasma.

To verify the validity of our numerical shearing co-rotating framework, we consider a pair plasma with β = 0.05. The external magnetic field is defined by the Alfvén velocity v_A = 0.05c and the angular frequency α by two different values of X=11,  33. 
The grid resolution of the numerical simulation is set to Δ x = 0.007 c/ω_pe with 1000 particles per cell. 
To force the excitation of a specific MRI wavelength in our simulation, we seed the plasma with a velocity profile v_seed/c = 1/20  v_A/c sin(2π z/L)ŷ, where L is the size of the domain, scanned over the values L=0.37, 0.4, 0.44, 0.5, 0.625, 0.8, 1.0, 1.25, 1.6, 2.0. 

Figure <ref> shows very good agreement between Eq. (<ref>) and the numerical results obtained with the one-dimensional version of the shearing co-rotating version of the PIC code OSIRIS that we developed for this work.

aasjournal


