On the Stability of IMEX Upwind gSBP Schemes for 1D Linear Advection-Diffusion Equations

A fully discrete energy stability analysis is carried out for linear advection-diffusion problems discretized by generalized upwind summation-by-parts~(upwind gSBP) schemes in space and implicit-explicit Runge-Kutta~(IMEX-RK) schemes in time. Hereby, advection terms are discretized explicitly while diffusion terms are solved implicitly. In this context, specific combinations of space and time discretizations enjoy enhanced stability properties. In fact, if the first and second-derivative upwind gSBP operators fulfill a compatibility condition, the allowable time step size is independent of grid refinement, although the advective terms are discretized explicitly. In one space dimension it is shown that upwind gSBP schemes represent a general framework including standard discontinuous Galerkin~(DG) schemes on a global level. While previous work for DG schemes has demonstrated that the combination of upwind advection fluxes and the central-type first Bassi-Rebay~(BR1) scheme for diffusion does not allow for grid-independent stable time steps, the current work shows that central advection fluxes are compatible with BR1 regarding enhanced stability of IMEX time stepping. Furthermore, unlike previous discrete energy stability investigations for DG schemes, the present analysis is based on the discrete energy provided by the corresponding SBP norm matrix and yields time step restrictions independent of the discretization order in space since no finite-element-type inverse constants are involved. Numerical experiments are provided confirming these theoretical findings.


Introduction
Advection-diffusion problems describe a wide range of phenomena arising due to the combination of two different physical processes, i.e., advective transport of a quantity due to a given velocity field and movement of particles from areas of higher to lower concentration.For instance, advection-diffusion equations arise as a model for the time evolution of chemical or biological species in a flowing medium.If the resulting partial differential equations are discretized on refined spatial grids, the presence of diffusion fluxes usually generates increasingly stiff problems.Particularly, explicit time discretization of semi-discrete advection-diffusion problems often results in an increasingly severe grid-dependent time step restriction since the time step Δt scales with the grid length scale Δx as Δt = O(Δx 2 ) .However, purely implicit time discretization requires the solution of large nonlinear systems of equations.Therefore, hybrid time integration schemes such as implicit-explicit (IMEX) methods have frequently been considered.Hereby, a common approach is to discretize the advective terms explicitly, while the diffusion terms are treated implicitly, see, e.g., [4,5,9,37].The additional complexity of a partially implicit treatment is justified by the observation that even though many applications contain nonlinear advection terms, the diffusion terms are often linear, only demanding the iterative solution of linear algebraic systems which are positive definite, symmetric, and sparse.In addition, the nonlinear systems resulting from potentially nonlinear diffusion terms are generally more efficient to solve than those arising from nonlinear advection terms.
While this approach alleviates the extremely restrictive time step scaling for the diffusion terms, it does not a priori dispose of the classical CFL-type time step restriction of the form Δt = O(Δx) due to the explicit time discretization of the advection terms.How- ever, for specific combinations of space and time discretizations, enhanced stability properties may be achieved.For instance, Calvo et al. [9] designed specific IMEX Runge-Kutta (IMEX-RK) methods which guarantee grid-independent time step restrictions of the form Δt = O(c∕a 2 ) for Fourier spatial discretizations, where a and c denote the advection and diffusion coefficients, respectively.Wang et al. [38] proved a similar unconditional L 2 -sta- bility result in the context of spatial discretization by the discontinuous Galerkin (DG) [12] method applied to linear advection-diffusion equations, with an extension to nonlinear problems in [39] and to the multi-dimensional case in [40].Hereby, diffusion terms were discretized by the local discontinuous Galerkin (LDG) scheme [11].In [27], an extension was provided for DG methods in one space dimension incorporating a larger class of diffusion schemes, namely, the ( , ) family [20,21] of diffusion discretizations which also includes the Bassi-Rebay schemes [6,7].Furthermore, Fu and Shu [14] proved the grid-independent L 2 -stability for the discretization of the diffusion terms by the embedded DG method and in [41], the direct DG method was considered.A related approach which provides the unconditional stability even if part of the problem is discretized explicitly was taken in [32,33] designing unconditionally stable IMEX linear multistep schemes, where both the IMEX splitting and the employed multistep schemes fulfill specific properties.
In this work, upwind generalized summation-by-parts (gSBP) space discretizations for advection-diffusion problems will be considered with respect to the combination with IMEX time integration and regarding stable time step choices.These upwind gSBP schemes do not necessarily provide a replacement for the DG space discretization, but rather a broader framework.
This general framework of SBP schemes provides structure-preserving spatial discretizations also including DG type methods.SBP operators fulfill a mimetic property, in particular mimicking the analytical concept of integration-by-parts. Originally, discrete derivative operators with an SBP property have been introduced for finite difference schemes for various types of PDEs including advection-diffusion equations and the linearized compressible Navier-Stokes equations, see [10,19,25,31,35,36]. Extensions of the SBP methodology have been provided for tensor-product grids on curvilinear elements [29] as well as for the general multidimensional case [18] including simplex elements.Via L 2 energy estimates, the spatial SBP operators automatically yield sta- ble semi-discrete schemes for periodic solutions of a broad class of linear equations.In addition, they have substantially profited from a combination with weakly enforced boundary conditions, most prominently using simultaneous approximation terms (SATs) which were first developed in [10].More recent investigations also focus on SBP operators within various popular classes of numerical schemes, e.g., finite volume schemes on unstructured dual grids [24], DG schemes with Legendre-Gauss-Lobatto nodes [15], DG schemes on Legendre-Gauss nodes both in one space dimension and on tensor-product grids [26], and flux-reconstruction schemes [28].In fact, the definition in [30] allows to derive gSBP properties for nodal DG schemes in one space dimension on arbitrary nodal sets, as long as their corresponding quadrature rule is sufficiently exact to mimic integration-by-parts at the discrete level.An extension of the SBP framework departing from their central character was originally proposed by Mattsson [23], where upwind SBP schemes of finite difference type were constructed.These schemes consist of dualpair SBP operators with non-central difference stencils.The incorporated upwind differencing leads to a built-in artificial dissipation in the semi-discrete energy analysis of periodic problems and has been used to construct upwind SBP schemes based on flux splitting for the shallow water equations in [22] and for general scalar hyperbolic conservation laws in [34].
Connections between nodal DG schemes and SBP methods have mainly fostered the discretization of conservation laws in the skew-symmetric form to preserve secondary quantities and to move towards semi-discrete energy or entropy stability results.In this context, the SBP properties of cellwise DG discretizations at the local level of one element represent the key aspect.In this work, emphasis is put on the global formulation of DG schemes in one space dimension which may be regarded as upwind gSBP schemes for a family of numerical flux functions ranging from upwind to central ones.In addition, we will take advantage of the connection to upwind gSBP schemes to improve the fully discrete stability estimates pertaining to DG and upwind gSBP schemes combined with the IMEX time integration.
While discretizing advection terms with a DG scheme is straightforward, various approaches to discretize diffusion terms have been introduced in the literature, since the discretization of higher order spatial derivatives is less natural.Either specially designed penalty terms are introduced as in [2,8] or the advection-diffusion equation is rewritten into a system of first-order equations using auxillary variables for the solution derivatives as in [3,6,7,11].The first method by Bassi and Rebay [6], termed the BR1 scheme, has actually been the first extension of the DG scheme to the compressible Navier-Stokes equations.It is based on rewriting the viscous terms into a larger, extended first-order degenerate system of PDEs with the gradient as a new unknown.After this reformulation, the standard DG approach is applied to the extended system which necessitates to prescribe two types of numerical fluxes.Hereby, the BR1 scheme is the simplest approach, using arithmetic means for both types of fluxes.Cockburn and Shu [11] analyzed various methods based on the reformulation into a first-order PDE and derived conditions on the numerical fluxes to guarantee the stability, the convergence, and a suboptimal error estimate of order N when using an approximation space of polynomial degree N. The analysis shows the sub-optimal convergence of the BR1 scheme for odd N, while the choice of alternating numerical fluxes usually associated with the LDG scheme by Cockburn and Shu [11] leads to the optimal convergence of order N + 1 .Being simple to code, parameter-free, and generic for nonlinear viscous fluxes and arbitrary grids, the BR1 scheme is still very popular and satisfies certain structural properties.In fact, in [17], the neutral behaviour of BR1 with respect to artificial dissipation over element interfaces and the resulting stability for the compressible Navier-Stokes equations has been proven.
The energy stability analysis in [27] for linear advection-diffusion problems discretized by the DG scheme in space using a ( , ) diffusion discretization has shown that the sim- ple BR1 approach coupled with upwind fluxes for advection does not inherit the enhanced stability properties in the case of IMEX time integration.In this regard, this current work provides a remedy, since the DG scheme based on BR1 fluxes may be regarded as a specific second-derivative upwind gSBP scheme compatible with central fluxes for advection.Furthermore, unlike previous discrete energy stability investigations for DG schemes, the present analysis is based on the discrete energy provided by the corresponding SBP norm matrix and yields time step restrictions independent of the discretization order in space, since no FE type inverse constants are involved.
This paper is organized as follows.Section 2 introduces the properties of first-and second-derivative upwind gSBP operators and reviews their connection to DG schemes on the global level.Furthermore, some properties and estimates required for the fully discrete stability analysis are provided.In Sect.3, an energy stability analysis based on the discrete energy defined by the associated upwind gSBP norm matrix is carried out for first-and second-order IMEX-RK time integrators combined with upwind gSBP schemes of arbitrary order in space.Numerical results verifying the conducted theoretical analysis are presented in Sect. 4. Finally, Sect. 5 contains a conclusion and an outlook on future work.

Discontinuous Space Discretization of the Linear advection-diffusion Equation
Throughout this work, we consider the linear advection-diffusion equation: with the advective velocity a > 0 and the diffusion coefficient c > 0 , supplemented by the periodic initial condition u(x, 0) = u 0 (x) in L 2 ( ) and periodic boundary conditions.For space discretization, the spatial domain Δx ⩾ ∈ ℝ + for Δx = max i Δx i under grid refinement Δx → 0 .On the closure Ī of the reference interval I = (−1, 1) , a nodal set of solu- tion points is specified and transferred to each grid cell, not necessarily including the cell boundaries in the nodal set.Furthermore, the discontinuities space discretization allows for a discontinuous representation of the semi-discrete solution at cell interfaces. (1)

Discontinuous Galerkin (DG) Formulation
DG schemes allow for piece-wise smooth approximate solutions with potential discontinuities across element interfaces.Due to its flexibility and generality, the DG scheme is a popular numerical method in a variety of applications.The main advantages are its local conservation property, an arbitrarily high order of accuracy, and superconvergence capabilities.
Dispensing with global continuity of the approximate solution, the DG approach results in relatively compact stencils which greatly facilitates both hp-adaptivity of the method and its implementation in parallel hardware environment.The DG scheme is built upon a variational formulation of the advection-diffusion equation.The basis functions and test functions used to define the DG scheme in one space dimension are taken from the finite element space , where P N (I i ) denotes the space of polynomial functions on I i of degree at most N.As usual in DG schemes, the functions in V h may be discontinuous across element boundaries.
For the discretization of diffusion terms, the PDE (1) will be reformulated as a first-order system of equations.Thus, the structure of the corresponding DG scheme is similar to the case of purely hyperbolic problems.The DG formulation for the linear diffusion equation will then be given in Sect.2.2.2.To formulate the DG scheme in matrix-vector notation similar to SBP schemes, it is sufficient to consider a scalar hyperbolic conservation law in one space dimension given by The DG scheme constructs an approximation u h of u which is piecewise continuous of the form using polynomial basis functions L i j ∈ P N I i .Here, nodal DG schemes are considered, where the basis functions are the Lagrange polynomials corresponding to a set of solution points   ∈ Ī,  = 1, ⋯ , N + 1 which are the nodes k of a suitable quadrature rule on the closure of the reference cell (−1, 1) .We assume that the quadrature rule is sym- metric and exact for polynomials of degree 2N − 1 and denote the corresponding weights by , = 1, ⋯ , N + 1 .More precisely, we consider the Lagrange polynomials defined by L i j (Λ i ( )) = j , where j denotes the Kronecker delta and Λ i transforms the reference cell I to the specific sub-interval I i , i.e., Λ i ( ) = 2 .The DG scheme in a weak form, obtained by multiplying (2) by a test function L i k , integrating in space over a cell I i , and applying partial integration, is given by denoting the values of a consistent numerical flux function f * evaluated at the left-and right-hand side limits of u h at the cell boundaries x i and the flux polynomial f h (x, t)� I i = ∑ N+1 j=1 f (u i j )L i j (x) , based on the pointwise values u i j = u h Λ i ( j ), t . (2) The integrals in (4), are then solved numerically by the chosen quadrature rule which exactly evaluates the last integral due to the degree of exactness.
The DG scheme in a strong from is obtained by a second partial integration of (4) resulting in the variational formulation Using the matrix-vector notation as shown, e.g., in [1,26], (5) rewrites in a simpler form when transformed to the reference cell I.For this purpose, we define the matrices M and Ŝ by where L k ∶ Ī → ℝ are the Lagrange polynomials corresponding to the quadrature nodes k .From the Lagrange interpolation property L j ( l ) = jl , we directly obtain hence M is diagonal.Using the solution vector u i and the vector of flux values f i , as well as the transformation of f h (x, t)| I i to the reference cell by f i h ( , t) = ∑ N+1 j=1 f (u i j )L j ( ) , the abbreviations f * ,i (1) = f * i+1 and f * ,i (−1) = f * i and the collection of the basis functions in the vector valued function L( ) = (L 1 ( ), ⋯ , L N+1 ( )) T , ( 5) is equivalent to the matrix- vector formulation see [1,26].By definition, the matrix M is invertible.Multiplying by the inverse of M , we obtain the final matrix-vector formulation of the DG scheme where the entries of D = M−1 Ŝ can be derived from (6) and are given by Djk = L � k ( j ).
Remark 1 It is worth noting that the matrix D provides a local first-derivative gSBP opera- tor on a specific DG cell, see, e.g., [26], but outsources the interface fluxes to be dealt with similar to the weak boundary treatment of SATs.In contrast, the next section investigates the upwind gSBP structure of the global DG scheme including the presence of numerical fluxes on cell interfaces. (5)

Upwind gSBP Space Discretization
Upwind SBP schemes, first derived in [23], are a class of finite difference methods which combine the structure-preserving properties of SBP derivative operators with artificial dissipation by choosing upwind directions.These schemes consist of dual-pair SBP operators with non-central difference stencils which lead to a built-in artificial dissipation when combined with flux splitting.While traditionally, SBP schemes were set up on equidistant grid points, we will consider generalized versions in this work on non-uniform nodal sets not necessarily containing boundary points which have been termed gSBP schemes.These schemes may be further modified by including upwind stencils, resulting in upwind gSBP schemes which will be considered in this work.The common goal of schemes within the SBP framework is to transfer continuous energy estimates to the semi-discrete level.For instance, considering solutions u(x, t) of the linear advection equation: subject to the homogeneous boundary condition u(x a , t) = 0 , the continuous energy estimate can be proven.By the SBP property, this estimate may be transferred to the spatial discretization.The corresponding operators approximate the first derivative x in (8) and are based on a set of grid points {x i } 1⩽i⩽M ⊂ [, ] called a nodal set.On this nodal set, approximate solutions of ( 8) are given by time-dependent solution vectors and u x in ( 8) is replaced by the term D − u , utilizing an upwind gSBP operator D − as defined below.
Definition 1 A pair of difference operators denoted by D + and D − , both approximating the first derivative x on an interval ( , ) , are called upwind gSBP operators of degree q if i. the matrices D + and D − provide accurate approximations to x of degree q, i.e., where T is the representation of the monomials x k on the grid points; ii.there exists a positive definite symmetric matrix M , also called a norm matrix, such that where B is a boundary operator of the form B = t t T − t t T , with boundary interpo- lation provided by t T x l = l and t T x l = l for 0 ⩽ l ⩽ r with r ⩾ q; iii.integration by parts is mimicked by the property iv.as an additional stability constraint, the symmetric matrix (8) Remark 2 Some remarks on the modifications with respect to classical SBP operators and on naming conventions are in order.
i. Due to the generalization provided by gSBP schemes beyond the traditional case of equidistant grid points including the boundary points and , the boundary operators of gSBP schemes extend the traditional boundary operator B = diag(−1, 0, ⋯ , 0, 1) to the given nodal set.ii.In the case C = 0 , we have D − = D + which is denoted a classical gSBP operator, since it does not feature any upwind directions.iii.Upwind gSBP operators based on diagonal norm matrices M are denoted diagonal- norm upwind SBP operators.
The motivation for using upwind gSBP operators is the additional artificial damping introduced by enforcing the negative semi-definiteness of the matrix C in (12).We will illustrate this property of an upwind gSBP scheme in comparison with classical gSBP schemes in the following example.
Considering the spatial discretization of the linear advection equation ( 8) via an upwind gSBP scheme, we use the operator D − , whereas for a < 0 , the operator D + would be applied.We obtain where the right-hand side of the above equation contains the required simultaneous-approximation-term (SAT), dealing with the boundary condition at the left boundary.Thereby, the corresponding SAT parameter ∈ ℝ is specified, such that a discrete energy estimate is fulfilled which mimics the continuous case, e.g., as in estimate (15).
Multiplying (13) from the left by u T M and adding the transpose now yields Using the upwind gSBP properties ( 11) and ( 12), we have and, therefore, setting u = t T u, For ⩽ − a 2 , the time evolution of the discrete energy may now be estimated by where the term containing the negative semi-definite matrix C introduces additional arti- ficial dissipation compared to the continuous estimate (9), and this additional dissipation vanishes for classical gSBP operators without upwind character, since C = 0 in this case.The upwind gSBP framework also enables the construction of second-derivative operators of the form which approximate the second derivative

Global DG Discretization of Advective Terms in Upwind SBP Framework
While the elementwise DG discretization of the linear advection equation results in a classical gSBP scheme on each grid cell, the corresponding global DG formulation over all grid cells results in its upwind variant.The classical SBP properties of element-wise 1D-DG schemes on general nodal sets are by now well-known, see, e.g., [15,26,28,31].In fact, the cellwise formulation (7) represents a gSBP scheme with the first-derivative gSBP operator D = M−1 S which approximates the continuous derivative to degree q = N and the degree of B is r = q = N , see [26, Lemma 1].
The boundary operator B is given by , hence, for a single cell (x i , x i+1 ) , we set = x i , = x i+1 and obtain t = L(−1) and t = L(1).
Next, we consider the global DG formulation over all grid cells.According to (7), on interior cells I i , i = 2, ⋯ , K − 1 , the cell-wise DG scheme in matrix-vector formulation applied to the linear advection equation (8), with the positive advection velocity a > 0 , is given by Hereby, a suitable numerical flux function (au) * ,i is employed, where (au) * ,i (−1) = (au) * i and (au) * ,i (1) = (au) * i+1 with Therefore, the choice of numerical flux functions ranges between the central flux for = 0 and the upwind flux for = 1 2 .Using the SAT boundary treatment with the SAT parameter , the scheme on the left boundary cell is given by while on the rightmost DG cell with outgoing information through the right boundary, we have ( 16) We now rewrite the global DG discretization on the computational domain as an upwind gSBP scheme which includes the interactions of degrees of freedom on adjacent cells.Under this assumption, we may rewrite (17) as with the global nodal DG solution vector given by u = (u 1 , ⋯ , u K ) T and the extended SATs Now, the global DG derivative operator D − is a block tridiagonal matrix with blocks corresponding to each DG element and its left and right adjacent cells.More precisely, we have where the repeating blocks corresponding to interior cells and the left and right boundary blocks are given by where D denotes the cellwise first-derivative gSBP operator of the DG scheme (7).
Proof The evaluation of the numerical flux function at an interface given in (18) may be rewritten as and the evaluation of the flux function au i h at cell boundaries is obtained as The desired block structure is obtained by inserting these contributions into (17) for interior DG cells and into (19) and (20) for boundary DG cells.Hereby, we account for the dependencies of the degrees of freedom on cell i only on u i , u i−1 , and u i+1 and set up the blocks accordingly.
The Proof We define the matrices Q − and Q + in (10) by Next, we show that Q − and Q + satisfy the SBP property (11), i.e., that the dual operator D + is chosen properly.We have where Q + + (Q − ) T = 0 , since the block representation given in Proposition 1 yields It remains to show the stability constraint (12), i.e., to show that is negative semi-definite.We have with Hence, we obtain i.e., C is negative semi-definite.
The following remark provides additional insight into the upwind gSBP character of the global DG formulation with numerical fluxes ranging from upwind to central fluxes.(

Periodic first-derivative upwind gSBP operators
Since the energy stability analysis in Sect. 3 regarding fully discrete IMEX upwind gSBP schemes is carried out for periodic advection-diffusion problems, the corresponding semi-discretization of (1) by upwind gSBP operators is devoid of SATs and periodic upwind gSBP operators are required which are of a slightly modified form, given by The construction of these periodic global discrete derivative operators is only based on the interior cell discretization (17).Therefore, they use the same building blocks given in Proposition 1 except for the absence of the boundary blocks A lb ( ), A rb ( ).

Rewriting LDG and BR1 Schemes as Second-Derivative Upwind SBP Operators
The upwind gSBP operators resulting from the DG discretization of the linear advection equation in one space dimension can be used to construct specific second-derivative operators.In particular, the construction (16) of second-derivative upwind gSBP operators by combining D − ( ) and D + ( ) is directly related the DG technique of rewriting diffusion equations into a first-order system using an auxiliary variable for the solution derivative as in [3,6,7,11].For the linear heat equation this procedure provides the system using an auxiliary variable q.For periodic boundary conditions, only the interior cell discretization ( 7) is relevant and the element-wise strong form DG semi-discretization yields Regarding the specification of the numerical fluxes, at each element boundary, the lefthand side and right-hand side values of a piecewise continuous function v are denoted by v − and v + , respectively.The corresponding jump at element interfaces is denoted by [v] = v + − v − and the arithmetic mean is given by {v} = 1 2 (v + + v − ) .The simplest choice for q * and u * is the BR1 scheme [6] given by the arithmetic means Furthermore, the original LDG [11] yields a parameter-dependent family of diffusion fluxes which, in one space dimension, contains the BR1 approach as a specific case with c 11 = c 12 = 0 .Commonly, LDG fluxes refer to the choice of alternating fluxes with c 12 = ±1, c 11 = 0 , which offers two different variants.One implementation is thus given by the second variant is specified by For the LDG a , the LDG b , and the BR1 schemes, the numerical fluxes q * and u * can now be written as for LDG b and q = u = 0 for BR1.Due to the assumed periodic boundary conditions, we now use the periodic global upwind gSBP operators (24), (25), and (26) developed in Sect.2.2.1.
Thus, the LDG a scheme can be rewritten as where u = u 1 , ⋯ , u K T , q = q 1 , ⋯ , q K T denote the global nodal solution vectors and the periodic upwind gSBP operator D − = D − = 1 2 is taken from (24).Hence, for the LDG a scheme, we have Analogously, for the LDG b variant we obtain Both of these formulations correspond to true upwind gSBP operators approximating the second derivative as defined in (16).In contrast, for the BR1 scheme, we have (29) q * ,BR1 = {q}, u * ,BR1 = {u}.
where D = D − ( = 0) = D + ( = 0) is a classical first-derivative gSBP operator, as already mentioned in Remark 3. Thus, the BR1 scheme represents a second-derivative gSBP operator which is obtained by applying the first-derivative gSBP operator twice.

Preliminary Estimates for First-and Second-Derivative Upwind gSBP Operators
The fully discrete stability analysis for IMEX upwind gSBP schemes in Sect.3.1 and Sect.3.2 will be restricted to periodic problems and is based on the discrete energy norm induced by the diagonal positive definite norm matrix M .Defining the inner product a discrete energy norm is given by Furthermore, the Cauchy-Schwarz inequality will be useful.As already mentioned, we assume periodic boundary conditions and the semi-discretization of the linear advection-diffusion equation ( 1) results in based on the periodic upwind gSBP operators given in ( 24) and (26).Using this notation, we have the following estimates and relations between the operators D ± and D 2 , which will be required for the subsequent discrete energy stability analysis.
Proposition 2 For the periodic upwind gSBP operators defined in (24), (25), and (26), and arbitrary nodal values u, v we have i.dissipativity of −D − , i.e., ii.dissipativity of D 2 , more precisely iii. a Cauchy-Schwarz-type inequality Proof First, making use of the periodicity, the upwind gSBP property (11) reduces to the equation: leading to (D − ) T M = −MD + which we will frequently use in the following derivations.Due to the symmetry of the inner product, we have and the first assertion is thus proven by Regarding the second assertion, we obtain (34) by and ( 35) directly follows.For (36), we may rewrite and the third assertion follows from the Cauchy-Schwarz inequality for the inner product (⋅, ⋅) M .
We will, furthermore, use the Young's inequality in the form for any constant C > 0 and any a, b ∈ ℝ.

IMEX Schemes for Advection-Diffusion Problems
Implicit-explicit one-step time integration methods for advection-diffusion equations are generally based on partitioned Runge-Kutta schemes applied to a partitioned system of ordinary differential equations of the form Hereby, different RK schemes, i.e., an explicit and an implicit one, are applied to the components F 1 and F 2 of the split right-hand side. (36) As in [4,9] we consider (s + 1)-stage explicit Runge-Kutta methods coupled with implicit s-stage DIRK schemes, where the absissae c 1 , c 2 of the explicit and implicit schemes, respec- tively, fulfill c 1 = 0 c 2 .The DIRK scheme is then recast into an (s + 1)-stage scheme as well by padding the first row and first column with zeros.We thus obtain the Butcher table Applied to the system of ODEs (38), the IMEX-RK scheme then has the form where u n , u n+1 denote the approximations at times t n and t n+1 = t n + Δt and u (i) are the inter- mediate stage values corresponding to the intermediate times t n,i given by t n,i = t n + c i Δt In [9], Calvo et al. derived necessary conditions on the coefficients of IMEX schemes of the above form to obtain a time step restriction of Δt = O(c∕a 2 ) for advection-diffusion equations discretized by the Fourier method.These conditions are similar to L-stability and include a 2 s+1,j = b 2 j , j = 2, ⋯ , s + 1 as a requirement for the implicit scheme.In this work, for the L 2 -stability analysis of IMEX-DG schemes, we consider the follow- ing first-and second-order IMEX schemes also used in [38]. First-order: Second-order:

fulfilling the identity
In addition, some of the numerical experiments carried out in Sect. 4 utilize the third-order IMEX scheme taken from [9] which is given by with

Stability Analysis for the First-Order IMEX Scheme
In principle, the energy stability analysis for IMEX upwind gSBP schemes is similar to the L 2 -stability analysis carried out for LDG schemes in [38] and for the ( , )-family in [27].The most significant difference of the current approach is the direct use of the discrete energy generated by the upwind gSBP norm matrix M instead of the usual L 2 -norm of the approximate solution.This leads to a very transparent energy estimate which is devoid of finite-element-type inverse constants and holds for upwind gSBP schemes of arbitrary order in space.In particular, the provided time step restriction does not depend on the polynomial degree in space.For the applicability of the following stability analysis, compatibility of the first-derivative and second-derivative upwind gSBP operators in the sense D 2 = D − D + is crucial.Discretizing the linear advection-diffusion equation ( 1) with a > 0 by first-and sec- ond-derivative upwind gSBP operators D − and D 2 in space and the first-order IMEX scheme (40) as time integrator yields the fully discrete IMEX upwind gSBP scheme Applying (⋅, ⋅) M on both sides of (43) with u n+1 as the second argument, we get Denoting u = u n+1 − u n and using identities analogously to the derivation in [38] com- bined with the estimate (33), we obtain Inserting this into (44) and employing (35), (36), and finally the Young's inequality (37), we have (42) Summarizing the above findings we have the following theorem.It is worth noting that the assertion of Theorem 2 is quite general, since it also holds for the LDG diffusion scheme coupled with upwind fluxes for advection, i.e., the scheme analyzed in [38].In addition, the constant in the above time step restriction is specifically determined and is independent of the polynomial degree of the space discretization, since no detour through finite-element spaces for u and u x is required.

Stability Analysis for the Second-Order IMEX Scheme
Discretizing the linear advection-diffusion equation ( 1) by upwind gSBP operators in space as in Sect.3.1 and using the second-order IMEX time integrator (41) yields the fully discrete IMEX upwind gSBP scheme Employing the discrete inner product (⋅, ⋅) M on both sides of (45) with u (n,1) as the second argument and on both sides of (46) with u n+1 as the second argument, we have where we have used − = 1 .Setting u ,1 = u (n,1) − u n and u ,2 = u n+1 − u (n,1) and adding (47) and (48), we obtain

Numerical Results
In this section, the stability and accuracy of IMEX upwind gSBP schemes are compared with respect to different combinations of the first-derivative and second-derivative upwind gSBP operators.All upwind gSBP operators in this section are based on uniform grids and the Legendre-Gauss-Lobatto quadrature rule.c Δt max , where Δt max is the maximum time step ensuring a non-increasing discrete energy (0, 0) First-order IMEX time integration

Numerical Study on Time Step Restrictions
We consider the exact solution of the linear advection-diffusion equation (1) in the interval (x a , x b ) = (−π, π) .The advec- tion-diffusion problem is discretized in space by different upwind gSBP operators on uniform grids using a number of N + 1 Legendre-Gauss-Lobatto nodes on each grid cell.The implemented upwind gSBP operators are specified by the parameters and , respectively, where we choose D − = D − ( = ) and D 2 = D − ( )D + ( ) .Thus, for = , the first-derivative and second-derivative operators are compatible in the sense of Theorems 2 and 3.For the choice represents upwind advection fluxes paired with LDG diffusion fluxes, while , = (0, 0) denotes the pairing of central advection fluxes with the BR1 diffusion scheme.The firstorder IMEX scheme (40) and the second-order IMEX scheme (41) are then used to discretize advection terms explicitly and diffusion terms implicitly in time.From the theoretical analysis, we expect the schemes to be stable for time steps Δt ⩽ c a 2 , where is some constant independent of grid refinement.
For different advection and diffusion parameters a and c, Table 1 lists the maximum stable time steps obtained by the IMEX schemes of first-and second-order coupled with upwind gSBP operators using N ⩽ 3 .In the numerical setup, we successively double the number of cells K and compute the numerical solution until the final time T = 1 000 is reached.The maximum stable time step Δt max is determined as the maximum time step for (50) u(x, t) = e −ct sin(x − at) Table 2 Values of = a 2 c Δt max for IMEX upwind gSBP schemes using the third-order IMEX method (42) which the discrete energy norm (u, u) M of the numerical solution is non-increasing and the corresponding values of = a 2 c Δt max are listed.An entry of + indicates that arbitrarily large time steps may be chosen without increasing the discrete energy in time.
As predicted by the theoretical analysis, for all schemes using compatible first-derivative and second-derivative upwind gSBP operators with = , the maximum admissible time step is bounded from below under grid refinement.The precise bound varies for different advection-diffusion parameters a, c but fulfills the theoretically predicted time step restriction as long as the compatibility condition is fulfilled.Moreover, the different schemes possess nearly the same time step bounds under the same setup with respect to the parameters a, c and the chosen IMEX scheme.Regarding the enhanced IMEX stability property, the time step bound is lower for second-order IMEX time integration as predicted by the theoretical analysis.
On the other hand, as already noted in [27], the choice of , = 1 2 , 0 repre- senting the BR1 diffusion scheme paired upwind advection fluxes does not admit gridindependent time step sizes, in accordance with the theoretical analysis with only covers the case = as already stated.Instead, we observe admissible time steps of = O(Δx) , analogous to the time step restrictions for explicitly discretized advection equa- tions.Replacing the upwind advection fluxes by central fluxes is compatible with the BR1 scheme and yields the enhanced IMEX stability which can also be observed in the results of Table 1.
Regarding higher order time discretization, Table 2 lists the corresponding maximum time steps ensuring a non-increasing discrete energy when combining the previously considered spatial upwind gSBP schemes for N ⩽ 3 with the third-order IMEX scheme (42).Although the theoretical analysis in Sect. 3 only covers the first-order and second-order IMEX schemes (40) and (41), respectively, an analogous stability behavior can be observed for the third-order IMEX scheme.In addition, for this setup, compatible first-and secondderivative upwind gSBP operators yield the enhanced IMEX stability, while for the combination of upwind advection fluxes with BR1 diffusion, i.e., , = 1 2 , 0 , the clas- sical CFL condition for the pure advection equation needs to be enforced.

Numerical Study on Accuracy of IMEX Upwind gSBP Schemes
In this section, we determine the experimental order of convergence (EOC) of specific IMEX upwind gSBP schemes, particularly for coarse time grids.Regarding the designed spatial and temporal order of the schemes, we study the cases N ⩽ 3 and utilize the IMEX schemes of the second order (41) and the third order (42).
First, we reconsider the decreasing exact solution (50) of the linear advection-diffusion equation (1).Table 3 lists the L 2 -errors and the derived EOC for the various schemes applied to this problem.Here, we compute the numeral solution for two different parameter sets of a = c = 0.1 and a = 1, c = 0.1 until the final time T = 10 and use time steps Δt = Δx for different constants as indicated in the table.In most cases, the predicted convergence rate is experimentally confirmed as the minimum of the expected orders in space and time.One exception is the poor performance in the case of , = 1 2 , 0 , i.e., upwind advection fluxes paired with the BR1 diffusion scheme as this choice suffers from the predicted loss of enhanced IMEX stability for larger time steps, particularly in the third-order case.A second exception is the case of , = (0, 0) , which yields a stable scheme but is subject to order reduction regarding the formally second-order scheme.As a matter of fact, this type of order reduction is a known phenomenon which is prominent in DG schemes employing numerical fluxes of central type in the case of odd polynomial degrees, also observed for the central-type kinetic energy preserving numerical fluxes for the compressible fluid flow, e.g., in [16,26].Furthermore, as in [38], we now supplement the linear advection-diffusion equation (1) by a source term, i.e., we consider the modified advection-diffusion problem The source term is set to g(x, t) = e ct (2c sin x + cos x) to provide the exact solution u(x, t) = e ct sin x , which is exponentially growing in time.
Computations are carried out for a diffusion parameter of c = 0.1 until the final time T = 10 with time steps Δt = Δx with as indicated in Table 4. Again, the anticipated order of convergence is reached in all setups with two exceptions.Analogously to the exponentially decreasing solution, order reduction occurs for the central scheme , = (0, 0) for N = 1 .Furthermore, the pairing of upwind advection with BR1 diffusion replicates the loss of the enhanced IMEX stability for all but the first time step choice in Table 4.

Numerical Results for the Viscous Burgers' Equation
Finally, we study the behavior of IMEX upwind gSBP schemes with respect to the viscous Burgers' equation as a prototype containing nonlinear convective terms, supplemented by the initial condition and periodic boundary conditions.Figure 1 shows the initial solution as well as the numerical solution for c = 0.1 at a given output time T. Two different upwind gSBP schemes with N = 2 are employed within the following semi-discrete formulation of the viscous Burgers' equation, with the vector of flux values f = 1 2 u 2 and the matrix C as in (12).This semi-discrete form is similar to the one used in [34] for the inviscid Burgers' equation but is based on the equation in a conservative form instead of the skew-symmetric formulation employed in [34].If D − is constructed from a DG scheme using upwind advection fluxes, the semi-dis- cretization (52) is equivalent to the use of global Lax-Friedrichs (LF) fluxes for the nonlinear convective term and will be labeled accordingly.If D − is constructed from central fluxes, we have C = 0 and thus obtain a central discretization of the inviscid fluxes.For time integration, the second-order IMEX scheme (41) is applied with a time step size of Δt = 0.1 .This same time step is used on two different grids with K = 50 and K = 100 ele- ments, respectively.First, we compute approximate solutions with the upwind gSBP scheme consisting of BR1 diffusion paired with global LF fluxes, i.e., the first of the above mentioned variants for D − .Since for linear advection, global LF fluxes reduce to upwind fluxes, the global LF fluxes can be viewed as a generalization of upwind fluxes to the nonlinear case.For this setup, the simulations needs to be stopped at T = 1.8 and T = 1 due to the expected lack of the stability, as shown in the top row of Fig. 1.These results indicate that also for nonlinear problems, a pairing of the BR1 scheme with upwind-type fluxes (51) might be less favorable in terms of the stability in the context of IMEX advection-diffusion splitting.A remedy is provided in the bottom row of Fig. 1 using central fluxes together with BR1.The corresponding stable numerical solution is shown at output time T = 2 .This solution is visually indistinguishable from a reference solution computed on a fine grid of K = 1 000 elements using the upwind gSBP scheme defined by and N = 3 combined with third-order IMEX time integration using a small time step of Δt = 0.01 .In addition, refining the grid does not necessitate choosing smaller time steps for this combination of central first-derivative and second-derivative gSBP operators.

Conclusion and Outlook
The current fully discrete energy stability analysis for linear advection-diffusion problems discretized by upwind gSBP schemes in space and IMEX-RK schemes in time based on advection-diffusion IMEX splitting demonstrates that if the employed upwind gSBP operators are chosen to fulfill D 2 = D − D + , the allowable time step size is inde- pendent of grid refinement, although the advective terms are discretized explicitly.In one space dimension, upwind gSBP schemes represent a general framework including standard DG schemes, thus the present results also hold for the LDG scheme combined with upwind advection fluxes.Furthermore, the analysis in this work is based on the discrete energy provided by the corresponding SBP norm matrix and does not involve inverse constants depending on the discretization order in space; therefore, the resulting time step bound is independent of the polynomial degree of the approximate solution.
While previous work for DG schemes has demonstrated that the combination of upwind advection fluxes with the BR1 diffusion scheme does not provide the unconditional stability in the sense of grid-independent time step restrictions, the current work shows that central advection fluxes are compatible with BR1.However, for this kind of overall central scheme, additional artificial dissipation might be necessary for discontinuous initial conditions and in the nonlinear convection-dominated case.These additional artificial dissipation terms could again be discretized either implicitly or explicitly within an IMEX approach.The numerical results for the nonlinear viscous Burgers' equation in this work indicate that a compatible pairing of first-derivative and second-derivative upwind gSBP operators is also required for nonlinear conservation laws with dissipative terms, if the enhanced IMEX stability is desired.This should be investigated more closely with respect to a discrete energy analysis which might require skew-symmetric formulations as in [13,34] instead of the divergence form of the given conservation law.
Although the focus of this work is restricted to one-dimensional advection-diffusion problems, the presented theoretical results on the enhanced IMEX stability can be carried over to spatial discretizations in higher dimensions as long as the upwind gSBP property is fulfilled and the corresponding first-derivative and second-derivative operators are compatible.This extension relies on the availability of such upwind gSBP operators on structured or unstructured meshes.While structured meshes can be dealt with using tensor-product grids, multi-dimensional generalized SBP operators on simplex meshes have been constructed in [18].These gSBP operators fulfill a multi-dimensional analog of the properties in Definition 1 with C = 0 , i.e., without upwind dissipa- tion.Considering for instance a two-dimensional linear advection-diffusion problem of the form subject to periodic boundary conditions, space discretization may be carried out by firstderivative upwind gSBP operators D − x and D − y discretizing x and y , respectively, and compatible second-order derivative operators D 2,x = D − x D + x and D 2,y = D − y D + y .Due to the upwind gSBP property (11), Proposition 2 then extends to the two-dimensional case and the extension of the stability results of Theorems 2 and 3 is only technical, since it involves an additive decomposition into the terms in x-and y-directions.The construction of different multi-dimensional upwind gSBP operators with C ≠ 0 and their application and com- parison in the context of practically relevant multi-dimensional advection-diffusion problems seems to be a promising approach for future exploration.This investigation would also profit from a comparison of the upwind gSBP approach to the multi-dimensional IMEX-LDG schemes considered in [40].

Theorem 1
specific form of D − allows us to show that the dual pair D − , D + with D − = D − ( ) and a suitably constructed dual operator D + satisfy the upwind gSBP properties specified in Definition 1 with = x a and = x b .The dual pair of discrete derivative operators is a dual pair of diagonal-norm upwind gSBP operators with respect to the global diagonal norm matrix and the generalized boundary operator with B l = −L(−1)L(−1) T and B r = L(1)L(1) T .

Remark 3
Since D + ( ) = D − (− ) , we note that i. the dual operator D + arises naturally from the DG scheme applied to the linear advec- tion equation with the negative advection velocity a < 0 , and ii. for = 0 , the global DG operator D = D − (0) = D + (0) is a classical diagonal-norm gSBP operator with global central character.

Theorem 2
If {D − , D + } is any dual pair of first-derivative upwind SBP operators with the norm matrix M and the second-derivative upwind SBP operator is chosen as D 2 = D − D + , then the fully discrete solution of the IMEX-upwind SBP scheme (43) satisfies the discrete energy stability of the form ‖u n+1 ‖ M ⩽ ‖u n ‖ M if the time step is bounded by

can be proven positive definite for = 1 − √ 2 2Theorem 3
and we have Using(33),(36), we can, furthermore, bound R 1 by since for the parameter = 1 − Inserting the estimates for R 1 and R 2 into (49), we obtain the following theorem.If {D − , D + } is any dual pair of first-derivative upwind SBP operators with the norm matrix M and the second-derivative upwind SBP operator is chosen as D 2 = D − D + , then the fully discrete solution of the IMEX-upwind SBP schemes (45), (46) satisfies the discrete energy stability of the form ‖u n+1 ‖ M ⩽ ‖u n ‖ M if the time step is bounded by

Table 4
Accuracy comparison of IMEX upwind gSBP schemes for the exponentially growing solution (51) with c = 0.1 .Computations carried out until the final time T = 10