The integration of the two presented systems requires some specific temporal integration methods. Indeed, these systems can be seen as convective-diffusive systems with stiff source terms [1]}, [2]}, [3]}, as presented in Appendix . Some of these source terms (in particular those associated with collisions and ionization/recombination) are numerically stiff since they are associated with physical processes where timescales often are much smaller than the convective timescales. Therefore, these terms are subject to strong CFL constraints leading to small timesteps and drastically increased computational costs. To avoid this, our strategy is to consider a first order Lie splitting approach in order to split the temporal integration on different subsystems. First, the convective terms have been integrated following a third-order predictor-corrector Hyman [4]}. Alternately, it is possible to consider a third-order explicit Runge-Kutta method. Then, the diffusive term (i.e.,  the Spitzer term) has been integrated by solving a non-Fickian (hyperbolic) diffusion equation as described by [5]}. Note that, in Ebysus, it is possible to integrate this term following an implicit multi-grid method, as described in [6]}. Finally, the stiff source terms (see \(S^{\text{stiff}}\)  in Appendix ) are integrated following an Ordinary Differential Equation (ODE) solver based on a fifth-order implicit Runge Kutta method also known as Radau IIA method [7]}, [8]}. In order to maintain the solenoidal constraint on the magnetic field, a Hodge projection method has been used [9]}.
