# High order asymptotic preserving and classical semi-implicit RK schemes for the euler-poisson system in the quasineutral limit

K.R. Arun, Nicolas Crouseilles, Saurav Samataray

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

## To cite this version:

K.R. Arun, Nicolas Crouseilles, Saurav Samataray. High order asymptotic preserving and classical semi-implicit RK schemes for the euler-poisson system in the quasineutral limit. 2022. hal-03783065

HAL Id: hal-03783065

https://hal.science/hal-03783065

Preprint submitted on 21 Sep 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.
HIGH ORDER ASYMPTOTIC PRESERVING AND CLASSICAL
SEMI-IMPLICIT RK SCHEMES FOR THE EULER-POISSON SYSTEM IN THE
QUASINEUTRAL LIMIT

K. R. ARUN, N. CROUSEILLES, AND S. SAMANTARAY

ABSTRACT. In this paper, the design and analysis of high order accurate IMEX finite volume schemes for the compressible Euler-Poisson (EP) equations in the quasineutral limit is presented. As the quasineutral limit is singular for the governing equations, the time discretisation is tantamount to achieving an accurate numerical method. To this end, the EP system is viewed as a differential algebraic equation system (DAEs) via the method of lines. As a consequence of this vantage point, high order linearly semi-implicit (SI) time discretisation are realised by employing a novel combination of the direct approach used for implicit discretisation of DAEs and, two different classes of IMEX-RK schemes: the additive and the multiplicative. For both the time discretisation strategies, in order to account for rapid plasma oscillations in quasineutral regimes, the nonlinear Euler fluxes are split into two different combinations of stiff and non-stiff components. The high order scheme resulting from the additive approach is designated as a classical scheme while the one generated by the multiplicative approach possesses the asymptotic preserving (AP) property. Time discretisations for the classical and the AP schemes are performed by standard IMEX-RK and SI-IMEX-RK methods, respectively so that the stiff terms are treated implicitly and the non-stiff ones explicitly. In order to discretise in space a Rusanov-type central flux is used for the non-stiff part, and simple central differencing for the stiff part. AP property is also established for the space-time fully-discrete scheme obtained using the multiplicative approach. Results of numerical experiments are presented, which confirm that the high order schemes based on the SI-IMEX-RK time discretisation achieve uniform second order convergence with respect to the Debye length and are AP in the quasineutral limit.

# 1. INTRODUCTION

Physics classifies the matter in the universe broadly into four different states: solids, liquids, gases and plasma. Even though the fourth state is what constitutes almost the entire universe, it is found in rarity in our proximity. At the same time matter in the form of plasma has far reaching scientific and industrial relevance. Mathematically, a plasma is primarily studied based on two different kinds of models: kinetic models and fluid models. Both these fields are indeed very active areas of research and with time new frontiers are being explored and established within them. Fluid models have the density, the temperature and the mean velocity as the base components under observation. Whereas kinetic models incorporate the time evolution of a distribution function which gives the probability of particles being in a given state in the six-dimensional phase space at a given time [5]. There are undoubtedly situations where kinetic models are inevitable for the application under consideration. But, at the same time there are many situations in which fluid models are adequate enough. Owing to their simplicity in dimensional requirements, fluids models are much easier to

<footnote>Date: September 21, 2022.</footnote>

<footnote n="2010">2010 Mathematics Subject Classification. Primary 35L45, 35L60, 35L65, 35L67; Secondary 65M06, 65M08.</footnote>

<footnote>Key words and phrases. Compressible Euler-Poisson system, Incompressible Euler system, Quasineutral limit, SI-IMEX-RK schemes, Asymptotic preserving, Finite volume method.</footnote>

<footnote>K. R. A. acknowledges Core Research Grant - CRG/2021/004078 from Science and Engineering Research Board, Department of Science & Technology, Government of India.</footnote>
numerically approximate and simulate. We refer the reader to [16, 38, 41, 52] for more details on the physics of plasma, their mathematical models and analysis.

In this paper we consider the one-fluid Euler-Poisson (EP) system [25, 26, 31] which is a macroscopic hydrodynamic model that simulates the dynamics of electrons in a plasma with constant positive ion background. It is a rather simple fluid model, albeit comprising of most of the numerical pathologies offered by more complex models. Therefore, it serves as a perfect test bed for the design and analysis of a new class of higher order methods. The one-fluid EP system of equations, is in itself a simplification of a more physically relevant plasma model, the two-fluid EP plasma model [16, 48]. The two-fluid EP system consists of equations for the conservation of masses and momentum for each of electrons and protons, which are coupled with a Poisson equation for the electrostatic potential. This constitutes a hyperbolic elliptic coupling. The literature on the analysis of the EP system leading to the existence and uniqueness of solution is quite vast; we refer the interested reader to [4, 13, 35, 54, 55, 56, 57, 58, 59] and the references cited therein, for an exposure of the same. For a detailed analytical study of the EP system, the reader is referred to, e.g. [21, 69]. The two fluid model gives rise to two important physical scales: the Debye length and the electron plasma period; see [16, 52] for details. The Debye length $\lambda_D$ is a measure of the typical length scale of charge imbalances in plasma, and the electron plasma period $\tau_p$ is the period of oscillations which occur (due to electrostatic restoring force) when such imbalances of charges take place. In the quasineutral region, the electric charges vanish and simultaneously, the electron plasma period reduces drastically. As a consequence, high frequency plasma oscillations are generated. Classical discretisation procedures fail to capture these oscillations unless their space and time mesh sizes resolve the Debye length and the electron plasma period. In other words, the space and time steps, $\Delta x$ and $\Delta t$, have to be severely restricted by

$$ (1.1) \qquad \Delta x \le \lambda_D, \quad \Delta t \le \tau_p. $$

Explicit discretisation methods where the discretisation parameters are constrained as in (1.1) are naturally computationally expensive and are more prone to developing instabilities.

There is an abundance of literature on methods which attempt to overcome the mesh size restrictions imposed by the presence of quasineutral regimes in the flow. These methods are developed for either kinetic models or fluid models. The literature in the case of kinetic models is quite vast; see [12, 17, 53, 60, 61, 62, 63, 70, 71, 76] and the references therein. However, the literature on fluid models is less plentifully rich. For some of the noted developments focusing on fluid models, we refer to [18, 31, 42, 72, 74]. However, in many practical problems of high relevance, one encounters the co-existence of quasineutral and non-quasineutral regions within the same problem. Therefore flow characteristics are multiscaled in nature. When subjected to appropriate scalings the multiscale features of plasma flows is often characterised by one or more singular perturbation parameters in the governing equations, in a non-dimensional form [20, 30, 36, 37, 44, 47, 49, 77]. The quasineutral limit is a singular limit of the EP system. Owing to this multiscaled nature, numerical simulation of plasmas is still an eminent challenge for the scientific community. It plagues numerical methods not only with outrageous stability restrictions but also with inaccuracies in the asymptotic regimes, especially close to the singular limit. This pathology yearns for the need to address the numerical approximations of plasma models in the asymptotic preserving (AP) framework [25, 45, 46]. An AP scheme has two main characteristic features, exactly aiming to remedy the above noted challenges associated with the numerical approximation of plasma flows, which are:

* its stability requirements are independent of the multiscale or a singular parameter;

* in the limit of the singular parameter, the numerical scheme transforms itself to a scheme for the limit system.
It has been reported in the literature that the ability of a scheme to be AP gets primarily dictated by the choice of the time integrator, rather than the spatial approximation techniques. That being said, it doesn't mean that caution shouldn't be exercised while choosing the spatial discretisation. To this end, semi-implicit Runge-Kutta(RK) or Implicit Explicit (IMEX)-RK [67, 68] stiff ordinary differential equation (ODE) integrators have played a dominating role in dispensing the tasking duty of tackling the stiffness in the system, and the limiting singular nature of the governing equations in hand. Recently a new class of IMEX-RK schemes called semi-implicit (SI)-IMEX-RK schemes for stiff ODEs has been developed in [9]. The difference between the new class i.e. the SI-IMEX-RK schemes and the earlier (additively partitioned) IMEX-RK schemes [67, 68] is that in case of the SI-IMEX-RK schemes there is no need to identify an additive stiff and non-stiff part in the ODE. The stiff variable can be determined implicitly. This gives a massive advantage in systems where the stiff variables cannot be optimally split into summable stiff-non-stiff parts. Moreover, in many cases they aid in obtaining methods which are only linearly implicit.

AP schemes based on IMEX time discretisation methods have found widespread relevance in a wide array of problems with singular parameters, such as: the low Mach number Euler equation [1, 6, 7, 10, 19, 27, 28, 39, 50, 51, 64, 66, 73, 75], to name but a few. First order IMEX schemes for plasma flow problems have been developed and studied in the literature, we refer the reader to [22, 23, 25, 26, 65] for a few developments in this direction. However, the methods proposed in these references are all first order accurate in time and space. In this paper we are interested to develop a general platform to design high-order IMEX schemes for plasma dynamics problems using the one-fluid EP model as a prototype system. In [29] authors have developed higher order semi-implicit schemes based on additively partitioned IMEX-RK schemes for hyperbolic kinetic models of plasma flows. Recently, SI-IMEX-RK based high order schemes for hyperbolic-elliptic kinetic models of plasma have been developed and analysed in [32, 33, 34]. In this paper our endeavour is to lay out a path way which facilitates the design of high order semi-implicit schemes for hyperbolic-elliptic coupled systems. To this end, via method of lines we perceive the coupled system as an index-1 differential algebraic system (DAEs). The time discretisation is designed by combining the direct approach [40] to design implicit RK schemes for DAEs with both the additively partitioned IMEX-RK schemes and the SI-IMEX schemes. Therefore this paper has also an underlying focus beyond the development of high-order schemes for the EP system i.e. this also showcases an approach to design high-order semi-implicit schemes for stiff DAEs as was done in [67, 68] for stiff ODEs, which is a novelty in itself.

The two main challenges in getting higher order accuracy are the following. First, the right hand side of the momentum equation contains a term which is non-linear and only one of the components in this non-linear product survives in the quasineutral limit. So, any AP discretisation should treat that component implicitly. In order to reduce the implicitness, the other component in the non-linear product should be explicit. In other words, we have a product of explicit and implicit terms, necessitating the use of SI-IMEX schemes instead of IMEX. Second, the Poisson equation has to be enforced at every stage of the RK time stepping. Since the Poisson equation is time independent, it has to be treated as an algebraic relation as far as the time discretisation is concerned. Formal asymptotic analysis of the schemes developed on both the additively partitioned IMEX-RK and the SI-IMEX-RK platforms reveals that the former is not asymptotically consistent, but the later is. This is in agreement with the analysis presented in [25] regarding the first-order classical and AP schemes, of which the schemes designed in this paper turn out to be higher order extensions.

The rest of this paper is organised as follows. In Section 2, we present the governing equations which consider the movement of the electrons within a constant non evolving positively charged ion background, leading to the so-called one-fluid EP system [25, 26]. Furthermore results pertaining to asymptotic analysis of the one-fluid model from the literature [25], which aid in the analysis of
the numerical schemes, are also presented in this section. Particularly, equivalent reformulations of the EP system which aid in obtaining the AP property of the schemes are furnished. In Section 3, a linearised version of the one dimensional EP system is derived. Upon Fourier transformation in space the linearised system turns out to be an index-1 DAEs [40]. Motivated by this the EP system is cast into a system of index-1 DAEs via method of lines. A brief discussion on the direct approach to design implicit RK schemes for index-1 DAEs, additively partition IMEX-RK schemes [67, 68] and SI-IMEX-RK schemes [9], is also included in this section. In Section 4 we combine the direct approach for index-1 DAEs with the additively partitioned IMEX-RK schemes and SI-IMEX-RK schemes to obtain two classes of high-order semi-implicit schemes for the index-1 DAEs, in in-turn for the EP system. Section 5 is devoted to the asymptotic analysis of the developed high-order schemes. Therein the additively partitioned IMEX-RK schemes are shown to be not consistent with the asymptotic limit system, whereas the SI-IMEX-RK scheme is shown to be AP being consistent with the incompressible limit system. In Section 6, the space discretisation technique is presented based on the finite volume framework in which a Godunov-type scheme for the explicit fluxes and a simple central differencing for the implicit gradient terms are used. The fully discrete scheme is also proved to be asymptotically consistent. In Section 7 some numerical case studies are carried out to test and display the high-order accuracy, AP property and the superior performance of the SI-IMEX-RK schemes to the IMEX-RK schemes. Finally, the paper is closed in Section 8 with some conclusions and future plans for further extension of the schemes to hyperbolic elliptic systems.

## 2. ONE-FLUID EULER-POISSON SYSTEM AND ITS QUASINEUTRAL LIMIT

The one-fluid EP system is a hydrodynamic model of plasma electrons under the action of an electrostatic force. The governing equations of this model are the compressible Euler equations for the conservation of mass and momentum, coupled with a Poisson equation for the electric potential. The scaled non-dimensionalised system of equations reads:

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

$$ (2.2) \qquad \partial_t q + \nabla \cdot \left( \frac{q \otimes q}{\rho} \right) + \nabla p(\rho) = \rho \nabla \phi, $$

$$ (2.3) \qquad \varepsilon^2 \Delta \phi = \rho - 1. $$

Here, the independent variables $t \ge 0$ and $x \in \mathbb{R}^d$ are, respectively, the time and space, and the dependent variables $\rho(t, x) \ge 0$, $q(t, x) = \rho u(t, x) \in \mathbb{R}^d$ and $\phi(t, x) \in \mathbb{R}$ are, respectively, the electron density, electron momentum and electric potential. The parameter $\varepsilon := \lambda_D/x_{\text{ref}}$ is a scaled Debye length; cf. [16, 52]. The pressure law in the momentum equation (2.2) is taken to be isentropic, i.e. $p(\rho) = \rho^\gamma$, with $\gamma > 1$. Finally, the spatial operators $\nabla, \nabla \cdot$ and $\Delta$ are respectively, the gradient, divergence and the Laplacian. For a complete derivation of the above system; see [25].

The one-fluid EP model considered here only takes into account the negatively charged electrons with the non-dimensionalised scaled negative charge being equal to $-1$. In the Poisson equation (2.3), the factor 1 on the right hand side term $\rho - 1$ signifies that the plasma consists of a uniform ion background with a constant ion density 1 in terms of the scaled dimensionless variables. The microscopic quantity $\varepsilon$ plays the role of a singular perturbation parameter; see also [26, 31] for more details.

### 2.1. Quasineutral Limit.
The rest of this section is devoted to a brief summary of some analytical results pertaining to the quasineutral limit of the EP system (2.1)-(2.3), deemed relevant for the analysis to be carried out later, along the lines of [25]. As discussed in the introduction, the multi-scale nature of the EP system is elucidated through the parameter $\varepsilon$, when there is a disparity in
the scale of the physical Debye length compared to the macroscopic length scale. The quasineutral limit model is obtained by letting $\epsilon \to 0$ in the system (2.1)-(2.3). Further, when $\epsilon \ll 1$, the flow regime can be considered to be close to the quasineutral limit.

Let us assume that the dependent variables $\rho, q$ and $\phi$ yield the following limits as $\epsilon \to 0$:

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

Passing to the limit $\epsilon \to 0$ in (2.1)-(2.3), cf. [25], we obtain the ensuing incompressible Euler system as the quasineutral limit of the Euler-Poisson system;:

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

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

$$ (2.7) \qquad \rho_{(0)} = 1. $$

Based on the above result we put together the following definition, which also concerns itself with the choice of the appropriate set of initial conditions for the EP system (2.1)-(2.3) and, reveals the multi-scale nature of the solution in accordance with the formal limits considered in (2.4).

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

$$ (2.8) \qquad \rho = \rho_{(0)} + \epsilon^2 \rho_{(2)}, \quad q = q_{(0)} + \epsilon^2 q_{(2)} $$

$$ (2.9) \qquad \rho_{(0)} = 1, \quad \nabla \cdot q_{(0)} = 0. $$

In the incompressible Euler limit system (2.5)-(2.7), the hydrodynamic pressure is identified as the negative of the limiting electric potential $\phi_{(0)}$. It is a mixed elliptic-hyperbolic system for the unknowns $u_{(0)}$ and $\phi_{(0)}$. The electric potential $\phi_{(0)}$ plays the role a Lagrange multiplier for the incompressibility constraint (2.5). Taking divergence of the momentum equation (2.6), and using the divergence-free condition (2.5) on the velocity $u_{(0)}$, see also [15]; yields the following elliptic equation for $\phi_{(0)}$:

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

The equation (2.10) can be used to eliminate the divergence constraint (2.5). In fact, the divergence-free condition (2.5) needs to be satisfied only at time $t = 0$. We summarise a few results from [25] as follows.

**Proposition 2.2.** The following reformulated incompressible Euler system:

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

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

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

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

is equivalent to the incompressible Euler system (2.5)-(2.7) for smooth solutions; in the sense that the solution of one system satisfies the other system and vice versa.

**Remark 2.3.** The divergence-free condition (2.14) on the initial velocity field in the reformulated Euler-Poisson system plays a vital role in obtaining the incompressibility condition for the limiting Euler system (2.5)-(2.7).

**Remark 2.4.** Even though the reformulated system doesn't serve much of a purpose in the analysis of solutions of the Euler-Poisson system or the incompressible Euler system, it finds its relevance in the analysis of the numerical schemes to be presented in later sections.
As is shown in the case of the limiting incompressible system (2.11)-(2.14), similarly, the scaled Euler-Poisson system (2.1)-(2.3) can be also shown to have a reformulated counterpart; see [25] for a formal derivation of the same.

**Proposition 2.5.** The reformulated Euler-Poisson system comprises of the mass conservation equation (2.1), the momentum conservation equation (2.2) with a reformulated poisson equation and two constraints on the initial conditions given as below:

$$ (2.15) \qquad \nabla \cdot ((\varepsilon^2 \partial_t^2 + \rho)\nabla\phi) = \nabla^2: \left( \frac{q \otimes q}{\rho} + pI \right), $$

$$ (2.16) \qquad \varepsilon^2 \Delta\phi = \rho - 1, \quad \text{at } t = 0, $$

$$ (2.17) \qquad \varepsilon^2 \Delta \partial_t \phi = \nabla \cdot q, \quad \text{at } t = 0. $$

The reformulated system is equivalent to the Euler-Poisson system (2.1)-(2.3) as long as their solutions are smooth.

**Remark 2.6.** As in the case of the reformulated incompressible system with respect to the divergence condition on the velocity, here also for the compressible system the reformulation has to satisfy the Poisson equation and its time derivative, only initially.

### 3. DAE FORMULATION AND RK SCHEMES

The challenges faced in the numerical approximation of the EP system (2.1)-(2.3) are multifold. First, the non-evolutionary Poisson equation (2.3) has to be discretised along with the evolutionary hyperbolic conservation equations (2.1)-(2.2). Second, in order to address the stiffness arising in the quasineutral regime, some implicitness is to be introduced whilst maintaining the maximal possible explicitness to achieve computational optimality. Third, the discretisation should allow high order extensions in space and time. Fourth, and most importantly, the resulting scheme should be AP in the quasineutral limit.

In this section, we address the first challenge mentioned above by drawing an analogy between the EP system (2.1)-(2.3) and a class of index-1 DAEs. We refer the reader to, e.g. [14] for a comprehensive treatment of the theory and numerics of DAEs. Subsequently, a brief summary of the so-called indirect approach [40] to derive implicit RK schemes for DAEs is put on display. Drawing inspiration from this approach, and to tackle the second concern, two different IMEX-RK techniques for stiff systems of ODEs are succinctly overviewed in order to design semi-implicit schemes for the EP system, possessing the required amount of implicitness.

#### 3.1. Linear Subsystem and Index-1 DAEs.
To elucidate the parallel between the EP system and index-1 DAEs, we first linearise the EP system (2.1)-(2.3), following the lines of [26], about a stationary homogeneous state $\bar{\rho} = 1$, $\bar{q} = 0$, $\partial_x \bar{\phi} = 0$. Let the perturbations from the considered stationary state be represented by an infinitesimal $\delta \ll 1$. An asymptotic expansion of the unknowns around the constant state with respect to $\delta$ is given by;

$$ (3.1) \qquad \rho = 1 + \delta\rho' + \mathcal{O}(\delta^2), \quad q = \delta u' + \mathcal{O}(\delta^2), \quad \phi = \delta\phi' + \mathcal{O}(\delta^2). $$

Considering only the first order terms in $\delta$, and dropping the primes for the sake of convenience, we obtain the following linearised EP system:

$$ (3.2) \qquad \partial_t \rho + \partial_x u = 0, $$

$$ (3.3) \qquad \partial_t u + c_s^2 \partial_x \rho = \partial_x \phi, $$

$$ (3.4) \qquad \varepsilon^2 \partial_{xx} \phi = \rho, $$
where $c_s = \sqrt{p'(1)}$ is the speed of sound. We take the partial Fourier transform, in space of the equations (3.2)-(3.4) to obtain the following system of equations:

$$ (3.5) \qquad \partial_t \hat{\rho} + i\xi \hat{u} = 0, $$

$$ (3.6) \qquad \partial_t \hat{u} + i\xi c_s^2 \hat{\rho} = i\xi \hat{\phi}, $$

$$ (3.7) \qquad -\varepsilon^2 \xi^2 \hat{\phi} = \hat{\rho}. $$

Here, $\hat{\rho}, \hat{u}$ and $\hat{\phi}$ are the Fourier transformed variables corresponding to $\rho, u$, and $\phi$, respectively. Furthermore, to comprehend the analogy between the Fourier transformed linearised EP system (3.5)-(3.7), the Fourier space variable $\xi$ is treated as a fixed parameter.

**Proposition 3.1.** The set of equations (3.5)-(3.7) constitutes an index-1 DAE system for the unknown vector $V = (\hat{\rho}, \hat{u}, \hat{\phi})$ with respect to the time variable $t$.

*Proof.* To begin the proof, we first make the observation that (3.7) is a purely algebraic equation free of any differential terms. In order to substantiate the claim in the proposition, i.e. the equations (3.5)-(3.7) form a DAE system, let us consider the functions $F, F_1, F_2$ and $F_3$, defined by

$$ (3.8) \qquad \mathbb{F}(t, V, V') := \begin{pmatrix} F_1(t, V, V') \\ F_2(t, V, V') \\ F_3(t, V, V') \end{pmatrix} := \begin{pmatrix} \frac{\partial \hat{\rho}}{\partial t} + i\xi \hat{u} \\ \frac{\partial \hat{u}}{\partial t} + i\xi c_s^2 \hat{\rho} - i\xi \hat{\phi} \\ -\varepsilon^2 \xi^2 \hat{\phi} - \hat{\rho} \end{pmatrix}. $$

Here we have denoted $V' = \partial_t V$. Evaluating the Jacobian $\frac{\partial \mathbb{F}}{\partial V'}$ yields

$$ (3.9) \qquad \frac{\partial \mathbb{F}}{\partial V'} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix}, $$

which is a non-invertible/singular matrix. Therefore, using the definition of DAEs, cf. [14], we conclude that the equations (3.5)-(3.7) constitute a system of DAEs.

Next, we address the claim that the DAE system (3.5)-(3.7) is of index-1. Towards this interest, we differentiate (3.7) with respect to $t$ and substitute for $\partial_t \hat{\rho}$ from (3.5), to get a differential equation for $\hat{\phi}$. At the end, the DAE system (3.5)-(3.7) gets transformed into

$$ (3.10) \qquad \partial_t \hat{\rho} + i\xi \hat{u} = 0, $$

$$ (3.11) \qquad \partial_t \hat{u} + i\xi c_s^2 \hat{\rho} - i\xi \hat{\phi} = 0, $$

$$ (3.12) \qquad -\varepsilon^2 \xi^2 \partial_t \hat{\phi} + i\xi \hat{u} = 0. $$

It can easily be verified that the set of equations (3.10)-(3.12) is purely differential, and not DAEs, by taking the Jacobian with respect to $V'$, as done above. Since only one differentiation was required to transform the DAEs (3.5)-(3.7) into the ODEs (3.10)-(3.12), we conclude that the index of the DAEs is one. $\square$

Motivated by the conclusion of the above proposition for the linearised EP system (3.2)-(3.4), we rewrite the non-linear one-fluid EP system (2.1)-(2.3) as a set of DAEs via the following definition.

**Definition 3.2.** The DAE formulation of the scaled EP system (2.1)-(2.3) is given by

$$ (3.13) \qquad \partial_t U + H(U, \phi) = 0, $$

$$ (3.14) \qquad G(U, \phi) = 0, $$
where

$$ (3.15) \qquad U := \begin{pmatrix} \rho \\ q \end{pmatrix}, \quad H(U, \phi) := F(U) - S(U, \phi), \quad G(U, \phi) := \varepsilon^2 \Delta\phi - \rho + 1, $$

and

$$ (3.16) \qquad F(U) := \begin{pmatrix} \nabla \cdot q \\ \nabla \cdot \left( \frac{q \otimes q}{\rho} \right) + \nabla p(\rho) \end{pmatrix}, \quad S(U, \phi) := \begin{pmatrix} 0 \\ \rho \nabla \phi \end{pmatrix}. $$

We underline the fact that the method of lines approach is used to write the system of PDEs (2.1)-(2.3) in a DAE form given by (3.13)-(3.14). Analogous treatment of PDEs using the method of lines philosophy to obtain a system of ODEs is a regular practice in the literature; see, e.g. [1, 3, 68] and the references therein for examples.

*Remark 3.3.* Note that unlike the yield typically obtained using the method of lines, i.e. a system of ODEs, in this paper, the semblance is drawn with a system of DAEs.

## 3.2. Indirect Approach for RK Methods for DAEs of Index 1.
Having recast the EP system as a set of DAEs in (3.13)-(3.14), we now turn towards the numerical approximation. It has to be noted that even though implicit RK schemes were primarily designed for stiff systems of ODEs, these methods can be systematically extended to even broader classes of equations, such as DAEs. What follows is a brief exposure to the so-called indirect approach [40] to design implicit RK schemes for index-1 DAEs. Dedicated to this aim, we consider an arbitrary set DAEs of index 1 in the form

$$ (3.17) \qquad y' = f(y, z), $$

$$ (3.18) \qquad 0 = g(y, z), $$

where $f$ and $g$ are sufficiently smooth functions, and $y$ and $z$ are vectors of appropriate dimensions with respect to the functions $f$ and $g$.

We now consider an $s$-stage implicit RK method defined by the triple $(A, c, \omega)$, where $A = (a_{k,l})$ are the coefficients, $c = (c_k)$ are the intermediate times and $\omega = (\omega_k)$ are the weights. The indirect approach [40] defines an RK scheme for (3.17)-(3.18) with intermediate stages $(Y^k, Z^k)$ in the following way:

$$ (3.19) \qquad Y^k = y^n + \Delta t \sum_{l=1}^{s} a_{k,l} f(Y^l, Z^l), \quad k = 1, 2, \dots, s, $$

$$ (3.20) \qquad 0 = g(Y^k, Z^k), \quad k = 1, 2, \dots, s, $$

where the last row of the coefficient matrix $A$ is taken to be equal to the vector $\omega$; in other words, the implicit RK scheme is stiffly accurate (SA); see [40]. Hence, the numerical solution $(y^{n+1}, z^{n+1})$ is defined by $y^{n+1} = Y^s$ and

$$ (3.21) \qquad g(y^{n+1}, z^{n+1}) = g(Y^s, Z^s) = 0. $$

Note that the indirect approach enables us to ensure that the manifold

$$ (3.22) \qquad 0 = g(y, z) $$

is preserved by all the intermediate stages and automatically by the update stage. For a detailed discussion on the derivation, validity and stability analysis of the direct approach of implicit RK schemes applied to index-1 DAEs, interested readers may refer to [40].

*Remark 3.4.* The manifold preserving property (3.21) of the indirect approach plays a crucial role in obtaining a consistent discretisation of non-evolutionary equations, such as the Poisson equation in the EP system (2.1)-(2.3).
The indirect approach presented in [40], and briefed above, in its core, takes implicit RK schemes for ODEs and derives implicit schemes for DAEs. Yes indeed an implicit scheme does fair well mathematically in the sense that it yields the asymptotic properties, such as stability and consistency which leads to the AP property. But, as a drawback, because of the fully implicit nature of the method, it becomes computationally very expensive when tackling problems, especially with nonlinearities, such as the EP system (3.13)-(3.14). In order to reduce the overdose of implicitness offered by fully implicit RK schemes via the indirect approach platform (3.19)-(3.21), in the following, we consider two semi-implicit RK approaches. The ultimate aim is to develop RK time steppings for index-1 DAEs and in turn for the EP system (3.13)-(3.14) as well, with an optimal portion of implicitness.

**3.3. IMEX-RK Schemes.** Semi-implicit or IMEX schemes trade unconditional stability of fully implicit schemes in order to draw some gains in the computational efficiency front. IMEX-RK schemes are a staple for the CFD community; see, e.g. [1, 3, 7, 8, 9, 43, 66, 68, 78] for IMEX schemes applied to tackle problems plagued with stiffness. Introduction of an adequate dosage of implicitness through IMEX schemes provides an additional amount of stability which is usually required for methods resolving singular parameters, while maintaining a pardonable level of computational performance. In this work we consider two types of IMEX-RK schemes, namely

* additively partitioned IMEX-RK schemes [2, 67];

* semi-implicit IMEX-RK schemes [9].

The former one, the class of additively partitioned IMEX-RK schemes, is primarily designed for stiff ODEs in which the identifiable stiff and non-stiff parts are summed together. To this end, consider the following system of additively partitioned stiff ODEs:

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

where $y = y(t) \in \mathbb{R}^n$, $f, g: \mathbb{R} \times \mathbb{R}^n \to \mathbb{R}^n$, and $\varepsilon > 0$ is a small parameter known as the stiffness parameter. Here, $f$ and $g$ are referred to as, respectively, the non-stiff and stiff parts.

**Definition 3.5.** Consider the stiff system (3.23) in additive form. Let $y^n$ be the numerical solution of (3.23) 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 following $s$ intermediate stages trailed by the update stage:

$$ (3.24) \qquad Y^k = y^n + \Delta t \sum_{\ell=1}^{k-1} \tilde{a}_{k,\ell} f(t^n + \tilde{c}_\ell \Delta t, Y^\ell) + \Delta t \frac{1}{\varepsilon} \sum_{\ell=1}^{s} a_{k,\ell} g(t^n + c_\ell \Delta t, Y^\ell), \quad k = 1, 2, \dots, s, $$

$$ (3.25) \qquad y^{n+1} = y^n + \Delta t \sum_{k=1}^{s} \tilde{\omega}_k f(t^n + \tilde{c}_k \Delta t, Y^k) + \Delta t \frac{1}{\varepsilon} \sum_{k=1}^{s} \omega_k g(t^n + c_k \Delta t, Y^k). $$

The IMEX-RK scheme (3.24)-(3.25) is characterised by the intermediate step sizes $\tilde{c} = (\tilde{c}_1, \tilde{c}_2, \dots, \tilde{c}_s)$ and $c = (c_1, c_2, \dots, c_s)$, the RK coefficients $\tilde{A} = (\tilde{a}_{i,j})$ and $A = (a_{i,j})$ and the weights $\tilde{\omega} = (\tilde{\omega}_1, \tilde{\omega}_2, \dots, \tilde{\omega}_s)$ and $\omega = (\omega_1, \omega_2, \dots, \omega_s)$. The matrices $\tilde{A}$ and $A$ are $s \times s$ matrices, where $\tilde{a}_{i,j} = 0$ for $j \ge i$ for the scheme to be explicit in $f$ and we further impose $a_{i,j} = 0$ for $j > i$, to get a diagonally Implicit RK (DIRK) scheme. The IMEX-RK schemes are usually represented in a compact form using a double Butcher tableaux

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

<caption>FIGURE 1. Double Butcher tableaux of an IMEX-RK scheme.</caption>
There exist a large section of stiff problems amongst ODEs wherein identifying an additive partition of the right hand side into stiff and non-stiff subparts may not be possible. There are even systems where it is possible to find an additive splitting of the stiff terms, however, stiffness can be associated with a nonlinear term. Consequently, in this scenario the implementation of additive IMEX-RK schemes may lead to a need for solving highly nonlinear equations, which in-turn requires the usage of expensive, iterative Newton solvers.

In the light of the aforestated predicaments the second type of schemes that we consider here, i.e. the semi-implicit IMEX-RK schemes, have been recently presented in [9]. The authors have devised an efficient platform to derive IMEX-RK schemes for problems where a multiplicative stiff-nonstiff splitting can be recognised. Further, the stiffness can be associated to some part of the unknown vector, and rest of it might be non-stiff. In such scenarios, these schemes have an advantage that either a linear solver can be employed or a Newton iteration for a problem with reduced non-linearity can be engaged; see [10, 11, 33] for some related works.

Consider a general stiff-ODE in the form

$$ (3.26) \qquad y' = H(y, y), $$

where in $H$ the first argument represents the non-stiff part of $y$ and the second argument the stiff-part. Note that the right hand side is assumed to have a stiff dependence only through the last argument. In order to distinguish these, we use subscripts and write

$$ (3.27) \qquad H(y, y) = H(y_E, y_I). $$

It has to be noted that both $y_E$ and $y_I$ represent the same variable $y$ itself. The subscripts are used to only distinguish between the stiff and non-stiff components. The class of problems characterised by the system (3.26)-(3.27) is said to be in a generalised partitioned form. The additively partitioned stiff-problems considered above also fall under this general class. In what follows, we iterate a semi-implicit time discretisation technique to numerically approximate the solution of (3.26)-(3.27) as presented in [9]. The discretisation is done in such a way that the variable $y$ without a stiff dependence is treated explicitly, and the other part implicitly. The formal definition of the scheme is given below, where we consider only a DIRK variant; the interested reader is referred to [9] for other variants and a detailed exposure to the scheme.

**Definition 3.6.** Suppose that the RK coefficients $\tilde{A} = (\tilde{a}_{k,l})$, $A = (a_{k,l})$ and the weights $\omega = (\omega_k)$ are given. A semi-implicit IMEX-RK (SI-IMEX-RK) scheme for (3.26) updates the numerical solution $y^n$ at time $t^n$ to the solution $y^{n+1}$ at time $t^{n+1} = t^n + \Delta t$ through the following $s$ intermediate steps:

$$ (3.28) \qquad Y_E^k = y^n + \Delta t \sum_{\ell=1}^{k-1} \tilde{a}_{k,\ell} H(Y_E^\ell, Y_I^\ell), \quad \text{for } k = 1, 2, \dots, s, $$

$$ (3.29) \qquad Y_I^k = y^n + \Delta t \sum_{l=1}^{k} a_{k,l} H(Y_E^l, Y_I^l), \quad \text{for } k = 1, 2, \dots, s, $$

and the final update step:

$$ (3.30) \qquad y^{n+1} = y^n + \Delta t \sum_{k=1}^{s} \omega_k H(Y_E^k, Y_I^k). $$

Here, the vectors $Y_I$ and $Y_E$ correspond to the implicit and explicit parts of the intermediate solutions, and the corresponding update formulae are implicit and explicit, respectively.
Remark 3.7. In order to facilitate the design of IMEX schemes for the EP system (2.1)-(2.3), based on the indirect approach discussed above, henceforth, we consider only SI-IMEX-RK schemes with stiffly accurate implicit solvers.

As a consequence of the above remark, the numerical solution of the SI-IMEX-RK schemes under consideration in this paper boils down to the solution of the last implicit intermediate stage s, i.e.,

$$ (3.31) \qquad y^{n+1} = Y_I^s $$

The order conditions of SI-IMEX schemes and additively partitioned IMEX-RK schemes are more or less the same [9]. Moreover the SI-IMEX schemes are also represented using a double Butcher tableaux, cf. Figure 1.

## 4. TIME SEMI-DISCRETE SCHEMES FOR THE EP SYSTEM

The goal of this section is to address the third and fourth challenges, outlined in Section 3, faced by numerical approximations of the EP system in the quasineutral limit: achieving both high orders of accuracy and the AP property. To this end, we develop two different classes of semi-implicit schemes for the EP system (2.1)-(2.3). Both the schemes are based on the indirect approach for DAEs as discussed in the preceding section. The main difference between the two schemes is that one is designed on the additively partitioned IMEX-RK platform and the other on the SI-IMEX-RK platform. The first scheme, despite being linearly implicit and hence computationally efficient, when subjected to an asymptotic analysis reveals that it is not consistent with the quasineutral limit. As a remedy to the paltry performance of the additively partitioned IMEX-RK schemes in the asymptotic limit, we develop the second class of schemes based on the SI-IMEX framework. A formal consistency analysis is then presented in Section 5 to establish its AP property.

### 4.1. Additively Partitioned IMEX-RK-DAE Schemes.
In order to design RK schemes for the EP system (2.1)-(2.3), henceforth we consider only its DAE formulation (3.13)-(3.14). As a first step, we identify an additive stiff-non-stiff splitting of the function H defined in (3.13). Note that the equation (3.15) already gives a natural splitting of H as the sum of the nonstiff hyperbolic flux F and the stiff source term S, cf. the asymptotic analysis presented in Section 2.1. This splitting corresponds to the so-called classical splitting reported in the literature; see [25, 31] for further details. Making use of the indirect approach from Section 3 with the above choice of F and S, we define an additively partitioned IMEX-RK-DAE scheme as follows.

**Definition 4.1.** An s-stage additively partitioned IMEX-RK-DAE scheme for the EP system (3.13)-(3.14) to update the solution $(U^n, \phi^n)$ at time $t^n$ to $(U^{n+1}, \phi^{n+1})$ at time $t^{n+1}$ uses the following s intermediate stages:

$$ (4.1) \qquad U^k = U^n - \Delta t \sum_{\ell=1}^{k-1} \tilde{a}_{k,\ell} F(U^\ell) + \Delta t \sum_{\ell=1}^{k} a_{k,\ell} S(U^\ell, \phi^\ell), \quad k = 1, 2, \dots, s, $$

$$ (4.2) \qquad G(U^k, \phi^k) = 0, \quad k = 1, 2, \dots, s, $$

and the final update,

$$ (4.3) \qquad (U^{n+1}, \phi^{n+1}) = (U^s, \phi^s). $$

Remark 4.2. Note that the above IMEX-RK schemes are stiffly accurate, and only DIRK variants are considered for the sake of maintaining simplicity and optimal computational performance.

The numerical solution is computed in the following chronology. For each $k = 1, 2, \dots, s$, first, the density $\rho^k$ is explicitly updated using the mass conservation update in the system (4.1). Using the
updated value of $\rho^k$ thus obtained, the linear elliptic problem (4.2) is solved for $\phi^k$. The momentum $q^k$ can then be obtained using the computed values of $\rho^k$ and $\phi^k$ from (4.1), which is now an explicit update.

Assuming that the data $\rho^n, q^n$, and $\phi^n$ at time $t^n$ for the above scheme (4.1)-(4.3) are well-prepared in the sense of Definition 2.1, for each $k = 1, 2, \dots, s$, taking the limit $\varepsilon \to 0$ in the Poisson equation (4.2) we obtain

$$ (4.4) \qquad \rho_{(0)}^k = 1, \text{ for } k = 1, 2, \dots, s. $$

Using the above equation in the mass updates for each $k = 1, 2, \dots, s$ we obtain

$$ (4.5) \qquad \nabla \cdot q_{(0)}^{k-1} = 0, \text{ for } k = 1, 2, \dots, s. $$

Therefore, in the limit $\varepsilon \to 0$, the additively partitioned IMEX-RK-DAE scheme (4.1)-(4.3) boils down to the following limiting scheme:

$$ (4.6) \qquad \rho_{(0)}^{n+1} = 1, $$

$$ (4.7) \qquad u_{(0)}^{n+1} = u_{(0)}^n - \Delta t \sum_{\ell=1}^{k-1} \tilde{a}_{k,\ell} \nabla \cdot (u_{(0)}^\ell \otimes u_{(0)}^\ell) + \Delta t \sum_{\ell=1}^{k} a_{k,\ell} \nabla \phi_{(0)}^\ell, $$

$$ (4.8) \qquad \nabla \cdot u_{(0)}^k = 0, \text{ for } k = 1, 2, \dots, s-1. $$

Note that the above limiting scheme doesn't provide a way for consistently caring forward the divergence-free condition to the numerical solution at time $t^{n+1}$. Precisely speaking, the IMEX-RK-DAE schemes (4.1)-(4.3) only offer the divergence-free condition up-till the $(s-1)^{th}$ stage via (4.8), and fail to provide it for $u_{(0)}^s = u_{(0)}^{n+1}$. A closer look at the limiting scheme (4.6)-(4.8) also reveals that, it lacks a way to obtain all the unknowns. Being fettered by the divergence-free condition on the leading order velocity, the IMEX-RK-DAE schemes fail to be consistent with the incompressible limit system (2.5)-(2.7) when $\varepsilon \to 0$ and also do not provide valid recursion to update all the variable at time $t^{n+1}$. Therefore, the IMEX-RK-DAE scheme cannot be AP.

A similar analysis of a first order version of the additively partitioned IMEX-RK-DAE scheme (4.1)-(4.3) can be found in [25]. In [25], the first order scheme is termed as the classical scheme which can be obtained in the present setting by using the IMEX-RK coefficients given in Figure 2. Consequently, the IMEX-RK-DAE schemes can be conceived as natural higher order extensions of the classical scheme. In [26], based on knowledge from plasma physics, it has been reported

<table><caption>FIGURE 2. Double Butcher tableaux of DIRK(1,1,1) scheme.</caption><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td></tr><tr><td></td><td>1</td><td>0</td><td></td><td>0</td><td>1</td></tr></table>

that the stability constraints for the classical scheme are extremely prohibitive in the quasineutral regime $\varepsilon \ll 1$. Therefore, combining this observation and the above presented proof of inconsistency, one can conclude that the IMEX-RK-DAE schemes derived from additively partitioned IMEX-RK schemes fail to qualify to being asymptotic preserving.

**Remark 4.3.** Note that it is possible to introduce other flux-splittings in the additive framework and develop IMEX-RK-DAE schemes possessing the AP property: e.g. the mass flux $q$ can be taken in $S$ instead of $F$. But, as a consequence of such a choice a non-linear elliptic problem has to be computationally tackled using iterative Newton solvers. In this paper the endeavour is to stay consistent with the continuous system which only has a linear elliptic problem.
**4.2. SI-IMEX-RK-DAE Schemes.** As a first step towards designing an SI-IMEX-RK scheme based on the indirect approach for the EP system (3.13)-(3.14), we specify the stiff and non-stiff terms. Consequently, it leads to the identification of terms that are to be treated implicitly or explicitly in $H$ and $G$. In the following, tracing the lines of SI-IMEX schemes presentation in Section 3, we subscript variables with $I$ and $E$, for those which are to be treated implicitly and explicitly, respectively. Based on the asymptotic analysis presented in Section 2, we recast the DAE form (3.13)-(3.14) of the EP system via the stiff and non-stiff conservative variables $U_I := (\rho_I, q_I)^T$ and $U_E := (\rho_E, q_E)^T$, and the stiff potential $\phi_I$ into,

$$ (4.9) \qquad \partial_t U + \mathbb{H}(U_E, U_I, \phi_I) = 0, $$

$$ (4.10) \qquad \mathbb{G}(U_I, \phi_I) = 0, $$

where

$$ (4.11) \qquad \mathbb{H}(U_E, U_I, \phi_I) := \begin{pmatrix} \nabla \cdot q_I \\ \nabla \cdot \left( \frac{q_E \otimes q_E}{\rho_E} \right) + \nabla p(\rho_E) - \rho_E \phi_I \end{pmatrix}, \\ \mathbb{G}(U_I, \phi_I) := -\epsilon^2 \Delta \phi_I - \rho_I + 1. $$

Using the indirect approach, we define the SI-IMEX-RK-DAE scheme for the above DAE formulation of the EP system as follows.

**Definition 4.4.** An $s$-stage SI-IMEX-RK-DAE scheme for the Euler-Poisson system (3.13)-(3.14) to update the solution $(U^n, \phi^n)$ at time $t^n$ to $(U^{n+1}, \phi^{n+1})$ at time $t^{n+1}$ uses the following $s$ intermediate steps:

$$ (4.12) \qquad U_E^k = U^n - \Delta t \sum_{\ell=1}^{k-1} \tilde{a}_{k,\ell} \mathbb{H} \left( U_E^\ell, U_I^\ell, \phi_I^\ell \right), \quad k = 1, 2, \dots, s, $$

$$ (4.13) \qquad U_I^k = U^n - \Delta t \sum_{\ell=1}^{k} a_{k,\ell} \mathbb{H} \left( U_E^\ell, U_I^\ell, \phi_I^\ell \right), \quad k = 1, 2, \dots, s, $$

$$ (4.14) \qquad \mathbb{G}(U_I^k, \phi_I^k) = 0, \quad k = 1, 2, \dots, s, $$

and the update stage

$$ (4.15) \qquad U^{n+1} = U_I^s, \quad \phi^{n+1} = \phi_I^s. $$

Even though the above scheme seems to have more implicit operations to be performed than the one based on additively partitioned approach, we will show that it can be equivalently reformulated so that it has the need to carry out same number of implicit operations as the scheme given in Definition 4.1.

*Remark 4.5.* Note that owing to Remark 3.7, the final update (4.15) of the SI-IMEX-RK-DAE scheme is same as the solution of the last implicit step.

**Definition 4.6.** The reformulation of the $s$-stage SI-IMEX-RK-DAE scheme (4.12)-(4.15) to update the solution $(\rho^n, q^n, \phi^n)$ at time $t^n$ to $(\rho^{n+1}, q^{n+1}, \phi^{n+1})$ at time $t^{n+1}$ uses $s$ intermediate stages and
an update stage as given below. For $k = 1, 2, \dots, s$, the fully explicit update of the $k^{th}$ stage reads

$$ (4.16) \qquad \rho_E^k = \rho^n - \Delta t \sum_{\ell=1}^{k-1} \tilde{a}_{k,\ell} \nabla \cdot q_I^\ell, $$

$$ (4.17) \qquad q_E^k = q^n - \Delta t \sum_{\ell=1}^{k-1} \tilde{a}_{k,\ell} \nabla \cdot \left( \frac{q_E^\ell \otimes q_E^\ell}{\rho_E^\ell} \right) - \Delta t \sum_{\ell=1}^{k-1} \tilde{a}_{k,\ell} \nabla p(\rho_E^\ell) + \Delta t \sum_{\ell=1}^{k-1} \tilde{a}_{k,\ell} \rho_E^\ell \nabla \phi_I^\ell. $$

The explicit part of the $k^{th}$ implicit stage is computed as

$$ (4.18) \qquad \hat{\rho}_I^k = \rho^n - \Delta t \sum_{\ell=1}^{k-1} a_{k,\ell} \nabla \cdot q_I^\ell, $$

$$ (4.19) \qquad \hat{q}_I^k = q^n - \Delta t \sum_{\ell=1}^{k} a_{k,\ell} \nabla \cdot \left( \frac{q_E^\ell \otimes q_E^\ell}{\rho_E^\ell} \right) - \Delta t \sum_{\ell=1}^{k} a_{k,\ell} \nabla p(\rho_E^\ell) + \Delta t \sum_{\ell=1}^{k-1} a_{k,\ell} \rho_E^\ell \nabla \phi_I^\ell. $$

The implicit solution is obtained by first solving the linear elliptic problem:

$$ (4.20) \qquad \nabla \cdot \left( \left( \varepsilon^2 + \Delta t^2 a_{k,k}^2 \rho_E^k \right) \nabla \phi_I^k \right) = \hat{\rho}_I^k - \Delta t a_{k,k} \nabla \cdot \hat{q}_I^k - 1, $$

and then performing the explicit updates:

$$ (4.21) \qquad q_I^k = \hat{q}_I^k + \Delta t a_{k,k} \rho_E^k \nabla \phi_I^k, $$

$$ (4.22) \qquad \rho_I^k = \hat{\rho}_I^k - \Delta t a_{k,k} \nabla \cdot q_I^k, $$

chronologically. Finally, the numerical solution is given by $(\rho^{n+1}, q^{n+1}, \phi^{n+1}) = (\rho_I^s, q_I^s, \phi_I^s)$.

**Proposition 4.7.** *The SI-IMEX-RK-DAE scheme (4.12)-(4.15) and the reformulated time semi-discrete scheme defined in Definition 4.6, are equivalent.*

*Proof.* We prove the proposition by using induction on the number of stages of the SI-IMEX-RK-DAE scheme. To start with, we consider the first stage, i.e. when $k = 1$, which corresponds to a trivial update for $\rho_E^1$ and $q_E^1$ given by

$$ (4.23) \qquad \rho_E^1 = \rho^n, \quad q_E^1 = q^n, $$

and $\rho_I^1, q_I^1$ and $\phi_I^1$ are obtained by solving the first implicit intermediate step:

$$ (4.24) \qquad \rho_I^1 = \rho^n - \Delta t a_{1,1} \nabla \cdot q_I^1, $$

$$ (4.25) \qquad q_I^1 = q^n - \Delta t a_{1,1} \nabla \cdot \left( \frac{q_E^1 \otimes q_E^1}{\rho_E^1} \right) - \Delta t a_{1,1} \nabla p(\rho_E^1) + \Delta t a_{1,1} \rho_E^1 \nabla \phi_I^1, $$

$$ (4.26) \qquad \varepsilon^2 \Delta \phi_I^1 = \rho_I^1 - 1. $$

Instead of solving the above linearly implicit system for $\rho_I^1, q_I^1$ and $\phi_I^1$, we use an elliptic reformulation. To this end we take the divergence of the momentum update (4.25) to obtain $\nabla \cdot q_I^1$ and, follow it up by substituting the expression thus obtained in place of $\nabla \cdot q_I^1$ in (4.24) to derive,

$$ (4.27) \qquad \rho_I^1 = \rho^n - \Delta t a_{1,1} \nabla \cdot q^n + \Delta t^2 a_{1,1}^2 \nabla^2: \left( \frac{q_E^1 \otimes q_E^1}{\rho_E^1} + p(\rho_E^1)I \right) - \Delta t^2 a_{1,1}^2 \nabla \cdot (\rho_E^1 \nabla \phi_I^1). $$

Next, we eliminate $\rho_I^1$ between (4.26) and (4.27) to get the following reformulated elliptic equation for $\phi_I^1$:

$$ (4.28) \qquad \nabla \cdot \left( \left( \varepsilon^2 + \Delta t^2 a_{1,1}^2 \rho_E^1 \right) \nabla \phi_I^1 \right) = \rho^n - 1 - \Delta t a_{1,1} \nabla \cdot q^n + \Delta t^2 a_{1,1}^2 \nabla^2: \left( \frac{q_E^1 \otimes q_E^1}{\rho_E^1} + p(\rho_E^1)I \right). $$
The above elliptic problem is solved to obtain $\phi_I^1$. The values of $q_I^1$ and $\rho_I^1$ are then computed using the anteriorly in-hand value of $\phi_I^1$, via the momentum update (4.25) and the mass update (4.24), chronologically. Note that these two updates performed following the prescribed sequence, are explicit evaluations for $q_I^1$ and $\rho_I^1$. Now for each $k = 2$ to $s$, we compute $\rho_E^k$ and $q_E^k$ explicitly using (4.12). Subsequently, using the evaluated values of $\rho_E^k$ and $q_E^k$, the explicit parts $\hat{\rho}_I^k$ and $\hat{q}_I^k$, of the implicit step (4.13) are calculated using (4.18)-(4.19). Finally, the implicit step (4.13) can be written in terms of $\hat{\rho}_I^k$ and $\hat{q}_I^k$, in the following compact form:

$$ (4.29) \qquad \rho_I^k = \hat{\rho}_I^k - \Delta t a_{k,k} \nabla \cdot q_I^k, $$

$$ (4.30) \qquad q_I^k = \hat{q}_I^k + \Delta t a_{k,k} \rho_E^k \nabla \phi_I^k, $$

$$ (4.31) \qquad \varepsilon^2 \Delta \phi_I^k = \rho_I^k - 1. $$

In order to solve (4.29)-(4.31) for $\rho_I^k$, $q_I^k$ and $\phi_I^k$, we use an elliptic reformulation, similar to the one performed for $k = 1$, followed by two explicit evaluations. Eliminating $\rho_I^k$ and $q_I^k$ amongst (4.29)-(4.31), yields the elliptic equation (4.20) for $\phi_I^k$. Once $\phi_I^k$ is known, $q_I^k$ and $\rho_I^k$ are computed using explicit evaluations of the momentum update (4.30) and the mass update (4.29), successively. Ultimately, the numerical solution $(\rho^{n+1}, q^{n+1}, \phi^{n+1})$ is updated as

$$ (4.32) \qquad \rho^{n+1} = \rho_I^s, \quad q^{n+1} = q_I^s, \quad \phi^{n+1} = \phi_I^s. $$

Therefore a solution to the SI-IMEX-RK-DAE scheme is a solution of the reformulated scheme in Definition 4.6. The converse can be proved in an analogous way, similar to [25]. □

We note that the above proof also serves as a step by step algorithm to sequentially obtain the numerical solution of the reformulated SI-IMEX-RK-DAE (4.16)-(4.22) scheme. We see that, each stage of the SI-IMEX-RK-DAE scheme comprises of three sub-steps: first, an explicit evaluation, second, solution of an elliptic equation, third and last, another explicit evaluation to update the solution of an implicit step. In other words, the numerical solution that we compute in this paper is not the solution of the actual SI-IMEX-RK-DAE scheme (4.12)-(4.15); rather it is the solution of a reformulated Poisson equation (4.20) coupled with the mass (4.22) and momentum (4.21) updates of the original numerical scheme (4.12)-(4.15) which are now explicit in nature.

**Remark 4.8.** Although the SI-IMEX scheme (4.12)-(4.15) for the non-linear EP system (2.1)-(2.3) is formulated as a semi-implicit scheme, the reformulation stated in Definition 4.6 shows that it is only linearly implicit.

As observed for the additively partitioned IMEX-RK-DAE scheme here also we see that by making a suitable choice of the RK tableaux; see Figure 3, the first order AP scheme presented in [25, 26] can be obtained as:

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

$$ (4.34) \qquad q^{n+1} = q^n - \Delta t \nabla \cdot \left( \frac{q^n \otimes q^n}{\rho^n} + p(\rho^n) I \right) + \Delta t \rho^n \nabla \phi^{n+1}, $$

$$ (4.35) \qquad \varepsilon^2 \Delta \phi^{n+1} = \rho^{n+1} - 1. $$

That being the case we come to the understanding that the SI-IMEX-RK-DAE scheme (4.12)-(4.15) serves as a natural higher order extension to the first order scheme (4.33)-(4.35). We refer the reader to [25, 26] and the references therein for a detailed study of the first order scheme; particularly, its computational cost compared to the so-called classical scheme; cf. [25, 26]. As stated there it is imperative to understand that the computational time expense typically associated with any amount of implicitness doesn't affect our scheme adversely.
$$ \begin{array}{c|c} 0 & 0 \\ \hline 1 & 1 \end{array} \quad \begin{array}{c|c} 1 & 1 \\ \hline 1 & 1 \end{array}. $$

<caption>FIGURE 3. Double Butcher tableau of SI-IMEX-Euler (1,1,1) scheme.</caption>

## 5. ANALYSIS OF THE TIME SEMI-DISCRETE SI-IMEX-RK-DAE SCHEME

In this section we study the AP property of the time semi-discrete scheme introduced in Section 4.2. It has been detailed, e.g. in [22, 23, 26], that standard (classical) discretisation methods applied to the Euler-Poisson system do not yield the correct quasineutral limit when subjected to very small Debye length $\epsilon$. In order to overcome this shortcoming of classical discretisation techniques, an AP scheme which is first order accurate is proposed and analysed in the above references.

We have observed in the preceding section that the AP scheme presented in the above mentioned references is a first order variant of the proposed SI-IMEX-RK-DAE scheme. Under the AP framework these schemes are designed to perform for a wide range of $\epsilon$: from non-stiff regimes $\epsilon = O(1)$ to stiff regimes where $\epsilon \to 0$. A scheme's attribute of being AP reflects that it performs uniformly with respect to $\epsilon$. The performance referred here pertains to the unprejudiced stability - asymptotic stability, and consistency in the limit of the singular parameter - asymptotic consistency of the method. In what follows we provide a formal proof for the asymptotic consistency of the time semi-discrete SI-IMEX-RK-DAE scheme with the quasineutral limit system and also post a short discussion regarding the asymptotic stability of the scheme; see also [26] for related discussions.

### 5.1. Asymptotic Consistency.

**Theorem 5.1.** Suppose that the data at time $t^n$ are well-prepared, i.e. $\rho^n$ and $q^n$ admit the decomposition (2.8)-(2.9). Then, for a SA-DIRK scheme, if the intermediate solutions $(\rho_E^k, q_E^k, \phi_E^k)$, $(\rho_I^k, q_I^k, \phi_I^k)$ admit the ansatz (2.8), then it follows that

$$ (5.1) \qquad \lim_{\epsilon \to 0} \rho_E^k = \rho_{E,(0)}^k = 1, \quad \lim_{\epsilon \to 0} \rho_I^k = \rho_{I,(0)}^k = 1, $$

$$ (5.2) \qquad \lim_{\epsilon \to 0} \nabla \cdot q_E^k = \nabla \cdot q_{E,(0)}^k = 0, \quad \lim_{\epsilon \to 0} \nabla \cdot q_I^k = \nabla \cdot q_{I,(0)}^k = 0, $$

i.e. they are well-prepared. As a result, the numerical solution $(\rho^{n+1}, q^{n+1}, \phi^{n+1})^T$ after one timestep is consistent with the reformulated incompressible limit (2.14)-(2.11); in the sense that an SI-IMEX-RK-DAE scheme based on an SA-DIRK scheme transforms to an s-stage scheme for the incompressible Euler system in the quasineutral limit $\epsilon \to 0$. In other words, the SI-IMEX-RK-DAE scheme (4.12)-(4.14) is asymptotically consistent.

*Proof.* We use induction on the number stages (s) to prove the theorem. Owing to the equivalence between the SI-IMEX-RK-DAE scheme (4.12)-(4.15) and its reformulation (4.16)-(4.22) cf. Proposition 4.7, we prove the consistency for the reformulated scheme instead.

As a starting step, consider the first stage, i.e. for $k=1$. We have from (4.16) and (4.17)

$$ (5.3) \qquad \rho_{E,(0)}^1 = \lim_{\epsilon \to 0} \rho_E^1 = \rho_{(0)}^n = 1, $$

$$ (5.4) \qquad \nabla \cdot q_{E,(0)}^1 = \lim_{\epsilon \to 0} \nabla \cdot q_E^1 = \nabla \cdot q_{(0)}^n = 0. $$

As a consequence of the above two equations, we obtain the divergence condition, $\nabla \cdot u_{E,(0)}^1 = 0$. Next, we take equation (4.20) into consideration. In the limit of $\epsilon \to 0$, it furnishes the following elliptic equation, for the limiting electric potential $\phi_{I,(0)}^1$:

$$ (5.5) \qquad \Delta \phi_{I,(0)}^1 = \nabla^2: (u_{E,(0)}^1 \otimes u_{E,(0)}^1). $$
We proceed further by letting $\epsilon \to 0$ in the implicit mass and momentum updates, (4.22) and (4.21), respectively, to yield

$$ (5.6) \qquad \rho_{I,(0)}^1 = 1 - \Delta t a_{1,1} \nabla \cdot q_{I,(0)}^1, $$

$$ (5.7) \qquad q_{I,(0)}^1 = u_{(0)}^n - \Delta t a_{1,1} \nabla \cdot (u_{E,(0)}^1 \otimes u_{E,(0)}^1) + \Delta t a_{1,1} \nabla \phi_{I,(0)}^1. $$

Taking divergence of the above equation (5.7) and using the value of $\Delta\phi_{I,(0)}^1$ from (5.5) in the expression thus obtain, we get the following:

$$ (5.8) \qquad \nabla \cdot q_{I,(0)}^1 = 0, $$

Now, using the above divergence-free expression for the limit momentum in the mass update (5.6) forthright yields the quasineutrality condition:

$$ (5.9) \qquad \rho_{I,(0)}^1 = 1, $$

for the limiting density. Summarising the above steps, for the first stage of the limiting scheme, we gather the following limiting scheme for $k=1$.

For the explicit part:

$$ (5.10) \qquad \rho_{E,(0)}^1 = \rho_{(0)}^n = 1, $$

$$ (5.11) \qquad \nabla \cdot q_{E,(0)}^1 = \nabla \cdot q_{(0)}^n = 0. $$

For the implicit part:

$$ (5.12) \qquad \rho_{I,(0)}^1 = 1, $$

$$ (5.13) \qquad u_{I,(0)}^1 = u_{(0)}^n - \Delta t a_{1,1} \nabla \cdot (u_{E,(0)}^1 \otimes u_{E,(0)}^1) + \Delta t a_{1,1} \nabla \phi_{I,(0)}^1, $$

$$ (5.14) \qquad \Delta \phi_{I,(0)}^1 = \nabla^2 : (u_{E,(0)}^1 \otimes u_{E,(0)}^1), $$

which provisions the divergence free condition, $\nabla \cdot q_{I,(0)}^1 = \nabla \cdot u_{I,(0)}^1 = 0$. Clearly, (5.12)-(5.14) is a consistent discretisation of the reformulated incompressible Euler system (2.5)-(2.7).

For $k=2$ onwards, proceeding as done in the case of $k=1$ we can obtain

For the explicit part:

$$ (5.15) \qquad \rho_{E,(0)}^k = \rho_{(0)}^n = 1, $$

$$ (5.16) \qquad u_{E,(0)}^k = u_{(0)}^n - \Delta t \sum_{\ell=1}^{k-1} \tilde{a}_{k,\ell} \nabla \cdot (u_{E,(0)}^\ell \otimes u_{E,(0)}^\ell) + \Delta t \sum_{\ell=1}^{k-1} \tilde{a}_{k,\ell} \nabla \phi_{I,(0)}^\ell, $$

$$ \text{and consequently, } \nabla \cdot q_{E,(0)}^k = \nabla \cdot u_{E,(0)}^k = 0. $$

For the implicit part:

$$ (5.17) \qquad \rho_{I,(0)}^k = 1, $$

$$ (5.18) \qquad u_{I,(0)}^k = u_{(0)}^n - \Delta t \sum_{\ell=1}^{k} a_{k,\ell} \nabla \cdot (u_{E,(0)}^\ell \otimes u_{E,(0)}^\ell) + \Delta t \sum_{\ell=1}^{k} a_{k,\ell} \nabla \phi_{I,(0)}^\ell, $$

$$ (5.19) \qquad \Delta \phi_{I,(0)}^k = \nabla^2 : (u_{E,(0)}^k \otimes u_{E,(0)}^k), $$

and ergo, the divergence condition:

$$ (5.20) \qquad \nabla \cdot q_{I,(0)}^k = \nabla \cdot u_{I,(0)}^k = 0. $$
Since the SI-IMEX-RK scheme under consideration is GSA, the numerical solution at $t^{n+1}$ is same as the solution of the implicit step for $k = s$. Therefore we have

$$ (5.21) \qquad \rho_{(0)}^{n+1} = 1, $$

$$ (5.22) \qquad u_{(0)}^{n+1} = u_{(0)}^{n} - \Delta t \sum_{\ell=1}^{s} a_{s,\ell} \nabla \cdot (u_{E,(0)}^{\ell} \otimes u_{E,(0)}^{\ell}) + \Delta t \sum_{\ell=1}^{s} a_{s,\ell} \nabla \phi_{I,(0)}^{\ell}, $$

$$ (5.23) \qquad \Delta \phi_{(0)}^{n+1} = \nabla^2 : (u_{E,(0)}^s \otimes u_{E,(0)}^s), $$

and consequently the divergence condition $\nabla \cdot q_{I,(0)}^{n+1} = \nabla \cdot u_{I,(0)}^{n+1} = 0$, i.e. the numerical solution $(\rho^{n+1}, q^{n+1}, \phi^{n+1})$ is also well prepared. Hence, we conclude that the limiting scheme is a consistent approximation of the reformulated incompressible system. $\square$

**Remark 5.2.** Note that even though the data at time level $t^n$ is non-wellprepared, still the SI-IMEX-RK scheme will be able to convert the following iterates into well-prepared data. For example, if the momentum $q^n$ is not well-prepared but the density $\rho^n$ is, then owing to the elliptic problem for $\phi^{n+1}$ leads to the divergence free condition for $q^{n+1}$ in the limit of $\varepsilon \to 0$. Similarly, if none of the initial data were well-prepared then it takes two iterates i.e. at time level $t^{n+2}$ the data is made well-prepared by the SI-IMEX-RK scheme. This is in alignment with the remark in [25] regarding the first order AP scheme.

## 5.2. Asymptotic Stability.
In [26] the authors have carried out a thorough linearised $L^2$-stability analysis of a semi-implicit time stepping scheme which corresponds to the first order variant of the SI-IMEX-RK-DAE scheme given in Figure 3. Through a Fourier analysis which captures the effects of a centred space discretisation it has been discovered that the above scheme fails to be uniformly stable. As a remedy to this some amount of artificial diffusions have been introduced into the mass and the momentum updates to retain stability; see also [23]. Though the numerical viscosity added is minimal, as a consequence, the transformed system corresponds to an upwind space discretisation. Furthermore, it has been shown the stability characteristics of the corrected scheme do not degrade with $\varepsilon$ in the asymptotic range and elsewhere as well; see also [24] for a similar analysis towards proving the stability of a first order SI-IMEX-RK-DAE scheme.

As is concluded in [26] and reiterated in [24], it is imperative to add some amount of artificial viscosity to the hydrodynamic equations, for the numerical scheme to be stable in the fully discrete setup. Adhering to this we introduce the necessary viscosity terms to the SI-IMEX-RK-DAE scheme while designing the space-time fully discrete scheme to be presented in Section 6. A detailed stability analysis for the second order version of the SI-IMEX-RK-DAE scheme is quite a tasking job and is a subject matter for some future work.

# 6. SPACE-TIME FULLY-DISCRETE SCHEME AND ITS ANALYSIS

This section is devoted to the construction of a space-time fully-discrete scheme corresponding to the SI-IMEX-RK-DAE scheme introduced in Section 3, and its asymptotic analysis. To this end, we avail a similar approach as is employed in [1], where, a fully-discrete scheme similar in its spirit to the one to be defined in this paper, for the low Mach number limit of the compressible Euler system was obtained using a coalescence of standard Rusanov-type flux for the explicit fluxes and second-order central differencing for the implicit derivative terms. A first-order fully discrete schemes and its asymptotic analysis can be found in [25, 26]. Subsequently, in this section, we also briefly discuss the asymptotic preserving properties of the fully discrete scheme, thus to be defined.
**6.1. Space-Time Fully-Discrete Scheme.** For ease of presentation, we introduce the following notations for the flux functions and the source term:

$$ (6.1) \qquad F_m(U) = \begin{pmatrix} F_{\rho,m}(U) \\ F_{q,m}(U) \end{pmatrix} = \begin{pmatrix} q_m \\ \frac{q_m}{\rho}q + pe_m \end{pmatrix}, \quad S_m(\phi) = \begin{pmatrix} S_{\rho,m}(\phi) \\ S_{q,m}(\phi) \end{pmatrix} = \begin{pmatrix} 0 \\ \phi e_m \end{pmatrix}, $$

where $e_m$ denotes the $m^{th}$ unit vector in $\mathbb{R}^3$, for $m = 1, 2, 3$. We rewrite the reformulated time semi-discrete scheme (4.16)-(4.22) in terms of (6.1) as follows. For the explicit stages,

$$ (6.2) \qquad \rho_E^k = \rho^n - \Delta t \sum_{\ell=1}^{k-1} \tilde{a}_{k,\ell} \sum_{m=1}^{3} \partial_{x_m} F_{\rho,m}(U_I^\ell), $$

$$ (6.3) \qquad q_E^k = q^n - \Delta t \sum_{\ell=1}^{k-1} \tilde{a}_{k,\ell} \sum_{m=1}^{3} \partial_{x_m} F_{q,m}(U_E^\ell) + \Delta t \sum_{\ell=1}^{k-1} \tilde{a}_{k,\ell} \rho_E^\ell \sum_{m=1}^{3} \partial_{x_m} S_{q,m}(\phi_I^\ell). $$

For the explicit part of the implicit stages,

$$ (6.4) \qquad \hat{\rho}_I^k = \rho^n - \Delta t \sum_{\ell=1}^{k-1} a_{k,\ell} \sum_{m=1}^{3} \partial_{x_m} F_{\rho,m}(U_I^\ell), $$

$$ (6.5) \qquad \hat{q}_I^k = q^n - \Delta t \sum_{\ell=1}^{k} a_{k,\ell} \sum_{m=1}^{3} \partial_{x_m} F_{q,m}(U_E^\ell) + \Delta t \sum_{\ell=1}^{k-1} a_{k,\ell} \rho_E^\ell \sum_{m=1}^{3} \partial_{x_m} S_{q,m}(\phi_I^\ell). $$

Finally, the implicit solution is obtained by solving,

$$ (6.6) \qquad \sum_{m=1}^{3} \partial_{x_m} \left( \left( \varepsilon^2 + \Delta t^2 a_{k,k}^2 \rho_E^k \right) \partial_{x_m} \phi_I^k \right) = \hat{\rho}_I^k - \Delta t a_{k,k} \sum_{m=1}^{3} \partial_{x_m} \hat{q}_{I,m}^k - 1, $$

and following up with the updates,

$$ (6.7) \qquad q_{I,m}^k = \hat{q}_{I,m}^k + \Delta t a_{k,k} \rho_E^k \sum_{m=1}^{3} \partial_{x_m} \phi_I^k, $$

$$ (6.8) \qquad \rho_I^k = \hat{\rho}_I^k - \Delta t a_{k,k} \sum_{m=1}^{3} \partial_{x_m} q_{I,m}^k. $$

In order to present the fully-discrete scheme, we introduce the vector $i = (i_1, i_2, i_3)$ where each $i_m$ for $m = 1, 2, 3$ represent $x_m$ the space direction. We further introduce the following finite difference and averaging operators: e.g. in the $x_m$-direction

$$ (6.9) \qquad \delta_{x_m} w_i = w_{i+\frac{1}{2}e_m} - w_{i-\frac{1}{2}e_m}, \quad \mu_{x_m} w_i = \frac{w_{i+\frac{1}{2}e_m} + w_{i-\frac{1}{2}e_m}}{2}, $$

for any grid function $w_i$.

**Definition 6.1.** The $k^{th}$ stage of an $s$-stage fully-discrete SI-IMEX-RK-DAE scheme for the Euler-Poisson system is defined as follows. The explicit solution is obtained by

$$ (6.10) \qquad \rho_{E,i}^k = \rho_i^n - \sum_{\ell=1}^{k-1} \tilde{a}_{k,\ell} \sum_{m=1}^{3} \nu_m \delta_{x_m} \mathcal{F}_{\rho,m,i}^{\ell} $$

$$ (6.11) \qquad q_{E,i}^k = q_i^n - \sum_{\ell=1}^{k-1} \tilde{a}_{k,\ell} \sum_{m=1}^{3} \nu_m \delta_{x_m} \mathcal{F}_{q,m,i}^{\ell} + \sum_{\ell=1}^{k-1} \tilde{a}_{k,\ell} \sum_{m=1}^{3} \nu_m \rho_{E,i}^{\ell} \delta_{x_m} S_{q,m,i}^{\ell}. $$
In an analogous manner, the explicit part of the implicit stage is computed via,

$$ (6.12) \qquad \hat{\rho}_{I,i}^k = \rho_i^n - \sum_{\ell=1}^{k-1} a_{k,\ell} \sum_{m=1}^3 \nu_m \delta_{x_m} \mathcal{F}_{\rho,m,i}^{\ell}, $$

$$ (6.13) \qquad \hat{q}_{I,i}^k = q_i^n - \sum_{\ell=1}^{k} a_{k,\ell} \sum_{m=1}^{3} \nu_m \delta_{x_m} \mathcal{F}_{q,m,i}^{\ell} + \sum_{\ell=1}^{k-1} a_{k,\ell} \sum_{m=1}^{3} \nu_m \rho_{E,i}^{\ell} \delta_{x_m} \mathcal{S}_{q,m,i}^{\ell}. $$

Then the implicit solution is obtained by first solving the discrete elliptic problem which is a linear system for {$\phi_{I,i}^k$} given by,

$$ (6.14) \qquad \sum_{m=1}^{3} \frac{\delta_{x_m}}{\Delta x_m} \left( \left( \varepsilon^2 + \rho_{E,i}^k \right) \frac{\delta_{x_m}}{\Delta x_m} \phi_{I,i}^k \right) = \hat{\rho}_{I,i}^k - a_{k,k} \sum_{m=1}^{3} \nu_m \delta_{x_m} \hat{q}_{I,m,i}^k, $$

followed by the explicit evalutions,

$$ (6.15) \qquad q_{I,i}^k = \hat{q}_{I,i}^k + a_{k,k} \sum_{m=1}^{3} \nu_m \rho_{E,i}^k \delta_{x_m} \phi_{I,i}^k, $$

$$ (6.16) \qquad \rho_{I,i}^k = \hat{\rho}_{I,i}^k - a_{k,k} \sum_{m=1}^{3} \nu_m \delta_{x_m} q_{I,m,i}^k + a_{k,k} \sum_{m=1}^{3} \nu_m \delta_{x_m} \mathcal{G}_{I,m,i}^k. $$

Here, the repeated index $m$ takes values in $\{1,2,3\}$, $\nu_m := \frac{\Delta t}{\Delta x_m}$ denote the mesh ratios and $\mathcal{F}_{\rho,m}$, $\mathcal{F}_{q,m}$ and $\mathcal{S}_m$ are, respectively, the numerical fluxes and numerical source term used to approximate the physical fluxes and physical source term $F_{\rho,m}$, $F_{q,m}$ and $\mathcal{S}_m$. Analogous to [1], in this section we have used a simple Rusanov-type flux for $\mathcal{F}_{q,m}$ and central flux for $\mathcal{F}_{\rho,m}$ and $\mathcal{S}_m$. The term $\mathcal{G}_{I,m}$ denotes the artificial viscosity added in the mass equation to adhere to the stability requirements discussed in Section 5. The numerical fluxes, source terms and the viscosity terms are defined as
$$ (6.17) \qquad \mathcal{F}_{\rho,m,i+\frac{1}{2}e_m}^k = \frac{1}{2} \left( F_{\rho,m} \left( U_{I,i+e_m}^k \right) + F_{\rho,m} \left( U_{I,i}^k \right) \right), $$

$$ \mathcal{F}_{\rho,m,i+\frac{1}{2}e_m}^k = \frac{1}{2} \left( F_{\rho,m} \left( U_{E,i+\frac{1}{2}e_m}^{k,+} \right) + F_{\rho,m} \left( U_{E,i+\frac{1}{2}e_m}^{k,-} \right) \right) - \frac{\alpha_{m,i+\frac{1}{2}e_m}}{2} \left( q_{m,E,i+\frac{1}{2}e_m}^{k,+} - q_{m,E,i+\frac{1}{2}e_m}^{k,-} \right), $$

$$ \mathcal{F}_{q,m,i+\frac{1}{2}e_m}^k = \frac{1}{2} \left( F_{q,m} \left( U_{E,i+\frac{1}{2}e_m}^{k,+} \right) + F_{q,m} \left( U_{E,i+\frac{1}{2}e_m}^{k,-} \right) \right) - \frac{\alpha_{m,i+\frac{1}{2}e_m}}{2} \left( q_{m,E,i+\frac{1}{2}e_m}^{k,+} - q_{m,E,i+\frac{1}{2}e_m}^{k,-} \right), $$

$$ \mathcal{S}_{q,m,i+\frac{1}{2}e_m}^k = \frac{1}{2} \left( S_{q,m} \left( \phi_{I,i+e_m}^k \right) + S_{q,m} \left( \phi_{I,i}^k \right) \right), \quad \mathcal{G}_{I,m,i+\frac{1}{2}e_m}^k = \frac{1}{\nu_m} \left( \rho_{E,i+e_m}^k - \rho_{E,i}^k \right). $$

Here, for any conservative variable $w$, we have denoted by $w_{i+\frac{1}{2}e_m}^{\pm}$, the interpolated states obtained using the piecewise linear reconstructions. The wave-speeds are computed as, e.g. in the $x_1$-direction

$$ (6.18) \qquad \alpha_{1,i+\frac{1}{2}e_1} := 2 \max \left( \left| \frac{q_{1,E,i+\frac{1}{2}e_1}^{k,-}}{\rho_{E,i+\frac{1}{2}e_1}^{k,-}} \right|, \left| \frac{q_{1,E,i+\frac{1}{2}e_1}^{k,+}}{\rho_{E,i+\frac{1}{2}e_1}^{k,+}} \right| \right). $$

*Remark 6.2.* Note that we have only introduced additional viscosity into the mass updates but not to the momentum equations. This choice is in accordance with the Rusanov-type flux used for the explicit momenta flux which introduces the required viscosity by default.
Remark 6.3. The space discretisation of the time semi-discrete high order classical scheme (4.1)-(4.3) is given as

$$ (6.19) \qquad U_i^{n+1} = U_i^n - \sum_{m=1}^{3} \nu_m \delta_{x_m} \mathcal{F}_{m,i}^n + \sum_{m=1}^{2} \nu_m \rho_i^{n+1} \delta_{x_m} \mathcal{S}_{m,i}^{n+1} $$

$$ (6.20) \qquad \varepsilon^2 \sum_{m=1}^{3} \left( \frac{\delta_{x_m}}{\Delta x_m} \right)^2 \phi_i^{n+1} = \rho_i^{n+1} - 1, $$

Where $\mathcal{F}_m$ in (6.19) is approximated using Rusanov fluxes similar to (6.17).

The eigenvalues of the Jacobians of $F_m(U)$ with respect to $U$ can be obtained as $\lambda_{m,1} = \lambda_{m,2} = \frac{q_m}{\rho}$, $\lambda_{m,3} = \frac{q_m}{\rho} - \sqrt{p'(\rho)}$ and $\lambda_{m,4} = \frac{q_m}{\rho} + \sqrt{p'(\rho)}$. Since the fluxes $F_m(U)$ for the classical scheme are entirely approximated by a Rusanov-type flux, the timestep $\Delta t$ at time $t^n$ is computed by the CFL condition:

$$ (6.21) \qquad \Delta t \max_i \max_m \left( \frac{|\lambda_{m,4}(U_i^n)|}{\Delta x_m} \right) = \nu, $$

where $\nu < 1$ is the given CFL number. For the SI-IMEX-RK-DAE scheme the mass flux is approximated using simple central differences. Hence, the eigenvalues of the Jacobians of the part of the flux which is approximated by Rusanov-type flux can be obtained as $\lambda_{m,1} = 0, \lambda_{m,2} = \lambda_{m,3} = \frac{q_m}{\rho}$ and $\lambda_{m,4} = 2\frac{q_m}{\rho}$ and the CFL condition is given by (6.21).

A detailed stability analysis of the first order AP scheme is carried in [26], where it has been shown that the scheme is stable under a CFL-like condition independent of $\varepsilon$. However, the classical scheme (6.19)-(6.20), as it involves the solution of an elliptic problem containing the stiffness parameter $\varepsilon$, should additionally satisfy the stability condition [31] given by

$$ (6.22) \qquad \Delta t \le \varepsilon. $$

**6.2. Asymptotic Preserving Property.** Here, we follow the approach of [1] to analyse the high order space-time fully discrete scheme, in-turn demonstrating that the space discretisation complements the designed time-discretisation optimally. The analysis of a first-order space-time fully discrete scheme for the Euler-Poisson system can be found in [23, 26]. We prove that the space discretisation coupled with the time-discretisation is indeed asymptotically consistent with the quasineutral limit, i.e. the scheme boils down to a consistent discretisation of the incompressible Euler system when $\varepsilon \to 0$. The stability properties of a first order version of the SI-IMEX-RK-DAE scheme have been studied in [26]. We do not intend to provide the stability considerations here.

**Definition 6.4 (Well-prepared data).** A numerical solution $(\rho_i, q_i, \phi_i)$ of the Euler-Poisson system is called well prepared if

$$ (6.23) \qquad \rho_i = \rho_{(0),i} + \varepsilon^2 \rho_{(2),i}, \qquad q_i = q_{(0),i} + \varepsilon^2 q_{(2),i} $$

$$ (6.24) \qquad \text{and,} \qquad \rho_{(0),i} = 1, \qquad \hat{\nabla} \cdot q_{(0),i} = 0. $$

Here, $\hat{\nabla}$ is the discrete gradient operator introduced by the central numerical fluxes, i.e. by replacing the derivatives by central differences.

**Theorem 6.5.** Suppose that the data at time $t^n$ are well-prepared, i.e. $\rho^n$ and $q^n$ admit the decompositions in (6.23). Then for an SA-DIRK scheme, the intermediate solutions $(\rho_I^k, q_I^k, \phi_I^k)$ satisfy

$$ (6.25) \qquad \lim_{\varepsilon \to 0} \rho_{I,i}^k = \rho_{I,(0),i}^k = 1 $$

$$ (6.26) \qquad \lim_{\varepsilon \to 0} \hat{\nabla} \cdot q_{I,i}^k = \hat{\nabla} \cdot q_{I,(0),i}^k = 0. $$
Moreover, the numerical solution $(\rho_i^{n+1}, q_i^{n+1}, \phi_i^{n+1})$ after one time-step is consistent with the reformulated incompressible limit (2.11)-(2.14). Therefore, the scheme transforms to an s-stage scheme for the limit system or in other words, the fully discrete scheme in Definition 6.1 is asymptotically consistent.

*Proof.* The proof follows from a combination of the proof of the AP property for the time semi-discrete system i.e. Theorem 5.1 and the use of mimetic operators as done in [1]. □

## 7. NUMERICAL CASE STUDIES

In this section we present the results of numerical case studies carried out with both the high-order schemes developed in this paper. For all the test cases we impose a combination of periodic and homogeneous Dirichlet boundary conditions for the hydrodynamic variables and the electric potential, respectively. The results presented here focus to numerically corroborate the following three properties as claimed throughout the paper:

(i) asymptotic consistency with the quasineutral limit system;

(ii) asymptotic convergence, i.e. second order convergence to the quasineutral limit solution;

(iii) superior performance of the SI-IMEX-RK-DAE scheme over high order classical schemes.

To this end we consider three test problems each one demonstrating some of the above stated properties of the SI-IMEX-RK-DAE scheme. Note that we use the LSDIRK(2,2,2) [9] variant of the SI-IMEX-RK scheme and the ARS(2,2,2) [67] variant of the additively partitioned IMEX-RK scheme. Note that both these variants are stiffly accurate and L-stable and their Butcher tableaux are given in Figure 4.

<table><caption>FIGURE 4. Double Butcher tableaux of the IMEX schemes. Left: ARS (2,2,2), and Right: LSDIRK(2,2,2). Here, $\gamma = 1 - \sqrt{2}/2$, $\sigma = 1/2\gamma$ and $\delta = 1 - \sigma$.</caption><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><td>0</td><td>0</td><td>0</td><td>0</td><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><td>γ</td><td>0</td><td>γ</td><td>0</td><td>γ</td><td>0</td><td>γ</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><td>1</td><td>0</td><td>1 - γ</td><td>γ</td><td>1</td><td>0</td><td>1 - γ</td><td>γ</td><td>1</td><td>0</td><td>1 - γ</td><td>γ</td></tr><tr><td colspan="2">δ</td><td>1 - δ</td><td>0</td><td colspan="2">0</td><td>1 - γ</td><td>γ</td><td colspan="2">0</td><td>1 - γ</td><td>γ</td><td colspan="2">0</td><td>1 - γ</td><td>γ</td><td>1</td><td>0</td><td>1 - γ</td><td>γ</td></tr></table>

### 7.1. A Small Perturbation of a Uniform Quasineutral Plasma.
We consider a quasineutral state where the density $\rho$ is a constant equal to one and the plasma is assumed to move in the horizontal direction with a uniform constant velocity one; see [23]. The constant value of density combined with the Poisson equation implies that the electric potential vanishes, i.e. $\phi = 0$. The aim of this case study is to test the AP scheme's ability to recover a quasineutral background state independent of the mesh-size: resolving or not resolving the parameter $\varepsilon$. To this end, a small cosine perturbation of magnitude $\delta = \varepsilon^2$ is added to the constant horizontal velocity $u_1$. The initial conditions read

$$ (7.1) \qquad \rho(0, x_1) = 1.0, \quad u_1(0, x_1) = 1.0 + \delta \cos 2\pi x_1, \quad \phi(0, x_1) = 0.0. $$

Note that the perturbation chosen here pertains to a set of well-prepared initial data cf. Definition 2.1. The parameter $\varepsilon$ is chosen as $\varepsilon = 10^{-4}$ and hence the plasma frequency $\omega = 1/\varepsilon = 10^4$. The plasma is considered in the domain $[0, 1]$. Throughout this test case the AP scheme is subjected to the advective CFL condition 6.21. The CFL number is set to be 0.45.

Figure 5 shows the density, velocity and electric potential at time $t = 0.01$ (100 cycles) computed on a mesh with $\Delta x = 10^{-4}$, resolving the Debye length. Similarly, Figure 6 shows the above mentioned components at time $t = 0.1$ (1000 cycles) computed on the same resolved mesh. The
time steps for the classical scheme is not only restricted by the CFL condition given by (6.21) but also by the severe stability restriction (6.22). It can be seen that both the classical and the AP schemes perform quite well to maintain the well-preparedness of the data over time for each of the components.

<figure><img src="image_2.png" /><figcaption>FIGURE 5. Small perturbation of a quasineutral state with AP and Classical Schemes. $\epsilon = 10^{-4}$, $N = 10^4$ (Resolved mesh). Left: Density Profile, Center: Velocity Profile, Right: Electric Potential Profiles, at time $t = 10^{-2}$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 6. Small perturbation of a quasineutral state with AP and Classical Schemes. $\epsilon = 10^{-4}$, $N = 10^4$ (Resolved mesh). Left: Density Profile, Center: Velocity Profile, Right: Electric Potential Profiles, at time $t = 10^{-1}$.</figcaption></figure>

Figure 7 shows the density, velocity and electric potential at time $t = 0.01$ computed on a coarse mesh with $\Delta x = 10^{-2}$, not resolving the Debye length. Figure 8 shows the corresponding components at time $t = 0.1$ on the same coarse mesh. As was done in the case of the fine mesh resolving the Debye length here also the classical scheme is subjected to the CFL conditions in (6.21) and (6.22). Both the figures convey that the classical as well as the AP scheme perform well and retain the well-preparedness aspect of the solution observed with fine mesh computations.

*Remark 7.1.* Note that in the numerical case studies presented until now the classical scheme is restricted to a rather stringent stability restriction and its performance which is seen to be similar to that of the AP scheme in case of the coarse mesh can be attributed to this strict time-stepping criteria.

Owing to the above remark we subject the classical scheme to satisfy the CFL condition (6.21) only, and therefore $\Delta t > \epsilon$ on a coarse mesh $\Delta x = 10^{-2}$. Figure 9 shows that the numerical solution blows up after a very small time $t = 6.1 \times 10^{-3}$ and hence inferring that the classical scheme is unstable for $\Delta t > \epsilon$. Therefore, the AP scheme is able to capture a solution close to the stationary
<figure><img src="image_4.png" /><figcaption>FIGURE 7. Small perturbation of a quasineutral state with AP ($\Delta t > \varepsilon$) and Classical Schemes ($\Delta t < \varepsilon$). $\varepsilon = 10^{-4}$, $N = 10^2$ (Unresolved mesh). Left: Density Profile, Center: Velocity Profile, Right: Electric Potential Profiles, at time $t = 10^{-2}$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 8. Small perturbation of a quasineutral state with AP ($\Delta t > \varepsilon$) and Classical Schemes ($\Delta t < \varepsilon$). $\varepsilon = 10^{-4}$, $N = 10^2$ (Unresolved mesh). Left: Density Profile, Center: Velocity Profile, Right: Electric Potential Profiles, at time $t = 10^{-1}$.</figcaption></figure>

state with a coarse mesh, whereas the classical scheme needs a very small time step to be stable, otherwise it blows up very rapidly. This is in agreement to the results provided in [23] for the first order classical scheme for the two-fluid EP system.

## 7.2. Slight Perturbation of a Maxwellian.
Here we consider a non-well prepared initial data. To this end we consider a test case from [65] where a small perturbation of amplitude $\delta = 10^{-2}$ to a Maxwellian was considered. The initial density profile reads:

$$ (7.2) \qquad \rho(0, x_1) = 1.0 + \delta \sin(\kappa\pi x_1), $$

where the frequency $\kappa = 2220$. The value of $\kappa$ is chosen such that $\kappa \sim \varepsilon^{-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 a stable equilibrium. The expectation from a numerical scheme is to be able to damp out the oscillations so that the electric field converges in time towards zero, i.e. the equilibrium is recovered in the long-time asymptotic.

Here we test both the AP and the classical schemes on a coarse mesh of 100 mesh points i.e. the mesh size $\Delta x$ not resolving $\varepsilon = 10^{-4}$ the Debye length. The AP scheme is subjected only to the advective CFL condition (6.21). But to adhere to the extra stability requirements of classical schemes, in addition to (6.21) the classical scheme is also subjected to the stability criteria (6.22). Figure 10 shows the density and the electric potential profiles for the AP scheme and the classical scheme at time $t = 0.1$. From this figure we conclude that the classical scheme even after being put under a far restrictive time stepping criteria is not able to converge asymptotically to
<figure><img src="image_6.png" /><figcaption>FIGURE 9. Small perturbation of a quasineutral state with Classical scheme ($\Delta t > \varepsilon$). $\varepsilon = 10^{-4}$, $N = 10^2$ (unresolved mesh). Electric Potential Profiles, at time $t = 0.0061$</figcaption></figure>

the equilibrium state on a coarse mesh. Whereas the AP scheme performs exceedingly well under
a very relaxed CFL condition to obtain a solution close to the quasineutral state.

<figure><img src="image_7.png" /><figcaption>FIGURE 10. Slight perturbation of a Maxwellian with AP and Classical schemes. $\varepsilon = 10^{-4}$, $N = 10^2$ (Unresolved mesh). Left:Density Profile and Right:Electric Potential Profile, at time = 0.1</figcaption></figure>

*Remark 7.2.* Owing to the Remark 1.1 in [25], we make the observation that even if the data is non-well prepared the AP scheme is able to damp out the oscillation on a coarse mesh, thereby leading towards a robust and a computationally efficient numerical method.

7.3. **Asymptotic Order of Convergence.** The goal of this experiment is to numerically validate the second order asymptotic convergence of the numerical solution to the incompressible limit solution. Let $v_{\varepsilon}(t^n)$ and $v_0(t^n)$ be the solution of the EP system (2.1)-(2.3) and the incompressible limit system (2.5)-(2.7), respectively, at time $t^n$. Let $v_{\varepsilon}^n$ be the numerical solution obtained using a
<table><caption>TABLE 1. AOC for SI-IMEX-RK-DAE schemes for $ε = 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>80</td><td>1.4641E-01</td><td></td><td>1.4816E-02</td><td></td><td>9.7052E-04</td><td></td></tr><tr><td>160</td><td>1.0083E-02</td><td>3.85</td><td>1.0301E-03</td><td>3.84</td><td>6.3394E-05</td><td>3.93</td></tr><tr><td>320</td><td>2.0090E-03</td><td>2.32</td><td>2.0869E-04</td><td>2.30</td><td>1.3126E-05</td><td>2.27</td></tr><tr><td>640</td><td>2.6880E-04</td><td>2.90</td><td>2.6880E-05</td><td>2.52</td><td>2.1804E-06</td><td>2.58</td></tr></tbody></table>

second order accurate AP scheme for the EP system (2.1)-(2.3) at the same time. Then we expect the following estimates:

$$ (7.3) \qquad v_{\varepsilon}^{n} - v_{0}(t^{n}) = v_{\varepsilon}^{n} - v_{\varepsilon}(t^{n}) + v_{\varepsilon}(t^{n}) - v_{0}(t^{n}) + v_{0}(t^{n}) - v_{0}(t^{n}) $$

$$ (7.4) \qquad = \mathcal{O}(\Delta t^2) + \mathcal{O}(\varepsilon) + \mathcal{O}(\Delta t^2). $$

Therefore from the above we get for $ε → 0$ we have:

$$ (7.5) \qquad v_{\varepsilon}^{n} - v_{0}(t^{n}) = \mathcal{O}(\Delta t^2) $$

Consequently, we compute the asymptotic order of convergence, as is defined in [1], as the experimental rate of convergence of a numerical scheme with respect to the asymptotic limit solution. Therefore the exact solution, used as a reference, is the solution of the incompressible limit (2.5)-(2.7).

To this end we take the initial conditions as considered in Section 7.1 with a modification in the simulation domain and the value of the perturbation parameter. Here we expand the plasma domain so that the initial perturbation can travel within the domain before encountering the boundary conditions at the domain walls. Consequently, the spatial domain is set to be the interval $[0, 10]$, and the perturbation parameter is increased to $δ = 10^{-2}$.

Here we have considered three different values of $ε ∈ \{10^{-4}, 10^{-5}, 10^{-6}\}$, all in the asymptotic range. Table 1 displays the second order convergence of the SI-IMEX-RK-DAE scheme in the quasineutral regime. From Table 1 we observe that even though the error is stuck around $\mathcal{O}(ε)$ still the AP scheme is able to recover the second order convergence.

## 8. CONCLUDING REMARKS

In this paper, we have presented a novel procedure which combines two semi-implicit time discretisation approaches, namely the additively partitioned IMEX-RK schemes and the SI-IMEX-RK schemes combined with the direct approach for implicit RK schemes for DAEs, in-order to develop high-order time semi-discrete schemes for the one-fluid Euler-Poisson system. The design is leveraged on the fact that via the method of lines, a linearised version of the governing equations can be perceived as a system of DAEs. Subsequently, we prove that the time semi-discrete SI-IMEX-RK-DAE schemes, which are developed on the SI-IMEX-RK platform are asymptotically consistent with the quasineutral limiting system. Whereas the schemes based on the additively partitioned IMEX platform fail to be AP, which is not surprising owing to the analysis of a first order scheme in [25]. Furthermore, we derive a space-time fully discrete scheme by approximating the spatial derivatives in a finite volume framework. Therein we chose a combination of Rusanov-type flux and central fluxes for the deemed stiff and the non-stiff terms, respectively. We show that the space discretisation complements the time semi-discrete scheme by preserving the AP property at the fully discrete level. Finally, the results of numerical case studies are provided to showcase the second order accuracy and the performance of the AP scheme in the quasineutral and non-quasineutral regimes.
The numerical case studies convey that the AP schemes perform equally well while resolving or not resolving the Debye length. Whereas the classical schemes even being put under much stricter stability criterion fail to replicate the fine mesh results on a coarse mesh. The SI-IMEX-RK scheme clearly demonstrates the second-order convergence of the numerical solution to the quasineutral limit solution.

Additionally, we make a note that the schemes developed here are designed for the EP system, on the account of its simplicity to aid the analysis and implementation of the same. But, the scheme can be extended to any stiff system of equations with hyperbolic and elliptic coupling, in particular for the two-fluid EP system [23], and the work is under process and will be reported in the future.

## REFERENCES

[1] K. R. Arun and S. Samantaray. Asymptotic preserving low Mach number accurate IMEX finite volume schemes for the isentropic Euler equations. *J. Sci. Comput.*, 82(2):Art. 35, 32, 2020.

[2] U. M. Ascher, S. J. Ruuth, and R. J. Spiteri. Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. *Appl. Numer. Math.*, 25(2-3):151-167, 1997. Special issue on time integration (Amsterdam, 1996).

[3] U. M. Ascher, S. J. Ruuth, and B. T. R. Wetton. Implicit-explicit methods for time-dependent partial differential equations. *SIAM J. Numer. Anal.*, 32(3):797-823, 1995.

[4] M. Bézard. Existence locale de solutions pour les équations d'Euler-Poisson. *Japan J. Indust. Appl. Math.*, 10(3):431-450, 1993.

[5] C. Birdsall and A. Langdon. *Plasma physics via computer simulation, Inst. of Phys*, volume 32. 1991.

[6] G. Bispen, K. R. Arun, M. Lukáčová-Medvidová, and S. Noelle. IMEX large time step finite volume methods for low Froude number shallow water flows. *Commun. Comput. Phys.*, 16(2):307-347, 2014.

[7] G. Bispen, M. Lukáčová-Medvidová, and L. Yelash. Asymptotic preserving IMEX finite volume schemes for low Mach number Euler equations with gravitation. *J. Comput. Phys.*, 335:222-248, 2017.

[8] S. Boscarino. Error analysis of IMEX Runge-Kutta methods derived from differential-algebraic systems. *SIAM J. Numer. Anal.*, 45(4):1600-1621, 2007.

[9] S. Boscarino, F. Filbet, and G. Russo. High order semi-implicit schemes for time dependent partial differential equations. *J. Sci. Comput.*, 68(3):975-1001, 2016.

[10] S. Boscarino, J.-M. Qiu, G. Russo, and T. Xiong. A high order semi-implicit IMEX WENO scheme for the all-Mach isentropic Euler system. *J. Comput. Phys.*, 392:594-618, 2019.

[11] S. Boscarino, G. Russo, and L. Scandurra. All Mach number second order semi-implicit scheme for the Euler equations of gas dynamics. *J. Sci. Comput.*, 77(2):850-884, 2018.

[12] J. U. Brackbill and D. W. Forslund. An implicit method for electromagnetic plasma simulation in two dimensions. *J. Comput. Phys.*, 46(2):271-308, 1982.

[13] U. Brauer and L. Karp. Local existence of solutions to the Euler-Poisson system, including densities without compact support. *J. Differential Equations*, 264(2):755-785, 2018.

[14] K. E. Brenan, S. L. Campbell, and L. R. Petzold. *Numerical solution of initial-value problems in differential-algebraic equations*, volume 14 of *Classics in Applied Mathematics*. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1996. Revised and corrected reprint of the 1989 original.

[15] Y. Brenier. Convergence of the Vlasov-Poisson system to the incompressible Euler equations. *Comm. Partial Differential Equations*, 25(3-4):737-754, 2000.

[16] F. Chen. *Introduction to Plasma Physics and Controlled Fusion: Volume 1: Plasma Physics*. Springer US, 2013.

[17] B. L. Cohen, A. B. Langdon, and A. Friedman. Implicit time integration for plasma simulation. *J. Comput. Phys.*, 46(1):15-38, 1982.

[18] P. Colella, M. R. Dorr, and D. D. Wake. A conservative finite difference method for the numerical solution of plasma fluid equations. *J. Comput. Phys.*, 149(1):168-193, 1999.

[19] F. Cordier, P. Degond, and A. Kumbaro. An asymptotic-preserving all-speed scheme for the Euler and Navier-Stokes equations. *J. Comput. Phys.*, 231(17):5685-5704, 2012.

[20] S. Cordier and E. Grenier. Quasineutral limit of an Euler-Poisson system arising from plasma physics. *Comm. Partial Differential Equations*, 25(5-6):1099-1113, 2000.

[21] S. Cordier and Y.-J. Peng. Système Euler-Poisson non linéaire. Existence globale de solutions faibles entropiques. *RAIRO Modél. Math. Anal. Numér.*, 32(1):1-23, 1998.
[22] P. Crispel, P. Degond, and M.-H. Vignal. An asymptotically stable discretization for the Euler-Poisson system in the quasi-neutral limit. *C. R. Math. Acad. Sci. Paris*, 341(5):323-328, 2005.

[23] P. Crispel, P. Degond, and M.-H. Vignal. An asymptotic preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit. *J. Comput. Phys.*, 223(1):208-234, 2007.

[24] N. Crouseilles, G. Dimarco, and M.-H. Vignal. 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):65-95, 2016.

[25] P. Degond. Asymptotic-preserving schemes for fluid models of plasmas. In *Numerical models for fusion*, volume 39/40 of *Panor. Synthèses*, pages 1-90. Soc. Math. France, Paris, 2013.

[26] P. Degond, J.-G. Liu, and M.-H. Vignal. Analysis of an asymptotic preserving scheme for the Euler-Poisson system in the quasineutral limit. *SIAM J. Numer. Anal.*, 46(3):1298-1322, 2008.

[27] P. Degond and M. Tang. All speed scheme for the low Mach number limit of the isentropic Euler equations. *Commun. Comput. Phys.*, 10(1):1-31, 2011.

[28] G. Dimarco, R. Loubère, and M.-H. Vignal. Study of a new asymptotic preserving scheme for the Euler system in the low Mach number limit. *SIAM J. Sci. Comput.*, 39(5):A2099-A2128, 2017.

[29] G. Dimarco and L. Pareschi. Asymptotic preserving implicit-explicit Runge-Kutta methods for nonlinear kinetic equations. *SIAM J. Numer. Anal.*, 51(2):1064-1087, 2013.

[30] D. Donatelli and P. Marcati. A quasineutral type limit for the Navier-Stokes-Poisson system with large data. *Nonlinearity*, 21(1):135-148, 2008.

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

[32] F. Filbet and L. M. Rodrigues. Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field. *SIAM J. Numer. Anal.*, 54(2):1120-1146, 2016.

[33] F. Filbet and L. M. Rodrigues. Asymptotically preserving particle-in-cell methods for inhomogeneous strongly magnetized plasmas. *SIAM J. Numer. Anal.*, 55(5):2416-2443, 2017.

[34] F. Filbet, L. M. Rodrigues, and H. Zakerzadeh. Convergence analysis of asymptotic preserving schemes for strongly magnetized plasmas. *Numer. Math.*, 149(3):549-593, 2021.

[35] P. Gamblin. Solution régulière à temps petit pour l'équation d'Euler-Poisson. *Comm. Partial Differential Equations*, 18(5-6):731-745, 1993.

[36] I. Gasser, C. D. Levermore, P. A. Markowich, and C. Schmeiser. The initial time layer problem and the quasineutral limit in the semiconductor drift-diffusion model. *European J. Appl. Math.*, 12(4):497-512, 2001.

[37] D. Gérard-Varet, D. Han-Kwan, and F. Rousset. Quasineutral limit of the Euler-Poisson system for ions in a domain with boundaries. *Indiana Univ. Math. J.*, 62(2):359-402, 2013.

[38] H. Goedbloed, R. Keppens, and S. Poedts. *Magnetohydrodynamics: Of Laboratory and Astrophysical Plasmas*. Cambridge University Press, 2019.

[39] J. Haack, S. Jin, and J.-G. Liu. An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations. *Commun. Comput. Phys.*, 12(4):955-980, 2012.

[40] E. Hairer and G. Wanner. *Solving ordinary differential equations. II*, volume 14 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition, 1996. Stiff and differential-algebraic problems.

[41] R. Hazeltine and J. Meiss. *Plasma Confinement*. Dover Books on Physics. Dover Publications, 2013.

[42] D. W. Hewett and C. W. Nielson. A multidimensional quasineutral plasma simulation model. *J. Comput. Phys.*, 29(2):219-236, 1978.

[43] I. Higueras, J. M. Mantas, and T. Roldán. Design and implementation of predictors for additive semi-implicit Runge-Kutta methods. *SIAM J. Sci. Comput.*, 31(3):2131-2150, 2009.

[44] F. Huang and Y. Li. Large time behavior and quasineutral limit of solutions to a bipolar hydrodynamic model with large data and vacuum. *Discrete Contin. Dyn. Syst.*, 24(2):455-470, 2009.

[45] S. Jin. Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations. *SIAM J. Sci. Comput.*, 21(2):441-454, 1999.

[46] S. Jin. Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review. *Riv. Math. Univ. Parma (N.S.)*, 3(2):177-216, 2012.

[47] Q. Ju, F. Li, and H. Li. The quasineutral limit of compressible Navier-Stokes-Poisson system with heat conductivity and general initial data. *J. Differential Equations*, 247(1):203-224, 2009.

[48] A. Jüngel. *Quasi-hydrodynamic semiconductor equations*, volume 41 of *Progress in Nonlinear Differential Equations and their Applications*. Birkhäuser Verlag, Basel, 2001.

[49] A. Jüngel and I. Violet. The quasineutral limit in the quantum drift-diffusion equations. *Asymptot. Anal.*, 53(3):139-157, 2007.
[50] R. Klein. Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics. I. One-dimensional flow. *J. Comput. Phys.*, 121(2):213-237, 1995.

[51] R. Klein, N. Botta, T. Schneider, C. D. Munz, S. Roller, A. Meister, L. Hoffmann, and T. Sonar. Asymptotic adaptive methods for multi-scale problems in fluid mechanics. *J. Engrg. Math.*, 39(1-4):261-343, 2001. Special issue on practical asymptotics.

[52] N. Krall and A. Trivelpiece. *Principles of Plasma Physics*. International series in pure and applied physics. San Francisco Press, 1986.

[53] A. B. Langdon, B. I. Cohen, and A. Friedman. Direct implicit large time-step particle simulation of plasmas. *J. Comput. Phys.*, 51(1):107-138, 1983.

[54] T. Makino. On a local existence theorem for the evolution equation of gaseous stars. In *Patterns and waves*, volume 18 of *Stud. Math. Appl.*, pages 459-479. North-Holland, Amsterdam, 1986.

[55] T. Makino. Blowing up solutions of the Euler-Poisson equation for the evolution of gaseous stars. In *Proceedings of the Fourth International Workshop on Mathematical Aspects of Fluid and Plasma Dynamics (Kyoto, 1991)*, volume 21, pages 615-624, 1992.

[56] T. Makino. Recent progress of the study of the Euler-Poisson equation for the evolution of gaseous stars. Number 824, pages 151-161. 1993. Mathematical analysis of phenomena in fluid and plasma dynamics (Japanese) (Kyoto, 1992).

[57] T. Makino and B. Perthame. Sur les solutions à symétrie sphérique de l'équation d'Euler-Poisson pour l'évolution d'étoiles gazeuses. *Japan J. Appl. Math.*, 7(1):165-170, 1990.

[58] T. Makino and S. Ukai. Sur l'existence des solutions locales de l'équation d'Euler-Poisson pour l'évolution d'étoiles gazeuses. *J. Math. Kyoto Univ.*, 27(3):387-399, 1987.

[59] P. Marcati and R. Natalini. Weak solutions to a hydrodynamic model for semiconductors and relaxation to the drift-diffusion equation. *Arch. Rational Mech. Anal.*, 129(2):129-145, 1995.

[60] R. J. Mason. Implicit moment particle simulation of plasmas. *J. Comput. Phys.*, 41(2):233-244, 1981.

[61] R. J. Mason. Implicit moment pic-hybrid simulation of collisional plasmas. *J. Comput. Phys.*, 51(3):484-501, 1983.

[62] R. J. Mason. Hybrid and collisional implicit plasma simulation models. In *Multiple time scales*, volume 3 of *Comput. Tech.*, pages 233-270. Academic Press, Orlando, FL, 1985.

[63] R. J. Mason. An electromagnetic field algorithm for 2D implicit plasma simulation. *J. Comput. Phys.*, 71(2):429-473, 1987.

[64] C.-D. Munz, S. Roller, R. Klein, and K. J. Geratz. The extension of incompressible flow solvers to the weakly compressible regime. *Comput. & Fluids*, 32(2):173-196, 2003.

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

[66] S. Noelle, G. Bispen, K. R. Arun, M. Lukáčová-Medvidová, and C.-D. Munz. A weakly asymptotic preserving low Mach number scheme for the Euler equations of gas dynamics. *SIAM J. Sci. Comput.*, 36(6):B989-B1024, 2014.

[67] L. Pareschi and G. Russo. Implicit-explicit Runge-Kutta schemes for stiff systems of differential equations. In *Recent trends in numerical analysis*, volume 3 of *Adv. Theory Comput. Math.*, pages 269-288. Nova Sci. Publ., Huntington, NY, 2001.

[68] L. Pareschi and G. Russo. Implicit-Explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation. *J. Sci. Comput.*, 25(1-2):129-155, 2005.

[69] F. Poupaud, M. Rascle, and J.-P. Vila. Global solutions to the isothermal Euler-Poisson system with arbitrarily large data. *J. Differential Equations*, 123(1):93-121, 1995.

[70] P. W. Rambo and J. Denavit. Fluid and field algorithms for time-implicit plasma simulation. *J. Comput. Phys.*, 92(1):185-212, 1991.

[71] P. W. Rambo and J. Denavit. Time-implicit fluid simulation of collisional plasmas. *J. Comput. Phys.*, 98(2):317-331, 1992.

[72] R. Schneider and C.-D. Munz. The approximation of two-fluid plasma flow with explicit upwind schemes. *Int. J. Numer. Model. 8*, page 399, 2005.

[73] T. Schneider, N. Botta, K. J. Geratz, and R. Klein. Extension of finite volume compressible flow solvers to multi-dimensional, variable density zero Mach number flows. *J. Comput. Phys.*, 155(2):248-286, 1999.

[74] U. Shumlak and J. Loverich. Approximate riemann solver for the two-fluid plasma model. *J. Comput. Phys.*, 187:620, 2003.

[75] M. Tang. Second order all speed method for the isentropic Euler equations. *Kinet. Relat. Models*, 5(1):155-184, 2012.
[76] J. M. Wallace, J. U. Brackbill, and D. W. Forslund. An implicit moment electromagnetic plasma simulation in cylindrical coordinates. *J. Comput. Phys.*, 63:434, 1986.

[77] S. Wang. Quasineutral limit of Euler-Poisson system with and without viscosity. *Comm. Partial Differential Equations*, 29(3-4):419–456, 2004.

[78] X. Zhong. Additive semi-implicit Runge-Kutta methods for computing high-speed nonequilibrium reactive flows. *J. Comput. Phys.*, 128(1):19–31, 1996.

SCHOOL OF MATHEMATICS, INDIAN INSTITUTE OF SCIENCE EDUCATION AND RESEARCH THIRUVANANTHAPURAM,
THIRUVANANTHAPURAM - 695551, INDIA.

Email address: arun@iisertvm.ac.in

UNIV RENNES, INRIA BRETAGNE ATLANTIQUE (MINGUS) & IRMAR UMR 6625 & ENS RENNES, FRANCE

Email address: nicolas.crouseilles@inria.fr

CENTRE FOR APPLICABLE MATHEMATICS, TATA INSTITUTE OF FUNDAMENTAL RESEARCH, BANGALORE - 560065,
INDIA.

Email address: saurav22@tifrbng.res.in