We present two-dimensional particle-in-cell (PIC) simulations of the fully kinetic collisionless magnetorotational instability (MRI) in weakly magnetized (high $\beta$) 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\rho<1$ and $k^{-3}$ for sub-Larmor scale ($k\rho>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 $\epsilon^{-2}$ in the high-energy region.