Fully Kinetic Large-scale Simulations of the Collisionless Magnetorotational Instability

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 (Pringle 1981). In particular, since the gravitational force decreases as the distance from the central object increases, the angular velocity of the plasma is therefore lower farther from the central object. This property of accretion disks induces the development of the so-called magnetorotational instability (MRI; Chandrasekhar et al. 1958; Balbus & Hawley 1991), through the action of which an initial seed (weak) magnetic field is exponentially amplified on a timescale 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) (Hawley et al. 1995; Balbus & Hawley 1998).

The current understanding of the MRI stems largely from MHD theory and simulation (Balbus & Hawley 1998). Goodman & Xu (1994) 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 (Goodman & Xu 1994; Latter et al. 2009; Pessah & Goodman 2009).

In radiatively inefficient accretion flow models for accretion onto compact objects, the accretion proceeds via a hot, low-density plasma with a proton temperature larger than the electron temperature (see Narayan et al. 1998; Quataert 2003 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., Quataert et al. 2002; Sharma et al. 2003, 2006, 2007; Krolik & Zweibel 2006) 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, so the magnetic moment $\mu ={{mv}}_{\perp }^{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, e.g., the mirror and the firehose (Rosenbluth & Garvin 1956; Parker 1958; Hasegawa 1969; Gary et al. 1993, 1997; Yoon et al. 1993; Hellinger & Matsumoto 2000; Pokhotelov et al. 2000, 2004; Gary & Karimabadi 2006; Kunz et al. 2014, 2015). 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 (Kunz et al. 2014, 2016; Mogavero & Schekochihin 2014; Melville et al. 2016).

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 (Gary & Karimabadi 2006). 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 and fully kinetic simulations—as does a 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 particle-in-cell (PIC) numerical studies, using a shearing co-rotating framework, Riquelme et al. showed that the MRI generates pressure anisotropies in both low (Riquelme et al. 2012) and high β (Riquelme et al. 2015) 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 (Hoshino 2013) and three-dimensional (Hoshino 2015) 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, Kunz et al. (2016) 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 and magnetic reconnection of the channel flows appears to be the key ingredient for triggering a turbulent regime.

This paper is organized as follows. In Section 2 we describe the shearing co-rotating framework implemented in our numerical code OSIRIS (Fonseca et al. 2002, 2013). The nonlinear evolution of the MRI is analyzed in Section 3 by describing the generation of the DKI in our larger simulations and how this influences the generation of the subsequent turbulent regime in the plasma. Section 4 focuses on the generation of pressure anisotropy and the consequent mirror instability. In Section 5 we investigate the turbulent saturation regime; the detailed analysis of the energy distribution and particle acceleration in this stage are reported in Section 6. Section 7 summarizes the results obtained in the paper and discusses future work.

Our numerical calculations employ the shearing co-rotating frame developed by Riquelme et al. (2012). For completeness, here we summarize 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 is 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₀ rotates with angular velocity ${\boldsymbol{\Omega }}={\rm{\Omega }}({x}_{0})\hat{z}\propto 1/{x}_{0}^{2}$. In our simulation, we study the evolution of a local portion of the disk, centered around the equilibrium position x₀, and require that the simulation domain in the radial direction, L_x, be small compared to the equilibrium position L_x ≪ x₀. In this approximation, the shearing velocity of the Keplerian disk, v₀ = ${\boldsymbol{\Omega }}$ ×  x₀, can be linearized to v₀ = $-3/2\alpha x\hat{y}$, where $\alpha =(d{\rm{\Omega }}/{dx}){| }_{{x}_{0}}$.

In order to include a differential shearing velocity v₀, the standard approach, used both in MHD simulations (see for example Hawley et al. 1995) and kinetic simulations (Hoshino 2013, 2015; Kunz et al. 2014), consists of the implementation of shearing periodic boundary conditions along the radial direction. An alternative method has been proposed by Riquelme et al. (2012), whereby one performs a Galilean transformation of Equations (1)–(5), implementing shearing coordinates.

To conduct our numerical simulations, we follow this last approach, modifying the PIC code OSIRIS (Fonseca et al. 2002, 2013) to include a local shearing co-rotating framework. In this particular frame, the shearing term v₀ 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 (Schiff 1939)

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{E}}=4\pi \rho +\displaystyle \frac{2{\boldsymbol{\alpha }}\cdot {\boldsymbol{B}}}{c}-\displaystyle \frac{{{\boldsymbol{v}}}_{0}}{c}\cdot {\rm{\nabla }}\times {\boldsymbol{B}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{B}}=0,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}=-c{\rm{\nabla }}\times {\boldsymbol{E}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial {\boldsymbol{E}}}{\partial t} & = & c{\rm{\nabla }}\times {\boldsymbol{B}}-4\pi {\boldsymbol{J}}+\displaystyle \frac{{{\boldsymbol{v}}}_{0}}{c}\times \displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}\\ & & -{\rm{\nabla }}\times \left({{\boldsymbol{v}}}_{0}\times \left({\boldsymbol{E}}-\displaystyle \frac{{{\boldsymbol{v}}}_{0}}{c}\times {\boldsymbol{B}}\right)\right),\end{eqnarray} \tag{ 4 }$$

where ${\boldsymbol{\alpha }}=\alpha \hat{z}$ 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 (1) and all terms proportional to v₀ in Equation (4). 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₀, provided that the equilibrium position x₀ is chosen to be large enough. With the non-relativistic approximation, the Maxwell's Equations (1)–(4) 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

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{p}}}{{dt}}=q\left({\boldsymbol{E}}+\displaystyle \frac{{\boldsymbol{v}}\times {\boldsymbol{B}}}{c}\right)-2{\boldsymbol{\alpha }}\times {\boldsymbol{p}},\end{eqnarray} \tag{ 5 }$$

where ${\boldsymbol{p}}$ and ${\boldsymbol{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 $| {\boldsymbol{u}}|$ is small compared to the shear velocity, $| {\boldsymbol{u}}| \ll | {{\boldsymbol{v}}}_{0}|$.

The equation of motion (5) is valid in any non-relativistic co-rotating frame. The non-relativistic restriction imposed above for the simplification of the Maxwell's Equations (1)–(4) and momentum (5) equation 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 also remains valid in the presence of relativistic particles, as long as the fluid motion of the plasma remains non-relativistic. For further discussion of the neglected relativistic effects in a generic co-rotating frame, we refer the reader to Riquelme et al. (2012).

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 Riquelme et al. (2012). The shearing, co-rotating Maxwell's Equations become

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{B}}=0,\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{E}}=4\pi \rho ,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}=-c{\rm{\nabla }}\times {\boldsymbol{E}}-\displaystyle \frac{3}{2}\alpha {B}_{x}\hat{y},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{E}}}{\partial t}=c{\rm{\nabla }}\times {\boldsymbol{B}}-4\pi {\boldsymbol{J}}-\displaystyle \frac{3}{2}\alpha {E}_{x}\hat{y},\end{eqnarray} \tag{ 9 }$$

and the equation of motion transforms to

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{p}}}{{dt}}=q\left({\boldsymbol{E}}+\displaystyle \frac{{\boldsymbol{v}}\times {\boldsymbol{B}}}{c}\right)-2{\boldsymbol{\alpha }}\times {\boldsymbol{p}}+\displaystyle \frac{3}{2}\alpha {p}_{x}\hat{y},\end{eqnarray} \tag{ 10 }$$

where we neglected the terms proportional to ∂/∂y in Equations (6)–(9), 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 (Krolik & Zweibel 2006), adapted for pair plasmas; this is reported in the Appendix.

2.1. Simulation Setup

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

Figure 1 shows the time evolution of the domain-averaged magnetic energy fluctuation $\langle \delta {B}^{2}\rangle$ 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.

Zoom In Zoom Out Reset image size

For simulation domains larger than 4λ₀, 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λ₀.

Figure 2 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 (Riquelme et al. 2012). The conclusions for the small box case are then the same as those of Riquelme et al. (2012); 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 the 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 (1)–(4) compared to the terms that we keep, and found that, on average, those ratios are less than 1% for the run D (16 × 16 box).

Zoom In Zoom Out Reset image size

We will focus the analysis on our largest simulation, corresponding to a simulation domain of 16λ₀ × 16λ₀, run D in Table 1.

Figure 3 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.1P₀. Interestingly, however, the value of β does not remain around unity and increases up to 10 during the saturation regime of the instability.⁴

Zoom In Zoom Out Reset image size

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

In Figure 4 we display snapshots of the average plasma density n = (n₊ + n₋)/2, the module of the in-plane magnetic field ${B}_{p}=\sqrt{{B}_{x}^{2}+{B}_{z}^{2}}$, the transverse magnetic field B_y, and the module of the in-plane current density ${j}_{p}=\sqrt{{j}_{x}^{2}+{j}_{z}^{2}}$ for run D at different times. T_orbit = 2P₀ 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₀. This regime is characterized by the formation of coherent structures—channel flows (Goodman & Xu 1994)—which, like the linear mode that they evolve from (see Figure 10 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.

Zoom In Zoom Out Reset image size

Channel flows are an exact solution of the nonlinear MHD equations (Goodman & Xu 1994); as such, their amplitudes will 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 Figure 4. From time T_orbit = 2P₀ to T_orbit = 2.8P₀, 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 Figure 5, right panel, which is a magnification of the red box identified in the plot of j_p at time T_orbit = 2.8P₀ in Figure 4. For contrast, the left panel of Figure 5 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.

Zoom In Zoom Out Reset image size

For times T_orbit > 2P₀, 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 (Pritchett et al. 1996; Daughton 1998, 1999a, 1999b) in the x–z plane of the simulation, along the x direction.⁵

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)λ₀ ∼ 20 c/ω_p—see Figure 5. 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 DKI, where k_DK δ ∼ 1 is expected (Daughton 1999b; Zenitani & Hoshino 2005).

Importantly, observe that λ_DK > λ₀, which partially accounts for the need to have large simulation domains. In particular, we note 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.

Figure 6 shows the spatial distribution of the pressure anisotropy at time T_orbit = 2P₀, 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 Figures 6(b) and (c). These results are consistent with previous numerical studies of the mirror instability (Riquelme et al. 2012; Hoshino 2013; Kunz et al. 2014).

Zoom In Zoom Out Reset image size

Figure 7 shows the distribution of the pressure anisotropy p_⊥/p_∥ as a function of β_∥ = 8πp_∥/B², where B is the total magnetic field, at different times. The solid line marks the stability threshold of the mirror instability, ${p}_{\perp }/{p}_{\parallel }\,-\,1/2\,\gt \sqrt{1/4+1/2{\beta }_{\parallel }}$ (Pokhotelov et al. 2000). The dashed line denotes the stability threshold of the firehose instability, ${p}_{\perp }/{p}_{\parallel }\,-\,1\lt 2/{\beta }_{\parallel }$ (Yoon et al. 1993; Hellinger & Matsumoto 2000). At time T_orbit = 2P₀, 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₀. At this time, there is thus a balance between pressure-anisotropy generation by the MRI, and its destruction by the mirror instability (Kunz et al. 2014; Riquelme et al. 2015; Melville et al. 2016).

Zoom In Zoom Out Reset image size

Note that in a realistic disk, one expects α/Ω₀ ≪ 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 Ω₀, 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.1P₀, 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.4P₀, this peak has grown, and we observe strong violation of the mirror stability boundary at β_∥ ∼ 3. As discussed in Section 3, 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 (Zenitani & Hoshino 2001; Sironi & Spitkovsky 2014; Dahlin et al. 2016) 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 is 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 (Lyubarsky 2005), where the factor 0.1 represents the characteristic relativistic reconnection rate (Liu et al. 2015; Cassak et al. 2017). With this assumption, we can estimate the velocity gain along y to be Δ(γv_y)/τ_A ∼ e/mE_y, where γ is the Lorentz factor. This yields a final proper velocity of u_y = γv_y ∼ 0.1v_Al/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 ≈ 10v_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 (Lyubarsky 2005). 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 Pokhotelov et al. (2000) τ_M ∼ (1 + (l/10d)²)τ_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 those of Kunz et al. (2016), 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.

Figure 8 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 $\langle {B}_{j}^{2}\rangle =\int d{{\rm{\Omega }}}_{k}{(k/2\pi )}^{2}| {B}_{j}{| }^{2}$, where $k=\sqrt{{k}_{x}^{2}+{k}_{z}^{2}}$ and Ω_k = tan⁻¹(k_x/k_z). During the linear regime (subplot a) at time T_orbit = 2.0P₀, 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 the Appendix). 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 (Pokhotelov et al. 2000). This coexistence of the MRI and the mirror instability is clearly visible in Figure 6.

Zoom In Zoom Out Reset image size

At T_orbit = 2.8P₀, 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 3 that β ≈ 1 at this time, so d ≈ ρ and thus also kd ≪ 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 (Kunz et al. 2014), 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 (Loureiro & Boldyrev 2017; Mallet et al. 2017; 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 4, second row.

At time T_orbit = 3.1P₀, 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.4P₀, the plasma is in a fully turbulent state, as shown in the last row of Figure 4. The last subplot in Figure 8 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 (Howes et al. 2008; Schekochihin et al. 2009; Chen et al. 2010; Boldyrev & Perez 2012; Passot & Sulem 2015), and also with reconnection-mediated turbulence  (Cerri & Califano 2017; Franci et al. 2017; Loureiro & Boldyrev 2017).

In Figure 9 we plot the evolution of the kinetic energy distribution at different times. At T_orbit = 2P₀, corresponding to the late linear evolution of the MRI, the energy distribution is still thermal. At times T_orbit = 2.8P₀ and T_orbit = 3P₀, 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 4 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 (Werner et al. 2016) 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²/(4πn(γm_ec² + 2.5T)), where T is the plasma temperature (Melzani et al. 2013). Using the values obtained in our simulations close to current sheets for these parameters ($l\sim 100-200c/{\omega }_{p}$ and σ ∼ 3−5), the predicted slope can vary between −2 and −2.5 (Werner et al. 2016), in good agreement with our simulations.

Zoom In Zoom Out Reset image size

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

In this work, we numerically analyzed the fully kinetic nonlinear evolution of a two-dimensional MRI in a collisionless, high -β pair plasma. Our main results are as follows. (i) The amplitude of the channel flows generated in the early nonlinear regime is limited by the onset of the 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.⁶ (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 also activated at larger scales 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, and is 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, the generation of pressure-anisotropy 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 electron and ion scales, which 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, such as for example, 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 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 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). N.F.L. was supported by US Department of Energy grant No. DE-FG02-91ER54109 and NSF CAREER award No. 1654168. G.I. acknowledges M. W. Kunz, M. Hoshino and C. Ruyer for valuable discussions.

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 a pure vertical initial magnetic field,⁷ and use the analytical results as a benchmark of our numerical code.

In the non-relativistic limit of the shearing co-rotating frame, the Faraday's and Ampére's equations are (see Equations (8)–(9))

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}=-c{\rm{\nabla }}\times {\boldsymbol{E}}-\displaystyle \frac{3}{2}\alpha {B}_{x}\hat{y}-\displaystyle \frac{3}{2}\alpha {ct}\displaystyle \frac{\partial {\boldsymbol{E}}}{\partial y}\times \hat{x},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{E}}}{\partial t}=c{\rm{\nabla }}\times {\boldsymbol{B}}-4\pi {\boldsymbol{J}}-\displaystyle \frac{3}{2}\alpha {E}_{x}\hat{y}-\displaystyle \frac{3}{2}\alpha {ct}\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial y}\times \hat{x}.\end{eqnarray} \tag{ 12 }$$

In the limit in which 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

$$\begin{eqnarray}{m}_{j}{n}_{j}\displaystyle \frac{\partial {{\boldsymbol{v}}}_{j}}{\partial t} & = & 2\alpha {m}_{j}{n}_{j}{v}_{z,j}\hat{y}-\displaystyle \frac{1}{2}\alpha {m}_{j}{n}_{j}{v}_{y,j}\hat{z}\\ & & +{q}_{j}{n}_{j}\left({\boldsymbol{E}}+\displaystyle \frac{{{\boldsymbol{v}}}_{{\boldsymbol{j}}}}{c}\times {\boldsymbol{B}}\right),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {n}_{j}}{\partial t}+{\rm{\nabla }}\cdot ({n}_{j}{{\boldsymbol{v}}}_{{\boldsymbol{j}}})=0,\end{eqnarray} \tag{ 14 }$$

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

$$\begin{eqnarray}&&{\boldsymbol{J}}=\displaystyle \sum _{j}{q}_{j}{n}_{j}{{\boldsymbol{v}}}_{j}.\end{eqnarray} \tag{ 15 }$$

We consider an external magnetic field ${\boldsymbol{B}}={B}_{0}\hat{z}$ and a quasi-neutral equilibrium (n_0,e ≈ n_0,i ≡ n₀). Linearizing Equations (11)–(15) and seeking solutions of the form exp(γt + ikz), we obtain the one-dimensional linear dispersion relation of MRI in the limit of a weak magnetic field (${v}_{{\rm{A}}}\to 0$), low rotational frequency (X ≫ 1), cold (β = 0), and pair plasma (R = 1),

$$\begin{eqnarray}&&4{\nu }^{4}+4(1+{\zeta }^{2}){\nu }^{2}+{\zeta }^{4}-6{\zeta }^{2}=0,\end{eqnarray} \tag{ 16 }$$

where we have adopted the following normalization:

$$\begin{eqnarray}&&\nu =\displaystyle \frac{\gamma }{\alpha },X=\displaystyle \frac{{\omega }_{{ci}}}{\alpha },R=\displaystyle \frac{{m}_{i}}{{m}_{e}},{v}_{{\rm{A}}}=\displaystyle \frac{{\omega }_{{ci}}}{{\omega }_{{pi}}},\zeta =\displaystyle \frac{{{kv}}_{{\rm{A}}}}{\alpha }.\end{eqnarray} \tag{ 17 }$$

Equation (16) is the same as that obtained in Krolik & Zweibel (2006) 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.007c/ω_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 ${{\boldsymbol{v}}}_{\mathrm{seed}}/c=1/20\,{v}_{{\rm{A}}}/c\,\sin (2\pi z/L)\hat{y}$, 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 10 shows very good agreement between Equation (16) 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.

Zoom In Zoom Out Reset image size