We describe a novel Godunov-type numerical method for solving the equations of resistive relativistic magnetohydrodynamics. In the proposed approach, the spatial components of both magnetic and electric fields are located at zone interfaces and are evolved using the constrained transport formalism. Direct application of Stokes' theorem to Faraday's and Ampere's laws ensures that the resulting discretization is divergence-free for the magnetic field and charge-conserving for the electric field. Hydrodynamic variables retain, instead, the usual zone-centred representation commonly adopted in finite-volume schemes. Temporal discretization is based on Runge-Kutta implicit-explicit (IMEX) schemes in order to resolve the temporal scale disparity introduced by the stiff source term in Ampere's law. The implicit step is accomplished by means of an improved and more efficient Newton-Broyden multidimensional root-finding algorithm. The explicit step relies on a multidimensional Riemann solver to compute the line-averaged electric and magnetic fields at zone edges and it employs a one-dimensional Riemann solver at zone interfaces to update zone-centred hydrodynamic quantities. For the latter, we introduce a five-wave solver based on the frozen limit of the relaxation system whereby the solution to the Riemann problem can be decomposed into an outer Maxwell solver and an inner hydrodynamic solver. A number of numerical benchmarks demonstrate that our method is superior in stability and robustness to the more popular charge-conserving divergence cleaning approach where both primary electric and magnetic fields are zone-centered. In addition, the employment of a less diffusive Riemann solver noticeably improves the accuracy of the computations.

We present a linear mode analysis of the relativistic MHD equations in the presence of finite electrical conductivity. Starting from the fully relativistic covariant formulation, we derive the dispersion relation in the limit of small linear perturbations. It is found that the system supports ten wave modes which can be easily identified in the limits of small or large conductivities. In the resistive limit, matter and electromagnetic fields decouple and solution modes approach pairs of light and acoustic waves as well as a number of purely damped (non-propagating) modes. In the opposite (ideal) limit, the frozen-in condition applies and the modes of propagation coincide with a pair of fast magnetosonic, a pair of slow and Alfv\'en modes, as expected. In addition, the contact mode is always present and it is unaffected by the conductivity. For finite values of the conductivity, the dispersion relation gives rise to either pairs of opposite complex conjugate roots or purely imaginary (damped) modes. In all cases, the system is dissipative and also dispersive as the phase velocity depends nonlineary on the wavenumber. Occasionally, the group velocity may exceed the speed of light although this does not lead to superluminal signal propagation.

We describe an implementation of a particle physics module available for the PLUTO code, appropriate for the dynamical evolution of a plasma consisting of a thermal fluid and a non-thermal component represented by relativistic charged particles, or cosmic rays (CR). While the fluid is approached using standard numerical schemes for magnetohydrodynamics, CR particles are treated kinetically using conventional Particle-In-Cell (PIC) techniques. The module can be used to describe either test particles motion in the fluid electromagnetic field or to solve the fully coupled MHD-PIC system of equations with particle backreaction on the fluid as originally introduced by \cite{Bai_etal.2015}. Particle backreaction on the fluid is included in the form of momentum-energy feedback and by introducing the CR-induced Hall term in Ohm's law. The hybrid MHD-PIC module can be employed to study CR kinetic effects on scales larger than the (ion) skin depth provided the Larmor gyration scale is properly resolved. When applicable, this formulation avoids to resolve microscopic scales offering a substantial computational saving with respect to PIC simulations. We present a fully-conservative formulation which is second-order accurate in time and space and extends to either Runge-Kutta (RK) or corner-transport-upwind (CTU) time-stepping schemes (for the fluid) while a standard Boris integrator is employed for the particles. For highly-energetic relativistic CRs and in order to overcome the time step restriction a novel sub-cycling strategy that retains second-order accuracy in time is presented. Numerical benchmarks and applications including Bell instability, diffusive shock acceleration and test particle acceleration in reconnecting layers are discussed.