An equilibrium-preserving discretization for the nonlinear Rosenbluth-Fokker-Planck operator in arbitrary multi-dimensional geometry

The Fokker–Planck collision operator is an advection-diffusion operator which describe dynamical systems such as weakly coupled plasmas [1,2], photonics in high temperature environment [3,4], biological [5], and even social systems [6]. For plasmas in the continuum, the Fokker–Planck collision operator supports such important physical properties as conservation of number, momentum, and energy, as well as positivity. It also obeys the Boltzmann’s H-theorem [7–11], i.e., the operator increases the system entropy while simultaneously driving the distribution function towards a Maxwellian. In the discrete, when these properties are not ensured, numerical simulations can either fail catastrophically or suffer from significant numerical pollution [12,13]. There is strong emphasis in the literature on developing numerical techniques to solve the Fokker–Planck equation while preserving these properties [12–24]. In this short note, we focus on the analytical equilibrium preserving property, meaning that the Fokker–Planck collision operator vanishes when acting on an analytical Maxwellian distribution function. The equilibrium preservation property is especially important, for example, when one is attempting to capture subtle transport physics. Since transport arises from small O ( ) corrections to the equilibrium [25] (where is a small expansion parameter), numerical truncation error present in the equilibrium solution may dominate, overwhelming transport dynamics. Chang and Cooper [15] first proposed a discrete equilibrium preservation scheme for the Fokker–Planck operator in a 1D isotropic linear system. Larsen et al. [14] developed an equilibrium preserving scheme for Compton scattering in 1D by using the analytical expression for the transport coefficients. Buet and Le Thanh [19,20] developed a mass and energy conserving, positivity and equilibrium preserving scheme in the Landau (integral) formulation for a 1D isotropic Fokker–Planck operator by leveraging the natural (integral) structure of the Fokker–Planck collision operator. All of the approaches above