Summary: Liouville’s equation describes light propagation through an optical system. It governs the evolution of an energy distribution, the basic luminance, on phase space. The basic luminance is discontinuous across optical interfaces, as is described by non-local boundary conditions at each optical interface. Curved optical interfaces manifest themselves as moving boundaries on phase space. A common situation is that an optical system is described by a piecewise constant refractive index field. Away from optical interfaces, the characteristics of Liouville’s equation reduce to straight lines. This property is exploited in the novel solver developed in this paper by employing a semi-Lagrangian discontinuous Galerkin (SLDG) scheme away from optical interfaces. Close to optical interfaces we apply an ADER discontinuous Galerkin (ADER-DG) method on a moving mesh. The ADER-DG method is a fully discrete explicit scheme which must obey a CFL condition that restricts the stepsize, whereas the SLDG scheme can be CFL-free. A moving mesh is used to align optical interfaces with the mesh. Very small elements cannot always be avoided, even when applying mesh refinement. Local time stepping is introduced in the solver to ensure these very small elements only have a local impact on the stepsize. By construction we allow elements of SLDG type to run at a stepsize independent of these small elements. The proposed SLDG scheme uses the exact evolution of the solution, as is described by the characteristics, to update the numerical solution. We impose the condition that no characteristic emanating from an SLDG element can cross an optical interface. In the novel hybrid SLDG and ADER-DG solver this condition is used to naturally divide the spatio-temporal domain into different regions, describing where the SLDG scheme and where the ADER-DG scheme need to be used. An intermediate element is introduced to efficiently couple an SLDG region with an ADER-DG region. Numerical experiments validate that the hybrid solver obeys energy conservation up to machine precision and numerical convergence results show the expected order of convergence. The performance of the hybrid solver is compared to a pure ADER-DG scheme with global time stepping to show the efficiency of the hybrid solver. In particular, a speed-up of 1.6 to 10, in favour of the hybrid solver, for computation times up to 4 minutes was seen in an example. Finally, the hybrid solver and the pure ADER-DG scheme are compared to quasi-Monte Carlo ray tracing. In the examples considered, amongst the three methods the hybrid solver is shown to converge the fastest to high accuracies.