  Kinetic equations model distributions of particles in position-velocity phase
space. Often, one is interested in studying the long-time behavior of particles
in high-collisional regimes in which an approximate (advection)-diffusion model
holds. In this paper we consider the diffusive scaling. Classical
particle-based techniques suffer from a strict time-step restriction in this
limit, to maintain stability. Asymptotic-preserving schemes avoid this problem,
but introduce an additional time discretization error, possibly resulting in an
unacceptably large bias for larger time steps. Here, we present and analyze a
multilevel Monte Carlo scheme that reduces this bias by combining estimates
using a hierarchy of different time step sizes. We demonstrate how to correlate
trajectories from this scheme, using different time steps. We also present a
strategy for selecting the levels in the multilevel scheme. Our approach
significantly reduces the computation required to perform accurate simulations
of the considered kinetic equations, compared to classical Monte Carlo
approaches.
