# Spectral scheme for an energetic Fokker-Planck equation with $\kappa$-distribution steady states.

Hugo Parada, Claudia Negulescu

► To cite this version:

Hugo Parada, Claudia Negulescu. Spectral scheme for an energetic Fokker-Planck equation with $\kappa$-distribution steady states.. 2025. hal-05240353

HAL Id: hal-05240353

https://hal.science/hal-05240353v1

Preprint submitted on 4 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# SPECTRAL SCHEME FOR AN ENERGETIC FOKKER-PLANCK EQUATION WITH $\kappa$-DISTRIBUTION STEADY STATES.

CLAUDIA NEGULESCU, HUGO PARADA

**ABSTRACT.** The concern of the present paper is the design of efficient numerical schemes for a specific Fokker-Planck equation describing the dynamics of energetic particles occurring in thermonuclear fusion plasmas (runaway electrons for example). In the long-time limit, the velocity distribution function of these particles tends towards a thermal non-equilibrium $\kappa$-distribution function which is a steady-state of the considered Fokker-Planck equation. These $\kappa$-distribution functions have the particularity of being only algebraically decaying for large velocities, thus describing very well suprathermal particle populations. Our aim is to present two efficient spectral methods for the simulation of such energetic particle dynamics. The first method will be based on rational Chebyshev basis functions, rather than on Hermite basis sets, which are the basis of choice for Maxwellian steady states. The second method is based on a different polynomial basis set, constructed via the Gram-Schmidt orthogonalisation process. These two new spectral schemes, specifically adapted to the here considered physical context, shall permit to cope with the long-time asymptotics without significant numerical costs.

## CONTENTS

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>1.1. Physical motivation of the energetic Fokker-Planck collision operator</td><td>4</td></tr><tr><td>1.2. Aim of this paper and outline</td><td>6</td></tr><tr><td>2. Fokker-Planck collision operator</td><td>8</td></tr><tr><td>2.1. Standard Fokker-Planck equation with $f_{eq} = M$.</td><td>9</td></tr><tr><td>2.2. Energetic Fokker-Planck equation for $f_{eq} = f_{\kappa}$.</td><td>10</td></tr><tr><td>2.3. Regularized energetic Fokker-Planck operator for $f_{eq} = f_{\kappa} e^{-a v^2/2}$</td><td>11</td></tr><tr><td>3. Spectral scheme based on rational Chebyshev functions</td><td>12</td></tr><tr><td>3.1. Rational Chebyshev polynomials</td><td>12</td></tr><tr><td>3.2. RC-spectral method</td><td>14</td></tr><tr><td>4. Spectral scheme based on non-standard orthogonal polynomials</td><td>16</td></tr><tr><td>5. Numerical simulations</td><td>19</td></tr><tr><td>5.1. Numerical simulations with a standard FD-method</td><td>20</td></tr><tr><td>5.2. Numerical results obtained with our spectral methods</td><td>23</td></tr><tr><td>6. Concluding remarks and perspectives</td><td>30</td></tr><tr><td>References</td><td>31</td></tr></table>

*Date:* September 4, 2025.

*Key words and phrases.* **Keywords:** Fusion plasma modelling, Fokker-Planck kinetic equation, energetic particles, kappa-distribution function, (thermal) non-equilibrium steady-states, asymptotic analysis, spectral methods, rational Chebyshev polynomials.
# 1. INTRODUCTION

Modelling the dynamics of particles which are not in thermodynamic equilibrium is a very interesting and intricate problem. Plasmas encountered in astrophysics [6, 51, 58] or fusion devices (tokamak) [39] are usually collisionless and do not attain thermal equilibrium, described by a Maxwellian velocity distribution function. Instead, such plasmas contain highly energetic particles, which are well-represented by so-called "$\kappa$-distribution functions", i.e. power-law distributions with heavy tails. The departure from thermal equilibrium is a consequence of an out-of-equilibrium phenomenon, where the driving (diffusive) force and the drag (dissipation) force do not satisfy any more the fluctuation-dissipation theorem [46], leading thus to velocity distributions with energetic tails, meaning an over-population of particles with very high kinetic energies.

The physical reason for the departure from Maxwellian equilibria is linked to the scattering cross-section dependence on the velocity, $\sigma(v)$, which is proportional to $v^{-4}$, meaning that while low-energy particles are often collisional, faster ones are generally not anymore, their slowing-down through Coulombian collisions becoming thus problematic at high speeds. An anomalous acceleration process, arising for example through turbulences or Dreicer electric fields, is then sufficient to accelerate them to infinite speeds (non-relativistic picture). In a relativistic framework, the frictional force does not vanish for $v \to \infty$, but increases slowly with higher energies, permitting finally to bound the particle speed by the speed of light. For more physical details about all these phenomena we refer the interested reader to [37, 49, 64].

Motivated by this rich physical context, we shall be concerned in this work with the mathematical description and numerical simulation of the energetic particle dynamics via a specific energetic Fokker-Planck equation. An accurate description of the energetic particle dynamics is of crucial importance for fusion reactor performances, in particular for confinement/heating reasons (α-particles) or stability reasons (runaway electrons). Runaway electrons are the main focus of this paper, their occurrence having an essential impact on the global stability of the fusion thermal plasma, such that their control is of fundamental importance.

$\kappa$-distributions were introduced for the first time in space physics by Olbert [55] and Vasyliunas [67] to describe suprathermal populations. Such distributions have a core of Maxwellian type and a power-law (thus not exponential) decline in the large velocity ranges (see Fig. 1), meaning they describe an over-population of particles at high speeds, when compared to Maxwellian distributions. Several definitions of $\kappa$-distributions exist in literature, characterizing the power-law nature of the energetic tails in different manners. The $\kappa$-distribution of first kind we shall consider here, is given by the formula [6]

$$f_{\kappa}(t, \mathbf{x}, \mathbf{v}) = A_{\kappa} n(t, \mathbf{x}) \left(1 + \frac{|\mathbf{v}|^2}{2 \kappa v_{th}^2}\right)^{-\kappa}, \quad A_{\kappa} := \frac{1}{(2 \pi \kappa v_{th}^2)^{3/2}} \frac{\Gamma(\kappa)}{\Gamma(\kappa - 3/2)}, \quad (1.1)$$

where the thermal speed $v_{th} := \sqrt{\frac{k_B T}{m}}$ and the particle density $n$ are fixed by the associated Maxwellian $M$, given by

$$M(t, \mathbf{x}, \mathbf{v}) := \frac{n(t, \mathbf{x})}{(2 \pi v_{th}^2)^{3/2}} e^{-\frac{|\mathbf{v}|^2}{2 v_{th}^2}}, \quad n := \int_{\mathbb{R}^3} M d\mathbf{v} = \int_{\mathbb{R}^3} f_{\kappa} d\mathbf{v}, \quad \frac{3}{2} n k_B T := \frac{m}{2} \int_{\mathbb{R}^3} |\mathbf{v}|^2 M d\mathbf{v}.$$

Let us recall here for simplicity the Euler-Gamma function $\Gamma$, defined for $r \in \mathbb{R}_*^+$ by

$$\Gamma(r) := \int_{0}^{\infty} s^{r-1} e^{-s} ds, \quad \Gamma(1) = 1, \quad \Gamma(1/2) = \sqrt{\pi}, \quad \Gamma(r+1) = r \Gamma(r),$$
and satisfying $\Gamma(r + \alpha) \sim_{r \to \infty} \Gamma(r) r^\alpha$ for all $\alpha \in \mathbb{R}^+$, where $\sim$ stands for the asymptotic equivalence. Remark that $f_\kappa$ contains three parameters, $v_{th}$ and $n$ which are linked to the associated Maxwellian, and $3/2 < \kappa < \infty$ which is the only free parameter and determines the distance away from the Maxwellian, in particular the lower the $\kappa$, the more pronounced are the tails. In the limit of large $\kappa \gg 1$ one recovers the standard Maxwellian distribution, i.e. one can easily prove that

$$ \lim_{\kappa \to \infty} f_{\kappa}(\mathbf{v}) = \mathcal{M}(\mathbf{v}), \quad \forall \mathbf{v} \in \mathbb{R}^3, $$

using

$$ e^{-\xi} = \lim_{\kappa \to \infty} \left(1 + \frac{\xi}{\kappa}\right)^{-\kappa}, \quad \lim_{\kappa \to \infty} A_{\kappa} = \frac{1}{(2\pi v_{th}^2)^{3/2}}. $$

One major problem with $\kappa$-distributions is that not all the moments are bounded, contrary to Maxwellian distributions, and this is due to the power-law decrease in the velocity-variable. Indeed, one has

$$ \int_{\mathbb{R}^3} |\mathbf{v}|^l f_{\kappa}(\mathbf{v}) d\mathbf{v} = \frac{2n}{\sqrt{\pi}} (2\kappa v_{th}^2)^{l/2} \frac{\Gamma\left(\frac{l+3}{2}\right) \Gamma\left(\kappa - \frac{l+3}{2}\right)}{\Gamma\left(\kappa - \frac{3}{2}\right)}, \quad \forall 0 \le l < 2\kappa - 3, $$

higher moments being divergent, a fact which alters a characterization of the $\kappa$-distribution function at the macroscopic level (for ex. via temperature, pressure, heat flux, etc.). To avoid this deficiency, one idea is to regularize the $\kappa$-distribution function by multiplying it with some exponentially decaying function, introducing for example $f_{\kappa,a}(\mathbf{v}) := f_{\kappa}(\mathbf{v}) e^{-a \frac{|\mathbf{v}|^2}{2v_{th}^2}}$, with $a \ll 1$ a small parameter permitting to regulate the cut-off of the $\kappa$-distribution at sufficiently high speeds. The regularized $\kappa$-distribution has now a power-law behaviour at low and intermediate velocities, and an exponential cut-off at higher velocities. Besides the fact that this regularized distribution has now finite moments at all orders, this strategy is somehow also very reasonable from a physical point of view, as it can be seen as a saturation in the velocity variable, avoiding the occurrence of particles with speeds higher than the speed of the light in vacuum.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. 1D first kind $\kappa$-distribution functions for several $\kappa$-values and associated Maxwellian distribution, recovered for $\kappa \to \infty$.</figcaption></figure>
Let us restrict in the following, for simplicity reasons, to the one-dimensional case $v \in \mathbb{R}$ and note that $\kappa$-distribution functions are steady state solutions of specific, energetic Fokker-Planck equations, where the drift and diffusion coefficients do not satisfy any more the fluctuation-dissipation theorem.

**1.1. Physical motivation of the energetic Fokker-Planck collision operator.** Let us explain in more details in this subsection how $\kappa$-distributions arise. The basic physical origin is the runaway phenomenon, which is a consequence of the fact that due to the strong energy-dependence of the Coulombian cross-section $\sigma(v)$, varying as $v^{-4}$, the dynamical friction (resistance due to long-range Coulomb interactions) of the particles cannot balance any more the turbulent acceleration forces (or some other acceleration forces not considered here and due to the electric field), resulting in a continuous acceleration of the particles if nothing else slows them down. To translate this in mathematics, let us proceed step by step.

From a microscopic point of view, the dynamics of a particle of mass $m$ is described via Langevin's equations

$$mv'(t) = -m\gamma(v)v + m\sqrt{2D(v)}\eta_t, \quad \forall t > 0,$$

where $\gamma(v) = v\sigma(v) \sim_{|v|\ge 1} v^{-3}$ is the friction coefficient and $D(v)$ the diffusion coefficient in the velocity variable. The second term on the right-hand side stands for the random force, assumed to be a Gaussian process (white noise) with infinitely short correlation times. The corresponding mesoscopic picture of Langevin's equation is the kinetic Fokker-Planck equation

$$\partial_t f = \partial_v [\gamma(v) v f + D(v) \partial_v f], \quad \forall (t, v) \in \mathbb{R}^+ \times \mathbb{R}.$$

At thermal equilibrium, the friction coefficient is related to the diffusion coefficient via $\frac{D(v)}{\gamma(v)} = \frac{k_B T}{m}$, which is the so-called fluctuation-dissipation theorem [46] and states that the friction is finally the mean by which an external work is dissipated into microscopic thermal energy. This relation leads then to the well-known form of the linear Fokker-Planck equation

$$\partial_t f = \partial_v [\gamma(v) (v f + v_{th}^2 \partial_v f)], \quad \text{with} \quad v_{th}^2 = \frac{k_B T}{m}, \quad \forall (t, v) \in \mathbb{R}^+ \times \mathbb{R}. \quad (1.2)$$

Remark that a transport term $v \partial_x f + E \partial_v f$ shall be introduced in this equation on the left-hand side, in order to get closer to a real thermonuclear fusion simulation with convection, for the moment we shall however only focus on the Fokker-Planck collision operator. Equation (1.2) can be recast into the more general one-dimensional Fokker-Planck form

$$\partial_t f = \nu \partial_v [D(v) (U'(v) f + \partial_v f)], \quad \forall (t, v) \in \mathbb{R}^+ \times \mathbb{R}, \quad (1.3)$$

where $\nu > 0$ is the collisional frequency, $D(v) \ge 0$ is the velocity dependent diffusion coefficient and $U(v)$ a potential function describing the drift mechanism. The stationary solution of (1.3) is given by

$$f_{\infty}(v) := \alpha e^{-U(v)}, \quad \forall v \in \mathbb{R},$$

with $\alpha > 0$ determined by the initial condition $f_{in}$ (conservation of mass).

The standard case (1.2) is recovered from (1.3) via $D(v) = \gamma(v) v_{th}^2$ and a quadratic potential $U(v) := \frac{v^2}{2v_{th}^2}$, leading to the usual thermal equilibrium $M(v) = \frac{n}{\sqrt{2\pi v_{th}^2}} e^{-\frac{v^2}{2v_{th}^2}}$ (Maxwellian distribution).
The fluctuation-dissipation theorem gets however violated in some situations, for example when the friction term cannot be any more balanced by the accelerating diffusion term, leading to an out-of-equilibrium situation characterized for example by heavy-tail $\kappa$-distributions. In mathematical terms, these $\kappa$-distributions arise by introducing in (1.2) an additional diffusive mechanism in the Fokker-Planck operator, such as

$$ \partial_t f = \nu \partial_v \{ \gamma(v) [v f + (v_{th}^2 + \bar{\mathcal{K}}_{turb}(v)) \partial_v f] \}, \quad \forall (t, v) \in \mathbb{R}^+ \times \mathbb{R}, \quad (1.4) $$

which includes now additionally to the Coulomb collisions with a dense plasma background, an empirical turbulent acceleration mechanism, described by the coefficient $\bar{\mathcal{K}}_{turb}(v)$. The idea is that the supplementary diffusion term injects energy into the system at a rate $\bar{\mathcal{K}}_{turb}(v)$. This mechanism of energy injection, combined with the friction coefficient $\gamma(v)$, determines together the shape of the equilibrium distribution function at large velocities.

One can now put this new Fokker-Planck equation under the form (1.3) with the diffusion coefficient $D(v) = \gamma(v) (v_{th}^2 + \bar{\mathcal{K}}_{turb}(v))$ and $U'(v) = v (v_{th}^2 + \bar{\mathcal{K}}_{turb}(v))^{-1}$. For large $\bar{\mathcal{K}}_{turb}(v)$, for example if $\bar{\mathcal{K}}_{turb}(v) = \bar{\mathcal{K}}_* v^2$ with $\bar{\mathcal{K}}_* > 0$, one gets a potential of logarithmic type $U(v) = \kappa \ln(1 + \frac{|v|^2}{2\kappa v_{th}^2})$, and the following $\kappa$-distribution functions as an equilibrium

$$ f_{\kappa}(v) = C_{\kappa} \left( 1 + \frac{v^2}{2 \kappa v_{th}^2} \right)^{-\kappa}, \quad \kappa := \frac{1}{2 \bar{\mathcal{K}}_{*}}, \quad v_{th} := \sqrt{\frac{k_B T}{m}}, \quad (1.5) $$

with constant $C_\kappa > 0$.

The problem now is that this is a non-relativistic picture, the particles being able to reach infinite speeds, which is from a physical point of view not so justified. To cope with this new problem, one has to introduce a saturation term, in order to restrain the speeds to speeds lower than the speed of light. From a physical point of view this arises naturally, as for very high electron speeds the friction coefficient starts to increase slowly, which is due mainly to the so-called Bremsstrahlung-Effekt [17, 64]. Starting from some threshold energy, the particles start to emit electromagnetic waves, such that their energy is dissipated, and their velocity cannot attain the speed of the light. This phenomenon can be described mathematically by starting from the previously introduced turbulent Fokker-Planck equation (1.4), rewritten here as

$$ \partial_t f = \nu \partial_v \left\{ D(v) \left[ (v_{th}^2 + \bar{\mathcal{K}}_{turb}(v))^{-1} v f + \partial_v f \right] \right\}, \quad (1.6) $$

and adding a new term $a v f$ in the friction part. This leads to

$$ \partial_t f = \nu \partial_v \left\{ D(v) \left[ \left( (v_{th}^2 + \bar{\mathcal{K}}_{turb}(v))^{-1} + a \right) v f + \partial_v f \right] \right\}, \quad (1.7) $$

recalling that $D(v) = \gamma(v) (v_{th}^2 + \bar{\mathcal{K}}_{turb}(v))$. Thus $U'(v) = v \left[ (v_{th}^2 + \bar{\mathcal{K}}_{turb}(v))^{-1} + a \right]$, and choosing again $\bar{\mathcal{K}}_{turb}(v) = \bar{\mathcal{K}}_* v^2$ leads to the potential $U(v) = \kappa \ln(1 + \frac{|v|^2}{2\kappa v_{th}^2}) + a \frac{v^2}{2}$ and the corresponding equilibria

$$ f_{\kappa,a}(v) := N_{\kappa} \left( 1 + \frac{v^2}{2 \kappa v_{th}^2} \right)^{-\kappa} e^{-a \frac{v^2}{2}}, \quad \kappa := \frac{1}{2 \bar{\mathcal{K}}_{*}}, $$

with a constant $N_\kappa > 0$. Equation (1.7) can be rewritten via these equilibria as

$$ \partial_t f = \nu \partial_v \left\{ D(v) f_{\kappa,a}(v) \partial_v \left( \frac{f}{f_{\kappa,a}} \right) \right\}, \quad \forall v \in \mathbb{R}. \quad (1.8) $$
The new friction force reads now $F_{fric} = -m\gamma(v)v[1 + a(v_{th}^2 + \bar{\kappa}_{turb}(v))]$, which means that the friction does not vanish any more for $|v| \to \infty$, but saturates at a constant value, as illustrated in Fig. 2, where we recall that the Coulombian friction satisfies $\gamma(v) \sim_{|v| \gg 1} v^{-3}$, whereas the turbulent diffusion has the form $\bar{\kappa}_{turb}(v) = \bar{\kappa}_{\star}v^2$.

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Friction force $F_{fric} = -mv\eta(v)$ with $\eta(v) := \gamma(v)[1 + a(v_{th}^2 + \bar{\kappa}_{turb}(v))]$ and $\gamma(v) = G(v)/v$, where G is the Chandrasekhar function [49]. Coulombian friction for $a = 0$ (blue dotted curve, vanishing at infinity); turbulent friction for $a = 0.012$ (red full curve).</figcaption></figure>

## 1.2. Aim of this paper and outline.
Having said some words about the physical motivation of this work, let us mention now some words about the mathematical goal of the present paper. The main objective is the efficient numerical resolution of the following 1D evolution problem

$$ \partial_t f = \partial_v \left[ f_{eq}(v) \partial_v \left( \frac{f(t, v)}{f_{eq}(v)} \right) \right], \quad \forall (t, v) \in \mathbb{R}_+ \times \mathbb{R}, \qquad (1.9) $$

associated with some initial condition $f_{in}$. Three cases shall be considered and compared, associated to the following three equilibria:

(i) Standard Fokker-Planck case: $f_{eq} = \mathcal{M}(v)$;

(ii) Energetic Fokker-Planck case: $f_{eq} = f_{\kappa}(v)$;

(iii) Regularized energetic Fokker-Planck case: $f_{eq} = f_{\kappa,a}(v) := f_{\kappa}(v)e^{-a\frac{v^2}{2}}$, where $a \ll 1$ is a small cut-off parameter, permitting to avoid supra-luminal particles.

Let us remark here that, contrary to the physical part, we set $D(v) \equiv 1$ for simplicity reasons (compare with (1.8)), the main difficulties being still present in this simplified framework. The mathematical study of the general Fokker-Planck equation has been investigated in [49]. Furthermore, let us recall the Maxwellian equilibrium, which is defined as

$$ \mathcal{M}(v) := \frac{1}{\sqrt{2\pi}} e^{-v^2/2}, \quad \forall v \in \mathbb{R}, \qquad (1.10) $$

whereas the $\kappa$-equilibrium is given for $\kappa > 1/2$ and $v \in \mathbb{R}$ by the formula

$$ f_{\kappa}(v) := \frac{c_{\kappa}}{\left(1 + \frac{v^2}{2\kappa}\right)^{\kappa}}, \quad \text{where} \quad c_{\kappa} := \frac{1}{\sqrt{2\pi\kappa}} \frac{\Gamma(\kappa)}{\Gamma(\kappa - 1/2)}, \qquad (1.11) $$
such that $\int_{\mathbb{R}} f_{\kappa}(v) dv = 1$. The cut-off parameter $0 < a \ll 1$ shall be chosen in such a manner to minimize the unrealistic effect of supra-luminal particles, or in other words to restrain the particle velocities to speeds lower than the speed of the light.

One remarks that

$$f_{\kappa,a} \xrightarrow{a \to 0} f_{\kappa} \xrightarrow{\kappa \to \infty} \mathcal{M},$$

and we recall, for later use, the moments of the $\kappa$-distribution function

$$\int_{\mathbb{R}} |v|^l f_{\kappa}(v) dv = \frac{1}{\sqrt{\pi}} (2\kappa)^{l/2} \frac{\Gamma(\frac{l+1}{2}) \Gamma(\kappa - \frac{l+1}{2})}{\Gamma(\kappa - \frac{1}{2})}, \quad \forall 0 \le l < 2\kappa - 1,$$

as well as those of the regularized $\kappa$-distribution, this time for all $l \in \mathbb{N}$

$$\begin{align}
\int_{\mathbb{R}} |v|^l f_{\kappa,a}(v) dv &= \frac{2^{\frac{l+1}{2}} a^{-\frac{l+1}{2}}}{\Gamma(\kappa)} E(\kappa; b; \frac{l+1}{2}; b; a\kappa) \tag{1.12} \\
&= 2^{\frac{l+1}{2}} \Gamma\left(\frac{l+1}{2}\right) a^{\frac{\kappa}{2} - \frac{l+1}{4} - \frac{1}{2}} \kappa^{\frac{\kappa}{2} + \frac{l+1}{4} - \frac{1}{2}} e^{\frac{a\kappa}{2}} W_{\frac{1}{2} - \frac{\kappa}{2} - \frac{l+1}{4}, \frac{\kappa}{2} - \frac{l+1}{4}}(a\kappa), \nonumber
\end{align}$$

where $E(\cdot)$ is the Mac Robert's E-function with $b \in \mathbb{R}$ an arbitrary value, and $W(\cdot)$ is the Whittaker function [72].

The first immediately visible difficulty in a numerical resolution of (1.9) is related to the unboundedness of the velocity-domain, which leads to challenging numerical complications. Several possibilities are available to treat such unbounded problems, for example spectral methods, domain truncation, mapping techniques (see [9] for more details). We shall concentrate in the present paper on spectral methods, based on basis functions intrinsic to unbounded domains, such as for example Hermite-basis functions, Sinc-functions, rational Chebyshev functions, etc. Spectral methods are so-called global methods and, when well-designed, are superieur in performance to local methods, such as standard finite-difference or finite-element methods, due to their exponential accuracy (greater precision with smaller number of points or modes). This exponential accuracy is however only achieved when the basis functions are well-chosen, an inadequate choice can have a drastic impact on the convergence rate.

The mathematical and numerical treatment of the three problems (i)-(iii) in view of the construction of a spectral scheme, is completely different, and this mainly due to the fact that the spectrum of the standard case (i) is discrete, the eigenvectors being well-known (Hermite basis functions), whereas case (ii) admits a continuous spectrum, requiring a continuous spectral representation with generalized eigenvectors not belonging to the functional space of interest. A basis not constituted of eigenvectors is then of rescue in this case (ii). The equilibrium $\kappa$-distributions are however not exponentially decaying in the velocity variable, but rather algebraically, such that the usual Hermite polynomials (useful in case (i)) will not be efficient here. When a function $f(t, v)$ decays as an inverse power of the velocity-variable, then rational Chebyshev polynomials are clearly more adequate than the exponentially decaying Hermite polynomials, such that our spectral method for case (ii) shall be based on the introduction of a rational Chebyshev basis set.

Finally, case (iii) admits again (due to the cut-off in the high velocity regimes) a discrete spectrum, however with non-standard eigenvectors. The regularization allows however now to employ the function $f_{\kappa,a}$ as a valid weight function, to start the construction (via Gram-Schmidt orthogonalisation process) of an adequate polynomial basis set, permitting by this manner to introduce an efficient spectral method for case (iii).
The outline of this paper is the following. In Section 2 we state the mathematical problem and present some of its properties (spectrum, asymptotic decay). Section 3 is concerned with the construction of a spectral scheme for the resolution of the Fokker-Planck equation in case (ii), scheme based on rational Chebyshev polynomials (so-called RC-scheme). Section 4 focuses on the regularized Fokker-Planck equation of case (iii) and proposes a spectral method based on a polynomial basis set, constructed via the weight function $f_{\kappa,a}$ (so-called GS-scheme). Finally, Section 5 presents the numerical results obtained with both spectral schemes and compares them with a standard finite-difference method (FD-scheme).

## 2. FOKKER-PLANCK COLLISION OPERATOR

The main objective of this section is the mathematical study of the following 1D evolution problem

$$ (FP) \quad \left\{ \begin{array}{l} \partial_t f = \partial_v \left[ f_{eq}(v) \partial_v \left( \frac{f}{f_{eq}} \right) \right], \quad \forall (t,v) \in \mathbb{R}_+ \times \mathbb{R}, \\ f(0, \cdot) = f_{in}(\cdot), \end{array} \right. \qquad (2.1) $$

where for the equilibrium distribution function we choose either $f_{eq} = \mathcal{M}(v)$ defined in (1.10), $f_{eq} = f_{\kappa}(v)$ defined in (1.11) or the regularized version $f_{eq} = f_{\kappa,a}(v)$. In the long-time limit $t \to \infty$ the solution to this equation will tend towards the stationary solution given by $f_{\infty}(v) := \frac{\langle f_{in} \rangle}{\langle f_{eq} \rangle} f_{eq}(v)$, where in the following we shall use the simplified notation for the average in the velocity space

$$ \langle f \rangle := \int_{\mathbb{R}} f(v) \, dv. $$

Let us denote the second order differential operator occurring in (2.1) by

$$ \mathcal{L}_{eq}(f) := -\partial_v \left[ f_{eq} \partial_v \left( \frac{f}{f_{eq}} \right) \right]. $$

Then one can show that $\mathcal{L}_{eq}: \mathcal{D}(\mathcal{L}_{eq}) \subset L_{f_{eq}^{-1}}^2 \to L_{f_{eq}^{-1}}^2$ is a linear, unbounded, self-adjoint and positive operator on the Hilbert-space

$$ L_{f_{eq}^{-1}}^2 := \left\{ f : \mathbb{R} \to \mathbb{R} \text{ measurable}, \quad \int_{\mathbb{R}} |f|^2 f_{eq}^{-1} \, dv < \infty \right\}, \quad (f,g)_{f_{eq}^{-1}} := \int_{\mathbb{R}} fg f_{eq}^{-1} \, dv, \qquad (2.2) $$

and with domain

$$ \mathcal{D}(\mathcal{L}_{eq}) = \{f \in L_{f_{eq}^{-1}}^2, \quad \mathcal{L}_{eq}(f) \in L_{f_{eq}^{-1}}^2\}. \qquad (2.3) $$

Introducing the Hilbert space and the evolution triplet

$$ H_{f_{eq}^{-1}}^{1} := \{f \in L_{f_{eq}^{-1}}^{2} / \partial_{v} f \in L_{f_{eq}^{-1}}^{2}\}, \quad H_{f_{eq}^{-1}}^{1} \hookrightarrow L_{f_{eq}^{-1}}^{2} \hookrightarrow (H_{f_{eq}^{-1}}^{1})^{*} $$

one can restrict $\mathcal{L}_{eq}$ to the linear and bounded operator

$$ \mathcal{L}_{eq}: H_{f_{eq}^{-1}}^{1} \to (H_{f_{eq}^{-1}}^{1})^{*}. $$

Standard functional analysis (Lions theorem, Hille-Yosida theorem) leads then to the following existence/uniqueness result.

**Proposition 1. (Existence/uniqueness of a solution)** [49] For each $f_{in} \in \mathcal{D}(\mathcal{L}_{eq})$ there exists a unique solution $f \in C^1([0, \infty); L_{f_{eq}^{-1}}^2) \cap C([0, \infty); \mathcal{D}(\mathcal{L}_{eq}))$ of the Fokker-Planck equation (2.1). For less regular initial conditions $f_{in} \in L_{f_{eq}^{-1}}^2$ and each $T > 0$, one has nonetheless a unique weak solution $f \in W_2^1(0, T; H_{f_{eq}^{-1}}^1, L_{f_{eq}^{-1}}^2) \subset C([0, T]; L_{f_{eq}^{-1}}^2)$.
Before going on, let us remark here, that one has the inclusions

$$L_{\mathcal{M}^{-1}}^2 \subset L_{f_{\kappa,a}^{-1}}^2 \subset L_{f_{\kappa}^{-1}}^2 \subset L^2(\mathbb{R}).$$

**2.1. Standard Fokker-Planck equation with $f_{eq} = \mathcal{M}$**. This case is very well documented, we recall here the results only for the sake of completeness. The equation we shall consider here is simply

$$(FP)_{\mathcal{M}} \quad \left\{ \begin{array}{l} \partial_t f = \partial_v \left[ \mathcal{M} \partial_v \left( \frac{f}{\mathcal{M}} \right) \right] = \partial_v [v f + \partial_v f], \quad \forall (t,v) \in \mathbb{R}_+ \times \mathbb{R}, \\ f(0, \cdot) = f_{in}(\cdot), \end{array} \right. \quad (2.4)$$

with the operator

$$\mathcal{L}_{\mathcal{M}}(f) := -\partial_v \left[ \mathcal{M} \partial_v \left( \frac{f}{\mathcal{M}} \right) \right]. \quad (2.5)$$

**Proposition 2.** [49, 65] (Hermite basis functions) The operator (2.5) is self-adjoint, positive and with compact resolvent, such that its spectrum is discrete, real, positive and consists of a sequence of eigenvalues $(\lambda_k)_{k \in \mathbb{N}} \subset \mathbb{R}$ satisfying $\lambda_k \to \infty$ as $k \to \infty$. In particular one has for all $k \in \mathbb{N}$

$$\mathcal{L}_{\mathcal{M}} \psi_k(v) = \lambda_k \psi_k(v), \quad \forall v \in \mathbb{R}, \quad \text{with eigenvalue } \lambda_k := k,$$

and the associated eigenvectors are the Hermite functions, defined recursively via $\psi_0 \equiv \mathcal{M}$, $\psi_1 \equiv v \mathcal{M}$ and for all $k \ge 1$ by the formulae

$$\sqrt{k+1} \psi_{k+1}(v) = v \psi_k(v) - \sqrt{k} \psi_{k-1}(v), \quad \forall v \in \mathbb{R}. \quad (2.6)$$

Remark that $\psi'_k(v) = -\sqrt{k+1} \psi_{k+1}(v)$ and that $\{\psi_k\}_{k \in \mathbb{N}}$ form a complete orthonormal basis set of the space $L_{\mathcal{M}^{-1}}^2$.

This proposition is useful both from an analytic point of view, permitting to study the decay rate towards the equilibrium, but also from a numerical perspective, permitting the construction of a spectral method. This is illustrated in the two following propositions.

**Proposition 3.** [49, 65] (Time decay) Let $f$ be the unique solution of the evolution problem (2.4), with initial condition $f_{in} \in L_{\mathcal{M}^{-1}}^2$. Then the following exponential time-decay of the solution towards its stationary state holds

$$\|f(t) - f_{\infty}\|_{\mathcal{M}^{-1}} \le \|f_{in} - \langle f_{in} \rangle \mathcal{M}\|_{\mathcal{M}^{-1}} e^{-t}, \quad \forall t \ge 0, \quad (2.7)$$

where the stationary solution is given by $f_{\infty}(v) = \langle f_{in} \rangle \mathcal{M}(v)$.

Now, one can make use of the orthonormal basis of eigenvectors $\{\psi_k\}_{k \in \mathbb{N}} \subset L_{\mathcal{M}^{-1}}^2$ to expand the solution of the evolution problem (2.4).

**Proposition 4. (Spectral representation)** The solution of (2.4) can be expressed in terms of the Hermite functions introduced in Proposition 2 as follows

$$f(t, v) = \sum_{k=0}^{\infty} \alpha_k(t) \psi_k(v) \quad \text{with} \quad \alpha_k(t) = (f(t), \psi_k)_{\mathcal{M}^{-1}} = \int_{\mathbb{R}} f(t, v) \psi_k(v) \mathcal{M}^{-1}(v) \, dv. \quad (2.8)$$

Inserting this expression in (2.4) yields an equation to be solved for the expansion coefficients $\alpha_k(t)$, leading for all $t \ge 0$ and $k \in \mathbb{N}$ to

$$\alpha'_{k}(t) + \lambda_{k} \alpha_{k}(t) = 0 \quad \Rightarrow \quad \alpha_{k}(t) = e^{-kt} \alpha_{\text{in},k}, \quad \alpha_{\text{in},k} := \int_{\mathbb{R}} f_{\text{in}}(v) \psi_{k}(v) \mathcal{M}^{-1}(v) \, dv. \quad (2.9)$$
This proposition is the starting point of numerical spectral methods for standard Fokker-Planck equations.

Let us remark firstly that, as the basis functions $\psi_k$ are eigenfunctions of the Fokker-Planck operator $\mathcal{L}_M$, the Fokker-Planck equation (2.4) can be diagonalized via the spectral decomposition, leading to the simple ODE system (2.9). Secondly, one observes immediately that $\alpha_0 \equiv \langle f_{in} \rangle$ and that the estimate (2.7) is a simple consequence of the exponential decay of the Hermite coefficients (2.9). Indeed, one has

$$ \|f(t) - f_{\infty}\|_{M^{-1}}^2 = \| \sum_{k=1}^{\infty} \alpha_k(t) \psi_k(v) \|_{M^{-1}}^2 = \sum_{k=1}^{\infty} |\alpha_k(t)|^2 = \sum_{k=1}^{\infty} e^{-2kt} |\alpha_{in,k}|^2 \le e^{-2t} \sum_{k=1}^{\infty} |\alpha_{in,k}|^2, $$

where the last sum is nothing else than $\|f_{in} - \langle f_{in} \rangle \mathcal{M}\|_{M^{-1}}^2$.

Finally, one can similarly estimate the truncation error when using only the first terms in the Hermite expansion (2.8), yielding

$$ \|f^N(t) - f_\infty\|_{M^{-1}}^2 \le e^{-2t} \|f_{in}^N - \langle f_{in} \rangle \mathcal{M}\|_{M^{-1}}^2, $$

where we denoted the truncation by $f^N(t, v) = \sum_{k=0}^N \alpha_k(t) \psi_k(v)$.

## 2.2. Energetic Fokker-Planck equation for $f_{eq} = f_\kappa$. We are now coming to the investigation of the energetic Fokker-Planck equation, which writes

$$ (FP)_\kappa \quad \left\{ \begin{array}{l} \partial_t f = \partial_v \left[ f_\kappa(v) \partial_v \left( \frac{f}{f_\kappa(v)} \right) \right] = \partial_v \left[ \frac{v}{1 + \frac{v^2}{2\kappa}} f + \partial_v f \right], \quad \forall (t,v) \in \mathbb{R}^+ \times \mathbb{R}, \\ f(0, \cdot) = f_{in}(\cdot), \end{array} \right. \qquad (2.10) $$

where the equilibrium is this time $f_\kappa(v) = c_\kappa (1 + \frac{v^2}{2\kappa})^{-\kappa}$, with $c_\kappa > 0$ the normalization constant. We shall work now in the Hilbert-space $L_{f_\kappa}^{2, -1}$ and shall take in this study some fixed $\kappa \in \mathbb{R}^+$.

The linear, unbounded operator

$$ \mathcal{L}_{\kappa}(f) := -\partial_v \left[ f_{\kappa}(v) \partial_v \left( \frac{f}{f_{\kappa}} \right) \right], \quad \mathcal{L}_{\kappa} : \mathcal{D}(\mathcal{L}_{\kappa}) \subset L_{f_{\kappa}}^{2, -1} \to L_{f_{\kappa}}^{2, -1}, \qquad (2.11) $$

has no more a discrete spectrum, in contrast to the standard Fokker-Planck case (see Proposition 2). Nevertheless, one can now follow the same strategy as in the standard Fokker-Planck case, by searching for generalized eigenfunctions and a continuous spectral representation, as has been done in [49]. In the present paper however, we shall follow a different strategy, by considering the expansion of the $\kappa$-distribution function in a basis set, not constituted of (generalized) eigenfunctions of the corresponding Fokker-Planck operator $\mathcal{L}_\kappa$, but constituted of a different basis set of the functional Hilbert-space $L_{f_\kappa}^{2, -1}$, namely of rational Chebyshev polynomials introduced in the next Section. For the moment, let us finish this Section, by recalling a time-decay result of the solution to (2.10) towards its stationary solution. Due to the fact that there is no spectral gap, the time-decay rate will be no more exponential as in the standard Fokker-Planck case, but algebraic.

**Theorem 1.** [10, 49] (Time decay) Let $f$ be a solution of the evolution problem (2.10), with an initial condition $f_{in} f_\kappa^{-1} \in L^\infty(\mathbb{R}_v)$. Then for all $0 < p < 2\kappa - 1$, with $\kappa > 1$, the
following estimate holds

$$ \|f(t) - f_{\infty}\|_{f_{\kappa}^{-1}}^2 \le \left[ \|f_{in} - \langle f_{in} \rangle f_{\kappa}\|_{f_{\kappa}^{-1}}^{-4/p} + K_{p,\kappa} \frac{4t}{p} \right]^{-p/2}, \quad \forall t \ge 0, \quad (2.12) $$

with $f_{\infty} = \langle f \rangle f_{\kappa} := \left(\int_{\mathbb{R}_{\nu}} f \, dv\right) f_{\kappa} = \langle f_{in} \rangle f_{\kappa}$ and $K_{p,\kappa} > 0$ a constant depending on $p$ and $\kappa$.

2.3. Regularized energetic Fokker-Planck operator for $f_{eq} = f_{\kappa} e^{-a v^2/2}$. Finally, let us come now to the investigation of the regularized energetic Fokker-Planck equation, which writes

$$ (FP)_{\kappa,a} \left\{ \begin{array}{l} \partial_t f = \partial_v \left[ f_{\kappa,a}(v) \partial_v \left( \frac{f}{f_{\kappa,a}} \right) \right] = \partial_v \left[ \left( \frac{1}{1+\frac{v^2}{2\kappa}} + a \right) v f + \partial_v f \right], \quad \forall (t,v) \in \mathbb{R}^+ \times \mathbb{R}, \\ f(0, \cdot) = f_{in}(\cdot), \end{array} \right. \quad (2.13) $$

where the equilibrium is this time $f_{\kappa,a}(v) = c_{\kappa} (1+\frac{v^2}{2\kappa})^{-\kappa} e^{-a \frac{v^2}{2}}$, with $c_{\kappa} > 0$ the normalization constant of $f_{\kappa}$. We shall work now in the space $L_{f_{\kappa,a}^{-1}}^2$ and shall take in this study some fixed $\kappa \in \mathbb{R}^+$ and $0 < a \ll 1$.

To study the spectrum of this Fokker-Planck equation, the Liouville transformation $f = g \sqrt{f_{\kappa,a}}$ is very useful, permitting to switch from the eigenvalue problem $\mathcal{L}_{\kappa,a}(f) = \lambda f$ to the Schrödinger problem (sharing the same spectrum with the corresponding Fokker-Planck operator), i.e.

$$ -\partial_{vv}g + Q_a(v)g = \lambda g, $$

where the occurring potential $Q_a$ is a confining potential given by

$$ Q_a(v) := \frac{v^2}{4} \left[ \frac{1}{1 + \frac{v^2}{2\kappa}} + a \right]^2 - \frac{1}{2} \left[ \frac{1 - \frac{v^2}{2\kappa}}{(1 + \frac{v^2}{2\kappa})^2} + a \right] \sim_{|v| \gg 1} \frac{1}{4} a^2 v^2 - \frac{a}{2}. $$

Spectral theorems of quantum mechanics [57] permit now to show that the spectrum of this Schrödinger equation is discrete for $a > 0$ (it is continuous for $a = 0$), and this essentially due to confining properties of the potential, i.e. $Q_a(v) \rightarrow_{|v| \rightarrow \infty} \infty$. We recall that for $a = 0$ one has $Q_0(v) \sim_{|v| \gg 1} \frac{\kappa(\kappa+1)}{v^2}$ [49], hence not confining. The main goal of the regularization is thus to have a discrete spectrum. Furthermore, the Fokker-Planck operator $\mathcal{L}_{\kappa,a}$ being self-adjoint and positive, one has $\sigma(\mathcal{L}_{\kappa,a}) = \{\lambda_k\}_{k \in \mathbb{N}} \subset [0, \infty)$, with $\lambda_0 = 0$ the first eigenvalue, corresponding to the eigenfunction $f_{\kappa,a}$ of $\mathcal{L}_{\kappa,a}$.

To investigate now the decay in time of the solution to (2.13) towards the stationary state $f_{\infty} = \frac{\langle f_{in} \rangle}{\langle f_{\kappa,a} \rangle} f_{\kappa,a}$, it is more suitable to filter out the equilibrium by setting $f = h f_{\kappa,a}$ and to consider the so-called *Ornstein-Uhlenbeck* evolution equation for the unknown $h$ in the weighted Hilbert-space $L_{f_{\kappa,a}}^2$

$$ (OU)_{\kappa,a} \left\{ \begin{array}{l} \partial_t h = \frac{1}{f_{\kappa,a}} \partial_v [f_{\kappa,a}(v) \partial_v h], \quad \forall (t,v) \in \mathbb{R}^+ \times \mathbb{R}, \\ h(0, \cdot) = h_{in}(\cdot) = f_{in}(\cdot) f_{\kappa,a}^{-1}. \end{array} \right. \quad (2.14) $$

This equation has the same spectrum as the original Fokker-Planck equation $(FP)_{\kappa,a}$, exhibiting thus a spectral gap which leads to a Poincaré inequality of the type

$$ \|h - \bar{h}\|_{L_{f_{\kappa,a}}^2}^2 \le C_P \|\partial_v h\|_{L_{f_{\kappa,a}}^2}^2, \quad \forall h \in H_{f_{\kappa,a}}^1, $$
with $C_P > 0$ the Poincaré constant and the average defined as $\bar{h} := \int_{\mathbb{R}} h f_{\kappa,a} dv$. Hence, one has the following exponential time-decay property.

**Theorem 2. (Time decay)** Let $f$ be a solution of the evolution problem (2.13), with an initial condition $f_{in} \in L^2_{f_{\kappa,a}}$. Then the following estimate holds, for some constant $C > 0$

$$ \|f(t) - f_{\infty}\|_{f_{\kappa,a}^{-1}} \le C \|f_{in} - \frac{\langle f_{in} \rangle}{\langle f_{\kappa,a} \rangle} f_{\kappa,a}\|_{f_{\kappa,a}^{-1}} e^{-\lambda_{1,a} t}, \quad \forall t \ge 0, \quad (2.15) $$

where the stationary solution is given by $f_{\infty} = \frac{\langle f_{in} \rangle}{\langle f_{\kappa,a} \rangle} f_{\kappa,a}$ and $\lambda_{1,a} > 0$ is the first positive eigenvalue of $\mathcal{L}_{\kappa,a}$, satisfying $\lambda_{1,a} \xrightarrow{a \to 0} 0$.

### 3. SPECTRAL SCHEME BASED ON RATIONAL CHEBYSHEV FUNCTIONS

We are coming now to the design of an adequate spectral method for our energetic particle Fokker-Planck equation (2.10), we recall here for clarity reasons

$$ (FP)_{\kappa} \left\{ \begin{array}{l} \partial_t f = \partial_v \left[ f_{\kappa}(v) \partial_v \left( \frac{f}{f_{\kappa}(v)} \right) \right] = \partial_v \left[ \frac{v}{1 + \frac{v^2}{2\kappa}} f + \partial_v f \right], \quad \forall (t,v) \in \mathbb{R}^{+} \times \mathbb{R}, \\ f(0, \cdot) = f_{in}(\cdot). \end{array} \right. \quad (3.1) $$

In the previous section we evoked the Hermite spectral method for the resolution of the standard Fokker-Planck equation (2.4). Hermite functions however are exponentially decaying for large $v$, and hence their use is not natural for the approximation of distribution functions which are slowly decaying for large $v$, such as the equilibrium of our problem (3.1) described by $\kappa$-distributions. New basis sets are thus required and have to be introduced. We shall not focus on basis functions composed of generalized eigenfunctions (continuous spectral representation) as has been done in [49]. Using instead transformations that map a finite domain into an infinite domain, allows one to generate a broad class of new basis functions for the unbounded domain, that are images of, for example, Chebyshev polynomials. We shall focus on this strategy in this section.

**3.1. Rational Chebyshev polynomials.** At the ground of our subsequent polynomial basis construction are the Chebyshev polynomials, which are nothing else than a change-of-variable in the Fourier cosine functions. Indeed, the map $s \in [0, \pi] \mapsto \xi = \cos(s) \in [-1, 1]$ transforms the cosine functions into the Chebyshev polynomials $T_n(\xi)$ given by

$$ T_n(\xi) := \cos(n s) \quad \text{where} \quad s = \arccos(\xi), \quad \forall \xi \in [-1, 1]. \quad (3.2) $$

The thus defined set $\{T_n\}_{n \in \mathbb{N}}$ forms an orthogonal basis of the space $L^2((-1, 1); \rho(\xi)d\xi)$ with the weight $\rho(\xi) := \frac{1}{\sqrt{1-\xi^2}}$, in particular one has the orthogonality relations

$$ \int_{-1}^{1} \frac{T_n(\xi) T_m(\xi)}{\sqrt{1 - \xi^2}} d\xi = \begin{cases} 0, & \text{if } m \neq n \\ \pi, & \text{if } m = n = 0, \\ \frac{\pi}{2}, & \text{if } m = n \neq 0, \end{cases} \quad \forall n, m \in \mathbb{N}. $$

At this point, using now the algebraic map

$$ \xi \in (-1, 1) \mapsto v := \frac{L\xi}{\sqrt{1 - \xi^2}} \in (-\infty, \infty), \quad \text{or equiv.} \quad \xi := \frac{v}{\sqrt{L^2 + v^2}}, \quad (3.3) $$
respectively

$$s \in (0, \pi) \mapsto v := L \cot(s) \in (-\infty, \infty), \quad \text{or equiv. } s = \arccot(v/L), \qquad (3.4)$$

permits to transform the finite interval into an infinite interval, and one defines the rational Chebyshev polynomials as

$$C_n^L(v) := T_n(\xi) = \cos(n s) \quad \text{where} \quad s \in (0, \pi), \xi \in (-1, 1), v \in (-\infty, \infty). \qquad (3.5)$$

One can equally start from the sinus functions in order to define

$$S_n^L(v) := S_n(\xi) = \sin((n+1)s) \quad \text{where} \quad s \in (0, \pi), \xi \in (-1, 1), v \in (-\infty, \infty). \qquad (3.6)$$

Let us remark that $L > 0$ is a mapping parameter, which has to be chosen adequately.

The two sets $\{C_n^L\}_{n \in \mathbb{N}}$ resp. $\{S_n^L\}_{n \in \mathbb{N}}$ form two independent orthogonal basis sets of the Hilbert-space $L_{\sigma_L}^2$, with weight $\sigma_L(v) := \frac{L}{L^2+v^2}$. In particular one has the orthogonality relations

$$\int_{-\infty}^{\infty} C_n^L(v) C_m^L(v) \frac{L}{L^2 + v^2} dv = \begin{cases} 0, & \text{if } m \neq n \\ \pi, & \text{if } m = n = 0, \\ \frac{\pi}{2}, & \text{if } m = n \neq 0, \end{cases} \quad \forall n, m \in \mathbb{N},$$

and equally for $\{S_n^L\}_{n \in \mathbb{N}}$. For an efficient computation of these rational Chebyshev basis functions, one has at hand the following recurrence formulae

$$\frac{2v/L}{\sqrt{1+(v/L)^2}} C_n^L(v) = C_{n+1}^L(v) + C_{n-1}^L(v), \quad (C_n^L)'(v) = \frac{n}{2L} \frac{1}{\sqrt{1+(v/L)^2}} [C_{n-1}^L - C_{n+1}^L],$$

as well as

$$\frac{4(v/L)^2}{1+(v/L)^2} C_n^L(v) = C_{n+2}^L + 2C_n^L + C_{n-2}^L, \quad \frac{4}{1+(v/L)^2} C_n^L(v) = -C_{n+2}^L + 2C_n^L - C_{n-2}^L.$$

Using all these formulae permits to compute the first polynomials, which are given via

$$C_0^L \equiv 1, \quad C_1^L(v) = \frac{v/L}{\sqrt{1+(v/L)^2}}, \quad C_2^L(v) = \frac{(v/L)^2 - 1}{1+(v/L)^2}, \quad C_3^L(v) = \frac{(v/L)[(v/L)^2 - 3]}{[1+(v/L)^2]^{3/2}},$$

and

$$S_0^L(v) = \frac{1}{\sqrt{1+(v/L)^2}}, \quad S_1^L(v) = \frac{2(v/L)}{1+(v/L)^2}, \quad S_2^L(v) = \frac{3(v/L)^2 - 1}{[1+(v/L)^2]^{3/2}}.$$

To give an example, let us consider the simple rational (algebraically decaying) functions $[1+(v/L)^2]^{-1}$ resp. $[1+(v/L)^2]^{-3/2}$ which are rather harmless, however approximating them with Hermite polynomials fails drastically. In this case rational Chebyshev polynomials are more powerful, as can be observed through the relations

$$\frac{1}{1+(v/L)^2} = \frac{1}{2} [C_0^L(v) - C_2^L(v)], \quad \frac{1}{[1+(v/L)^2]^{3/2}} = \frac{1}{4} [3S_0^L(v) - S_2^L(v)].$$

More generally, one has for all $n \in \mathbb{N}$ the finite expansions

$$\frac{1}{[1+(v/L)^2]^n} = \alpha_0 C_0^L + \alpha_2 C_2^L + \cdots + \alpha_{2n} C_{2n}^L, \quad \frac{1}{[1+(v/L)^2]^{n+1/2}} = \beta_0 S_0^L + \beta_2 S_2^L + \cdots + \beta_{2n} S_{2n}^L. \qquad (3.7)$$

The next theorem describes the way of convergence of the rational Chebyshev expansion series.
**Theorem 3.** [9] (Convergence of the rational Chebyshev expansion)

Let $g \in L_{\sigma_L}^2$ be a function free of singularities on $\mathbb{R}$. Then one has the following properties:

(i) If $g$ has an asymptotic power series for $|v| \gg 1$ that contains only even, non-negative powers of $1/v$, the expansion series of $g$ in the basis $\{C_n^L(v)\}_{n \in \mathbb{N}}$ is exponentially convergent.

(ii) If $g$ has an asymptotic power series for $|v| \gg 1$ that contains only odd, positive powers of $1/v$, the expansion series of $g$ in the basis $\{S_n^L(v)\}_{n \in \mathbb{N}}$ is exponentially convergent.

(iii) If $g$ has an asymptotic power series for $|v| \gg 1$ that contains only integral, non-negative powers of $1/v$, the expansion series of $g$ in the full basis $\{C_n^L(v), S_n^L(v)\}_{n \in \mathbb{N}}$ is exponentially convergent.

(iv) If $g$ has an asymptotic power series for $|v| \gg 1$ that contains non-integral powers of $1/v$, the expansion series of $g$ in one of the basis $\{C_n^L(v)\}_{n \in \mathbb{N}}$ or/and $\{S_n^L(v)\}_{n \in \mathbb{N}}$ is only algebraically convergent.

**Remark 1.** Algebraic respectively exponential convergence for an expansion

$$f(v) = \sum_{n=0}^{\infty} \alpha_n \varphi_n(v), \quad \forall v \in \mathbb{R},$$

means there exist some constants $p > 0$, $\mu > 0$ and $\iota$, $N_* \in \mathbb{N}$, such that one of the two items is satisfied

* $|\alpha_n| \le p n^{-\iota}$, $\forall n \ge N_*$ (alg. conv.),

* $|\alpha_n| \le p e^{-\mu n}$, $\forall n \ge N_*$ (exp. conv.).

3.2. **RC-spectral method.** We shall set $L = \sqrt{2\kappa}$ and take $\kappa \in \mathbb{N}$ in the following. The question now is to find a basis set for our weighted Hilbert-space $L_{f_\kappa}^{2, -1}$, based on the just introduced rational Chebyshev basis sets of the space $L_{\sigma_\kappa}^2$, where we recall that

$$f_{\kappa}(v) = \frac{c_{\kappa}}{(1 + \frac{v^2}{2\kappa})^{\kappa}}, \quad c_{\kappa} := \frac{1}{\sqrt{2\pi\kappa}} \frac{\Gamma(\kappa)}{\Gamma(\kappa - 1/2)}, \quad \sigma_{\kappa}(v) := \frac{1}{\sqrt{2\kappa}} \frac{1}{1 + \frac{v^2}{2\kappa}}, \quad \forall \kappa \in \mathbb{N}.$$

Scaling the rational Chebyshev polynomials $\{C_n^L, S_n^L\}_{n \in \mathbb{N}}$ as

$$\Theta_n(v) := C_n^L(v) \Upsilon_\kappa(v), \quad \Xi_n(v) := S_n^L(v) \Upsilon_\kappa(v),$$

with

$$\Upsilon_{\kappa}(v) := [\sigma_{\kappa}(v) f_{\kappa}(v)]^{1/2} = \frac{(c_{\kappa}/\sqrt{2\kappa})^{1/2}}{(1 + \frac{v^2}{2\kappa})^{\frac{\kappa+1}{2}}}, \quad \forall v \in \mathbb{R},$$

permits to introduce two independent orthogonal basis sets $\{\Theta_n, \Xi_n\}_{n \in \mathbb{N}}$ for the functional space we are considering, namely $L_{f_\kappa}^2$. The recurrence formulae remain the same as above, the only modification is in the derivative, which becomes

$$\Theta'_{n}(v) = \frac{n}{2} \frac{1}{\sqrt{2\kappa}} \frac{1}{\sqrt{1 + \frac{v^2}{2\kappa}}} [\Theta_{n-1}(v) - \Theta_{n+1}(v)] - \frac{\kappa+1}{2\kappa} \frac{v}{1 + \frac{v^2}{2\kappa}} \Theta_{n}(v), \quad \forall v \in \mathbb{R}.$$

Taking an odd $\kappa \in \mathbb{N}$ prescribes an expansion in $\{\Theta_n\}_{n \in \mathbb{N}}$ (see Thm. 3 (i)) and symmetry properties permit furthermore to show that even functions $f(t, |v|)$ possess an expansion using solely even rational basis functions. Altogether we shall thus search for a solution to (3.1) under the form

$$f(t, v) = \sum_{n=0}^{\infty} \alpha_{2n}(t) \Theta_{2n}(v), \quad \forall (t, v) \in \mathbb{R}^{+} \times \mathbb{R}.$$
Let us now introduce this spectral expansion in the Fokker-Planck equation (3.1), in order to find an ODE system for the spectral coefficients {$\alpha_{2n}(t)$}$_{n\in\mathbb{N}}$. Straightforward computations yield firstly that

$$
\begin{align}
-\mathcal{L}_{\kappa}(\Theta_0) &= c_0 \Theta_0 + d_0 \Theta_2 + e_0 \Theta_4, \nonumber \\
-\mathcal{L}_{\kappa}(\Theta_2) &= b_2 \Theta_0 + c_2 \Theta_2 + d_2 \Theta_4 + e_2 \Theta_6, \tag{3.8} \\
-\mathcal{L}_{\kappa}(\Theta_{2n}) &= a_{2n} \Theta_{2n-4} + b_{2n} \Theta_{2n-2} + c_{2n} \Theta_{2n} + d_{2n} \Theta_{2n+2} + e_{2n} \Theta_{2n+4}, \quad \forall n \ge 2, \nonumber
\end{align}
$$

with

$$
c_0 = -\frac{(\kappa - 1)^2}{16\kappa}, \quad d_0 = -\frac{\kappa - 1}{4\kappa}, \quad e_0 = \frac{(\kappa - 1)(\kappa + 3)}{16\kappa},
$$

$$
b_2 = -\frac{\kappa - 1}{8\kappa}, \quad c_2 = -\frac{25 + \kappa(\kappa - 6)}{32\kappa}, \quad d_2 = -\frac{\kappa - 9}{8\kappa}, \quad e_2 = \frac{(\kappa - 3)(\kappa + 5)}{32\kappa},
$$

and

$$
a_n = \frac{(\kappa - 1 + n)(\kappa + 3 - n)}{32\kappa}, \quad b_n = -\frac{\kappa - (n-1)^2}{8\kappa}, \quad c_n = \frac{-3n^2 - (\kappa - 1)^2}{16\kappa}
$$

$$
d_n = -\frac{\kappa - (n+1)^2}{8\kappa}, \quad e_n = \frac{(\kappa - 1 - n)(\kappa + 3 + n)}{32\kappa}, \quad \forall n \ge 4.
$$

Secondly, taking now the weighted scalar-product (in $L_{f\kappa}^{2}$) of the Fokker-Planck equation
$\partial_t f = -\mathcal{L}_\kappa f$ with $\Theta_{2j}$, yields

$$
\left\{
\begin{aligned}
\alpha'_{0}(t) &= c_{0}\alpha_{0}(t) + b_{2}\alpha_{2}(t) + a_{4}\alpha_{4}(t), \\
&\vdots \\
\alpha'_{2j}(t) &= e_{2j-4}\alpha_{2j-4}(t) + d_{2j-2}\alpha_{2j-2}(t) + c_{2j}\alpha_{2j}(t) + b_{2j+2}\alpha_{2j+2}(t) + a_{2j+4}\alpha_{2j+4}(t).
\end{aligned}
\right.
$$

This leads to the following penta-diagonal ODE system to be solved for the computation
of the spectral coefficients $\mathcal{X}(t) := (\alpha_0(t), \alpha_2(t), \alpha_4(t), \dots)^t$, i.e.

$$
\mathcal{X}'(t) = M \mathcal{X}(t), \tag{3.9}
$$

with

$$
M = \begin{pmatrix}
c_0 & b_2 & a_4 & 0 & \cdots & 0 & 0 & 0 \\
d_0 & c_2 & b_4 & a_6 & \cdots & 0 & 0 & 0 \\
e_0 & d_2 & c_4 & b_6 & \cdots & 0 & 0 & 0 \\
\vdots & \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\
0 & 0 & 0 & \cdots & d_{2N-4} & c_{2N-2} & b_{2N} & 0 \\
0 & 0 & 0 & \cdots & e_{2N-4} & d_{2N-2} & c_{2N} & 0
\end{pmatrix} .
$$

**Remark 2.** The just introduced spectral method is very efficient for physical $\kappa$-values arising in astrophysics or fusion plasmas, such as $\kappa \in (1, 10)$, however it turns out to be problematic when $\kappa$ becomes larger and larger, preventing the investigation of the limit $\kappa \to \infty$. Indeed, one can immediately observe that some entries of the matrix $M$ become infinite for $\kappa \to \infty$ and the rational basis functions $\{C_n^L\}_{n \in \mathbb{N}}$ get also more and more similar in this limit,
rendering the RC-spectral method not adequate to treat cases with $\kappa \gg 1$. Furthermore, the expansion of the asymptotic equilibrium $f_\kappa(v)$ requires the use of $\kappa$ Chebyshev basis terms, as

$$f_{\kappa}(v) = \sum_{n=0}^{N} \alpha_{2n}(t) C_{2n}(v) \Upsilon_{\kappa}(v) \Leftrightarrow f_{\kappa}(v) \Upsilon_{\kappa}^{-1}(v) = \frac{\sqrt{c_{\kappa}} (2\kappa)^{1/4}}{(1 + \frac{v^2}{2\kappa})^{\frac{\kappa-1}{2}}} = \sum_{n=0}^{N} \alpha_{2n}(t) C_{2n}(v),$$

such that with (3.7) one notices that this is possible for $N = \kappa - 1$. Thus, with $\kappa \to \infty$ the expansion requires more and more terms, leading to an inefficient Chebyshev-spectral method. There are however some advantages for small $\kappa$ values, we shall discuss in the numerical section, such that we find that this scheme has its right to exist.

### 4. SPECTRAL SCHEME BASED ON NON-STANDARD ORTHOGONAL POLYNOMIALS

The aim of this section will be thus to construct a different spectral scheme, which shall degenerate into the Hermite spectral method of section 2.1, when $\kappa$ tends towards infinity. The idea of this new scheme is based on the Gram-Schmidt orthogonalisation procedure, which requires the existence of the moments of the weight function ($\kappa$-distribution), such that we shall be forced to consider in this section rather the regularized Fokker-Planck equation (2.13), i.e.

$$(FP)_{\kappa,a} \left\{ \begin{array}{ll} \partial_t f = -\mathcal{L}_{\kappa,a}(f), & \forall (t,v) \in \mathbb{R}^+ \times \mathbb{R}, \\ f(0, \cdot) = f_{in}(\cdot), & \end{array} \right. \qquad (4.1)$$

with the collision operator $\mathcal{L}_{\kappa,a} : \mathcal{D}(\mathcal{L}_{\kappa,a}) \subset L_{f_{\kappa,a}}^2 \to L_{f_{\kappa,a}}^2$ given by

$$\mathcal{L}_{\kappa,a}(f) := -\partial_v \left[ \left( \frac{1}{1 + \frac{v^2}{2\kappa}} + a \right) v f + \partial_v f \right] = -\partial_v \left[ f_{\kappa,a} \partial_v (f/f_{\kappa,a}) \right]. \qquad (4.2)$$

Even if this operator has a discrete spectrum as in the standard Fokker-Planck case, the eigenfunctions are non-standard and have no analytical expression, such that one has to search for them numerically. In this section however we shall follow a distinct strategy, and shall construct a different orthogonal basis set $\{\zeta_k(v)\}_{k \in \mathbb{N}}$ of $L_{f_{\kappa,a}}^2$, based on the Gram-Schmidt orthogonalization process, enabling then the expansion of the solution to (4.1) as

$$f(t, v) = \sum_{k=0}^{\infty} \alpha_k(t) \zeta_k(v), \quad \forall (t, v) \in \mathbb{R}^{+} \times \mathbb{R}, \qquad (4.3)$$

with $\{\alpha_k(t)\}_{k \in \mathbb{N}}$ to be determined.

To simplify the presentation, let us filter out the equilibrium from the solution to (4.1) by setting $f = h f_{\kappa,a}$, where the unknown $h \in L_{f_{\kappa,a}}^2$ is this time solution of the Ornstein-Uhlenbeck equation

$$(OU)_{\kappa,a} \left\{ \begin{array}{ll} \partial_t h = \frac{1}{f_{\kappa,a}} \partial_v [f_{\kappa,a} \partial_v h], & \forall (t,v) \in \mathbb{R}^+ \times \mathbb{R}, \\ h(0, \cdot) = h_{in}(\cdot) = f_{in} f_{\kappa,a}^{-1}. & \end{array} \right. \qquad (4.4)$$

We shall construct now an orthogonal polynomial basis set $\{p_k(v)\}_{k \in \mathbb{N}}$ of $L_{f_{\kappa,a}}^2$, considering $\omega(v) := f_{\kappa,a}(v)$ as weight function. The regularization is crucial here, as this construction requires the existence of all moments of the weight function. The Gram-Schmidt orthogonalization process yields the following recursive formula for the construction of a monic
polynomial set

$$p_{-1} \equiv 0, \quad p_0 \equiv 1, \quad p_1(v) = v, \qquad (4.5)$$
$$p_{k+1}(v) = vp_k(v) - \beta_k p_{k-1}(v), \quad \forall k \ge 1,$$

where

$$\beta_k := \frac{(v p_k, p_{k-1})_\omega}{\|p_{k-1}\|_{\omega}^2} = \frac{\|p_k\|_{\omega}^2}{\|p_{k-1}\|_{\omega}^2}, \quad \forall k \ge 1.$$

Classical weights, such as Maxwellian ones, give rise to classical families of orthogonal polynomials (Hermite polynomials for ex.). For such families explicit formulae are available for the recursion coefficients. Outside the classical framework however, numerical strategies are required for the computation of these coefficients, and the difficult part in the construction of the polynomial basis set is the efficient (accurate and stable) computation of these coefficients {$\beta_k$}$_{k \in \mathbb{N}}$. One has to take into account that we deal with polynomials on unbounded domains, and small (round-off or truncation) errors are immediately visible. In particular, polynomial spectral methods are more sensitive to errors than for ex. Fourier spectral methods. We shall present an efficient computational strategy for the recursion coefficients below, for the moment however, let us remark that the first coefficients write

$$\beta_1 = \frac{m_2}{m_0}, \quad \beta_2 = \frac{m_0 m_4 - m_2^2}{m_0 m_2}, \qquad (4.6)$$

where {$m_l$}$_{l \in \mathbb{N}}$ are the moments of order $l$ of the weight function $\omega(v)$, i.e.

$$m_l := \int_{\mathbb{R}} v^l \omega(v) dv, \quad \forall l \in \mathbb{N},$$

computed via Whittaker formula (1.12). Supposing now that we have constructed in this manner an orthogonal basis set {$p_k$}$_{k \in \mathbb{N}}$ of $L_{f_{\kappa,a}}^2$, the solution to (4.4) writes

$$h(t, v) = \sum_{k=0}^{\infty} \alpha_k(t) p_k(v), \quad \forall (t, v) \in \mathbb{R}^{+} \times \mathbb{R}. \qquad (4.7)$$

Inserting this decomposition in the Ornstein-Uhlenbeck evolution equation (4.4) and taking the $L_\omega^2$ scalar product with $p_l(v)$, permits to get a system for the computation of the spectral coefficients {$\alpha_k(t)$}$_{k \in \mathbb{N}}$. Indeed, one gets

$$\sum_{k=0}^{\infty} \alpha'_{k}(t) \int_{\mathbb{R}} p_{k}(v) p_{l}(v) \omega(v) dv = - \sum_{k=0}^{\infty} \alpha_{k}(t) \int_{\mathbb{R}} p'_{k}(v) p'_{l}(v) \omega(v) dv, \quad \forall l \in \mathbb{N}.$$

Remarkning that $(p_k, p_l)_\omega = \gamma_l \delta_{k,l}$, where we denote the weighted norm by $\gamma_l := ||p_l||_\omega^2 = \int_{\mathbb{R}} |p_l(v)|^2 \omega(v) dv$, one obtains the linear ODE system

$$\begin{pmatrix} \alpha_0'(t) \gamma_0 \\ \alpha_1'(t) \gamma_1 \\ \vdots \\ \alpha_l'(t) \gamma_l \\ \vdots \end{pmatrix} = - \begin{pmatrix} \theta_{00} & \cdots & \theta_{0k} & \cdots \\ \vdots & \ddots & \vdots & \ddots \\ \theta_{l0} & \cdots & \theta_{lk} & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix} \begin{pmatrix} \alpha_0(t) \\ \alpha_1(t) \\ \vdots \\ \alpha_k(t) \\ \vdots \end{pmatrix}, \qquad (4.8)$$

where the entries of the matrix are

$$\theta_{kl} = \theta_{lk} := \int_{\mathbb{R}} p'_{k}(v) p'_{l}(v) \omega(v) dv. \qquad (4.9)$$
It remains now to find an efficient manner to compute the integrals {$\theta_{lk}$}$_{l,k \in \mathbb{N}}$ as well as the recurrence coefficients {$\beta_k$}$_{k \in \mathbb{N}}$ and the weighted norms {$\gamma_k$}$_{k \in \mathbb{N}}$.

Several methods have been proposed in literature for an efficient computation of the re-
cursion coefficients, see for example [31]. The classical approach in terms of moments (see
(4.6)), computed via the Hankel determinants, is numerically problematic, being very ill-
conditioned. In this work we shall use the so called *Modified Chebyshev algorithm*, which is a
more stable and a better conditioned strategy. The idea is to work with "modified moments",
corresponding to a different measure (weight function) $d\mu := \mu(v) dv$, to be chosen close to
the original measure $dw := \omega(v) dv$, and with the additional requirement of possessing ex-
plicit recurrence coefficients. This shall permit to bypass the stability issues and generate in
a stable manner a non-classical orthogonal polynomial basis set.

To be more precise, taking for instance as weight $\mu(v) = \mathcal{M}(v) e^{-a \frac{v^2}{2}}$, which is close to our
original weight function $\omega_{\kappa}(v) = f_{\kappa} e^{-a \frac{v^2}{2}}$ in the sense that

$$
\lim_{\kappa \to \infty} \omega_{\kappa} = \mu, \tag{4.10}
$$

we construct the corresponding orthogonal monic polynomial basis set {$q_k$}$_{k \in \mathbb{N}}$ of $L_\mu^2(\mathbb{R})$, via
the three-term recurrence relation

$$
\begin{gather}
q_{-1} \equiv 0, \quad q_0 \equiv 1, \quad q_1(v) = v, \tag{4.11} \\
q_{k+1}(v) = vq_k(v) - b_k q_{k-1}(v), \quad k \ge 1. \nonumber
\end{gather}
$$

For such Maxwellian weights $\mu(v)$, the recursion coefficients $\{b_k\}_{k \in \mathbb{N}}$ are explicitly known and
$\{q_k\}_{k \in \mathbb{N}}$ are Hermite-like polynomials, in particular one has

$$
b_k = \frac{k}{1+a}, \quad \|q_k\|_{\mu}^{2} = \frac{k!}{(1+a)^{\frac{2k+1}{2}}}, \quad k \ge 1.
$$

Based on this new basis set, let us introduce the mixed moments

$$
\sigma_{k,l} = \int_{\mathbb{R}} p_k(v) q_l(v) \omega(v) dv, \quad \forall k, l \ge -1,
$$

which shall help to construct the original basis set {$p_k$}$_{k \in \mathbb{N}}$.

By orthogonality, we have $\sigma_{k,l} = 0$ for $k > l$. Moreover, as $p_k$ and $q_k$ are monic polynomials of degree $k$, we can write

$$
p_k(v) = q_k(v) + r(v), \quad \text{with } \deg(r) < k,
$$

such that using the orthogonality of {$p_k$}$_{k \in \mathbb{N}}$, one has

$$
\gamma_k = \|p_k\|_{\omega}^2 = \int_{\mathbb{R}} p_k(v) p_k(v) \omega(v) dv = \int_{\mathbb{R}} p_k(v) q_k(v) \omega(v) dv = \sigma_{k,k}. \quad (4.12)
$$

Thus

$$
\beta_k = \frac{\|p_k\|_{\omega}^2}{\|p_{k-1}\|_{\omega}^2} = \frac{\sigma_{k,k}}{\sigma_{k-1,k-1}}, \quad \forall k \ge 1. \quad (4.13)
$$

Using the recurrence formulae (4.5) for $p_k$ respectively (4.11) for $q_k$, we derive

$$
\sigma_{k,l} = \sigma_{k-1,l+1} - \beta_{k-1} \sigma_{k-2,l} + b_l \sigma_{k-1,l-1}, \quad \forall k, l \ge 1. \quad (4.14)
$$

Hence, starting from the values of $\sigma_{-1,l} = 0$ and $\{\sigma_{0,l}\}_{l \in \mathbb{N}} = \int_{\mathbb{R}} q_l(v) \omega(v) dv$, which are noth-
ing more than linear combinations of the moments $\{m_l\}_{l \in \mathbb{N}}$, computed via the Whittaker
formulae (1.12), one can compute step by step the values of $\{\sigma_{k,l}\}_{k,l \in \mathbb{N}}$, $\{\beta_k\}_{k \in \mathbb{N}}$ and $\{\gamma_k\}_{k \in \mathbb{N}}$
via (4.12)-(4.14).

For the resolution of our Fokker-Planck problem (4.4), thus of (4.8), one needs now to compute the entries of the matrix, in particular $\theta_{kl} := \int_{\mathbb{R}} p'_{k}(v) p'_{l}(v) \omega(v) dv$. This shall be done by introducing the following quantities

$$ \chi_{k,l} := \int_{\mathbb{R}} p_k(v) p'_l(v) \omega(v) dv, \quad \xi_{k,l} := \int_{\mathbb{R}} p'_k(v) p_l(v) \omega(v) dv. $$

Using now the recursion formulae (4.5)-(4.11) permits to obtain for all $k, l \in \mathbb{N}$

$$ \theta_{k,l} = \theta_{k-1,l+1} + \beta_l \theta_{k-1,l-1} + \chi_{k-1,l} - \beta_{k-1} \theta_{k-2,l} - \xi_{k-1,l}, $$

$$ \chi_{k,l} = \chi_{k-1,l+1} + \beta_l \chi_{k-1,l-1} - \gamma_l \delta_{k-1,l} - \beta_{k-1} \chi_{k-2,l}, $$

$$ \xi_{k,l} = \xi_{k-1,l+1} + \beta_l \xi_{k-1,l-1} + \gamma_l \delta_{k-1,l} - \beta_{k-1} \xi_{k-2,l}. $$

Thus starting from the values $\theta_{0,l} = 0$, $\chi_{-1,l} = 0$, $\xi_{0,l} = 0$ for all $l \in \mathbb{N}$, as well as from the values one has to compute as an initial step

$$ \theta_{1,l} := \int_{\mathbb{R}} p'_{l}(v) \omega(v) dv, \quad \chi_{0,l} := \int_{\mathbb{R}} p'_{l}(v) \omega(v) dv, \quad \xi_{1,l} := \int_{\mathbb{R}} p_{l}(v) \omega(v) dv, $$

one gets the desired values of $\{\theta_{k,l}\}_{k,l \in \mathbb{N}}$. The computation of $\theta_{1,l}$ shall be done via Whittaker's formulae (1.12), for accuracy reasons. Remark also that $\xi_{1,l} = \gamma_0 \delta_{0,l}$ due to the orthogonality of the polynomials $\{p_k\}_{k \in \mathbb{N}}$, and that $\chi_{k,l} = \xi_{l,k}$ for all $k, l \in \mathbb{N}$. Furthermore, as one has the decomposition

$$ p'_{l}(v) = \sum_{k=0}^{l-1} c_{k}^{(l)} p_{k}(v), \quad \text{with coefficients } c_{k}^{(l)} = \frac{\xi_{l,k}}{\gamma_{k}}, $$

one has immediately that $\chi_{k,l} = 0$ for all $k > l - 1$ as well as $\xi_{k,l} = 0$ for all $l > k - 1$. Finally this decomposition permits also to observe that

$$ \chi_{0,l} = c_{0}^{(l)} \gamma_{0} = \theta_{1,l}, \quad c_{0}^{(l)} := (p'_{l}, p_{0})_{\omega}. $$

**Remark 3.** Let us mention that the Gram-Schmidt spectral method is very performant for large $\kappa$-values, $\kappa \gg 1$, as in this case the weight $\omega_\kappa$ is close to the Maxwellian weight $\mu$, whereas for small values of $\kappa$ the use of Whittaker's formulae (1.12) becomes problematic due to the fact, that we deal with polynomials on unbounded domains, rendering the utility of the GS-scheme very awkward. There is a delicate interplay between the regularization parameter $\alpha$ (thus somehow the truncation of the domain) and the $\kappa$-parameter (meaning the hot tails of the distribution function) one has to care of, in order to avoid instability problems. Let us also observe that the asymptotic equilibrium $f_{\kappa,a}$ requires only one term in the GS-expansion, which is a huge advantage of the scheme, permitting the use of very few basis terms. Finally, in the limit $\kappa \to \infty$, the GS-scheme turns into the Hermite-spectral scheme, and this due our specific construction, in particular see (4.5)-(4.10)-(4.11).

## 5. NUMERICAL SIMULATIONS

In this section we present the numerical simulations obtained with the spectral methods constructed in this work, namely the Chebyshev spectral method (RC-scheme) resp. the Gram-Schmidt spectral method (GS-scheme). The aim is to study the efficiency of these schemes for the resolution of the energetic Fokker-Planck equation, comparing them with some standard finite-difference methods (FD-scheme). We are particularly interested in the
accuracy of these methods, simulation times as well as the asymptotic convergence rates towards the equilibrium, for different values of the parameter $\kappa$ (smaller and larger ones).

### 5.1. Numerical simulations with a standard FD-method.
Let us start with presenting the test case and the numerical results obtained with a standard finite difference scheme. We shall be concerned in the following with the resolution of the spatially homogeneous Fokker-Planck equation

$$ (FP)_{\kappa,a} \left\{ \begin{array}{l} \partial_t f = \partial_v \left[ \left( \frac{v}{1 + \frac{v^2}{2\kappa}} + av \right) f + \partial_v f \right], \quad \forall (t,v) \in \mathbb{R}^+ \times \mathbb{R}, \\ f(0, \cdot) = f_{in}(\cdot), \end{array} \right. \quad (5.1) $$

where we shall set either $a = 0$ (for the FD- and RC-scheme) or $a = 10^{-3}$ (regularized version for the GS-scheme), play with different values of $\kappa \in (1, 50)$, and choose an initial condition of two-bump form, composed of two shifted $\kappa$-distribution functions

$$ f_{in}(v) := \frac{1}{2} [f_{\kappa}(v+u) + f_{\kappa}(v-u)] \quad \text{with} \quad u=2. \quad (5.2) $$

In the long-time limit $t \to \infty$ the solution to this problem tends towards the stationary solution, given by

$$ f_{\infty}(v) := n f_{\kappa,a}(v), \quad n := \frac{\langle f_{in} \rangle}{\langle f_{\kappa,a} \rangle}, $$

where we recall the time-decay Theorems 1 and 2 as well as the notation

$$ \langle g \rangle := \int_{\mathbb{R}} g(v) dv. $$

Let us remark that the initial condition enters the asymptotic limit function $f_{\infty}$ only via the mass-constant $n$, which has thus to be computed precisely. The aim is now to study the performance of both spectral schemes in describing the relaxation of the distribution function $f$ towards the equilibrium, when compared to a standard FD-scheme we shall present now.

Set $\kappa = 3$ in this section and denote for simplicity reasons

$$ \eta(v) := \frac{1}{1 + \frac{v^2}{2\kappa}} + a. $$

We shall now truncate the velocity domain as $[-v_{max}, v_{max}]$ with $v_{max}$ to be defined later on, impose homogeneous Dirichlet boundary conditions on $f$, and consider a (not necessarily uniform) velocity mesh $\{v_j\}_{j=0}^{N_v}$, defining the mid-points $v_{j+1/2} := (v_j + v_{j+1})/2$ on which equation (5.1) will be approximated via a finite volume scheme. Integrating (5.1) on $[v_{j-1/2}, v_{j+1/2}]$ and using the finite volume notation $f_j(t) := \frac{1}{\Delta v_j} \int_{v_{j-1/2}}^{v_{j+1/2}} f(t, v) dv$ with the cell
length $\Delta v_j := v_{j+1/2} - v_{j-1/2}$, leads to

$$
\begin{align}
\frac{d}{dt} f_j(t) &= \frac{1}{\Delta v_j} \left[ \left( \eta_{j+1/2} v_{j+1/2} f(t, v_{j+1/2}) + (\partial_v f)(t, v_{j+1/2}) \right) \right. \nonumber \\
&\qquad \left. - \left( \eta_{j-1/2} v_{j-1/2} f(t, v_{j-1/2}) + (\partial_v f)(t, v_{j-1/2}) \right) \right] \nonumber \\
&\approx \frac{1}{\Delta v_j} \left[ \left( \eta_{j+1/2} v_{j+1/2} \frac{f_j(t) + f_{j+1}(t)}{2} + \frac{f_{j+1}(t) - f_j(t)}{v_{j+1} - v_j} \right) \right. \nonumber \\
&\qquad \left. - \left( \eta_{j-1/2} v_{j-1/2} \frac{f_j(t) + f_{j-1}(t)}{2} + \frac{f_j(t) - f_{j-1}(t)}{v_j - v_{j-1}} \right) \right] \nonumber \\
&= \frac{1}{\Delta v_j} \left\{ f_{j-1}(t) \left( -\eta_{j-1/2} \frac{v_{j-1/2}}{2} + \frac{1}{v_j - v_{j-1}} \right) \right. \nonumber \\
&\qquad \left. + f_j(t) \left[ \left( \eta_{j+1/2} \frac{v_{j+1/2}}{2} - \frac{1}{v_{j+1} - v_j} \right) - \left( \eta_{j-1/2} \frac{v_{j-1/2}}{2} + \frac{1}{v_j - v_{j-1}} \right) \right] \right. \nonumber \\
&\qquad \left. + f_{j+1}(t) \left( \eta_{j+1/2} \frac{v_{j+1/2}}{2} + \frac{1}{v_{j+1} - v_j} \right) \right\}, \tag{5.3}
\end{align}
$$

for all $j = 1, \dots, N_v - 1$, whereas we set for the boundary values $f_0 = f_{N_v} \equiv 0$.

In what concerns the time discretization, an implicit Euler scheme is used to solve the thus obtained ODE system, leading to

$$ (\mathbb{I} - \Delta t A) f^{n+1} = f^n, $$

with A a tridiagonal matrix of size $(N_v - 1) \times (N_v - 1)$, given by

$$ A_{j,j-1} := \frac{1}{\Delta v_j} \left( - \eta_{j-1/2} \frac{v_{j-1/2}}{2} + \frac{1}{v_j - v_{j-1}} \right), $$

$$ A_{j,j} := \frac{1}{\Delta v_j} \left[ \left( \eta_{j+1/2} \frac{v_{j+1/2}}{2} - \frac{1}{v_{j+1} - v_j} \right) - \left( \eta_{j-1/2} \frac{v_{j-1/2}}{2} + \frac{1}{v_j - v_{j-1}} \right) \right], $$

$$ A_{j,j+1} := \frac{1}{\Delta v_j} \left( \eta_{j+1/2} \frac{v_{j+1/2}}{2} + \frac{1}{v_{j+1} - v_j} \right). $$

We plotted on the left of Figure 3 the initial condition $f_{in}(v)$ corresponding to the two-bump distribution function (5.2), as well as the time evolution of this initial velocity distribution function towards the asymptotic stationary state $f_\infty(v)$. In order to evaluate the convergence towards this equilibrium, we plotted then on Figure 4 the time evolution of the $L_{f_\kappa}^2$ error between $f(t, \cdot)$ and $f_\infty$. One remarks that after an initial algebraic (quasi exponential) decay, the error starts to increase or saturates, as observed from the right plot of Fig. 4. The slope of this $log-log$ error-curve for large $t \gg 1$ (before saturation) is approximately equal to $-5$. Examining the error estimate in Theorem 1, we have that $-p/4 = -5$, meaning $p = 20$, which indicates that the scheme is performing better than it should, after Theorem 1. However, let us underline here that it is very intricate to verify numerically the theoretical error estimate, due to several reasons. Firstly, the result of this theorem is not necessarily optimal, for example the explicit expression of the constant $K_{p,\kappa}$ is not known. Secondly, our numerical simulations are performed with a specific initial condition, and we went only to final times of the order of $T \sim 30$, which is not necessary enough for comparison with the long-time asymptotics of the theoretical error estimate. And finally, an essential argument is that we
truncated the domain in our numerical simulations, a fact that changes an algebraic decay into an exponential decay.

Furthermore, we also observe that the mass is not accurately conserved numerically (see right of Fig. 3). The reason for these shortcomings comes from the fact, that we truncated the velocity domain at $v_{max} = 15$ and set homogeneous Dirichlet boundary conditions $f(t, \pm v_{max}) = 0$ for all times $t \ge 0$. However, since we are dealing with energetic particles, thus distribution functions with heavy tails, truncating the domain necessarily leads to mass losses as well as to a distribution function that gets away after a while from the asymptotic state $f_\infty$. Truncating the domain more far away reduces this error, however increases naturally the computational time. It is because of such difficulties in the numerical discretization of the energetic particle dynamics, that we proposed in this paper the two spectral schemes, the RC- and the GS-scheme, which we shall confront in detail in the next subsection.

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Left: Time evolution of the velocity distribution function $f(t, v)$ starting from the two-bump initial condition (5.2) with $\kappa = 3$. Right: Plot of the mass $\int_{\mathbb{R}} f(t, v) dv$ for three different domain truncations $[-v_{max}, v_{max}]$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Error evolution in time, $\kappa = 3$. Left: $\|f(t) - f_\infty\|_{L^2_{f_\kappa^{-1}}}$ for $v_{max} = 30$. Right: $\|f(t) - f_\infty\|_{L^2_{f_\kappa^{-1}}}$ in log-log scale and for two different domain truncations $v_{max} = 15$ and $v_{max} = 30$.</figcaption></figure>
5.2. **Numerical results obtained with our spectral methods.** Let us investigate now the numerical results obtained for the resolution of the previously introduced problem (5.1)-(5.2) with both spectral methods, namely the RC-scheme presented in Section 3 as well as the GS-scheme designed in Section 4. The solution $f(t, v)$ of this problem is approximated by one of the two following truncated series (corresponding to the two different basis-sets, introduced previously)

$$f_{RC}^{(N)}(t, v) := \sum_{k=0}^{\tilde{N}_{RC}} \mathbf{a}_{2k}(t) \Theta_{2k}(v), \quad f_{GS}^{(N)}(t, v) := \sum_{k=0}^{N_{GS}} \alpha_k(t) p_k(v) f_{\kappa,a}(v),$$

with $\tilde{N}_{RC} = N_{RC}/2$, $\mathbf{a}_{2k}(t)$ resp. $\alpha_k(t)$ being solutions of the linear systems (3.9) for the RC-scheme resp. system (4.8) for the GS-scheme, solved via a first-order implicit Euler-scheme. A second-order Crank-Nicolson scheme could be used to improve the accuracy. Let us remark that in the asymptotic limit $t \to \infty$ the steady state needs for the RC-scheme $N_{RC} + 1 = \kappa$ terms (rightly only half of them), whereas the GS-scheme requires only one term, thus $N_{GS} = 0$. The parameters used for all the following simulations are summarized in Table 1. The initial condition we choose here is again the two-bump distribution defined in (5.2), its time-evolution being similar to the one plotted on the left of Fig. 3.

<table><caption>TABLE 1. Numerical parameters used for the confrontation of the different schemes.</caption><thead><tr><th></th><th>DF<sub>ref</sub>-scheme</th><th>DF-scheme</th><th>RC-scheme</th><th>GS-scheme</th></tr></thead><tbody><tr><td>a</td><td>0</td><td>0</td><td>0</td><td>10<sup>-3</sup></td></tr><tr><td>κ</td><td>3 or 31</td><td>3 or 31</td><td>3</td><td>31</td></tr><tr><td>N<sub>v</sub>; v<sub>max</sub></td><td>10 001; 1 000</td><td>[101 ... 1 001]; 15</td><td>Θ</td><td>Θ</td></tr><tr><td>N<sub>RC/GS</sub></td><td>Θ</td><td>Θ</td><td>8, 10, 16</td><td>6, 8, 10, 16</td></tr><tr><td>T; N<sub>t</sub></td><td>10; 10 000</td><td>10; [100 ... 4 000]</td><td>10; [100 ... 4 000]</td><td>10; [100 ... 4 000]</td></tr></tbody></table>

To start the procedure, we need to expand the initial condition in the corresponding orthogonal basis set, which is done via

$$\int_{\mathbb{R}} f_{in}(v) \Theta_{2l}(v) f_{\kappa}^{-1}(v) dv = \mathbf{a}_{2l}(0) \mathbf{c}_l, \quad \int_{\mathbb{R}} f_{in}(v) p_l(v) dv = \alpha_l(0) \gamma_l, \quad \forall l \ge 0,$$

where $\mathbf{c}_0 = \pi$, $\mathbf{c}_l = \pi/2$ for $l \ge 1$ and we recall that $\gamma_l$ is the weighted norm $\gamma_l = \|p_l\|_{f_{\kappa,a}}^2$. The computation of these integrals is an important step in our spectral methods, a bad (not-accurate or too time-consuming) computation of these projections would lead to a loss of all the advantages of the spectral schemes. Exact computations or Gaussian quadratures could lead to very precise starting points, we shall use in this paper however standard integration methods and we illustrated, for validation, in Figures 7 resp. 8 the reconstruction of the initial condition (with the corresponding basis set). From these Figures one can identify the number of basis-functions required in order to have a precise simulation in the initial time-phase, for example the RC-scheme requires $N_{RC} = 16$ for a precision of $1.6 \cdot 10^{-3}$ whereas the GS-scheme requires $N_{GS} = 12$ basis functions for a precision of $2.6 \cdot 10^{-3}$. The initial condition (and thus the reconstruction) is not any more important in the asymptotic regime, only the approximation of the mass enters in the steady state.

A few simulations are carried out now in the aim to compare the proposed schemes with respect to convergence, accuracy, computational efficiency and the applicability to the specific physical situation we are interested in. For these comparison studies a reference solution
$f_{ref}(t, v)$ of (5.1)-(5.2) has to be constructed, and this is done here via the standard FD-scheme on a very large domain $v_{max} = 1000$ and with a very fine mesh of $N_t = 10000$ and $N_v = 10001$. This reference solution can be considered somehow as an exact solution. The discrete $L_v^2$-norm we shall use in the following for the error investigations, is defined as

$$ ||f_{ref} - f_{num}||_{L_v^2(t_*)} := \frac{2v_{max}}{N_v} \sum_{j=0}^{N_v-1} |f_{ref}(t_*, v_j) - f_{num}(t_*, v_j)|^2. $$

Let us start by presenting in Figures 5 and 6 the curves corresponding to the reference distribution function $f_{ref}(t, v)$ as well as the three distributions $f_{num}(t, v)$ computed with the FD-, RC- resp. GS-schemes at two different times $t_* = 0.2$ and $t_* = 2$ as well as two different $\kappa$-values, namely $\kappa = 3$ and $\kappa = 31$. What can be observed from these plots is that for small times, the accurate projection of the initial condition on the corresponding basis-set is of great importance. In particular the RC-scheme has some problems for large $\kappa$ and this due to the necessity to take a large number of basis-functions for the reconstruction, namely $N_{RC} > \kappa - 1$ as explained in Remark 2, whereas the GS-scheme has some problems for small $\kappa$-values, due to the hot tails of the distribution functions and the need to take a large (non-physical) regularization parameter $a > 0$ in order to overcome stability problems, as explained in Remark 3. For larger times, all numerical solutions get closer and closer to the equilibrium solution $f_\infty(v)$, and the initial condition does no longer play a so important role. Each scheme is thus better adapted for a particular situation, meaning a particular $\kappa$-regime. We shall hence take in the following $\kappa = 3$ for the RC-scheme and $\kappa = 31$ for the GS-scheme. In order to find out how much basis-functions are required for each scheme, to have a good precision especially in the small-time regime, we plotted on Figures 7 resp. 8 the reconstruction of the initial condition with the RC-scheme resp. the GS-scheme.

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Left: Velocity distribution functions $f_{ref}(t_*, v)$ resp. $f_{num}(t_*, v)$ for $t_* = 0.2$ and $\kappa = 3$. Right: Same velocity distribution functions at instant $t_* = 2$. Furthermore $N_{RC} = N_{GS} = 10$.</figcaption></figure>

To demonstrate the convergence of our spectral schemes, we depicted on Figure 9 resp. 10 the $L_v^2$-errors of our spectral methods with respect to the reference solution, i.e. $||f_{ref}(t_*, \cdot) - f_{RC/GS}(t_*, \cdot)||_{L_v^2}$, errors as a function of the time-discretization step $\Delta t$ and this at three different instants $t_*$ and for three different truncation indices $N_{RC/GS}$. For the Chebyshev scheme we used $\kappa = 3$, while for the Gram-Schmidt scheme we chose $\kappa = 31$ in accordance with the remarks made above. What can be observed is the convergence of both schemes,
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Left: Velocity distribution functions $f_{ref}(t_{\star}, v)$ resp. $f_{num}(t_{\star}, v)$ for $t_{\star} = 0.2$ and $\kappa = 31$. Right: Same velocity distribution functions at instant $t_{\star} = 2$. Again $N_{RC} = N_{GS} = 10$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 7. Reconstruction of the initial condition $f_{in}(v)$ via the RC basis-set ($\kappa = 3$), for several truncation indices $N_{RC}$. Left: Velocity distribution function. Right: Error plots $||f_{recon}^{(N)} - f_{in}||_{L^2}$.
For $N = 16$, the error is of $1.6 \cdot 10^{-3}$.</figcaption></figure>

the errors being linear in the discretization step $\Delta t$, which is completely natural, as we are discretizing the system in time via the implicit Euler scheme. Note that the slope becomes very small in the last plot of both Figures 9 resp. 10. This is simply due to saturation effects, the reference solution being not an exact solution, but a numerical solution, where we fixed the time-step to $\Delta t = 10^{-3}$. This artifact is clearly visible in the last plot of Figure 10, where the error saturates at some value of $1.212 \cdot 10^{-3}$. As a consequence, our RC- and GS-schemes become at large times too accurate, to be any more compared with this reference solution. Furthermore, one observes that in the GS-scheme the number of basis-functions does not any more play a role, after some transition-phase in time, the asymptotic regime needing indeed only one term in the expansion.

Naturally, increasing the domain via $v_{max}$, the number of grid points $N_v$ or the truncation index $N_{RC/GS}$ shrinks the error, but also increases the simulation time. Let us investigate now whether our spectral schemes have an advantage with respect to standard difference methods. To compare the simulation times of the three different methods, we performed
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 8. Reconstruction of the initial condition $f_{in}(v)$ via the GS basis set ($\kappa = 31$), for several truncation indices $N_{GS}$. Left: Velocity distribution function. Right: Error plots $||f_{recon}^{(N)} - f_{in}||_{L^2_v}$. For $N = 12$, the error is of $2.6 \cdot 10^{-3}$.</figcaption></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 9. Error curves $||f_{ref}(t_{\star}, \cdot) - f_{RC}(t_{\star}, \cdot)||_{L^2_v}$ for $\kappa = 3$, three time instants $t_{\star} = 0.2$; 2; 10 and three different truncation indices $N = 8$; 10; 16.</figcaption></figure>

several simulations for three different sets ($S_1$, $S_2$, $S_3$) of parameters, summarized in Table 2, and three different instants $t_{\star} = 0.2$; 2; 10. We then collected in Tables 3-5 the obtained
<figure><img src="image_14.png" /><figcaption>FIGURE 10. Error curves $||f_{ref}(t_{\star}, \cdot) - f_{GS}(t_{\star}, \cdot)||_{L^2}$ for $\kappa = 31$, three time instants $t_{\star} = 0.2$; 2; 10 and three different truncation indices $N = 8$; 10; 16.</figcaption></figure>

errors (with respect to the reference solution) in mass, denoted by $\mathbb{E}_m$, the overall error $||f_{ref}(t_{\star}, \cdot) - f_{num}(t_{\star}, \cdot)||_{L^2}$ denoted simply by $\mathbb{E}_f$, as well as the computational times. One observes that for short times ($t_{\star} = 0.2$) the DF-scheme is clearly the most advantageous method for both $\kappa$-values. For intermediate times ($t_{\star} = 2$) the GS-scheme starts to be competitive with respect to the FD-scheme, the RC-scheme still being behind in efficiency. In long-time regimes ($t_{\star} = 10$ and larger times) the RC- and the GS-spectral methods clearly take the overhand.

<table><caption>TABLE 2. Three parameter sets for the simulations of Tables 3-5.</caption><thead><tr><th>Δt = 0.01</th><th>N<sub>v</sub></th><th>v<sub>max</sub></th><th>N<sub>RC</sub></th><th>N<sub>GS</sub></th></tr></thead><tbody><tr><td>S<sub>1</sub></td><td>100</td><td>15</td><td>8</td><td>6</td></tr><tr><td>S<sub>2</sub></td><td>500</td><td>30</td><td>10</td><td>8</td></tr><tr><td>S<sub>3</sub></td><td>1000</td><td>30</td><td>16</td><td>10</td></tr></tbody></table>

Let us finish our comparison studies, by investigating some physical resp. mathematical properties of our spectral solutions. Firstly, let us remark that both spectral schemes preserve the mass equally accurate for all sets of parameters, a property that comes from our mass-computation procedure, synthesized here. Integrating the spectral expansions over the
<table><caption>TABLE 3. Simulation times for the FD-scheme, RC- resp. GS-scheme, at instant $t_* = 0.2$. Three sets of parameters have been chosen, summarized in Table 2. The discretization step is fixed at $\Delta t = 0.01$.</caption><thead><tr><th>t = 0.2</th><th>ε<sub>m</sub></th><th>ε<sub>f</sub></th><th>time</th><th>ε<sub>m</sub></th><th>ε<sub>f</sub></th><th>time</th><th>ε<sub>m</sub></th><th>ε<sub>f</sub></th><th>time</th></tr></thead><tbody><tr><td>DF (κ = 3)</td><td>9.21 · 10<sup>-5</sup></td><td>1.23 · 10<sup>-3</sup></td><td>0.0035</td><td>2.11 · 10<sup>-6</sup></td><td>5.76 · 10<sup>-4</sup></td><td>0.043</td><td>2.11 · 10<sup>-6</sup></td><td>5.7 · 10<sup>-4</sup></td><td>0.5</td></tr><tr><td>RC (κ = 3)</td><td>4.75 · 10<sup>-5</sup></td><td>3.2 · 10<sup>-2</sup></td><td>0.041</td><td>4.75 · 10<sup>-5</sup></td><td>3.1 · 10<sup>-2</sup></td><td>0.1</td><td>4.75 · 10<sup>-5</sup></td><td>2.7 · 10<sup>-2</sup></td><td>1.82</td></tr><tr><td>DF (κ = 31)</td><td>7.7 · 10<sup>-16</sup></td><td>1.43 · 10<sup>-3</sup></td><td>0.0033</td><td>0.0</td><td>8.53 · 10<sup>-4</sup></td><td>0.039</td><td>0.0</td><td>9.87 · 10<sup>-4</sup></td><td>0.49</td></tr><tr><td>GS (κ = 31)</td><td>4.44 · 10<sup>-16</sup></td><td>8.6 · 10<sup>-3</sup></td><td>0.022</td><td>4.44 · 10<sup>-16</sup></td><td>7.39 · 10<sup>-3</sup></td><td>0.03</td><td>4.44 · 10<sup>-16</sup></td><td>6.46 · 10<sup>-3</sup></td><td>0.042</td></tr></tbody></table>

<table><caption>TABLE 4. Simulation times for the FD-scheme, RC- resp. GS-scheme, at instant $t_* = 2$. Three sets of parameters have been chosen, summarized in Table 2. The discretization step is fixed at $\Delta t = 0.01$.</caption><thead><tr><th>t = 2</th><th>ε<sub>m</sub></th><th>ε<sub>f</sub></th><th>time</th><th>ε<sub>m</sub></th><th>ε<sub>f</sub></th><th>time</th><th>ε<sub>m</sub></th><th>ε<sub>f</sub></th><th>time</th></tr></thead><tbody><tr><td>DF (κ = 3)</td><td>9.6 · 10<sup>-4</sup></td><td>5.59 · 10<sup>-4</sup></td><td>0.012</td><td>2.5 · 10<sup>-5</sup></td><td>7.97 · 10<sup>-4</sup></td><td>0.21</td><td>2.5 · 10<sup>-5</sup></td><td>9.77 · 10<sup>-4</sup></td><td>4.0</td></tr><tr><td>RC (κ = 3)</td><td>4.75 · 10<sup>-5</sup></td><td>3.94 · 10<sup>-2</sup></td><td>0.044</td><td>4.75 · 10<sup>-5</sup></td><td>3.68 · 10<sup>-2</sup></td><td>0.1</td><td>4.75 · 10<sup>-5</sup></td><td>3.76 · 10<sup>-2</sup></td><td>1.82</td></tr><tr><td>DF (κ = 31)</td><td>1.34 · 10<sup>-14</sup></td><td>1.01 · 10<sup>-3</sup></td><td>0.01</td><td>1.77 · 10<sup>-15</sup></td><td>1.24 · 10<sup>-3</sup></td><td>0.21</td><td>4.66 · 10<sup>-15</sup></td><td>1.53 · 10<sup>-3</sup></td><td>4.05</td></tr><tr><td>GS (κ = 31)</td><td>4.44 · 10<sup>-16</sup></td><td>5.29 · 10<sup>-3</sup></td><td>0.022</td><td>4.44 · 10<sup>-16</sup></td><td>5.3 · 10<sup>-3</sup></td><td>0.03</td><td>4.44 · 10<sup>-16</sup></td><td>5.3 · 10<sup>-3</sup></td><td>0.042</td></tr></tbody></table>

<table><caption>TABLE 5. Simulation times for the FD-scheme, RC- resp. GS-scheme, at instant $t_* = 10$. Three sets of parameters have been chosen, summarized in Table 2. The discretization step is fixed at $\Delta t = 0.01$.</caption><thead><tr><th>t = 10</th><th>ε<sub>m</sub></th><th>ε<sub>f</sub></th><th>time</th><th>ε<sub>m</sub></th><th>ε<sub>f</sub></th><th>time</th><th>ε<sub>m</sub></th><th>ε<sub>f</sub></th><th>time</th></tr></thead><tbody><tr><td>DF (κ = 3)</td><td>2.3 · 10<sup>-3</sup></td><td>1.05 · 10<sup>-3</sup></td><td>0.048</td><td>8.08 · 10<sup>-5</sup></td><td>7.2 · 10<sup>-4</sup></td><td>1.025</td><td>8.11 · 10<sup>-5</sup></td><td>1.02 · 10<sup>-3</sup></td><td>20.5</td></tr><tr><td>RC (κ = 3)</td><td>4.75 · 10<sup>-5</sup></td><td>2.2 · 10<sup>-3</sup></td><td>0.045</td><td>4.75 · 10<sup>-5</sup></td><td>2.65 · 10<sup>-3</sup></td><td>0.1</td><td>4.75 · 10<sup>-5</sup></td><td>2.91 · 10<sup>-3</sup></td><td>1.82</td></tr><tr><td>DF (κ = 31)</td><td>7.5 · 10<sup>-14</sup></td><td>1.8 · 10<sup>-3</sup></td><td>0.047</td><td>1.42 · 10<sup>-14</sup></td><td>9.17 · 10<sup>-4</sup></td><td>1.05</td><td>3.06 · 10<sup>-14</sup></td><td>1.3 · 10<sup>-3</sup></td><td>20.5</td></tr><tr><td>GS (κ = 31)</td><td>4.44 · 10<sup>-16</sup></td><td>1.2 · 10<sup>-3</sup></td><td>0.024</td><td>4.44 · 10<sup>-16</sup></td><td>1.2 · 10<sup>-3</sup></td><td>0.033</td><td>4.44 · 10<sup>-16</sup></td><td>1.2 · 10<sup>-3</sup></td><td>0.044</td></tr></tbody></table>

velocity space, yields

$$ \langle f_{RC}^{(N)} \rangle = \sum_{k=0}^{\tilde{N}_{RC}} \mathfrak{a}_{2k}(t) \langle \Theta_{2k} \rangle, \quad \langle f_{GS}^{(N)} \rangle = \sum_{k=0}^{N_{GS}} \alpha_k(t) \langle p_k f_{\kappa,a} \rangle, $$

where we recall firstly that

$$ \langle p_k f_{\kappa,a} \rangle = (p_k, p_0)_{f_{\kappa,a}} = \gamma_0 \delta_{k,0} \quad \Rightarrow \quad \langle f_{GS}^{(N)} \rangle = \alpha_0(t) \gamma_0. $$

and secondly

$$ \Theta_{2k} = C_{2k} \Upsilon_{\kappa} = C_{2k} \frac{(c_{\kappa}/\sqrt{2\kappa})^{1/2}}{\left(1 + \frac{v^2}{2\kappa}\right)^{\frac{\kappa+1}{2}}} = (c_{\kappa}\sqrt{2\kappa})^{1/2} C_{2k} \frac{1}{\left(1 + \frac{v^2}{2\kappa}\right)^{\frac{\kappa-1}{2}}} \sigma_{\kappa}, $$

which, for $\kappa = 3$ and with (3.7), yields

$$ \langle f_{RC}^{(N)} \rangle = (c_{\kappa} \sqrt{2\kappa})^{1/2} \sum_{k=0}^{\tilde{N}_{RC}} \frac{\mathfrak{a}_{2k}(t)}{2} [(C_{2k}, C_0)_{\sigma_{\kappa}} - (C_{2k}, C_2)_{\sigma_{\kappa}}] = (c_{\kappa} \sqrt{2\kappa})^{1/2} \frac{\pi}{2} \left( \mathfrak{a}_0(t) - \frac{1}{2} \mathfrak{a}_2(t) \right). $$

The mass is also very precise (machine accuracy) for large $\kappa$-values, which comes from the fact that we have a more concentrated distribution function (close to a Maxwellien), such that no mass is lost with the FD-scheme. Compare this with the conservation of mass property of the FD-scheme, as shown in Figure 3. Secondly, the convergence rates of the solution towards the equilibrium is investigated in Figures 11 for both schemes. These results are better than expected from the (non-optimal) mathematical time-decay Theorems 1 and 2 for the energetic Fokker-Planck equation.
<figure><img src="image_15.png" /><figcaption>FIGURE 11. Relaxation of the distribution function $f(t, v)$ towards the equilibrium $f_\infty(v)$ in the long-time limit. Error plots $||f(t) - f_\infty||_{L^2_{f_\infty}}$ for the RC-scheme ($\kappa = 3$) on the left, and the GS-scheme ($\kappa = 31$) on the right.</figcaption></figure>

Interesting is also to investigate how the spectral coefficients, namely $a_{2k}(t)$ for the RC-spectral method resp. $\alpha_k(t)$ for the GS-scheme behave in time as well as with respect to the index $k$. The time-dependence of the first coefficients is very well illustrated on Figure 12. One observes that for the RC-scheme only two coefficients survive in the asymptotic limit $t \to \infty$, namely $a_0(t) \to_{t \to \infty} 0.4626$ and $a_2(t) \to_{t \to \infty} -0.4567$, which is consistent with the fact, that we took $\kappa = 3$. The GS-scheme requires in the asymptotic limit $t \to \infty$ only the first coefficient $\alpha_0(t) \to_{t \to \infty} 1.000525$. To be more efficient, it could be a good idea, to adapt the number of used coefficients with time, meaning using lesser and lesser coefficients as time goes on. Such a strategy has been used in [27] by one of the authors, is however not followed here. The magnitude of the different coefficients at a fixed time instant $t_* = 0.2$; 2; 10 is illustrated in Figure 13. These plots permit to see how much coefficients are needed at different time instants.

<figure><img src="image_16.png" /><figcaption>FIGURE 12. Time-evolution of the first spectral coefficients $a_k(t)$ resp. $\alpha_k(t)$. Left: RC-scheme, $\kappa = 3$. Right: GS-scheme, $\kappa = 31$. Furthermore, we took $N_{RC} = N_{GS} = 10$.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>FIGURE 13. Index-dependence of the spectral coefficients $\alpha_k(t_*)$ resp. $\alpha_k(t_*)$ for three different time instants $t_* = 0.2; 2; 10$. Left: RC-scheme, $\kappa = 3$ and $N_{RC} = 20$. Right: GS-scheme, $\kappa = 31$ and $N_{GS} = 10$.</figcaption></figure>

To conclude, the FD-scheme is relatively easy to implement, but is somehow inferior in accuracy and simulation time, in the long-time regime, as compared to our spectral methods. The RC-scheme is very efficient for small $\kappa$-values and long-time simulations, while the GS-scheme is very efficient for large $\kappa$-values and long-time simulations. In the limit $\kappa \to \infty$ the GS-scheme transforms in the very well-known Hermite spectral scheme, used in cases one has Maxwellian equilibria.

## 6. CONCLUDING REMARKS AND PERSPECTIVES

Let us conclude this paper by summarizing what was achieved in this work and what remains still to be done in future works. The main part of this work was concerned with the design of two efficient spectral schemes for the resolution of a specific Fokker-Planck equation, whose stationary states are given by $\kappa$-distribution functions, which are (thermal) non-equilibrium distributions and describe the energetic particle population in a fusion plasma gas or in space plasmas. At the heart of spectral methods is the fact that any smooth function $f(v)$ can be expanded in the form

$$f(v) = \sum_{k=0}^{\infty} \alpha_k \varphi_k(v), \quad \forall v \in \mathbb{R},$$

where $\{\varphi_k\}_{k \in \mathbb{N}}$ are global basis functions. An efficient scheme (accurate and stable) is achieved when the basis functions are adequately chosen. We presented in this paper two different polynomial basis sets for the resolution of the energetic Fokker-Planck equation, and compare these schemes with a standard finite difference scheme. Our conclusion is that spectral schemes (if well-designed) are superior in performance than local methods (FD), when long-time simulations are performed. For short-time simulations, standard FD-schemes are unbeatable. However, exponential accuracy is achieved for spectral schemes only when the basis functions are adapted to the specific problem, lacking such imperative can have a drastic impact on the convergence rate. In a future work we shall come closer to the physical reality, by including the diffusion coefficient $D(v)$ in the Fokker-Planck operator, as well as introducing a transport term (coupling with Poisson equation). The here presented schemes have naturally to be then adapted to the new problem.
**Acknowledgments.** This work has been carried out within the framework of the EURO-fusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 — EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. H. Parada was funded by the Agence Nationale de la Recherche, Labex CIMI under grant agreement ANR-11-LABX-0040.

REFERENCES

[1] P. Andries, K. Aoki, and B. Perthame *A consistent bgk-type model for gas mixtures*, Journal of Statistical Physics, **106** (2002), pp. 993–1018.

[2] B. Ayuso, J. A. Carrillo and C. W. Shu *Discontinuous Galerkin methods for the multi-dimensional Vlasov-Poisson problem*, Mathematical Models and Methods in Applied Sciences, **22**, (2021) pp. 1250042.

[3] D. Bakry, F. Barthe, P. Cattiaux, A. Guillin *A simple proof of the Poincaré inequality for a large class of probability measures*. Electronic Communications in Probability, 13(none) 60-66 2008.

[4] N. Ben Abdallah, A. Mellet, M. Puel, *Fractional diffusion limit for collisional kinetic equations: a Hilbert expansion approach*, Kinetic and Related Models **4** (2011), no. 4, Doi: 10.3934/krm.2011.4.873.

[5] S. Benzoni-Gavage *Calcul différentiel et équations différentielles*. Dunod, (2021).

[6] N. H. Bian, A. G. Emslie, D. J. Stackhouse, E. P. Kontar *The formation of kappa-distribution accelerated electron populations in solar flares*, The Astrophysical Journal **796** (2014).

[7] A. Blanchet, M. Bonforte, J. Dolbeault, and al. *Asymptotics of the Fast Diffusion Equation via Entropy Estimates*, Arch. Rational Mech. Anal. 191 (2009) 347–385

[8] A. Bobylev, M. Bisi, M. Groppi, G. Spiga, I. Potapenko *A general consistent BGK model for gas mixtures*, Kinet. Relat. Models **11** (2018), no. 6, 1377–1393.

[9] J. P. Boyd *Spectral methods using rational basis functions on an infinite interval*, JCP **69** (1987), 112–142.

[10] E. Bouin, J. Dolbeault, L. Ziviani L² *Hypocoercivity methods for kinetic Fokker-Planck equations with factorised Gibbs states*, (2023), preprint

[11] K.L. Cartwright, J. P. Verboncoeur and C. K. Birdsall *Nonlinear hybrid Boltzmann-particle-in-cell acceleration algorithm*, Physics of Plasmas (1994-present), 7.8 (2000): 3252-3264.

[12] L. Cesbron, A. Mellet, K. Trivisa, *Anomalous transport of particles in plasma physics*, Applied Mathematics Letters **25** (2012), no. 12, 2344–2348.

[13] F. F. Chen *Plasma Physics and controlled fusion*, Springer Verlag New York, (2006).

[14] Y. Chen, S. Parker *A gyrokinetic ion zero electron inertia fluid electron model for turbulence simulations*, Phys. Plasmas **8**, no. 2, 441–446 (2001)

[15] C. Cheverry, N. Raymond *A guide to spectral theory: applications and exercises*. Springer Nature, 2021.

[16] D. Clemence *On the Singular Behaviour of the Titchmarsh-Weyl m-Function for the Perturbed Hill's Equation on the Line*. Canadian Mathematical Bulletin. **40**, 416-421 (1997,12), Cambridge University Press

[17] E. S. Cramer, J. R. Dwyer, S. Arabshahi, I. B. Vodopiyanov, N. Liu, H. K. Rassoul *An analytical approach for calculating energy spectra of relativistic runaway electron avalanches in air*, J. Geophys. Res. Space Physics **119** (2014), 7794–7823.

[18] E.A. Coddington, N. Levinson *Theory of Ordinary Differential Equations*, Krieger Publishing, Malabar, FL. 1984.

[19] N. Crouseilles, C. Negulescu *Hybrid modelling of energetic particles and bulk plasma*, in preparation.

[20] P. Degond *Chapter 1 - Asymptotic Continuum Models for Plasmas and Disparate Mass Gaseous Binary Mixtures*, in Material Substructures in Complex Bodies, edited by Gianfranco CaprizPaolo Maria Mariano, Elsevier Science Ltd, Oxford, 2007, Pages 1-62.

[21] Y. Di, Y. Fan, Z. Kou, R. Li, Y. Wang *Filtered hyperbolic moment method for the Vlasov equation*, J. Sci. Comput. **79** (2019) 969–991.

[22] J. Dolbeault, G. Toscani *Fast diffusion equations: Matching large time asymptotics by relative entropy methods*, Kinetic and Related Models, **4** (2011) 701-716.
[23] J. Dominski, S. Brunner, S.K. Aghdam, T. Goerler, F. Jenko, D. Told Identifying the role of non-adiabatic passing electrons in ITG/TEM microturbulence by comparing fully kinetic and hybrid electron simulations, Journal of Physics: Conference Series **401** (2012).

[24] M. Eastham, C. Fulton, S. Pruess Using the SLEDGE package on Sturm-Liouville problems having nonempty essential spectra. *ACM Transactions On Mathematical Software.* **22**, 423-446 (1996,12,1).

[25] F. Filbet, T. Xiong *Conservative Discontinuous Galerkin/Hermite Spectral Method for the Vlasov-Poisson System*, Commun. Appl. Math. Comput., 2020

[26] F. Filbet, M. Bessemoulin-Chatard *On the stability of conservative discontinuous Galerkin/Hermite Spectral methods for the Vlasov-Poisson System*, J. Comp. Physics, 2022

[27] F. Filbet, C. Negulescu Fokker-Planck multi-species equations in the adiabatic asymptotics, Journal of Computational Physics **471** (2022).

[28] C. Fulton, S. Pruess *The computation of spectral density functions for singular Sturm-Liouville problems involving simple continuous spectra*. *ACM Transactions On Mathematical Software.* **24**, 107-129 (1998,3,1).

[29] G. Furioli, A. Pulvirenti, E. Terraneo, G. Toscnai, *One-Dimensional Fokker-Planck Equations and Functional Inequalities for Heavy Tailed Densities*, Milan J. Math. **90** (2022), 177-208.

[30] X. Garbet et al. *Global simulations of ion turbulence with magnetic shear reversal*, Physics of Plasmas (1994-present) 8.6 (2001): 2793-2803.

[31] W. Gautschi *Orthogonal Polynomials: Computation and Approximation*, Oxford University Press 2004.

[32] R. J. Goldston, P. H. Rutherford *Plasma Physics*, Taylor & Francis Group, (1995).

[33] H. Goedbloed, S. Poedts *Principles of Magnetohydrodynamics*, Cambridge University Press, Cambridge, (2004).

[34] H. Grad, *On the kinetic theory of rarefied gases*, Comm. Pure Appl. Math. **2** (1949), 331-407.

[35] J. Greene *Improved Bhatnagar-Gross-Krook model of electron-ion collisions*. Phys. Fluids **16**, 2022–2023 (1973)

[36] E. P. Gross and M. Krook *Model for collision processes in gases: Small-amplitude oscillations of charged two-component systems*, Phys. Rev., **102** (1956), pp. 593-604,

[37] A. V. Gurevich, J. A. Valdivia, G. M. Milikh, K. Papadopoulos *Runaway electrons in the atmosphere in the presence of a magnetic field*, Radio Science **31** (1996), no. 6, 1541–1554.

[38] B. B. Hamel *Kinetic model for binary gas mixtures*, Physics of Fluids, **8** (1965), pp. 418-425.

[39] A. Hasegawa, K. Mima, M. Duong-Van *Plasma distribution function in a suprathermal radiation field*, Phys. Reviex Letters **54** (1985), no. 24, 2608-2610.

[40] R.D. Hazeltine, J.D. Meiss *Plasma confinement*, Dover Publications, Inc. Mineola, New York (2003).

[41] M. Herda, L. M. Rodrigues *Anisotropic Boltzmann-Gibbs dynamics of strongly magnetized Vlasov-Fokker-Planck equations*, Kinet. Relat. Models **12** (2019), no. 3, 593-636.

[42] M. Klaus *On the variation-diminishing property of Schrödinger operators*. *Oscillations, Bifurcation And Chaos*. pp. 199-211 (1987)

[43] A. Kostenko, G. Teschl, J. H. Toloza *Dispersion Estimates for One-Dimensional Schrödinger Equations with Singular Potentials*, Annales Henri Poincaré **17** (2016), 3147-3176.

[44] H. Kovarík, F. Truc *Schrödinger Operators on a Half-Line with Inverse Square Potentials*. *Mathematical Modelling Of Natural Phenomena.* **9**, 170-176 (2014).

[45] C. Klingenberg, M. Pirner, G. Puppo *A consistent kinetic model for a two-component mixture with an application to plasma*, Kinetic and Related Models, **10** (2017), pp. 445-465.

[46] R. Kubo *The fluctuation-dissipation theorem*, Rep. Prog. Phys. **29** (1966), 255-284.

[47] D.T.K. Kwok *A hybrid Boltzmann electrons and PIC ions model for simulating transient state of partially ionized plasma*, Journal of Computational Physics 227.11 (2008): 5758-5777.

[48] R. Lavine *Absolute continuity of positive spectrum for Schrödinger operators with long-range potentials*, Journal of Functional Analysis, Volume 12, Issue 1, 1973, Pages 30-54,

[49] E. Lehman, C. Negulescu Fokker-Planck equation for energetic particles. *The kappa-distribution function*, submitted.

[50] R. Li, Y. Ren, Yinuo, Y. Wang *Hermite spectral method for Fokker-Planck-Landau equation modeling collisional plasma*, J. Comput. Phys. **434** (2021).

[51] G. Livadiotis *Kappa Distributions. Theory and Applications in Plasmas*, Elsevier Science Publishing Co Inc, (2017).

[52] A. Mellet *Diffusion Limit of a Non-Linear Kinetic Model without the Detailed Balance Principle*, Monatsh. Math. **134** (2002), 305-329.
[53] A. Mellet, S. Mischler, C. Mouhot, *Fractional diffusion limit for collisional kinetic equations*, Archive for Rational Mechanics and Analysis **199** (2011), 493–525.

[54] T. Miyazawa *Low-energy asymptotic expansion of the Green function for one-dimensional Fokker-Planck and Schrödinger equations*. *Journal Of Physics A: Mathematical And Theoretical.* **41**, 315304 (2008,7)

[55] S. Olbert *Summary of Experimental Results from M.I.T. Detector on IMP-1*, Physics of the Magnetosphere (1968), 641–659.

[56] F.W.J. Olver *Asymptotics and Special Functions*, A. K. Peters, Wellesley, MA, 1997.

[57] D. Pearson *Quantum scattering and spectral theory*. (Academic Press,1988)

[58] V. Pierrard, M. Lazar *Kappa distributions: Theory and applications in space plasmas*, Solar Phys. **267** (2010), 153–174.

[59] S. Pruess, C. Fulton *Mathematical software for Sturm-Liouville problems*. *ACM Transactions On Mathematical Software.* **19**, 360-376 (1993,9,1)

[60] M. Reed, B. Simon *I: Functional Analysis*. (Academic Press,1981,2)

[61] M. Röckner, F.Y. Wang *Weak Poincaré Inequalities and $L^2$-Convergence Rates of Markov Semigroups*, Journal of Functional Analysis **185** (2001), 564–603.

[62] J.W. Schumer, J.P. Holloway *Vlasov simulations using velocity-scaled Hermite representations*, J. Comput. Phys. **b**f 144 (2) (1998) 626–661.

[63] W. Schlag *Dispersive estimates for Schrödinger operators: a survey. Mathematical Aspects Of Nonlinear Dispersive Equations*. pp. 255-285 (2007), Publisher: Princeton Univ. Press

[64] A. Stahl *Momentum-space dynamics of runaway electrons in plasmas*, PhD thesis 2017, Department of Physics Chalmers University of Technology Gothenburg, Sweden.

[65] D. Stürzer, A. Arnold, *Spectral Analysis and Long-Time Behaviour of a Fokker-Planck Equation with a Non-Local Perturbation*, Rendiconti Lincei Cl. Sci. Fis. Mat. Natur **25** (2014), no. 1, 53–89.

[66] W. T. Taitano, B. D. Keenan, L. Chacón, S. E. Anderson, H. R. Hammer, A. N. Simakov *An Eulerian Vlasov-Fokker-Planck algorithm for spherical implosion simulations of inertial confinement fusion capsules* Comput. Phys. Commun. **263** (2021), Paper No. 107861, 42 pp.

[67] V.M. Vasyliunas *A survey of low-energy electrons in the evening sector of the magnetosphere with OGO 1 and OGO 3*, J. Geophys. Res. 73 (1968), 2839–2884.

[68] J. Vencels, G.L. Delzanno, A. Johnson, I.B. Peng, E. Laure and S. Markidis *Spectral solver for multi-scale plasma physics simulations with dynamically adaptive number of moments* Procedia Computer Science, **51**, (2015) pp.1148-1157.

[69] J. Wilkening, A. Cerfon *A spectral transform method for singular Sturm-Liouville problems with applications to energy diffusion in plasma physics*, SIAM J. Appl. Math., 75/2 (2015) 350-392.

[70] Z.-Q. Wang and B.-Y. Guo, *A rational approximation and its applications to nonlinear partial differential equations on the whole line*. *Journal of Mathematical Analysis and Applications*, **274(1)**, 2002, 374-403.

[71] D.R. Yafaev *The low energy scattering for slowly decreasing potentials*. *Communications in Mathematical Physics* **85** (1982): 177-196.

[72] D. Zwillinger *Table of Integrals, Series, and Products*, Academic Press, 8th edition, 2014.

UNIVERSITÉ DE TOULOUSE & CNRS, UPS, INSTITUT DE MATHÉMATIQUES DE TOULOUSE UMR 5219, F-31062 TOULOUSE, FRANCE.

Email address: claudia.negulescu@math.univ-toulouse.fr, hugo.parada@math.univ-toulouse.fr,