4-vector. We recall that $\mathcal{L}_{\boldsymbol{\beta}}$ denotes the Lie derivative along the shift vector and that $\epsilon^{ijk}$ is related to the four-dimensional Levi-Civitá tensor via $\epsilon^{\nu\kappa\lambda} =
\epsilon^{\mu\nu\kappa\lambda}n_{\mu}$ or alternatively $\epsilon^{ijk} =
\eta^{ijk}/\sqrt{\gamma}$, where $\gamma$ now is the determinant of the 3-metric. The scalar fields $\phi,~\psi$ measure the deviation from the constrained solution, with $\phi$ driving the solution of Eq.  toward the zero-divergence condition $\nabla_i B^i=0$, and $\psi$ driving the solution of Eq.  toward the condition $\nabla_i
E^i=q$. This driving is exponentially fast and over a time scale $1/\kappa$. This approach, named hyperbolic divergence cleaning in the context of IMHD, was proposed in Ref. [@Dedner:2002] as a simple way of solving the Maxwell equations and enforcing the conservation of the divergence-free condition for the magnetic field. This method has been extended to the resistive relativistic case in Refs. [@Komissarov2007; @Palenzuela:2008sf].

An obvious consequence of the Maxwell equations in RMHD is the conservation law associated with the electric charge $$\begin{aligned}
\nabla_{\mu} I^{\mu} = 0\,,\end{aligned}$$ which provides an evolution equation for the charge density $$\begin{aligned}
\label{charge}
(\partial_t - \mathcal{L}_{\beta}) q + \nabla_i (\alpha J^i) = \alpha\, K
q\,.\end{aligned}$$ Combining the MHD and Maxwell equations we obtain the following set of evolution equations, which 