The Gaussian Radial Basis Function Method for Plasma Kinetic Theory
A fundamental macroscopic description of a magnetized plasma is the Vlasov equation supplemented by the nonlinear inverse-square force Fokker-Planck collision operator [Rosenbluth et al., Phys. Rev., 107, 1957]. The Vlasov part describes advection in a six-dimensional phase space whereas the collision operator involves friction and diffusion coefficients that are weighted velocity-space integrals of the particle distribution function. The Fokker-Planck collision operator is an integro-differential, bilinear operator, and numerical discretization of the operator is far from trivial. In this letter, we describe a new approach to discretize the entire kinetic system based on an expansion in Gaussian Radial Basis functions (RBFs). This approach is particularly well-suited to treat the collision operator because the friction and diffusion coefficients can be analytically calculated. Although the RBF method is known to be a powerful scheme for the interpolation of scattered multidimensional data, Gaussian RBFs also have a deep physical interpretation as local thermodynamic equilibria. In this letter we outline the general theory, highlight the connection to plasma fluid theories, and also give 2D and 3D numerical solutions of the nonlinear Fokker-Planck equation. A broad spectrum of applications for the new method is anticipated in both astrophysical and laboratory plasmas.
Motivation A fundamental macroscopic description of a magnetized plasma is the Vlasov equation supplemented by the nonlinear inverse-square force Fokker-planck collision operator Rosenbluth et al. (1957)
where is the distribution of species with charge and mass . The electric and magnetic fields depend on the dsitribution through Maxwell’s equations. This model assumes a statistical description of Coulomb interaction in the limit of small-angle scattering, with the changes in due to collisions with species described by
The friction and diffusion coefficients are given by the expressions
where and is the Coulomb logarithm which respresents a physical cut-off for small-angle collisions. The Rosenbluth potentials appearing in the friction and diffusion coefficients are weighted integrals of the distribution function
and they satisfy the velocity-space Poisson equations and . Expressed in terms of the potential functions, the Fokker-Planck collision operator is
where . A common approach for numerical evaluation of follows a two-phase method where one first inverts the velocity-space Laplacian operators and then directly evaluates the collision operator. Boundary conditions for the Poisson equations can be obtained by limiting the solution to a sub-space and evaluating the expressions at the boundary using a multipole expansion of the potentials Dorf et al. (2014), or by imposing the free-space solution outside the sub-space Pataki and Greengard (2011). Another sophisticated approach is based on fast spectral decomposition as described in Ref. Pareschi et al. (2000). In this letter, we propose a new approach using a mesh-free shifted-Maxwellian representation which is intuitively appealing and straightforward to implement. The solution thus obtained is smooth, extends to , and allows compact representation of any interesting macroscopic quantity (number, momentum, energy density, and so on).
The Gaussian RBF Method To solve the kinetic equation, Eq. (1), we write the total distribution as a finite sum of shifted Maxwellians
where each Maxwellian, , is normalized to unity and the weights allowed to evolve in time. The width parameters and mean velocities can be arbitrary functions of position. The shifted Maxwellians are nothing other than Gaussian Radial Basis functions (RBFs) which have found numerous applications in applied mathematics – in particular for the construction of neural networks Park and Sandberg (1991). For compactness, in what follows we will retain the spatial dependence of all quantities but will not write the dependence explicitly. The potential functions then take the form
where the Gaussian RBF potentials and are defined in terms of the functions and , where is the error function. We thus find a simple bilinear expression for the complete nonlinear operator
where the Gaussian RBF collision tensor is
such that . As we have analytical expressions for , , and , the tensor is easy to implement and fast to evaluate at any point in velocity space. In problems with azimuthal symmetry, a 2D RBF scheme can be developed with axisymmetric ring-like RBF-basis:
where is the order-zero modified Bessel function of the first kind and are the cylindrical velocity-space coordinates. Although explicit expressions for axisymmetric RBF potentials and are not available in a closed form, they can be evaluated numerically to machine precision at any requested point.
Collocation Options We describe two different methods for obtaining an ordinary differential equation for the time evolution of weights: the Galerkin and the center-collocation projections. In the Galerkin method, the kinetic equation – already expanded in the RBF basis – is multiplied by each individual basis function and then integrated over the entire domain. Conversely, in the center-collocation method, the kinetic equation is evaluated at the center of each RBF. Both methods yield equations for the RBF weights, , and result in a differential equation for the weights that can be written in a matrix form. For the moment, let us illustrate the method by considering the spatially homogeneous case with no electromagnetic fields. Then, the matrix equation is
In the Galerkin projection (GP), the matrix is symmetric, typically diagonally dominant, and given by the expression
whereas in the center-collocation (CC) method, the matrix is no longer necessarily symmetric, but can still be dominated by the diagonal components:
On the right-hand-side of Eq. (2), the tensor becomes
for the Galerkin and center-collocation, respectively. Obtaining the center-collocation tensor is merely a matter of evaluating the RBF tensor at the collocation points. Evaluation of the tensor for the Galerkin projection is somewhat more complicated, although the result may be potentially be more accurate or robust. Nevertheless, to maintain simplicity, we focus hereafter on the center collocation-method and omit the CC subscript for brevity. In this case the RBF equations for the full nonlinear system become
where the operator
retains the familiar appearance of the Vlasov operator even though the velocity-space has been completely removed from the problem. Note that depends explicitly on species and implicitly on via the Maxwell equations. In we have defined the temperature-like parameter and also dropped some additional terms that arise if the parameters and depend on position. The RBF-kinetic equation (3) describes collisional fluid-like evolution of the weights in time and space. One physically appealing feature of the RBF expansion is that familiar expressions for macroscopic fluid moments retain their intuitive form:
Another is that Maxwell’s equations can also written compactly:
where for simplicity we have neglected the displacement current. In contrast to the standard fluid approach, where calculation of the higher-order velocity-space moments becomes increasingly cumbersome, the RBF-kinetic approach, however, offers a straightforward and intuitive alternative by describing the velocity space with Gaussian functions that naturally appear in the kinetic theory as thermodynamic equilibrium solutions. Moreover, so long as the RBF spacing is equal to about one e-folding length, the problem remains well-conditioned even for hundreds or thousands of basis functions.
It should also be evident that the Gaussian RBF method is well-suited to linearized models.The linearized collision operator in this case is just the sum of and ; namely, the first row and column of the nonlinear collision tensor. Because the method is mesh-free, it is also natually suited to multi-species plasmas, unlike contemporary algorithms that are based on a velocity mesh Xiong et al. (2008).
Nonlinear Simulations in 2D and 3D To assess the viability of the Gaussian RBF approach, we study a single species plasma using a uniform grid for the collocation points and a fixed value for all RBFs. We also normalize time with the collision time, i.e., . Our problem is thus to solve the nonlinear single species Fokker-Planck equation
We emphasize that this is a highly nontrivial numerical problem. First, we test how well the discretization preserves a (stationary) Maxwellian equilibrium state, where the temperature of the equilibrium is fixed and much warmer than the individual RBFs. For the analytic equilibrium we choose , with and arrange the collocation points into uniform rectilinear grids in spaces and for the full 3D and axisymmetric 2D methods respectively. The equilibrium is projected onto the RBF basis to obtain the numerical equivalent of the steady state distribution function. Then the collision operator in Eq. (4) is evaluated for an increasing number of collocation points and different values of . As an error metric, we choose , presented in Fig. 1. This test is sensitive to the error both in the region interior to the RBF collocation centers, as well as to the region beyond them. As is evident from Fig. 1, the maximum of the global value decreases when the number of collocation points is increased indicating that a numerical equivalent to the analytical steady-state can be reached.
Next, we follow the relaxation of a non-trivial distribution function to the equilibrium state and track the conservation of number, momentum, and energy densities. The initial state of the distribution function is set to , where , , for the components and respectively. The initial state is thus axisymmetric with respect to the axis and the the axisymmetric setting is thus chosen to be . The collocation points are chosen uniformly in the regions and , the initial state projected to the Gaussian RBF basis to get initial values for the weights, and the weights then evolved in time according to the Eq. (2) using a standard fourth-order Runge-Kutta method. From the time-evolution of the weights, we have calculated the evolution of the velocity-space moments, and recorded their maximal deviations from the initial values during the simulation. The results are illustrated in Fig. 2, where one observes good conservation of number (), momentum (velocity components , , and ), and energy () densities, and convergence of their relative error towards zero when the number of RBFs is increased for both full 3D and axisymmetric 2D solvers. Even better conservation properties can be reached by embedding the conservation laws into the time-evolution of the weights. For example, replacing one row in the matrix with ones (unity moments of the basis functions) and setting the corresponding element in to zero neglects information from one of the basis nodes but adds the density conservation. With this approach we observed, e.g., density conservation up to 12 digits without relevant increase in the computation time.
The time evolution of the distribution in the conservation study using 3D RBFs and axisymmetric RBFs is illustrated in Fig. 3 with slices in -plane for the full 3D solver and in -plane for the 2D axisymmetric solver. The time slices are chosen for the initial state () in Figs. a and e, for the beginning of the relaxation process () in Figs. b and f, for the merging phase towards the equilibrium () in Figs. c and g, and for the final state () where the distribution function is close to spherical symmetric and equilibrium in Figs. d and h. Both the 2D and 3D solvers qualitatively describe the same solution to the initially axisymmetric problem and preserve the initial symmetry. As a more quantitative measure, we also follow the time evolution of the -moment of the distribution function which is not a conserved quantity. From the analytical initial state one can calculate the value to be and for the corresponding analytical equilibrium state the value is . The time evolution of -moment from the numerical calculation is shown in Fig. 4. As is clear, the 2D and 3D solutions agree with each other and also confirm our claim that the solution in Figs. d and h is close to the analytical equilibrium state. In other words, we have solved the relaxation problem accurately in both 2D and 3D and demonstrated the capabilities of the RBFs for discretizing the nonlinear collision operator.
Summary In this letter, we have derived a completely new approach to collisional dynamics in plasmas. We have shown how to implement a Gaussian radial-basis-function ansatz to discretize the velocity space, and demonstrated its capabilities in non-linear Fokker-Planck calculations in both full 3D and axisymmetric 2D velocity-space. In a wider connection to kinetic theory we have provided also the discretization of the advection terms and given analytical ansatz expressions for all relevant fluid-moments. In future, our method could be used to solve various problems where non-equilibrium kinetic effects are important, as long as the Landau-approximation for the collision integral is valid.
- M. N. Rosenbluth, W. M. MacDonald, and D. L. Judd, Phys. Rev. 107, 1 (1957).
- M. A. Dorf, R. H. Cohen, M. Dorr, J. Hittinger, and T. D. Rognlien, Contrib. Plasma Phys. 54, 517 (2014).
- A. Pataki and L. Greengard, J. Comput. Phys. 230, 7840 (2011).
- L. Pareschi, G. Russo, and G. Toscani, J. Comput. Phys. 165, 216 (2000).
- J. Park and I. W. Sandberg, Neural Comput. 3, 246 (1991).
- Z. Xiong, R. Cohen, T. Rognlien, and X. Xu, J. Comput. Phys. 227, 7192 (2008).