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

HAL
open science

# Hybrid modelling of energetic α-particles interacting with the thermal bulk plasma

Nicolas Crouseilles, Claudia Negulescu

► To cite this version:

Nicolas Crouseilles, Claudia Negulescu. Hybrid modelling of energetic α-particles interacting with the thermal bulk plasma. 2024. hal-04844525

HAL Id: hal-04844525

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

Preprint submitted on 17 Dec 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# HYBRID MODELLING OF ENERGETIC α-PARTICLES INTERACTING WITH THE THERMAL BULK PLASMA

NICOLAS CROUSEILLES, CLAUDIA NEGULESCU

**ABSTRACT.** The main concern of the present paper is the mathematical modelling and numerical simulation of thermonuclear fusion plasmas, constituted of two different kinds of particles, namely the thermal electron/ion bulk and an energetic α-particle population, created by the fusion reaction at very high speeds (3.5 MeV). This α-particle population behaves differently than the thermal plasma bulk, and due to its high energy contain, can have a considerable impact on the stability and thus confinement of the plasma bulk. An adequate modelling of this α-population is therefore crucial, and the present paper intends to make use of kinetic Fokker-Planck models to describe the evolution of these energetic particles experiencing Coulomb collisions. The main mathematical and numerical difficulties encountered in this study are related to the multi-scale nature of the overall problem, involving different types of particles and their multiscale dynamics.

**Keywords:** Fusion plasma modelling, energetic particles, thermal plasma bulk, Fokker-Planck kinetic equations, collisional processes, hybrid modelling, multi-scale numerical scheme.

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>2. The α-particle collision operator</td><td>3</td></tr><tr><td>2.1. Maxwellian approximation of the plasma background</td><td>4</td></tr><tr><td>2.2. Rescaling of the α-particle collision operator and simplifications</td><td>8</td></tr><tr><td>3. Study of the multi-scale α-particle collision operator</td><td>10</td></tr><tr><td>3.1. Multiscale α-particle slowing-down</td><td>10</td></tr><tr><td>3.2. Mathematical study</td><td>11</td></tr><tr><td>4. Multi-scale numerical scheme</td><td>14</td></tr><tr><td>4.1. Scaling of the Fokker-Planck kinetic equation</td><td>15</td></tr><tr><td>4.2. Velocity discretization</td><td>16</td></tr><tr><td>4.3. Time discretisation</td><td>16</td></tr><tr><td>4.4. Transport approximation</td><td>17</td></tr><tr><td>5. Numerical results and interpretations</td><td>18</td></tr><tr><td>5.1. Spatially homogeneous case</td><td>18</td></tr><tr><td>5.2. Including space transport</td><td>19</td></tr><tr><td>5.3. Fokker-Planck-Poisson</td><td>22</td></tr><tr><td>6. Concluding remarks and perspectives</td><td>24</td></tr><tr><td>References</td><td>26</td></tr></table>

Date: December 14, 2024.
# 1. INTRODUCTION

In the core region of thermonuclear tokamak fusion plasmas, which are mainly constituted of thermal ions and electrons, characterized by Maxwellian distribution functions (and called in the following bulk plasma), one also observes the occurrence of highly energetic particles. These particles are characterized by the fact that their kinetic energy $E_h$ is much higher than the thermal energy $k_B T_c$ of the bulk plasma, and their density is much lower, namely one has

$$E_h = \frac{1}{2}m|v_h|^2 \gg k_B T_c, \quad n_h \ll n_c,$$

where the index "h" stands for the energetic (hot) population and the index "c" for the thermal (cold) population. The energetic particles we shall consider in the present work are ion populations, generated via different mechanisms, such as

* neutral beam injection;

* nuclear fusion reactions (α-particles);

* radio frequency waves or ion cyclotron resonant heating.

Let us say some words about all these mechanisms. At ignition, the fusion plasma is self-heated by energetic ions (α-particles) produced during the fusion reaction at energies of about 3.5 MeV and which transfer their energy to the bulk plasma through collisional processes. Additional heating systems are however required to bring the fusion plasma to ignition. These are based either on the injection of highly energetic neutral particles, which become ionized as they penetrate into the dense bulk plasma and thermalize there through multiple Coulomb collisions. Another heating mechanism is based on the injection of electromagnetic wave energy into the plasma, in particular accelerating minority species of ions to high energies by radio-frequency waves, using a frequency equal to the cyclotron frequency of the minority ions (resonance). Altogether, the energetic particles generated via the mechanisms mentioned above, interact with the bulk plasma and have various impacts on the overall behaviour of the fusion plasma, for example

* driving instabilities, such as the TAE (toroidal Alfvén eigenmodes), fishbone modes;

* stabilizing sawteeth (internal kink modes);

* heating the thermal bulk via Coulomb collisional processes.

The success of magnetically confined fusion reactors relies upon a proper confinement of the energetic fusion products, and this for sufficiently long times such that they can transfer sufficient energy to the fuel ions in order to permit the fusion reaction to take place in a self-sustained manner. Therefore a detailed understanding of the impact of the energetic ion population on the confinement properties, on the global stability as well as on the heating of the fusion bulk plasma is of crucial importance for fusion reactor performances. More details about this physical context can be found in [9, 17, 20].

Widely used strategies to study the energetic particle populations are hybrid approaches, which combine a macroscopic description of the bulk plasma via a MHD model, and a mesoscopic description of the energetic particles (EP) via a kinetic equation [8, 26, 33]. Both models are coupled via current-coupling [7, 35] or pressure-coupling strategies [10]. These models have been introduced in the aim to account more precisely for the energetic particle population (which requires a kinetic description), evolving in a thermalized plasma bulk. However all the above mentioned models suffer from an inadequate description of the collisions between the energetic particles and the electron/ion thermal bulk. Very often the collision terms are even neglected due to computational reasons, and as a consequence, all
these models do not incorporate momentum and energy exchange terms between the energetic species and the thermal bulk, effects which are fundamental in particular for the heating of the fusion plasma by the α-particles. These existing hybrid models are hence valid only for time-scales smaller than the characteristic collisional times of the energetic particle species.

Several works exist in astrophysics [2, 34], modelling the transport of highly energetic particles accelerated during solar or stellar flares and releasing energy by interactions with the ambient atmosphere, until an eventual thermalization regime is reached. Such works are based on kinetic descriptions of the energetic particles, including several physical phenomena, such as Coulomb collisions described via Rosenbluth operators, relativistic particles, reconnection of magnetic field lines, radiation bursts etc. The aim of the present paper was to introduce more physics in the previously mentioned tokamak current-coupling resp. pressure-coupling MHD models, by proposing a more accurate, however still simple description of the energetic α-particle population (emerging during the fusion reactions), permitting thus to investigate the longer-time dynamics of these α-particle species and their influence on the thermal tokamak bulk plasma. Some physical effects from the astrophysical literature are included in our model, however we kept this model as simple as possible in order to concentrate on the multi-scale nature of the mathematical problem and propose an efficient numerical scheme. This paper concentrates firstly only on one dimension in the velocity space (focusing on the slowing-down of the energetic particles, drag and energy scattering), and prepares the ground for a future paper, taking more velocity dimensions into account, hence also the deflection scattering.

Our first goal is to construct a Fokker-Planck operator for the energetic ions (α-particles in the following) and to study it. The background plasma is assumed to be composed of Maxwellian electrons and ions, and we shall suppose that the density $n_α$ of the energetic ions is much lower than the density of the background plasma ($n_e \sim n_i$). The energetic ions are supposed to be injected or created at an initial velocity $v_{α,*}$ inbetween the thermal ion and electron velocities, thus satisfying $v_{th,i} \ll v_{α,*} \ll v_{th,e}$, which is a physically confirmed hypothesis [16, 32]. In a first step, the result of the collisions with the plasma bulk will be mainly frictional drag (which will cause the slowing-down of the energetic ions) and angular scattering (which will cause the deflection of the energetic ions from their original direction). By this manner, the α-particles will transfer their energy to the surrounding bulk plasma and will leave the energetic velocity regime, without becoming Maxwellians, as the inter-species α-particle collisions are rare. In a second, final step however, the energetic α-particle population becomes also thermalized and can be accurately described by a Maxwellian in the long-time limit. To describe this multi-scale α-particle evolution as accurately as possible, a kinetic framework is required for the energetic species and an adequate coupling with the thermal bulk plasma, via well-chosen collision operators. This is the goal of the present paper. The main mathematical and numerical difficulties of the obtained Fokker-Planck equation are linked to the multi-scale nature of the problem, thus requiring a delicate asymptotic analysis.

The outline of this paper is the following. Section 2 presents the physical framework and introduces an α-particle Fokker-Planck collision operator, describing the collisions of the energetic species with the thermal bulk. Section 3 is concerned with the mathematical study of a simplified multi-scale α-particle Fokker-Planck equation (in one dimension and omitting the transport term), focusing on the existence and uniqueness results, as well as on the $t \to \infty$ asymptotics. Section 4 proposes an efficient multi-scale numerical scheme for the resolution of this α-particle Fokker-Planck equation, the difficulty being to cope properly
with the different time and velocity scales. Finally, we test our scheme in a physical fusion test case in Section 5 and conclude the paper with some remarks and prespectives.

## 2. THE α-PARTICLE COLLISION OPERATOR

Let us introduce in this section the collision operator corresponding to the interaction of the α-particles with the bulk, constituted of thermalized electrons and ions. The evolution of the fast α-particles is described via a distribution function $f_{\alpha}(t, \mathbf{x}, \mathbf{v})$ which is solution of a kinetic equation of the form

$$ (2.1) \qquad \partial_t f_{\alpha} + \mathbf{v} \cdot \nabla_{\mathbf{x}} f_{\alpha} + \frac{e_{\alpha}}{m_{\alpha}} (\mathbf{E} + \mathbf{v} \times \mathbf{B}) \cdot \nabla_{\mathbf{v}} f_{\alpha} = Q_{\alpha\alpha} + Q_{\alpha e} + Q_{\alpha i}, $$

defined for times $t > 0$ and evolving on the phase-space $(\mathbf{x}, \mathbf{v}) \in \mathbb{T}^3 \times \mathbb{R}^3$ ($\mathbb{T}^3$ being the three-dimensional torus), whereas $m_{\alpha}$ is the mass and $e_{\alpha}$ the elementary charge of the considered α-particles. The electro-magnetic fields $(\mathbf{E}, \mathbf{B})$ are computed usually via Maxwell's equations. The intra-species collisions between α-particles $Q_{\alpha\alpha}$ will be neglected in the following due to the assumption of low α-particle densities, i.e. $n_{\alpha} \ll n_e \sim n_i$, meaning that α-particles do not encounter very often, such that collisions within their population are negligible. However the collisions with the thermalized plasma bulk are important, relaxing the α-particle population and permitting the energy-transfer to the bulk.

Our starting point for the inter-species collisions is the Rosenbluth collision operator [20, 21], which is an alternative expression of the Landau collision operator and is given by

$$ (2.2) \qquad Q_{\alpha s}(\mathbf{v}) := 4\pi \ln \Lambda_{\alpha s} \left( \frac{e_{\alpha} e_{s}}{m_{\alpha} \varepsilon_{0}} \right)^{2} \nabla_{\mathbf{v}} \cdot \left[ -\frac{m_{\alpha}}{m_{s}} \nabla_{\mathbf{v}} \varphi_{s} f_{\alpha} + \frac{1}{2} H_{\psi_{s}} \nabla_{\mathbf{v}} f_{\alpha} \right], \quad \forall \mathbf{v} \in \mathbb{R}^{3}, $$

where $s$ stands for either the thermalized electron species ($s = e$) or the thermalized ion species ($s = i$), $\ln \Lambda_{\alpha s}$ is the Coulomb logarithm (of reasonable size $\sim 20-40$), $m_s$ resp. $m_{\alpha}$ are the masses of the corresponding particles, $e_s$ resp. $e_{\alpha}$ their elementary charge and finally $\varepsilon_0$ is the vacuum permittivity. We denoted by $H_{\chi}$ the Hesse-matrix corresponding to a function $\chi$ and the functions $\varphi_s$ and $\psi_s$ are the so-called Rosenbluth potentials, defined through the particle distribution functions $f_s$ as follows

$$ (2.3) \qquad \varphi_s(\mathbf{v}) := \int_{\mathbb{R}^3} \frac{1}{|\mathbf{v} - \mathbf{v}'|} f_s(\mathbf{v}') \, d^3\mathbf{v}', \quad \psi_s(\mathbf{v}) := \int_{\mathbb{R}^3} |\mathbf{v} - \mathbf{v}'| f_s(\mathbf{v}') \, d^3\mathbf{v}'. $$

Remark that one has the following Poisson equations for the computation of the Rosenbluth potentials

$$ (2.4) \qquad -\Delta_{\mathbf{v}}\varphi_s = 4\pi f_s, \quad -\Delta_{\mathbf{v}}\psi_s = -2\varphi_s. $$

These potentials describe the cumulative contribution of the background plasma bulk to the drag and diffusion terms in the α-particle collision operator.

Denoting, for simplicity reasons, the integration in the velocity variable by $\langle \xi \rangle = \int_{\mathbb{R}^3} \xi(\mathbf{v}) \, d^3\mathbf{v}$, one can show that the inter-species Rosenbluth collision operators $Q_{\alpha s}$ (2.2) preserve certain physical properties, such as the mass $\langle Q_{\alpha s} \rangle = 0$, the momentum

$$ (2.5) \qquad m_{\alpha} \langle Q_{\alpha s} \mathbf{v} \rangle + m_{s} \langle Q_{s\alpha} \mathbf{v} \rangle = 0, \quad \forall s \in \{e, i\}, $$

the total energy

$$ (2.6) \qquad m_{\alpha} \langle Q_{\alpha s} \frac{|\mathbf{v}|^2}{2} \rangle + m_{s} \langle Q_{s\alpha} \frac{|\mathbf{v}|^2}{2} \rangle = 0, \quad \forall s \in \{e, i\}, $$
and satisfy the entropy decay relation

$$m_{\alpha} \langle Q_{\alpha s} \ln(f_{\alpha}) \rangle + m_{s} \langle Q_{s\alpha} \ln(f_{s}) \rangle \le 0, \quad \forall s \in \{e, i\}.$$

These relations are useful for the coupling of the α-particle dynamics with the thermal bulk, described usually via fluid (MHD) models.

The Rosenbluth form is rather cumbersome for numerical simulations, such that several approximations shall be considered in the following in order to simplify it. These simplifications will be based on the following assumptions:

* thermalized plasma bulk, meaning we shall assume $f_e = M_{n_e, u_e, T_e}$ and $f_i = M_{n_i, u_i, T_i}$;

* thermal bulk characteristics: $T_e \sim T_i$ and $n_e \sim n_i$; $m_e \ll m_i = m_\alpha$;

* $n_\alpha \ll n_e \sim n_i$, leading to the disregard of the intra-species collision operator $Q_{\alpha\alpha}$;

* $v_{th,i} \ll v_{\alpha,\star} \ll v_{th,e}$, hence α-particles are created with a velocity inbetween the thermal ion and thermal electron velocities, where in a general manner we define the thermal velocity via $v_{th,s} := \sqrt{\frac{2k_B T_s}{m}}$, for $s \in \{e, i\}$.

## 2.1. Maxwellian approximation of the plasma background.
In the case the species s is described by a Maxwellian distribution function, one can compute analytically the Rosenbluth potentials and simplify the collision operators. Let us first start with the simpler case of zero mean Maxwellians, and generalize then the expressions in subsection 2.1.2.

### 2.1.1. Maxwellian background with zero mean velocities.
Let us suppose in this subsection that the bulk particles, namely the thermal electrons and ions, have zero mean velocities $\mathbf{u}_e = \mathbf{u}_i = 0$, and are described by $f_s := M_{n_s, 0, T_s}$. In this case, due to the fact that $f_s$ is spherically symmetric, the Rosenbluth potentials are likewise spherically symmetric and can be written (see [21, 34] for more details and setting $v := |\mathbf{v}|$) as

$$\varphi_s(v) = \frac{n_s}{v} \phi(x_s), \quad x_s := \frac{v}{v_{th,s}}, \quad \forall v \in \mathbb{R},$$

$$\psi_s(v) = n_s v_{th,s} \left[ \frac{1}{2} \phi'(x_s) + \left( \frac{1}{2x_s} + x_s \right) \phi(x_s) \right], \quad v_{th,s} := \sqrt{\frac{2k_B T_s}{m_s}},$$

where the error function $\phi \in C^\infty(\mathbb{R})$ and the Chandrasekhar function $G \in C^\infty(\mathbb{R})$ are defined as (see Figure 1)

$$(2.7) \qquad \phi(x) := \frac{2}{\sqrt{\pi}} \int_0^x e^{-y^2} dy; \qquad G(x) := \frac{\phi(x) - x \phi'(x)}{2x^2}, \quad G(0) := 0.$$

This permits to simplify the Rosenbluth collision operator (2.2), due to

$$\varphi'_{s}(v) = -2 \frac{n_s}{v_{th,s}^2} G(x_s), \quad \nabla_{\mathbf{v}} \varphi_{s}(v) = \varphi'_{s}(v) \frac{\mathbf{v}}{v}, \quad \psi'_{s}(v) = n_s [\phi(x_s) - G(x_s)],$$

$$\psi''_{s}(v) = 2 \frac{n_s}{v} G(x_s), \quad H_{\psi_s} = \left(\mathbb{I} - \frac{\mathbf{v} \otimes \mathbf{v}}{v^2}\right) \frac{1}{v} \psi'_{s}(v) + \frac{\mathbf{v} \otimes \mathbf{v}}{v^2} \psi''_{s}(v),$$

with $\mathbb{I}$ the identity matrix. Remark that $H_{\psi_s}$ contains two parts, a projection on the plane perpendicular to the velocity $\mathbf{v}$ and a projection in the velocity vector direction. To render the ensuing expressions more readable, let us introduce the following notation

$$v_0^{\alpha s} := c_{\alpha s} \frac{n_s}{v_{th,s}^3}, \quad c_{\alpha s} := 4\pi \ln \Lambda_{\alpha s} \left( \frac{e_\alpha e_s}{m_\alpha \varepsilon_0} \right)^2,$$
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Left: Plot of the odd error function $\phi(x)$. Right: Plot of the odd Chandrasekhar function $G(x)$.</figcaption></figure>

as well as the deflection, slowing down and the energy scattering collisional frequencies

$$ (2.8) \qquad \begin{aligned} \nu_D^{\alpha s} &:= \nu_0^{\alpha s} \frac{\phi(x_s) - G(x_s)}{x_s^3}, & x_s &:= v/v_{th,s}, \\ \nu_s^{\alpha s} &:= 2 \nu_0^{\alpha s} \left(1 + \frac{m_\alpha}{m_s}\right) \frac{G(x_s)}{x_s}, \\ \nu_{\|\|}^{\alpha s} &:= 2 \nu_0^{\alpha s} \frac{G(x_s)}{x_s^3}. \end{aligned} $$

Passing now in (2.2) to spherical coordinates $(v, \theta, \varphi) \in \mathbb{R}^+ \times [0, \pi] \times [0, 2\pi)$ in the velocity variable, yields the following formulae for the collision operator [21, 34] (2.9)

$$ Q_{\alpha s} = \nu_D^{\alpha s} \mathcal{L}(f_\alpha) + \frac{1}{v^2} \partial_v [v^2 \Gamma_v^{\alpha s}(f_\alpha)], \qquad \Gamma_v^{\alpha s}(f_\alpha) := v \left[ \nu_s^{\alpha s} \left( \frac{m_\alpha}{m_\alpha + m_s} \right) f_\alpha + \frac{1}{2} \nu_{\|\|}^{\alpha s} v \partial_v f_\alpha \right], $$

where $\Gamma_v^{\alpha s}(f_\alpha)$ is the collisional flux due to the drag and energy scattering and $\mathcal{L}(f_\alpha)$ is the Lorentz scattering operator, defined in spherical resp. pitch angle coordinates $(\xi := \cos\theta)$, as (2.10)

$$ \mathcal{L}(f_\alpha) := \frac{1}{2} \left( \frac{1}{\sin\theta} \partial_\theta[\sin\theta \partial_\theta f_\alpha] + \frac{1}{\sin^2\theta} \partial_\varphi^2 f_\alpha \right) = \frac{1}{2} \partial_\xi [(1 - \xi^2) \partial_\xi f_\alpha] + \frac{1}{2} \frac{1}{(1 - \xi^2)} \partial_\varphi^2 f_\alpha. $$

Remark that $\mathcal{L}$ corresponds to the pitch-angle diffusion operator (no energy scattering), coming from

$$ \nabla_{\mathbf{v}} \cdot [W(\mathbf{v}) \nabla_{\mathbf{v}} f_{\alpha}] = \frac{2}{v^3} \mathcal{L}(f_{\alpha}), \qquad W(\mathbf{v}) := \frac{1}{v} \left( \mathbb{I} - \frac{\mathbf{v} \otimes \mathbf{v}}{v^2} \right), $$

and tends to make the $\alpha$-particle distribution function isotropic in $(\theta, \varphi)$. The deflection collision frequency $\nu_D^{\alpha s}$ determines how quickly the direction of the velocity vector of the $\alpha$-particle changes during a collision with species $s$, or in other words it gives the rate at which an $\alpha$-particle is diffusing in the plane perpendicular to the velocity $\mathbf{v}$, due to collisions with the bulk. The slowing-down frequency $\nu_s^{\alpha s}$ describes the rate at which an $\alpha$-particle is decelerated by a collision with species $s$ (drag) and finally the energy scattering collisional frequency $\nu_{\|\|}^{\alpha s}$ gives the parallel velocity diffusion frequency.
The collision operator (2.9) can be also rewritten, making use of only $\nu_0^{\alpha s}$, as follows

$$ (2.11) \quad Q_{\alpha s} = \nu_0^{\alpha s} \frac{\phi(x_s) - G(x_s)}{x_s^3} \mathcal{L}(f_\alpha) + 2 \nu_0^{\alpha s} \frac{1}{v^2} \partial_v \left[ v^2 \frac{G(x_s)}{x_s} \left( \frac{m_\alpha}{m_s} v f_\alpha + \frac{1}{2} v_{th,s}^2 \partial_v f_\alpha \right) \right]. $$

No approximation has been performed to obtain this collision operator from (2.9), we only assumed that species s is described by a Maxwellian distribution function with zero mean, such that the computations of the Rosenbluth potentials could be done analytically. Hence, all the conservation laws (see (2.5)-(2.6)) are still valid in this case. Later on, we shall however make use of some assumptions, such as $v_{th,i} \ll v_{\alpha,\star} \ll v_{th,e}$ as well as $m_e \ll m_i$, in order to render this operator more tractable for numerical simulations. The asymptotic developments of the error and Chandrasekhar functions, for small as well as large arguments, are then required and summarized in Table 1.

<table><caption>TABLE 1. Asymptotic developments of the error and Chandrasekhar functions [37].</caption><thead><tr><th></th><th>0 &lt; x &lt;&lt; 1</th><th>x &gt;&gt; 1</th></tr></thead><tbody><tr><td>φ(x) ~</td><td><sup>2x</sup>&frasl;<sub>&radic;π</sub> (1 - <sup>x<sup>2</sup></sup>&frasl;<sub>3</sub> + <sup>x<sup>4</sup></sup>&frasl;<sub>10</sub> ...)</td><td>1 - <sup>e<sup>-x<sup>2</sup></sup></sup>&frasl;<sub>&radic;πx</sub> (1 - <sup>1</sup>&frasl;<sub>2x<sup>2</sup></sub> + <sup>3</sup>&frasl;<sub>4x<sup>4</sup></sub> ...)</td></tr><tr><td>G(x) ~</td><td><sup>2x</sup>&frasl;<sub>3&radic;π</sub> - <sup>2x<sup>3</sup></sup>&frasl;<sub>5&radic;π</sub> ...</td><td><sup>1</sup>&frasl;<sub>2x<sup>2</sup></sub> - <sup>e<sup>-x<sup>2</sup></sup></sup>&frasl;<sub>&radic;πx</sub> (1 + <sup>1</sup>&frasl;<sub>2x<sup>2</sup></sub> - <sup>1</sup>&frasl;<sub>4x<sup>4</sup></sub> ...)</td></tr></tbody></table>

Remark that these developments permit to picture out the shape of the coefficient $\frac{G(x_s)}{x_s}$, occurring in the second term of the collision operator (2.11) (the rather standard Fokker-Planck term), and plotted on Figure 2. It behaves as $1/x^3$ for large values of $x$, and is bounded for small $x$-values, which signifies that with larger and larger velocities, the $\alpha$-particles are lesser and lesser slowed down (specificity of energetic particles).

<figure><img src="image_4.png" /><figcaption>FIGURE 2. Plot of the coefficient term $G(x)/x$, occurring in (2.11).</figcaption></figure>

### 2.1.2. Shifted Maxwellian background.
Let us now expose the differences when supposing non-vanishing mean velocities for the bulk ion/electron plasma. Starting again from (2.2) and inserting the shifted Maxwellians $f_s := M_{n_s, u_s, T_s}$ into the Rosenbluth potentials (2.3), yields

$$ (2.12) \quad Q_{\alpha s} := 4\pi \ln \Lambda_{\alpha s} \left( \frac{e_{\alpha} e_{s}}{m_{\alpha} \epsilon_{0}} \right)^{2} \nabla_{\mathbf{v}} \cdot \left[ -\frac{m_{\alpha}}{m_{s}} \nabla_{\mathbf{v}} \tilde{\varphi}_{s} f_{\alpha} + \frac{1}{2} H_{\tilde{\psi}_{s}} \nabla_{\mathbf{v}} f_{\alpha} \right], $$
where this time

$$ \nabla_{\mathbf{v}} \tilde{\varphi}_s(\mathbf{v}) = -\frac{2n_s}{v_{th,s}^2} G\left(\frac{|\mathbf{v} - \mathbf{u}_s|}{v_{th,s}}\right) \frac{\mathbf{v} - \mathbf{u}_s}{|\mathbf{v} - \mathbf{u}_s|}, \quad x_s := \frac{|\mathbf{v} - \mathbf{u}_s|}{v_{th,s}}, $$

$$ H_{\tilde{\psi}_s}^{\alpha}(\mathbf{v}) = \frac{n_s}{|\mathbf{v} - \mathbf{u}_s|} \left\{ \mathbb{I} - \frac{(\mathbf{v} - \mathbf{u}_s) \otimes (\mathbf{v} - \mathbf{u}_s)}{|\mathbf{v} - \mathbf{u}_s|^2} \right\} [\phi(x_s) - G(x_s)] \\ + \frac{2 n_s}{|\mathbf{v} - \mathbf{u}_s|} \frac{(\mathbf{v} - \mathbf{u}_s) \otimes (\mathbf{v} - \mathbf{u}_s)}{|\mathbf{v} - \mathbf{u}_s|^2} G(x_s). $$

Defining now the collisional frequencies as

$$ \nu_0^{\alpha s} := c_{\alpha s} \frac{n_s}{v_{th,s}^3}, \quad c_{\alpha s} := 4 \pi \left( \frac{e_{\alpha} e_s}{m_{\alpha} \varepsilon_0} \right)^2 \ln \Lambda_{\alpha s} $$

$$ \nu_D^{\alpha s} := \nu_0^{\alpha s} \frac{v_{th,s}^3}{|\mathbf{v} - \mathbf{u}_s|^3} [\phi(x_s) - G(x_s)], $$

$$ \nu_s^{\alpha s} := 2 \nu_0^{\alpha s} \frac{v_{th,s}}{|\mathbf{v} - \mathbf{u}_s|} \left(1 + \frac{m_{\alpha}}{m_s}\right) G(x_s), $$

$$ \nu_{\parallel}^{\alpha s} := 2 \nu_0^{\alpha s} \frac{v_{th,s}^3}{|\mathbf{v} - \mathbf{u}_s|^3} G(x_s), $$

permits to rewrite the collision operator as

$$ Q_{\alpha s}(\mathbf{v}) := \nabla_{\mathbf{v}} \cdot \left[ \nu_s^{\alpha s} (\mathbf{v} - \mathbf{u}_s) f_{\alpha}(\mathbf{v}) - 2 \nu_0^{\alpha s} v_{th,s} G(x_s) \frac{\mathbf{v} - \mathbf{u}_s}{|\mathbf{v} - \mathbf{u}_s|} f_{\alpha}(\mathbf{v}) \right] \\ + \frac{1}{2} \nabla_{\mathbf{v}} \cdot \left[ \nu_D^{\alpha s} |\mathbf{v} - \mathbf{u}_s|^2 \left( \mathbb{I} - \frac{(\mathbf{v} - \mathbf{u}_s) \otimes (\mathbf{v} - \mathbf{u}_s)}{|\mathbf{v} - \mathbf{u}_s|^2} \right) \cdot \nabla_{\mathbf{v}} f_{\alpha} \right] \\ + \frac{1}{2} \nabla_{\mathbf{v}} \cdot \left[ \nu_{\parallel}^{\alpha s} (\mathbf{v} - \mathbf{u}_s) \otimes (\mathbf{v} - \mathbf{u}_s) \cdot \nabla_{\mathbf{v}} f_{\alpha} \right]. $$

One has now to introduce the new shifted velocity variables $\mathbf{w} := \mathbf{v} - \mathbf{u}_s$ and to pass then to spherical variable $(w, \theta_w, \varphi_w)$ in $\mathbf{w}$ to recover thus the collision operator

$$ Q_{\alpha s}(w) = \nu_0^{\alpha s} \frac{\phi(x_s) - G(x_s)}{x_s^3} \mathcal{L}(f_{\alpha}) + 2 \nu_0^{\alpha s} \frac{1}{w^2} \partial_w \left[ w^2 \frac{G(x_s)}{x_s} \left( \frac{m_{\alpha}}{m_s} w f_{\alpha} + \frac{1}{2} v_{th,s}^2 \partial_w f_{\alpha} \right) \right], $$

which has the same form as (2.11).

## 2.2. Rescaling of the α-particle collision operator and simplifications.
To identify the dominant mechanisms, let us start from the case of zero mean velocities for the bulk plasma (Subsection 2.1.1) and rescale the electron resp. ion part of the α-particle collision operator (2.11), by making use of the assumptions $v_{th,i} \ll v_{\alpha,\star} \ll v_{th,e}$, $m_e/m_i \ll 1$, $n_e \sim n_i \gg n_\alpha$ and $T_e \sim T_i$. We recall here, that $v_{\alpha,\star}$ is the creation velocity of the α-particles, these ones being afterwards slowed-down by collisions towards smaller and smaller velocities. Let us introduce some small parameters describing the situation we are interested to describe. Firstly we define

$$ \varepsilon := \frac{v_{th,i}}{v_{\alpha,\star}} = \frac{v_{\alpha,\star}}{v_{th,e}} \ll 1, \quad \delta^2 := \frac{m_e}{m_i} \ll 1. $$

Remarking that $v_{th,s} := \sqrt{\frac{2k_B T_s}{m_s}}$, one has

$$ \frac{v_{th,i}}{v_{th,e}} = \frac{v_{th,i}}{v_{\alpha,\star}} \frac{v_{\alpha,\star}}{v_{th,e}} = \varepsilon^2 \quad \text{but also} \quad \frac{v_{th,i}}{v_{th,e}} = \sqrt{\frac{T_i m_e}{T_e m_i}} \sim \delta, $$
such that we will fix $\delta := \epsilon^2$ in the following and shall continue with only one fixed, but small parameter $\epsilon \in (0, 1)$, having thus a mass ratio of $\frac{m_e}{m_i} = \epsilon^4$.

Having introduced the characteristic scales of our problem, let us rescale now the $\alpha$-particle collision operator (2.11), performing the following change of variable

$$v = v_{\alpha, \star} v', \quad f_{\alpha}(v) = \bar{f}_{\alpha} f'_{\alpha}(v'), \quad v' \in [0, 1],$$

which leads to

$$Q_{\alpha s}(v') = \nu_0^{\alpha s} \bar{f}_{\alpha} \frac{\phi(x_s) - G(x_s)}{x_s^3} \mathcal{L}(f'_{\alpha}) + 2 \nu_0^{\alpha s} \bar{f}_{\alpha} \frac{1}{(v')^2} \partial_{v'} \left[ (v')^2 \frac{G(x_s)}{x_s} \left( \frac{m_{\alpha}}{m_s} v' f'_{\alpha} + \frac{1}{2} \frac{v_{th,s}^2}{v_{\alpha, \star}^2} \partial_{v'} f'_{\alpha} \right) \right].$$

where we recall that $\nu_0^{\alpha s} = c_{\alpha s} \frac{n_s}{v_{th,s}^3}$ and $x_s = \frac{v_{\alpha, \star}}{v_{th,s}} v'$. Remarking that $c_{\alpha i} \sim c_{\alpha e}$, let us rewrite the collision frequency as

$$\nu_0^{\alpha s} = c_{\alpha s} \frac{n_{\alpha}}{v_{\alpha, \star}^3} \frac{n_s}{n_{\alpha}} \frac{v_{\alpha, \star}^3}{v_{th,s}^3} = \bar{\nu}_0 \frac{n_s}{n_{\alpha}} \frac{v_{\alpha, \star}^3}{v_{th,s}^3}, \quad \text{setting} \quad \bar{\nu}_0 := c_{\alpha s} \frac{n_{\alpha}}{v_{\alpha, \star}^3},$$

and introducing furthermore the characteristic value $\bar{Q} := \bar{\nu}_0 \bar{f}_{\alpha}$. The $\alpha$-particle density being small as compared to the electron/ion bulk density, we introduce finally another small parameter $\beta \ll 1$ defined as $\beta := \frac{\bar{\nu}_{\alpha}}{\bar{n}_s}$, yielding (2.13)

$$Q_{\alpha s}(v') = \frac{\bar{Q}}{\beta} \frac{v_{\alpha, \star}^3}{v_{th,s}^3} \left\{ \frac{\phi(x_s) - G(x_s)}{x_s^3} \mathcal{L}(f'_{\alpha}) + \frac{2}{(v')^2} \partial_{v'} \left[ (v')^2 \frac{G(x_s)}{x_s} \left( \frac{m_{\alpha}}{m_s} v' f'_{\alpha} + \frac{1}{2} \frac{v_{th,s}^2}{v_{\alpha, \star}^2} \partial_{v'} f'_{\alpha} \right) \right] \right\}.$$

This rescaled operator (2.13) will be the starting point of the next sections. But before passing to the mathematical and numerical investigations, let us comment at the end of this section a little bit on this operator, to get a feeling of its action, in particular focusing on its asymptotic, simplified forms in the regime of small $\epsilon \ll 1$. For this we summarized in Table 2 the asymptotic behaviours of some of its coefficients.

<table><caption>TABLE 2. Asymptotic developments of the coefficients occurring in (2.13).</caption><thead><tr><th></th><th>0 &lt; x &lt;&lt; 1</th><th>x &gt;&gt; 1</th></tr></thead><tbody><tr><td>$\frac{\phi(x)-G(x)}{x^3} \sim \frac{1}{x^3} \left[ \frac{4x}{3\sqrt{\pi}} - \frac{4x^3}{15\sqrt{\pi}} + \mathcal{O}(x^5) \right]$</td><td>$\frac{1}{x^3} \left[ 1 - \frac{1}{2x^2} + \frac{1}{x^3\sqrt{\pi}} e^{-x^2} + \mathcal{O}(x^{-5}) \right]$</td></tr><tr><td>$\frac{G(x)}{x} \sim \frac{1}{x} \left[ \frac{2x}{3\sqrt{\pi}} - \frac{2x^3}{5\sqrt{\pi}} + \mathcal{O}(x^5) \right]$</td><td>$\frac{1}{x} \left[ \frac{1}{2x^2} - \frac{e^{-x^2}}{x\sqrt{\pi}} \left( 1 + \frac{1}{2x^2} + \mathcal{O}(x^{-4}) \right) \right]$</td></tr></tbody></table>

For the electrons the situation is rather simple, as we always have $x_e = \frac{v_{\alpha, \star}}{v_{th,e}} v' = \epsilon v' \ll 1$, such that skipping the primes in (2.13) for notational simplicity leads to the rescaled and simplified (suprathermal) electron collision operator (in the new variables)

$$(2.14) \qquad Q_{\alpha e}^{ST} = \frac{4}{3\sqrt{\pi}} \frac{\bar{Q}}{\epsilon\beta} \left\{ \frac{\epsilon^2}{v^2} \mathcal{L}(f_{\alpha}) + \frac{1}{v^2} \partial_v \left[ v^2 \left( v f_{\alpha} + \frac{\epsilon^2}{2} \partial_v f_{\alpha} \right) \right] \right\}.$$

For the ions, we have to consider two cases. Immediately after the $\alpha$-particle creation, we are in the situation $v_{th,i} \ll v_{\alpha, \star} v' \ll v_{th,e}$, such that $x_i \gg 1$, which leads to the rescaled and
simplified (suprathermal) collision operator

$$ (2.15) \qquad Q_{\alpha i}^{ST} = \frac{\bar{Q}}{\beta} \left\{ \frac{1}{v^3} \mathcal{L}(f_{\alpha}) + \frac{1}{v^2} \partial_v \left[ v^2 \frac{1}{v^3} \left( v f_{\alpha} + \frac{\varepsilon^2}{2} \partial_v f_{\alpha} \right) \right] \right\}. $$

However, when the $\alpha$-particles are so slowed down, that their velocity is of the order of $v_{th,i}$ or even lower, the approximation (2.15) for the collision operator $Q_{\alpha i}$ is no more valid, and one has to go back to the more general form (2.13). In the extreme case of $x_i \ll 1$ one can simplify this operator to the thermal form

$$ (2.16) \qquad Q_{\alpha i}^{T} = \frac{4}{3\sqrt{\pi}} \frac{\bar{Q}}{\varepsilon^3 \beta} \left\{ \frac{\varepsilon^2}{v^2} \mathcal{L}(f_{\alpha}) + \frac{1}{v^2} \partial_v \left[ v^2 \left( v f_{\alpha} + \frac{\varepsilon^2}{2} \partial_v f_{\alpha} \right) \right] \right\}. $$

The rescaled versions (2.14)-(2.15) permit to identify the dominant collisional processes during the primary suprathermal $\alpha$-particle evolution process. One remarks that once created, the $\alpha$-particles are firstly slowed down by the Coulomb collisions with the electrons (dominant term $\sim 1/\varepsilon$). The deflection collision frequency and the energy-scattering frequency are negligible for electrons, and this due to the small electron-to-ion mass ratio. Then when slowed-down up to a critical velocity $v_c$, the thermal ions take the overhand and dominate the final suprathermal relaxation process of the $\alpha$-particles (slowing-down and deflection are now comparable). The critical velocity $v_c$, called also crossover velocity [16] as it delimitates the main influence of each operator, can be computed by equating the electron to the ion frictional drags in (2.14)-(2.15), namely putting

$$ \frac{4}{3\sqrt{\pi}} n_e c_{\alpha e} \frac{1}{\varepsilon} = n_i c_{\alpha i} \frac{1}{v_c^3}, $$

leading to

$$ v_c := \left( \frac{3\sqrt{\pi}}{4} \frac{n_i \ln \Lambda_{\alpha i}}{n_e \ln \Lambda_{\alpha e}} \varepsilon \right)^{1/3}. $$

Finally, when the $\alpha$-particle velocity becomes close to the thermal ion bulk velocity $v_{th,i}$ or even lower, the thermal collision operator (2.16) takes the overhand. We are thus in a typical situation of a multiscale process, and to be able to describe accurately all the divers phenomena, we shall separate in Section 4 the velocity domain in distinct parts, another good idea being the separation of the $\alpha$-particle population in a suprathermal and a thermal population.

### 3. STUDY OF THE MULTI-SCALE $\alpha$-PARTICLE COLLISION OPERATOR

The aim of this section is to investigate from a more mathematical point of view the one-dimensional $v$-part of the collision operator (2.11) (corresponding to the drag and energy scattering), in particular to underline its multiscale nature, study its mathematical properties and prepare thus the design of an adequate multiscale discretization strategy. The investigation of the multi-dimensional problem, with the inclusion of the Lorentz scattering operator $\mathcal{L}(f_{\alpha})$ (2.10), will be the aim of a forthcoming paper.

3.1. Multiscale $\alpha$-particle slowing-down. Goal of this subsection is to investigate in more details the influence of each part of the one-dimensional $\alpha$-particle collision operators, $Q_{\alpha e}$ resp. $Q_{\alpha i}$, on the whole particle dynamics. Starting from (2.2)-(2.3) with $f_s = M_{n_s,0,T_s}$ leads in one dimension to the operator

$$ (3.1) \qquad Q_{\alpha s}(v) = 2 \nu_0^{\alpha s} \partial_v \left[ \frac{G(v/v_{th,s})}{v/v_{th,s}} \left( \frac{m_\alpha}{m_s} v f_\alpha + \frac{v_{th,s}^2}{2} \partial_v f_\alpha \right) \right], \quad \forall v \in \mathbb{R}, $$
where $\nu_0^{\alpha s} := c_{\alpha s} \frac{n_s}{v_{th,s}^3}$ and $c_{\alpha s} := 4 \pi \left( \frac{e_\alpha e_s}{m_\alpha \varepsilon_0} \right)^2 \ln \Lambda_{\alpha s}$. Let us introduce as previously the collision frequency $\bar{\nu}_0 := c_{\alpha s} \frac{n_\alpha}{v_{\alpha,\star}^3}$, depending only on the $\alpha$-particle characteristics.

To compute the dynamics associated with the $\alpha$-particles, we consider the following equation

$$ \partial_t f_\alpha = Q_{\alpha e} + Q_{\alpha i}, \quad \forall (t,v) \in \mathbb{R}^+ \times \mathbb{R}, $$

forgetting for the moment about the transport term. Under an appropriate scaling (discussed in subsection 2.2) in the velocity variable, namely $v = v_{\alpha,\star} v'$ with $v_{th,i} \ll v_{\alpha,\star} \ll v_{th,e}$, and rescaling also the time variable as $t = \bar{t} t'$, yields

$$ (3.2) \qquad \partial_t f_\alpha = 2 \underbrace{\frac{\bar{\nu}_0 \bar{t}}{\varepsilon \beta} \partial_v \left[ \frac{G(\varepsilon v)}{\varepsilon v} \left( v f_\alpha + \frac{\varepsilon^2}{2} \partial_v f_\alpha \right) \right]}_{Q_{\alpha e}} + 2 \underbrace{\frac{\bar{\nu}_0 \bar{t}}{\varepsilon^3 \beta} \partial_v \left[ \frac{G(v/\varepsilon)}{v/\varepsilon} \left( v f_\alpha + \frac{\varepsilon^2}{2} \partial_v f_\alpha \right) \right]}_{Q_{\alpha i}}. $$

Defining now for simplicity reasons $\eta := 2 \frac{\bar{\nu}_0 \bar{t}}{\varepsilon^2 \beta}$ ($\varepsilon$-dependent), we can rewrite equation (3.2) as

$$ (3.3) \qquad \partial_t f_\alpha = \eta \partial_v \left[ D_\varepsilon(v) \left( v f_\alpha + \frac{\varepsilon^2}{2} \partial_v f_\alpha \right) \right], \quad D_\varepsilon(v) := \frac{1}{v} [G(\varepsilon v) + G(v/\varepsilon)], \quad \forall v \in \mathbb{R}. $$

The goal of the next subsection will be to analyze this equation from a mathematical point of view. However, to get a first insight into its main features we plotted on Figure 3 the different rescaled Chandrasekhar functions $G(v)$, $G(\varepsilon v)$ and $G(v/\varepsilon)$ as well as the corresponding diffusion coefficient $D_\varepsilon(v)$ for $\varepsilon = 0.1$. What one can remark from these Figures is the following:

(a) the diffusion coefficient is positive, strictly decreasing, with $D_\varepsilon(0) = \frac{2}{3\sqrt{\pi}}(\varepsilon + \frac{1}{\varepsilon})$, which is approximately equal to 3.8 for $\varepsilon = 0.1$, and $D_\varepsilon(v) \to 0$ as $v \to \infty$;

(b) the two functions $G(\varepsilon v)$ and $G(v/\varepsilon)$ are nothing else than a rescaling of the original Chandrasekhar function $G(v)$, and they intersect at some value $v_c = \left( \frac{3\sqrt{\pi}}{4} \varepsilon \right)^{1/3} \approx 0.51$ for $\varepsilon = 0.1$.

<figure><img src="image_5.png" /><figcaption>FIGURE 3. Left: Plot of the Chandrasekhar functions $G(v)$, $G(\varepsilon v)$ and $G(v/\varepsilon)$ for $\varepsilon = 0.1$. Right: Plot of the corresponding diffusion coefficient $D_\varepsilon(v)$.</figcaption></figure>

From these two Remarks and Figures one can now understand some features of the kinetic equation (3.3). Firstly, the two functions $G(\varepsilon v)$ and $G(v/\varepsilon)$ act in different velocity regions, in particular $G(\varepsilon v)$, which comes from the electron/$\alpha$-particle collision operator $Q_{\alpha e}$, acts
at large velocities, thus in the initial phase of the α-particle injection, and $G(v/ε)$, which comes from the ion/α-particle collision operator $Q_{αi}$ acts at smaller velocities, thus after the α-particles have been slowed-down. The influence of these two operators becomes equal at the crossing velocity $v_c$ approximately. Secondly, the diffusion coefficient $D_ε(v)$ has a very fine boundary layer at small velocities, which gets finer and finer with smaller $ε$, as one can observe comparing Figure 3 ($ε = 0.1$) with Figure 4 ($ε = 0.05$). Thus, in order to precisely take into account for this diffusion coefficient, one has to adopt a grid which is finer in the small velocity ranges and larger in the rest of the velocity domain, leading naturally to a choice of an adaptive grid or a domain-decomposition method (see Section 4).

<figure><img src="image_6.png" /><figcaption>FIGURE 4. Plot of the Chandrasekhar functions $G(v)$, $G(εv)$, $G(v/ε)$ as well as the corresponding diffusion coefficient $D_ε(v)$ for $ε = 0.05$.</figcaption></figure>

## 3.2. Mathematical study.
Let us come now to the more theoretical investigations of the kinetic equation

$$ (3.4) \qquad \begin{cases} \partial_t f_\alpha = \eta \partial_v \left[ D_\epsilon(v) \left( v f_\alpha + \frac{\epsilon^2}{2} \partial_v f_\alpha \right) \right], & D_\epsilon(v) := \frac{1}{v} [G(\epsilon v) + G(v/\epsilon)], \\ f_\alpha(0, v) = f_{\alpha,in}(v), & \forall v \in \mathbb{R}, \end{cases} $$

for fixed, independent parameters $ε \in (0, 1)$ and $η \in \mathbb{R}^+$. We are interested in particular in the existence and uniqueness of a solution to (3.4), which is a standard result based on Lions' theorem [36]. Furthermore, we aim to investigate the time decay of this solution towards the corresponding stationary state $f_{\alpha,\infty}(v) = \frac{n_\alpha}{\sqrt{\pi} \epsilon} e^{-\frac{v^2}{\epsilon^2}}$, with $n_\alpha := \int_{\mathbb{R}} f_{\alpha,in}(v) dv$ obtained from the conservation of mass property, namely $\partial_t \int_{\mathbb{R}} f_\alpha(t, v) dv = 0$.

**Remark 1.** The factor $η$ in front of the collision operator can be eliminated with a simple change of variable in time. This shall be done in this more theoretical section (setting hence $η = 1$), however for the physical and numerical part, we shall further keep this factor to accelerate the simulations. At the same time, the $ε$-dependence of the diffusion coefficient $D_ε(v)$ is very important and shall be kept in the whole paper.

For all this theoretical investigations an adequate functional framework is required, and we shall adopt for theoretical reasons the functional transformation $f_\alpha = h f_{eq}$, where $f_{eq}(v) := \frac{1}{\sqrt{\pi} \epsilon} e^{-\frac{v^2}{\epsilon^2}}$ is an equilibrium of the Fokker-Planck equation. This leads finally to the following
evolution problem for the unknown distribution function $h(t, v)$

$$ (3.5) \qquad \left\{ \begin{array}{ll} \partial_t h = \frac{1}{f_{eq}(v)} \partial_v [D_\varepsilon(v) f_{eq}(v) \partial_v h], & \forall (t,v) \in \mathbb{R}^+ \times \mathbb{R}, \\ h(0,v) = h_{in}(v), & \forall v \in \mathbb{R}, \end{array} \right. $$

associated with an initial condition $h_{in}$ and with the corresponding stationary state $h_\infty \equiv \bar{h}_{in} = \int_{\mathbb{R}} h_{in}(v) f_{eq}(v) dv$, which is obtained from the fact that $\partial_t \bar{h} = \partial_t \int_{\mathbb{R}} h(t, v) f_{eq}(v) dv = 0$.

Standard arguments, such as Hille Yosida's theorem respectively Lions' theorem, permit to show the following existence and uniqueness result of a solution to the Fokker-Planck equation (3.5).

**Proposition 1. (Existence/uniqueness of a solution)** Let us denote the second order differential operator occurring in (3.5) by

$$ \mathcal{L}_{D,eq}(h) := -\frac{1}{f_{eq}} \partial_v [D_\varepsilon f_{eq} \partial_v h], \quad D_\varepsilon(v) := \frac{1}{v} [G(\varepsilon v) + G(v/\varepsilon)]. $$

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

$$ (3.6) \quad L_{eq}^2 := \left\{ h : \mathbb{R} \to \mathbb{R} \text{ measurable}, \quad \int_{\mathbb{R}} |h|^2 f_{eq} \, dv < \infty \right\}, \quad (h, g)_{L_{eq}^2} := \int_{\mathbb{R}} h g f_{eq} \, dv, $$

with definition domain

$$ (3.7) \qquad \mathcal{D}(\mathcal{L}_{D,eq}) = \{h \in L_{eq}^2, \mathcal{L}_{D,eq}(h) \in L_{eq}^2\}. $$

Furthermore, for each $h_{in} \in \mathcal{D}(\mathcal{L}_{D,eq})$ there exists a unique solution $h \in C^1([0, \infty); L_{eq}^2) \cap C([0, \infty); \mathcal{D}(\mathcal{L}_{D,eq}))$ of the Fokker-Planck equation (3.5). For less regular initial conditions $h_{in} \in L_{eq}^2$ and arbitrary $T > 0$, one has nonetheless a unique weak solution $h \in W_2^1(0, T; H_{eq}^1, L_{eq}^2) \subset C([0, T]; L_{eq}^2)$, where

$$ (3.8) \quad H_{D,eq}^1 := \{h \in L_{eq}^2 / \sqrt{D_\varepsilon(v)} \partial_v h \in L_{eq}^2\}, \quad \|h\|_{H_{D,eq}^1}^2 = \|h\|_{L_{eq}^2}^2 + \|\sqrt{D_\varepsilon} \partial_v h\|_{L_{eq}^2}^2. $$

Finally one has also a maximum principle, meaning that if furthermore $h_{in} \in L^\infty(\mathbb{R}_v)$, then

$$ \|h(t, \cdot)\|_{L^\infty(\mathbb{R}_v)} \le \|h_{in}\|_{L^\infty(\mathbb{R}_v)}, \quad \forall t \ge 0. $$

Let us pass now to the decay of the solution $h$ of (3.5) towards the stationary solution $h_\infty$.

**Theorem 1. (Time decay)** Let $h$ be a solution of the evolution problem (3.5), with an initial condition $h_{in} \in L^\infty(\mathbb{R}_v) \subset L_{eq}^2$. Then the following estimate holds

$$ (3.9) \qquad \|h(t) - h_\infty\|_{L_{eq}^2}^2 \le \left[ \|h_{in} - \bar{h}_{in}\|_{L_{eq}^2}^{-2/r} + K_{\varepsilon,r} \frac{t}{r} \right]^{-r}, \quad \forall t \ge 0, \forall r > 0, $$

where $h_\infty = \bar{h} = \int_{\mathbb{R}} h f_{eq} \, dv = \bar{h}_{in}$ and $K_{\varepsilon,r} > 0$ is a constant explicit in the proof.

We have only a "super-algebraic decay" rate, and not an exponential one as in the standard Fokker-Planck case, because although the functional space we are working with enjoys a Poincaré-Wirtinger inequality, the operator on the right-hand side of (3.5) features a vanishing diffusion coefficient, such that we need a Nash-type inequality. The proof of this theorem is based on the following Poincaré inequality.
**Lemma 1. (Poincaré-Wirtinger inequality)** Let $D_\varepsilon(v)$ be the diffusion coefficient given in Proposition 1. Then, the following Poincaré inequality holds

$$ (3.10) \quad \int_{\mathbb{R}} |h - \tilde{h}|^2 D_{\varepsilon}(v) f_{eq} \, dv \le C_{P,\varepsilon} \int_{\mathbb{R}} |\partial_v h(v)|^2 D_{\varepsilon}(v) f_{eq} \, dv, \quad \forall h \in H^1(\mathbb{R}; f_{eq} \, dv), $$

where $\tilde{h} := \int_{\mathbb{R}} h D_{\varepsilon}(v) f_{eq} \, dv$ and $C_{P,\varepsilon} > 0$ is the Poincaré constant, dependent on $\varepsilon \in (0,1)$.

*Proof.* The proof of this Lemma is based on Corollary 1.6 of [4], with $V(v) := \frac{v^2}{\varepsilon^2} + \ln(\sqrt{\pi} \varepsilon) - \ln(D_\varepsilon(v))$. $\square$

*Proof of Theorem 1.* To prove the time-decay of the solution, let us integrate equation (3.5) against $h(t, v) f_{eq} \, dv$, to get

$$ (3.11) \quad \frac{1}{2} \frac{d}{dt} \|h(t) - \bar{h}\|_{L_{\varepsilon q}^2}^2 = - \int_{\mathbb{R}} D_{\varepsilon}(v) (\partial_v h)^2 f_{eq} \, dv \le - \frac{1}{C_{P,\varepsilon}} \int_{\mathbb{R}} |h - \tilde{h}|^2 D_{\varepsilon}(v) f_{eq} \, dv, $$

where we used the Poincaré inequality (3.10). Concluding at this stage is problematic, due to the occurrence of the diffusion term $D_\varepsilon(v)$ on the right-hand side. To overcome this difficulty, let us proceed as follows

$$ \begin{aligned} \int_{\mathbb{R}} |h(t) - \bar{h}|^2 f_{eq} \, dv &= \inf_{c \in \mathbb{R}} \int_{\mathbb{R}} |h - c|^2 f_{eq} \, dv, & \bar{h} &:= \int_{\mathbb{R}} h f_{eq} \, dv \\ &\le \int_{\mathbb{R}} |h - \tilde{h}|^2 f_{eq} \, dv, & \tilde{h}(t) &:= \int_{\mathbb{R}} h(t) D_{\varepsilon}(v) f_{eq} \, dv \\ &\stackrel{\text{Hölder}}{\le} \left( \int_{\mathbb{R}} |h - \tilde{h}|^2 D_{\varepsilon}(v) f_{eq} \, dv \right)^{\frac{1}{p}} \left( \int_{\mathbb{R}} |h - \tilde{h}|^2 D_{\varepsilon}(v)^{-\frac{1}{p-1}} f_{eq} \, dv \right)^{\frac{p-1}{p}}. \end{aligned} $$

Now, since the Maxwellian has bounded moments, we have for all $p > 1$

$$ (3.12) \qquad M_{\varepsilon,p} := \int_{\mathbb{R}} D_{\varepsilon}(v)^{-\frac{1}{p-1}} f_{eq} \, dv < \infty. $$

And secondly, the maximum principle permits to show that

$$ |h - \tilde{h}|^2(t, v) \le ||h_{in}||_{L^{\infty}(\mathbb{R}_v)}^2 \left(1 + \int_{\mathbb{R}} D_{\varepsilon}(v) f_{eq} \, dv\right)^2, \quad \forall (t, v) \in \mathbb{R}^+ \times \mathbb{R}_v, $$

such that one gets

$$ \int_{\mathbb{R}} |h(t) - \bar{h}|^2 f_{eq} \, dv \le \left( \int_{\mathbb{R}} |h - \tilde{h}|^2 D_{\varepsilon}(v) f_{eq} \, dv \right)^{\frac{1}{p}} \mathcal{R}_{\varepsilon,p}^{\frac{p-1}{p}}, $$

where

$$ \mathcal{R}_{\varepsilon,p} := M_{\varepsilon,p} \left( 1 + \int_{\mathbb{R}} D_{\varepsilon}(v) f_{eq} \, dv \right)^2 ||h_{in}||_{\infty}^2 < \infty. $$

Introducing the constant $K_{\varepsilon,p} := \frac{2}{C_{P,\varepsilon} \mathcal{R}_{\varepsilon,p}^{p-1}}$, permits thus to obtain altogether

$$ \frac{d}{dt} \|h - \bar{h}\|_{L_{\varepsilon q}^2}^2 \le -K_{\varepsilon,p} \left( \int_{\mathbb{R}} |h - \bar{h}|^2 f_{eq} \, dv \right)^p, \quad \forall p > 1. $$

We can conclude finally at this stage with Gronwall's lemma, leading for all $p > 1$ to

$$ \|h(t) - h_{\infty}\|_{L_{\varepsilon q}^2}^2 \le \left[ \|h_{in} - \bar{h}_{in}\|_{L_{\varepsilon q}^2}^{-2(p-1)} + (p-1) K_{\varepsilon,p} t \right]^{-\frac{1}{p-1}}, \quad \forall t \ge 0. $$

Setting then $r := \frac{1}{p-1} \in (0, \infty)$ yields (3.9).

$\square$
## 4. MULTI-SCALE NUMERICAL SCHEME

Let us say some words about the difficulties encountered when trying to solve the multiscale problem (3.4) with a standard method. Concerning the velocity grid, we mentioned earlier that a non-homogeneous grid will be necessary in order to get an accurate resolution of the stiff gradients of the diffusion coefficient $D_\epsilon(v)$, arising in the small velocity region $v \in (0, \epsilon^{1/3})$. What about the time-discretization? An explicit discretization is linked with a CFL-type (or parabolic-type) stability condition, which restricts the time steps $\Delta t > 0$. This one will be driven by the small cells in $v$, localized in the velocity region $v \in (0, \epsilon^{1/3})$, and will lead immediately, by forcing an overall small time-step, to very long simulation times. One way to overcome this problem is a domain decomposition strategy, separating the velocity domain in the boundary layer region $v \in (0, \epsilon^{1/3})$ (fine grid in $(t, v)$) and the outer region $v \in (\epsilon^{1/3}, 1)$ (coarse grid in $(t, v)$). The coupling between both regions can be done either via adequate boundary conditions (Dirichlet-to-Neumann for example) or by separating the $\alpha$-particle population in supra-thermal and thermal particles, evolving separately on the two distinct velocity regions (with distinct discretizations), and permitting an exchange between both populations. In this paper we shall remain in one dimension in the velocity space, and shall rather use an implicit scheme for the resolution of the Fokker-Planck kinetic equation, which is unconditionally stable, thus not requiring any restriction on the time-step. Implicit schemes require however the resolution of a linear system, which is not so drastic in one dimension. When dealing however with the multi-dimensional case, including also the Lorentz collision operator, more performant schemes are required, based for example on the above mentioned ideas of decomposition, and will be the aim of some forthcoming papers.

### 4.1. Scaling of the Fokker-Planck kinetic equation.
Let us consider in the following the full $\alpha$-particle Fokker-Planck kinetic equation

$$ \partial_t f_\alpha + v \partial_x f_\alpha + \frac{e_\alpha E}{m_\alpha} \partial_v f_\alpha = Q_{\alpha i} + Q_{\alpha e}, $$

with $Q_{\alpha s}$ defined in (3.1) (or rescaled version (3.2)). To propose a scaling of this equation, we introduce characteristic values for all the occurring quantities and perform a change of variables, which leads to

$$ \partial_t f_\alpha + \frac{\bar{v} \bar{t}}{\bar{x}} v \partial_x f_\alpha + \frac{e_\alpha \bar{E} \bar{t}}{m_\alpha \bar{v}} E \partial_v f_\alpha = \frac{\bar{Q} \bar{t}}{\bar{f}} Q(f). $$

The characteristic scales are summarized here:

* Microscopic velocity scale (Parameter: $\epsilon$)

$$ \bar{v} := v_{\alpha, \star}, \quad \epsilon := \frac{v_{th,i}}{v_{\alpha, \star}} = \frac{v_{\alpha, \star}}{v_{th,e}} \ll 1, \quad v_{th,s} := \sqrt{\frac{2 k_B \bar{T}_s}{m_s}}; $$

* Disparate masses (Parameter: $\epsilon$) and same bulk temperatures

$$ \epsilon^4 = \frac{m_e}{m_i} \ll 1, \quad \bar{T}_i = \bar{T}_e. $$

* Macroscopic time and length scale (Parameter: $\alpha$)

$$ \bar{x}, \bar{t}: \text{ observation time and distance, } \alpha := \frac{\bar{v} \bar{t}}{\bar{x}}; $$
* Electric energy versus kinetic energy (Parameter: $\gamma$)

$$ \gamma := \frac{e \bar{\phi}}{m_{\alpha} \bar{v}^2}; $$

* Energetic versus thermal particle densities (Parameter: $\beta$)

$$ \beta := \frac{\bar{n}_{\alpha}}{\bar{n}_{s}}, \quad \bar{n}_{i} = \bar{n}_{e}; $$

* Collisional frequency (Parameter: $\eta$)

$$ \eta := 2 \frac{\bar{\nu}_0 \bar{t}}{\varepsilon^2 \beta}, \quad \bar{Q} = \bar{\nu}_0 \bar{f}_{\alpha}. $$

Introducing these parameters permits to obtain the rescaled kinetic equation (4.1)

$$ \partial_t f_{\alpha} + \frac{1}{\alpha} v \partial_x f_{\alpha} + \frac{\gamma}{\alpha} E v f_{\alpha} = \eta \partial_v \left[ D_{\varepsilon}(v) \left( v f_{\alpha} + \frac{\varepsilon^2}{2} \partial_v f_{\alpha} \right) \right], \quad D_{\varepsilon}(v) := \frac{1}{v} [G(\varepsilon v) + G(v/\varepsilon)]. $$

Several asymptotic regimes can now be investigated. As we are interested in the present work in an efficient discretization of the multi-scale $\alpha$-particle collision operator, we shall set $\alpha = \gamma = 1$ and keep only $\varepsilon, \eta \in (0, 1)$ as parameters.

### 4.2. Velocity discretization.
To propose a numerical scheme for (4.1), we shall start by considering first the following homogeneous equation

$$ (4.2) \qquad \left\{ \begin{array}{ll} \partial_t f = \eta \partial_v \left[ D_{\varepsilon}(v) \left( v f + \frac{\varepsilon^2}{2} \partial_v f \right) \right], & \forall (t, v) \in (0, t_{\max}) \times (-v_{\max}, v_{\max}), \\ f(0, v) = f_{in}(v), & \forall v \in [-v_{\max}, v_{\max}], \end{array} \right. $$

associated with homogeneous Dirichlet boundary conditions. For the discretization of this evolution problem a non-uniform grid is required in the velocity variable. Indeed, on one hand, in order to capture correctly the stationary state $f_{\infty} = \frac{n}{\sqrt{\pi \varepsilon}} e^{-\frac{v^2}{\varepsilon^2}}$, $n = \int_{-\infty}^{\infty} f_{in} dv$, a very fine mesh is required close to $v = 0$, namely $\Delta v \approx \varepsilon$. On the other hand, to take into account for the rapid $\alpha$-particle dynamics a large velocity domain $[-v_{\max}, v_{\max}]$ is mandatory. Considering both constraints necessarily leads, if a uniform velocity grid is employed, to a very large number of points.

For this reason, we shall rather consider a non-uniform velocity mesh $\{v_j\}_{j=0}^{N_v}$ and define the mid-points $v_{j+1/2} := (v_j + v_{j+1})/2$ on which equation (4.2) will be approximated via a finite volume scheme. Integrating (4.2) 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) &= \eta \frac{1}{\Delta v_j} \left[ D_{\varepsilon}(v_{j+1/2}) \left( v_{j+1/2} f(t, v_{j+1/2}) + \frac{\varepsilon^2}{2} (\partial_v f)(t, v_{j+1/2}) \right) \right. \nonumber \\
&\qquad \left. - D_{\varepsilon}(v_{j-1/2}) \left( v_{j-1/2} f(t, v_{j-1/2}) + \frac{\varepsilon^2}{2} (\partial_v f)(t, v_{j-1/2}) \right) \right] \nonumber \\
&\approx \eta \frac{1}{\Delta v_j} \left[ D_{\varepsilon}(v_{j+1/2}) \left( v_{j+1/2} \frac{f_j(t) + f_{j+1}(t)}{2} + \frac{\varepsilon^2}{2} \frac{f_{j+1}(t) - f_j(t)}{v_{j+1} - v_j} \right) \right. \nonumber \\
&\qquad \left. - D_{\varepsilon}(v_{j-1/2}) \left( v_{j-1/2} \frac{f_j(t) + f_{j-1}(t)}{2} + \frac{\varepsilon^2}{2} \frac{f_j(t) - f_{j-1}(t)}{v_j - v_{j-1}} \right) \right] \nonumber \\
&= \eta \frac{1}{\Delta v_j} \left[ f_{j-1}(t) D_{\varepsilon}(v_{j-1/2}) \left( -\frac{v_{j-1/2}}{2} + \frac{\varepsilon^2}{2(v_j - v_{j-1})} \right) \right. \nonumber \\
&\qquad \left. + f_j(t) \left( \frac{D_{\varepsilon}(v_{j+1/2}) v_{j+1/2} - D_{\varepsilon}(v_{j-1/2}) v_{j-1/2}}{2} - \frac{\varepsilon^2}{2} \left( \frac{D_{\varepsilon}(v_{j-1/2})}{v_j - v_{j-1}} + \frac{D_{\varepsilon}(v_{j+1/2})}{v_{j+1} - v_j} \right) \right) \right. \nonumber \\
&\qquad \left. + f_{j+1}(t) D_{\varepsilon}(v_{j+1/2}) \left( \frac{v_{j+1/2}}{2} + \frac{\varepsilon^2}{2(v_{j+1} - v_j)} \right) \right], \quad \forall j = 1, \dots, N_v - 1, \tag{4.3}
\end{align}
$$

whereas we set for the boundary values $f_0 = f_{N_v} \equiv 0$. Having discretized the equation in the velocity variable, one has still to generate the non-uniform mesh. The first criteria for the mesh is to have a velocity step size of order $\varepsilon$ close to the origin, in order to capture the asymptotic state $f_\infty$. Our strategy for constructing a non-uniform velocity mesh is the following: for a given $v_{\max}$, large enough to take into account for the support of the distribution function, one defines firstly $w_{\max} := v_{\max}^{1/\text{pow}}$ and a uniform mesh $\{w_j\}_{j=0}^{N_v}$ as $w_j := -w_{\max} + j\Delta w$, with $\Delta w := 2w_{\max}/N_v$ and we deduce then the non-uniform mesh $\{v_j\}_{j=0}^{N_v}$ from the relation $v_j := w_j^{\text{pow}}$ for all indices $j$. This simple non-uniform mesh construction permits to obtain refined cells close to the origin $v = 0$, since $v_j = O(\Delta v^{\text{pow}})$ whereas far from the origin, $\Delta v = O(1)$. In practice, we choose pow=2 or pow=3. Figure 5 illustrates the mesh obtained with a cube rule (pow=3).

<figure><img src="image_7.png" /><figcaption>FIGURE 5. Example of a non-uniform mesh with a 'cube' rule (pow=3). Here $v_{\max} = 10$.</figcaption></figure>

4.3. **Time discretisation.** In what concerns the time discretization, one has to care about the fact that, due to a very refined velocity mesh close to the origin $v = 0$, any explicit
time integrator will suffer from a very stringent CFL parabolic condition $\Delta t = O(\Delta v^2)$, with $\Delta t > 0$ the time step. Moreover, the presence of the factor $\eta$ (which may be very large) in front of the collision operator leads to $\Delta t = O(1/\eta)$. Hence, we shall prefer an implicit time discretization which requires however the solution of a linear system induced by (4.3). Introducing the time discretization $t^n = n\Delta t, n \in \mathbb{N}$, the numerical unknown is denoted by $f_j^n \approx \frac{1}{\Delta v_j} \int_{v_{j-1/2}}^{v_{j+1/2}} f(t^n, v)dv$ and the vector $f^n \in \mathbb{R}^{N_v-1}$ (with components $\{f_j^n\}_{j=1}^{N_v-1}$ and Dirichlet boundary conditions $f_0^n = f_{N_v}^n = 0$) satisfies the following linear system $Af^{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{\eta}{\Delta v_j} D_{\varepsilon}(v_{j-1/2}) \left( -\frac{v_{j-1/2}}{2} + \frac{\varepsilon^2}{2(v_j - v_{j-1})} \right),$$

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

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

### 4.4. Transport approximation.
When one is interested in the coupling between the transport and collision operator, splitting methods are very popular so that dedicated methods can be used for each specific operator (transport or diffusion). In this work we shall use a 3-step splitting procedure, consisting of:

*   **Advection in the x-direction:** $\partial_t f + v \partial_x f = 0$. Here a Fourier spectral method is employed since we are supposed to be periodic in the space-variable $x$, the equation becomes then an ODE in the Fourier space. Denoting by $\hat{f}_k(t, v)$ the $k$-mode of $f(t, x, v)$, the solution of this transport equation writes $\hat{f}_k(t, v) = e^{-ikv(t-t^n)} \hat{f}_k(t^n, v)$.

*   **Poisson equation:** $\partial_x E(t, x) = \int_{\mathbb{R}} f(t, x, v)dv - 1$. The electric field is approximated equally using Fourier spectral method and we impose zero space-average, to fix the constant.

*   **Advection in the v-direction:** $\partial_t f + E \partial_v f = 0$. For this step, let us recall that a non-uniform velocity mesh has to be taken into account. For a fixed position $x$, we have to deal with a linear transport equation since $E$ does not evolve during one time-step. Hence, we use a semi-Lagrangian method, approach based on the fact that the function $f$ is constant along the characteristics: $f(t, x, v) = f(t^n, x, v - (t - t^n)E(t^n, x))$. Since the foot $v - (t - t^n)E$ does not coincide with a grid point $v_j$, an interpolation is performed. The update of $f$ can then be written as $f(t, x, v) = \mathcal{I}[f(t^n, x, \cdot)](v - (t - t^n)E(t^n, x))$ with $\mathcal{I}$ an interpolation operator (a cubic spline interpolation is used). Let us remark that as we are using a non-uniform mesh, the localisation of the foot of the characteristics $v - (t - t^n)E(t^n, x)$ is not as straightforward as in the uniform mesh case. We refer to [1, 6] for related approaches.

*   **Collision part:** $\partial_t f = \eta Q^\varepsilon(f)$. For this part, a time-implicit method, combined with a finite volume method on a non-uniform velocity mesh, is used, as described in the two previous subsections.

These steps are solved successively to get a first order in time approximation of the collisional Vlasov-Poisson equation. Of course, a Strang splitting extension enables to get second order accuracy. The so-obtained scheme, which is halfway between an Eulerian and a Lagrangian method, does not require any CFL condition due to the implicit treatment of the collisional
part and the semi-Lagrangian method for the transport parts. It allows larger time-steps $\Delta t > 0$ compared to Eulerian codes, while avoiding numerical noise, standard in Lagrangian methods.

## 5. NUMERICAL RESULTS AND INTERPRETATIONS

In this section we shall present some numerical results for the resolution of the Fokker-Planck Vlasov-Poisson problem in different physical situations, results obtained with the previously introduced numerical scheme. Firstly, a study in the spatially homogeneous case is provided. Then, some tests will be given, including at the same time the collision operator and the free flow transport term, before investigating the fully coupled Fokker-Planck-Poisson problem. Aim is to validate the numerical scheme, in particular the discretization of the Fokker-Planck collision operator, to observe the different relaxation mechanisms in the energetic $\alpha$-particle dynamics (slowing-down as a consequence of the collisions with the thermal bulk) and to verify the mathematical time-decay Theorem 1, predicting a super-algebraic decay of the solution towards the steady state $f_\infty(v) = \frac{n}{\sqrt{\pi\epsilon}}e^{-v^2/\epsilon^2}$.

5.1. **Spatially homogeneous case.** Let us illustrate in this subsection the dynamics of the solution of (4.2), obtained with the numerical scheme presented in Section 4. The initial condition is chosen as

$$f_{in}(v) := \frac{1}{\sqrt{2\pi T}} e^{-\frac{(v-u)^2}{2T}}, \quad \text{with } u = 5, T = 0.2 \text{ and } \forall v \in (-v_{max}, v_{max}).$$

This initial condition describes $\alpha$-particles generated by fusion reactions at a certain temperature and with a certain mean velocity, our goal being to study their slowing-down and relaxation through the collisions with the thermal bulk. In our context, the dynamics of the $\alpha$-particles is supposed to converge in the long-time limit $t \to \infty$ towards the asymptotic state given by $f_\infty(v) = \frac{n}{\sqrt{\pi\epsilon}}e^{-v^2/\epsilon^2}$, with the constant density $n = \int_{-\infty}^\infty f_{in}(v) dv \equiv 1$.

First, we consider $\epsilon = 0.1$ (which is close to the relevant physical value) and $\eta = 1/\epsilon^2$ in (4.2), with a mesh generated from $N_v = 300$, $v_{max} = 10$ and pow= 3, whereas the time step is chosen as $\Delta t = 10^{-2}$ (the final time being $t_{max} = 2$). The results are given in Figure 6: on the left plot are represented the $\alpha$-particle velocity distributions $f(t, v)$ for different times, whereas on the right plot we displayed the difference between the final numerical solution (at $t_{max} = 2$) and the stationary state $f_\infty(v)$. We can observe the two stages in the dynamics: firstly a slowing-down of the $\alpha$-particles can be observed between $t = 0$ and $t = 0.4$ (due to collisions with the light bulk electrons) and secondly a thermalization occurs in the long-time limit $t \to \infty$ (due to collisions with the heavy bulk ions). After $t = 1$, the solution becomes very close to the stationary state $f_\infty$. Finally, let us remark that the discrete total mass, defined as $m^n = \sum_{j=1}^{N_v-1} (v_{j+1/2} - v_{j-1/2})f_j^n$, is preserved with machine accuracy.

The same test is performed with a smaller value of $\epsilon = 0.01$ ($\eta = 1/\epsilon^2$) and the results are displayed on Figure 7. The two scales in the $\alpha$-particle dynamics can be equally observed here (slowing-down and thermalization), but at a much faster rate, due to the presence of a larger $\eta = 1/\epsilon^2$ in front of the collision operator $Q^\epsilon$. Comparing the two absolute-error plots on the right of Figures 6 and 7 permits to remark that the error is more localized for smaller $\epsilon$-values, however it is larger due to the fact that with smaller $\epsilon$-values the steady-state $f_\infty(v) = \frac{n}{\sqrt{\pi\epsilon}}e^{-v^2/\epsilon^2}$ becomes higher, the relative error being nevertheless similar for both $\epsilon$-values.
The goal of a second test case was to illustrate the time-decay Theorem 1, stated in subsection 3.2, by computing the quantity $H(t) := \|h(t)-h_\infty\|_{L^2_{\varepsilon,\infty}} = (\int_{\mathbb{R}} |f(t,v)-f_\infty(v)|^2 f_{eq}^{-1}(v) dv)^{1/2}$. In particular, we are interested in the $\varepsilon$-dependence of the constant $K_{\varepsilon,r} > 0$, which was not explicited in the theoretical part. To do so, we run a simulation with a close to equilibrium initial condition

$$f_{in}(v) := \frac{1}{\sqrt{\pi} \varepsilon} e^{-\frac{(v-u)^2}{\varepsilon^2}}, \quad \text{with} \quad u = 0.01.$$

Indeed, due to the presence of $f_{eq}^{-1} = \varepsilon \sqrt{\pi} e^{v^2/\varepsilon^2}$ in the entropy $H(t)$ and the fact that $\varepsilon \ll 1$, the numerical approximation of $H(t)$ is very delicate and to be able to evaluate it correctly, one has to start close to $f_\infty$. For different values of $\varepsilon = k \times 10^{-1}$, $k = 1, \dots, 10$ and setting here $\eta = 1$, we consider the time evolution of the entropy $H(t)$ in semi-log scale, and compute the decay rate for each value of $\varepsilon$. On the right plot of Figure 8, the time evolution of $H(t)$ (in semi-log scale) is illustrated for $\varepsilon = 0.4$ and one observes clearly an exponential decay towards the equilibrium, whereas the theory predicts a super-algebraic one (Theorem 1). The decay rates (or slopes) are supposed to be an approximation of $K_{\varepsilon,r}$ for large $r$. On the left plot we displayed the decay rates of $H(t)$ for different $\varepsilon$-values. One notices that this constant depends on $\varepsilon$, in particular when $\varepsilon$ becomes small, the constant seems to increase, this being linked to the $\varepsilon$-dependent velocity diffusion-part in the collision operator.

<figure><img src="image_8.png" /><figcaption>FIGURE 6. Case $\varepsilon = 0.1$. Left: Relaxation of $f(t,v)$ towards $f_\infty(v) = \frac{1}{\sqrt{\pi} \varepsilon} e^{-\frac{v^2}{\varepsilon^2}}$. Right: Distance to the stationary distribution function at final time $t_{max} = 2$, i.e. $f(t = 2, v) - f_\infty(v)$. Here $N_v = 200, \Delta t = 0.01$.</figcaption></figure>

5.2. **Including space transport.** Let us include now the free-flow transport term and consider the following equation

$$\left\{ \begin{array}{ll} \partial_t f + v \partial_x f = \eta Q^\varepsilon(f), & \forall (t, x, v) \in (0, t_{\max}) \times (0, L) \times (-v_{\max}, v_{\max}), \\ f(0, x, v) = f_{in}(x, v), & \forall (x, v) \in [0, L] \times [-v_{\max}, v_{\max}], \end{array} \right.$$

associated with periodic boundary conditions in $x$ and homogeneous Dirichlet boundary conditions in $v$. The splitting method briefly described in Section 4 is used for its resolution. The chosen initial condition is

$$f_{in}(x, v) := \frac{1}{\sqrt{2\pi}} e^{-\frac{|v-u|^2}{2}} (1 + 0.5 \cos(kx)), \quad \text{with} \quad k = 0.5, x \in [0, 2\pi/k], u \in \{1.3, 5\},$$
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Case $\epsilon = 0.01$. Left: Relaxation of $f(t, v)$ towards $f_\infty(v) = \frac{1}{\sqrt{\pi}\epsilon}e^{-\frac{v^2}{\epsilon^2}}$. Right: Distance to the stationary distribution function at final time $t_{\max} = 2$, i.e. $f(t = 0.2, v) - f_\infty(v)$. Here $N_v = 500$, $\Delta t = 0.01$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 8. Illustration of the entropy $H(t)$-decay from Theorem 1. Left: exponential decay rates as a function of $\epsilon$. Right: Time evolution (in semi-log scale) of the entropy $H(t)$ for $\epsilon = 0.4$.</figcaption></figure>

which corresponds to a perturbation of a Maxwellian distribution in velocity by a “sinusoidal” density in space. The aim is to investigate the decay towards the stationary state, given by the global Maxwellian $f_\infty(v) := \frac{n}{\sqrt{\pi}\epsilon}e^{-\frac{v^2}{\epsilon^2}}$ with the density $n$ independent on $(t, x)$, fact which is a consequence of the competition between the collision operator $Q^\epsilon$ and the transport operator. The parameter $\epsilon$ as well as the initial $\alpha$-particle mean velocity $u$ are playing here an important role in particular in the determination of the different relaxation time-spans. The collision operator, acting only on the velocity variable, is responsible for the relaxation towards a local Maxwellian, the signification of “local” being that the macroscopic variables are still $(t, x)$-dependent. Simultaneously, the transport operator mixes space and velocity variables, leading in the long-time limit $t \to \infty$ to the global (space-independent) equilibrium. This is the well-known *hypocoercivity* theory [12].

In the following simulations, the number of points in space is $N_x = 64$ (with periodic boundary conditions) whereas $N_v = 512$ points are considered in velocity for the velocity
domain [-9, 9] and with pow=2. We choose a time-step of $\Delta t = 0.5$.

Figures 9 and 10 investigate the local/global relaxation phases for two different initial $\alpha$-particle mean velocities, namely $u = 1.3$ resp. $u = 5$, however for fixed $\varepsilon = 0.1$ and $\eta = 1$. On the left plots of both Figures we plot the space-integrated distribution function $v \rightarrow \int f(t, x, v)dx$ for different times ($t = 0, 20, 100$). As in the spatially homogeneous case (see Figures 6 resp. 7), we observe the two regimes corresponding to collisions with the bulk electrons first and then with the bulk ions, leading to a slowing-down of the $\alpha$-particles. This slowing-down is visible also on the right plots of Figures 9 and 10, the plotted mean velocity is indeed fastly decreasing towards zero with time. On the middle and right plots of these two Figures, we explored the relaxation of the macroscopic quantities (depending on ($t, x$)), meaning we are interested in the relaxation of the distribution function towards the global equilibrium. One can observe clearly (comparing Figures 9 and 10) that a different initial mean velocity can induce different relaxation times, and this is due to the fact that the two relaxations towards a local (microscopic relaxation) resp. a global equilibrium (macroscopic relaxation) are the result of the interplay between the collision operator (acting only on the velocity variable) and the transport term (mixing $x$ and $v$ variables).

At first glance, it seems to be quite surprising to observe that mixing occurs at different time-scales according to the choice of the initial mean velocity $u$. Indeed, for $u = 1.3$ one can see that the solution still enjoys some space dependency at $t = 100$ (Figure 9, middle) whereas for $u = 5$, still at the same time $t = 100$, the solution has already reached the homogeneous state $f_\infty$ (Figure 10, middle). One possible explanation could come from the fact that for a larger $u$, the interaction between the transport and collision terms is stronger so that the global relaxation is faster. For smaller $u$ however, the distribution function fastly relaxes towards a space dependent equilibrium centered around $v = 0$ and with a small temperature $\varepsilon = 0.1$. Hence the transport operator $v\partial_x$ is not very efficient to mix the $x$ and $v$ variables, in order to relax towards the global equilibrium and this global relaxation takes more time. Some tests with larger values of $\varepsilon$ confirm this point. A theoretical study would enable to understand this phenomena more in depth and is the aim of a future work.

Figure 11 studies now the influence of the $\varepsilon$-parameter in the relaxation process towards the global equilibrium, by plotting for two different $\varepsilon$ values ($\varepsilon = 1$ and $0.1$) the evolution $t \rightarrow \|f(t) - f_\infty\|_{L^2_{xv}}$ for $\eta = 1$ and fixed $u = 1.3$. When $\varepsilon = 1$, the distribution function becomes homogeneous in space in short times since at $t \approx 20$ one has $\|f - f_\infty\|_{L^2_{xv}} \approx 10^{-3}$ (see Figure 11, left). Considering a smaller value of $\varepsilon$ (see Figure 11, right, for $\varepsilon = 0.1$) reduces the effect of the velocity-diffusion part in the collision operator $Q^\varepsilon$, so that the space homogenization occurs much later, at times of about $t = 2000$. Indeed, even if the mean velocity $U(t)$ is very small at time $t = 100$ (see Figure 9 (right)), $f$ still enjoys important space variations and becomes close to the homogeneous equilibrium $f_\infty$ later on (around $t = 2000$).

### 5.3. Fokker-Planck-Poisson.
To get more insight into the behaviour of the full non-linear problem, let us consider in this part the following collisional Fokker-Planck-Poisson equation
<figure><img src="image_11.png" /><figcaption>FIGURE 9. Space transport case for $\epsilon = 0.1$, $\eta = 1$ and $u = 1.3$. Left: Time evolution of the space integrated distribution function $v \to \int f(t, x, v)dx$. Middle: Time evolution of the density $n(t, x) = \int f(t, x, v)dv$. Right: Time evolution of the mean velocity $U(t) = \int v f(t, x, v)dv$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 10. Space transport case for $\epsilon = 0.1$, $\eta = 1$ and $u = 5$. Left: Time evolution of the space integrated distribution function $v \to \int f(t, x, v)dx$. Middle: Time evolution of the density $n(t, x) = \int f(t, x, v)dv$. Right: Time evolution of the mean velocity $U(t) = \int v f(t, x, v)dv$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 11. Space transport case. Time evolution of $\|f(t) - f_\infty\|_{L^2}$. Left: $\eta = 1, \epsilon = 1$. Right: $\eta = 1, \epsilon = 0.1$.</figcaption></figure>

satisfied by $(f, E)$, with $E(t, x)$ the self-consistent electric field, i.e.

$$ egin{cases} 
\partial_t f + v \partial_x f + E \partial_v f = \eta Q^\epsilon(f), \quad \forall (t, x, v) \in (0, t_{\max}) \times (0, L) \times (-v_{\max}, v_{\max}), \\ 
\partial_x E = \int_{\mathbb{R}} f dv - 1, \quad \forall (t, x) \in (0, t_{\max}) \times (0, L), \\ 
f(0, x, v) = f_{\text{in}}(x, v), \quad \forall (x, v) \in [0, L]_{23} \times [-v_{\max}, v_{\max}]. 
\end{cases} $$
A 3-stage splitting procedure (described in subsection 5.1) is used to approximate this collisional Vlasov-Poisson equation. In the sequel, we will use the following usual notation for a Maxwellian $M_{n,u,T}$

$$M_{n,u,T}(v) := \frac{n}{\sqrt{2\pi T}} e^{\frac{-|v-u|^2}{2T}}, \quad \forall v \in (-v_{\max}, v_{\max}),$$

and we will be interested in the time evolution of the electric energy through the computation of the $L^2$-norm $\|E(t, \cdot)\|_{L_x^2}$. The goal is here to evaluate the influence of the collision operator on the well-known Landau damping as well as on the bump-on-tail instability.

The first test is a Landau damping test-case for which the initial condition is a small perturbation of a Maxwellian steady state

$$f_{in}(x, v) := M_{1,0,1}(v)(1 + a \cos(kx)), \quad x \in [0, 2\pi/k], \quad a = 10^{-3}, \quad k = 0.5,$$

and whose decay towards a steady state can be interpreted as a stability effect. In Figure 12, we plot the time evolution of the electric energy for $\eta = 0$ (which corresponds to the standard collisionless Vlasov-Poisson equation), for $\eta = 0.1, \varepsilon = 1$ (weakly collisional regime) and for $\eta = 10^3, \varepsilon = 1$, in order to investigate the effect of the collision operator $Q^\varepsilon$ on the Landau damping phenomena. First, when $\eta = 0$, we observe the usual Landau damping with the expected rate of 0.15 obtained by analytical computations [31]. When $\eta = 0.1$, we can observe that the collisional effects are combined with the mixing effects of the transport part, leading to stronger Landau damping rates. Finally, when $\eta = 10^3$, the system is close to the fluid regime resulting in essentially no damping on the electric energy, the reason being that the Landau damping is a pure kinetic effect. This has already been observed for the Vlasov-Poisson-BGK in [14].

We then consider a second test-case from [29], for which the initial condition is of the form of a perturbed two-bump Maxwellian

$$f_{in}(x, v) := (M_{\iota,0,1}(v) + M_{(1-\iota)/2,4,0.2}(v) + M_{(1-\iota)/2,-4,0.2}(v))(1 + a \cos(2\pi x/L)), \quad x \in [0, L].$$

This initial condition is known to generate a two stream instability in the collisionless case. The parameters are chosen as in [29], namely $a = 0.00056, L = 22, \iota = 0.97$. Here, we investigate the effect of the collision operator $Q^\varepsilon$ on the Vlasov-Poisson instability. In Figure 13, we plot the time evolution of the electric energy $\|E(t, \cdot)\|_{L_x^2}$ for $\eta = 0$ (which corresponds to the collisionless Vlasov-Poisson equation), for $\eta = \varepsilon = 0.1$ and for $\eta = 0.1, \varepsilon = 1$. For $\eta = 0$, the instability is observed with the expected growth rate in the linear stage of the evolution, rate of about 0.076 (see [29]). The nonlinear evolution-phase starts at $t \approx 100$, where the amplitude of the electric energy oscillates around a saturation value, fact which corresponds to trapped particles in phase-space vortex structures. For $\eta = 0.1$, we can observe that the linear growth phase of the instability is modified, since the growth rate does not correspond any more to the collisionless one (kept in red in the middle and right Figures 13). For large times, we can even observe that the electric energy decreases towards zero, meaning the non-linear evolution-phase is affected by collisions, these one reducing the trapping vortex structures and pushing the distribution function towards a Maxwellian. This is more clear from Figures 14 where $v \to \int f(t, x, v)dx$ is plotted at different times for $\eta = 0$ (left), $\eta = \varepsilon = 0.1$ (middle) and $\eta = 0.1, \varepsilon = 1$ (right). When $\eta = 0$, one can see trapped particles around $v = \pm 4$ in the final distribution function at $t = 400$ (left of Figure 14). When $\eta = \varepsilon = 0.1$, the trapping does not occur and we observe instead the distribution function converging towards the Maxwellian equilibrium $f_\infty$ (middle of Figure 14). For the
same value of $\eta$ but with $\varepsilon = 1$, the equilibrium is almost reached at $t = 400$ (right of Figure 14). Similar comments are available for Figure 15 where the phase-space contour-plot of the final distribution function is displayed for $\eta = 0$, $\eta = \varepsilon = 0.1$ and $\eta = 0.1, \varepsilon = 1$. The vortices around $v = 4$ can be observed in the collisionless case (typical signature of the particle trapping), whereas Maxwellian profiles are obtained in the collisional case.

<figure><img src="image_14.png" /><figcaption>FIGURE 12. Vlasov-Poisson case: Time evolution of the electric energy $\|E(t)\|_{L^2}$. Left: $\eta = 0$ (no collision). Middle: $\eta = 0.1, \varepsilon = 1$. Right: $\eta = 10^3, \varepsilon = 1$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIGURE 13. Valentini case: Time evolution of the electric energy $\|E(t)\|_{L^2}$. Left: $\eta = 0$ (no collision). Middle: $\eta = \varepsilon = 0.1$. Right: $\eta = 0.1, \varepsilon = 1$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 14. Valentini case: Time evolution of the electric energy $\|E(t)\|_{L^2}$. Left: $\eta = 0$ (no collision). Middle: $\eta = \varepsilon = 0.1$. Right: $\eta = 0.1, \varepsilon = 1$.</figcaption></figure>

## 6. CONCLUDING REMARKS AND PERSPECTIVES

Let us conclude this paper by summarizing what was done, and what remains to be done in some forthcoming papers. Energetic particle collisions with the thermal bulk can have
<figure><img src="image_17.png" /><figcaption>FIGURE 15. Distribution function at the final time $t = 400$ (Valentini case). Left: $\eta = 0$ (no collision). Middle: $\eta = \varepsilon = 0.1$. Right: $\eta = 0.1, \varepsilon = 1$.</figcaption></figure>

significant effects on plasma instabilities and plasma heating, such that an accurate description is crucial for plasma simulations. The main goal of the present paper was to model mathematically in a precise manner energetic $\alpha$-particles generated through nuclear fusion reactions. Several dynamics have been identified in the dynamics of these energetic particles, linked firstly to the collisions with the light thermal electrons, and secondly with the heavy thermal ions. In the long-time limit, after having been slowed-down, the $\alpha$-particles become thermalized, their distribution being thus accurately described by a Maxwellian. The overall particle dynamics being a multi-scale problem, an adequate and efficient numerical simulation of these energetic particle population is essential. We proposed in this first work a simple implicit time-discretization, together with a non-uniform velocity grid, in order to get in contact with the main numerical difficulties. In more dimensions, such scheme could become very time-consuming, such that we shall propose in a future work a more performant method, based on domain-decomposition strategies, taking into account also for the Lorentz collision operator (2.10), describing the pitch-angle diffusion operator (no energy scattering).

**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.

## REFERENCES

[1] B. Afeyan, F. Casas, N. Crouseilles, A. Dodhy, E. Faou, M. Mehrenberger, E. Sonnendrücker *Simulations of kinetic electrostatic electron nonlinear (KEEN) waves with variable velocity resolution grids and high-order time-splitting*, Eur. Phys. J. B, **68** (2014), p. 295.

[2] J. C. Allred, M. Alaoui, A. F. Kowalski, G. S. Kerr, *Modeling the Transport of Nonthermal Particles in Flares Using Fokker-Planck Kinetic Theory*, The Astrophysical Journal **902** (2020), no. 1,

[3] 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.

[4] 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** (2008), pp. 60–66.

[5] 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).
[6] E. Bourne, Y. Munschy, V. Grandgirard, M. Mehrenberger, P. Ghendrih, *Non-uniform splines for semi-Lagrangian kinetic simulations of the plasma sheath*, J. Comput. Phys. **488** (2023), 112229.

[7] S. Briguglio, G. Vlad, F. Zonca, C. Kar, *Hybrid magnetohydrodynamic-gyrokinetic simulation of toroidal Alfvén modes*, Phys. Plasmas 2 (10), 3711–3723 (1995).

[8] J.W. Burby, and C. Tronci, *Variational approach to low-frequency kinetic-MHD in the current coupling scheme*, Plasma Phys. Control. Fusion 59 (4), 045013 (2017).

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

[10] C.Z. Cheng, *A kinetic-magnetohydrodynamic model for low-frequency phenomena*, J. Geophys. Res. Space Phys. 96 (A12), 21159–21171.

[11] A. Chrisment, P. Loiseau, J.-L. Feugeas, P.-E. Masson-Laborde, J. Mathi- aud, V. Tikhonchuk, and Ph. Nicolai, *Analysis of a kinetic model for electron heat transport in inertial confinement fusion plasmas*, Physics of Plasmas 29(6) (2022).

[12] L. Desvillettes, C. Villani *On the trend to global equilibrium for spatially inhomogeneous kinetic systems: The Boltzmann equation*, Inventiones mathematicae **159** (2005), 245–316.

[13] 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 elec-tron simulations*, Journal of Physics: Conference Series **401** (2012).

[14] A. Crestetto, N. Crouseilles, M. Lemou, *Kinetic/Fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equations using particles*, Kin. Rel. Models **5** (2012), pp. 787-816.

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

[16] J. D. Gaffey Jr., *Energetic ion distribution resulting from neutral beam injection in tokamaks*, Journal of Plasma Physics **16** (1976), no. 2, 149–169.

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

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

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

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

[21] F.L. Hinton *Collisional transport in plasma*, Handbook of Plasma Physics, Eds. R.N. Sudan, Netherlands: North-Holland, 1983.

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

[23] W. Manheimer, D. Colombant, and A. Schmitt, *Calculations of nonlocal electron energy transport in laser produced plasmas in one and two dimensions using the velocity dependent Krook model*, Phys. Plasmas 19(5) (2012).

[24] A. Mellet *Diffusion Limit of a Non-Linear Kinetic Model without the Detailed Balance Principle*, Monatsh. Math. **134** (2002), 305–329.

[25] A. Mellet, S. Mischler, C. Mouhot, *Fractional diffusion limit for collisional kinetic equations*, Archive for Rational Mechanics and Analysis **199** (2011), 493–525.

[26] P.J. Morrisson, E. Tassi, C. Tronci, *Energy stability analysis for a hybrid fluid-kinetic plasma model*, In Nonlinear Physical Systems-Spectral Analysis, Stability and Bifurcations (ed. O. Kirillov & D. Pelinovsky), pp. 311–329. Wiley (2014).

[27] Ph. Nicolai, J.-L. Feugeas, and G. Schurtz, *A practical nonlocal model for heat transport in magnetized laser plasmas*, Phys. Plasmas 13(3) (2006).

[28] B. Peigney, O. Larroche, V. Tikhonchuk, *Fokker-Planck kinetic modeling of suprathermal α-particles in a fusion plasma*, J. Comput. Phys. 278, pp. 416-444 (2014).

[29] O. Pezzi, F. Valentini, P. Veltri, *Nonlinear regime of electrostatic waves propagation in presence of electron-electron collisions*, Phys. Plasmas, **22**, 042112 (2015).

[30] M.N. Rosenbluth, and W.M. MacDonald, and D.L. Judd, *Fokker-Planck Equation for an Inverse-Square Force*, Phys. Rev. 107(1), pp. 1-6, 1957.

[31] E. Sonnendrücker, *Numerical methods for Vlasov equations*, Lecture notes.

[32] V. Tikhonchuk, *Progress and opportunities for inertial fusion energy in Europe*, Phil. Trans. R. Soc. A 378.2184 (2020).

[33] C. Tronci, E. Tassi, E. Camporeale, P.J. MORRISON, *Hybrid vlasov-MHD models: Hamiltonian vs. non-Hamiltonian*, Plasma Phys. Control. Fusion 56 (9), 095008 (2014).
[34] B. A. Trubnikov, *Particle Interactions in a Fully Ionized Plasma*, Reviews of Plasma Physics **1** (1965), 105-140.

[35] J. Zhu, Z.W. Ma, S. Wang, *Hybrid simulations of Alfvén modes driven by energetic particles*, Phys. Plasmas **23** (12), 122506 (2016).

[36] E. Zeidler, *Nonlinear Functional Analysis and Its Applications, II/A: Linear Monotone Operators*, Springer Verlag New York Inc., 1990.

[37] D. Zwillinger, *Table of Integrals, Series, and Products*, Academic Press Inc. 8th Edition, 2014.

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

Email address: nicolas.crouseilles@inria.fr, claudia.negulescu@math.univ-toulouse.fr