Concerning the spatial discretization, a sixth-order differential operator has been considered for the convective fluxes, as in Bifrost. The conservative quantities are not co-located in space and are considered on a staggered grid, thus, a fifth-order interpolation is used to relocate variables as needed. Since a high-order finite difference scheme has been considered, it is convenient to numerically stabilize shocks and discontinuities with artificial numerical terms. In the literature, several type of artificial terms have been introduced, as described by [1]}, [2]} and references therein. By inheritance of the Bifrost  code, the structure of the artificial terms are the hyperdiffusive terms described by [3]}, [4]}. These terms are considered in order to capture shocks and discontinuities while treating correctly turbulence in a high Reynolds number regime. In Ebysus  the hyperdiffusive terms have been expanded for the multi-fluid case from those considered in Bifrost  since they depend on the velocity and speeds of sound for each species and the Alfven speed of all ions. The constant coefficients introduced in these terms have been chosen based on fits to a wide variety of shock and turbulence test problems.
