Abstract: For sufficiently strong electric fields, electrons may break away from thermal equilibrium and approach relativistic speeds. Such “runaway” electrons, common in tokamaks, traverse orbits at much faster time scales than collisional ones, while dynamics of interest saturate on time scales much longer than these. In this study, we propose a 2D-2P semi-Lagrangian scheme to efficiently bridge these temporal scales. The approach reformulates the Vlasov equation as an integro-differential operator using Green’s functions along electron orbits (defined by their conserved relativistic energy and magnetic moment), and employs operator splitting to decouple integrals over the relativistic collisional source [1]. We consider 2D magnetic fields (e.g., nested flux surfaces), but the formulation generalizes to arbitrary 3D magnetic fields. Our treatment of the relativistic collisional step is fully implicit, highly scalable, and conservative [2]. The proposed 2D-2P treatment is formally first-order accurate in time, but (i) preserves asymptotic properties associated with stiff Vlasov term, (ii) is uniformly accurate in Δt/ε, where Δt is the timestep and ε is the ratio of advection to collisional time scales, and (iii) is optimal (i.e., scalable with the total number of mesh points in the domain). We will demonstrate the algorithm with applications of interest in realistic tokamak geometries [3]. [1] Chacón, Daniel, Taitano, JCP, 449 (2022)  [2] Daniel, Taitano, Chacón, CPC, 254 (2020)  [3] McDevitt and Tang, EPL, 127 (2019)