Higher-order Time Integration of Coulomb Collisions in a Plasma Using Langevin Equations

We examine the extension of Langevin-equation Monte-Carlo algorithms for Coulomb collisions from the conventional Euler $O$(\textit{$\Delta $t}$^{1/2})$-order time integration to the next higher [Milstein-$O$(\textit{$\Delta $t})] order. In one common Langevin-equation approach, the angular scattering step is treated with a combination of near-Cartesian stochastic velocity-direction kicks, in a unit-vector frame that is rotated so that at the beginning of each timestep, one axis is aligned with the velocity direction. We find that in such schemes, the angular component of the collisional scattering cannot be extended beyond the Euler order. Instead, the extension to higher order proceeds through a formulation of the angular scattering directly as stochastic differential equations in the two fixed-frame spherical coordinates. Such an algorithm involves generation of random numbers that sample the joint distribution function of both the (Gaussian) random coordinate displacements and of double stochastic ``area integrals.'' The sampling of the area integrals can be made using simple but highly accurate approximations to results on ``Levy-area'' processes. Implications for particle simulation of Coulomb collisions in plasmas are discussed.