Alexandre Ern

Jean-Luc Guermond

Invariant-domain preserving high-order time stepping: II. IMEX schemes *

Keywords: Time integration, implicit-Explicit time integration methods, conservation equations, hyperbolic systems, invariant-domains, high-order method. MSC. 35L65, 65M60, 65M12, 65N30

We consider high-order discretizations of a Cauchy problem where the evolution operator comprises a hyperbolic part and a parabolic part with diffusion and stiff relaxation terms. We propose a technique that makes every implicit-explicit (IMEX) time stepping scheme invariantdomain preserving and mass conservative. Following the ideas introduced in Part I on explicit Runge-Kutta schemes, the IMEX scheme is written in incremental form. At each stage, we first combine a low-order and a high-order hyperbolic update using a limiting operator, then we combine a low-order and a high-order parabolic update using another limiting operator. The proposed technique, which is agnostic to the space discretization, allows to optimize the time step restrictions induced by the hyperbolic sub-step. To illustrate the proposed methodology, we derive four novel IMEX methods with optimal efficiency. All the implicit schemes are singly diagonal. One of them is A-stable and the other three are L-stable. The novel IMEX schemes are evaluated numerically on a stiff ODE system and a scalar nonlinear conservation equation.

Introduction

This work is the second part of a project started in [START_REF] Ern | Invariant-domain-preserving high-order time stepping: I. explicit runge-kutta schemes[END_REF] whose objective is to develop Runge-Kutta time stepping schemes that are invariant-domain preserving (IDP) and conservative. The scope of the present work lies in the approximation of the following Cauchy problem posed on the space domain D ⊂ R d and the time interval J := (0, T ) with T > 0:

∂ t u + f (u) + g(u, ∇u) = 0 in D × J, u(0) = u 0 in D, (1)
supplemented with appropriate boundary conditions. The dependent variable u takes values in R m with m ≥ 1. The operator f : A → R m represents the hyperbolic part of the problem, and the operator g : A × R m×d → R m represents the parabolic part, typically associated with diffusion and (stiff) relaxation processes. Here, A is the domain of f and A × R m×d is the domain of g. In the applications we have in mind, these operators have the following structure:

f (u) = ∇•f(u), g(u, ∇u) = ∇•d(u, ∇u) + r(u), (2)
with the hyperbolic flux f : A → R m×d , the diffusive flux d : A × R m×d → R m×d , and the relaxation operator r : A → R m . As it is out of the scope of this paper to discuss the existence and uniqueness of solutions to (1), we assume that this problem admits a reasonable class of solutions. We also assume that the set A ⊂ R m is an invariant domain for this solution class. This means that if u 0 (x) ∈ A for a.e.

x in D (and up to perturbations resulting from boundary conditions which go beyond the scope of this paper), then any admissible solution to [START_REF] Ascher | Implicit-explicit methods for time-dependent partial differential equations[END_REF] takes values in A at a.e. x in D at a.e. time t ∈ [0, T ]. The set A may depend on u 0 . A simple example is the scalar convection-diffusion equation (i.e., m = 1), in which case the interval A := [ess inf x∈D u 0 (x), ess sup x∈D u 0 (x)] ⊂ R is an invariant domain. Two more elaborate examples are the compressible Euler equations and the compressible Navier-Stokes equations. For these equations, the conserved variable u takes values in R d+2 and its components are the density, the momentum, and the total mechanical energy (i.e., m = d + 2). An invariant domain for the compressible Euler equations is the set A composed of those states with positive density, positive internal energy, and specific entropy s(u) larger than ess inf x∈D s(u 0 (x)). An invariant domain for the compressible Navier-Stokes equations is the set A composed of those states with positive density and positive internal energy. Another important property of (1) is conservation. Letting (r i ) i∈{1:m} be the components of r, we assume that there is an index subset C ⊂ {1:m} := {1, . . . , m} such that r p (u) = 0 for all p ∈ C. This means that the relaxation process does not affect the dependent variables indexed in the subset C. Then, again in the absence of perturbations due to the boundary conditions, the following conservation property holds:

D u p (t, •) dx = D u 0,p dx, ∀t ∈ J, ∀p ∈ C. (3)
The objective of this work is to construct high-order discretizations in space and time that are conservative and leave the set A invariant. Such methods are called invariant-domain preserving for A, or (IDP) for short. To stay general, our starting point is a system of ordinary differential equations (ODEs) obtained after discretization in space of the conservation equation [START_REF] Ascher | Implicit-explicit methods for time-dependent partial differential equations[END_REF]. We mainly focus in this paper on the time discretization. The time discretization methods we are going to present can be combined with various space discretization techniques (e.g., discontinuous and continuous finite elements, finite volumes, finite differences, etc.). We assume that the ODE system takes the following generic form:

M∂ t U = F(U) + G(U), ∀t ∈ J, U(0) = U 0 . (4)
The mass matrix M is induced by the space discretization (it is the Gram matrix associated with the global shape functions of the space approximation). The dependent variable U(t) takes values in (R m ) I where I ≥ 1 is the number of degrees of freedom (dofs) employed in the space discretization. We set V := {1:I} := {1, . . . , I} and write U(t) := (U i (t)) i∈V . For all i ∈ V, the local state vector U i (t) = (U p,i (t)) p∈{1:m} is viewed as an approximation of the exact solution u(t, •) at some point in D, say x i . The nonlinear mappings F ∈ C 0 (A I ; (R m ) I ) and G ∈ C 0 (A I ; (R m ) I ) result from the space discretization of the operators -f and -g in [START_REF] Ascher | Implicit-explicit methods for time-dependent partial differential equations[END_REF], respectively, and U 0 ∈ A I is an appropriate approximation in space of the initial datum u 0 . We loosely refer to the mappings F and G as the hyperbolic flux and the parabolic flux, respectively. Assuming that U i (0) ∈ A for all i ∈ V, a natural requirement for the space approximation is that it is invariant-domain preserving, i.e., U(t) ∈ A I , ∀t ∈ J.

A second requirement is that conservation holds:

i∈V m i U p,i (t) = i∈V m i U p,i (0), ∀t ∈ J, ∀p ∈ C. (6)
Using an implicit scheme to discretize in time the ODE system (4) is often too expensive owing to the nonlinearities involved in the fluxes F and G, whereas using an explicit scheme results in a severe restriction on the time step owing to stiffness induced by the parabolic flux G. A well known remedy to this conundrum is to resort to implicit-explicit (IMEX) schemes, where the numerical flux F is treated explicitly, and the numerical flux G is treated implicitly. The origin of IMEX schemes can be traced back to Crouzeix [START_REF] Crouzeix | Une méthode multipas implicite-explicite pour l'approximation des équations d'évolution paraboliques[END_REF], Varah [START_REF] Varah | Stability restrictions on second order, three level finite difference schemes for parabolic equations[END_REF]. We also refer the reader to Ascher et al. [START_REF] Ascher | Implicit-explicit methods for time-dependent partial differential equations[END_REF][START_REF] Ascher | Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations[END_REF], Burman and Ern [START_REF] Burman | Implicit-explicit Runge-Kutta schemes and finite elements with symmetric stabilization for advection-diffusion equations[END_REF], Kennedy and Carpenter [START_REF] Kennedy | Additive Runge-Kutta schemes for convection-diffusionreaction equations[END_REF], Pareschi and Russo [START_REF] Pareschi | Implicit-explicit Runge-Kutta schemes for stiff systems of differential equations[END_REF][START_REF] Pareschi | Implicit-Explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation[END_REF], Zhong [START_REF] Zhong | Additive semi-implicit Runge-Kutta methods for computing high-speed nonequilibrium reactive flows[END_REF] for other developments. Despite these advances, a crucial question that still remains open is how to reconcile the use of an IMEX time stepping scheme with the above invariant-domain property, while at the same time ensuring conservation. Building on [START_REF] Ern | Invariant-domain-preserving high-order time stepping: I. explicit runge-kutta schemes[END_REF], we propose an answer to this question in this paper. More precisely, we introduce a technique that makes every IMEX Runge-Kutta (RK) time stepping method invariant-domain preserving (IDP) and conservative. The resulting schemes are called "IDP-IMEX" schemes.

This work is organized as follows. In Section 2, we outline the discrete setting in space and time, we identify the key assumptions underlying this work, and we exemplify these notions for the Euler IDP-IMEX scheme. In Section 3, we extend these ideas and build higher-order IDP-IMEX schemes. We introduce a generic IDP-IMEX algorithm composed of the steps (51) to (61) whose properties are stated in Theorem 3.3. In Section 4, we review some examples of higher-order IMEX schemes and we derive some novel examples with optimal efficiency. Finally, in Section 5, we present numerical illustrations on stiff ODEs and a nonlinear scalar conservation equation.

Preliminaries

The goal of this section is threefold: (i) introduce the discrete setting in space and time; (ii) identify the key ideas and assumptions underlying this work; (iii) exemplify these notions for the Euler IDP-IMEX scheme. All this material is used in Section 3 where we introduce the novel higher-order IDP-IMEX schemes.

Time discretization and quasi-linearization

Let t n ∈ [0, T ] be the current time with n ∈ {0:N }, t 0 := 0, and t N := T . Let τ n be the current time step and let t n+1 := t n +τ n . To simplify the notation, we henceforth write τ instead of τ n . Let U n be the approximation of the solution to (4) at the discrete time t n . The key invariant-domain property we want to achieve is the following:

(U n ∈ A I ) =⇒ (U n+1 ∈ A I ), ∀n ≥ 0. (7)
Moreover, we want to achieve this goal while maintaining conservation. The notion of conservation will be made more precise below, but for the time being, conservation is expressed at the global level by requiring that

i∈V m i U n+1 p,i = i∈V m i U n p,i , ∀p ∈ C. (8)
To avoid solving a nonlinear problem at each time-step when the parabolic fluxes are made implicit, we introduce a quasi-linearization process (the way this is done is made more precise in §2.2 and §2.3). We consider a quasi-linearized parabolic flux

G lin ∈ C 0 (A I ×(R m ) I ; (R m ) I ) that is consistent with G, i.e., G lin (W; W) = G(W), ∀W ∈ A I . (9)
We assume that this flux is such that for all W ∈ A I , the problem consisting of seeking

U ∈ (R m ) I so that MU -τ G lin (W; U) = MW (10)
is well-posed and easy to solve. For instance, this problem could only involve linear solves. Notice that this does not mean that the mapping U → G(W; U) is linear; see §5.1 for an example. Owing to the above quasi-linearization process, we reformulate (4) over the time interval

J n := [t n , t n+1 ] as follows: Find U ∈ C 1 (J n ; (R m ) I ) so that U(t n ) = U n and for all t ∈ J n , M∂ t U = F(U) + G(U) -G lin (U n ; U) explicit + G lin (U n ; U) implicit . ( 11
)

Space discretization and conservation structure

Let us now give details on the space discretization. We consider two space discretizations. The first one is low-order accurate and referred to with the superscript L . The second one is highorder accurate and referred to with the superscript H . The low-order scheme is based on a loworder invertible mass matrix M L ∈ R I×I and low-order fluxes F L , G L : A I → (R m ) I . The highorder scheme is based on a high-order invertible mass matrix M H ∈ R I×I and high-order fluxes

F H , G H : A I → (R m ) I .
We assume that M H is symmetric positive-definite with entries (m ij ) i,j∈V and M L is diagonal with entries (δ ij m i ) i,j∈V . For all i ∈ V, we introduce the subset I(i) V such that m ij = 0 for all j ∈ I(i). We call I(i) stencil at i. The notion of stencil is symmetric, i.e., j ∈ I(i) if and only if i ∈ I(j) because m ij = m ji . Finally, we assume that

m i = j∈V m ij = j∈V m ji , ∀i ∈ V, (12)
In the finite element terminology, this means that the low-order mass matrix M L is the lumped version of the high-order mass matrix M H . For every matrix M ∈ R I×I and every vector V ∈ (R m ) I with components V p,i , with p ∈ {1:m} and i ∈ V, the components of the vector MV ∈ (R m ) I are defined to be (MV) p,i := j∈V m ij V p,j for all p ∈ {1:m} and i ∈ V.

The components of the low-order and high-order hyperbolic fluxes are denoted

F L i (V) ∈ R m and F H i (V) ∈ R m
for all i ∈ V and all V ∈ A I . To account for the fact that the hyperbolic fluxes are associated with a conservation principle, we assume that these fluxes admit the following stencil-based decomposition:

F L i (V) = j∈I(i) F L ij (V), F H i (V) = j∈I(i) F H ij (V), ∀V ∈ A I , (13)
where

F L ij , F H ij ∈ C 0 (A I ; R m
), and we assume the following skew-symmetry property:

F L ij (V) = -F L ji (V), F H ij (V) = -F H ji (V), ∀i ∈ V, ∀j ∈ I(i). (14)
The same structure is assumed for the parabolic fluxes, namely (for brevity, we only write the statements for the high-order fluxes)

G H i (V) = j∈I(i) D H ij (V) + R H i (V), D H ij (V) = -D H ji (V), ∀i ∈ V, ∀j ∈ I(i). (15)
Consistently with our assumption that r p (u) = 0 for all p ∈ C, we assume that R H p,i (V) = 0 for all p ∈ C and all i ∈ V.

The quasi-linearization process mentioned in [START_REF] Guermond | Invariant domains and first-order continuous finite element approximation for hyperbolic systems[END_REF] is performed for both the low-order and high-order parabolic fluxes. This leads to quasi-linearized parabolic fluxes G L,lin , G H,lin ∈ C 0 (A I × (R m ) I ; (R m ) I ) which we assume satisfy the following decompositions and properties:

G L,lin i (W; V) = j∈I(i) D L,lin ij (W; V) + R L,lin i (W; V), (16)
G H,lin i (W; V) = j∈I(i) D H,lin ij (W; V) + R H,lin i (W; V), (17)
D L,lin ij (W; V) = -D L,lin ji (W; V), D H,lin ij (W; V) = -D H,lin ji (W; V), (18)
R L,lin p,i = R H,lin p,i = 0, ∀p ∈ C, ∀i ∈ V. (19)
In conclusion, we are going to consider two versions of the ODE system (11) over the time interval

J n = [t n , t n+1 ].
One corresponds to the low-order space discretization:

M L ∂ t U L = F L (U L ) explicit + G L,lin (U n ; U L ) implicit . ( 20
)
The other one corresponds to the high-order space discretization:

M H ∂ t U H = F H (U H ) + G H (U H ) -G H,lin (U n ; U H ) explicit + G H,lin (U n ; U H ) implicit . (21)
Consistently with our assumption on [START_REF] Ern | Invariant-domain-preserving high-order time stepping: I. explicit runge-kutta schemes[END_REF], we assume that for all V ∈ (R m ) I and all W ∈ A I , the problems consisting of seeking U L , U H ∈ (R m ) I so that

M L U L -τ G L,lin (W; U L ) = M L V, M H U H -τ G H,lin (W; U H ) = M H V, (22)
are well-posed and easy to solve.

Structural IDP assumptions

To gently introduce our ideas, let us consider the well-known IMEX method consisting of combining the forward and the backward Euler time steppings. We call this method Euler IMEX. Let us apply the method to the low-order ODE system [START_REF] Kuzmin | Flux-Corrected Transport: Principles, Algorithms, and Applications[END_REF]. Let U n ∈ A I . The first step is explicit and consists of computing the hyperbolic prediction

W L,n := I + τ (M L ) -1 F L (U n ). (23)
The second step is implicit and consists of computing the final state U L,n+1 by solving the quasilinear problem

I -τ (M L ) -1 G L,lin (W L,n ; •) (U L,n+1 ) = W L,n . (24)
Altogether, we have

M L U L,n+1 = M L U n + τ F L (U n ) + τ G L,lin (W L,n ; U L,n+1 ). (25)
This leads us to formulate the following two key structural assumptions on the low-order fluxes:

Assumption 2.1 (Low-order fluxes). There exists τ * > 0 s.t. for all τ ∈ (0, τ * ],

(i) the low-order hyperbolic flux satisfies the following property:

V ∈ A I =⇒ I + τ (M L ) -1 F L (V) ∈ A I . (26)
(ii) the low-order quasi-linearized parabolic flux satisfies the following property:

V ∈ A I =⇒ I -τ (M L ) -1 G L,lin (V; •) -1 (V) ∈ A I . ( 27
)
The following result prefigures what we are aiming at. We omit the proof since it is somewhat standard.

Lemma 2.2 (Low-order Euler IDP-IMEX scheme). Assume that U n ∈ A I and τ ∈ (0, τ * ]. Then, the low-order Euler IMEX scheme (25) is well-defined. Moreover, it is IDP and conservative, i.e., U L,n+1 ∈ A I and i∈V m i U L,n+1

p,i = i∈V m i U n p,i for all p ∈ C.
Remark 2.3 (Time step restriction). In many situations, the time step restriction τ ∈ (0, τ * ] is only required for the invariant-domain property of the hyperbolic step [START_REF] Sanders | A third-order accurate variation nonexpansive difference scheme for single nonlinear conservation laws[END_REF]. The invariant-domain property of the parabolic step [START_REF] Shu | Efficient implementation of essentially non-oscillatory shockcapturing schemes[END_REF] can often be shown to hold for every time step τ > 0.

Of course, the above result is of little interest since what we actually want is to use a highorder approximation in space. The Euler IMEX scheme applied to the high-order ODE system (21) consists of seeking U H,n+1 ∈ (R m ) I so that

M H U H,n+1 = M H U n + τ F H (U n ) + τ G H,lin (U n ; U H,n+1 ). (28)
Similarly to [START_REF] Pareschi | Implicit-Explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation[END_REF], the method ( 28) is composed of two steps. The first one consists of computing the forward Euler prediction W H,n :=

I+τ (M H ) -1 F H )(U n ).
The second one consists of computing the parabolic update U H,n+1 by solving the quasi-linear problem

I-τ (M H ) -1 G H,lin (U n ; •) (U H,n+1 ) = W H,n .
Unfortunately, even though we are using forward and backward Euler time stepping, there is no guarantee that U H,n+1 belongs to A I , i.e., the high-order counterpart of Lemma 2.2 does not hold true in general. This problem is solved in the literature by using nonlinear limiting operators; see Boris and Book [START_REF] Boris | Flux-corrected transport. I. SHASTA, a fluid transport algorithm that works[END_REF], Harten [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF], Osher and Chakravarthy [START_REF] Osher | High resolution schemes and the entropy condition[END_REF], Zalesak [START_REF] Zalesak | Fully multidimensional flux-corrected transport algorithms for fluids[END_REF]. Limiting is realized in the discontinuous Galerkin and finite volume settings by squeezing the high-order approximation towards the piecewise constant approximation over each mesh cell (see Sanders [26,Thm. 2 [START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF]Thm. 2.5]). A wellknown limiting method for scalar conservation equations is the so-called flux-corrected transport (FCT) technique of Boris and Book [START_REF] Boris | Flux-corrected transport. I. SHASTA, a fluid transport algorithm that works[END_REF] and Zalesak [START_REF] Zalesak | Fully multidimensional flux-corrected transport algorithms for fluids[END_REF]. The reader is also referred to Kuzmin and Turek [START_REF] Kuzmin | Flux correction tools for finite elements[END_REF] and Kuzmin et al. [START_REF] Kuzmin | Flux-Corrected Transport: Principles, Algorithms, and Applications[END_REF] for other extensions on this method in the context of finite elements. When the bounds to be enforced are non-affine (as is often the case for hyperbolic systems), one has to use nonlinear methods like in, e.g., [START_REF] Sanders | A third-order accurate variation nonexpansive difference scheme for single nonlinear conservation laws[END_REF]Lem. 3.3] [START_REF] Guermond | Second-order invariant domain preserving approximation of the Euler equations using convex limiting[END_REF][START_REF] Guermond | Invariant domain preserving discretizationindependent schemes and convex limiting for hyperbolic systems[END_REF]). The key idea common to all the above techniques is the decomposition of the flux over the stencils in skew-symmetric components as in ( 13), [START_REF] Guermond | Second-order invariant domain preserving approximation of the compressible Navier-Stokes equations[END_REF], and [START_REF] Hairer | Solving Ordinary Differential Equations. II. Stiff and Differentialalgebraic Problems[END_REF].

In the present work, we are going to consider two limiters: one to compute the hyperbolic prediction and another to compute the parabolic update. We now introduce the corresponding notation. Let L be the collection of the sparse symmetric matrices with coefficients in [0, 1] and with the sparsity pattern induced by the stencils (I(i)) i∈V . We also let M be the collection of the skew-symmetric matrices in (R m ) I×I with the block-sparsity pattern induced by the stencils (I(i)) i∈V . Finally, we define

B := {z := (z 1 , . . . , z m ) ∈ R m | z p = 0, ∀p ∈ C}. Definition 2.4 (Conservative hyperbolic limiter). Let V i + τ mi j∈I(i) A ij i∈V be an hyperbolic prediction, with V := (V i ) i∈V ∈ A I and A := (A ij ) i∈V,j∈I(i) ∈ M.
We call conservative hyperbolic limiter any operator hyp : A I ×M (V, A) → ( ij ) i∈V,j∈I(i) ∈ L s.t. the following holds:

V i + τ m i j∈I(i) ij A ij ∈ A, ∀i ∈ V. ( 29
)
For brevity, the state

(V i + m -1 i τ j∈I(i) ij A ij ) i∈V ∈ A I is denoted hyp (V, A). Definition 2.5 (Conservative parabolic limiter). Let V i + τ mi j∈I(i) A ij + τ mi B i i∈V be a parabolic update with V := (V i ) i∈V ∈ A I , A := (A ij ) i∈V,j∈I(i) ∈ M, and B := (B i ) i∈V ∈ B I .
We call conservative parabolic limiter any operator par :

A I ×M×B I (V, A, B) → ( a ij ) i∈V,j∈I(i) , ( b i ) i∈V ∈ L×[0, 1] I s.t.

the following holds:

V i + τ m i j∈I(i) a ij A ij + τ m i b i B i ∈ A, ∀i ∈ V. ( 30
)
For brevity, the state

V i + τ mi j∈I(i) a ij A ij + τ mi b i B i is denoted par (V, A, B).
The existence of limiters is guaranteed since the trivial limiters hyp (V, A) = V (i.e., ij = 0 for all i ∈ V and all j ∈ I(i)) and par (V, A, B) = V (i.e., a ij = b i = 0 for all i ∈ V and all j ∈ I(i)) are always admissible because V ∈ A I . Of course, the trivial limiters are inefficient. The goal of limiters is to construct the limiting coefficients ij , a ij and b i as close to 1 as possible. Regardless of the values taken by the limiters, an important property of the limiters is conservativity. Lemma 2.6 (Conservation). For all (V, A, B) ∈ A I ×M×B I , all p ∈ C, we have

i∈V m i hyp (V, A) p,i = i∈V m i V p,i , i∈V m i par (V, A, B) p,i = i∈V m i V p,i . (31)
Proof. For the parabolic limiter, we have

i∈V m i par (V, A, B) p,i = i∈V m i V i +τ i∈V j∈I(i) a ij A ij for all p ∈ C because B ∈ B I .
But the symmetry property a ij = a ji and the skew-symmetry property A ij = -A ji imply that i∈V j∈I(i) a ij A ij = 0, whence the assertion. The proof for the hyperbolic limiter is similar.

Euler IDP-IMEX scheme

All the ingredients are now in place to make the Euler IMEX scheme invariant-domain preserving with high-order space discretization. Given U n ∈ A I , the scheme is decomposed into the following four steps:

U n (1) -→ (W L,n+1 , W H,n+1 ) (2)
-→

hyperbolic step W n+1 (3) -→ (U L,n+1 , U H,n+1 ) (4) -→ parabolic step U n+1 . ( 32
)
Let us now give the details of the four steps. We essentially follow Zalesak's limiting strategy [29, Eq. ( 4)] for both the hyperbolic and the parabolic steps.

Hyperbolic steps ( 1) and ( 2)

Step (1) consists of computing the low-order and high-order hyperbolic updates

M L W L,n+1 := M L U n + τ F L (U n ), (33)
M H W H,n+1 := M H U n + τ F H (U n ). (34)
In Step (2), we apply the hyperbolic limiting operator. Subtracting (33) from (34) and using ( 13) and ( 14), elementary manipulations show that for all i ∈ V,

W H,n+1 i = W L,n+1 i + τ m i j∈I(i) A n ij , (35)
with

A n ij := F H ij (U n ) -F L ij (U n ) + m i δ ij -m ij τ (W H,n+1 j -U n j -W H,n+1 i + U n i ). (36)
Notice that A n is indeed skew-symmetric. Then using Definition 2.4, the conservative IDP hyperbolic high-order update is obtained by setting

W n+1 := hyp (W L,n+1 , A n ). (37)
Parabolic steps ( 3) and ( 4)

Step (3) consists of computing the low-order and high-order parabolic updates by solving

M L U L,n+1 -τ G L,lin (W n+1 ; U L,n+1 ) := M L W n+1 , (38)
M H U H,n+1 -τ G H,lin (U n ; U H,n+1 ) := M H W n+1 . ( 39
)
The quasi-linearization in (38) is based on W n+1 to invoke Assumption [START_REF] Shu | Efficient implementation of essentially non-oscillatory shockcapturing schemes[END_REF]. The quasi-linearization in (39) is based on U n to be consistent with the higher-order case to be explained in the next section (see also Remark 3.1). In Step (4), we apply the parabolic limiting operator. Subtracting (38) from (39) and using ( 16), [START_REF] Kennedy | Additive Runge-Kutta schemes for convection-diffusionreaction equations[END_REF], and ( 18), elementary manipulations show that for all i ∈ V,

U H,n+1 i = U L,n+1 i + τ m i j∈I(i) A n ij + τ m i B n i , (40)
with

A n ij := D H,lin ij (U n ; U H,n+1 ) -D L,lin ij (W n+1 ; U L,n+1 ) + τ -1 (m i δ ij -m ij )(U H,n+1 j -W n+1 j -U H,n+1 i + W n+1 i ), (41)
B n i := R H,lin i (U n ; U H,n+1 ) -R L,lin i (W n+1 ; U L,n+1 ). (42)
Notice that A n is indeed skew-symmetric and B n ∈ B I . Then using Definition 2.5, the conservative IDP parabolic high-order update is obtained by setting

U n+1 := par (U L,n+1 , A n , B n ). (43)

Conclusion

The main result for the above construction is the following.

Lemma 2.7 (High-order Euler IDP-IMEX scheme). Let Assumption 2.1 hold and assume τ ∈ (0, τ * ]. Assume that the limiters match Definitions 2.4-2.5. Let U n ∈ A I . Then, the above Euler IMEX schemes (33)-( 43) satisfies the following properties:

(i) It is well defined;

(ii) It is IDP, i.e., it satisfies (7);

(iii) It is conservative, i.e., it satisfies (8).

Proof.

(1) Let us first prove the assertions (i) and (ii). Since U n ∈ A I , invoking Assumption 2.1(i) gives W L,n+1 ∈ A I . The definition of the hyperbolic limiter then implies that W n+1 ∈ A I . Invoking Assumption 2.1(ii) then shows that U L,n+1 is well-defined and U L,n+1 ∈ A I . Finally, the definition of the parabolic limiter implies that U n+1 ∈ A I .

(2) Conservativity follows from the following identities for all p ∈ C:

i∈V m i U n+1 p,i = i∈V m i U L,n+1 p,i = i∈V m i W n+1 p,i = i∈V m i W L,n+1 p,i = i∈V m i U n p,i ,
where the first and third equalities follow from Lemma 2.6, whereas the second and fourth equalities follow from the skew-symmetry assumption on the low-order fluxes.

3 High-order IDP-IMEX schemes

In this section, we extend the construction described in §2.4 to any IMEX scheme combining an explicit Runge-Kutta (ERK) scheme with a diagonally implicit RK scheme whose first stage is fully explicit (EDIRK). We assume that both schemes consist of s stages, s ≥ 2. The main original ideas of this paper are in this section.

Butcher tableaux

The ERK and EDIRK schemes are described by their respective Butcher tableau which we assume to have the following form:

0 0 c 2 a e 2,1 0 c 3 a e 3,1 a e 3,2 0 . . . . . . . . . . . . . . . c s a e s,1 a e s,2 • • • a e s,s-1 0 b e 1 b e 2 • • • b e s-1 b e s 0 0 c 2 a i 2,1 a i 2,2 c 3 a i 3,1 a i 3,2 a i 3,3 . . . . . . . . . . . . . . . c s a i s,1 a i s,2 • • • a i s,s-1 a i s,s b i 1 b i 2 • • • b i s-1 b i s (44)
Notice that a e l,k := 0 for all k ≥ l, a i l,k := 0 for all k > l. Here, the superscript e refers to the ERK scheme and the superscript i refers to the EDIRK scheme. Recall that the coefficients c j define the intermediate time steps t n,j := t n +c j τ . Notice that both schemes share the same set of coefficients (c j ) j∈{1:s} . For convenience, we set c s+1 := 1. We assume that c 1 = 0 and c j ≥ 0 for all j ∈ {2:s}. To simplify some expressions, we set a e s+1,j := b e j and a i s+1,j := b i j for all j ∈ {1:s}. Whenever a i i,1 = 0 for all i ∈ {1:s + 1}, the EDIRK scheme is called zero-padded. If all the diagonal entries a i i,i are equal (except the first one which is zero), we speak of singly diagonal EDIRK scheme (ESDIRK). We are going to assume that l∈{1:j} a e j,l = l∈{1:j}

a i j,l = c j , ∀j ∈ {1:s}. ( 45
)
This is one of Butcher's simplifying assumptions for each RK scheme. This assumption implies that a e 1,1 = a i 1,1 = 0. Moreover, since consistency requires that j∈{1:s} b e j = j∈{1:s} b i j = 1, the identity (45) also holds true for j = s + 1. Finally, we use the convention that a i s+1,s+1 := 0.

High-order IMEX scheme in incremental form

Following the ideas developed in [START_REF] Ern | Invariant-domain-preserving high-order time stepping: I. explicit runge-kutta schemes[END_REF] for ERK schemes, we write the IMEX scheme in incremental form. To this purpose, for all l ∈ {2:s + 1}, we define the stage index l (l) to be the largest index in {1:l -1} so that c l -c l (l) is the minimal value of c l -c k for all k ∈ {1:l -1} so that c l -c k ≥ 0:

l (l) := min{k ∈ {1:l -1} | c l -c k ≥ 0}, ∀l ∈ {2:s + 1}. ( 46
)
Owing to the assumption c l ≥ 0 = c 1 for all l ∈ {2:s+1}, we infer that

1 ∈ {k ∈ {1:l-1} | c l -c k ≥ 0}, which means that the set {k ∈ {1:l -1} | c l -c k ≥ 0}
is nonempty and the above definition makes sense. The definition of l (l) remains meaningful for so-called confluent RK methods for which several c l 's take the same value. If the sequence (c l ) l∈{1:s} is nondecreasing, then l (l) = l -1 for all l ∈ {2:s + 1}. The reason for looking for the smallest difference c l -c l (l) is to minimize the CFL restriction on the time step (see Assumption 2.1(i)). For further reference, we define

∆c max := max 2≤l≤s+1 c l -c l (l) . (47)
Notice that ∆c max ≥ 1 s and ∆c max = 1 s whenever all the stages of the ERK method are equidistributed, i.e., c l = l-1 s , l ∈ {1:s + 1}. In the rest of this paper, we simply write l instead of l (l) to simplify the notation.

We can now approximate in time the high-order ODE system (21) by using the IMEX method defined by the two Butcher tableaux in (44). We first set U n,1 := U n . Then, for all l ∈ {2:s + 1}, the l-th stage of the IMEX scheme consists of computing the following high-order update:

M H U n,l := M H U n,l + τ k∈{1:l-1} (a e l,k -a e l ,k ) F H (U n,k ) + G H (U n,k ) -G H,lin (U n ; U n,k ) + τ a i l,l G H,lin (U n ; U n,l ) + τ k∈{1:l-1} (a i l,k -a i l ,k )G H,lin (U n ; U n,k ).
This incremental form of the IMEX scheme is obtained by subtracting the equation defining the IMEX update at stage l from the equation defining the update at stage l. We decompose the above problem into one hyperbolic prediction followed by one parabolic update as follows:

M H W n,l := M H U n,l + τ k∈{1:l-1} (a e l,k -a e l ,k )F H (U n,k ), (48)
M H U n,l -τ a i l,l G H,lin (U n ; U n,l ) := M H W n,l + τ k∈{1:l-1} (a i l,k -a i l ,k )G H,lin (U n ; U n,k ) + τ k∈{1:l-1} (a e l,k -a e l ,k ) G H (U n,k ) -G H,lin (U n ; U n,k ) . ( 49
)
Notice that, following [START_REF] Guermond | Invariant domains and first-order continuous finite element approximation for hyperbolic systems[END_REF], the quasi-linearization process uses the initial state U n at all stages. We now explain how to make the method (48)-(49) IDP.

IDP-IMEX scheme

We proceed as in §2.4 to make the scheme (48)-(49) invariant-domain preserving. Given U n ∈ A I , we set U n,1 := U n and we decompose each stage l ∈ {2:s + 1} into the following four steps:

U n,l (1) -→ (W L,l , W H,l ) (2)
-→

hyperbolic step (48) W n,l (3) -→ (U L,l , U H,l ) (4)
-→

parabolic step (49) U n,l . (50)
Hyperbolic steps ( 1) and ( 2)

The IDP realization of the hyperbolic update (48) is done as in [START_REF] Ern | Invariant-domain-preserving high-order time stepping: I. explicit runge-kutta schemes[END_REF]. One first computes the low-order and high-order hyperbolic updates defined by

M L W L,l := M L U n,l + τ (c l -c l )F L (U n,l ), (51)
M H W H,l := M H U n,l + τ k∈{1:l-1} (a e l,k -a e l ,k )F H (U n,k ). ( 52
)
By proceeding as in ( 35)-(36), we have

W H,l i = W L,l i + τ m i j∈I(i) A n,l ij , ∀i ∈ V, (53)
with A n,l ij :=

k∈{1:l-1} (a e l,k -a e l ,k )F H ij (U n,k ) -(c l -c l )F L ij (U n,l ) (54)
-τ -1 (m ij -m i δ ij )(W H,l j -U n,l j -W H,l i + U n,l i ).
Notice that A n,l is skew-symmetric in compliance with Definition 2.4. Using the hyperbolic limiter, we then set W n,l := hyp (W L,l , A n,l ).

Parabolic steps (3) and ( 4)

We now compute the low-order and high-order parabolic updates defined by solving the following two problems:

M L U L,l -τ (c l -c l )G L,lin (W n,l ; U L,l ) := M L W n,l , (56)
M H U H,l -τ a i l,l G H,lin (U n ; U H,l ) := M H W n,l (57)
+ k∈{1:l-1} τ (a i l,k -a i l ,k )G H,lin (U n ; U n,k ) + k∈{1:l-1} τ (a e l,k -a e l ,k ) G H (U n,k ) -G H,lin (U n ; U n,k ) .
The update U n,l is obtained by employing the conservative parabolic limiter and by proceeding as for the Euler IMEX scheme. Subtracting (56) from (57) yields

U H,l i = U L,l i + τ m i j∈I(i) A n,l ij + τ m i B n,l i , ∀i ∈ V, ( 58
) with A n,l ij := k∈{1:l-1} (a e l,k -a e l ,k ) D H ij (U n,k ) -D H,lin ij (U n ; U n,k ) -(c l -c l )D L,lin ij (W n,l ; U L,l ) + k∈{1:l} (a i l,k -a i l ,k )D H,lin ij (U n ; U H,l ) (59) -τ -1 (m ij -m i δ ij ) U H,l j -W n,l j -U H,l i + W n,l i , and B n,l i := k∈{1:l-1} (a e l,k -a e l ,k ) R H i (U n,k ) -R H,lin i (U n ; U n,k ) -(c l -c l )R L,lin i (W n,l ; U L,l ) + k∈{1:l} (a i l,k -a i l ,k )R H,lin i (U n ; U H,l ). (60)
Notice that A n,l is skew-symmetric and B n,l ∈ B I , in compliance with Definition 2.5. Using the parabolic limiter, we finally set

U n,l := par (U L,l , A n,l , B n,l ). ( 61
)
At the end of the loop, the final update is obtained by setting U n+1 := U n,s+1 .

Remark 3.1 (Quasi-linearization). We observe that the high-order update (21) involves a quasilinearization based on U n . It is essential that the quasi-linearization for the high-order update be the same for all the stages of the IMEX scheme to preserve the high-order accuracy in time of the method. But, the quasi-linearization for the low-order update at each stage l ∈ {1:s + 1} is based on W n,l ; this allows us to invoke the invariant-domain property stated in Assumption [START_REF] Shu | Efficient implementation of essentially non-oscillatory shockcapturing schemes[END_REF].

Remark 3.2 (Complexity). The low-order update (56) requires solving a quasi-linear system for all l ∈ {2:s + 1}. The high-order update (57) requires solving a quasi-linear system for all l ∈ {2:s} and amounts to an explicit update for l = s + 1 because a i s+1,s+1 = 0. Thus, the above method requires solving (2s -1) quasi-linear systems over each time interval.

Conclusion

The main result motivating the construction introduced above is the following assertion. (62)

Assume that the limiters hyp and par match Definitions 2.4-2.5. Let U n ∈ A I . Consider the s-stage IMEX scheme composed of the steps (51)-(61) for l ∈ {2:s + 1}. This scheme satisfies the following properties:

(i) It is well defined;

(ii) It is IDP, i.e., it satisfies (7);

(iii) It is conservative, i.e., it satisfies (8).

Proof. Assume (62) and U n ∈ A I . We are going to show by induction that all the intermediate updates (U n,l ) l∈{1:s+1} are well defined and are in A I . The definition U n,1 := U n implies that the assumption holds true for l = 1. Let now l ∈ {2:s + 1}. We make the following observations: The low-order hyperbolic update (51) has the same structure as (33); the high-order hyperbolic update (53)-(54) has the same structure as (35)-(36); the low-order parabolic update (56) has the same structure as (38); the high-order parabolic update (58)-(61) has the same structure as (40)-( 43).

Hence we can apply Lemma 2.7 to the scheme (51)-(61) provided the effective time step (c l -c l )τ used in the low-order hyperbolic and parabolic stages (51) and ( 56) is in the interval (0, τ * ]. But this is the case owing to the assumption (62). Then Lemma 2.7 implies that U n,l is well defined and is in A I . This establishes (i) and (ii). Notice that Lemma 2.7 also asserts that

i∈V m i U n,l p,i = i∈V m i U n,l p,i , ∀p ∈ C.
An induction argument readily gives i∈V m i U n,l p,i = i∈V m i U n p,i for all l ∈ {1:s + 1}, and the conservation property (8) follows from U n+1 := U n,s+1 . This proves (iii).

Following Shu and Osher [START_REF] Shu | Efficient implementation of essentially non-oscillatory shockcapturing schemes[END_REF] (see also [10, Def. 2.2]), the quantity

c eff := 1 s∆c max , (63)
is called efficiency ratio of the s-stage IMEX scheme. Recall that ∆c max is defined in (47). By construction, we have c eff ≤ 1. Theorem 3.3 shows that the IMEX scheme is IDP for all τ ∈ (0, c eff sτ * ]. Hence it desirable to have an efficiency ratio as large as possible for computational efficiency. In particular, the largest time allowed is s×τ * when c eff = 1. The optimal value c eff = 1 is attained when the coefficients c j are equi-distributed, i.e., c j := j-1 s for all j ∈ {1:s}.

Examples of IMEX schemes

In this section, we review some examples of IMEX schemes and introduce some novel schemes. We only consider schemes with p-order accuracy where p ∈ {2, 3, 4}. Recall that the IMEX schemes under consideration combine an ERK scheme and an EDIRK scheme. Both schemes consist of s ≥ p stages and are described by the Butcher tableaux introduced in (44). In what follows, we use the terminology IMEX(s, p; c eff ) for an IMEX scheme with s stages, order p, and efficiency ratio c eff . Four new schemes with optimal efficiency and the following characteristics are introduced in this section:

(1) IMEX(3, 3; 1), singly diagonal, A-stable implicit part, (81);

(2) IMEX(4, 3; 1), singly diagonal, L-stable implicit part, (82);

(3) IMEX(5, 4; 1), singly diagonal, L-stable implicit part, (83); (4) IMEX(6, 4; 1), singly diagonal, L-stable implicit part, (84).

Main properties of IMEX schemes

Three important notions for IMEX schemes are the consistency order, the stability of the implicit scheme, and the efficiency ratio. We now briefly discuss these three properties.

For simplicity, we focus on IMEX schemes for which we have b e i = b i i =: b i for all i ∈ {1:s}. We denote by B the row vector in R s having components (b i ) i∈{1:s} . We denote by C the column vector in R s having components (c j ) j∈{1:s} . We also use the notation C p , p ≥ 0, for the column vector in R s having components (c p j ) j∈{1:s} . To be coherent with the literature we set U := C 0 = (1, . . . , 1) T and use the symbol C instead of C 1 . We denote by B C the row vector in R s having components (b j c j ) j∈{1:s} . We denote by A e (resp., A i ) the square matrix of order s with entries (a e i,j ) i,j∈{1:s} (resp., (a i i,j ) i,j∈{1:s} ). Notice that A e is strictly lower triangular, whereas A i is lower triangular.

We adopt the following terminology often used in the literature. The identity matrix of order s is denoted I s .

Consistency order

Recall that necessary consistency conditions for the explicit and the implicit methods to be separately of order p are BA r-1 C q-1 = (q -1)! (q -1 + r)! , ∀r ∈ {1:p}, q ∈ {1:p -r + 1}.

These conditions are sufficient for p ≤ 2. They are also sufficient for all p ≥ 2 if the ODE systems are autonomous and linear. Additional nonlinear conditions must be enforced for nonlinear autonomous systems when p ≥ 2. Coupling conditions must be added for IMEX schemes to be of order p ≥ 2.

The consistency properties of IMEX methods are reviewed in Pareschi and Russo [25, §2.1] and Kennedy and Carpenter [18, §2.2]. The analysis therein is based on the following simplifying assumption (see ( 45)), which we systematically enforce:

A e U = C, A i U = C. (65)
The (linear order) conditions to achieve second-order are

BU = 1, BC = 1 2 , (66)
while the conditions BA e U = BA i U = 1 2 follow from (65) and (66). The conditions to achieve third-order accuracy are (65)-(66) together with the following (linear order) conditions

BC 2 = 1 3 , BA e C = BA i C = 1 6 , (67)
while the conditions B(A e ) 2 U = B(A i ) 2 U = 1 6 follow from (65) and (67). The conditions to achieve fourth-order accuracy are (65), (66), (67), together with the (linear order) conditions

BC 3 = 1 4 , BA e C 2 = BA i C 2 = 1 12 , B(A e ) 2 C = B(A i ) 2 C = 1 24 , (68)
the (nonlinear order) condition

(B C)A e C = (B C)A i C = 1 8 , (69)
and the coupling condition

BA e A i C = BA i A e C = 1 24 . (70)
The conditions B(A e ) 3 U = B(A i ) 3 U = 1 24 follow from (65) and (68). Finally, we are also going to make use of the fifth-order linear order conditions for a six-stage, fourth-order method

BC 4 = 1 5 , BA e C 3 = BA i C 3 = 1 20 , B(A e ) 2 C 2 = B(A i ) 2 C 2 = 1 60 , B(A e ) 3 C = B(A i ) 3 C = 1 120 . (71)
The conditions B(A e ) 4 U = B(A i ) 4 U = 1 120 follow from (65) and (71).

Stability

The amplification function associated with a DIRK scheme is ). For DIRK schemes, A i is invertible if all the diagonal entries of A i are nonzero. In this case, L-stability amounts to B(A i ) -1 U = 1. However, for EDIRK schemes, the first diagonal entry of A i is zero. In this case, one considers the block decompositions

R(z) := 1 + zB(I s -zA i ) -1 U, z ∈ C. (72)
A i = 0 0 α Ã , B = (β, B), (73)
with α ∈ R s-1 (column vector), Ã ∈ R s-1,s-1 , β ∈ R, and B ∈ R s-1 (row vector). Then, the amplification function defined in (72) can be rewritten as

R(z) = 1 + zβ + z B(I s-1 -z Ã) -1 ( Ũ + zα), (74)
where Ũ ∈ R s-1 is the column vector having all entries equal to one. Assuming that à is invertible, one readily verifies that the EDIRK scheme is L-stable if it is A-stable and if the following holds:

β = B Ã-1 α, B Ã-1 Ũ + B Ã-2 α = 1. (75)
Notice that the first condition in (75) implies that lim t→-∞ R(t) = 1 -B Ã-1 Ũ -B Ã-2 α, and the second condition then implies that lim t→-∞ R(t) = 0.

Second-order IMEX schemes

A first possibility to obtain a two-stage, second-order IMEX method consists of combining Heun's second-order scheme with the Crank-Nicolson (A-stable) scheme. The Butcher tableaux are

0 0 1 1 0 1 1 2 1 2 0 0 1 1 2 1 2 1 1 2 1 2 (76)
We have l (l) = l -1 for all l ∈ {2:3}, and the efficiency ratio is c eff = 1 2 . We call this method IMEX(2, 2; 1 2 ). A second possibility consists of combining the explicit and implicit (A-stable) midpoint rules. The corresponding Butcher tableaux are 0 0

1 2 1 2 0 1 0 1 0 0 1 2 0 1 2 1 0 1 (77)
We have l (l) = l -1 for all l ∈ {2:3}, and in this case the efficiency ratio reaches the optimal value c eff = 1. We call this method IMEX(2, 2; 1). The amplification function is R(z) = 2+z 2-z for the Crank-Nicolson scheme and the midpoint rule. It is remarkable that the amplification function is the same for both schemes. However, the efficiency of the Crank-Nicolson scheme is only 1 2 , whereas that of the midpoint rule is 1.

A third possibility (see Ascher et al. [2,Sec. 2.5]) is to consider a three-stage, second-order scheme in which the implicit scheme is an L-stable, zero-padded, two-stage ESDIRK scheme. The Butcher tableaux are 0 0

γ γ 0 1 δ 1 -δ 0 1 0 1 -γ γ 0 0 γ 0 γ 1 0 1 -γ γ 1 0 1 -γ γ (78)
with γ := 1 -1 √ 2 ≈ 0.29289 and δ is an adjustable parameter for which the value δ = -2 3 √ 2 is recommended. We have l (l) = l -1 for all l ∈ {2:4}, but the efficiency ratio is only c eff = 1 3 (1 -γ) ≈ 0.24. We call this method IMEX(3, 2; 0.24). Remark 4.1 (Strang's splitting). Strang's splitting can be rewritten as an IMEX scheme. Consider for instance that the explicit (resp., implicit) midpoint rule is used for the explicit (resp., implicit) steps. One can verify that the whole process can be rewritten as a five-stage IMEX scheme with the following Butcher tableaux 0 0

1 4 1 4 0 1 2 0 1 2 0 1 2 0 1 2 0 0 3 4 0 1 2 0 1 4 0 1 0 1 2 0 0 1 2 0 0 1 4 0 0 1 2 0 0 1 2 1 2 0 0 1 0 3 4 0 0 1 0 0 1 0 0 1 0 0
As expected, there is only one implicit substep (the third one). We have l (l) = (1, 2, 3, 4, 5) for l ∈ {2:6} with max l∈{2:6} (c l -c l ) = 1 4 . Notice that the fourth substep does not involve extra flux computations with respect to the third substep. Hence, the efficiency ratio is c eff = 4s -1 with s = 4 (rather than c eff = 4s -1 with s = 5), i.e., the method has optimal efficiency. A variant of this method is implemented in [START_REF] Guermond | Second-order invariant domain preserving approximation of the compressible Navier-Stokes equations[END_REF] to solve the compressible Navier-Stokes equations (the method SPPRK(3, 3) is used therein instead of the midpoint rule though).

Third-order IMEX schemes

In this section, we consider third-order IMEX schemes composed of three or four stages. Threestage schemes in which the implicit scheme is A-stable are available in the literature, but none of these methods has optimal efficiency. We derive here a three-stage IMEX scheme achieving optimality and whose implicit tableau is A-stable. We also construct a four-stage scheme with optimal efficiency whose implicit tableau is L-stable.

Three-stage schemes

A first possibility to obtain a three-stage, third-order IMEX method consists of using the two-stage, third-order, zero-padded ESDIRK scheme (Crouzeix [8], Nørsett [START_REF] Nørsett | Semi explicit Runge-Kutta methods[END_REF]) for the implicit scheme and combining it with the three-stage, third-order ERK scheme sharing the same coefficients c j and b j . The corresponding Butcher tableaux are (see Ascher et al. [2,Sec. 2.4])

0 0 γ γ 0 1 -γ γ -1 2 -2γ 0 1 0 1 2 1 2 0 0 γ 0 γ 1 -γ 0 1 -2γ γ 1 0 1 2 1 2 (79)
This ESDIRK scheme is L-stable.

The companion ERK method that shares the same set of coefficients c j and b j has already been proposed in Ern and Guermond [START_REF] Ern | Invariant-domain-preserving high-order time stepping: I. explicit runge-kutta schemes[END_REF]. The six coefficients of the matrix A e are obtained by enforcing the fourth-order linear consistency conditions (three linear equations from (65), one linear equation form (67), one linear and one nonlinear equation from ( 68)) This is the only four-stage, third-order ERK method with optimally distributed coefficients that is also fourth-order accurate on linear problems. Its Butcher tableau is as follows:

0 0 1 4 1 4 0 1 2 0 1 2 0 3 4 0 1 4 1 2 0 1 0 2 3 -1 3 2 3 (82b)
Near the origin along the imaginary axis, we have |R e (i )| = 1+ρ e 6 6 +O( 8) with ρ e 6 = -2B(A e ) 4 C+ 2B(A e ) 3 C -B(A e ) 2 C + 1 36 = -1 72 . Remark 4.2 (Other four-and five-stage methods). A third-order IMEX method combining a four-stage ERK method with a four-stage L-stable DIRK method is described in Ascher et al. [2, Sec. 2.7] (the DIRK method has actually three stages since the third and fourth stages are identical). The efficiency ratio is close to c eff = 0.46. A variant of the implicit four-stage scheme is studied in Calvo et al. [START_REF] Calvo | Linearly implicit Runge-Kutta methods for advectionreaction-diffusion equations[END_REF]. A method combining a five-stage ERK method with a five-stage L-stable DIRK method is described in [2, Sec. 2.8] (the DIRK method has actually four stages since the fourth and fifth stages are identical). The efficiency ratio is only c eff = 1 8 .

Fourth-order IMEX schemes

Some fourth-order IMEX methods composed of five and six stages are discussed in Kennedy and Carpenter [18, Sec. 3.2 & 3.3], but the efficiency ratio of these methods is far from being optimal.

Here, we devise five-and six-stage, fourth-order IMEX schemes with optimal efficiency; the implicit scheme is L-stable in both cases. We call these methods IMEX(5, 4; 1) and IMEX(6, 4; 1).

Five-stage scheme

We set c l := l-1 5 for all l ∈ {1:5} to achieve optimal efficiency. There are 19 coefficients to be determined: 14 entries in the matrix A i and the five components of the vector b. We enforce Butcher's simplifying assumption (65) (4 equations), the (linear) order conditions (66), (67), and (68) (7 equations), the (nonlinear) order condition (69) (one equation), and the two conditions in (75) which are necessary for L-stability. This gives 14 equations. We additionally require that the scheme be singly diagonal, giving the three additional equations a i 2,2 = a i 3,3 = a i 4,4 = a i 5,5 , and that a i 5,1 = 0, giving one additional equation. The resulting under-determined set of nonlinear equations (18 equations, 19 unknowns) is solved using julia with 2×10 -16 tolerance. The following solution is found:

This ESDIRK scheme is L-stable. Along the imaginary axis near the origin, we have

|R i (i )| = 1 + ρ i 6 6 + O( 8 ) with ρ i 6 ≈ -0.0846, where ρ i 6 = -2B(A i ) 4 C + 2B(A i ) 3 C -1 72
. We now devise the companion ERK scheme that shares the same set of coefficients c j and b j . There are 10 unknowns (the entries of the strictly lower triangular matrix A e ). We enforce Butcher's simplifying assumption (65) (4 equations) and the (linear) order conditions (67) involving the matrix A e (three equations, since the remaining order conditions have already been accounted for in the design of the ESDIRK scheme above), the (nonlinear) order condition (69) (one equation), and the two coupling conditions (70). This gives 10 equations. The resulting set of nonlinear equations (10 equations, 10 unknowns) is solved using julia with 2×10 -16 tolerance. The following solution is found:

0 0 1 5 0.2 0 2 5
0.26075582269554909 0.13924417730445096 0

We have |R e (i )| = 1 + ρ e 6 6 + O( 8) with ρ e 6 ≈ -0.0148 along the imaginary axis near the origin, where ρ e 6 = 2B(A e ) 3 C -1 72 (notice that 2B(A e ) 4 C = 0).

Six-stage schemes

A six-stage, fourth-order method with L-stable (actually, stiffly accurate) implicit scheme is designed in Calvo et al. [START_REF] Calvo | Linearly implicit Runge-Kutta methods for advectionreaction-diffusion equations[END_REF] using a L-stable, six-stage (actually five distinct stages) fourth-order SDIRK method from Hairer and Wanner [START_REF] Hairer | Solving Ordinary Differential Equations. II. Stiff and Differentialalgebraic Problems[END_REF]. More precisely, the explicit tableau is given by Equation ( 14) in [START_REF] Calvo | Linearly implicit Runge-Kutta methods for advectionreaction-diffusion equations[END_REF] and the implicit tableau is given by Equation (IV.6.16) in [START_REF] Hairer | Solving Ordinary Differential Equations. II. Stiff and Differentialalgebraic Problems[END_REF] (see also Table IV.6.5). The efficiency ratio of this method is only c eff = 1 12 ≈ 0.08. We call this method IMEX(6, 4; 0.08).

We now propose a six-stage, fourth-order IMEX method with optimal efficiency. We set c l := l-1 6 and l (l) = l -1 for all l ∈ {1:6} to achieve optimal efficiency. There are 26 coefficients to be determined for the EDIRK scheme: 20 entries in the matrix A i and the six components of the vector b. We enforce Butcher's simplifying assumption (65) (5 equations), the (linear) order conditions (66), (67), (68) (7 equations), the (nonlinear) order condition (69) (one equation), and the two conditions in (75) which are necessary for L-stability. We also enforce the fifth-order linear order conditions (71) (4 equations). This gives 19 equations. We additionally require that the scheme be singly diagonal, giving the four equations a i 2,2 = a i 3,3 = a i 4,4 = a i 5,5 = a i 6,6 . We also set b 4 = 0.47, giving six additional equations. The resulting set of nonlinear equations (24 equations, 26 unknowns) is solved using julia with 4×10 -16 tolerance. The following solution is found:

0 0 1 6 -0.1113871744697862 0.2780538411364528 2
We have |R i (i )| = 1+ρ i 6 6 +O( 8) along the imaginary axis near the origin, with ρ i 6 ≈ -1.06×10 -3 , where we recall that ρ i

6 := -2B(A i ) 4 C + 2B(A i ) 3 C -1 72 = -2B(A i ) 4 C + 1 360
. We now proceed to find a companion ERK scheme sharing the same set of coefficients c j and b j . There are 15 unknowns (the entries of the strictly lower triangular matrix A e ). We enforce Butcher's simplifying assumption (65) (5 equations), the (linear) order conditions (67), (68) involving the matrix A e (three equations), the (nonlinear) order condition (69) (one equation), and the two coupling conditions (70). This gives 11 equations. We also enforce the fifth-order linear order conditions (71) (4 equations). In total we have 14 equations. The resulting under-determined set of nonlinear equations (14 equations, 15 unknowns) is solved using julia. The following solution is found with 4×10 -16 tolerance:

0 0 0 0 0 1 6
0.1666666666666667 0 0 0

Numerical illustrations

We illustrate the IMEX methods proposed in the paper. We start with convergence tests on an ODE system. Then we solve a scalar nonlinear conservation equation with hyperbolic and parabolic fluxes.

Convergence tests

We test the convergence properties of the new IMEX methods proposed in this paper and compare them to the other published methods listed in Section 4. Following Kennedy and Carpenter [18, §5.1], we consider the 2×2 ODE system

∂ t y 1 (t) = -2y 1 + -1 (y 2 2 -y 1 ), ∂ t y 2 (t) = y 1 -y 2 -y 2 2 , (85)
with > 0 and initial condition y 1 (0) = y 2 (0) = 1. The solution is y 1 (t) = y 2 2 (t), y 1 (t) = e -t . As → 0, the above problem degenerates into the index-1 differential algebraic equation ∂ t y 2 (t) = y 1 -y 2 -y 2 2 with y 1 = y 2 2 . We denote by U := (u 1 , u 2 ) T the approximate solution produced by the IMEX methods. Referring to [START_REF] Guermond | Invariant domains and first-order continuous finite element approximation for hyperbolic systems[END_REF] for the notation, we set M := I 2 , where I 2 is the 2×2 identity matrix, and

F(U) := (-2u 1 , u 1 -u 2 -u 2 2 ) T , G(U) := ( -1 (u 2 2 -u 1 ), 0) T , G lin (W, U) := G(U). (86)
Notice that G lin is not linear, but solving the problem (10) is simple:

(I -τ G lin (W, •)) -1 (W) = ( 1 +τ ( w 1 + τ w 2 2 ), w 2 ) T . ( 87
)
We test all the methods mentioned in Section 4 by solving the above problem over the time interval [0, T ] with T = 4, the initial data (u 1 (0), u 2 (0)) = (1, 1), and for ∈ {1, 10 -6 }. For each method, we compute the errors |u 1 (T )-y 1 (T )|/|y 1 (T )+y 2 (T )| and |u 2 (T )-y 2 (T )|/|y 1 (T )+y 2 (T )|. The results are reported in Figures 1,2,3. The symbols "y1,e0", "y1,e-6", "y2,e0", and "y2,e-6" in the legend refer to the error on the variable y 1 with = 10 0 , the error on the variable y 1 with = 10 -6 , the error on the variable y 2 with = 1, and the error on the variable y 2 with = 10 -6 , respectively. 2 ), IMEX(2, 2; 1), IMEX(3, 2; 0.24) (with the Butcher tableaux (76), (77), (78) from left to right).

The three panels in Figure 1 show the results for the second-order methods IMEX(2, 2; 1 2 ), IMEX(2, 2; 1) and IMEX(3, 2; 0.24) (from left to right and with the Butcher tableaux (76), (77), (78)). We observe that they all deliver second-order accuracy uniformly with respect to for both y 1 and y 2 .

The results for the methods IMEX(3, 3; 0.26), IMEX(3, 3; 1), IMEX(4, 3; 1) (see (79), (81), and (82) respectively) are shown in Figure 2. We observe third-order accuracy for the three methods on both variables when = 1, but as expected, the convergence on y 1 reduces to second-order when = 10 -6 . This order reduction in the pre-asymptotic range (i.e., < τ ) is well documented in the literature and we refer the reader to Boscarino and Pareschi [START_REF] Boscarino | On the asymptotic properties of IMEX Runge-Kutta schemes for hyperbolic balance laws[END_REF] for an analysis of this phenomenon. We observe that the two new methods introduced in this paper (i.e., IMEX(3, 3; 1), IMEX(4, 3; 1)) perform as expected. Recall that IMEX(3, 3; 0.26) and IMEX(3, 3; 1) are both three-stage, thirdoder methods. The efficiency of IMEX(3, 3; 0.26) is only 0.26 whereas that of IMEX(3, 3; 1) is optimal. The implicit components of these methods are A-stable. They also have the same amplification functions. The method IMEX(4, 3; 1) is composed of four stages, has efficiency 1, and its implicit component is L-stable.

Finally, the results for the three fourth-order methods IMEX(5; 4; 1) (see (83)), IMEX(6; 4; 0.08), and IMEX(6, 4; 1) (see (84)) are shown in Figure 3. We observe fourth-order accuracy for the three methods on both variables when = 1. The convergence rate on y 2 is still 4 when = 10 -6 , but the convergence rate for y 1 reduces to second order when = 10 -6 . We observe that the two new methods introduced in this paper (i.e., IMEX(5; 4; 1) and IMEX(6; 4; 1)) perform as expected. Recall that IMEX(5; 4; 1) has five stages with efficiency c eff = 1, IMEX(6; 4; 0.08) has six stages with efficiency c eff = 0.08, and IMEX(6; 4; 1) has six stages with efficiency c eff = 1. The implicit components of these three methods are L-stable.

Nonlinear scalar conservation equation

In this section, we illustrate the method on the following scalar nonlinear conservation equation:

∂ t u + ∇•f (u) -∆u = 0, x ∈ D ∞ , t > 0, (88)
posed in the two-dimensional domain D ∞ := R×(0, 1). The flux is defined by f (u) := (u(1-u), 0) T . With the notation x := (x, y), the initial data is

u 0 (x) := µ + δ tanh δ (x -x 0 ) , µ := 1 2 (u L + u R ), δ := 1 2 (u R -u L ). (89)
Assuming homogeneous Neuman boundary conditions on the top and bottom parts of the domain, the solution to this Cauchy problem is a wave moving at speed s := 1 -2µ: v(x, t) = u 0 (x -st) with s := (s, 0).

The method described in this paper is implemented using continuous finite elements. The tests are done with continuous P 1 and P 3 finite elements. The low-order solution method for the hyperbolic subproblem is fully described in [START_REF] Guermond | Invariant domains and first-order continuous finite element approximation for hyperbolic systems[END_REF]. The high-order method and the limiting are described in [START_REF] Guermond | Second-order invariant domain preserving approximation of the Euler equations using convex limiting[END_REF][START_REF] Guermond | Invariant domain preserving discretizationindependent schemes and convex limiting for hyperbolic systems[END_REF]. We use FCT to perform the limiting as the problem is scalar-valued. Local bounds are used at every grid point. Relaxation of the bounds guaranteeing high-order convergence is done as explained in [12, §4.7.1] and [13, §7.6].

We set u L := -1 and u R := 1, so that the solution to (88) is a wave moving at speed s = 1. The numerical simulations are done in the truncated computational domain D := (x L , x R )×(y B , y T ) with x L = y B := 0, x R := 1, y T := Neumann conditions are enforced on ∂D par N for the parabolic subproblems. The enforcement of the boundary condition for the hyperbolic subproblems is done at the end of each stage of the IMEX step.

In all the tests, the time step is computed by using the expression

τ := CFL×s×τ * , (91)
where CFL > 0 is a fixed parameter, s is the number of stages of the IMEX method, and τ * is the maximum time step for which the low-order hyperbolic update is IDP; see Assumption 2.1(i). This definition guarantees that for a given simulation time T , the total number of flux evaluations (which is a measure of the algorithmic cost) is approximately s× T τ = s× T CFL×s×τ * = T CFL×τ * . Hence, for a given mesh, a given final time T , and a given CFL number, the algorithmic cost of two IMEX methods with different number of stages is approximately identical. The simulations are done up to T = 1 2 for = 2×10 -n , n ∈ {2, 3, 4}. We use unstructured Delaunay meshes. We test the following five methods: IMEX(2, 2; 1); IMEX(3, 3; 1); IMEX(4, 3; 1); IMEX(5, 4; 1); and IMEX(6, 4; 1). All the errors are evaluated at T and are relative.

We show the errors and the convergence rates for continuous P 1 elements in Table 5.2. We observe that the methods deliver second-order accuracy when the mesh size is small enough to capture the viscous layer of size . The accuracy is limited to second-order due to our using P 1 elements. We also notice that all the methods deliver first-order accuracy when the mesh size cannot capture the viscous layer. First-order accuracy is optimal in this case.

We show the errors and the convergence rates for continuous P 3 elements in Table 5.2. The methods deliver optimal accuracy when the mesh size is small enough to capture the viscous layer, that is, second-order for IMEX(2, 2; 1) and fourth-order for the other methods. There seems to be some super-convergence effect for the third-order methods IMEX(3, 3; 1) and IMEX(4; 4; 1). Here again, all the methods deliver first-order accuracy when the mesh size cannot capture the viscous layer.

Table 1: Problem (88) for = 2×10 -n , n ∈ {2, 3, 4}. P 1 finite elements. Error in the L 1 -norm. First row: IMEX(2, 2; 1). Second row: IMEX(3, 3; 1) and IMEX(4, 3; 1). Third row IMEX(5, 4; 1) and IMEX(6, 4; 1).

Conclusions

A new time stepping technique making every IMEX method invariant-domain preserving has been introduced. New IMEX methods with optimal efficiency have been constructed. The numerical experiments done in §5.1 and §5.2 demonstrate that the new IMEX methods proposed in this paper behave as predicted by the theory. All the methods tested are invariant-domain preserving and deliver the expected accuracy. A natural perspective of this work is to show how present methodology can be used to solve nonlinear systems of conservation equations like the compressible Navier-Stokes equations.

theorem 3 . 3

 33 (s-stage IDP-IMEX). Let Assumption 2.1 hold and let τ ∈ 1, τ * ∆c max .

  05144383172900784 1.17898889035791732 -0.90036112111104449 . . . 1 -0.10511678454691901 0.87880047152100838 -0.58903404061484477 0.46213380485434047 . . .

4 5 .

 45 . . 0.57281606248213512 1 . . . 0.35321654878641495

  0.51534223099602405 -0.81662794199265554 0.88505294668159373 . . . 1 -0.10511678454691901 0.87880047152100838 -0.58903404061484477 0.46213380485434047 . . .

6 6

 6 have |R e (i )| = 1+ρ e +O( 8 ) along the imaginary axis near the origin, with ρ e 6 ≈ -9.67×10 -5 , where ρ e 6 := -2B(A e ) 4 C + 2B(A e ) 3 C -1 72 = -2B(A e ) 4 C + 1 360 . The value of the coefficient b 4 is adjusted to the value 0.47 to make ρ e 6 negative.

Figure 1 :

 1 Figure1: Convergence test on problem (85) for the second-order methods IMEX(2, 2; 1 2 ), IMEX(2, 2; 1), IMEX(3, 2; 0.24) (with the Butcher tableaux (76), (77), (78) from left to right).

Figure 2 :

 2 Figure 2: Convergence test on problem (85) for the third-order methods IMEX(3, 3; 0.26), IMEX(3, 3; 1); IMEX(4, 3; 1) (with the Butcher tableaux (79), (81), (82) from left to right).

1 4 .

 4 Let ∂D hyp D = ∂D par D := {x L , x R }×(y B , y T ), := ({x L } ∪ {x R })×(y B , y T ), and ∂D par N := (x L , x R )×{y B , y T } (it happens here that ∂D hyp D = ∂D par D because 1 -2u R < 0). At each stage of the IMEX method, Dirichlet boundary conditions are enforced at ∂D hyp D for the hyperbolic subproblems and at ∂D par D for the parabolic subproblems. Homogeneous

Liu and Osher [21, Thm. 1], and Zhang and Shu

	.1], Coquel and
	LeFloch [7, Thm .4.3],

  , [7, Thm .4.3], [21, Thm. 2], or Zhang and Shu [31, Lem. 2.4], or other nonlinear variants like convex limiting (see, e.g., Guermond et al.

  62E-02 1.20 1.59E-02 1.22 7.65E-03 1.16 7.42E-03 1.18 2.93E-03 1.47 2.81E-03 1.49 7.22E-04 2.06 6.95E-04 2.05 1.33E-04 2.47 1.28E-04 2.46 = 10 -4 (5,4;1) rate (6,4;1) rate 3.60E-02 -3.59E-02 -1.52E-02 1.41 1.50E-02 1.42 7.65E-03 1.07 7.56E-03 1.06 4.01E-03 0.99 3.98E-03 0.99 1.99E-03 1.03 1.96E-03 1.04 9.86E-04 1.02 9.63E-04 1.04

		= 10 -2		= 10 -3		= 10 -4
	I	IMEX(2,1;1) rate	(2,1;1) rate	(2,1;1) rate
	106	1.98E-02	-	3.36E-02 -	3.60E-02 -
	360	4.13E-03 2.56	1.61E-02 1.20	1.52E-02 1.41
	1309	8.12E-04 2.52	7.60E-03 1.16	7.64E-03 1.07
	4825	2.03E-04 2.13	2.88E-03 1.49	4.02E-03 0.98
	18846 4.99E-05 2.06	7.01E-04 2.07	1.99E-03 1.03
	74510 1.25E-05 2.01	1.29E-04 2.46	9.83E-04 1.03
			= 10 -2		= 10 -3	= 10 -4
	I	IMEX(3,3;1) rate IMEX(4,3;1) rate	(3,3;1) rate (4,3;1) rate	(3,3;1) rate (4,3;1) rate
	106	1.97E-02	-	1.97E-02	-	3.36E-02 -3.36E-02 -	3.60E-02 -3.60E-02 -
	360	4.13E-03 2.56 4.13E-03 2.56	1.62E-02 1.20 1.61E-02 1.20	1.52E-02 1.41 1.52E-02 1.41
	1309	8.26E-04 2.49 8.27E-04 2.49	7.62E-03 1.17 7.60E-03 1.17	7.66E-03 1.06 7.63E-03 1.06
	4825	2.04E-04 2.15 2.04E-04 2.15	2.88E-03 1.49 2.88E-03 1.49	4.03E-03 0.98 4.01E-03 0.98
	18846 5.00E-05 2.06 5.00E-05 2.06	7.03E-04 2.07 7.04E-04 2.07	2.00E-03 1.03 1.99E-03 1.03
	74510 1.25E-05 2.02 1.25E-05 2.02	1.32E-04 2.44 1.32E-04 2.44	9.84E-04 1.03 9.82E-04 1.03
			= 10 -2		= 10 -3
	I	IMEX(5,4;1) rate IMEX(6,4;1) rate	(5,4;1) rate (6,4;1) rate
	106	1.98E-02	-	1.97E-02	-	3.36E-02 -3.35E-02 -
	360	4.13E-03 2.56 4.10E-03 2.57	1.
	1309	8.26E-04 2.49 8.11E-04 2.51
	4825	2.03E-04 2.15 2.03E-04 2.13
	18846 4.99E-05 2.06 4.99E-05 2.06
	74510 1.25E-05 2.02 1.25E-05 2.02

Table 2 :

 2 Problem (88) for = 2×10 -n , n ∈ {2, 3, 4}. P 3 finite elements. Error in the L 1 -norm. First row: IMEX(2, 2; 1). Second row: IMEX(3, 3; 1) and IMEX(4, 3; 1). Third row IMEX(5, 4; 1) and IMEX(6, 4; 1). 73E-09 4.09 3.73E-09 4.09 665302 2.65E-10 3.83 2.65E-10 3.83 = 10 -3 (5,4;1) rate (6,4;1) rate 1.80E-02 -1.80E-02 -6.79E-03 1.48 6.79E-03 1.48 2.07E-03 1.77 2.06E-03 1.77 2.95E-04 2.92 2.94E-04 2.93 6.66E-06 5.50 6.65E-06 5.50 2.29E-07 4.88 2.29E-07 4.87 = 10 -4 (5,4;1) rate (6,4;1) rate 1.93E-02 -1.93E-02 -8.16E-03 1.31 8.16E-03 1.31 4.04E-03 1.05 4.04E-03 1.05 2.00E-03 1.06 2.00E-03 1.06 9.09E-04 1.14 9.09E-04 1.14 3.20E-04 1.51 3.20E-04 1.51

		= 10 -2		= 10 -3			= 10 -4
	I	IMEX (2,1;1) rate	(2,1;1) rate	(2,1;1) rate
	778	1.48E-03	-	1.80E-02 -		1.94E-02 -
	2902	2.10E-05	6.46	6.79E-03 1.48	8.16E-03 1.31
	11113	1.18E-06	4.29	2.07E-03 1.77	4.04E-03 1.05
	42097	1.45E-07	3.14	2.95E-04 2.92	2.00E-03 1.06
	166966	3.02E-08	2.28	6.82E-06 5.47	9.09E-04 1.14
	665302	7.44E-09	2.03	2.62E-07 4.72	3.20E-04 1.51
			= 10 -2			= 10 -3	= 10 -4
	I	IMEX (3,3;1) rate IMEX (4,3;1) rate	(3,3;1) rate (4,3;1) rate	(3,3;1) rate (4,3;1) rate
	778	1.48E-03	-	1.48E-03		-	1.80E-02 -1.80E-02 -	1.93E-02 -1.93E-02 -
	2902	2.01E-05	6.53	2.00E-05	6.53	6.79E-03 1.48 6.79E-03 1.48	8.16E-03 1.31 8.16E-03 1.31
	11113	9.46E-07	4.55	9.45E-07	4.55	2.06E-03 1.77 2.07E-03 1.77	4.04E-03 1.05 4.04E-03 1.05
	42097	6.28E-08	4.07	6.25E-08	4.08	2.94E-04 2.92 2.95E-04 2.92	2.00E-03 1.06 2.00E-03 1.06
	166966	3.79E-09	4.08	3.73E-09	4.09	6.67E-06 5.50 6.67E-06 5.50	9.09E-04 1.14 9.09E-04 1.14
	665302	2.79E-10	3.77	2.67E-10	3.82	2.29E-07 4.88 2.29E-07 4.88	3.20E-04 1.51 3.20E-04 1.51
			= 10 -2
	I	IMEX(5,4;1) rate IMEX(6,4;1) rate
	778	1.48E-03	-	1.47E-03	-
	2902	2.00E-05 6.53 2.00E-05 6.53
	11113	9.45E-07 4.55 9.46E-07 4.55
	42097	6.26E-08 4.08 6.27E-08 4.08
	166966 3.

Acknowledgments

This material is based upon work supported in part by the National Science Foundation via grants DMS2110868; the Air Force Office of Scientific Research, USAF, under grant/contract number FA9550-18-1-0397; by the Army Research Office under grant/contract number W911NF-19-1-0431; and the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contracts B640889. The support of INRIA through the International Chair program is also acknowledged.

* This material is based upon work supported in part by the National Science Foundation grant DMS2110868, the Air Force Office of Scientific Research, USAF, under grant/contract number FA9550-18-1-0397, the Army Research Office, under grant number W911NF-19-1-0431, and the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contracts B640889. The support of INRIA through the International Chair program is acknowledged.

with γ := 1 2 + 1 2 √ 3 ≈ 0.78867 (i.e., γ 2 = γ - 1 6 ). The amplification function is

The zero-padded ESDIRK scheme is A-stable, but not L-stable because we only have lim t→-∞ R(t) = 1 -√ 3 ≈ -0.73205. Finally, we observe that the values for l are (1, 1, 2), and the efficiency ratio is only c eff = 1 3 γ ≈ 0.26. We call this method IMEX (3, 3; 0.26). We now propose a three-stage, third-order IMEX method with optimal efficiency. We call this method IMEX(3, 3; 1). We use the third-order Heun method for the ERK part, and we design the corresponding three-stage, third-order EDIRK scheme. To this purpose, we first request that the EDIRK scheme has the same set of coefficients c j and b j as Heun's method, so that there remains to determine the matrix A i . Since this matrix is lower triangular, of order three, and a i 1,1 = 0, this leaves five entries to be determined. Four equations can be enforced: two from Butcher's simplifying assumption (65) (there are two equations corresponding to the rows i ∈ {2, 3} in (65) since the row corresponding to i = 1 is trivial), one is the (linear order) condition BA i C = 1 6 stated in (67) (the remaining linear order conditions are already satisfied), and one is the first stability condition in (75). One can show that it is not possible to enforce the second equality in (75) (there would be no solution). The fifth condition we use to close the system consists of minimizing lim t→-∞ R(t). Solving this problem leads to an A-stable method with lim t→-∞ R(t) = 1 -√ 3. Incidentally, the implicit scheme turns out to be singly diagonal, although this property has not been enforced explicitly. The Butcher arrays of the ERK and ESDIRK scheme are as follows:

(81) with (again) γ := 1 2 + 1 2 √ 3 ≈ 0.78867. We have l (l) = l -1 for all l ∈ {2:4}, and the efficiency ratio reaches the optimal value c eff = 1. Quite remarkably, the amplification function for the above ESDIRK scheme is still given by (80). Hence, the amplification functions of the methods described by the Butcher tableaux (79) and (81) are identical, but the efficiency of (79) is only 1 3 γ ≈ 0.26, whereas that of the new method (81) is 1.

Four-stage schemes

It is possible to devise a four-stage, third-order IMEX method with optimal efficiency in which the implicit part is an ESDIRK L-stable scheme. We call this method IMEX(4, 3; 1). We set c l := l-1 4 for all l ∈ {1:4} to achieve optimal efficiency. There are 13 coefficients to be determined: 9 entries in the matrix A i and the four components of the vector b. We enforce Butcher's simplifying assumption (65) (3 equations), the (linear) order conditions (66) and (67) (4 equations), and the two conditions in (75) which are necessary for L-stability. This gives 9 equations. We additionally require that the scheme be singly diagonal, giving the two additional equations a i 2,2 = a i 3,3 = a i 4,4 , and that BC 3 = 1 4 (this is the first of the fourth-order (linear) conditions in (68)). The resulting under-determined set of nonlinear equations (12 equations, 13 unknowns) is solved using julia with 10 -15 tolerance. The following solution is found: (82a)