In the presence of a sufficiently strong electric field, electrons may break away from thermal equilibrium and approach relativistic speeds. Such “runaway” electrons, common in tokamaks, traverse orbits (trapped or passing) at much faster time scales than collisional ones, while dynamics of interest saturate on time scales much longer than these. Therefore, accurate and efficient modeling of orbit dynamics beyond collisional timescales is essential to model runaway dynamics in tokamaks. Common strategies to deal with this scale separation are based on bounce averaging, but these are brittle when separation between collisional and bounce timescales is lost [1], and are not easily generalizable to arbitrary 3D magnetic field configurations. In this study, we propose a 2D-2P semi-Lagrangian scheme to efficiently bridge these 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 [2], making the method tractable (similarly to what was done in previous work for anisotropic diffusion [3]). While for the time being our implementation considers 2D magnetic fields (e.g., nested flux surfaces), the formulation can be straightforwardly generalized to arbitrary 3D magnetic fields. Our 2P relativistic Fokker-Planck collisional treatment [4] is fully implicit, highly scalable algorithmically and in parallel, and ensures positivity and discrete energy and momentum conservation for the electron-electron interactions. Using a multigrid-preconditioned Anderson Acceleration nonlinear solver, we demonstrate parallel scalability up to 4096 cores and algorithmic scalability with time steps approximately 1000 times larger than explicit ones. The proposed 2D-2P treatment is formally first-order accurate in time, but promises to (i) preserve asymptotic properties associated with stiff Vlasov term, (ii) be uniformly accurate in $\epsilon$, where $\epsilon$ is the ratio of advection to collisional time scales, and (iii) be optimal (i.e. scalable with the total number of mesh points in the domain). We will demonstrate the algorithm with realistic applications of interest [1]. [1] McDevitt and Tang, EPL, 127 (2019) [2] Braams and Karney, PRL, 59 (1987) [3] Chacón et al., JCP, 272 (2014) [4] Daniel, Taitano, Chacón, CPC, submitted (2019)