  Supersonic magnetohydrodynamic (MHD) turbulence is a ubiquitous state for
many astrophysical plasmas. However, even the basic statistics for this type of
turbulence remains uncertain. We present results from supersonic MHD turbulence
simulations at unparalleled resolutions, with plasma Reynolds numbers of over a
million. In the kinetic energy spectrum we find a break between the scales that
are dominated by kinetic energy, with spectral index $-2$, and those that
become strongly magnetized, with spectral index $-3/2$. By analyzing the
Helmholtz decomposed kinetic energy spectrum, we find that the compressible
modes are not passively mixed through the cascade of the incompressible modes.
At high magnetic Reynolds number, above $10^5$, we find a power law in the
magnetic energy spectrum with spectral index $-9/5$. On the strongly
magnetized, subsonic scales the plasma tends to self-organize into locally
relaxed regions, where there is strong alignment between the current density,
magnetic field, velocity field and vorticity field, depleting both the
nonlinearities and magnetic terms in the MHD equations, which we attribute to
plasma relaxation on scales where the magnetic fluctuations evolve on shorter
timescales than the velocity fluctuations. This process constrains the cascade
to inhomogenous, volume-poor, fractal surfaces between relaxed regions, which
has significant repercussions for understanding the nature of magnetized
turbulence in astrophysical plasmas and the saturation of the fluctuation
dynamo.
