  Radiative transfer plays a major role in high-energy astrophysics. In
multiple scenarios and in a broad range of energy scales, the coupling between
matter and radiation is essential to understand the interplay between theory,
observations and numerical simulations. In this paper, we present a novel
scheme for solving the equations of radiation relativistic magnetohydrodynamics
within the parallel code L\'ostrego. These equations, which are formulated
taking successive moments of the Boltzmann radiative transfer equation, are
solved under the gray-body approximation and the M1 closure using an IMEX time
integration scheme. The main novelty of our scheme is that we introduce for the
first time in the context of radiation magnetohydrodynamics a family of
Jacobian-free Riemann solvers based on internal approximations to the
Polynomial Viscosity Matrix, which were demonstrated to be robust and accurate
for non-radiative applications. The robustness and the limitations of the new
algorithms are tested by solving a collection of one-dimensional and
multi-dimensional test problems, both in the free-streaming and in the
diffusion radiation transport limits. Due to its stable performance, the
applicability of the scheme presented in this paper to real astrophysical
scenarios in high-energy astrophysics is promising. In future simulations, we
expect to be able to explore the dynamical relevance of photon-matter
interactions in the context of relativistic jets and accretion discs, from
microquasars and AGN to gamma-ray bursts.
