  Particle distribution functions evolving under the Lorentz operator can be
simulated with the Langevin equation for pitch angle scattering. This approach
is frequently used in particle based Monte-Carlo simulations of plasma
collisions, among others. However, most numerical treatments do not guarantee
energy conservation, which may lead to unphysical artifacts such as numerical
heating and spectra distortions. We present a novel structure-preserving
numerical algorithm for the Langevin equation for pitch angle scattering.
Similar to the well-known Boris algorithm, the proposed numerical scheme takes
advantage of the structure-preserving properties of the Cayley transform when
calculating the velocity-space rotations. The resulting algorithm is explicitly
solvable, while preserving the norm of velocities down to machine precision. We
demonstrate that the method has the same order of numerical convergence as the
traditional stochastic Euler-Maruyama method. The numerical scheme is
benchmarked by simulating the pitch angle scattering of a particle beam, and
comparing with the analytical solution. Benchmark results show excellent
agreement with theoretical predictions, showcasing the remarkable long time
accuracy of the proposed algorithm.
