2D Monte Carlo operators describing the Coulomb collisions in pitch angle and radius are studied, which are applicable to axisymmetric toroidal plasmas. The coupling between the spatial and velocity coordinates in toroidal plasmas requires the spatial Jacobian to be included in the collision operator, which is essential to solve the Fokker–Planck equation correctly. The sharp variation of diffusion at the trapped–passing boundary for a tokamak gives rise to differences in the de-trapping probability of the particles into co- and counter-passing orbits, which also produces asymmetry in the distribution function. This problem is solved by using standard and non-standard drift terms, and by symmetrizing the transport across the trapped–passing boundary. Collision operators that relax the distribution function to a prescribed density profile have been developed for simplified models. To obtain converged results, different models are developed and tested that are applicable to diffusion problems with discontinuous diffusion coefficients.