  In this study, we consider a numerical implementation of the nonlinear
Rosenbluth-Trubnikov collision operator for particle simulations in plasma
physics in the framework of the finite element method (FEM). The relevant
particle evolution equations are formulated as stochastic differential
equations, both in the Stratonovich and It\^o forms, and are then solved with
advanced high-order stochastic numerical schemes. Due to its formulation as a
stochastic differential equation, both the drift and diffusion components of
the collision operator are treated on an equal footing. Our investigation
focuses on assessing the accuracy of these schemes. Previous studies on this
subject have used the Euler-Maruyama scheme, which, although popular, is of low
order, and requires small time steps to achieve satisfactory accuracy. In this
work, we compare the performance of the Euler-Maruyama method to other
high-order stochastic methods known in the stochastic differential equations
literature. Our study reveals advantageous features of these high-order
schemes, such as better accuracy and improved conservation properties of the
numerical solution. The main test case used in the numerical experiments is the
thermalization of isotropic and anisotropic particle distributions.
