# Hyperbolic reduced model for Vlasov-Poisson equation with Fokker-Planck collision

Emmanuel Franck, Ibtissem Lannabi, Youssouf Nasseri, Laurent Navoret, Giuseppe Parasiliti Rantone, Guillaume Steimer

► To cite this version:

Emmanuel Franck, Ibtissem Lannabi, Youssouf Nasseri, Laurent Navoret, Giuseppe Parasiliti Rantone, et al.. Hyperbolic reduced model for Vlasov-Poisson equation with Fokker-Planck collision. 2024. hal-04099697v3

HAL Id: hal-04099697

https://hal.science/hal-04099697v3

Preprint submitted on 15 Apr 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Hyperbolic reduced model for Vlasov-Poisson equation with Fokker-Planck collision

Emmanuel Franck, Ibtissem Lannabi, Youssouf Nasseri, Laurent Navoret, Giuseppe Parasiliti Rantone, Guillaume Steimer

► To cite this version:

Emmanuel Franck, Ibtissem Lannabi, Youssouf Nasseri, Laurent Navoret, Giuseppe Parasiliti Rantone, et al.. Hyperbolic reduced model for Vlasov-Poisson equation with Fokker-Planck collision. 2024. hal-04099697v3

HAL Id: hal-04099697

https://hal.science/hal-04099697v3

Preprint submitted on 15 Apr 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Hyperbolic reduced model for Vlasov-Poisson equation
with Fokker-Planck collision

Emmanuel Franck*,¹, Ibtissem Lannabi†,², Youssouf Nasseri‡,¹, Laurent Navoret§,¹,
Giuseppe Parasiliti Rantone¶,³, and Guillaume Steimer||,¹

¹Institut de Recherche Mathématique Avancée, UMR 7501, Université de Strasbourg et CNRS, 7 rue René Descartes, 67000 Strasbourg, France & INRIA Nancy-Grand Est, TONUS Project, Strasbourg, France

²CNRS/Univ Pau & Pays Adour/E2S UPPA, Laboratoire de Mathématiques et de leurs Applications de Pau - Fédération IPRA, UMR5142 64000, Pau, France & Cagire team, Inria Bordeaux Sud-Ouest, France

³ Sorbonne Université, CNRS UMR 7190, Institut Jean Le Rond d'Alembert, 75005 Paris, France

April 15, 2024

Abstract

This paper proposes a reduced model to simulate the one-dimensional Vlasov-Poisson equation with the non-linear Fokker-Planck operator. The model provides the space-time dynamics of a few macroscopic quantities constructed following the Reduced Order Method (ROM) in the velocity variable: the compression is thus applied to the semi-discretization of the Vlasov equation. To gain efficiency, a Discrete Empirical Interpolation Method (DEIM) is applied to the compressed non-linear Fokker-Planck operator. The size of the resulting reduced model is chosen empirically according to the Knudsen number. Furthermore, we propose a correction to the reduced collision operator that ensures the reduced moments to satisfy an Euler-type system. Numerical simulations of the reduced model show that the model can capture the plasma dynamics in different collisional regimes and initial conditions at a low cost.

# 1 Introduction

The Vlasov-Poisson-Fokker-Planck equation is a model for the transport of the distribution function of charged particles in the six-dimensional position-velocity phase space. The non-linear Fokker-Planck operator describes the short-range binary interactions between charged particles, called collisions. The weight of collisions in the dynamics is measured by the dimensionless Knudsen number $\epsilon$, the scaled mean-free path between collisions: a small Knudsen number corresponds to a collisional regime.

Simulations of such dynamics are very computationally demanding. Indeed, since phase space is of dimension six, simulations require a lot of memory and CPU resources. In addition, capturing collisional dynamics leads to numerical constraints on the time step and/or phase-space discretization. Indeed, the Fokker-Planck operator is a diffusive operator in the velocity variable: stability conditions for explicit numerical schemes have very stringent stability conditions when the Knudsen number $\epsilon$ is small. We refer to [BDS01, CA22] and references therein.

<footnote>*emmanuel.franck@inria.fr</footnote>

<footnote>†ibtissem.lannabi@univ-pau.fr</footnote>

<footnote>‡youssouf.nasseri@inria.fr</footnote>

<footnote>§laurent.navoret@math.unistra.fr</footnote>

<footnote>¶giuseppe.parasiliti_rantone@sorbonne-universite.fr</footnote>

<footnote>lguillaume.steimer@inria.fr</footnote>
Hyperbolic reduced model for Vlasov-Poisson equation
with Fokker-Planck collision

Emmanuel Franck*,¹, Ibtissem Lannabi†,², Youssouf Nasseri‡,¹, Laurent Navoret§,¹,
Giuseppe Parasiliti Rantone¶,³, and Guillaume Steimer||,¹

¹Institut de Recherche Mathématique Avancée, UMR 7501, Université de Strasbourg et CNRS, 7 rue René Descartes, 67000 Strasbourg, France & INRIA Nancy-Grand Est, TONUS Project, Strasbourg, France

²CNRS/Univ Pau & Pays Adour/E2S UPPA, Laboratoire de Mathématiques et de leurs Applications de Pau - Fédération IPRA, UMR5142 64000, Pau, France & Cagire team, Inria Bordeaux Sud-Ouest, France

³ Sorbonne Université, CNRS UMR 7190, Institut Jean Le Rond d'Alembert, 75005 Paris, France

April 15, 2024

Abstract

This paper proposes a reduced model to simulate the one-dimensional Vlasov-Poisson equation with the non-linear Fokker-Planck operator. The model provides the space-time dynamics of a few macroscopic quantities constructed following the Reduced Order Method (ROM) in the velocity variable: the compression is thus applied to the semi-discretization of the Vlasov equation. To gain efficiency, a Discrete Empirical Interpolation Method (DEIM) is applied to the compressed non-linear Fokker-Planck operator. The size of the resulting reduced model is chosen empirically according to the Knudsen number. Furthermore, we propose a correction to the reduced collision operator that ensures the reduced moments to satisfy an Euler-type system. Numerical simulations of the reduced model show that the model can capture the plasma dynamics in different collisional regimes and initial conditions at a low cost.

# 1 Introduction

The Vlasov-Poisson-Fokker-Planck equation is a model for the transport of the distribution function of charged particles in the six-dimensional position-velocity phase space. The non-linear Fokker-Planck operator describes the short-range binary interactions between charged particles, called collisions. The weight of collisions in the dynamics is measured by the dimensionless Knudsen number $\epsilon$, the scaled mean-free path between collisions: a small Knudsen number corresponds to a collisional regime.

Simulations of such dynamics are very computationally demanding. Indeed, since phase space is of dimension six, simulations require a lot of memory and CPU resources. In addition, capturing collisional dynamics leads to numerical constraints on the time step and/or phase-space discretization. Indeed, the Fokker-Planck operator is a diffusive operator in the velocity variable: stability conditions for explicit numerical schemes have very stringent stability conditions when the Knudsen number $\epsilon$ is small. We refer to [BDS01, CA22] and references therein.

<footnote>*emmanuel.franck@inria.fr</footnote>

<footnote>†ibtissem.lannabi@univ-pau.fr</footnote>

<footnote>‡youssouf.nasseri@inria.fr</footnote>

<footnote>§laurent.navoret@math.unistra.fr</footnote>

<footnote>¶giuseppe.parasiliti_rantone@sorbonne-universite.fr</footnote>

<footnote>lguillaume.steimer@inria.fr</footnote>
In order to avoid full model simulations, a classical strategy is to design reduced models that are valid in some parameter regimes. Typically, these reduced models provide the space-time dynamics of macroscopic quantities, which are integrated quantities over the velocity variables. For instance, the Euler system, satisfied by the density, the momentum, and the energy, can be obtained with the moment method and is valid in the strong collisional regime ($\epsilon \ll 1$). The Navier-Stokes equations can be considered for lower Knudsen numbers with $\epsilon < 10^{-2}$. Models with more moments have been designed in order to be valid at higher Knudsen number regimes.

Another possible approach is to consider reduced models based on data as proposed by the Reduced Order Modeling (ROM). The method consists in constructing an adapted basis in which the solution can be well approximated with few components. The basis is computed through a Proper Orthogonal decomposition applied to samples of the unknown in the considered physical regime. The reduced model is then obtained with a projection Galerkin method. We refer to [HPR22b] for more references. This method has been applied for particle discretization of the Vlasov-Poisson system [HPR22a] with a reduction in both space and velocity. In order to efficiently tackle Eulerian discretization in the six-dimensional phase space, low-rank tensor bases have been proposed [FL22, VD17], as well as their dynamical version [GEQ23].

Here we propose a mixed method where the ROM approach only compresses the dynamics in the velocity variable. Like moment models, the resulting reduced model provides the spatial dynamics of the reduced quantities. Therefore, the method provides a model independent of the spatial discretization of the computational domain. In order to provide a first assessment of the method, this study focuses on the one-dimensional dynamics.

The non-linear collision operator leads to a reduced operator that would require to alternate between reduced and full dimensional data. To avoid these costly computations, we propose to use the Discrete Empirical Interpolation Method (DEIM) [CS10]. We also propose a correction of the reduced collision operator to ensure that it preserves reduced moments (like mass, momentum, and energy). Therefore, the reduced moments associated with the reduced data will satisfy an Euler-type equation.

This article is structured as follows. Section 2 introduces the Vlasov-Poisson-Fokker-Planck model, its semi-discretization in velocity, and the reduced model obtained after using the ROM approach in velocity. In Section 3, we recall the DEIM strategy to reduce the computational complexity for the non-linear reduced collision operator. Then we present a correction to it such that it preserves the moments. Finally, in Section 4, we assess the capability of the reduced model in capturing non-linear and linear dynamics. To this aim, we perform a Landau damping test case with different Knudsen numbers and perturbation amplitudes.

# 2 Reduced Order Modeling in velocity of the one-dimensional Vlasov-Poisson Fokker-Planck model

## 2.1 Vlasov-Poisson Fokker-Planck model

We consider the one-dimensional Vlasov-Poisson-Fokker-Planck dynamics that describes the evolution of a statistical distribution of charged particles in position-velocity phase-space, with long-range interactions through the self-induced electric field and short-range interactions. We define $f(t, x, v)$ as the distribution of particles at time $t \ge 0$ in the phase space $(x, v) \in [0, L] \times \mathbb{R}$ and the system reads:

$$ \partial_t f(t, x, v) + v \partial_x f(t, x, v) + \partial_x \phi(t, x) \partial_v f(t, x, v) = \frac{1}{\epsilon} Q(f)(t, x, v), \quad (1) $$

$$ -\partial_x^2 \phi(t, x) = \frac{1}{L} \int_{[0,L] \times \mathbb{R}} f(t, y, v) dy dv - \int_{\mathbb{R}} f(t, x, v) dv, \quad (2) $$

where $\phi(t, x)$ denotes the electric potential. The electric field is defined by the relation: $E(t, x) = -\partial_x \phi(t, x)$.

The right-hand side of the equation models collision interactions: $Q(f)$ models the short-range interactions, and the positive parameter $\epsilon$ represents the collision timescale. This collision operator is local in space and acts only on the velocity variable. Here, we consider a non-linear Fokker-Planck collision operator given by:

$$ Q(f) = \partial_v \left( (v - u_f) f + T_f \partial_v f \right), $$
In order to avoid full model simulations, a classical strategy is to design reduced models that are valid in some parameter regimes. Typically, these reduced models provide the space-time dynamics of macroscopic quantities, which are integrated quantities over the velocity variables. For instance, the Euler system, satisfied by the density, the momentum, and the energy, can be obtained with the moment method and is valid in the strong collisional regime ($\epsilon \ll 1$). The Navier-Stokes equations can be considered for lower Knudsen numbers with $\epsilon < 10^{-2}$. Models with more moments have been designed in order to be valid at higher Knudsen number regimes.

Another possible approach is to consider reduced models based on data as proposed by the Reduced Order Modeling (ROM). The method consists in constructing an adapted basis in which the solution can be well approximated with few components. The basis is computed through a Proper Orthogonal decomposition applied to samples of the unknown in the considered physical regime. The reduced model is then obtained with a projection Galerkin method. We refer to [HPR22b] for more references. This method has been applied for particle discretization of the Vlasov-Poisson system [HPR22a] with a reduction in both space and velocity. In order to efficiently tackle Eulerian discretization in the six-dimensional phase space, low-rank tensor bases have been proposed [FL22, VD17], as well as their dynamical version [GEQ23].

Here we propose a mixed method where the ROM approach only compresses the dynamics in the velocity variable. Like moment models, the resulting reduced model provides the spatial dynamics of the reduced quantities. Therefore, the method provides a model independent of the spatial discretization of the computational domain. In order to provide a first assessment of the method, this study focuses on the one-dimensional dynamics.

The non-linear collision operator leads to a reduced operator that would require to alternate between reduced and full dimensional data. To avoid these costly computations, we propose to use the Discrete Empirical Interpolation Method (DEIM) [CS10]. We also propose a correction of the reduced collision operator to ensure that it preserves reduced moments (like mass, momentum, and energy). Therefore, the reduced moments associated with the reduced data will satisfy an Euler-type equation.

This article is structured as follows. Section 2 introduces the Vlasov-Poisson-Fokker-Planck model, its semi-discretization in velocity, and the reduced model obtained after using the ROM approach in velocity. In Section 3, we recall the DEIM strategy to reduce the computational complexity for the non-linear reduced collision operator. Then we present a correction to it such that it preserves the moments. Finally, in Section 4, we assess the capability of the reduced model in capturing non-linear and linear dynamics. To this aim, we perform a Landau damping test case with different Knudsen numbers and perturbation amplitudes.

# 2 Reduced Order Modeling in velocity of the one-dimensional Vlasov-Poisson Fokker-Planck model

## 2.1 Vlasov-Poisson Fokker-Planck model

We consider the one-dimensional Vlasov-Poisson-Fokker-Planck dynamics that describes the evolution of a statistical distribution of charged particles in position-velocity phase-space, with long-range interactions through the self-induced electric field and short-range interactions. We define $f(t, x, v)$ as the distribution of particles at time $t \ge 0$ in the phase space $(x, v) \in [0, L] \times \mathbb{R}$ and the system reads:

$$ \partial_t f(t, x, v) + v \partial_x f(t, x, v) + \partial_x \phi(t, x) \partial_v f(t, x, v) = \frac{1}{\epsilon} Q(f)(t, x, v), \quad (1) $$

$$ -\partial_x^2 \phi(t, x) = \frac{1}{L} \int_{[0,L] \times \mathbb{R}} f(t, y, v) dy dv - \int_{\mathbb{R}} f(t, x, v) dv, \quad (2) $$

where $\phi(t, x)$ denotes the electric potential. The electric field is defined by the relation: $E(t, x) = -\partial_x \phi(t, x)$.

The right-hand side of the equation models collision interactions: $Q(f)$ models the short-range interactions, and the positive parameter $\epsilon$ represents the collision timescale. This collision operator is local in space and acts only on the velocity variable. Here, we consider a non-linear Fokker-Planck collision operator given by:

$$ Q(f) = \partial_v \left( (v - u_f) f + T_f \partial_v f \right), $$
where $\rho_f(t, x)$, $u_f(t, x)$ and $T_f(t, x)$ denote the particle density (in space), the velocity, and the temperature, respectively, and are defined by:

$$ \rho_f(t, x) = \int_{\mathbb{R}} f(t, x, v) dv, \quad \rho_f(t, x) u_f(t, x) = \int_{\mathbb{R}} v f(t, x, v) dv, \qquad (3) $$

and

$$ p_f(t, x) = \rho_f(t, x) T_f(t, x) = \int_{\mathbb{R}} (v - u_f(t, x))^2 f(t, x, v) dv. \qquad (4) $$

The collision operator can be rewritten as:

$$ Q(f) = T_f \partial_v \left( M_f \partial_v \left( \frac{f}{M_f} \right) \right), \qquad (5) $$

where $M_f(v) = M_{\rho_f, u_f, T_f}(v)$ is the so-called Maxwellian distribution defined by:

$$ M_{\rho,u,T}(v) = \frac{\rho}{\sqrt{2\pi T}} e^{-\frac{(v-u)^2}{2T}}. \qquad (6) $$

This expression of $Q$ shows that the equilibrium of the collision operator (satisfying $Q(f) = 0$) is precisely the Maxwellian distribution.

To reduce the complexity of model (1)-(2), we want to project the velocity distribution function on a basis adapted to the desired physical parameters. For this purpose, we start by considering a finite-dimensional approximation of the dynamics.

## 2.2 Semi-discretized model in velocity

We perform a semi-discretization in velocity of the Vlasov-Poisson equations (1)-(2) using a finite difference method. We take a velocity interval $[-v_{\max}, v_{\max}]$ and an odd number of nodes $N_v$: $v_i = (i - (N_v - 1)/2)\Delta v$ with $\Delta v = 2v_{\max}/(N_v - 1)$ and $i \in \{0, N_v - 1\}$. The unknown vector $f(t, x) = (f(t, x, v_0), \dots, f(t, x, v_{N_v-1}))$ satisfies the following semi-discretized Vlasov equation:

$$ \partial_t f(t, x) + A \partial_x f(t, x) + (\partial_x \phi(t, x)) D f(t, x) = \frac{1}{\varepsilon} Q_{\Delta v}(f(t, x)), \qquad (7) $$

$$ - \Delta \phi(t, x) = \frac{1}{L} \int_{[0,L]} \Delta v \mathbf{1}^T \mathbf{f}(t, x) \, dx - \Delta v \mathbf{1}^T \mathbf{f}(t, x), \qquad (8) $$

where $A = \text{Diag}(\mathbf{v})$ denotes the diagonal matrix of size $N_v$ with the vector $\mathbf{v} = (v_0, \dots, v_{N_v-1})$ on the diagonal. Matrix $D$ refers to the matrix associated with the centered finite difference of the velocity derivative with Dirichlet boundary conditions:

$$ D = \frac{1}{2\Delta v} \begin{pmatrix} 0 & 1 & 0 & \cdots & 0 \\ -1 & \ddots & \ddots & \ddots & \vdots \\ 0 & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & 1 \\ 0 & \cdots & 0 & -1 & 0 \end{pmatrix}. $$

Symbol **1** denotes the vector $(1, \dots, 1)^T \in \mathbb{R}^{N_v}$. Equation (7) is a hyperbolic system with a source term.

For the discretization of the collision operator, we consider the scheme proposed in [BDS01] that is explicit, preserves the density, momentum, and energy, and dissipates entropy. Based on the following expression of the non-linear Fokker-Planck collision operator,

$$ Q(f) = T \partial_v \left( f \partial_v \log \left( \frac{f}{M_f} \right) \right), $$

the discretization can be written as:

$$ Q_{\Delta v}(\mathbf{f})_j = \frac{\mathbf{f}_{j+\frac{1}{2}} \left( \log \left( \frac{\mathbf{f}_{j+1}}{(\mathbf{M}_{\mathbf{f}})_{j+1}} \right) - \log \left( \frac{\mathbf{f}_i}{(\mathbf{M}_{\mathbf{f}})_j} \right) \right) - \mathbf{f}_{j-\frac{1}{2}} \left( \log \left( \frac{\mathbf{f}_j}{(\mathbf{M}_{\mathbf{f}})_j} \right) - \log \left( \frac{\mathbf{f}_{j-1}}{(\mathbf{M}_{\mathbf{f}})_{j-1}} \right) \right)}{\Delta v^2}, $$
where $\rho_f(t, x)$, $u_f(t, x)$ and $T_f(t, x)$ denote the particle density (in space), the velocity, and the temperature, respectively, and are defined by:

$$ \rho_f(t, x) = \int_{\mathbb{R}} f(t, x, v) dv, \quad \rho_f(t, x) u_f(t, x) = \int_{\mathbb{R}} v f(t, x, v) dv, \qquad (3) $$

and

$$ p_f(t, x) = \rho_f(t, x) T_f(t, x) = \int_{\mathbb{R}} (v - u_f(t, x))^2 f(t, x, v) dv. \qquad (4) $$

The collision operator can be rewritten as:

$$ Q(f) = T_f \partial_v \left( M_f \partial_v \left( \frac{f}{M_f} \right) \right), \qquad (5) $$

where $M_f(v) = M_{\rho_f, u_f, T_f}(v)$ is the so-called Maxwellian distribution defined by:

$$ M_{\rho,u,T}(v) = \frac{\rho}{\sqrt{2\pi T}} e^{-\frac{(v-u)^2}{2T}}. \qquad (6) $$

This expression of $Q$ shows that the equilibrium of the collision operator (satisfying $Q(f) = 0$) is precisely the Maxwellian distribution.

To reduce the complexity of model (1)-(2), we want to project the velocity distribution function on a basis adapted to the desired physical parameters. For this purpose, we start by considering a finite-dimensional approximation of the dynamics.

## 2.2 Semi-discretized model in velocity

We perform a semi-discretization in velocity of the Vlasov-Poisson equations (1)-(2) using a finite difference method. We take a velocity interval $[-v_{\max}, v_{\max}]$ and an odd number of nodes $N_v$: $v_i = (i - (N_v - 1)/2)\Delta v$ with $\Delta v = 2v_{\max}/(N_v - 1)$ and $i \in \{0, N_v - 1\}$. The unknown vector $\mathbf{f}(t, x) = (f(t, x, v_0), \dots, f(t, x, v_{N_v-1}))$ satisfies the following semi-discretized Vlasov equation:

$$ \partial_t \mathbf{f}(t, x) + A \partial_x \mathbf{f}(t, x) + (\partial_x \phi(t, x)) D \mathbf{f}(t, x) = \frac{1}{\varepsilon} Q_{\Delta v}(\mathbf{f}(t, x)), \qquad (7) $$

$$ - \Delta \phi(t, x) = \frac{1}{L} \int_{[0,L]} \Delta v \mathbf{1}^T \mathbf{f}(t, x) dx - \Delta v \mathbf{1}^T \mathbf{f}(t, x), \qquad (8) $$

where $A = \text{Diag}(\mathbf{v})$ denotes the diagonal matrix of size $N_v$ with the vector $\mathbf{v} = (v_0, \dots, v_{N_v-1})$ on the diagonal. Matrix $D$ refers to the matrix associated with the centered finite difference of the velocity derivative with Dirichlet boundary conditions:

$$ D = \frac{1}{2\Delta v} \begin{pmatrix} 0 & 1 & 0 & \cdots & 0 \\ -1 & \ddots & \ddots & \ddots & \vdots \\ 0 & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & 1 \\ 0 & \cdots & 0 & -1 & 0 \end{pmatrix}. $$

Symbol $\mathbf{1}$ denotes the vector $(1, \dots, 1)^T \in \mathbb{R}^{N_v}$. Equation (7) is a hyperbolic system with a source term.

For the discretization of the collision operator, we consider the scheme proposed in [BDS01] that is explicit, preserves the density, momentum, and energy, and dissipates entropy. Based on the following expression of the non-linear Fokker-Planck collision operator,

$$ Q(f) = T \partial_v \left( f \partial_v \log \left( \frac{f}{M_f} \right) \right), $$

the discretization can be written as:

$$ Q_{\Delta v}(\mathbf{f})_j = \frac{\mathbf{f}_{j+\frac{1}{2}} \left( \log \left( \frac{\mathbf{f}_{j+1}}{(\mathbf{M}_{\mathbf{f}})_{j+1}} \right) - \log \left( \frac{\mathbf{f}_i}{(\mathbf{M}_{\mathbf{f}})_j} \right) \right) - \mathbf{f}_{j-\frac{1}{2}} \left( \log \left( \frac{\mathbf{f}_j}{(\mathbf{M}_{\mathbf{f}})_j} \right) - \log \left( \frac{\mathbf{f}_{j-1}}{(\mathbf{M}_{\mathbf{f}})_{j-1}} \right) \right)}{\Delta v^2}, $$
where $\mathbf{f}_{\frac{1}{2}} = \mathbf{f}_{N_v - \frac{1}{2}} = 0$ and $\mathbf{f}_{j + \frac{1}{2}}$ is the entropic average:

$$ \text{for } j = 1, \dots, N_v - 2, \quad \mathbf{f}_{j + \frac{1}{2}} = \begin{cases} \frac{\mathbf{f}_{j+1} - \mathbf{f}_j}{\log(\mathbf{f}_{j+1}) - \log(\mathbf{f}_j)}, & \text{if } \mathbf{f}_{j+1} \neq \mathbf{f}_j, \\ \mathbf{f}_j, & \text{otherwise.} \end{cases} \qquad (9) $$

The discretization involves the discrete Maxwellian, which equals the evaluation of the Maxwellian at the velocity grid point, $\mathbf{M}_{\mathbf{f}} = (M_{\rho_{\mathbf{f}}, \tilde{u}_{\mathbf{f}}, (\tilde{\rho}_{\mathbf{f}}, \tilde{u}_{\mathbf{f}}, \tilde{T}_{\mathbf{f}})/\rho_{\mathbf{f}}}(v_j))_j$, where $\tilde{\rho}_{\mathbf{f}}, \tilde{u}_{\mathbf{f}}, \tilde{T}_{\mathbf{f}}$ refer to the modified discrete moments:

$$ \tilde{\rho}_{\mathbf{f}} = \sum_{j=1}^{N_v-2} \Delta v \mathbf{f}_{j+\frac{1}{2}}, $$

$$ \tilde{\rho}_{\mathbf{f}} \tilde{u}_{\mathbf{f}} = \sum_{j=1}^{N_v-2} \Delta v v_{j+\frac{1}{2}} \mathbf{f}_{j+\frac{1}{2}} + T_{\mathbf{f}}(\mathbf{f}_{N_v} - \mathbf{f}_1), $$

$$ \tilde{\rho}_{\mathbf{f}} \tilde{T}_{\mathbf{f}} = \sum_{j=1}^{N_v-2} \Delta v (v_{j+\frac{1}{2}} - \tilde{u}_{\mathbf{f}})^2 \mathbf{f}_{j+\frac{1}{2}} + \tilde{\rho}_{\mathbf{f}} T_{\mathbf{f}} (\mathbf{f}_{N_v} (v_{N_v+\frac{1}{2}} - \tilde{u}_{\mathbf{f}}) + \mathbf{f}_1 (v_{\frac{1}{2}} - \tilde{u}_{\mathbf{f}})), $$

and where $\rho_{\mathbf{f}}, u_{\mathbf{f}}, T_{\mathbf{f}}$ refer to the classical discrete moments associated with the discrete distribution $\mathbf{f}$ defined by:

$$ \rho_{\mathbf{f}} = \Delta v \mathbf{1}^T \mathbf{f}, $$

$$ \rho_{\mathbf{f}} u_{\mathbf{f}} = \Delta v \mathbf{1}^T \mathrm{Diag}(\mathbf{v}) \mathbf{f}, $$

$$ \rho_{\mathbf{f}} u_{\mathbf{f}}^2 + \rho_{\mathbf{f}} T_{\mathbf{f}} = \Delta v \mathbf{1}^T \mathrm{Diag}(\mathbf{v})^2 \mathbf{f}. $$

We also define the moment matrix $m$:

$$ m = \begin{bmatrix} m_{\rho} \\ m_{\rho u} \\ m_{w} \end{bmatrix} = \begin{bmatrix} \Delta v \mathbf{1}^T \\ \Delta v \mathbf{1}^T \mathrm{Diag}(\mathbf{v}) \\ \Delta v \mathbf{1}^T \mathrm{Diag}(\mathbf{v})^2 \end{bmatrix} \in M_{3,N_v}(\mathbb{R}). $$

As already said, the discrete moments $\rho_{\mathbf{f}}$ and $\rho_{\mathbf{f}} u_{\mathbf{f}}$ are conserved by the scheme. We also note that this discrete collision operator vanishes on the discrete Maxwellian with the modified mean velocity and temperature.

## 2.3 Reduced model

We now apply the ROM methodology. Using a Proper Orthogonal Decomposition (POD), we first generate a reduced basis of size $K \ll N_v$ from samples of the distribution function and define the linear decompression operator $\Phi$, that we apply to a reduced vector $\hat{f}$ to get the full vector $\mathbf{f}$, and the compression operator $\Phi^T$, that we apply to the full vector $\mathbf{f}$ to get a reduced vector $\hat{f}$:

$$ \mathbf{f} \in \mathbb{R}^{N_v} \xrightarrow{\Phi^T} \hat{\mathbf{f}} \in \mathbb{R}^K \xrightarrow{\Phi} \tilde{\mathbf{f}} \in \mathbb{R}^{N_v}. $$

Then we perform a Galerkin projection to obtain the reduced model on the reduced quantity $\hat{\mathbf{f}}$.

**Reduced basis and compression/decompression operators.** We first consider $N_s$ samples of the distribution in time and space obtained from the numerical simulations of the model of the previous section. Spatial discretization is done using finite volume schemes, and the time discretization is done using an explicit first-order scheme. We refer to Appendix A for more details. The sample matrix is then given by:

$$ X = [\mathbf{f}_1(t_1, x_1), \dots, \mathbf{f}_n(t_n, x_n)] \in M_{N_v, N_s}(\mathbb{R}), $$

where $(\mathbf{f}_i)_i$ are obtained from $N_s$ possible different initial data and different Knudsen numbers and $(t_i, x_i)$ are chosen randomly. The linear decompression operator $\Phi \in M_{N_v, k}(\mathbb{R})$ is defined such that the compression-decompression error on the samples is minimized:

$$ \min_{\substack{\Phi \in M_{N_v, k}(\mathbb{R}) \\ \Phi^T \Phi = \mathrm{Id}}} \|X - \Phi \Phi^T X\|_F^2, \qquad (10) $$
where $\mathbf{f}_{\frac{1}{2}} = \mathbf{f}_{N_v - \frac{1}{2}} = 0$ and $\mathbf{f}_{j + \frac{1}{2}}$ is the entropic average:

$$ \text{for } j = 1, \dots, N_v - 2, \quad \mathbf{f}_{j + \frac{1}{2}} = \begin{cases} \frac{\mathbf{f}_{j+1} - \mathbf{f}_j}{\log(\mathbf{f}_{j+1}) - \log(\mathbf{f}_j)}, & \text{if } \mathbf{f}_{j+1} \neq \mathbf{f}_j, \\ \mathbf{f}_j, & \text{otherwise.} \end{cases} \qquad (9) $$

The discretization involves the discrete Maxwellian, which equals the evaluation of the Maxwellian at the velocity grid point, $\mathbf{M}_f = (M_{\rho_f, \tilde{u}_f, (\tilde{\rho}_f \tilde{T}_f) / \rho_f}(v_j))_j$, where $\tilde{\rho}_f, \tilde{u}_f, \tilde{T}_f$ refer to the modified discrete moments:

$$ \tilde{\rho}_f = \sum_{j=1}^{N_v-2} \Delta v \mathbf{f}_{j+\frac{1}{2}}, $$

$$ \tilde{\rho}_f \tilde{u}_f = \sum_{j=1}^{N_v-2} \Delta v v_{j+\frac{1}{2}} \mathbf{f}_{j+\frac{1}{2}} + T_f(\mathbf{f}_{N_v} - \mathbf{f}_1), $$

$$ \tilde{\rho}_f \tilde{T}_f = \sum_{j=1}^{N_v-2} \Delta v (v_{j+\frac{1}{2}} - \tilde{u}_f)^2 \mathbf{f}_{j+\frac{1}{2}} + \tilde{\rho}_f T_f (\mathbf{f}_{N_v} (v_{N_v+\frac{1}{2}} - \tilde{u}_f) + \mathbf{f}_1 (v_{\frac{1}{2}} - \tilde{u}_f)), $$

and where $\rho_f, u_f, T_f$ refer to the classical discrete moments associated with the discrete distribution $f$ defined by:

$$ \rho_f = \Delta v \mathbf{1}^T \mathbf{f}, $$

$$ \rho_f u_f = \Delta v \mathbf{1}^T \mathrm{Diag}(\mathbf{v}) \mathbf{f}, $$

$$ \rho_f u_f^2 + \rho_f T_f = \Delta v \mathbf{1}^T \mathrm{Diag}(\mathbf{v})^2 \mathbf{f}. $$

We also define the moment matrix $m$:

$$ m = \begin{bmatrix} m_{\rho} \\ m_{\rho u} \\ m_{w} \end{bmatrix} = \begin{bmatrix} \Delta v \mathbf{1}^T \\ \Delta v \mathbf{1}^T \mathrm{Diag}(\mathbf{v}) \\ \Delta v \mathbf{1}^T \mathrm{Diag}(\mathbf{v})^2 \end{bmatrix} \in M_{3,N_v}(\mathbb{R}). $$

As already said, the discrete moments $\rho_f$ and $\rho_f u_f$ are conserved by the scheme. We also note that this discrete collision operator vanishes on the discrete Maxwellian with the modified mean velocity and temperature.

## 2.3 Reduced model

We now apply the ROM methodology. Using a Proper Orthogonal Decomposition (POD), we first generate a reduced basis of size $K \ll N_v$ from samples of the distribution function and define the linear decompression operator $\Phi$, that we apply to a reduced vector $\hat{f}$ to get the full vector $f$, and the compression operator $\Phi^T$, that we apply to the full vector $f$ to get a reduced vector $\hat{f}$:

$$ \mathbf{f} \in \mathbb{R}^{N_v} \xrightarrow{\Phi^T} \hat{\mathbf{f}} \in \mathbb{R}^K \xrightarrow{\Phi} \tilde{\mathbf{f}} \in \mathbb{R}^{N_v}. $$

Then we perform a Galerkin projection to obtain the reduced model on the reduced quantity $\hat{f}$.

**Reduced basis and compression/decompression operators.** We first consider $N_s$ samples of the distribution in time and space obtained from the numerical simulations of the model of the previous section. Spatial discretization is done using finite volume schemes, and the time discretization is done using an explicit first-order scheme. We refer to Appendix A for more details. The sample matrix is then given by:

$$ X = [\mathbf{f}_1(t_1, x_1), \dots, \mathbf{f}_n(t_n, x_n)] \in M_{N_v, N_s}(\mathbb{R}), $$

where $(\mathbf{f}_i)_i$ are obtained from $N_s$ possible different initial data and different Knudsen numbers and $(t_i, x_i)$ are chosen randomly. The linear decompression operator $\Phi \in M_{N_v, k}(\mathbb{R})$ is defined such that the compression-decompression error on the samples is minimized:

$$ \min_{\substack{\Phi \in M_{N_v, k}(\mathbb{R}) \\ \Phi^T \Phi = \mathrm{Id}}} \|X - \Phi \Phi^T X\|_F^2, \qquad (10) $$
where $\|A\|_F = \sqrt{\sum_{ij} A_{ij}^2}$ denotes the Frobenius norm for matrices. In practice, the columns of $\Phi$ are given by the first $K$ eigenvectors of the matrix $X^T X \in M_{N_v}(\mathbb{R})$. These first $K$ eigenvectors define the reduced basis.

**Reduced model by Galerkin projection.** Applying the Galerkin projection method to equations (7)-(8) consists in inserting the ansatz $\mathbf{f} = \Phi \hat{\mathbf{f}}$ into the equation, applying the projection $\Phi^T$ and then using the identity $\Phi^T \Phi = \text{Id}$. We get the following model on the reduced variable $\hat{\mathbf{f}}$:

$$ \partial_t \hat{\mathbf{f}}(t, x) + \hat{A} \partial_x \hat{\mathbf{f}}(t, x) + (\partial_x \phi(t, x)) \hat{D} \hat{\mathbf{f}}(t, x) = \frac{1}{\varepsilon} \hat{Q}_{\Delta v}(\hat{\mathbf{f}}), \quad (11) $$

$$ - \Delta \phi(t, x) = \hat{m}_{\rho} \hat{\mathbf{f}}(t, x), \quad (12) $$

where matrices $\hat{A}, \hat{D}, m_{\rho}$ are respectively of size $K \times K$, $K \times K$ and $1 \times K$ and given by:

$$ \hat{A} = \Phi^T \text{Diag}(\mathbf{v}) \Phi, \quad \hat{D} = \Phi^T D \Phi, \quad \hat{m}_{\rho} = \Delta v \mathbf{1}^T \Phi, $$

and the reduced collision operator is defined by:

$$ \hat{Q}_{\Delta v}(\hat{\mathbf{f}}) = \Phi^T Q_{\Delta v}(\Phi \hat{\mathbf{f}}(t, x)). \quad (13) $$

This system is a hyperbolic system of size $K$ with two source terms: the first term is linear and comes from the transport in the velocity variable, while the other one is non-linear and comes from the collision operator.

# 3 Hyper-reduction and corrections of the reduced collision operator

The non-linear collision operator (13) requires to evaluate the collision operator $Q_{\Delta v}$ on the full vector $\Phi \hat{\mathbf{f}}$ of size $N_v$. Since this could be quite numerically expensive, several techniques called hyper-reduction have been developed. In the present work, we propose to use the classical Discrete Empirical Interpolation Method [CS10].

The original discrete collision operator has been constructed to preserve mass, momentum, and energy. Such properties are crucial to capture the appropriate physical dynamics. We aim to ensure the same properties on the reduced collision operator. We thus define reduced moments and an associated corrected collision operator that preserves them.

## 3.1 DEIM hyper-reduction

Let us briefly describe the DEIM method. Using the Proper Orthogonal Decomposition on samples of the non-linear collision operator $Q_{\Delta v}(\mathbf{f})$ (without reduction), the method first generates a basis made of $K_Q$ vectors of size $N_v$: let $\Phi_Q$ the matrix of size $N_v \times K_Q$ that gathers these vectors. Then for any $\hat{\mathbf{f}}$, we are looking at coefficients $c(\hat{\mathbf{f}}) \in \mathbb{R}^{K_Q}$ such that:

$$ Q_{\Delta v}(\Phi \hat{\mathbf{f}}) \approx \Phi_Q c(\hat{\mathbf{f}}), \quad (14) $$

The coefficients are chosen such that $K_Q$ rows are indeed equalities (hence the name interpolation):

$$ P_Q^T Q_{\Delta v}(\Phi \hat{\mathbf{f}}) = P_Q^T \Phi_Q c(\hat{\mathbf{f}}), \quad (15) $$

where $P_Q$ is a selection matrix of size $N_v \times K_Q$, whose columns are made of canonical basis vectors. The selection of these $K_Q$ rows is made iteratively with a greedy-like algorithm. In particular, $P_Q$ is determined such that the matrix $P_Q^T \Phi_Q$ is invertible. Hence plugging back (15) into (14), we get the following approximation:

$$ Q_{\Delta v}(\Phi \hat{\mathbf{f}}) \approx \Phi_Q (P_Q^T \Phi_Q)^{-1} P_Q^T Q_{\Delta v}(\Phi \hat{\mathbf{f}}), $$

and finally the operator $\Phi^T Q_{\Delta v}(\Phi \hat{\mathbf{f}})$ in expression (13) is replaced by

$$ \hat{Q}_{\Delta v}^{\text{DEIM}}(\hat{\mathbf{f}}) = \Phi^T \Phi_Q (P_Q^T \Phi_Q)^{-1} P_Q^T Q_{\Delta v}(\Phi \hat{\mathbf{f}}), $$

We refer to [CS10] for more details. The computational gain of this method comes from the fact that only the $K_Q$ selected rows (among $N_v$) of $Q_{\Delta v}(\Phi \hat{\mathbf{f}})$ are actually computed.
where $\|A\|_F = \sqrt{\sum_{ij} A_{ij}^2}$ denotes the Frobenius norm for matrices. In practice, the columns of $\Phi$ are given by the first $K$ eigenvectors of the matrix $X^T X \in M_{N_v}(\mathbb{R})$. These first $K$ eigenvectors define the reduced basis.

**Reduced model by Galerkin projection.** Applying the Galerkin projection method to equations (7)-(8) consists in inserting the ansatz $\mathbf{f} = \Phi \hat{\mathbf{f}}$ into the equation, applying the projection $\Phi^T$ and then using the identity $\Phi^T \Phi = \text{Id}$. We get the following model on the reduced variable $\hat{\mathbf{f}}$:

$$ \partial_t \hat{\mathbf{f}}(t, x) + \hat{A} \partial_x \hat{\mathbf{f}}(t, x) + (\partial_x \phi(t, x)) \hat{D} \hat{\mathbf{f}}(t, x) = \frac{1}{\varepsilon} \hat{Q}_{\Delta v}(\hat{\mathbf{f}}), \quad (11) $$

$$ - \Delta \phi(t, x) = \hat{m}_{\rho} \hat{\mathbf{f}}(t, x), \quad (12) $$

where matrices $\hat{A}, \hat{D}, m_{\rho}$ are respectively of size $K \times K$, $K \times K$ and $1 \times K$ and given by:

$$ \hat{A} = \Phi^T \text{Diag}(\mathbf{v}) \Phi, \quad \hat{D} = \Phi^T D \Phi, \quad \hat{m}_{\rho} = \Delta v \mathbf{1}^T \Phi, $$

and the reduced collision operator is defined by:

$$ \hat{Q}_{\Delta v}(\hat{\mathbf{f}}) = \Phi^T Q_{\Delta v}(\Phi \hat{\mathbf{f}}(t, x)). \quad (13) $$

This system is a hyperbolic system of size $K$ with two source terms: the first term is linear and comes from the transport in the velocity variable, while the other one is non-linear and comes from the collision operator.

# 3 Hyper-reduction and corrections of the reduced collision operator

The non-linear collision operator (13) requires to evaluate the collision operator $Q_{\Delta v}$ on the full vector $\Phi \hat{\mathbf{f}}$ of size $N_v$. Since this could be quite numerically expensive, several techniques called hyper-reduction have been developed. In the present work, we propose to use the classical Discrete Empirical Interpolation Method [CS10].

The original discrete collision operator has been constructed to preserve mass, momentum, and energy. Such properties are crucial to capture the appropriate physical dynamics. We aim to ensure the same properties on the reduced collision operator. We thus define reduced moments and an associated corrected collision operator that preserves them.

## 3.1 DEIM hyper-reduction

Let us briefly describe the DEIM method. Using the Proper Orthogonal Decomposition on samples of the non-linear collision operator $Q_{\Delta v}(\mathbf{f})$ (without reduction), the method first generates a basis made of $K_Q$ vectors of size $N_v$: let $\Phi_Q$ the matrix of size $N_v \times K_Q$ that gathers these vectors. Then for any $\hat{\mathbf{f}}$, we are looking at coefficients $c(\hat{\mathbf{f}}) \in \mathbb{R}^{K_Q}$ such that:

$$ Q_{\Delta v}(\Phi \hat{\mathbf{f}}) \approx \Phi_Q c(\hat{\mathbf{f}}), \quad (14) $$

The coefficients are chosen such that $K_Q$ rows are indeed equalities (hence the name interpolation):

$$ P_Q^T Q_{\Delta v}(\Phi \hat{\mathbf{f}}) = P_Q^T \Phi_Q c(\hat{\mathbf{f}}), \quad (15) $$

where $P_Q$ is a selection matrix of size $N_v \times K_Q$, whose columns are made of canonical basis vectors. The selection of these $K_Q$ rows is made iteratively with a greedy-like algorithm. In particular, $P_Q$ is determined such that the matrix $P_Q^T \Phi_Q$ is invertible. Hence plugging back (15) into (14), we get the following approximation:

$$ Q_{\Delta v}(\Phi \hat{\mathbf{f}}) \approx \Phi_Q (P_Q^T \Phi_Q)^{-1} P_Q^T Q_{\Delta v}(\Phi \hat{\mathbf{f}}), $$

and finally the operator $\Phi^T Q_{\Delta v}(\Phi \hat{\mathbf{f}})$ in expression (13) is replaced by

$$ \hat{Q}_{\Delta v}^{\text{DEIM}}(\hat{\mathbf{f}}) = \Phi^T \Phi_Q (P_Q^T \Phi_Q)^{-1} P_Q^T Q_{\Delta v}(\Phi \hat{\mathbf{f}}), $$

We refer to [CS10] for more details. The computational gain of this method comes from the fact that only the $K_Q$ selected rows (among $N_v$) of $Q_{\Delta v}(\Phi \hat{\mathbf{f}})$ are actually computed.
## 3.2 Preservation of reduced moments

The reduced collision operator does not a priori satisfy the preservation of the reduced moments nor vanishes on discrete Maxwellians. We propose to slightly change the discrete collision operator to enforce these two properties, as already proposed in [Rif20]. Similar corrections have been proposed in the context of spectral methods [GT09], discontinuous Galerkin methods [ZG18] or discrete-velocity methods [DL13].

We first introduced the reduced moment operator as follows:

$$ \hat{m} = \begin{bmatrix} \hat{m}_{\rho} \\ \hat{m}_{\rho u} \\ \hat{m}_{w} \end{bmatrix} = \begin{bmatrix} \Delta v \mathbf{1}^T \\ \Delta v \mathbf{1}^T \Phi (\Phi^T \text{Diag}(\mathbf{v}) \Phi) \\ \frac{1}{2} \Delta v \mathbf{1}^T \Phi (\Phi^T \text{Diag}(\mathbf{v}) \Phi)^2 \end{bmatrix} \in M_{3,K}(\mathbb{R}). \quad (16) $$

Note that this definition differs from the moment operator $m\Phi$, which would be the natural definition since it first decompresses up to the full size and then applies the classical moment operator. Instead, the considered reduced moment operator first multiplies the reduced discrete distribution $\hat{\mathbf{f}} \in \mathbb{R}^K$ by the reduced velocity operator $(\Phi^T \text{Diag}(\mathbf{v}) \Phi)$, then decompresses and finally integrates. With such a definition, we have the relations:

$$ \hat{m}_{\rho u} = \hat{m}_{\rho} \hat{A}, \quad \hat{m}_{w} = \frac{1}{2} \hat{m}_{\rho u} \hat{A}. \quad (17) $$

which will be useful in obtaining a good structure of the reduced moments equation.

We modify the reduced collision operator $\hat{Q}_{\Delta v}$ by taking the closest vector $\hat{Q}$ such that $\hat{m}\hat{Q} = 0$. We thus define the corrected collision operator $\hat{Q}_{\Delta v}^c$ as the solution to the minimizing problem:

$$ \hat{Q}_{\Delta v}^c (\Phi \hat{\mathbf{f}}) = \underset{\hat{Q} \in \mathbb{R}^K \text{ s.t. } \hat{m} \hat{Q} = 0}{\operatorname{argmin}} \| \Phi^T \hat{Q} - Q_{\Delta v} (\Phi \hat{\mathbf{f}}) \|^2. $$

The solution to this minimization problem is given by:

$$ \hat{Q}_{\Delta v}^c (\Phi \hat{\mathbf{f}}) = \Phi^T Q_{\Delta v} (\Phi \hat{\mathbf{f}}) - \hat{m}^T (\hat{m} \hat{m}^T)^{-1} \hat{m} \Phi^T Q_{\Delta v} (\Phi \hat{\mathbf{f}}). \quad (18) $$

The proof can be found in Appendix B. Naturally, we can apply the DEIM strategy presented in the previous section to this corrected collision operator. Namely,

$$ \hat{Q}_{\Delta v}^{\text{DEIM}} (\Phi \hat{\mathbf{f}}) = \hat{Q}_{\Delta v}^{\text{DEIM}} (\hat{\mathbf{f}}) - \hat{m}^T (\hat{m} \hat{m}^T)^{-1} \hat{m} \hat{Q}_{\Delta v}^{\text{DEIM}} (\hat{\mathbf{f}}). \quad (19) $$

## 3.3 Reduced Maxwellian distributions

Given a set of moments $\mu \in \mathbb{R}^3$, we denote $\mathbf{M}(\mu) = \mathbf{M}_{\rho,u,T}$ the discrete Maxwellian distribution, where $\rho, u, T$ are defined such that $\mu = (\rho, \rho u, \rho u^2/2 + \rho T/2)^T$. We then defined the reduced Maxwellian $\hat{\mathbf{M}}(\mu)$ associated with reduced moments $\mu$ as the closest to the full Maxwellian $\mathbf{M}(\mu)$ after decompression:

$$ \hat{\mathbf{M}}(\mu) = \underset{\hat{\mathbf{f}} \in \mathbb{R}^K \text{ s.t. } \hat{m}\hat{\mathbf{f}} = \mu}{\operatorname{argmin}} \|\Phi\hat{\mathbf{f}} - \mathbf{M}(\mu)\|^2. \quad (20) $$

As in the previous section, the solution to this minimization problem denoted is given by:

$$ \hat{\mathbf{M}}(\mu) = \Phi^T \mathbf{M}(\mu) + \hat{m}^T (\hat{m} \hat{m}^T)^{-1} (\mu - \hat{m} \Phi^T \mathbf{M}(\mu)). $$

If $\mu$ is associated with the reduced moments of $\hat{\mathbf{f}}$, then we get the following expression:

$$ \hat{\mathbf{M}}(\hat{m}\hat{\mathbf{f}}) = \Phi^T \mathbf{M}(\hat{m}\hat{\mathbf{f}}) + \hat{m}^T (\hat{m}\hat{m}^T)^{-1} (\hat{m}\hat{\mathbf{f}} - \hat{m}\Phi^T \mathbf{M}(\hat{m}\hat{\mathbf{f}})). \quad (21) $$

The reduced collision operator is not guaranteed to vanish for reduced Maxwellians. We thus propose to consider the following second correction:

$$ \hat{Q}_{\Delta v}^{c2} (\Phi \hat{\mathbf{f}}) = \left( \frac{\| \hat{\mathbf{f}} - \hat{\mathbf{M}}(\hat{m}\hat{\mathbf{f}}) \|}{\| \hat{\mathbf{f}} - \hat{\mathbf{M}}(\hat{m}\hat{\mathbf{f}}) \| + \delta} \right) \hat{Q}_{\Delta v}^c (\Phi \hat{\mathbf{f}}), \quad (22) $$
## 3.2 Preservation of reduced moments

The reduced collision operator does not a priori satisfy the preservation of the reduced moments nor vanishes on discrete Maxwellians. We propose to slightly change the discrete collision operator to enforce these two properties, as already proposed in [Rif20]. Similar corrections have been proposed in the context of spectral methods [GT09], discontinuous Galerkin methods [ZG18] or discrete-velocity methods [DL13].

We first introduced the reduced moment operator as follows:

$$ \hat{m} = \begin{bmatrix} \hat{m}_{\rho} \\ \hat{m}_{\rho u} \\ \hat{m}_{w} \end{bmatrix} = \begin{bmatrix} \Delta v \mathbf{1}^T \\ \Delta v \mathbf{1}^T \Phi (\Phi^T \text{Diag}(\mathbf{v}) \Phi) \\ \frac{1}{2} \Delta v \mathbf{1}^T \Phi (\Phi^T \text{Diag}(\mathbf{v}) \Phi)^2 \end{bmatrix} \in M_{3,K}(\mathbb{R}). \quad (16) $$

Note that this definition differs from the moment operator $m\Phi$, which would be the natural definition since it first decompresses up to the full size and then applies the classical moment operator. Instead, the considered reduced moment operator first multiplies the reduced discrete distribution $\hat{\mathbf{f}} \in \mathbb{R}^K$ by the reduced velocity operator $(\Phi^T \text{Diag}(\mathbf{v}) \Phi)$, then decompresses and finally integrates. With such a definition, we have the relations:

$$ \hat{m}_{\rho u} = \hat{m}_{\rho} \hat{A}, \quad \hat{m}_{w} = \frac{1}{2} \hat{m}_{\rho u} \hat{A}. \quad (17) $$

which will be useful in obtaining a good structure of the reduced moments equation.

We modify the reduced collision operator $\hat{Q}_{\Delta v}$ by taking the closest vector $\hat{Q}$ such that $\hat{m}\hat{Q} = 0$. We thus define the corrected collision operator $\hat{Q}_{\Delta v}^c$ as the solution to the minimizing problem:

$$ \hat{Q}_{\Delta v}^c (\Phi \hat{\mathbf{f}}) = \underset{\hat{Q} \in \mathbb{R}^K \text{ s.t. } \hat{m} \hat{Q} = 0}{\operatorname{argmin}} \| \Phi^T \hat{Q} - Q_{\Delta v} (\Phi \hat{\mathbf{f}}) \|^2. $$

The solution to this minimization problem is given by:

$$ \hat{Q}_{\Delta v}^c (\Phi \hat{\mathbf{f}}) = \Phi^T Q_{\Delta v} (\Phi \hat{\mathbf{f}}) - \hat{m}^T (\hat{m} \hat{m}^T)^{-1} \hat{m} \Phi^T Q_{\Delta v} (\Phi \hat{\mathbf{f}}). \quad (18) $$

The proof can be found in Appendix B. Naturally, we can apply the DEIM strategy presented in the previous section to this corrected collision operator. Namely,

$$ \hat{Q}_{\Delta v}^{\text{DEIM}} (\Phi \hat{\mathbf{f}}) = \hat{Q}_{\Delta v}^{\text{DEIM}} (\hat{\mathbf{f}}) - \hat{m}^T (\hat{m} \hat{m}^T)^{-1} \hat{m} \hat{Q}_{\Delta v}^{\text{DEIM}} (\hat{\mathbf{f}}). \quad (19) $$

## 3.3 Reduced Maxwellian distributions

Given a set of moments $\mu \in \mathbb{R}^3$, we denote $\mathbf{M}(\mu) = \mathbf{M}_{\rho,u,T}$ the discrete Maxwellian distribution, where $\rho, u, T$ are defined such that $\mu = (\rho, \rho u, \rho u^2/2 + \rho T/2)^T$. We then defined the reduced Maxwellian $\hat{\mathbf{M}}(\mu)$ associated with reduced moments $\mu$ as the closest to the full Maxwellian $\mathbf{M}(\mu)$ after decompression:

$$ \hat{\mathbf{M}}(\mu) = \underset{\hat{\mathbf{f}} \in \mathbb{R}^K \text{ s.t. } \hat{m}\hat{\mathbf{f}} = \mu}{\operatorname{argmin}} \|\Phi\hat{\mathbf{f}} - \mathbf{M}(\mu)\|^2. \quad (20) $$

As in the previous section, the solution to this minimization problem denoted is given by:

$$ \hat{\mathbf{M}}(\mu) = \Phi^T \mathbf{M}(\mu) + \hat{m}^T (\hat{m} \hat{m}^T)^{-1} (\mu - \hat{m} \Phi^T \mathbf{M}(\mu)). $$

If $\mu$ is associated with the reduced moments of $\hat{\mathbf{f}}$, then we get the following expression:

$$ \hat{\mathbf{M}}(\hat{m}\hat{\mathbf{f}}) = \Phi^T \mathbf{M}(\hat{m}\hat{\mathbf{f}}) + \hat{m}^T (\hat{m}\hat{m}^T)^{-1} (\hat{m}\hat{\mathbf{f}} - \hat{m}\Phi^T \mathbf{M}(\hat{m}\hat{\mathbf{f}})). \quad (21) $$

The reduced collision operator is not guaranteed to vanish for reduced Maxwellians. We thus propose to consider the following second correction:

$$ \hat{Q}_{\Delta v}^{c2} (\Phi \hat{\mathbf{f}}) = \left( \frac{\| \hat{\mathbf{f}} - \hat{\mathbf{M}}(\hat{m}\hat{\mathbf{f}}) \|}{\| \hat{\mathbf{f}} - \hat{\mathbf{M}}(\hat{m}\hat{\mathbf{f}}) \| + \delta} \right) \hat{Q}_{\Delta v}^c (\Phi \hat{\mathbf{f}}), \quad (22) $$
with $\delta > 0$.

As for the collision operator, we can apply the DEIM method to avoid the computation of $M(\mu)$ appearing in (21) and which is evaluated in the full space. We consider samples of the Maxwellian distribution and construct a reduced basis of size $K_M$ and gather it into the matrix $\Phi_M$ of dimension $N_v \times K_M$. Then we consider the following reduced Maxwellian:

$$ \widehat{M}^{\text{DEIM}}(\hat{m}\hat{\mathbf{f}}) = \Phi^T (\Phi_M (P_M^T \Phi_M)^{-1} P_M^T \mathbf{M}(\hat{m}\hat{\mathbf{f}})) \\ + \hat{m}^T (\hat{m}\hat{m}^T)^{-1} (\hat{m}\hat{\mathbf{f}} - \hat{m}\Phi^T (\Phi_M (P_M^T \Phi_M)^{-1} P_M^T \mathbf{M}(\hat{m}\hat{\mathbf{f}}))) $$

where $P_M$ is a matrix of size $N_v \times K_M$ which selects $K_M$ rows of $M$. Then combining the above definition of the reduced Maxwellian with the second corrected reduced collisional operator (23), the reduced collision operator writes:

$$ \hat{Q}_{\Delta v}^{\text{c2DEIM}}(\Phi \hat{\mathbf{f}}) = \left( \frac{\|\hat{\mathbf{f}} - \widehat{M}^{\text{DEIM}}(\hat{m}\hat{\mathbf{f}})\|}{\|\hat{\mathbf{f}} - \widehat{M}^{\text{DEIM}}(\hat{m}\hat{\mathbf{f}})\| + \delta} \right) \hat{Q}_{\Delta v}^{\text{cDEIM}}(\Phi \hat{\mathbf{f}}). \quad (23) $$

## 3.4 Reduced moment equations

In this part, we show that the equation on the reduced moments is an Euler-type system. Indeed, let us consider the semi-discretized reduced equation as

$$ \partial_t \hat{\mathbf{f}} + \hat{A} \partial_x \hat{\mathbf{f}} + (\partial_x \phi) \hat{D} \hat{\mathbf{f}} = \frac{1}{\varepsilon} \hat{Q}_{\Delta v}^c (\Phi \hat{\mathbf{f}}), \quad (24) $$

where $\hat{Q}^{(2)}$ is the corrected reduced collision kernel defined in (23). As this operator has vanishing reduced moments, we have:

$$ \partial_t(\hat{m}\hat{\mathbf{f}}) + \partial_x(\hat{m}\hat{A}\hat{\mathbf{f}}) + (\partial_x\phi)\hat{m}\hat{D}\hat{\mathbf{f}} = 0, $$

which is equivalent to the system of equations:

$$ \left\{ \begin{aligned} \partial_t(\hat{m}_\rho \hat{\mathbf{f}}) + \partial_x(\hat{m}_\rho \hat{A}\hat{\mathbf{f}}) + (\partial_x\phi)\hat{m}_\rho \hat{D}\hat{\mathbf{f}} &= 0, \\ \partial_t(\hat{m}_{\rho u} \hat{\mathbf{f}}) + \partial_x(\hat{m}_{\rho u} \hat{A}\hat{\mathbf{f}}) + (\partial_x\phi)\hat{m}_{\rho u} \hat{D}\hat{\mathbf{f}} &= 0, \\ \partial_t(\hat{m}_w \hat{\mathbf{f}}) + \partial_x(\hat{m}_w \hat{A}\hat{\mathbf{f}}) + (\partial_x\phi)\hat{m}_w \hat{D}\hat{\mathbf{f}} &= 0, \end{aligned} \right. \quad (25) $$

Next, thanks to relations (17), we get the following property:

**Proposition 1.** Denoting $(\hat{\rho}, \hat{\rho}\hat{u}, \hat{w})^T = \hat{m}\hat{\mathbf{f}}$, system (25) becomes the following Euler-Poisson type system:

$$ \left\{ \begin{aligned} \partial_t \hat{\rho} + \partial_x (\hat{\rho}\hat{u}) &= -(\partial_x\phi)\hat{m}_\rho \hat{D}\hat{\mathbf{f}}, \\ \partial_t (\hat{\rho}\hat{u}) + \partial_x (\hat{\rho}\hat{u}^2 + \hat{p}) &= -(\partial_x\phi)\hat{m}_{\rho u} \hat{D}\hat{\mathbf{f}}, \\ \partial_t \hat{w} + \partial_x (\hat{w}\hat{u} + \hat{p}\hat{u} + \hat{q}) &= -(\partial_x\phi)\hat{m}_w \hat{D}\hat{\mathbf{f}}, \end{aligned} \right. \quad (26) $$

where $\hat{p}$ is such that $\hat{w} = \hat{\rho}\hat{u}^2/2 + \hat{p}/2$ and

$$ \hat{q} = \frac{1}{2} \Delta v \mathbf{1}^T \Phi (\hat{A} - \hat{u} \text{Id})^3 \hat{\mathbf{f}}. \quad (27) $$

*Proof.* Indeed, using relations (17), the fluxes of the reduced density and momentum write:

$$ \begin{aligned} \hat{m}_\rho \hat{A} \hat{\mathbf{f}} &= \hat{m}_{\rho u} \hat{A} \hat{\mathbf{f}} = \hat{\rho} \hat{u}, \\ \hat{m}_{\rho u} \hat{A} \hat{\mathbf{f}} &= 2 \hat{m}_w \hat{\mathbf{f}} = 2 \hat{w} = \hat{\rho} \hat{u}^2 + \hat{p}. \end{aligned} $$
with $\delta > 0$.

As for the collision operator, we can apply the DEIM method to avoid the computation of $\mathbf{M}(\mu)$ appearing in (21) and which is evaluated in the full space. We consider samples of the Maxwellian distribution and construct a reduced basis of size $K_M$ and gather it into the matrix $\Phi_M$ of dimension $N_v \times K_M$. Then we consider the following reduced Maxwellian:

$$ \widehat{\mathbf{M}}^{\text{DEIM}}(\hat{m}\hat{\mathbf{f}}) = \Phi^T (\Phi_M (P_M^T \Phi_M)^{-1} P_M^T \mathbf{M}(\hat{m}\hat{\mathbf{f}})) \\ + \hat{m}^T (\hat{m}\hat{m}^T)^{-1} (\hat{m}\hat{\mathbf{f}} - \hat{m}\Phi^T (\Phi_M (P_M^T \Phi_M)^{-1} P_M^T \mathbf{M}(\hat{m}\hat{\mathbf{f}}))) $$

where $P_M$ is a matrix of size $N_v \times K_M$ which selects $K_M$ rows of $\mathbf{M}$. Then combining the above definition of the reduced Maxwellian with the second corrected reduced collisional operator (23), the reduced collision operator writes:

$$ \hat{Q}_{\Delta v}^{\text{c2DEIM}}(\Phi \hat{\mathbf{f}}) = \left( \frac{\|\hat{\mathbf{f}} - \widehat{\mathbf{M}}^{\text{DEIM}}(\hat{m}\hat{\mathbf{f}})\|}{\|\hat{\mathbf{f}} - \widehat{\mathbf{M}}^{\text{DEIM}}(\hat{m}\hat{\mathbf{f}})\| + \delta} \right) \hat{Q}_{\Delta v}^{\text{cDEIM}}(\Phi \hat{\mathbf{f}}). \quad (23) $$

## 3.4 Reduced moment equations

In this part, we show that the equation on the reduced moments is an Euler-type system. Indeed, let us consider the semi-discretized reduced equation as

$$ \partial_t \hat{\mathbf{f}} + \hat{A} \partial_x \hat{\mathbf{f}} + (\partial_x \phi) \hat{D} \hat{\mathbf{f}} = \frac{1}{\varepsilon} \hat{Q}_{\Delta v}^c (\Phi \hat{\mathbf{f}}), \quad (24) $$

where $\hat{Q}^{(2)}$ is the corrected reduced collision kernel defined in (23). As this operator has vanishing reduced moments, we have:

$$ \partial_t(\hat{m}\hat{\mathbf{f}}) + \partial_x(\hat{m}\hat{A}\hat{\mathbf{f}}) + (\partial_x\phi)\hat{m}\hat{D}\hat{\mathbf{f}} = 0, $$

which is equivalent to the system of equations:

$$ \left\{ \begin{aligned} \partial_t(\hat{m}_\rho \hat{\mathbf{f}}) + \partial_x(\hat{m}_\rho \hat{A} \hat{\mathbf{f}}) + (\partial_x \phi) \hat{m}_\rho \hat{D} \hat{\mathbf{f}} &= 0, \\ \partial_t(\hat{m}_{\rho u} \hat{\mathbf{f}}) + \partial_x(\hat{m}_{\rho u} \hat{A} \hat{\mathbf{f}}) + (\partial_x \phi) \hat{m}_{\rho u} \hat{D} \hat{\mathbf{f}} &= 0, \\ \partial_t(\hat{m}_w \hat{\mathbf{f}}) + \partial_x(\hat{m}_w \hat{A} \hat{\mathbf{f}}) + (\partial_x \phi) \hat{m}_w \hat{D} \hat{\mathbf{f}} &= 0, \end{aligned} \right. \quad (25) $$

Next, thanks to relations (17), we get the following property:

**Proposition 1.** Denoting $(\hat{\rho}, \hat{\rho}\hat{u}, \hat{w})^T = \hat{m}\hat{\mathbf{f}}$, system (25) becomes the following Euler-Poisson type system:

$$ \left\{ \begin{aligned} \partial_t \hat{\rho} + \partial_x (\hat{\rho}\hat{u}) &= -(\partial_x \phi) \hat{m}_\rho \hat{D} \hat{\mathbf{f}}, \\ \partial_t (\hat{\rho}\hat{u}) + \partial_x (\hat{\rho}\hat{u}^2 + \hat{p}) &= -(\partial_x \phi) \hat{m}_{\rho u} \hat{D} \hat{\mathbf{f}}, \\ \partial_t \hat{w} + \partial_x (\hat{w}\hat{u} + \hat{p}\hat{u} + \hat{q}) &= -(\partial_x \phi) \hat{m}_w \hat{D} \hat{\mathbf{f}}, \end{aligned} \right. \quad (26) $$

where $\hat{p}$ is such that $\hat{w} = \hat{\rho}\hat{u}^2/2 + \hat{p}/2$ and

$$ \hat{q} = \frac{1}{2} \Delta v \mathbf{1}^T \Phi (\hat{A} - \hat{u} \operatorname{Id})^3 \hat{\mathbf{f}}. \quad (27) $$

*Proof.* Indeed, using relations (17), the fluxes of the reduced density and momentum write:

$$ \begin{aligned} \hat{m}_{\rho} \hat{A} \hat{\mathbf{f}} &= \hat{m}_{\rho u} \hat{A} \hat{\mathbf{f}} = \hat{\rho} \hat{u}, \\ \hat{m}_{\rho u} \hat{A} \hat{\mathbf{f}} &= 2 \hat{m}_w \hat{\mathbf{f}} = 2 \hat{w} = \hat{\rho} \hat{u}^2 + \hat{p}. \end{aligned} $$
Then the flux of the reduced energy is given by:

$$
\begin{align*}
\hat{m}_w \hat{A} \hat{M} &= \frac{1}{2} \Delta v \mathbf{1}^T \Phi \hat{A}^3 \hat{M} \\
&= \frac{1}{2} \Delta v \mathbf{1}^T \Phi (\hat{A} - \hat{u} \mathrm{Id} + \hat{u} \mathrm{Id})^3 \hat{M} \\
&= \frac{1}{2} \Delta v \mathbf{1}^T \Phi \left[ (\hat{A} - \hat{u} \mathrm{Id})^3 + 3\hat{u}(\hat{A} - \hat{u} \mathrm{Id})^2 + 3\hat{u}^2(\hat{A} - \hat{u} \mathrm{Id}) + \hat{u}^3 \mathrm{Id} \right] \hat{M} \\
&= \hat{q} + \frac{3}{2} \hat{u} \hat{p} + 0 + \frac{1}{2} \hat{\rho} \hat{u}^3 \\
&= \hat{q} + \hat{u} \hat{p} + \hat{w} \hat{u},
\end{align*}
$$

where, in the second last equality, we used the definition of $\hat{q}$ given in (27) and the identities:

$$
\begin{align*}
\Delta v \mathbf{1}^T \Phi (\hat{A} - \hat{u} \mathrm{Id}) \hat{f} &= \hat{\rho} \hat{u} - \hat{\rho} \hat{u} = 0, \\
\Delta v \mathbf{1}^T \Phi (\hat{A} - \hat{u} \mathrm{Id})^2 \hat{f} &= \Delta v \mathbf{1}^T \Phi (\hat{A}^2 - 2\hat{u}\hat{A} + \hat{u}^2 \mathrm{Id}) \hat{M} \\
&= 2\hat{w} - 2\hat{\rho}\hat{u}^2 + \hat{\rho}\hat{u}^2 \\
&= \hat{p}.
\end{align*}
$$

Inserting these expressions of fluxes into (25), we obtain (26). □

Consequently, the particular choice of the reduced moments (16) ensures the recovery of an Euler-Poisson type system. Note, however, that the heat flux is a priori non-zero, and there is a priori no conservation of mass.

# 4 Numerical results

We test the reduced order model developed in this paper on the Landau damping test case as considered in [CCL75]. This test consists in considering the following initial distribution:

$$f(0, x, v) = \frac{1}{\sqrt{2\pi}} \exp(-v^2/2) (1 + \alpha \cos(kx)), \quad x \in [0, 2\pi/k], \quad v \in [-v_{\max}, v_{\max}],$$

with $k = 0.5$, $v_{\max} = 6$. We set $\alpha \in [0.01, 0.2]$ and $\varepsilon \in [0.01, 10]$ in order to assess the method on both non-linear Landau damping (large $\alpha$) and collisional regime (small $\varepsilon$). The Vlasov-Poisson-Fokker-Planck model is discretized with $N_x = N_v = 128$.

## 4.1 Reduction for given parameters

In this section, we consider the following parameters: $(\varepsilon, \alpha) \in \{0.01, 0.1, 1, 10\} \times \{0.01, 0.1, 0.2\}$. For each pair $(\varepsilon, \alpha)$, the reduced model is built on $N_s = 200$ uniform samples. Then we will discuss the reduced dimensions $K$, $K_M$, and $K_Q$ needed to recover the right dynamics. The final time is taken equal to $T = 40$ in this test and $\delta = 1 \times 10^{-10}$.

Figure 1 shows the singular value distributions corresponding to the samples matrices $N_s \times N_s$ of the distribution function $f$, the Maxwellian $M(mf)$ and the collision operator $Q(f)$ for the different choices of $(\varepsilon, \alpha)$. We remark that the singular values associated with the Maxwellian (in green dots) decrease fast down to $10^{-11}$ regardless $(\varepsilon, \alpha)$. In red dots, the singular values of the distribution function decrease slowly for large $\alpha$ and are similar to those of the Maxwellians for small $\varepsilon$. Similarly, the singular values for the collision operator (in blue dots) increase with both $\alpha$ and $\varepsilon$. For instance, with $\varepsilon = 10$, they reach a plateau of $10^{-5}$ for $\alpha = 0.01$, $10^{-3}$ for $\alpha = 0.1$ and $10^{-2}$ for $\alpha = 0.2$. Consequently, the model becomes stiffer with larger values of $\alpha$ and $\varepsilon$. Moreover, vertical lines in Figure 1 represent the threshold values in order to obtain correct damping results: only eigenvectors with associated singular values below the threshold are considered in the reduced model. These values are reported in Table 1.

Figure 2 shows the evolution of the $L^2$ norm of the electric field:

$$|E(t, \cdot)|_{L^2} = \sqrt{\int_0^{\frac{2\pi}{k}} (-\partial_x \phi(t, x))^2 dx},$$
Then the flux of the reduced energy is given by:

$$
\begin{align*}
\hat{m}_w \hat{A} \hat{M} &= \frac{1}{2} \Delta v \mathbf{1}^T \Phi \hat{A}^3 \hat{M} \\
&= \frac{1}{2} \Delta v \mathbf{1}^T \Phi (\hat{A} - \hat{u} \mathrm{Id} + \hat{u} \mathrm{Id})^3 \hat{M} \\
&= \frac{1}{2} \Delta v \mathbf{1}^T \Phi \left[ (\hat{A} - \hat{u} \mathrm{Id})^3 + 3 \hat{u} (\hat{A} - \hat{u} \mathrm{Id})^2 + 3 \hat{u}^2 (\hat{A} - \hat{u} \mathrm{Id}) + \hat{u}^3 \mathrm{Id} \right] \hat{M} \\
&= \hat{q} + \frac{3}{2} \hat{u} \hat{p} + 0 + \frac{1}{2} \hat{\rho} \hat{u}^3 \\
&= \hat{q} + \hat{u} \hat{p} + \hat{w} \hat{u},
\end{align*}
$$

where, in the second last equality, we used the definition of $\hat{q}$ given in (27) and the identities:

$$
\begin{align*}
\Delta v \mathbf{1}^T \Phi (\hat{A} - \hat{u} \mathrm{Id}) \hat{f} &= \hat{\rho} \hat{u} - \hat{\rho} \hat{u} = 0, \\
\Delta v \mathbf{1}^T \Phi (\hat{A} - \hat{u} \mathrm{Id})^2 \hat{f} &= \Delta v \mathbf{1}^T \Phi (\hat{A}^2 - 2\hat{u}\hat{A} + \hat{u}^2 \mathrm{Id}) \hat{M} \\
&= 2\hat{w} - 2\hat{\rho}\hat{u}^2 + \hat{\rho}\hat{u}^2 \\
&= \hat{p}.
\end{align*}
$$

Inserting these expressions of fluxes into (25), we obtain (26). □

Consequently, the particular choice of the reduced moments (16) ensures the recovery of an Euler-Poisson type system. Note, however, that the heat flux is a priori non-zero, and there is a priori no conservation of mass.

# 4 Numerical results

We test the reduced order model developed in this paper on the Landau damping test case as considered in [CCL75]. This test consists in considering the following initial distribution:

$$f(0, x, v) = \frac{1}{\sqrt{2\pi}} \exp(-v^2/2) (1 + \alpha \cos(kx)), \quad x \in [0, 2\pi/k], \quad v \in [-v_{\max}, v_{\max}],$$

with $k = 0.5$, $v_{\max} = 6$. We set $\alpha \in [0.01, 0.2]$ and $\varepsilon \in [0.01, 10]$ in order to assess the method on both non-linear Landau damping (large $\alpha$) and collisional regime (small $\varepsilon$). The Vlasov-Poisson-Fokker-Planck model is discretized with $N_x = N_v = 128$.

## 4.1 Reduction for given parameters

In this section, we consider the following parameters: $(\varepsilon, \alpha) \in \{0.01, 0.1, 1, 10\} \times \{0.01, 0.1, 0.2\}$. For each pair $(\varepsilon, \alpha)$, the reduced model is built on $N_s = 200$ uniform samples. Then we will discuss the reduced dimensions $K$, $K_M$, and $K_Q$ needed to recover the right dynamics. The final time is taken equal to $T = 40$ in this test and $\delta = 1 \times 10^{-10}$.

Figure 1 shows the singular value distributions corresponding to the samples matrices $N_s \times N_s$ of the distribution function $f$, the Maxwellian $M(mf)$ and the collision operator $Q(f)$ for the different choices of $(\varepsilon, \alpha)$. We remark that the singular values associated with the Maxwellian (in green dots) decrease fast down to $10^{-11}$ regardless $(\varepsilon, \alpha)$. In red dots, the singular values of the distribution function decrease slowly for large $\alpha$ and are similar to those of the Maxwellians for small $\varepsilon$. Similarly, the singular values for the collision operator (in blue dots) increase with both $\alpha$ and $\varepsilon$. For instance, with $\varepsilon = 10$, they reach a plateau of $10^{-5}$ for $\alpha = 0.01$, $10^{-3}$ for $\alpha = 0.1$ and $10^{-2}$ for $\alpha = 0.2$. Consequently, the model becomes stiffer with larger values of $\alpha$ and $\varepsilon$. Moreover, vertical lines in Figure 1 represent the threshold values in order to obtain correct damping results: only eigenvectors with associated singular values below the threshold are considered in the reduced model. These values are reported in Table 1.

Figure 2 shows the evolution of the $L^2$ norm of the electric field:

$$|E(t, \cdot)|_{L^2} = \sqrt{\int_0^{\frac{2\pi}{k}} (-\partial_x \phi(t, x))^2 dx},$$
for some pairs $(\epsilon, \alpha)$ obtained with the reduced model, in log scale, and obtained with both the reduced and the kinetic models. Comparing the results obtained with the reduced model to those obtained with the kinetic model, we observe that a good damping rate is recovered using the reduced model. Also, we check that the DEIM and corrected moments do not deteriorate the results.

In conclusion, the reduced model provides good results for various ranges of parameters $\epsilon$ and $\alpha$ and enables a drastic reduction of unknowns.

<figure><img src="image_2.png" /><figcaption>Figure 1: Singular values $(\sigma_i)_i$ of the samples matrix for the distribution function $f$, the Maxwellian $M$ and the collision operator $Q$.</figcaption></figure>
for some pairs $(\epsilon, \alpha)$ obtained with the reduced model, in log scale, and obtained with both the reduced and the kinetic models. Comparing the results obtained with the reduced model to those obtained with the kinetic model, we observe that a good damping rate is recovered using the reduced model. Also, we check that the DEIM and corrected moments do not deteriorate the results.

In conclusion, the reduced model provides good results for various ranges of parameters $\epsilon$ and $\alpha$ and enables a drastic reduction of unknowns.

<figure><img src="image_2.png" /><figcaption>Figure 1: Singular values $(\sigma_i)_i$ of the samples matrix for the distribution function $f$, the Maxwellian $M$ and the collision operator $Q$.</figcaption></figure>
<table><caption>Table 1: Number of eigenvectors kept for the solution space (K), the Maxwellian (K<sub>M</sub>) and the collision operator (K<sub>Q</sub>).</caption><thead><tr><th>α</th><th>0.01</th><th>0.1</th><th>0.2</th></tr></thead><tbody><tr><th>ε</th><td></td><td></td><td></td></tr><tr><th>0.01</th><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (6, 6, 6)</td><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (9, 9, 9)</td><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (15, 10, 10)</td></tr><tr><th>0.1</th><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (5, 3, 3)</td><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (9, 5, 5)</td><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (10, 7, 7)</td></tr><tr><th>1</th><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (7, 7, 7)</td><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (7, 7, 7)</td><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (7, 7, 7)</td></tr><tr><th>10</th><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (11, 11, 11)</td><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (11, 11, 11)</td><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (11, 11, 11)</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2: Logarithm of the L<sup>2</sup> norm of the electric field |E(t, .)|<sub>L<sup>2</sup></sub> as a function of time for the reference solution, reduced model, with DEIM and DEIM with corrected moments.</figcaption></figure>

## 4.2 Preservation of Maxwellian distributions

In Section 3.3, we proposed a second correction of the reduced collisional operator (see Eq. (23)), which forces it to vanish for reduced Maxwellians thanks to a prefactor. The correction is parameterized by δ ≥ 0. This correction is aimed at preserving Maxwellian distributions in time, up to the compression-decompression error. To assess the effectiveness of this correction, we consider a homogeneous test case with a Maxwellian initial distribution:

$$ \partial_t f(t, x, v) = \frac{1}{\varepsilon} Q(f)(t, x, v), \quad x \in [0, 2\pi/k], \ v \in [-v_{\max}, v_{\max}], \\ f(0, x, v) = \frac{1}{\sqrt{2\pi}} e^{-\frac{v^2}{2}}, $$

with the parameter ε ∈ {0.01, 0.1, 1, 10} on the time interval [0, T] with T = 20. A specific reduced model is built for each ε using N<sub>s</sub> = 200 uniform samples and reduced dimensions (K, K<sub>Q</sub>, K<sub>M</sub>) = (11, 11, 11). We then compare the solution obtained using the reduced models at time T with the exact solution f(T, x, v) = f(0, x, v), when using different values of the parameter δ involved in the correction: from δ = 0, which corresponds to no correction, to δ = 1 × 10<sup>-2</sup>. In Figure 3, we observe the L<sup>2</sup> error for the four different reduced models built for each ε. While for ε = 10, the correction improves the result by a factor 5 only, the error can be reduced by a factor 1000 when ε = 1. For a smaller ε = 0.1, the second correction is even more important as the reduced model without correction (δ = 0) leads to a large error of order 5 × 10<sup>-1</sup>, while the correction reduces the error to about 2 × 10<sup>-8</sup> for δ larger than 1 × 10<sup>-7</sup>. Additionally, for ε = 0.01, when δ values range between 0 and 10<sup>-5</sup>, the error remains constant and relatively high. Conversely, the error improves as delta exceeds 10<sup>-4</sup>.
<table><caption>Table 1: Number of eigenvectors kept for the solution space (K), the Maxwellian (K<sub>M</sub>) and the collision operator (K<sub>Q</sub>).</caption><thead><tr><th>α</th><th>0.01</th><th>0.1</th><th>0.2</th></tr></thead><tbody><tr><th>ε</th><td></td><td></td><td></td></tr><tr><th>0.01</th><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (6, 6, 6)</td><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (9, 9, 9)</td><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (15, 10, 10)</td></tr><tr><th>0.1</th><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (5, 3, 3)</td><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (9, 5, 5)</td><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (10, 7, 7)</td></tr><tr><th>1</th><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (7, 7, 7)</td><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (7, 7, 7)</td><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (7, 7, 7)</td></tr><tr><th>10</th><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (11, 11, 11)</td><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (11, 11, 11)</td><td>(K, K<sub>Q</sub>, K<sub>M</sub>) = (11, 11, 11)</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2: Logarithm of the L<sup>2</sup> norm of the electric field |E(t, .)|<sub>L<sup>2</sup></sub> as a function of time for the reference solution, reduced model, with DEIM and DEIM with corrected moments.</figcaption></figure>

## 4.2 Preservation of Maxwellian distributions

In Section 3.3, we proposed a second correction of the reduced collisional operator (see Eq. (23)), which forces it to vanish for reduced Maxwellians thanks to a prefactor. The correction is parameterized by δ ≥ 0. This correction is aimed at preserving Maxwellian distributions in time, up to the compression-decompression error. To assess the effectiveness of this correction, we consider a homogeneous test case with a Maxwellian initial distribution:

$$ \partial_t f(t, x, v) = \frac{1}{\varepsilon} Q(f)(t, x, v), \quad x \in [0, 2\pi/k], \ v \in [-v_{\max}, v_{\max}], \\ f(0, x, v) = \frac{1}{\sqrt{2\pi}} e^{-\frac{v^2}{2}}, $$

with the parameter ε ∈ {0.01, 0.1, 1, 10} on the time interval [0, T] with T = 20. A specific reduced model is built for each ε using N<sub>s</sub> = 200 uniform samples and reduced dimensions (K, K<sub>Q</sub>, K<sub>M</sub>) = (11, 11, 11). We then compare the solution obtained using the reduced models at time T with the exact solution f(T, x, v) = f(0, x, v), when using different values of the parameter δ involved in the correction: from δ = 0, which corresponds to no correction, to δ = 1 × 10<sup>-2</sup>. In Figure 3, we observe the L<sup>2</sup> error for the four different reduced models built for each ε. While for ε = 10, the correction improves the result by a factor 5 only, the error can be reduced by a factor 1000 when ε = 1. For a smaller ε = 0.1, the second correction is even more important as the reduced model without correction (δ = 0) leads to a large error of order 5 × 10<sup>-1</sup>, while the correction reduces the error to about 2 × 10<sup>-8</sup> for δ larger than 1 × 10<sup>-7</sup>. Additionally, for ε = 0.01, when δ values range between 0 and 10<sup>-5</sup>, the error remains constant and relatively high. Conversely, the error improves as delta exceeds 10<sup>-4</sup>.
<figure><img src="image_4.png" /><figcaption>Figure 3: $L^2$ error of the DEIM with second correction reduced model as a function of $\delta$</figcaption></figure>

## 4.3 Generalization to other parameters

In Section 4.1, a different reduced model has been constructed for each set of parameters. Here, we explore the ability to construct a reduced model valid in a full range of the parameters. This generalisation property is a key feature: the reduced model could then be used for parameters that are not involved in the construction of the reduced model, without the need to carry out the full kinetic simulations. Here, we would like to build a reduced model valid for all the values $(\epsilon, \alpha)$ in the domain $\mathcal{D} = [1, 10] \times [0.01, 0.1]$. Therefore, we consider 25 pairs $(\epsilon, \alpha)$, randomly chosen in $\mathcal{D}$ and the reduced models is built using 30 samples per pairs $(\epsilon, \alpha)$, resulting into $N_s = 750$ samples. The final time equals $T = 20$. We set $\delta = 1 \times 10^{-10}$.

According to singular value distributions of the sample matrices (see Figure 4), we set $K = K_M = K_Q = 10$. In Figure 4, on the right, is depicted the electric energy field for $(\epsilon, \alpha) = (5, 5 \times 10^{-2}) \in \mathcal{D}$. We observe that the dynamics are well recovered for each algorithm (with/without DEIM or corrected moments). We also compute the relative $L^2$ error between the solution of the reduced model after decompression $\Phi \hat{f}$ and the solution of the underlying kinetic model $f$ at the final time $T$. Furthermore we compute the $L^2$ norm of the electric field over the time interval $[0, T]$ in Figure 2. We observe that each algorithm performs well. Moreover, the error is larger for the corrected DEIM ROM method.

<figure><img src="image_5.png" /><figcaption>Figure 4: Singular values distributions (left) and logarithm of the norm of the electric field $|E(t, \cdot)|_{L^2}$ as a function of time for $(\epsilon, \alpha) = (5, 5 \times 10^{-2}) \in \mathcal{D}$.</figcaption></figure>

<table><caption>Table 2: relative $L^2$ errors, $(\epsilon, \alpha) = (5, 5 \times 10^{-2}) \in \mathcal{D}$.</caption><thead><tr><th></th><th>ROM</th><th>DEIM</th><th>corrected DEIM</th></tr></thead><tbody><tr><td>$\Phi \hat{f}(T)$</td><td>$5.36 \times 10^{-4}$</td><td>$5.62 \times 10^{-4}$</td><td>$5.66 \times 10^{-4}$</td></tr><tr><td>$|E(t, \cdot)|_{L^2}$</td><td>$1.73 \times 10^{-2}$</td><td>$1.85 \times 10^{-2}$</td><td>$2.63 \times 10^{-2}$</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>Figure 3: $L^2$ error of the DEIM with second correction reduced model as a function of $\delta$</figcaption></figure>

## 4.3 Generalization to other parameters

In Section 4.1, a different reduced model has been constructed for each set of parameters. Here, we explore the ability to construct a reduced model valid in a full range of the parameters. This generalisation property is a key feature: the reduced model could then be used for parameters that are not involved in the construction of the reduced model, without the need to carry out the full kinetic simulations. Here, we would like to build a reduced model valid for all the values $(\varepsilon, \alpha)$ in the domain $\mathcal{D} = [1, 10] \times [0.01, 0.1]$. Therefore, we consider 25 pairs $(\varepsilon, \alpha)$, randomly chosen in $\mathcal{D}$ and the reduced models is built using 30 samples per pairs $(\varepsilon, \alpha)$, resulting into $N_s = 750$ samples. The final time equals $T = 20$. We set $\delta = 1 \times 10^{-10}$.

According to singular value distributions of the sample matrices (see Figure 4), we set $K = K_M = K_Q = 10$. In Figure 4, on the right, is depicted the electric energy field for $(\varepsilon, \alpha) = (5, 5 \times 10^{-2}) \in \mathcal{D}$. We observe that the dynamics are well recovered for each algorithm (with/without DEIM or corrected moments). We also compute the relative $L^2$ error between the solution of the reduced model after decompression $\Phi \hat{f}$ and the solution of the underlying kinetic model $f$ at the final time $T$. Furthermore we compute the $L^2$ norm of the electric field over the time interval $[0, T]$ in Figure 2. We observe that each algorithm performs well. Moreover, the error is larger for the corrected DEIM ROM method.

<figure><img src="image_5.png" /><figcaption>Figure 4: Singular values distributions (left) and logarithm of the norm of the electric field $|E(t, \cdot)|_{L^2}$ as a function of time for $(\varepsilon, \alpha) = (5, 5 \times 10^{-2}) \in \mathcal{D}$.</figcaption></figure>

<table><caption>Table 2: relative $L^2$ errors, $(\varepsilon, \alpha) = (5, 5 \times 10^{-2}) \in \mathcal{D}$.</caption><thead><tr><th></th><th>ROM</th><th>DEIM</th><th>corrected DEIM</th></tr></thead><tbody><tr><td>$\Phi \hat{f}(T)$</td><td>$5.36 \times 10^{-4}$</td><td>$5.62 \times 10^{-4}$</td><td>$5.66 \times 10^{-4}$</td></tr><tr><td>$|E(t, \cdot)|_{L^2}$</td><td>$1.73 \times 10^{-2}$</td><td>$1.85 \times 10^{-2}$</td><td>$2.63 \times 10^{-2}$</td></tr></tbody></table>
In addition, we aim to test the reduced model outside its learned domain. We set a final time of simulation $T_{test} = 25 > T$ and consider $(\varepsilon, \alpha) \in D_{test} \subset D$. We choose $D_{test} = \{0.1, 0.5, 1, 10, 15, 20\} \times \{0.2, 0.3\}$ to have more non-linear damping with both collisional and non-collisional regime. The obtained relative $L^2$ errors on $\Phi\hat{f}$ at time $T_{test}$ and on the norm of the electric field on the time interval $[0, T_{test}]$ are given in Table 3. Overall, each ROM method provides good results with $(t, (\alpha, \varepsilon)) \in [0, T] \times D$ and shows good generalisation performance with $(t, (\alpha, \varepsilon)) \in [T, T_{test}] \times D_{test}$. For instance, even with $\varepsilon = 0.1$ and 20, solution errors remain bounded by $2 \times 10^{-1}$, $2 \times 10^{-2}$ respectively. Going outside the training set, the precision of the models decrease slowly and we still have correctness on both solutions and electric energies.

Going further, we notice that our corrected DEIM ROM performs better when $\varepsilon < 1$. That is due to the fact that its moments are corrected, such as to obtain an Euler-type system. For example, let us observe the energy damping with $(\varepsilon, \alpha) = (0.1, 0.3)$ on Figure 5a. While the ROM and DEIM ROM drift from the reference solution, the corrected DEIM ROM remains close to it. On Figure 5b, we observe that it is no more the case with larger $\varepsilon$.

In conclusion, we have shown that our corrected moments DEIM performs better for small Knudsen numbers $\varepsilon \le 1$ and is better to be used instead of the usual DEIM. Conversely, this correction is not effective with large $\varepsilon > 1$ and the latter is better to be used.

<table><thead><tr><th rowspan="2">α</th><th rowspan="2">Method</th><th colspan="6">ε</th></tr><tr><th>0.1</th><th>0.5</th><th>1</th><th>10</th><th>15</th><th>20</th></tr></thead><tbody><tr><td rowspan="3">0.2</td><td>ROM</td><td>1.41 × 10<sup>-1</sup></td><td>5.58 × 10<sup>-2</sup></td><td>3.55 × 10<sup>-2</sup></td><td>6.02 × 10<sup>-3</sup></td><td>4.44 × 10<sup>-1</sup></td><td>6.89 × 10<sup>-3</sup></td></tr><tr><td>DEIM</td><td>1.36 × 10<sup>-1</sup></td><td>5.56 × 10<sup>-2</sup></td><td>3.59 × 10<sup>-2</sup></td><td>6.26 × 10<sup>-3</sup></td><td>5.48 × 10<sup>-3</sup></td><td>8.20 × 10<sup>-3</sup></td></tr><tr><td>cDEIM</td><td>1.53 × 10<sup>-2</sup></td><td>1.42 × 10<sup>-2</sup></td><td>1.34 × 10<sup>-2</sup></td><td>6.09 × 10<sup>-3</sup></td><td>5.93 × 10<sup>-3</sup></td><td>7.30 × 10<sup>-3</sup></td></tr><tr><td rowspan="3">0.3</td><td>ROM</td><td>1.94 × 10<sup>-1</sup></td><td>1.18 × 10<sup>-1</sup></td><td>9.62 × 10<sup>-2</sup></td><td>2.15 × 10<sup>-2</sup></td><td>1.75 × 10<sup>-2</sup></td><td>1.62 × 10<sup>-2</sup></td></tr><tr><td>DEIM</td><td>1.98 × 10<sup>-1</sup></td><td>1.21 × 10<sup>-1</sup></td><td>9.84 × 10<sup>-2</sup></td><td>2.23 × 10<sup>-2</sup></td><td>1.97 × 10<sup>-2</sup></td><td>1.95 × 10<sup>-2</sup></td></tr><tr><td>cDEIM</td><td>3.99 × 10<sup>-2</sup></td><td>4.78 × 10<sup>-2</sup></td><td>3.93 × 10<sup>-2</sup></td><td>1.55 × 10<sup>-2</sup></td><td>1.57 × 10<sup>-2</sup></td><td>1.64 × 10<sup>-2</sup></td></tr></tbody></table>

(a) relative errors on $\Phi\hat{f}(T_{test})$.

<table><thead><tr><th rowspan="2">α</th><th rowspan="2">Method</th><th colspan="6">ε</th></tr><tr><th>0.1</th><th>0.5</th><th>1</th><th>10</th><th>15</th><th>20</th></tr></thead><tbody><tr><td rowspan="3">0.2</td><td>ROM</td><td>5.98 × 10<sup>-1</sup></td><td>3.30 × 10<sup>-1</sup></td><td>2.15 × 10<sup>-1</sup></td><td>2.68 × 10<sup>-2</sup></td><td>5.75 × 10<sup>-2</sup></td><td>6.96 × 10<sup>-2</sup></td></tr><tr><td>DEIM</td><td>5.79 × 10<sup>-1</sup></td><td>3.27 × 10<sup>-1</sup></td><td>2.12 × 10<sup>-1</sup></td><td>2.69 × 10<sup>-2</sup></td><td>6.52 × 10<sup>-2</sup></td><td>7.92 × 10<sup>-2</sup></td></tr><tr><td>cDEIM</td><td>5.42 × 10<sup>-2</sup></td><td>7.35 × 10<sup>-2</sup></td><td>7.83 × 10<sup>-2</sup></td><td>6.75 × 10<sup>-2</sup></td><td>8.92 × 10<sup>-2</sup></td><td>9.30 × 10<sup>-2</sup></td></tr><tr><td rowspan="3">0.3</td><td>ROM</td><td>6.37 × 10<sup>-1</sup></td><td>4.77 × 10<sup>-1</sup></td><td>3.65 × 10<sup>-1</sup></td><td>6.32 × 10<sup>-2</sup></td><td>5.88 × 10<sup>-2</sup></td><td>6.18 × 10<sup>-2</sup></td></tr><tr><td>DEIM</td><td>6.29 × 10<sup>-1</sup></td><td>4.80 × 10<sup>-1</sup></td><td>3.65 × 10<sup>-1</sup></td><td>6.14 × 10<sup>-2</sup></td><td>5.79 × 10<sup>-2</sup></td><td>6.36 × 10<sup>-2</sup></td></tr><tr><td>cDEIM</td><td>1.03 × 10<sup>-1</sup></td><td>1.49 × 10<sup>-1</sup></td><td>1.54 × 10<sup>-1</sup></td><td>1.05 × 10<sup>-1</sup></td><td>1.05 × 10<sup>-1</sup></td><td>1.03 × 10<sup>-1</sup></td></tr></tbody></table>

(b) relative errors on $|E(t, \cdot)|_{L^2}$ over $[0, T_{test}]$.

<caption>Table 3: relative $L^2$ errors, $(\alpha, \varepsilon) \in D_{test}$.</caption>
In addition, we aim to test the reduced model outside its learned domain. We set a final time of simulation $T_{test} = 25 > T$ and consider $(\varepsilon, \alpha) \in D_{test} \subset D$. We choose $D_{test} = \{0.1, 0.5, 1, 10, 15, 20\} \times \{0.2, 0.3\}$ to have more non-linear damping with both collisional and non-collisional regime. The obtained relative $L^2$ errors on $\Phi\hat{f}$ at time $T_{test}$ and on the norm of the electric field on the time interval $[0, T_{test}]$ are given in Table 3. Overall, each ROM method provides good results with $(t, (\alpha, \varepsilon)) \in [0, T] \times D$ and shows good generalisation performance with $(t, (\alpha, \varepsilon)) \in [T, T_{test}] \times D_{test}$. For instance, even with $\varepsilon = 0.1$ and 20, solution errors remain bounded by $2 \times 10^{-1}$, $2 \times 10^{-2}$ respectively. Going outside the training set, the precision of the models decrease slowly and we still have correctness on both solutions and electric energies.

Going further, we notice that our corrected DEIM ROM performs better when $\varepsilon < 1$. That is due to the fact that its moments are corrected, such as to obtain an Euler-type system. For example, let us observe the energy damping with $(\varepsilon, \alpha) = (0.1, 0.3)$ on Figure 5a. While the ROM and DEIM ROM drift from the reference solution, the corrected DEIM ROM remains close to it. On Figure 5b, we observe that it is no more the case with larger $\varepsilon$.

In conclusion, we have shown that our corrected moments DEIM performs better for small Knudsen numbers $\varepsilon \le 1$ and is better to be used instead of the usual DEIM. Conversely, this correction is not effective with large $\varepsilon > 1$ and the latter is better to be used.

<table><thead><tr><th rowspan="2">α</th><th rowspan="2">Method</th><th colspan="6">ε</th></tr><tr><th>0.1</th><th>0.5</th><th>1</th><th>10</th><th>15</th><th>20</th></tr></thead><tbody><tr><td rowspan="3">0.2</td><td>ROM</td><td>1.41 × 10<sup>-1</sup></td><td>5.58 × 10<sup>-2</sup></td><td>3.55 × 10<sup>-2</sup></td><td>6.02 × 10<sup>-3</sup></td><td>4.44 × 10<sup>-1</sup></td><td>6.89 × 10<sup>-3</sup></td></tr><tr><td>DEIM</td><td>1.36 × 10<sup>-1</sup></td><td>5.56 × 10<sup>-2</sup></td><td>3.59 × 10<sup>-2</sup></td><td>6.26 × 10<sup>-3</sup></td><td>5.48 × 10<sup>-3</sup></td><td>8.20 × 10<sup>-3</sup></td></tr><tr><td>cDEIM</td><td>1.53 × 10<sup>-2</sup></td><td>1.42 × 10<sup>-2</sup></td><td>1.34 × 10<sup>-2</sup></td><td>6.09 × 10<sup>-3</sup></td><td>5.93 × 10<sup>-3</sup></td><td>7.30 × 10<sup>-3</sup></td></tr><tr><td rowspan="3">0.3</td><td>ROM</td><td>1.94 × 10<sup>-1</sup></td><td>1.18 × 10<sup>-1</sup></td><td>9.62 × 10<sup>-2</sup></td><td>2.15 × 10<sup>-2</sup></td><td>1.75 × 10<sup>-2</sup></td><td>1.62 × 10<sup>-2</sup></td></tr><tr><td>DEIM</td><td>1.98 × 10<sup>-1</sup></td><td>1.21 × 10<sup>-1</sup></td><td>9.84 × 10<sup>-2</sup></td><td>2.23 × 10<sup>-2</sup></td><td>1.97 × 10<sup>-2</sup></td><td>1.95 × 10<sup>-2</sup></td></tr><tr><td>cDEIM</td><td>3.99 × 10<sup>-2</sup></td><td>4.78 × 10<sup>-2</sup></td><td>3.93 × 10<sup>-2</sup></td><td>1.55 × 10<sup>-2</sup></td><td>1.57 × 10<sup>-2</sup></td><td>1.64 × 10<sup>-2</sup></td></tr></tbody></table>

(a) relative errors on $\Phi\hat{f}(T_{test})$.

<table><thead><tr><th rowspan="2">α</th><th rowspan="2">Method</th><th colspan="6">ε</th></tr><tr><th>0.1</th><th>0.5</th><th>1</th><th>10</th><th>15</th><th>20</th></tr></thead><tbody><tr><td rowspan="3">0.2</td><td>ROM</td><td>5.98 × 10<sup>-1</sup></td><td>3.30 × 10<sup>-1</sup></td><td>2.15 × 10<sup>-1</sup></td><td>2.68 × 10<sup>-2</sup></td><td>5.75 × 10<sup>-2</sup></td><td>6.96 × 10<sup>-2</sup></td></tr><tr><td>DEIM</td><td>5.79 × 10<sup>-1</sup></td><td>3.27 × 10<sup>-1</sup></td><td>2.12 × 10<sup>-1</sup></td><td>2.69 × 10<sup>-2</sup></td><td>6.52 × 10<sup>-2</sup></td><td>7.92 × 10<sup>-2</sup></td></tr><tr><td>cDEIM</td><td>5.42 × 10<sup>-2</sup></td><td>7.35 × 10<sup>-2</sup></td><td>7.83 × 10<sup>-2</sup></td><td>6.75 × 10<sup>-2</sup></td><td>8.92 × 10<sup>-2</sup></td><td>9.30 × 10<sup>-2</sup></td></tr><tr><td rowspan="3">0.3</td><td>ROM</td><td>6.37 × 10<sup>-1</sup></td><td>4.77 × 10<sup>-1</sup></td><td>3.65 × 10<sup>-1</sup></td><td>6.32 × 10<sup>-2</sup></td><td>5.88 × 10<sup>-2</sup></td><td>6.18 × 10<sup>-2</sup></td></tr><tr><td>DEIM</td><td>6.29 × 10<sup>-1</sup></td><td>4.80 × 10<sup>-1</sup></td><td>3.65 × 10<sup>-1</sup></td><td>6.14 × 10<sup>-2</sup></td><td>5.79 × 10<sup>-2</sup></td><td>6.36 × 10<sup>-2</sup></td></tr><tr><td>cDEIM</td><td>1.03 × 10<sup>-1</sup></td><td>1.49 × 10<sup>-1</sup></td><td>1.54 × 10<sup>-1</sup></td><td>1.05 × 10<sup>-1</sup></td><td>1.05 × 10<sup>-1</sup></td><td>1.03 × 10<sup>-1</sup></td></tr></tbody></table>

(b) relative errors on $|E(t, \cdot)|_{L^2}$ over $[0, T_{test}]$.

<caption>Table 3: relative $L^2$ errors, $(\alpha, \varepsilon) \in D_{test}$.</caption>
<figure><img src="image_6.png" /><figcaption>Figure 5: Logarithm of the norm of the electric field as a function of time for several $(t, \mu) \in [0, T_{test}] \times D_{test}^{\mu}$.</figcaption></figure>

# 5 Conclusion

In this paper, we have proposed a reduced order modeling approach in velocity applied to the Vlasov-Poisson-Fokker-Planck equation. This results in a hyperbolic system with source terms approximating the dynamics for a given initial data and a given Knudsen number or a range of initial data and Knudsen numbers. The numerical results show that both collisional and non-collisional regimes can be approximated with a reduced system of size ≈ 10. We also note that we can consider a smaller reduced system in a collisional regime. Inversely, the system size should be larger to capture non-linear dynamics. We have introduced a corrected reduced collision operator that preserves moments. In addition to ensuring some mathematical structure of the reduced moment equations, it numerically provides better results in collisional regimes. Finally, to extend the generalization range of the method, we would have to turn to non-linear reduction methods like quadratic ones [RSK23] or auto-encoder reduction [KK20].

# References

[BDS01] Christophe Buet, Stéphane Dellacherie, and Rémi Sentis. Numerical solution of an ionic fokker-planck equation with electronic temperature. *SIAM Journal on Numerical Analysis*, 39(4):1219–1253, 2001.

[CA22] Nicolas Crouseilles and Ibrahim Almuslimani. Conservative stabilized runge-kutta methods for the vlasov-fokker-planck equation. *preprint. hal-03911417*, 2022.

[CCL75] Anaïs Crestetto, Nicolas Crouseilles, and Mohammed Lemou. Kinetic/fluid micro-macro numerical schemes for vlasov-poisson-bgk equation using particles. *Kinetic and Related Models*, AIMS, 5(4):787–816, 2012, ff10.3934/krm.2012.5.787ff. ffhal-00728875.

[CS10] Saifon Chaturantabut and Danny C. Sorensen. Nonlinear model reduction via Discrete Empirical Interpolation. *SIAM Journal on Scientific Computing*, 32(5):2737–2764, 2010.

[DL13] Giacomo Dimarco and Raphaël Loubere. Towards an ultra efficient kinetic scheme. part ii: The high order case. *Journal of computational physics*, 255:699–719, 2013.

[FL22] Cassini Fabio and Einkemmer Lukas. Efficient 6D Vlasov simulation using the dynamical low-rank framework Ensign. *Computer Physics Communications*, 280:108489, 2022.

[GEQ23] Wei Guo, Jannatul Ferdous Ema, and Jing-Mei Qiu. A local macroscopic conservative (lomac) low rank tensor method with the discontinuous galerkin method for the vlasov dynamics. *Communications on Applied Mathematics and Computation*, pages 1–26, 2023.
<figure><img src="image_6.png" /><figcaption>Figure 5: Logarithm of the norm of the electric field as a function of time for several $(t, \mu) \in [0, T_{test}] \times D_{test}^{\mu}$.</figcaption></figure>

# 5 Conclusion

In this paper, we have proposed a reduced order modeling approach in velocity applied to the Vlasov-Poisson-Fokker-Planck equation. This results in a hyperbolic system with source terms approximating the dynamics for a given initial data and a given Knudsen number or a range of initial data and Knudsen numbers. The numerical results show that both collisional and non-collisional regimes can be approximated with a reduced system of size ≈ 10. We also note that we can consider a smaller reduced system in a collisional regime. Inversely, the system size should be larger to capture non-linear dynamics. We have introduced a corrected reduced collision operator that preserves moments. In addition to ensuring some mathematical structure of the reduced moment equations, it numerically provides better results in collisional regimes. Finally, to extend the generalization range of the method, we would have to turn to non-linear reduction methods like quadratic ones [RSK23] or auto-encoder reduction [KK20].

# References

[BDS01] Christophe Buet, Stéphane Dellacherie, and Rémi Sentis. Numerical solution of an ionic fokker-planck equation with electronic temperature. *SIAM Journal on Numerical Analysis*, 39(4):1219–1253, 2001.

[CA22] Nicolas Crouseilles and Ibrahim Almuslimani. Conservative stabilized runge-kutta methods for the vlasov-fokker-planck equation. *preprint. hal-03911417*, 2022.

[CCL75] Anaïs Crestetto, Nicolas Crouseilles, and Mohammed Lemou. Kinetic/fluid micro-macro numerical schemes for vlasov-poisson-bgk equation using particles. *Kinetic and Related Models*, AIMS, 5(4):787–816, 2012, ff10.3934/krm.2012.5.787ff. ffhal-00728875.

[CS10] Saifon Chaturantabut and Danny C. Sorensen. Nonlinear model reduction via Discrete Empirical Interpolation. *SIAM Journal on Scientific Computing*, 32(5):2737–2764, 2010.

[DL13] Giacomo Dimarco and Raphaël Loubere. Towards an ultra efficient kinetic scheme. part ii: The high order case. *Journal of computational physics*, 255:699–719, 2013.

[FL22] Cassini Fabio and Einkemmer Lukas. Efficient 6D Vlasov simulation using the dynamical low-rank framework Ensign. *Computer Physics Communications*, 280:108489, 2022.

[GEQ23] Wei Guo, Jannatul Ferdous Ema, and Jing-Mei Qiu. A local macroscopic conservative (lomac) low rank tensor method with the discontinuous galerkin method for the vlasov dynamics. *Communications on Applied Mathematics and Computation*, pages 1–26, 2023.
[GT09] Irene M Gamba and Sri Harsha Tharkabhushanam. Spectral-lagrangian methods for collisional models of non-equilibrium statistical states. *Journal of Computational Physics*, 228(6):2012–2036, 2009.

[HPR22a] Jan S. Hesthaven, Cecilia Pagliantini, and Nicolò Ripamonti. Adaptive symplectic model order reduction of parametric particle-based Vlasov-Poisson equation, 2022.

[HPR22b] Jan S. Hesthaven, Cecilia Pagliantini, and Gianluigi Rozza. Reduced basis methods for time-dependent problems. *Acta Numerica*, 31:265–345, 2022.

[KK20] Lee Kookjin and Carlberg Kevin, T. Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders. *Journal of Computational Physics*, 404:108973, 2020.

[Rif20] Sébastien Riffaud. *Reduced-order models: convergence between scientific computing and data for fluid mechanics*. PhD thesis, Université de Bordeaux, 2020.

[RSK23] Geelen Rudy, Wright Stephen, and Willcox Karen. Operator inference for non-intrusive model reduction with quadratic manifolds. *Computer Methods in Applied Mechanics and Engineering*, 403:115717, jan 2023.

[VD17] Ehrlacher Virginie and Lombardi Damiano. A dynamical adaptive tensor method for the Vlasov-Poisson system. *Journal of Computational Physics*, 339:285–306, jun 2017.

[ZG18] Chenglong Zhang and Irene M Gamba. A conservative discontinuous galerkin solver for the space homogeneous boltzmann equation for binary interactions. *SIAM Journal on Numerical Analysis*, 56(5):3040–3070, 2018.

# A Numerical discretization details

## A.1 Discretization of the full model

Equation (7) is a hyperbolic system with a source term whose transport parts can be solved using a finite volume method, while the Poisson equation (8) is discretized with a finite difference method. The n-th iteration of the scheme reads:

$$-\frac{\phi_{i+1}^n - 2\phi_i^n - \phi_{i-1}^n}{\Delta x^2} = \frac{1}{L} \sum_{i=0}^{N_x-1} \Delta x \Delta v(1, \dots, 1)^T \mathbf{f}_i^n - \Delta v(1, \dots, 1)^T \mathbf{f}_i^n, \\ \frac{\mathbf{f}_i^{n+1} - \mathbf{f}_i^n}{\Delta t} + \frac{\mathbf{f}_{i+\frac{1}{2}}^n - \mathbf{f}_{i-\frac{1}{2}}^n}{\Delta x} + \left( \frac{\phi_{i+1}^n - \phi_{i-1}^n}{2\Delta x} \right) D \mathbf{f}_i^n = \frac{1}{\varepsilon} Q_{\Delta v}(\mathbf{f}_i^n)_i$$

where $\mathbf{f}_{i+\frac{1}{2}}^n$ is an upwind approximation with respect to the velocity given by:

$$\mathbf{f}_{i+\frac{1}{2}}^n = \text{Diag}(\mathbf{v}) \frac{\mathbf{f}_{i+1}^n + \mathbf{f}_i^n}{2} - \text{Diag}(|\mathbf{v}|) \frac{\mathbf{f}_{i+1}^n - \mathbf{f}_i^n}{2}.$$

To obtain better accuracy in space, we use a MUSCL method in the upwind flux. The transport scheme in space and velocity induces a CFL stability condition:

$$\Delta t \le \min \left( \frac{\Delta x}{v_{\max}}, \frac{\Delta v}{\max_i \left| \frac{\phi_{i+1}^n - \phi_{i-1}^n}{2\Delta x} \right|} \right).$$

The explicit treatment of the collision operator results in an additional constraint: $\Delta t \le \varepsilon C \Delta v^2$ with $C > 0$ related to the discrete Maxwellian distributions [BDS01]. In practice, we take $\Delta t$ sufficiently small in order to ensure the stability of the numerical simulations. Recently, high-order Runge Kutta methods have been proposed to relax the constraint [CA22].
[GT09] Irene M Gamba and Sri Harsha Tharkabhushanam. Spectral-lagrangian methods for collisional models of non-equilibrium statistical states. *Journal of Computational Physics*, 228(6):2012–2036, 2009.

[HPR22a] Jan S. Hesthaven, Cecilia Pagliantini, and Nicolò Ripamonti. Adaptive symplectic model order reduction of parametric particle-based Vlasov-Poisson equation, 2022.

[HPR22b] Jan S. Hesthaven, Cecilia Pagliantini, and Gianluigi Rozza. Reduced basis methods for time-dependent problems. *Acta Numerica*, 31:265–345, 2022.

[KK20] Lee Kookjin and Carlberg Kevin, T. Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders. *Journal of Computational Physics*, 404:108973, 2020.

[Rif20] Sébastien Riffaud. *Reduced-order models: convergence between scientific computing and data for fluid mechanics*. PhD thesis, Université de Bordeaux, 2020.

[RSK23] Geelen Rudy, Wright Stephen, and Willcox Karen. Operator inference for non-intrusive model reduction with quadratic manifolds. *Computer Methods in Applied Mechanics and Engineering*, 403:115717, jan 2023.

[VD17] Ehrlacher Virginie and Lombardi Damiano. A dynamical adaptive tensor method for the Vlasov-Poisson system. *Journal of Computational Physics*, 339:285–306, jun 2017.

[ZG18] Chenglong Zhang and Irene M Gamba. A conservative discontinuous galerkin solver for the space homogeneous boltzmann equation for binary interactions. *SIAM Journal on Numerical Analysis*, 56(5):3040–3070, 2018.

# A Numerical discretization details

## A.1 Discretization of the full model

Equation (7) is a hyperbolic system with a source term whose transport parts can be solved using a finite volume method, while the Poisson equation (8) is discretized with a finite difference method. The n-th iteration of the scheme reads:

$$-\frac{\phi_{i+1}^n - 2\phi_i^n - \phi_{i-1}^n}{\Delta x^2} = \frac{1}{L} \sum_{i=0}^{N_x-1} \Delta x \Delta v(1, \dots, 1)^T \mathbf{f}_i^n - \Delta v(1, \dots, 1)^T \mathbf{f}_i^n, \\ \frac{\mathbf{f}_i^{n+1} - \mathbf{f}_i^n}{\Delta t} + \frac{\mathbf{f}_{i+\frac{1}{2}}^n - \mathbf{f}_{i-\frac{1}{2}}^n}{\Delta x} + \left( \frac{\phi_{i+1}^n - \phi_{i-1}^n}{2\Delta x} \right) D \mathbf{f}_i^n = \frac{1}{\varepsilon} Q_{\Delta v}(\mathbf{f}_i^n)_i$$

where $\mathbf{f}_{i+\frac{1}{2}}^n$ is an upwind approximation with respect to the velocity given by:

$$\mathbf{f}_{i+\frac{1}{2}}^n = \text{Diag}(\mathbf{v}) \frac{\mathbf{f}_{i+1}^n + \mathbf{f}_i^n}{2} - \text{Diag}(|\mathbf{v}|) \frac{\mathbf{f}_{i+1}^n - \mathbf{f}_i^n}{2}.$$

To obtain better accuracy in space, we use a MUSCL method in the upwind flux. The transport scheme in space and velocity induces a CFL stability condition:

$$\Delta t \le \min \left( \frac{\Delta x}{v_{\max}}, \frac{\Delta v}{\max_i \left| \frac{\phi_{i+1}^n - \phi_{i-1}^n}{2\Delta x} \right|} \right).$$

The explicit treatment of the collision operator results in an additional constraint: $\Delta t \le \varepsilon C \Delta v^2$ with $C > 0$ related to the discrete Maxwellian distributions [BDS01]. In practice, we take $\Delta t$ sufficiently small in order to ensure the stability of the numerical simulations. Recently, high-order Runge Kutta methods have been proposed to relax the constraint [CA22].
## A.2 Discretization of the reduced model

Since the reduced model is a hyperbolic system, we use a finite volume scheme. Since the transport is linear, we consider the upwind scheme, whose fluxes write:

$$ \hat{f}_{i+\frac{1}{2}}^n = \hat{A} \frac{\hat{f}_{i+1}^n + \hat{f}_i^n}{2} - |\hat{A}| \frac{\hat{f}_{i+1}^n - \hat{f}_i^n}{2}. \quad (28) $$

As for the kinetic solver, we use a MUSCL strategy to obtain second order accuracy.

# B Solution to the minimization problem

**Proposition 2.** Let $\Phi \in M_{K,N_v}$, with $\Phi^T\Phi = \text{Id}$, $\hat{m} \in M_{3,K}$ of full rank. Then the unique solution to the least-square problem:

$$ \hat{N}^c = \underset{\hat{N} \in \mathbb{R}^K \text{ s.t. } \hat{m}\hat{N} = b}{\operatorname{argmin}} \|\Phi\hat{N} - N\|^2, $$

is given by:

$$ \hat{N}^c = \Phi^T N + \hat{m}^T (\hat{m}\hat{m}^T)^{-1} (b - \hat{m}\Phi^T N). $$

*Proof.* $\Phi$ being of full rank, the least-square problem under constraints has a unique solution denoted $\hat{N}^c$. Let $f(\hat{N}) = \|\Phi\hat{N} - N\|^2$. Its gradient is given by:

$$ \nabla f(\hat{N}) = 2(\Phi^T \Phi \hat{N} - \Phi^T N) = 2(\hat{N} - \Phi^T N). $$

Applying the Lagrange multiplier method, there exists $\lambda \in \mathbb{R}^K$ such that

$$ 2(\hat{N}^c - \Phi^T N) + \hat{m}^T \lambda = 0, \quad (29) $$

$$ \hat{m}\hat{N}^c = b. \quad (30) $$

Then multiplying (29) by $\hat{m}$ and using (30), we get the value of $\lambda$:

$$ \lambda = -2(\hat{m}\hat{m}^T)^{-1}(b - \hat{m}\Phi^T N). $$

Then, reporting this value into (29), we obtain the expected solution. $\square$
## A.2 Discretization of the reduced model

Since the reduced model is a hyperbolic system, we use a finite volume scheme. Since the transport is linear, we consider the upwind scheme, whose fluxes write:

$$ \hat{\mathbf{f}}_{i+\frac{1}{2}}^{n} = \hat{A} \frac{\hat{\mathbf{f}}_{i+1}^{n} + \hat{\mathbf{f}}_{i}^{n}}{2} - |\hat{A}| \frac{\hat{\mathbf{f}}_{i+1}^{n} - \hat{\mathbf{f}}_{i}^{n}}{2}. \quad (28) $$

As for the kinetic solver, we use a MUSCL strategy to obtain second order accuracy.

# B Solution to the minimization problem

**Proposition 2.** Let $\Phi \in M_{K,N_v}$, with $\Phi^T\Phi = \text{Id}$, $\hat{m} \in M_{3,K}$ of full rank. Then the unique solution to the least-square problem:

$$ \hat{N}^c = \underset{\hat{N} \in \mathbb{R}^K \text{ s.t. } \hat{m}\hat{N} = b}{\operatorname{argmin}} \|\Phi\hat{N} - N\|^2, $$

is given by:

$$ \hat{N}^c = \Phi^T N + \hat{m}^T (\hat{m}\hat{m}^T)^{-1} (b - \hat{m}\Phi^T N). $$

*Proof.* $\Phi$ being of full rank, the least-square problem under constraints has a unique solution denoted $\hat{N}^c$. Let $f(\hat{N}) = \|\Phi\hat{N} - N\|^2$. Its gradient is given by:

$$ \nabla f(\hat{N}) = 2(\Phi^T \Phi \hat{N} - \Phi^T N) = 2(\hat{N} - \Phi^T N). $$

Applying the Lagrange multiplier method, there exists $\lambda \in \mathbb{R}^K$ such that

$$ 2(\hat{N}^c - \Phi^T N) + \hat{m}^T \lambda = 0, \quad (29) $$

$$ \hat{m}\hat{N}^c = b. \quad (30) $$

Then multiplying (29) by $\hat{m}$ and using (30), we get the value of $\lambda$:

$$ \lambda = -2(\hat{m}\hat{m}^T)^{-1}(b - \hat{m}\Phi^T N). $$

Then, reporting this value into (29), we obtain the expected solution. □