  This paper proposes a novel numerical integrator for modeling multispecies
Coulomb collisions in kinetic plasmas. The proposed scheme provides an energy-,
momentum-, and positivity-preserving particle discretization of the nonlinear
Landau collision operator, extending the works of J.A. Carrillo et al., Journal
of Computational Physics, 7, 100066 (2020) and E. Hirvijoki, Plasma Physics and
Controlled Fusion, 63, 044003 (2021). The discrete-time conservation properties
are analyzed both algebraically and numerically, and an efficient,
GPU-parallelized implementation is validated against inhomogeneous temperature
relaxation, isotropization and thermalization examples. The results agree with
analytical estimates, confirming the method capable of reproducing physics.
