# High order Asymptotic-Preserving penalized numerical schemes for the Euler-Poisson system in the quasineutral limit

Nicolas Crouseilles, Giacomo Dimarco, Saurav Samantaray

## ► To cite this version:

Nicolas Crouseilles, Giacomo Dimarco, Saurav Samantaray. High order Asymptotic-Preserving penalized numerical schemes for the Euler-Poisson system in the quasineutral limit. 2024. hal-04844514

HAL Id: hal-04844514

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

Preprint submitted on 17 Dec 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# High order Asymptotic-Preserving penalized numerical schemes for the Euler-Poisson system in the quasineutral limit

Nicolas Crouseilles¹, Giacomo Dimarco², and Saurav Samantaray³

¹Université de Rennes, Inria Rennes (Mingus team) and IRMAR UMR CNRS 6625, F-35042 Rennes, France & ENS Rennes

²Department of Mathematics and Computer Science & Center for Computing, Modeling and Statistics (CMCS), University of Ferrara. Italy

³Department of Mathematics, Indian Institute Of Technology Madras, Chennai, India

## Abstract

In this work, we focus on developing a new class of numerical methods able to both handle quasineutrality and charge separation in plasmas. At large temporal and spatial scales, plasmas tend to be quasineutral, meaning that the local net charge density is nearly zero. However, when the scale at which one observes the plasma dynamics is smaller than the characteristic distance over which the electric field and charges are typically screened, then quasineutrality breaks down. In such regimes, standard numerical methods face severe stability constraints, rendering them practically unusable. To address this issue, in this work, we introduce and analyze a new class of finite volume penalized-IMEX Runge-Kutta methods for the Euler-Poisson (EP) system, specifically designed to handle the quasineutral limit. We show that, these proposed schemes are uniformly stable with respect to the Debye length and degenerate into high order methods as the quasineutral limit is approached. Several numerical tests confirm that this new class of methods exhibits the desired properties.

**Keywords**— Plasma simulation, quasineutral limit, Euler-Poisson system, Asymptotic-Preserving schemes, Debye length

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>One-fluid Euler-Poisson system and its quasineutral limit</td><td>3</td></tr><tr><td>3</td><td>Time discretization</td><td>4</td></tr><tr><td>3.1</td><td>A possible first order time discretization</td><td>5</td></tr><tr><td>3.2</td><td>An alternative first order time discretization for the Euler-Poisson system</td><td>7</td></tr><tr><td>3.3</td><td>An IMEX Runge-Kutta class of methods for handling quasineutrality in plasma</td><td>9</td></tr><tr><td>3.3.1</td><td>Generalities on IMEX Runge-Kutta methods</td><td>9</td></tr><tr><td>3.3.2</td><td>High order time semi-discretization for the Euler-Poisson system</td><td>10</td></tr><tr><td>4</td><td>Analysis of the penalized-IMEX schemes</td><td>11</td></tr><tr><td>4.1</td><td>Analysis of type-A IMEX Runge-Kutta schemes</td><td>12</td></tr><tr><td>4.2</td><td>Analysis of type-CK IMEX Runge-Kutta schemes</td><td>15</td></tr><tr><td>5</td><td>Space discretization</td><td>15</td></tr><tr><td>6</td><td>Numerical results</td><td>16</td></tr><tr><td>6.1</td><td>Asymptotic order of convergence</td><td>16</td></tr><tr><td>6.2</td><td>Perturbation of a quasineutral plasma in 1D with ρ = 1, u = 1</td><td>18</td></tr><tr><td>6.3</td><td>Perturbation of a quasineutral plasma in 1D with ρ = 1, u = 0</td><td>19</td></tr><tr><td>6.4</td><td>Perturbation of a quasineutral plasma in 2D with ρ = 1, u = (1, 1)</td><td>21</td></tr><tr><td>7</td><td>Conclusions</td><td>22</td></tr></table>
# 1 Introduction

The modeling and numerical simulation of plasma phenomena is a highly active field of research, driven by its significance in both applied technologies and fundamental scientific exploration [10]. Plasmas, often referred to as the fourth state of matter, are basically ionized gases that exhibit unique behaviors under the influence of electromagnetic fields. Their study is crucial for a broad range of technological applications, such as in the development of fusion energy, where confined plasmas hold the potential for generating nearly limitless clean energy [23, 25]. Plasmas are also essential in industries like semiconductor manufacturing [43], space propulsion systems [37], medical devices, and environmental control technologies [26], where ionized gases are used for applications such as waste treatment and air purification.

In this paper, we tackle the specific problem of developing suitable numerical methods for the description of plasmas [17, 18], focusing on the challenges associated with the coexistence of quasineutrality and neutrality breakdown. Our study concentrates on a macroscopic description, namely the one-fluid Euler-Poisson (EP) system, which models the dynamics of electrons in a plasma with a constant ion background. In this system, the density and momentum evolve according to a hyperbolic structure influenced by an electrostatic potential that solves an elliptic equation. Despite of its apparent simplicity, this model has been chosen since it contains most of the challenges related to the treatment of the quasineutrality in plasmas. Future investigations will be about the extension of the study here presented to more sophisticated descriptions of plasma phenomena.

An important physical parameter, which plays a fundamental role in this context is the so-called Debye length, $\lambda$, which measures the typical length scale of charge imbalances in a plasma. Indeed in quasineutral regions, the electric charges macroscopically effectively cancel out while plasma still carries a non zero electric current. To take into account this fact and to be able to describe such phenomena, a numerical scheme has to be able to capture the small scales corresponding to the Debye length, typically much smaller than the overall size of the system. The direct consequence is that, in such cases, using explicit time discretization schemes leads to prohibitively expensive simulations [17]. Furthermore, in many practical situations, such as sheath formation, plasma diode modeling, or arc formation on satellite solar cells, the coexistence of both quasineutral and non-quasineutral regions is present. These scenarios require consequently specialized numerical treatments and the design of new suitable numerical schemes. One promising approach in this direction involves the use of domain decomposition techniques, which have proven highly effective in managing multiple regimes in plasmas, see for instance [21, 12]. These methods offer a viable alternative to standard explicit numerical schemes that would otherwise necessitate resolving the smallest scale (the Debye length), leading to extremely small time steps and too high computational costs. However, to the best of our knowledge, the use of domain decomposition methods to connect the asymptotic quasineutral model with the Euler-Poisson system has not yet been explored. This connection presents unique challenges, as developing appropriate numerical methods and accurately identifying the interface between different models requires specific attention. Indeed, the task of ensuring a smooth transition between the models is complex and demands further research and developments.

Several alternative approaches have been proposed in the literature to address the challenges posed by the coexistence of different regimes in plasmas. Among the most successful in recent years, there are the so-called Asymptotic-Preserving (AP) schemes [32, 33]. These schemes enjoy the property of being uniformly stable with respect to small scale parameters. Moreover, when an asymptotic model can be identified from these small scale dynamics, they permit a smooth transition into a consistent discretization for this limit model. In other words, the limit of such AP schemes is both a consistent and accurate discretization for the corresponding limiting model. In the context of the quasineutral scaling for the Euler-Poisson system, notable works include the researches [13, 15, 17, 3, 1], as well as [24, 19], which provide comprehensive numerical analyses of the problem. Additional relevant studies on numerical methods for handling quasineutrality can be found in [6, 5, 14, 18, 16, 35]. More generally, Asymptotic-Preserving schemes have garnered significant interest in recent decades for their ability to approximate multiscale partial differential equations (PDEs), which are ubiquitous in plasma physics. For further exploration of this topic, see [32, 33, 38, 17]. These AP schemes are often based on Implicit-Explicit (IMEX) time integrators, which have a long history of development and have benefited from numerous contributions in the past. Important references on IMEX methods include [41, 29, 42, 34, 39, 7, 22].

The primary goal of this paper is to design a new class of high order schemes for the Euler-Poisson system that are uniformly stable with respect to the Debye length and degenerate into consistent and accurate discretizations of the quasineutral system in the limit when the Debye length approaches zero. Although significant progress has been made in developing Asymptotic-Preserving (AP) schemes in the quasineutral limit [13, 15, 17, 3, 1, 36, 31], few contributions have focused on the development of high order AP schemes which is the direction pursued in this work. Another key contribution of this research is the ability of the proposed schemes to handle very general initial conditions. In fact, most of the existing studies on the quasineutral limit, assume that the initial conditions are well-prepared, meaning they are compatible with the asymptotic model. However, this assumption can be quite restrictive in practice, as one would expect the numerical scheme to be capable of handling a broad range of initial conditions. For this reason, in this paper, we explore the challenges posed by the choice of general initial conditions, the implications they have for the numerical solution, and we propose a strategy that ensures stability and consistency with the limit system, even when arbitrarily chosen initial conditions are used. In particular, we
will show that in the cases in which the construction of the numerical method does not take into account the possibility of handling rather arbitrarily initial conditions, the solution may grow and explode in finite time.

The strategy employed to ensure the above discussed properties is primarily based on an Implicit-Explicit (IMEX) time integration approach combined with a suitable penalization technique. Following the ideas from [1] and several works that focused on the low Mach limit in fluid dynamics [4, 8, 9, 20], an implicit scheme for density $\rho$ can be derived by utilizing the divergence of the momentum equation and thanks to the introduction of a suitable penalization term containing the asymptotic model solution. This approach allows for the projection of the density onto the quasineutral state, enabling the Poisson equation to be solved in a stable manner. Following the above path, in this work, we first introduce such a penalization technique and then extend a particular first order time discretization through the use of IMEX Runge-Kutta methods. We successively perform a numerical analysis of the proposed class of IMEX schemes [22, 30, 2] and we prove that type-A IMEX Runge-Kutta methods are able to ensure the Asymptotic-Preserving (AP) property for a broad class of initial conditions when the penalization is introduced. This stands in contrast to other classes of methods, such as the type-CK IMEX Runge-Kutta schemes, which fail to handle initial data that are not consistent with the limit system. In particular, we will show that this latter class of methods exhibits instabilities in certain situations.

The rest of the paper is organized as follows. In the next section, the Euler-Poisson model and its formal quasineutral limit are recalled. Then, in Section 3, two different first order time discretizations are introduced and discussed and then a suitable extension to high order IMEX Runge-Kutta schemes is derived. In Section 4, a numerical analysis of the proposed schemes is performed showing, in particular, consistency of the methods with the asymptotic model for arbitrarily initial data. In the same section, we illustrate the limits of an alternative class of IMEX schemes by showing that they may not be able to provide the required Asymptotic-Preserving property and stability. Section 5 details the space approximation, and finally, in Section 6, some numerical results are proposed to illustrate the capability of this new class of schemes compared to other approaches from the literature. A final section permits us to draw some conclusions and suggest future investigations.

## 2 One-fluid Euler-Poisson system and its quasineutral limit

We model a plasma by considering electrons to be always in thermodynamic equilibrium while heavier ions to be in a state of a fixed background bath. Under the previous hypothesis, we have a macroscopic fluid of negative charged particles coupled with an ensemble of steady atoms with a positive net electric charge. Finally, if the fluid is also considered isentropic then this system of particles is well described by the so-called Euler-Poisson (EP) equations [17]

$$ 
\begin{aligned}
\partial_t \rho + \nabla \cdot q &= 0, \\
m \partial_t q + m \nabla \cdot (\rho u \otimes u + p I) &= e \rho \nabla \phi, \\
\Delta \phi &= \frac{e}{\varepsilon_0} (\rho - \rho_0), \quad x \in \mathbb{T}.
\end{aligned}
 $$

where $\rho(t, x)$ is the density of electrons, $\rho_0$ the neutralizing ion background, $q(t, x) := (\rho u)(t, x)$ their momentum ($u$ being the mean fluid velocity), $m$ is the electron mass, $\varepsilon_0$ the permittivity, $e$ the electron charge while the pressure obeys to the law $p(\rho) = C\rho^\beta$ with $\beta \ge 1$ with $I$ denoting the identity matrix. Finally, the spatial differential operators $\nabla$ and $\Delta$ used are the gradient and the Laplacian, respectively, while $\mathbb{T}$ is the domain. In the non dimensional form the same Euler-Poisson system reads as:

$$ \partial_t \rho + \nabla \cdot q = 0, \qquad (1) $$

$$ \partial_t q + \nabla \cdot (\rho u \otimes u + p(\rho)I) = \rho \nabla \phi, \qquad (2) $$

$$ \lambda^2 \Delta \phi = \rho - 1, \quad x \in \mathbb{T}, \qquad (3) $$

with now $p(\rho) = \rho^\beta$. Moreover, the electric field $E$ is obtained as the negative gradient of the electric potential $\phi(t, x)$, i.e. $E = -\nabla\phi$ under the hypothesis of zero magnetic field. The system is supplemented with initial conditions $(\rho, q)(t = 0, x)$, while periodic boundary conditions are considered in space if not otherwise stated. The uniform ion background is, in the non dimensional case, represented by the constant density 1 in the Poisson equation (3). In (3) it appears the scaled Debye length

$$ \lambda = \left( \frac{\varepsilon_0 k_B T_0}{x_0^2 e^2 \rho_0} \right)^{1/2} $$

with $x_0$ the typical length of the problem, $\rho_0$ and $T_0$ the typical density and temperature of the electrons, $k_B$ the Boltzmann constant. Note also that in the case of one-species model, the scaled plasma period equals the Debye length:

$$ \tau = \frac{\tau_p}{x_0/v_0} \left( \frac{\varepsilon_0 m}{e^2 \rho_0} \right)^{1/2} = \frac{1}{x_0} \left( \frac{\varepsilon_0 m v_0^2}{e^2 \rho_0} \right)^{1/2} = \lambda, \quad \text{since} \quad v_0 = (K_B T_0 / m)^{1/2}, $$
meaning that the electrons oscillate around the steady ions at a speed proportional to $1/\lambda$.

Let now observe that the Euler-Poisson system (1)-(2)-(3) becomes singular in the limit of $\lambda \to 0$, see [3, 17] for details. This is the so-called quasineutral limit in which macroscopically charge imbalances disappear. This can be formally observed by taking $\lambda \to 0$ in the Poisson equation (3) obtaining $\rho = 1$, i.e. a singular limit. In other words, the electrostatic potential becomes the Lagrange multiplier of the constraint $\rho = 1$ as in the incompressible Euler equations in which the pressure is the Lagrange multiplier of the divergence free constraint. In this limit, it is not possible to recover directly an explicit expression for the electric field, while it is well known that inside a plasma, even in the quasineutral regime, currents are present.

We then discuss a formal derivation of the asymptotic limit system obtained when $\lambda \to 0$ from the EP system (1)-(3) giving rise to an explicit equation for the potential $\phi$ and consequently to a non-zero electric field. We refer the readers to [13, 17] for more details, and to [40, 11, 27] for theoretical aspects and a rigorous derivation associated to the quasineutral limit. Formally the variables $\rho, q$ and $\phi$ yield the following limits as $\lambda \to 0$:

$$ \lim_{\lambda \to 0} \rho = \rho_{(0)} = 1, \quad \lim_{\lambda \to 0} q = q_{(0)} = \rho_{(0)} u_{(0)}, \quad \lim_{\lambda \to 0} \phi = \phi_{(0)}. $$

Hence, from (1)-(3) we obtain the incompressible Euler-Poisson system:

$$ \nabla \cdot u_{(0)} = 0, \qquad (4) $$

$$ \partial_t u_{(0)} + \nabla \cdot (u_{(0)} \otimes u_{(0)}) = \nabla \phi_{(0)}, \qquad (5) $$

$$ \rho_{(0)} = 1, \qquad (6) $$

where one can notice that the negative electric potential $-\phi_{(0)}$ plays the role of the incompressible pressure. Let us now remark that the divergence free constraint $\nabla \cdot u_{(0)} = 0$ can be reformulated. Indeed, taking the divergence of (5) and using (4) give

$$ \nabla^2 : (u_{(0)} \otimes u_{(0)}) = \Delta \phi_{(0)}, $$

so that the limit system (4)-(5) can be reformulated as

$$ \partial_t u_{(0)} + \nabla \cdot (u_{(0)} \otimes u_{(0)}) = \nabla \phi_{(0)}, \qquad (7) $$

$$ \Delta \phi_{(0)} = \nabla^2 : (u_{(0)} \otimes u_{(0)}), \qquad (8) $$

where $\nabla^2$ denotes the tensor of second order derivatives and $:$ the contracted product of two tensors. Note that the initial data has to be well-prepared in the following sense for the solution of the EP system (1)-(3) to converge to the solution of the quasineutral system (4)-(6).

**Definition 1.** A triple $(\rho, q, \phi)$ of solution of the EP system (1)-(3) is called well-prepared if it admits the following form:

$$ \begin{aligned} \rho &= \rho_{(0)} + \lambda^2 \rho_{(2)}, & q &= q_{(0)} + \lambda^2 q_{(2)} \\ \rho_{(0)} &= 1, & \nabla \cdot q_{(0)} &= 0. \end{aligned} \qquad (9) $$

The primary objective of this work is to develop appropriate time discretizations for (1)-(3) that are uniformly stable with respect to the Debye length $\lambda$ and degenerate, as $\lambda \to 0$, into a consistent discretization of the asymptotic quasineutral model (4)-(6). The numerical schemes proposed here are then designed to handle different regimes that may arise in plasma simulations, without imposing time step restrictions due to fast-scale dynamics. In other words, the main goal consists in having a tool able to approximate complex plasma phenomena independently on the fast and/or slow scales which may appear during the simulations. The above discussed property is often called in the literature Asymptotic-Preserving, where a formal well-established definition is

**Definition 2** (Strong Asymptotic-Preserving property). A consistent time discretization method for (1)-(3) of stepsize $\Delta t$ is Asymptotic-Preserving (AP) if, independently of the initial data and of the stepsize $\Delta t$, in the limit $\lambda \to 0$ it becomes a consistent time discretization method for the asymptotic system (4)-(6).

As opposite, we refer to as the weak AP property when the above definition is fulfilled only when the initial data are well-prepared as stated in Definition 1. Finally, let observe that these definitions of strong and weak AP schemes do not imply that the schemes preserve the order of accuracy in time in the stiff limit $\lambda \to 0$. In the latter case we will say that the schemes are also asymptotically accurate (AA).

## 3 Time discretization

We start our analysis of possible time discretizations for the Euler-Poisson model by presenting two different first order schemes. This will permit to highlight the main challenges and difficulties related to the quasineutrality. Indeed, a first order discretization permits to have an easy access into the design and analysis elements of the semi-discrete approach developed and to suggest the best strategy to pursue as it will be detailed later. In a second part of the section, we propose a new class of high order extension of the first order methods presented.

The main components of the time semi-discrete schemes are summarized in the following three steps:
* First, a suitable penalization for the electric potential is introduced into the momentum equation.

* Second, two discretizations of a reformulated Euler-Poisson system are proposed where some implicit terms appear.

* Third, equivalent reformulated time semi-discrete schemes are derived from the previous step which are shown to be fully explicit.

Let us now consider the first step of the above procedure. We reformulate the EP system (1)-(3) forming the basis for the time discretizations to follow. The penalization procedure consists in introducing in the momentum equation (2) the quasineutral electric potential solution, i.e.

$$ \partial_t q + \nabla \cdot (\rho u \otimes u + p(\rho)I) = (\rho - 1)\nabla\phi + \nabla\phi, \quad (10) $$

obtained by simply adding and subtracting $\nabla\phi$. With this formulation in hand, we introduce the following notations

$$ \begin{aligned} U := (\rho, q)^T, \quad F := (q \otimes q)/\rho + p(\rho)I, \quad G_{\text{ex}} := (0, \nabla \cdot F - (\rho - 1)\nabla\phi)^T, \quad G_{\text{im}} := (\nabla \cdot q, -\nabla\phi)^T \\ H(U, \phi) := -\lambda^2 \Delta\phi - \rho + 1, \end{aligned} \quad (11) $$

so that the Euler-Poisson system (1)-(3) becomes

$$ \begin{aligned} \partial_t U + G_{\text{ex}} + G_{\text{im}} &= 0, \\ H(U, \phi) &= 0. \end{aligned} \quad (12) $$

Based on this reformulation, two different Implicit-Explicit time discretizations [8, 16] are proposed by taking in both cases the mass flux and the gradient of the electric potential $\nabla\phi$ implicitly while the rest of the terms appearing the Euler-Poisson system explicitly. In the next paragraph, we discuss and analyse the schemes obtained based on the above choice, highlighting the main differences between the two.

## 3.1 A possible first order time discretization

The most natural first order time discretization of the Euler-Poisson system (1)-(3) through the decomposition (12) reads:

$$ \begin{aligned} \frac{U^{n+1} - U^n}{\Delta t} &= -G_{\text{ex}}^n - G_{\text{im}}^{n+1}, \\ \lambda^2 \Delta \phi^{n+1} &= \rho^{n+1} - 1, \end{aligned} \quad (13) $$

which, using the notations $U = (\rho, q)$ and $G_{\text{ex}}, G_{\text{im}}$ from (11), corresponds to

$$ \rho^{n+1} = \rho^n - \Delta t \nabla \cdot q^{n+1}, \quad (14) $$

$$ \begin{aligned} q^{n+1} &= q^n - \Delta t \nabla \cdot (\rho^n u^n \otimes u^n + p^n I) + \Delta t \rho^n \nabla \phi^n - \Delta t \nabla \phi^n + \Delta t \nabla \phi^{n+1}, \\ &= q^n - \Delta t \nabla \cdot F^n + \Delta t (\rho^n - 1) \nabla \phi^n + \Delta t \nabla \phi^{n+1}, \end{aligned} \quad (15) $$

$$ \lambda^2 \Delta \phi^{n+1} = \rho^{n+1} - 1. \quad (16) $$

Note the above scheme is semi-implicit in nature and both, the mass and momentum updates have implicit terms. In order to simplify the scheme with respect to the implicit quantities, we insert the divergence of the momentum equation (15) into the mass equation (14). This leads to:

$$ \rho^{n+1} = \rho^n - \Delta t \nabla \cdot q^{n+1} = \rho^n - \Delta t \nabla \cdot q^n + \Delta t^2 \nabla^2 : F^n + \Delta t^2 \nabla \cdot ((\rho^n - 1) \nabla \phi^n) - \Delta t^2 \Delta \phi^{n+1}. $$

Now, using the Poisson equation $\lambda^2 \Delta\phi^{n+1} = \rho^{n+1} - 1$ and subtracting 1 on both sides give:

$$ \rho^{n+1} - 1 = \rho^n - 1 - \Delta t \nabla \cdot q^{n+1} = \rho^n - 1 - \Delta t \nabla \cdot q^n + \Delta t^2 \nabla^2 : F^n + \Delta t^2 \nabla \cdot ((\rho^n - 1) \nabla \phi^n) - \frac{\Delta t^2}{\lambda^2} (\rho^{n+1} - 1), $$

from which we deduce $\rho^{n+1}$

$$ \begin{aligned} \rho^{n+1} &= 1 + \frac{\lambda^2}{\lambda^2 + \Delta t^2} [\rho^n - 1 - \Delta t \nabla \cdot q^n + \Delta t^2 \nabla^2 : F^n + \Delta t^2 \nabla \cdot ((\rho^n - 1) \nabla \phi^n)] \\ &= 1 + \frac{\lambda^2}{\lambda^2 + \Delta t^2} [\rho^n - 1 - \Delta t \nabla \cdot q^n + \Delta t^2 \nabla^2 : F^n + \Delta t^2 \nabla \rho^n \cdot \nabla \phi^n + \Delta t^2 (\rho^n - 1) \Delta \phi^n] \\ &= 1 + \frac{\lambda^2}{\lambda^2 + \Delta t^2} [\rho^n - 1 - \Delta t \nabla \cdot q^n + \Delta t^2 \nabla^2 : F^n + \Delta t^2 \nabla \rho^n \cdot \nabla \phi^n] + \frac{\Delta t^2}{\lambda^2 + \Delta t^2} (\rho^n - 1)^2. \end{aligned} \quad (17) $$
Now, under the assumption $\rho(t=0) = 1 + O(\lambda^2)$ on the initial conditions, $\rho^{n+1}$ is updated from (17) so that we can prove $\rho^{n+1} = 1 + O(\lambda^2)$. This particular choice on the initial data is referred to as well-prepared initial conditions following the Definition 1. So under the above condition, the Poisson equation can be now solved:

$$ 
\begin{aligned}
\Delta\phi^{n+1} &= \frac{\rho^{n+1} - 1}{\lambda^2} \\
&= \frac{1}{\lambda^2 + \Delta t^2} \left[ \rho^n - 1 - \Delta t \nabla \cdot q^n + \Delta t^2 \nabla^2 : F^n + \Delta t^2 \nabla \rho^n \cdot \nabla \phi^n \right] + \frac{\Delta t^2}{\lambda^2 + \Delta t^2} \frac{(\rho^n - 1)^2}{\lambda^2},
\end{aligned}
\quad (18) $$

and finally the momentum equation can be updated explicitly

$$ q^{n+1} = q^n - \Delta t \nabla \cdot F^n + \Delta t(\rho^n - 1)\nabla\phi^n + \Delta t \nabla\phi^{n+1}. \quad (19) $$

The above proposed first order in time discretization possesses the following property:

**Theorem 1.** Let us assume that the initial conditions

$$ \rho^0 = \rho(t=0), \quad q^0 = q(t=0), \quad \text{and} \quad \phi^0 = \phi(t=0), $$

are well-prepared in the sense of Definition 1 and $\lambda^2\phi^0 = \rho^0 - 1$. Then, the time semi-discrete scheme (17)-(19) initialised with the above initial conditions enjoys the following properties when $\lambda \to 0$

$$ \rho^n = 1 + O(\lambda^2), \quad \nabla \cdot u^n = O(\lambda^2), \quad \phi^n = O(1), \quad n \ge 1. $$

Moreover, under the same hypothesis on the initial conditions, the scheme (17)-(19) degenerates when $\lambda \to 0$ towards

$$ u^{n+1} = u^n - \Delta t \nabla \cdot (u^n \otimes u^n) + \Delta t \nabla \phi^{n+1}, \quad \Delta \phi^{n+1} = \nabla^2 : (u^n \otimes u^n), \quad \text{with} \quad \rho^{n+1} = 1, \nabla \cdot u^{n+1} = 0 $$

which is a first order discretization of the asymptotic limit model (4)-(6).

*Proof.* Under the well-prepared initial conditions assumption, it is clear from (17) that $\rho^{n+1} = 1 + O(\lambda^2)$. To prove $\nabla \cdot u^{n+1} = O(\lambda^2)$, let us take the divergence of (19) to get

$$ \nabla \cdot q^{n+1} = \nabla \cdot q^n - \Delta t \nabla^2 : F^n + \Delta t \nabla \cdot \left[ (\rho^n - 1) \nabla \phi^n \right] + \Delta t \Delta \phi^{n+1}. $$

Using well-prepared initial conditions and the expression for $F$ from (11) leads to

$$ \nabla \cdot q^{n+1} = O(\lambda^2) - \Delta t \nabla^2 : (u^n \otimes u^n) + \Delta t \Delta \phi^{n+1}. $$

Now, let us consider (18), which degenerates when $\lambda \to 0$ (still under well-prepared initial conditions) into

$$ \Delta \phi^{n+1} = \nabla^2 : (u^n \otimes u^n) + O(\lambda^2), $$

so that, substituting (3.1) into (20) leads to $\nabla \cdot q^{n+1} = \nabla \cdot u^{n+1} = O(\lambda^2)$. Finally, going back to (19) and using the previous limits gives

$$ u^{n+1} = u^n - \Delta t \nabla \cdot (u^n \otimes u^n) + \Delta t \nabla \phi^{n+1}. $$

□

In conclusion, under the assumption of well-prepared initial conditions, the first order semi-discrete time scheme proposed in (17)-(19) reduces to a scheme that is consistent with the asymptotic model (4)-(6). Conversely, if the assumption of well-prepared initial conditions is not satisfied, it becomes evident from (18) that the term $\Delta\phi^n$ exhibits stiffness. Specifically, $\Delta\phi^n$ typically scales as $O(1/\lambda^2)$ when $\rho^n$ is not well-prepared. This significantly complicates the computation of the solution for small but positive values of $\lambda$, and makes it infeasible when $\lambda$ approaches zero. In the following sections, we will first explore a different first order discretization and then we discuss possible extensions to high order schemes. In particular, we will demonstrate that a suitable combination of implicit and explicit discretization together with the penalization (10) can overcome the limitations imposed by the well-preparedness assumption. A key result of this analysis is the identification of type-A IMEX Runge-Kutta schemes [42, 39] as an effective approach together with a suitable penalization to stabilize the Euler-Poisson system in the quasineutral limit, even when the initial data are not well-prepared.

**Remark 1.** The method discussed above relies on the Butcher tableau of Figure 1. One can observe that this method belong to a particular class of Implicit-Explicit Runge-Kutta methods for which the matrix of the implicit tableau is singular. This choice, as it will be discussed later, will play a fundamental role for the possibility of handling different kind of initial conditions for the Euler-Poisson system.
<table><caption>Figure 1: Double Butcher tableau for the first order Explicit-Implicit IMEX-RK scheme (14)-(16).</caption><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>1</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td></tr><tr><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td></tr></table>

**Remark 2.** Let notice that some similar approaches have been proposed in the recent past [3, 13, 15]. However, the order in which the unknowns are updated in (14)-(16) is different from the strategy proposed in those works. Indeed, the scheme proposed in [3] and previous works relied on a reformulated Poisson equation where the electric potential was first computed so that the momentum could be updated explicitly. Finally, an explicit update for the density was suggested. This is opposite from what has been discussed in this paragraph: the density is first computed explicitly, then the potential by solving the usual Poisson equation and finally an explicit update for the momentum is performed. This difference causes the schemes here proposed to possess different properties as discussed in the numerical test section.

## 3.2 An alternative first order time discretization for the Euler-Poisson system

We discuss now a possible modification of the method discussed in the previous section. The new tableau is shown in Figure 2. We refer to this method as Implicit-Explicit Euler IMEX-RK scheme as opposite to the previous

<table><caption>Figure 2: Double Butcher tableau for the first order Implicit-Explicit IMEX-RK scheme.</caption><tr><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td></tr><tr><td>1</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td></tr><tr><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td></tr></table>

one classically named Explicit-Implicit Euler IMEX-RK scheme. This new method consists of two stages and it gives at the first stage, where the apexes indicate the stage values,

$$ \rho^{n,(1)} = \rho^n - \Delta t \nabla \cdot q^{n,(1)}, \quad (20) $$

$$ q^{n,(1)} = q^n + \Delta t \nabla \phi^{n,(1)}, \quad (21) $$

$$ \lambda^2 \Delta \phi^{n,(1)} = \rho^{n,(1)} - 1, \quad (22) $$

while at the second stage reads

$$ \rho^{n,(2)} = \rho^n - \Delta t \nabla \cdot q^{n,(2)}, \quad (23) $$

$$ q^{n,(2)} = q^n - \Delta t \nabla \cdot (\rho^{n,(1)} u^{n,(1)} \otimes u^{n,(1)} + p^{n,(1)} I) + \Delta t \rho^{n,(1)} \nabla \phi^{n,(1)} - \Delta t \nabla \phi^{n,(1)} + \Delta t \nabla \phi^{n,(2)}, \quad (24) $$

$$ \lambda^2 \Delta \phi^{n,(2)} = \rho^{n,(2)} - 1. \quad (25) $$

Finally the solution at time step $(n + 1)$ coincides with the solution provided by the second stage, namely:

$$ \rho^{n+1} = \rho^{n,(2)}, \quad q^{n+1} = q^{n,(2)}, \quad \phi^{n+1} = \phi^{n,(2)}. \quad (26) $$

We now proceed as in the previous section and we compute the density of the first stage by inserting the divergence of the momentum (21)

$$ \rho^{n,(1)} = \rho^n - \Delta t \nabla \cdot q^{n,(1)} = \rho^n - \Delta t \nabla \cdot q^n - \Delta t^2 \Delta \phi^{n,(1)}, $$

and by substituting the Poisson equation into the above expression, we get:

$$ \rho^{n,(1)} - 1 = \rho^n - 1 - \Delta t \nabla \cdot q^n - \frac{\Delta t^2}{\lambda^2} (\rho^{n,(1)} - 1) \Rightarrow \rho^{n,(1)} = 1 + \frac{\lambda^2}{\lambda^2 + \Delta t^2} (\rho^n - 1 - \Delta t \nabla \cdot q^n). $$

Let us observe that now, as opposite to the case treated in the previous paragraph, $\rho^{n,(1)} = 1 + O(\lambda^2)$ without any hypothesis on the initial datum for the density. Moreover, for what concerns the Poisson equation we have:

$$ \Delta \phi^{n,(1)} = \frac{\rho^{n,(1)} - 1}{\lambda^2} = \frac{1}{\lambda^2 + \Delta t^2} (\rho^n - 1 - \Delta t \nabla \cdot q^n), $$

while the momentum equation gives

$$ q^{n,(1)} = q^n - \Delta t \nabla \phi^{n,(1)}. $$
This concludes the computation of the first stage of the scheme where one can notice that the density after the first stage is distant from the constant value 1 by $O(\lambda^2)$. The second stage is very close to the Explicit-Implicit Euler method (14)-(16) discussed previously and gives for the density

$$ \rho^{n,(2)} - 1 = \rho^n - 1 - \Delta t \nabla \cdot q^{n,(2)} = \rho^n - 1 - \Delta t \nabla \cdot q^n + \Delta t^2 \nabla^2 : F^{n,(1)} + \Delta t^2 \nabla \cdot ((\rho^{n,(1)} - 1) \nabla \phi^{n,(1)}) - \frac{\Delta t^2}{\lambda^2} (\rho^{n,(2)} - 1), $$

and then

$$ \begin{aligned} \rho^{n,(2)} = & \ 1 + \frac{\lambda^2}{\lambda^2 + \Delta t^2} \left[ \rho^n - 1 - \Delta t \nabla \cdot q^n + \Delta t^2 \nabla^2 : F^{n,(1)} + \Delta t^2 \nabla \rho^{n,(1)} \cdot \nabla \phi^{n,(1)} \right] + \right. \\ & \ \left. + \frac{\Delta t^2 \lambda^2}{\lambda^2 + \Delta t^2} (\rho^{n,(1)} - 1) \Delta \phi^{n,(1)}. \end{aligned} \quad (27) $$

The Poisson equation instead becomes

$$ \begin{aligned} \Delta \phi^{n,(2)} = & \ \frac{\rho^{n,(2)} - 1}{\lambda^2} = \frac{1}{\lambda^2 + \Delta t^2} \left[ \rho^n - 1 - \Delta t \nabla \cdot q^n + \Delta t^2 \nabla^2 : F^{n,(1)} + \Delta t^2 \nabla \rho^{n,(1)} \cdot \nabla \phi^{n,(1)} \right] + \right. \\ & \ \left. + \frac{\Delta t^2 \lambda^2}{\lambda^2 + \Delta t^2} (\rho^{n,(1)} - 1) \Delta \phi^{n,(1)}, \end{aligned} \quad (28) $$

and the momentum equation can be updated explicitly

$$ q^{n,(2)} = q^n - \Delta t \nabla \cdot F^{n,(1)} + \Delta t (\rho^{n,(1)} - 1) \nabla \phi^{n,(1)} + \Delta t \nabla \phi^{n,(2)}. \quad (29) $$

Finally, the solutions at the next time step coincide with the solution produced by the second stage as stated in (26). We can now summarize the main property of this new time discretization in the following Theorem.

**Theorem 2.** For any initial data, the time semi-discrete scheme (20)-(26), enjoys the following properties when $\lambda \to 0$, after the first two time steps,

$$ \rho^n = 1 + \mathcal{O}(\lambda^2), \quad \nabla \cdot u^n = \mathcal{O}(\lambda^2), \quad \phi^n = \mathcal{O}(1), \quad n \ge 2. $$

Moreover, after two time steps, the scheme (20)-(26) degenerates when $\lambda \to 0$ towards

$$ u^{n+1} = u^n - \Delta t \nabla \cdot (u^n \otimes u^n) + \Delta t \nabla \phi^{n+1}, \quad \Delta \phi^{n+1} = \nabla^2 : (u^n \otimes u^n), \quad \text{with } \rho^{n+1} = 1, \nabla \cdot u^{n+1} = 0 $$

which is a first order discretization of the asymptotic limit model (4)-(6).

*Proof.* Let first of all observe that the property enounced in the above Theorem is not common in the context of semi-discrete methods. The scheme indeed does not project the solution into a consistent discretization of the quasineutral limit in one single time step but it needs two consecutive time steps to perform such projection. In particular, this is true for the constraint on the divergence of the velocity while the density is projected towards the quasineutral state after the first time step. To see that this property holds true, let us notice that from (27) one obtains $\rho^{n+1} = \rho^{n,(2)} = 1 + \mathcal{O}(\lambda^2)$, then by taking the divergence of (29) one gets

$$ \nabla \cdot q^{n,(2)} = \nabla \cdot q^n - \Delta t \nabla^2 : F^{n,(1)} + \Delta t \nabla \cdot ((\rho^{n,(1)} - 1) \nabla \phi^{n,(1)}) + \Delta t \Delta \phi^{n,(2)}. \quad (30) $$

From (28), one can see that when $\lambda \to 0$ the electric potential degenerates to

$$ \Delta \phi^{(2)} = \nabla^2 : (u^{(1)} \otimes u^{(1)}) + \frac{1}{\Delta t^2} (\rho^n - 1 - \Delta t \nabla \cdot q^n) + \mathcal{O}(\lambda^2), \quad (31) $$

so that, substituting (31) into (30) leads to

$$ \nabla \cdot q^{n,(2)} = \nabla \cdot q^n + \frac{1}{\Delta t} (\rho^n - 1 - \Delta t \nabla \cdot q^n) + \mathcal{O}(\lambda^2) = \frac{1}{\Delta t} (\rho^n - 1) + \mathcal{O}(\lambda^2), $$

and $\nabla \cdot q^{n+1} = \nabla \cdot q^{n,(2)}$. But if the density $\rho^n$ is far from 1 (which is typically the case for $n = 0$ if the initial density is not well-prepared), thus we can only conclude $\nabla \cdot q^{n+1} = \mathcal{O}(1)$. So, as stated before, while the constraint on the density is satisfied after one time step, this is not the case for the constraint on the divergence of the velocity. We now consider the next time step, namely $(n+2)$. The first stage of this second time step gives as before $\rho^{n+1,(1)} = 1 + \mathcal{O}(\lambda^2)$, while for the electric potential one has

$$ \Delta \phi^{n+1,(1)} = \frac{\rho^{n+1,(1)} - 1}{\lambda^2} = \frac{1}{\Delta t} \nabla \cdot q^{n+1} + \mathcal{O}(\lambda^2), \quad (32) $$

and for the momentum

$$ q^{n+1,(1)} = q^{n+1} + \Delta t \nabla \phi^{n+1,(1)}. \quad (33) $$
The second stage instead becomes $\rho^{n+2,(2)} = 1 + \mathcal{O}(\lambda^2)$,

$$ \nabla \cdot q^{n+1,(2)} = \nabla \cdot q^{n+1} - \Delta t \nabla^2 : F^{n+1,(1)} + \Delta t \nabla \cdot ((\rho^{n+1,(1)} - 1)\nabla\phi^{n+1,(1)}) + \Delta t \Delta\phi^{n+1,(2)}, \quad (34) $$

and

$$ \Delta\phi^{n+1,(2)} = \nabla^2 : (u^{n+1,(1)} \otimes u^{n+1,(1)}) - \frac{1}{\Delta t} (\nabla \cdot q^{n+1}) + \mathcal{O}(\lambda^2), \quad (35) $$

so that, substituting (35) into (34) gives $\nabla \cdot q^{n+1,(2)} = \mathcal{O}(\lambda^2)$. Finally, going back to (34) and using the previous limits gives

$$ u^{n+1,(2)} = u^{n+1} - \Delta t \nabla \cdot (u^{n+1,(1)} \otimes u^{n+1,(1)}) + \Delta t \nabla \phi^{n+1,(2)}. \quad (36) $$

Finally, observing that the values computed at second stage coincide with the ones of the numerical solution at time $(n+2)$ concludes the proof. □

In conclusion, the semi-discrete Implicit-Explicit first order Runge-Kutta method (20)-(26) reduces to a scheme that is consistent with the asymptotic model (4)-(6) in two time steps when initial data are far from equilibrium state. That is opposite to what happens with the Explicit-Implicit Euler method of the previous section (14)-(16) which needs consistency of the initial data with the asymptotic state. The price to pay is a more expensive numerical method, since it involves two stages instead of the one stage scheme proposed before. In the next section we discuss high order extensions of such methods.

## 3.3 An IMEX Runge-Kutta class of methods for handling quasineutrality in plasma

We try now to generalize the approaches presented in the previous sections to the high order case. The mathematical framework is the one discussed in [29, 42, 39]. The class of IMEX-RK schemes provides a robust and efficient basis to design schemes which are uniformly stable with respect to fast scale dynamics. In the recent past, they have been extensively used [32, 33, 39, 16, 22] for treating singular perturbation problems with the specific task of obtaining high order numerical methods consistent with the limit model.

### 3.3.1 Generalities on IMEX Runge-Kutta methods

From now on we restrict ourselves to a specific subclass of IMEX-RK schemes, namely we rely on diagonally implicit (DIRK) schemes. This is a typical choice in many different context where computational cost is one of the main request. In general, an s-stage IMEX-RK scheme is characterised by the two $s \times s$ lower triangular matrices $\tilde{A} = (\tilde{a}_{i,j})$, and $A = (a_{i,j})$, the coefficients $\tilde{c} = (\tilde{c}_1, \tilde{c}_2, \dots, \tilde{c}_s)$ and $c = (c_1, c_2, \dots, c_s)$, and the weights $\tilde{\omega} = (\tilde{\omega}_1, \tilde{\omega}_2, \dots, \tilde{\omega}_s)$ and $\omega = (\omega_1, \omega_2, \dots, \omega_s)$. Here, for the DIRK class, the entries of $\tilde{A}$ and $A$ satisfy the conditions $\tilde{a}_{i,j} = 0$ for $j \ge i$, and $a_{i,j} = 0$ for $j > i$. In order to introduce the method, let us consider the following general stiff system of ODEs in additive form:

$$ y' = f(t, y) + \frac{1}{\varepsilon} g(t, y), \quad (37) $$

where $0 < \varepsilon \ll 1$ is called the stiffness parameter. The functions $f$ and $g$ are known as, respectively, the non-stiff part and the stiff part of the system (37) (see, e.g. [29], for a comprehensive treatment of such systems). Let now $y^n$ be a numerical solution of (37) at time $t^n$ and let $\Delta t$ denote a fixed time step. An s-stage IMEX-RK scheme updates $y^n$ to $y^{n+1}$ through the computation of s intermediate stages:

$$ Y^{(i)} = y^n + \Delta t \sum_{j=1}^{i-1} \tilde{a}_{i,j} f(t^n + \tilde{c}_j \Delta t, Y^{(j)}) + \Delta t \sum_{j=1}^{s} a_{i,j} \frac{1}{\varepsilon} g(t^n + c_j \Delta t, Y^{(j)}), \quad 1 \le i \le s, \quad (38) $$

$$ y^{n+1} = y^n + \Delta t \sum_{i=1}^{s} \tilde{\omega}_i f(t^n + \tilde{c}_i \Delta t, Y^{(i)}) + \Delta t \sum_{i=1}^{s} \omega_i \frac{1}{\varepsilon} g(t^n + c_i \Delta t, Y^{(i)}), \quad (39) $$

where the coefficients $\tilde{a}_{i,j}$, $a_{i,j}$, $\tilde{c}_i$, $c_i$ and $\tilde{\omega}_i$, $\omega_i$ are components of the so-called Butcher tableaux, with $\tilde{c}, \tilde{b}, c, b \in \mathbb{R}^s$

$$ \begin{array}{c|c|c} \tilde{c}^T & \tilde{A} & c^T \\ \hline \tilde{\omega}^T & & \omega^T \end{array} $$

<caption>Figure 3: General double Butcher tableaux for an IMEX-RK scheme.</caption>

and $\tilde{A}, A \in \mathcal{M}_{s,s}(\mathbb{R})$. As mentioned above, we consider in the following two subcategories of IMEX-RK namely, the type-A and type-CK schemes which are defined below (see [8, 22, 34] for details).

**Definition 3.** An IMEX-RK scheme is said to be of
• type-A, if the matrix A is invertible;

• type-CK, if the matrix $A \in \mathbb{R}^{s \times s}$, $s \ge 2$, can be written as

$$A = \begin{pmatrix} 0 & 0 \\ \alpha & A_{s-1} \end{pmatrix},$$

where $\alpha \in \mathbb{R}^{s-1}$ and $A_{s-1} \in \mathbb{R}^{s-1 \times s-1}$ is invertible.

Another property which is crucial to obtain the AP property (from Definition 2) is that the IMEX-RK scheme is globally stiffly accurate.

**Definition 4.** An IMEX-RK scheme with the Butcher tableaux given in Figure 3 is said to be globally stiffly accurate (GSA), if

$$\tilde{a}_{s,j} = \tilde{\omega}_j, \quad a_{s,j} = \omega_j \text{ for all } j = 1, \dots, s. \tag{40}$$

This property is enough to guarantee that the solution obtained from the last stage of the IMEX-RK method coincides with the numerical solution $y^{n+1} = Y^{(s)}$. In the next paragraph, we make use of the scheme (38)-(39) for the Euler-Poisson system in the quasineutral scaling.

### 3.3.2 High order time semi-discretization for the Euler-Poisson system

Inspired by the scheme in (13) (using the notations (11)) and by the two first order discretizations discussed in the previous section, we introduce now the following class of numerical schemes

$$U^{(i)} = U^n - \Delta t \sum_{j=1}^{i-1} \left[ \tilde{a}_{i,j} G_{\text{ex}}^{(j)} + a_{i,j} G_{\text{im}}^{(j)} \right] - \Delta t a_{i,i} G_{\text{im}}^{(i)}, \tag{41}$$

$$\lambda^2 \Delta \phi^{(i)} = \rho^{(i)} - 1 \tag{42}$$

where the intermediate RK stages are, for $i = 1, \dots, s$, and the final update at time $t^{n+1} = t^n + \Delta t$ is:

$$U^{n+1} = U^n - \Delta t \sum_{j=1}^{s} \left[ \tilde{\omega}_j G_{\text{ex}}^{(j)} + \omega_j G_{\text{im}}^{(j)} \right], \tag{43}$$

$$\lambda^2 \Delta \phi^{n+1} = \rho^{n+1} - 1, \tag{44}$$

where the coefficients $\tilde{a}_{i,j}$, $a_{i,j}$ and $\tilde{\omega}_j$, $\omega_j$ are usually prescribed using the Butcher tableaux of Figure 3 involving the matrices $\tilde{A}$ and $A$ as discussed in the previously. As mentioned before in our analysis, we will only use DIRK-GSA schemes. A direct consequence of this latter choice is the last stage $s$ of (41)-(42) is the same as (43)-(44). Therefore, in the sequel, we will assume the stage $s$ corresponds to the iteration update i.e. $(U^{n+1}, \phi^{n+1}) = (U^{(s)}, \phi^{(s)})$. Now, using the components $U = (\rho, q)$, the scheme (41)-(42) can be written, for $i = 1, \dots, s$, as

$$\rho^{(i)} = \rho^n - \Delta t \sum_{j=1}^{i-1} a_{i,j} \nabla \cdot q^{(j)} - \Delta t a_{i,i} \nabla \cdot q^{(i)} = \hat{\rho}^{(i)} - \Delta t a_{i,i} \nabla \cdot q^{(i)}, \tag{45}$$

$$\begin{aligned} q^{(i)} &= q^n - \Delta t \sum_{j=1}^{i-1} \left[ \tilde{a}_{i,j} \left( \nabla \cdot F^{(j)} - (\rho^{(j)} - 1) \nabla \phi^{(j)} \right) - a_{i,j} \nabla \phi^{(j)} \right] + \Delta t a_{i,i} \nabla \phi^{(i)} = \\ &= \hat{q}^{(i)} + \Delta t a_{i,i} \nabla \phi^{(i)}, \end{aligned} \tag{46}$$

$$\lambda^2 \Delta \phi^{(i)} = \rho^{(i)} - 1, \tag{47}$$

where the notations $\hat{\rho}^{(i)}$ and $\hat{q}^{(i)}$ are defined by

$$\hat{\rho}^{(i)} = \rho^n - \Delta t \sum_{j=1}^{i-1} a_{i,j} \nabla \cdot q^{(j)}, \quad \hat{q}^{(i)} = q^n - \Delta t \sum_{j=1}^{i-1} \left[ \tilde{a}_{i,j} \left( \nabla \cdot F^{(j)} - (\rho^{(j)} - 1) \nabla \phi^{(j)} \right) - a_{i,j} \nabla \phi^{(j)} \right]. \tag{48}$$

Written under this form, the scheme looks like an implicit scheme, involving the resolution of non linear systems. However, in the sequel, as done for the first order schemes presented in the previous section, we manipulate the equations showing that the proposed class of schemes does not require any nonlinear solver. It instead exhibits the same complexity of a fully explicit scheme. To this end, taking the divergence of the momentum equation (46) gives

$$\begin{aligned} \nabla \cdot q^{(i)} &= \nabla \cdot \hat{q}^{(i)} + \Delta t a_{i,i} \Delta \phi^{(i)}, \\ &= \nabla \cdot q^n - \Delta t \sum_{j=1}^{i-1} \left[ \tilde{a}_{i,j} \left( \nabla^2 : F^{(j)} - \nabla \cdot \left( (\rho^{(j)} - 1) \nabla \phi^{(j)} \right) \right) - a_{i,j} \Delta \phi^{(j)} \right] + \Delta t a_{i,i} \Delta \phi^{(i)}, \end{aligned} \tag{49}$$
and inserting the above expression for $\nabla \cdot q^{(i)}$ in the mass update (45) gives,

$$
\begin{align*}
\rho^{(i)} &= \hat{\rho}^{(i)} - \Delta t a_{i,i} \left\{ \nabla \cdot q^n - \Delta t \sum_{j=1}^{i-1} \left[ \tilde{a}_{i,j} \left( \nabla^2 : F^{(j)} - \nabla \cdot \left( (\rho^{(j)} - 1) \nabla \phi^{(j)} \right) \right) - a_{i,j} \Delta \phi^{(j)} \right] + \Delta t a_{i,i} \Delta \phi^{(i)} \right\}, \\
&= \hat{\rho}^{(i)} - \Delta t a_{i,i} \nabla \cdot q^n \\
&\quad + \Delta t^2 a_{i,i} \sum_{j=1}^{i-1} \left\{ \tilde{a}_{i,j} \left( \nabla^2 : F^{(j)} - \nabla \cdot \left( (\rho^{(j)} - 1) \nabla \phi^{(j)} \right) \right) - a_{i,j} \Delta \phi^{(j)} \right\} - \Delta t^2 a_{i,i}^2 \Delta \phi^{(i)}.
\end{align*}
$$

Using the Poisson equation $\lambda^2 \Delta \phi^{(i)} = \rho^{(i)} - 1$ to substitute for $\Delta \phi^{(i)}$ and subtracting 1 from both sides in the above equation gives,

$$
\begin{equation}
\begin{split}
\rho^{(i)} = {}& 1 + \frac{\lambda^2}{\lambda^2 + \Delta t^2 a_{i,i}^2} \left\{ \rho^n - 1 - \Delta t \sum_{j=1}^{i-1} \tilde{a}_{i,j} \nabla \cdot q^{(j)} - \Delta t a_{i,i} \nabla \cdot q^n \right. \\
& \left. \qquad + \Delta t^2 a_{i,i} \sum_{j=1}^{i-1} \left\{ \tilde{a}_{i,j} \left( \nabla^2 : F^{(j)} - \nabla \cdot \left( (\rho^{(j)} - 1) \nabla \phi^{(j)} \right) \right) - a_{i,j} \Delta \phi^{(j)} \right\} \right\}.
\end{split}
\tag{50}
\end{equation}
$$

Now, using notations $\hat{\rho}^{(i)}$ and $\hat{q}^{(i)}$ from (48) the update for the density reads:

$$
\rho^{(i)} = 1 + \frac{\lambda^2}{\lambda^2 + \Delta t^2 a_{i,i}^2} \left( \hat{\rho}^{(i)} - 1 - \nabla \cdot \hat{q}^{(i)} \right), \quad (51)
$$

while for the electric potential

$$
\Delta \phi^{(i)} = \frac{1}{\lambda^2 + \Delta t^2 a_{i,i}^2} \left( \hat{\rho}^{(i)} - 1 - \nabla \cdot \hat{q}^{(i)} \right). \quad (52)
$$

**Remark 3.** Note that in the first stage the only possible coefficient which may be different from zero is $a_{1,1}$ and so it may be written separately from the rest of the stages. Indeed, for $i = 1$, (50) becomes (using $\hat{\rho}^{(i)} = \rho^n$ and $\hat{q}^{(i)} = q^n$)

$$
\rho^{(1)} = 1 + \frac{\lambda^2}{\lambda^2 + \Delta t^2 a_{1,1}^2} \left( \rho^n - 1 - a_{1,1} \Delta t \nabla \cdot q^n \right) \quad (53)
$$

$$
q^{(1)} = q^n + a_{1,1} \Delta t \nabla \phi^{(1)}, \quad (54)
$$

$$
\lambda^2 \Delta \phi^{(1)} = \rho^{(1)} - 1. \tag{55}
$$

It has to be noted that in general $a_{1,1}$ is different from zero only for type-A IMEX scheme while for type-CK IMEX method one has $a_{1,1} = 0$.

The schemes can be summarised in the following Algorithm 1.

<caption>**Algorithm 1** High order AP schemes for the Euler-Poisson system in the quasineutral scaling</caption>

1: **for** i = 1, 2, ..., s **do**

2: compute $\rho^{(i)}$ using (51)

3: compute $\phi^{(i)}$ using (52)

4: compute $q^{(i)}$ using (46)

5: **end for**

In the following section we discuss some relevant properties of the schemes presented.

4 Analysis of the penalized-IMEX schemes

The goal of this section is to investigate the behaviour of the penalized-IMEX schemes (41)-(44) when $\lambda \to 0$, for a time step $\Delta t > 0$ fixed. In particular, we will focus on type-A schemes, a short discussion on type-CK schemes will end the section. Let us notice that an example of type-A scheme has been given in Section 3.2 while an example of type-CK scheme has been given in Section 3.1.
## 4.1 Analysis of type-A IMEX Runge-Kutta schemes

In this part, we study the AP property of the penalized-IMEX schemes combined with type-A Runge-Kutta method in the quasineutral limit for the Euler-Poisson system. We divide the following analysis into two parts. First, we check that the constraints

$$ \rho^n - 1 = \mathcal{O}(\lambda^2), \quad \nabla \cdot u^n = \mathcal{O}(\lambda^2), \quad \phi^n = \mathcal{O}(1), \quad n \ge 2, \qquad (56) $$

are satisfied. This will ensure the compatibility of the numerical solution with the asymptotic model (4)-(6). Second, we prove that the penalized-IMEX scheme degenerates as $\lambda \to 0$ towards a high order time discretization of the asymptotic model (4)-(6).

The next result is about the projection of the unknowns of the Euler-Poisson system to the quasineutral equilibrium, i.e. $\rho = 1$ and $\nabla \cdot u = 0$ when general initial conditions are considered.

**Theorem 3.** The s-stage type-A IMEX scheme (51)-(52)-(46) is a consistent time approximation of (1)-(2)-(3). Moreover, for any initial conditions $\rho^0 = \rho(t=0)$, $q^0 = q(t=0)$ with $\lambda^2\Delta\phi^0 = \rho^0 - 1$, as $\lambda \to 0$, the scheme satisfies the following properties

$$ \rho^n - 1 = \mathcal{O}(\lambda^2), \quad \nabla \cdot u^n = \mathcal{O}(\lambda^2), \quad \phi^n = \mathcal{O}(1), \quad n \ge 2. $$

*Proof.* The proof is separated into two parts corresponding to the first two iterations $n = 0$ and $n = 1$. These two parts being also divided into two sub-parts corresponding to the first stage $i = 1$ and the others $i = 2, \dots, s$.

**1. First iteration (from $n = 0$ to $n = 1$).** We investigate the behaviour of the scheme after the first time iteration. For general initial conditions, we verify that the type-A Runge-Kutta schemes enables the projection of the numerical solution solution to the quasineutral manifold for what concerns the density $\rho$.

### First stage ($i = 1$)

For the first stage for prescribed $\rho, q$ we have:

$$ \begin{aligned} \rho^{(1)} &= \rho^n - \Delta t a_{1,1} \nabla \cdot q^{(1)} \\ q^{(1)} &= q^n + \Delta t a_{1,1} \nabla \phi^{(1)}. \end{aligned} \qquad (57) $$

Now, with the same manipulations performed in the previous section, we have for $\rho$ via (53):

$$ \rho^{(1)} = 1 + \frac{\lambda^2}{\lambda^2 + \Delta t^2 a_{1,1}^2} (\rho^n - 1 - \Delta t a_{1,1} \nabla \cdot q^n), \qquad (58) $$

which ensures $\rho^{(1)} = 1 + \mathcal{O}(\lambda^2)$ for any (smooth) initial conditions. Hence, the Poisson equation (52) becomes

$$ \Delta\phi^{(1)} = \frac{1}{\lambda^2 + \Delta t^2 a_{1,1}^2} (\rho^n - 1 - \Delta t a_{1,1} \nabla \cdot q^n), \qquad (59) $$

which ensures that $\phi^{(1)} = \mathcal{O}(1)$ for any initial conditions. Taking the divergence of (57) and using the expression (59) for $\Delta\phi^{(1)}$ gives

$$ \nabla \cdot q^{(1)} = \nabla \cdot q^n + \Delta t a_{1,1} \Delta\phi^{(1)} = \nabla \cdot q^n + \frac{\Delta t a_{1,1}}{\lambda^2 + \Delta t^2 a_{1,1}^2} (\rho^n - 1 - \Delta t a_{1,1} \nabla \cdot q^n), $$

which in turns gives

$$ \nabla \cdot q^{(1)} = \left( \frac{\lambda^2}{\lambda^2 + \Delta t^2 a_{1,1}^2} \right) \nabla \cdot q^n + \frac{\Delta t a_{1,1}}{\lambda^2 + \Delta t^2 a_{1,1}^2} (\rho^n - 1) = \mathcal{O}(\lambda^2) + \frac{\Delta t a_{1,1}}{\lambda^2 + \Delta t^2 a_{1,1}^2} (\rho^n - 1). \qquad (60) $$

From the above computations, one can see that the density is projected on the constant state $\rho(1) = 1$ so that the potential remains bounded, however, the divergence of the momentum is not projected to 0 at the end of the first stage if we start with a general (non well-prepared) initial condition.

We show now by induction that these three conditions, namely $\rho^{(i)} - 1 = \mathcal{O}(\lambda^2)$, $\phi^{(i)} = \mathcal{O}(1)$ and $\nabla \cdot q^{(i)} = \mathcal{O}(1)$ are propagated through all the stages of the first time iteration.
### Stages $i > 1$

Let us assume that for the stages $j = 1, \dots, i - 1$, the following conditions hold true

$$ \rho^{(j)} - 1 = \mathcal{O}(\lambda^2), \quad \phi^{(j)} = \mathcal{O}(1), \quad \nabla \cdot q^{(j)} = \mathcal{O}(1), \qquad (61) $$

then we want to prove that these conditions remain true for the stage $i$. The density update at stage $i$ follows equation (51) and consequently gives:

$$ \rho^{(i)} = \hat{\rho}^{(i)} - \Delta ta_{i,i} \nabla \cdot q^{(i)} = \hat{\rho}^{(i)} - \Delta ta_{i,i} \nabla \cdot \hat{q}^{(i)} - \frac{\Delta t^2 a_{i,i}^2}{\lambda^2} (\rho^{(i)} - 1), $$

so that

$$ \rho^{(i)} - 1 = \frac{\lambda^2}{\lambda^2 + \Delta t^2 a_{i,i}^2} \left( \hat{\rho}^{(i)} - 1 - \Delta ta_{i,i} \nabla \cdot \hat{q}^{(i)} \right). \qquad (62) $$

Now given the expressions of $\hat{\rho}^{(i)}$ and $\hat{q}^{(i)}$ from the definitions (48) and using the induction assumption (61), we conclude that $\rho^{(i)} - 1 = \mathcal{O}(\lambda^2)$ and then the elliptic problem (47) for $\phi^{(i)}$ give $\phi^{(i)} = \mathcal{O}(1)$. Let us now focus on $\nabla \cdot q^{(i)}$. We have that

$$ \nabla \cdot q^{(i)} = \nabla \cdot \hat{q}^{(i)} + \Delta ta_{i,i} \Delta \phi^{(i)}, $$

and thanks to the induction assumption (61), we easily deduce $\nabla \cdot q^{(i)} = \mathcal{O}(1)$.

## 2. Second iteration (from $n = 1$ to $n = 2$). After the first iteration, we have that the numerical solution ($\rho^1, q^1, \phi^1$), at time $t^1$ ($t^1 = t^0 + \Delta t$), satisfies the following conditions:

$$ \rho^1 - 1 = \rho^{(s)} - 1 = \mathcal{O}(\lambda^2), \quad \phi^1 = \phi^{(s)} = \mathcal{O}(1), \quad \nabla \cdot q^1 = \nabla \cdot q^{(s)} = \mathcal{O}(1), \qquad (63) $$

where $(s)$ is the last stage of the previous Runge-Kutta iteration corresponding to the numerical solution at time $n + 1 = 1$ since the schemes are globally stiffly accurate by hypothesis.

### First stage ($i = 1$)

We analyse now the first stage. Since the expressions obtained for $\rho^{(1)}$ in (58) and $\phi^{(1)}$ in (59) are valid for any iteration $n \ge 0$, we get

$$ \rho^{(1)} - 1 = \mathcal{O}(\lambda^2) \quad \text{and} \quad \phi^{(1)} = \mathcal{O}(1). $$

The divergence of the momentum instead verifies

$$ \nabla \cdot q^{(1)} = \mathcal{O}(\lambda^2) + \frac{\Delta ta_{1,1}}{\lambda^2 + \Delta t^2 a_{1,1}^2} (\rho^n - 1). $$

However now, since from the first iteration we have that $(\rho^n - 1) = \rho^1 - 1 = \mathcal{O}(\lambda^2)$, we get that $\nabla \cdot q^{(1)} = \mathcal{O}(\lambda^2)$ where we recall that the stage values indicated with the parenthesis correspond to the solutions computed between the time level $n = 1$ and $n = 2$. To summarise, we obtain after the first stage:

$$ \rho^{(1)} - 1 = \mathcal{O}(\lambda^2), \quad \phi^{(1)} = \mathcal{O}(1), \quad \text{and} \quad \nabla \cdot q^{(1)} = \mathcal{O}(\lambda^2). $$

### Stages $i > 1$

Let us assume that for the stages $j = 1, \dots, i - 1$ the following conditions hold true

$$ \rho^{(j)} - 1 = \mathcal{O}(\lambda^2), \quad \phi^{(j)} = \mathcal{O}(1), \quad \nabla \cdot q^{(j)} = \mathcal{O}(\lambda^2), \qquad (64) $$

then we want to prove, as done for the first time step, that these conditions remain true for the stage $i$. From the density update shown for the previous time step (62) which holds true also here we have that $\rho^{(i)} - 1 = \mathcal{O}(\lambda^2)$. In the same way, using this result for the density in the Poisson equation (47) we obtain $\phi^{(i)} = \mathcal{O}(1)$. Let us now focus on the proving last property i.e. $\nabla \cdot q^{(i)} = \mathcal{O}(\lambda^2)$. From the divergence expression (49) and the induction assumption (64), we get

$$ \begin{aligned} \nabla \cdot q^{(i)} &= \nabla \cdot \hat{q}^{(i)} + \Delta ta_{i,i} \Delta \phi^{(i)} \\ &= \nabla \cdot \hat{q}^{(i)} + \frac{\Delta ta_{i,i}}{\lambda^2} (\rho^{(i)} - 1) \\ &= \nabla \cdot \hat{q}^{(i)} + \frac{\Delta ta_{i,i}}{\lambda^2} (\hat{\rho}^{(i)} - 1 - a_{i,i} \Delta t \nabla \cdot q^{(i)}), \\ (1 + \frac{a_{i,i}^2 \Delta t^2}{\lambda^2}) \nabla \cdot q^{(i)} &= \nabla \cdot \hat{q}^{(i)} + \frac{\Delta ta_{i,i}}{\lambda^2} (\hat{\rho}^{(i)} - 1). \end{aligned} $$
Hence we have using (48) for $\hat{\rho}^{(i)}$

$$ 
\begin{aligned}
\nabla \cdot q^{(i)} &= \frac{\lambda^2}{\lambda^2 + a_{i,i}^2 \Delta t^2} \left[ \nabla \cdot \hat{q}^{(i)} + \frac{\Delta t a_{i,i}}{\lambda^2} (\hat{\rho}^{(i)} - 1) \right] \\
&= \frac{\lambda^2}{\lambda^2 + a_{i,i}^2 \Delta t^2} \nabla \cdot \hat{q}^{(i)} + \frac{a_{i,i} \Delta t}{\lambda^2 + a_{i,i}^2 \Delta t^2} (\rho^n - 1) - \frac{a_{i,i} \Delta t^2}{\lambda^2 + a_{i,i}^2 \Delta t^2} \left( \sum_{j=1}^{i-1} a_{i,j} \nabla \cdot q^{(j)} \right).
\end{aligned}
$$

Now, $\hat{q}^{(i)}$ is uniformly bounded and $\nabla \cdot \hat{q}^{(i)}$ also. Second, from (64), we have $\nabla \cdot q^{(j)} = \mathcal{O}(\lambda^2)$ for $j = 1, \dots, i-1$ and $\rho^n - 1 = \mathcal{O}(\lambda^2)$. Hence, we conclude that $\nabla \cdot q^{(i)} = \mathcal{O}(\lambda^2)$. $\square$

We discuss now the capability of the type-A IMEX Runge-Kutta schemes (50)-(47)-(46) to become a consistent and high order discretization of the limit model (4)-(6) when $\lambda \to 0$. The following result holds true

**Theorem 4.** For any initial conditions $\rho^0 = \rho(t=0)$, $q^0 = q(t=0)$ with $\lambda^2 \Delta \phi^0 = \rho^0 - 1$, as $\lambda \to 0$, the s-stage type-A IMEX Runge-Kutta schemes (50)-(47)-(46) degenerate towards

$$ u^{(i)} = u^n - \Delta t \sum_{j=1}^{i-1} \tilde{a}_{i,j} \nabla^2 : (u^{(j)} \otimes u^{(j)}) + \Delta t \sum_{j=1}^{i} a_{i,j} \nabla \phi^{(j)}, \quad \text{and} \quad \nabla \cdot u^{(i)} = 0 \quad \text{for } i=1, \dots, s, $$

which is an s-stage high order IMEX time discretization of the asymptotic model (4)-(6).

*Proof.* From Theorem 3, we have that $\rho^{(i)} - 1 = \mathcal{O}(\lambda^2)$ and

$$ q^{(i)} = \rho^{(i)} u^{(i)} = u^{(i)} + \mathcal{O}(\lambda^2), $$

while for the momentum flux the following result holds true (recalling $p(\rho) = C\rho^\beta$)

$$ F^{(i)} \equiv \rho^{(i)} u^{(i)} \otimes u^{(i)} + p^{(i)} I = u^{(i)} \otimes u^{(i)} + I + \mathcal{O}(\lambda^2). $$

Thus, using (46), the asymptotic scheme for the mean velocity $u$ becomes

$$ 
\begin{aligned}
u^{(i)} &= u^n - \Delta t \sum_{j=1}^{i-1} \left[ \tilde{a}_{i,j} \nabla \cdot (u^{(j)} \otimes u^{(j)}) - a_{i,j} \nabla \phi^{(j)} \right] + \Delta t a_{i,i} \nabla \phi^{(i)} + \mathcal{O}(\lambda^2), \\
&= u^n - \Delta t \sum_{j=1}^{i-1} \tilde{a}_{i,j} \nabla \cdot (u^{(j)} \otimes u^{(j)}) + \Delta t \sum_{j=1}^{i} a_{i,j} \nabla \phi^{(j)} + \mathcal{O}(\lambda^2),
\end{aligned}
\quad (65) $$

which is nothing else but an IMEX-RK scheme for the asymptotic quasineutral model (5) on $u$. Now, let us determine the equation for the electric potential $\phi$. By considering the divergence of the momentum equation (46) in combination with Theorem 3 we get

$$ 
\begin{aligned}
-\Delta t a_{i,i} \Delta \phi^{(i)} &= \nabla \cdot q^n - \Delta t \sum_{j=1}^{i-1} \left[ \tilde{a}_{i,j} \left( \nabla^2 : F^{(j)} - \nabla \cdot \left( (\rho^{(j)} - 1) \nabla \phi^{(j)} \right) \right) - a_{i,j} \Delta \phi^{(j)} \right] + \mathcal{O}(\lambda^2) \\
&= -\Delta t \sum_{j=1}^{i-1} \tilde{a}_{i,j} \nabla^2 : (u^{(j)} \otimes u^{(j)}) + \Delta t \sum_{j=1}^{i-1} a_{i,j} \Delta \phi^{(j)} + \mathcal{O}(\lambda^2).
\end{aligned}
\quad (66) $$

Then, we proved the asymptotic scheme (65)-(66) is a s-stage IMEX scheme for the asymptotic model (4)-(5). $\square$

**Remark 4.** Note that for $i=1$, i.e. for the first stage of each time iteration, the expression (66) becomes

$$ \Delta \phi^{(1)} = \mathcal{O}(\lambda^2), $$

however, for the next stages, $\phi^{(j)} = \mathcal{O}(1)$, as one can see for the case $i=2$,

$$ a_{2,2} \Delta \phi^{(2)} = \tilde{a}_{2,2} \nabla^2 : (u^{(1)} \otimes u^{(1)}) + \mathcal{O}(\lambda^2). $$
## 4.2 Analysis of type-CK IMEX Runge-Kutta schemes

In this part, we quickly discuss how type-CK IMEX schemes behave in the quasineutral limit for the Euler-Poisson system. One leading example of such class of schemes has been given in Section 3.1 as already stated. The only difference with type-A schemes lies in the first stage ($i = 1$) where type-CK IMEX schemes give

$$ \rho^{(1)} = \rho^n, \quad q^{(1)} = q^n, \quad \lambda^2 \Delta\phi^{(1)} = \rho^n - 1, $$

being all coefficients of the two tableaux equal to zero. Thus, one can immediately appreciate the difference with respect to the type-A schemes (studied in the previous subsection) for which the first stage presents always non zero coefficients for the implicit part. Consequently, with this class of type-A methods, one is able to project the density towards the quasineutral manifold during the first stage, i.e. $\rho^{(1)} - 1 = O(\lambda^2)$. This in turn ensures for the potential to remain bounded, i.e. $\phi^{(1)} = O(1)$. As opposite, when the initial data are not well-prepared as discussed in Definition 1, for the type-CK schemes, one has that $\phi^{(1)} = \phi^n = O(1/\lambda^2)$ which tends to infinity as $\lambda \to 0$. We have consequently the following general result.

**Theorem 5.** For well-prepared initial data definition: $\rho^0 - 1 = O(\lambda^2) = \nabla \cdot q^0$, $\phi^0 = O(1)$, the s-stage type-CK IMEX Runge-Kutta schemes (50)-(47)-(46) degenerate towards a s-stage type-CK IMEX scheme for the asymptotic model (4)-(5)-(6). Conversely, for general initial data, the s-stage type-CK IMEX Runge-Kutta schemes (50)-(46)-(47) satisfies $\phi^{(i)} = O(1/\lambda^2)$.

*Proof.* For the well-prepared case by following the path of Theorem 4, one has $\rho^{(i)} - 1 = O(\lambda^2)$, $q^{(i)} = u^{(i)} + O(\lambda^2)$, $\forall i$, $\Delta\phi^{(i)} = O(1)$ for stage value ($i = 1$) and then for ($i > 1$) $\Delta\phi^{(i)}$ can be computed from (66). This is enough to prove that $u^{(i)}$ is computed from (65) starting from the first time step $n = 1$ and thus to prove that (50)-(47)-(46) realize a consistent and high order discretization of the quasineutral limit (4)-(5). For general initial conditions, we have $\rho^{(1)} - 1 = O(1)$ so immediately the electric potential becomes $\phi^{(1)} = O(1/\lambda^2)$. By induction, one can easily show using (50) that these two estimates are propagated along the stages. □

The consequence of the above result is that for initial data far from the quasineutral regime, the use of type-CK schemes may bring to instabilities in the numerical solution when positive and sufficiently small values of the Debye length $\lambda$ are employed. In the numerical section we will show that this is indeed the case.

# 5 Space discretization

In this section, we describe the spatial discretization strategy coupled with the time semi-discrete penalized-IMEX scheme introduced in Section 3.3 and we present the space-time fully-discrete numerical scheme. We follow a similar path to the one used in [3], where, a second order fully-discrete scheme for the Euler-Poisson system was developed using a combination of Rusanov fluxes and simple central differencing. Higher order spatial discretization can be realised by combining standard high order shock-capturing schemes with WENO reconstructions and high order central differences. The final scope is to obtain a fully-discrete scheme which is compatible enough to maintain the properties of the time semi-discrete penalized-IMEX scheme presented in Section 3.3.

In order to present the fully-discrete scheme, we introduce the vector $k = (k_1, k_2, \dots, k_d)$ where each $k_m$ for $m = 1, 2, \dots, d$ represents the space mesh in direction $m$: $x_m = k_m \Delta x_m$ with $\Delta x_m = L_m/N_m$, $x_m \in [0, L_m]$ and $N_m$ the number of space points in each direction. We further introduce the following finite difference and averaging operators in the $x_m$-direction

$$ \delta_{x_m} w_k = w_{k+\frac{1}{2}e_m} - w_{k-\frac{1}{2}e_m}, \quad \mu_{x_m} w_k = \frac{w_{k+\frac{1}{2}e_m} + w_{k-\frac{1}{2}e_m}}{2}, \qquad (67) $$

for any grid function $w_k$ and where $e_m$ indicates the basis vector of $\mathbb{R}^d$. The full discretization is introduced in the following definition.

**Definition 5.** The $i^{th}$ stage of the s-stage fully-discrete penalized-IMEX scheme for the Euler-Poisson system is defined as follows. The explicit variables are computed as:

$$ \hat{\rho}_k^{(i)} = \rho_k^n - \sum_{j=1}^{i-1} a_{i,j} \sum_{m=1}^{d} \nu_m \delta_{x_m} \mathcal{G}_{m,k}^{(i)}, \qquad (68) $$

$$ \hat{q}_k^{(i)} = q_k^n - \sum_{j=1}^{i-1} \tilde{a}_{i,j} \left( \sum_{m=1}^{d} \nu_m \delta_{x_m} \mathcal{F}_{m,k}^{(j)} + (\rho_k^{(j)} - 1) \nu_m \delta_{x_m} \mu_{x_m} \phi_k^{(j)} e_m \right) + \sum_{j=1}^{i-1} a_{i,j} \nu_m \delta_{x_m} \mu_{x_m} \phi_k^{(j)} e_m, \qquad (69) $$
where $\nu_m := \frac{\Delta t}{\Delta x_m}$ denote the mesh ratios, the vector $G_m$ is an approximation of the component mass flux $q_m$ while the vector $F_m$ is an approximation of the $m^{th}$ element of the flux $F$ in (11) defined as:

$$ 
\begin{aligned}
G_{m,k+\frac{1}{2}e_m}^{(j)} &= \frac{1}{2}(q_{m,k+e_m}^{(j)} + q_{m,k}^{(j)}) \\
F_{m,k+\frac{1}{2}e_m}^{(j)} &= \frac{1}{2}(F_m(U_{k+\frac{1}{2}e_m}^{(j),+}) + F_m(U_{k+\frac{1}{2}e_m}^{(j),-})) - \frac{\alpha_{m,k+\frac{1}{2}e_m}^{(j)}}{2}(q_{k+\frac{1}{2}e_m}^{(j),+} - q_{k+\frac{1}{2}e_m}^{(j),-}),
\end{aligned}
\quad (70) $$

where $F_m$ denotes the $m^{th}$ element of the flux $F = (F_1, \dots, F_m)$ and where, for any conservative variable $w$, we denote by $w_{i+\frac{1}{2}e_m}^{\pm}$, the interpolated states obtained using the piecewise linear reconstructions. Following, the above update $\rho_k^{(i)}$ is then given by:

$$ \rho_k^{(i)} = 1 + \frac{\lambda^2}{\lambda^2 + \Delta t^2 a_{i,i}^2} \left( \hat{\rho}_k^{(i)} - 1 - \sum_{m=1}^{d} \nu_m \mu_{x_m} \delta_{x_m} \hat{q}_{m,k}^{(i)} \right). \quad (71) $$

Successively, the implicit solution is obtained by first solving the discrete elliptic problem which is a linear system for $\{\phi_i^{(k)}\}$ given by,

$$ \sum_{m=1}^{d} \frac{\delta_{x_m}}{\Delta x_m} \left( \frac{\delta_{x_m}}{\Delta x_m} \phi_k^{(i)} \right) = \rho_k^{(i)} - 1, \quad (72) $$

followed by the explicit evaluations of the momentum:

$$ q_{m,k}^{(i)} = \hat{q}_{m,k}^{(i)} - \Delta t a_{i,i} \nu_m \delta_{x_m} \phi_k^{(i)}, \quad m = 1, 2, \dots, d, \quad (73) $$

where the repeated index $m$ takes values in $\{1, 2, \dots, d\}$. Finally, the wave-speeds are computed as follows

$$ \alpha_{m,k+\frac{1}{2}e_m}^{(j)} := 2 \max \left( \left| \frac{q_{m,k+\frac{1}{2}e_m}^{(j),-}}{\rho_{k+\frac{1}{2}e_m}^{(j),-}} \right|, \left| \frac{q_{m,k+\frac{1}{2}e_m}^{(j),+}}{\rho_{k+\frac{1}{2}e_m}^{(j),+}} \right| \right). \quad (74) $$

Let us observe that while the momentum flux $F$ is approximated by a Rusanov-type scheme instead the mass flux is approximated using simple central differences. The eigenvalues of the Jacobians flux which is approximated by Rusanov-type flux are

$$ \lambda_{m,1} = 0, \quad \lambda_{m,2} = \lambda_{m,3} = \frac{q_m}{\rho}, \quad \lambda_{m,4} = 2 \frac{q_m}{\rho}, $$

while the CFL condition is given, at time $t^n$, by:

$$ \Delta t \max_k \max_m \left( \frac{|\lambda_{m,3,k}|}{\Delta x_m}, \frac{|\lambda_{m,4,k}|}{\Delta x_m} \right) = \nu, \quad (75) $$

where $\nu < 1$ is the given CFL number, we refer to [28] for details about the choice of the CFL conditions.

# 6 Numerical results

In this section, we present the results of several numerical test cases carried out to corroborate the theoretical results of the previous sections. We specifically consider the two butcher tableaux of Figure 4 corresponding to two type-A GSA IMEX-RK schemes from [22] and one type-CK GSA IMEX scheme from [39], see Figure 5, applied to the penalized-IMEX schemes (50)-(47)-(46) presented in this work. Comparisons with the SI-IMEX scheme from [3] are also performed. This SI-IMEX scheme (based on the semi-implicit approach proposed in [7]) enables to design high order AP scheme in the quasineutral limit, but requires the initial condition to be well-prepared (weak AP property). We first investigate 1D test with a pressure law $\beta = 1$ whereas a 2D test will end this section with $\beta = 2$.

## 6.1 Asymptotic order of convergence

We start by a numerical validation of the asymptotic order of convergence of the numerical solution produced by the penalized-IMEX schemes. The exact solution of the incompressible limit (4)-(6) is used as a reference solution. The initial conditions reads:

$$ \rho(0, x_1) = 1, \quad u_1(0, x_1) = 1 + \delta^2 \cos(2K\pi x_1) \quad \phi(0, x_1) = 0, \quad (76) $$

where $x_1 \in [0, 20]$, $K = 1$ and the perturbation parameter is $\delta = 10^{-2}$. We consider three different values of $\lambda \in \{10^{-4}, 10^{-5}, 10^{-6}\}$, all close to the asymptotic regime. The numerical asymptotic order of convergence (AOC) is reported in Table 1 for the type-A scheme DP2A2. The results show that the uniform second order convergence of the penalized-IMEX scheme in the QN regime is achieved.
<table><caption>Figure 4: Double Butcher tableaux for two type-A GSA additive IMEX schemes. Top: DP1A2 (for DP1-A(2, 4, 2)). Bottom: DP2A2 (for DP2-A(2, 4, 2)).</caption><tbody><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1/2</td><td>1/2</td><td>0</td><td>0</td><td>0</td></tr><tr><td>1/3</td><td>0</td><td>0</td><td>0</td><td>0</td><td>2/3</td><td>1/6</td><td>1/2</td><td>0</td><td>0</td></tr><tr><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>1/2</td><td>-1/2</td><td>1/2</td><td>1/2</td><td>0</td></tr><tr><td>1</td><td>1/2</td><td>0</td><td>1/2</td><td>0</td><td>1</td><td>3/2</td><td>-3/2</td><td>1/2</td><td>1/2</td></tr><tr><td></td><td>1/2</td><td>0</td><td>1/2</td><td>0</td><td></td><td>3/2</td><td>-3/2</td><td>1/2</td><td>1/2</td></tr><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>γ</td><td>γ</td><td>0</td><td>0</td><td>0</td></tr><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-γ</td><td>γ</td><td>0</td><td>0</td></tr><tr><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td><td>0</td><td>1 - γ</td><td>γ</td><td>0</td></tr><tr><td>1</td><td>0</td><td>1/2</td><td>1/2</td><td>0</td><td>1</td><td>0</td><td>1/2</td><td>1/2 - γ</td><td>γ</td></tr><tr><td>0</td><td>0</td><td>1/2</td><td>1/2</td><td>0</td><td></td><td>0</td><td>1/2</td><td>1/2 - γ</td><td>γ</td></tr></tbody></table>

<table><caption>Figure 5: Double Butcher tableaux for a type-CK GSA additive IMEX schemes: ARS (for ARS(2,2,2)). Here, $\gamma = 1 - \sqrt{2}/2$, $\sigma = 1/2\gamma$ and $\delta = 1 - \sigma$.</caption><tbody><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>γ</td><td>γ</td><td>0</td><td>0</td><td>γ</td><td>0</td><td>γ</td><td>0</td></tr><tr><td>1</td><td>δ</td><td>1 - δ</td><td>0</td><td>1</td><td>0</td><td>1 - γ</td><td>γ</td></tr><tr><td></td><td>δ</td><td>1 - δ</td><td>0</td><td></td><td>0</td><td>1 - γ</td><td>γ</td></tr></tbody></table>

<table><caption>Table 1: Asymptotic order of convergence for the penalized-IMEX (DP2A2) scheme for $\lambda = 10^{-4}, 10^{-5}, 10^{-6}$.</caption><thead><tr><th rowspan="2">N</th><th colspan="2">λ = 10<sup>-4</sup></th><th colspan="2">λ = 10<sup>-5</sup></th><th colspan="2">λ = 10<sup>-6</sup></th></tr><tr><th>L<sup>2</sup> Error in φ</th><th>AOC</th><th>L<sup>2</sup> Error in φ</th><th>AOC</th><th>L<sup>2</sup> Error in φ</th><th>AOC</th></tr></thead><tbody><tr><td>320</td><td>3.2851E-03</td><td></td><td>3.2896E-03</td><td></td><td>3.2886E-03</td><td></td></tr><tr><td>640</td><td>8.1925E-04</td><td>2.00</td><td>8.2002E-04</td><td>2.00</td><td>8.2648E-04</td><td>1.99</td></tr><tr><td>1280</td><td>2.4011E-04</td><td>1.77</td><td>2.4018E-05</td><td>1.77</td><td>2.3955E-05</td><td>1.78</td></tr><tr><td>2560</td><td>6.2858E-05</td><td>1.93</td><td>6.2700E-06</td><td>1.93</td><td>8.0036E-06</td><td>1.58</td></tr></tbody></table>
## 6.2 Perturbation of a quasineutral plasma in 1D with $\rho = 1, u = 1$

A quasineutral state, characterized by a constant density and constant velocity, is considered:

$$ \bar{\rho}(x_1) = 1, \quad \bar{u}(x_1) = 1, \quad \bar{\phi}(x_1) = 0. \qquad (77) $$

The constant density, combined with the Poisson equation, implies that the electric potential is identically zero. This quasineutral state is perturbed by introducing a small cosine perturbation of magnitude $\delta^2$ to the uniform velocity field (see [15] for details). All computations for this test case are performed in the 1D spatial domain $x_1 \in [0, 1]$, with periodic boundary conditions applied to the hydrodynamic variables $\rho$ and $q$, and homogeneous Dirichlet boundary conditions imposed on the electric potential $\phi$. The purpose of this test is to verify the Asymptotic-Preserving (AP) property of the schemes (penalized-IMEX and SI-IMEX from [3]), i.e., their ability to recover a quasineutral background state independently of the spatial mesh size relative to the Debye length $\lambda$. To this end, we consider two types of perturbations added to the quasineutral state (77), resulting in two types of initial data:

* Case 1: well-prepared initial conditions $\delta^2 = \lambda^2$.

* Case 2: non well-prepared initial conditions $\delta^2 = 10^{-2}$.

### Case 1

We choose the velocity perturbation equal to $\delta^2 = \lambda^2$ (according to Definition 1) so that the initial condition reads:

$$ \rho(0, x_1) = \bar{\rho}(x_1), \quad u_1(0, x_1) = \bar{u}_1(x_1) + \delta^2 \cos(2K\pi x_1) \quad \phi(0, x_1) = \bar{\phi}(x_1), \qquad (78) $$

where $K = 16$, $\bar{\rho}$, $\bar{u}$, $\bar{\phi}$ are given by (77). Figure 6 shows the density perturbation ($|\rho(t, \cdot) - 1|$), the divergence of the velocity field $|\partial_{x_1} u_1(t, x_1)|$ and the absolute value of the electric potential ($|\phi(t, \cdot)|$) in semi-log scale and at the final time $t = 10^{-1}$ for the two type-A (DP2A1 and DP2A2), the type-CK (ARS) variants of penalized-IMEX schemes and for the SI-IMEX scheme from [3] (denoted by LSDIRK). The computation was carried out on a resolved mesh with $N = 10^4$ points corresponding to $\Delta x = 1/N$ with a CFL = 0.45 (CFL = 0.1 for ARS) and for the scaled Debye length $\lambda = 10^{-4}$. The long time behaviour of all the penalized-IMEX schemes on a fine mesh

<figure><img src="image_2.png" /><figcaption>Figure 6: Case 1, resolved mesh ($N = 10^4$) with type-A penalized-IMEX and SI-IMEX schemes. $\lambda = 10^{-4}$. Left: perturbed density $|\rho(t, x_1) - 1|$. Center: divergence of the velocity $|\partial_{x_1} u(t, x_1)|$. Right: electric potential $|\phi(t, x_1)|$, at time $t = 10^{-1}$.</figcaption></figure>

in Figure 6 shows the ability of the penalized-IMEX schemes to maintain the quasineutrality for large times. In Figure 7, the same quantities as of Figure 6 are shown with a under resolved mesh: $N = 10^2$ corresponding to a size mesh $\Delta x = 1/N$, the CFL condition remaining the same of the previous situation, as well as the scaled Debye length value $\lambda = 10^{-4}$. The results show that on this chosen coarser mesh, all the tested schemes maintain the solution over the quasi-equilibrium manifold. One can also see DP2A2 and LSDIRK project more quickly towards to the equilibrium compared to DP2A1. The gain in computational cost can be made clear by observing that a standard explicit scheme would need a time step of the order of $\lambda^2 = 10^{-8}$ while in the two tested situations we have employed respectively $\Delta t \sim 0.00225$ for the coarse mesh and $\Delta t \sim 0.000022$ for the fine mesh which lead to 46 and 4446 time steps, respectively. Moreover, one has to observe that with standard explicit numerical method the solution of the Poisson equation would be very difficult with such value of the Debye length, the resulting matrix being ill conditioned.

### Case 2

We now test our schemes with a larger perturbation since $\delta^2 = 10^{-2}$ in (78). The initial data is consequently not well-prepared in the sense of Definition 1. Figure 8 shows the density perturbation, the divergence of velocity and the absolute value of the potential for the type-A penalized-IMEX schemes (DP2A1 and DP2AA2) on a coarse mesh $N = 10^2$ (not resolving the small scale $\lambda = 10^{-4}$) with a CFL = 0.25 at time $t = 10^{-1}$ on the top
<figure><img src="image_3.png" /><figcaption>Figure 7: Case 1, under-resolved mesh ($N = 10^2$) with type-A penalized-IMEX and SI-IMEX schemes. $\lambda = 10^{-4}$. Left: perturbed density $|\rho(t, x_1) - 1|$. Center: divergence of the velocity $|\partial_{x_1} u(t, x_1)|$. Right: electric potential $|\phi(t, x_1)|$, at time $t = 10^{-1}$.</figcaption></figure>

and at a later time $t = 1$ on the bottom. Looking at the results, we observe that DP2A2 scheme is able to achieve the equilibrium state earlier, whereas DP2A1 scheme takes longer to capture the equilibrium. This is due to the stronger damping property of DP2A2. The results presented corroborate the conclusion of the analysis presented

<figure><img src="image_4.png" /><figcaption>Figure 8: Case 2, under-resolved mesh $N = 10^2$ with type-A penalized-IMEX schemes. $\lambda = 10^{-4}$. Left: perturbed density $|\rho(t, x_1) - 1|$. Center: divergence of the velocity $|\partial_{x_1} u(t, x_1)|$. Right: electric potential $|\phi(t, x_1)|$. Top time $t = 10^{-1}$. Bottom $t = 1$.</figcaption></figure>

via Theorem 3: even for initial data far from the quasineutrality the type-A penalized-IMEX schemes bring the solution close to the asymptotic solution.

Figure 9 shows the results of the same test in terms of density, divergence of the velocity field and electric potential as that of Figure 8 but for the type-CK penalized-IMEX scheme ARS at time $t = 0.025$. We can clearly see, as a confirmation of the theoretical results shown in the previous sections, that the type-CK scheme is not able to maintain the solution close to the asymptotic regime i.e. $\rho - 1 = O(\lambda^2)$, $\nabla \cdot u = O(\lambda^2)$ and $\phi = O(1)$. In particular, the solution starts blowing up in a very small time.

## 6.3 Perturbation of a quasineutral plasma in 1D with $\rho = 1, u = 0$

In [3] an IMEX Runge-Kutta scheme for the Euler-Poisson system in the quasineutral regime has been presented, and it was shown that the schemes were plagued with stability issues when subjected to a coarse mesh without the time step restriction of $\Delta t < \lambda$ (let us remark that the SI-IMEX scheme does not have this issue). Thus, this test case is aimed at showing the superiority of the new high order penalized-IMEX schemes here studied over the high order schemes presented in [3]. The same Butcher tableaux i.e. DP2A2 for both the scheme (denoted
<figure><img src="image_5.png" /><figcaption>Figure 9: Case 2, under-resolved mesh $N = 10^2$ with type-CK penalized-IMEX Schemes. $\lambda = 10^{-4}$. Left: perturbed density $|\rho(t, x_1) - 1|$. Center: divergence of the velocity $|\partial_{x_1}u(t, x_1)|$. Right: electric potential $|\phi(t, x_1)|$ at time $t = 0.025$.</figcaption></figure>

by CL-DP2A2) of [3] and the type-A DP2A2) penalized-IMEX scheme (denoted by AP-DP2A2) discussed in this work is used as a comparison in what follows.

The test case is from [38], wherein a small perturbation of amplitude $\delta = 10^{-2}$ is added to the density so that the initial condition is

$$ \rho(0, x_1) = 1.0 + \delta \sin(K\pi x_1), \quad u(0, x_1) = 0, \qquad (79) $$

with $K = 2220$ and $x_1 \in [0, 1]$. The value of $K$ is chosen such that $K \sim \lambda^{-1}$ to ensure that the wavelength of the density perturbation is of the same order as the Debye length. This is a slight perturbation of the steady state equilibrium, so initially the plasma is at rest. Again for all the results related to this test case, a semi-log

<figure><img src="image_6.png" /><figcaption>Figure 10: Perturbation of the density with the type-A penalized-IMEX (AP-DP2A2) and with the classical scheme (CL-DP2A2 from [3]). $\lambda = 10^{-4}$ and $N = 10^2$ (under-resolved mesh). Left: perturbed density $|\rho(t, x_1) - 1|$. Right: electric potential $|\phi(t, x_1)|$. Time = 0.035</figcaption></figure>

scale is used. Figure 10 shows the results for the density perturbation ($|\rho - 1|$) and for the absolute value of electric potential computed with the type-A (DP2A2) penalized-IMEX (denoted by AP-DP2A2) on the left and with the (denoted by CL-DP2A2) IMEX scheme on the right. The number of points in space is $N = 10^2$ such that the small physical space dynamics identified by $\lambda$ is not resolved. For the scheme (CL-DP2A2) the time step is $\Delta t = 10^{-4} = \lambda$ which is required for stability and which resolves the time scale identified by $\lambda$. Whereas for the type-A penalized-IMEX scheme (AP-DP2A2) the time step are chosen as $\Delta t = 0.5\Delta x$, not resolving the time scale identified by $\lambda$. Clearly, even with much smaller time steps the scheme CL-DP2A2 does not perform well (the potential starts blowing up). On the other side, even with a much larger time step, the new penalized-IMEX scheme is able to maintain the equilibrium state quite accurately ($|\rho - 1|$ is about $\lambda^2 = 10^{-8}$). Finally, Figure 11 shows the same physical quantities (perturbed density and potential) but just for the penalized-IMEX scheme and at a later time $t = 10^{-1}$ where the other class of schemes blows up. We can clearly see that the numerical solution is now very close to the equilibrium state.
<figure><img src="image_7.png" /><figcaption>Figure 11: Perturbation of the density with the type-A penalized-AP scheme (AP-DP2A2). $\lambda = 10^{-4}$, $N = 10^2$ (under-resolved mesh). Left: perturbed density $|\rho(t, x_1) - 1|$. Right: electric potential $|\phi(t, x_1)|$. Time = $10^{-1}$.</figcaption></figure>

## 6.4 Perturbation of a quasineutral plasma in 2D with $\rho = 1, u = (1, 1)$

We finally perform a two-dimensional test case. This is the extension of the test considered in Section 6.2. The following initial conditions are chosen

$$ 
\begin{align} 
\rho(0, x_1, x_2) &= 1, \nonumber \\ 
u_1(0, x_1, x_2) &= 1 + \sin(K\pi(x_1 - x_2)) + \lambda \sin(K\pi(x_1 + x_2)), \tag{80} \\ 
u_2(0, x_1, x_2) &= 1 + \sin(K\pi(x_1 - x_2)) + \lambda \cos(K\pi(x_1 + x_2)), \nonumber \\ 
\phi(0, x_1, x_2) &= 0, \nonumber 
\end{align} 
$$

where $(x_1, x_2) \in [0, 1]^2$, $K = 16$ and the perturbation is considered on a (divergence free) velocity field $u = (u_1, u_2)$. The perturbation chosen here pertains to a set of non well-prepared initial data (according to Definition 1). The CFL number is set to be 0.45 ($N = ??$ in each space direction) while periodic boundary conditions are imposed for all the components and a mean-free condition on $\phi$ is added to ensure its uniqueness. The value of the isentropic parameter $\beta$ for the pressure is taken to be 2: $p(\rho) = \rho^2$ (whereas $\beta = 1$ in the previous tests). Figure 12 shows

<figure><img src="image_8.png" /><figcaption>Figure 12: 2D test: Small perturbation of a quasineutral state with the type-A penalized-IMEX (DP2A2) scheme. Left: $t \to ||\rho(t, \cdot) - 1||_{L^2}$. Right: $t \to ||\nabla \cdot u(t, \cdot)||_{L^2}$. $\lambda = 10^{-2}, 10^{-3}$ and $10^{-4}$.</figcaption></figure>

the time evolution of the quantities $||\rho(t, \cdot) - 1||_{L^2}$ and $||\nabla \cdot u(t, \cdot)||_{L^2}$ in semi-log scale, for the time domain $t \in [0, 0.5]$ and $\lambda$ in the range $\{10^{-2}, 10^{-3}, 10^{-4}\}$. For each of the values of $\lambda$ we see that the type-A (DP2A2) penalized-IMEX scheme is able to maintain the initial proximity of the solution to the quasineutral limit.
# 7 Conclusions

In this paper, we developed and analyzed a new class of high order penalized-IMplicit-Explicit (IMEX) schemes for the Euler-Poisson system in the quasineutral scaling. These methods are designed to be uniformly stable with respect to the Debye length and to degenerate into accurate discretizations of the quasineutral system as the Debye length approaches zero. This feature enables the schemes to efficiently manage the coexistence of quasineutral and non-quasineutral regions, a scenario commonly encountered in practical plasma applications.

One of the key achievements of this work is addressing a major limitation of previous studies, which required well-prepared initial conditions to ensure consistency with the asymptotic regime. In contrast, the schemes proposed here are capable of handling more general initial conditions without losing stability or consistency with the asymptotic limit. This significantly broadens the applicability of these numerical methods to more realistic plasma scenarios where such assumptions cannot always be guaranteed.

Our analysis and numerical experiments demonstrate that the proposed schemes exhibit the desired strong Asymptotic-Preserving (AP) property as well as the asymptotic accuracy. Additionally, we investigated the performance of different IMEX schemes, revealing that certain classes of IMEX schemes, namely type-CK IMEX Runge-Kutta, may fail consistency and lose stability under specific conditions when close to the quasineutral limit.

In conclusion, the penalized-IMEX schemes introduced in this paper provide a robust and computationally efficient approach for simulating plasma phenomena in the quasineutral limit, making them a valuable tool for various applications. Future work could extend these schemes to more general models, incorporating departures from thermodynamic equilibrium and addressing the challenges related to the different masses of the particles.

# References

[1] ALVAREZ LAGUNA, A., PICHARD, T., MAGIN, T., CHABERT, P., BOURDON, A., AND MASSOT, M. An asymptotic preserving well-balanced scheme for the isothermal fluid equations in low-temperature plasmas at low-pressure. *Journal of Computational Physics* 419 (2020), 109634.

[2] ANANDAN, A., BOUTIN, B., AND CROUSEILLES, N. High order asymptotic preserving scheme for linear kinetic equations with diffusive scaling. *M2AN* (2024).

[3] ARUN, K. R., CROUSEILLES, N., AND SAMANTARAY, S. High order asymptotic preserving and classical semi-implicit RK schemes for the Euler-Poisson system in the quasineutral limit. *J. Sci. Comput.* 100, 1 (2024), Paper No. 24, 34.

[4] ARUN, K. R., AND SAMANTARAY, S. Asymptotic preserving low Mach number accurate IMEX finite volume schemes for the isentropic Euler equations. *J. Sci. Comput.* 82, 2 (2020), Art. 35, 32.

[5] BELAOUAR, R., CROUSEILLES, N., DEGOND, P., AND SONNENDRÜCKER, E. An asymptotically stable semi-Lagrangian scheme in the quasi-neutral limit. *J. Sci. Comput.* 41, 3 (2009), 341–365.

[6] BESSEMOLIN-CHATARD, M., CHAINAIS-HILLAIRET, C., AND VIGNAL, M.-H. Study of a finite volume scheme for the drift-diffusion system. Asymptotic behavior in the quasi-neutral limit. *SIAM J. Numer. Anal.* 52, 4 (2014), 1666–1691.

[7] BOSCARINO, S., FILBET, F., AND RUSSO, G. High order semi-implicit schemes for time dependent partial differential equations. *J. Sci. Comput.* 68, 3 (2016), 975–1001.

[8] BOSCARINO, S., PARESCHI, L., AND RUSSO, G. Implicit-Explicit Runge-Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit. *SIAM J. Sci. Comput.* 35, 1 (2013), A22–A51.

[9] BOSCARINO, S., QIU, J.-M., RUSSO, G., AND XIONG, T. A high order semi-implicit IMEX WENO scheme for the all-Mach isentropic Euler system. *J. Comput. Phys.* 392 (2019), 594–618.

[10] CHEN, F. F. *Introduction to plasma physics*. New York: Plenum Press, 1974.

[11] CORDIER, S., AND GRENIER, E. Quasineutral limit of an Euler-Poisson system arising from plasma physics. *Comm. Partial Differential Equations* (2000), 1099–1113.

[12] CRESTETTO, A., DELUZET, F., AND NEGULESCU, C. A hybrid method for anisotropic elliptic problems based on the coupling of an asymptotic-preserving method with the asymptotic limit model. *SIAM J. Sci. Comput.* 38, 3 (2016), A1821–A1847.

[13] CRISPEL, P., DEGOND, P., AND VIGNAL, M.-H. An asymptotically stable discretization for the Euler-Poisson system in the quasi-neutral limit. *C. R. Math. Acad. Sci. Paris* 341, 5 (2005), 323–328.

[14] CRISPEL, P., DEGOND, P., AND VIGNAL, M.-H. Quasi-neutral fluid models for current-carrying plasmas. *J. Comput. Phys.* 205, 2 (2005), 408–438.

[15] CRISPEL, P., DEGOND, P., AND VIGNAL, M.-H. An asymptotic preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit. *J. Comput. Phys.* 223, 1 (2007), 208–234.
[16] CROUSEILLES, N., DIMARCO, G., AND VIGNAL, M.-H. Multiscale schemes for the BGK-Vlasov-Poisson system in the quasi-neutral and fluid limits. Stability analysis and first order schemes. *Multiscale Model. Simul.* 14, 1 (2016), 65–95.

[17] DEGOND, P. Asymptotic-preserving schemes for fluid models of plasmas. In *Numerical models for fusion*, volume 39/40 of Panoramas Synthèses (2013), 1–90.

[18] DEGOND, P., DELUZET, F., AND DOYEN, D. Asymptotic-preserving particle-in-cell methods for the Vlasov-Maxwell system in the quasi-neutral limit. *J. Comput. Phys.* 330 (2017), 467–492.

[19] DEGOND, P., LIU, J.-G., AND VIGNAL, M.-H. Analysis of an asymptotic preserving scheme for the Euler-Poisson system in the quasineutral limit. *SIAM J. Numer. Anal.* 46, 3 (2008), 1298–1322.

[20] DIMARCO, G., LOUBÈRE, R., AND VIGNAL, M.-H. Study of a new asymptotic preserving scheme for the Euler system in the low Mach number limit. *SIAM J. Sci. Comput.* 39, 5 (2017), A2099–A2128.

[21] DIMARCO, G., MIEUSSENS, L., AND RISPOLI, V. An asymptotic preserving automatic domain decomposition method for the Vlasov-Poisson-BGK system with applications to plasmas. *J. Comput. Phys.* 274 (2014), 122–139.

[22] DIMARCO, G., AND PARESCHI, L. Asymptotic Preserving Implicit-Explicit Runge-Kutta Methods for Non-linear Kinetic Equations. *SIAM J. Numer. Anal.* 51, 2 (2013), 1064–1087.

[23] DIMITS, A., BATEMAN, G., BEER, M., COHEN, B., DORLAND, W., HAMMETT, G., KIM, C., KINSEY, J., KOTSCHENREUTHER, M., KRITZ, A., LAO, L., MANDREKAS, J., NEVINS, W., PARKER, S., REDD, A., SHUMAKER, D., SYDORA, R., AND WEILAND, J. Comparisons and physics basis of tokamak transport models and turbulence simulations. *Physics of Plasmas* 7, 3 (2000), 969–983.

[24] FABRE, S. Stability analysis of the Euler-Poisson equations. *J. Comput. Phys.* 101, 2 (1992), 445–451.

[25] FREIDBERG, J. P. *Plasma Physics and Fusion Energy*. Cambridge University Press, 2007.

[26] FRIDMAN, A. *Plasma Chemistry in Energy Systems and Environmental Control*. Cambridge University Press, 2008, p. 755–847.

[27] GÉRARD-VARET, D., HAN-KWAN, D., AND ROUSSET, F. Quasineutral limit of the Euler-Poisson system for ions in a domain with boundaries. *Indiana Univ. Math. J.* 62, 2 (2013), 359–402.

[28] GUERMOND, J.-L., AND POPOV, B. Fast estimation from above of the maximum wave speed in the Riemann problem for the Euler equations. *J. Comput. Phys.* 321 (2016), 908–926.

[29] HAIRER, E., AND WANNER, G. Solving ordinary differential equations. II, second ed., vol. 14 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 1996. Stiff and differential-algebraic problems.

[30] JANG, J., LI, F., QIU, J., AND XIONG, T. High order asymptotic preserving dg-imex schemes for discrete-velocity kinetic equations in a diffusive scaling. *J. of Comput. Phys.* 281 (2015), 199–224.

[31] JI, L., YANG, Z., LI, Z., WU, D., JIN, S., AND XU, Z. An asymptotic-preserving and energy-conserving particle-in-cell method for Vlasov-Maxwell equations. *J. Math. Phys.* 64, 6 (2023), Paper No. 063503, 19.

[32] JIN, S. Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations. *SIAM J. Sci. Comput.* 21 (1999), 441–454.

[33] JIN, S. Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review. *Riv. Mat. Univ. Parma* (2012), 177–216.

[34] KENNEDY, C. A., AND CARPENTER, M. H. Additive runge-kutta schemes for convection-diffusion-reaction equations. *Applied Numerical Mathematics* 44, 1 (2003), 139–181.

[35] LABRUNIE, S., CARRILLO, J. A., AND BERTRAND, P. Numerical study on hydrodynamic and quasi-neutral approximations for collisionless two-species plasmas. *J. Comput. Phys.* 200, 1 (2004), 267–298.

[36] MA, C., AND JIN, S. An implicit, asymptotic-preserving and energy-charge-conserving method for the Vlasov-Maxwell system near quasi-neutrality. *Commun. Comput. Phys.* 35, 3 (2024), 724–760.

[37] MAZOUFFRE, S. Electric propulsion for satellites and spacecraft: established technologies and novel approaches. *Plasma Sources Science Technology* 25, 3 (June 2016).

[38] NEGULESCU, C. Asymptotic-preserving schemes. Modeling, simulation and mathematical analysis of magnetically confined plasmas. *Riv. Math. Univ. Parma* (N.S.) 4, 2 (2013), 265–343.

[39] PARESCHI, L., AND RUSSO, G. Implicit-explicit Runge-Kutta schemes for stiff systems of differential equations. In *Recent trends in numerical analysis*, vol. 3 of *Adv. Theory Comput. Math*. Nova Sci. Publ., Huntington, NY, 2001, pp. 269–288.

[40] SLEMROD, M., AND STERNBERG, N. Quasi-neutral limit for Euler-Poisson system. *J. Nonlinear Sci.* 11, 3 (2001), 193–209.
[41] U. M. ASCHER, S. J. R., AND WETTON, B. T. R. Implicit-explicit methods for time-dependent partial differential equations. SIAM J. Numer. Anal. (1995), 797-823.

[42] U.M. ASCHER, S.J. RUUTH, R. S. Implicit-explicit runge-kutta methods for time-dependent partial differential equations. Appl. Numer. Math. (1997), 151-167.

[43] VLADIMIROV, V. V., VOLKOV, A. F., AND MEILIKHOV, E. Z. Semiconductor plasmas. Moscow Atomizdat (1979).