domain: math info

text:
Martin Schreiber
email: mar-tin.schreiber@univ-grenoble-alpes.fr

Jed Brown

A generalized rational

Keywords: Exponential integrators, rational approximation, parallel-in-time, Cauchy contour, Butcher table, diagonalization AMS subject classifications

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1. Introduction. Time integration of IVPs is one of the most traditional tasks in scientific computing, having seen two centuries of research. The IVPs we are interested in are given entirely or partly by linear autonomous PDEs, which are ubiquitous in applications ranging from daily weather forecasting [START_REF]The Strength of a Common Goal: A Roadmap To 2025[END_REF] to full waveform inversion [START_REF] Virieux | An introduction to full waveform inversion[END_REF]. Integration of such systems is sequential in time using conventional methods such as explicit and diagonally implicit Runge-Kutta [START_REF] Runge | Über die numerische Auflösung von Differentialgleichungen[END_REF][START_REF] Kutta | Beitrag zur näherungsweisen integration totaler Differentialgleichungen[END_REF]: Without special structure [START_REF] Jackson | The potential for parallelism in RK methods. Part 1: RK formulas in standard form[END_REF], the state at each stage is necessary to compute the next stage, either explicitly or implicitly. The time step size is typically limited by stability and/or accuracy requirements and the method is purely sequential in the time dimension.

With the desire to solve PDEs with ever-higher resolutions, the demands on highperformance computers (HPC) have increased. The steady and ongoing increase in HPC performance is provided almost exclusively by increased parallelism; increasing resolution in space (spatial scalability) can be solved in the same amount of time per time step, but the wallclock time to simulate for a fixed physical duration increases due to the increasing number of time steps to satisfy the Courant-Friedrichs-Lewy (CFL) constraint [START_REF] Courant | Über die partiellen Differenzengleichungen der mathematischen Physik[END_REF] for transport phenomena. Consequently, refinement to increase accuracy on a transient physical problem is always a scaling challenge, and many applications are unable to increase spatial resolution without sacrificing external timelines such as IPCC assessment reports [START_REF] Calvin | Climate Change 2023: Synthesis Report. Contribution of Working Groups I, II and III to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF] or design/manufacturing timelines. Parallelism in the time dimension is seen as an opportunity to utilize greater parallelism to meet stringent simulation timelines. The Rational Approximation of Exponential Integration (REXI) family of methods, which we briefly explain next, are a promising candidate for hyperbolic PDEs. Consider a linear autonomous PDE given by ∂U (t) ∂t = LU(t) with U(t) the current state and L a linear differential operator. Discretizing state variable and operator in space leads to (1.1)

∂U (t) ∂t = LU (t)
with L the discrete linear operator and U (t) the discrete state variables at time t.

Solving such IVPs have been intensively studied over the last decades with various approaches, and one of the direct methods is the application of an exponential integration

U (t + ∆t) = exp(∆tL)U (t) (1.2)
with the solution U (t) at time t. We want to emphasize that no time discretization has been introduced and that the only approximations are related to space. The REXI method exploits the feature that exp(∆tL) only needs to be approximated within a spectrum related to time step size ∆t and the spectrum of L. In the present work, we can express a variety of different time integration methods by what we refer to as the "unified REXI" formulation given by

U (t + ∆t) ≈ γU (t) + N n=1 β n (∆tL -α n ) -1 U (t)
Parallelization (1.3) with the time step size given by ∆t, and (typically) complex valued REXI coefficients α n , β n and real-valued γ is a new generalization that we will use in the following sections. The remainder of this work investigates different ways to infer these REXI coefficients and their relation to a class of Runge-Kutta methods. Based on this, we will study their numerical properties in the linear and nonlinear context.

Related Work.

2.1. Exponential integration. Exponential integration methods are formulated for nonlinear systems factored as (2.1)

∂U (t) ∂t = LU (t) + N (U (t)),
where the linear part L is intended to capture the "fast" dynamics and N is the remaining nonlinear part. An exact ansatz for advancing this split equation over a

This manuscript is for review purposes only.

finite time interval is given by (2.2) U (t + ∆t) = exp(∆tL)U (t) + ∆t 0 exp((∆t -τ )L)N (U (t + τ ))dτ.

In this form, the linear parts are integrated precisely by an exponential function, hence overcoming potential stiffness challenges caused by the linear parts. Due to this advantageous property, the interest in these exponential integrators has steadily increased over the last decades (see, e.g., [START_REF] Moler | Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later[END_REF][START_REF] Hochbruck | Exponential integrators[END_REF]) where various approaches have been taken to approximate the integral of the nonlinearities. One of the most commonly known approximations of the integral is, e.g., given by (see [START_REF] Cox | Exponential time differencing for stiff systems[END_REF])

(2.3) U (t + ∆t) = ϕ 0 (∆tL)U (t) + ∆tϕ 1 (∆tL)N (U (t))

where we used the notations ϕ 0 (Z) = e Z and ϕ 1 (Z) = e Z -I Z . We skip further examples for discretized exponential integrator formulations and only like to point out the ϕ functions to be omnipresent in higher-order exponential integration methods, which is generally given, e.g., by

ϕ i+1 (Z) = (ϕ i (Z) -ϕ i (0)) Z -1 for i ≥ 0. (2.4)
An investigation of all different varieties of discretizations of exponential integrators incorporating the nonlinearities is beyond the scope of this work, and we continue on the linear parts. These linear parts can be either given by full linear PDEs or by time integrating only a part of linear PDEs where the underlying requirement of time integration results in problems of the form U (t+∆t) = ϕ 0 (∆tL)U (t) = exp(∆tL)U (t).

In contrast to state-of-the-art time integration methods, which are used in operational codes, exponential integrators for linear operators avoid any time-discretization errors. However, the computational complexity can be tremendous and triggered the development of various ways to tackle this challenge [START_REF] Moler | Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later[END_REF]. We will briefly summarize the ones recently researched, namely based on Krylov subspaces and REXI.

The exponential can be approximated using Krylov subspace solvers (see [START_REF] Niesen | A Krylov subspace algorithm for evaluating the phi-functions appearing in exponential integrators[END_REF][START_REF] Tokman | Efficient integration of large stiff systems of ODEs with exponential propagation iterative (EPI) methods[END_REF][START_REF] Tokman | A new class of exponential propagation iterative methods of Runge-Kutta type (EPIRK)[END_REF][START_REF] Clancy | On the use of exponential time integration methods in atmospheric models[END_REF]) where we see polynomial approximations (e.g., based on Chebyshev) as a subclass of them. The advantage of such methods is their simplicity -assuming the Krylov solver framework given -since only vector multiplications with the linear operator are required. However, the potential drawbacks of Krylov subspace solvers are their inherent property of sequential iterations over the Krylov subspace, hence not providing ways to exploit additional degrees of parallelization. An alternative is to use the REXI method, which will be discussed in the next section in further detail.

REXI.

The particular way to evaluate the ϕ functions in the present work is strongly related to Padé approximations, which can be used as a first instance to approximate the ϕ functions. This approximation is most naturally related to how all Runge-Kutta formulations, e.g., based on the Butcher table, can be formulated for linear autonomous operators. However, higher-order polynomials in the denominator also make the development of solvers for such Padé approximations more challenging, and a partial fraction decomposition can be used. This well-known decomposition transforms a higher-order Padé approximation into a sum of lower-order terms, which can eventually be used to develop solvers parallelized over all terms. Although not being explained in the context of Padé approximations, REXI methods using the partial fraction decomposed form have been developed in different contexts. They This manuscript is for review purposes only. can be interpreted as a Padé approximation, which is why they are mentioned in this context.

In what follows, we will provide an overview of different methods, which can all be phrased in REXI form. One of the earliest REXI formulations for hyperbolic PDE time integrators is related to the Laplace transformation (cf. [START_REF] Lynch | Filtering integration schemes based on the Laplace and Z transforms[END_REF][START_REF] Clancy | Laplace transform integration of the shallow-water equations. Part I: Eulerian formulation and Kelvin waves[END_REF]). Here, the PDE is transformed with the Laplace operator, where the backward transform is conducted with a Cauchy Contour integral. This transformation can be again related to an exponential integration scheme, namely to the Cauchy Contour method mentioned below, see also [START_REF] Trefethen | The Exponentially Convergent Trapezoidal Rule[END_REF]. More recently, time integration based on Laplace transformations with a circle-based Cauchy contour integration have been more intensively studied in [START_REF] Rojas Mendoza | Numerical Solution of Ordinary Differential equations using Laplace transform integration[END_REF] with ODEs, in particular filtering properties. However, it needed a more extensive (community) effort to develop other, e.g., higher-order methods around them, as has been extensively the case for exponential integration methods. This lack of advancements with the REXI Laplace transform is also why we only concentrate on exponential integration-based formulations. We want to point out that the same approaches could also be taken from the Laplace transform perspective.

Another way to infer REXI coefficients originates from the REXI method based on Gaussian basis functions originally developed in [START_REF] Haut | A high-order time-parallel scheme for solving wave propagation problems via the direct construction of an approximate time-evolution operator[END_REF], which only targets purely oscillatory problems (hence L has only imaginary or zero eigenvalues). This method also showed excellent properties regarding the wallclock-time vs. error for the linear shallow-water equations on the plane (see [START_REF] Schreiber | Beyond spatial scalability limitations with a massively parallel method for linear oscillatory problems[END_REF]) and on the rotating sphere [START_REF] Schreiber | A parallel time integrator for solving the linearized shallow water equations on the rotating sphere[END_REF][START_REF] Schreiber | Exponential integrators with parallel-in-time rat. approx. for the shallow-water equations on the rotating sphere[END_REF].

Although initially developed for analytical reasons, the Cauchy contour integration method can indeed be used for REXI time integration. As pointed out above, one of the first times this has been used as a REXI-like method was with the Laplace transformations. However, exponential formulations (see Eq. 2.2) provide a more direct and substantial established way to integrate in time. Here, the property of ϕ i being an analytical function plays a fundamental role in the Cauchy contour integration method as well as the exponentially fast converging trapezoidal rule to approximate the contour [START_REF] Trefethen | The Exponentially Convergent Trapezoidal Rule[END_REF]. This method has already been used in different works: The approximation of ϕ i (x) evaluations on scalar values has been used in various works to overcome singularities of ϕ i>0 singularities at the origin, see, e.g., [START_REF] Buvoli | A class of exponential integrators based on spectral deferred correction[END_REF]. It has been used mainly for parabolic problems [START_REF] Sheen | A parallel method for time-discretization of parabolic problems based on contour integral representation and quadrature[END_REF], also pointing out the potential of parallelization for the first time, as well as using a Carathéodory-Fejér method [START_REF] Schmelzer | Evaluating matrix functions for exponential integrators via carathéodory-fejér approximation and contour integrals[END_REF]. Regarding real applications, it was applied to nonlinear shallow-water equations on the rotating sphere [START_REF] Schreiber | Exponential integrators with parallel-in-time rat. approx. for the shallow-water equations on the rotating sphere[END_REF], providing improved wallclock time-to-solution by using an enlarged and shifted contour to avoid numerical cancellation errors.

Parallel-in-time.

Overcoming the wallclock time limitations of simulations, which cannot be accomplished by any further increase of parallelization in the spatial dimension, is the main focus of the parallel-in-time algorithms. Here, two different types of approaches exist: (a) minimally-invasive methods that take existing time integration methods and incorporate them into an iterative-in-time correction scheme (see, e.g., Parareal [START_REF] Lions | Résolution d'EDP par un schéma en temps "pararéel[END_REF] and PFASST [START_REF] Minion | Interweaving PFASST and Parallel Multigrid[END_REF]); and (b) invasive methods that replace an existing time stepping with one that works entirely differently. Very often, one likes to use a combination of these approaches to enhance the convergence speed of the correction scheme in time. REXI is an invasive parallel-in-time algorithm (see [START_REF] Schreiber | Beyond spatial scalability limitations with a massively parallel method for linear oscillatory problems[END_REF]) since it requires efficient complex-valued solvers for each REXI term.

Unified REXI formulation.

We start directly with the REXI formulation which will provide a standard fundament for the different variants to infer REXI coefficients. Given a discrete linear operator L, we can use an eigendecomposition This manuscript is for review purposes only. L = QΛQ -1 with the eigenvectors stored in the columns of Q and the eigenvalues placed correspondingly on the diagonal of Λ.

(3.1) ∂U (t) ∂t = LU (t) = QΛQ -1 U (t)
where Q and Λ are the matrices with the eigenvectors and eigenvalues on the diagonal, respectively. In terms of the characteristic variable u = Q -1 U and due to diagonalonly Λ, we get independent equations of the form

∂u i (t) ∂t = λ i u i (t)
with λ i the individual Eigenvalues on the diagonal of Λ. In characteristic variables, the unified REXI formulation (1.3) becomes

(3.2) u i (t + ∆t) ≈ γu i (t) + N n=1 β n (∆tλ i -Iα n ) -1 u i (t).
Since each component u i is decoupled, we can freely drop the subscript. For the purpose of time integration, the linear operator L is completely described by its eigenvalues λ, where imaginary components (λ) represent oscillation and negative real values (λ) < 0 describe a diffusive/damping behavior. Note that substituting

λ = 1, t = 0, ∆t = x, u(0) = 1 in (3.2) yields exp(x) = γ + n β n (x -α n ) -1 , which
provides intuition as a sum of rational functions.

Exploiting symmetry of coefficients.

We note that it is possible to reduce the workload by a factor of two for real-valued operators L when the poles α consist of complex conjugate pairs (see, e.g., [START_REF] Lynch | Filtering integration schemes based on the Laplace and Z transforms[END_REF][START_REF] Haut | A high-order time-parallel scheme for solving wave propagation problems via the direct construction of an approximate time-evolution operator[END_REF]). This optimization does not change the relative performance of the methods we consider here, so for simplicity, we do not apply it.

3.2. REXI-derived higher-order ϕ forms. Particular higher-order exponential time integrators such as (2.3) require evaluations of higher-order ϕ i|i>0 . REXI coefficients for these functions are so far computed with methods tailored to them, see [START_REF] Haut | A high-order time-parallel scheme for solving wave propagation problems via the direct construction of an approximate time-evolution operator[END_REF][START_REF] Schreiber | A parallel time integrator for solving the linearized shallow water equations on the rotating sphere[END_REF]. We briefly present an new alternative way to compute them which is easily applicable. Given REXI coefficients for

ϕ i (x) ≈ γ + n β n (x -α n ) -1
we can compute higher-order REXI approximations with

ϕ i+1 (x) = ϕ i (x) -ϕ i (0) x = γ + n βn x-αn -ϕ i (0) x = n β n α n (x -α n ) + 1 x      n - β n α n + γ -ϕ i (0) =0      = n βn αn
x -α n .

The cancellation of the terms is a consequence of the stationary modes which require n -βn αn + γ = ϕ i (0). Note that this leads to different coefficients compared to tailored computations.

This manuscript is for review purposes only. known that no off-the-shelf solvers such as GMRES and multigrid methods work in a highly-scalable way (see, e.g., [START_REF] Ernst | Why it is Difficult to Solve Helmholtz Problems with Classical Iterative Methods[END_REF]). This is ongoing research, and in the present work, we are using solvers developed in spherical harmonics formulations, hence a solver tailored particularly to one PDE problem.

4. REXI methods. These sections cover various ways to infer REXI coefficients, which represent, from our point of view, the most interesting cases. The goal is not to show all methods in great detail but their fundamental properties.

Although we present methods in characteristic form (3.2), the proposed methods also hold in system form (1.3). In the following, we will use the error

(4.1) e(z) = γ + n β n (z -α n ) -1 -exp(z)
to compute the deviation from ϕ 0 (z) = exp(z) with z = λ∆t denoting the point on the complex plane to evaluate. Since approximating diffusive problems is relatively straightforward, we focus on purely oscillatory problems with λ ∈ i R. The REXI methods we consider have complex-conjugate poles α, thus e(z) = e(z) and so we only plot errors for [z] ≥ 0.

4.1. B-REXI: Butcher/Bickart. A Butcher table [START_REF] Butcher | Implicit Runge-Kutta Processes[END_REF] provides a canonical representation of s-stage Runge-Kutta methods [START_REF] Runge | Über die numerische Auflösung von Differentialgleichungen[END_REF][START_REF] Kutta | Beitrag zur näherungsweisen integration totaler Differentialgleichungen[END_REF] in terms of a matrix A ∈ R s×s and completion vector b ∈ R s , with c = A1 determining the abscissa (which we will see are related to REXI poles and 1 is a column vector of ones). The coefficients are selected to achieve the desired order of accuracy and stability properties as well as solution procedure, such as explicit, diagonally implicit, and fully implicit.

For fully nonlinear and non-autonomous ODEs ∂u ∂t = f (t, u), a Runge-Kutta method in Butcher form requires solving a system of stage equations (4.2)

y s = u n + ∆t S j=1
A sj f (t + c j ∆t, y j ), i = 1, ..., S and evaluating the completion formula u n+1 = u n + ∆t S j=1 b j f (t + c j ∆t, y j ). Here, ∆t is the time step size, and y = {y j } S j=1 is the vector of stage solutions. For linear autonomous equations, we can choose characteristic variables, in which case f (t, u) = λu, and the stage equations (4.2) reduce to y = 1u + ∆tλAy and (4.3)

u n+1 = 1 + ∆tλb T (I -∆tλA) -1 1 R(∆tλ) u n ,
where we have identified the stability function R(z) ≈ exp(z).

Derivation.

We now show that unified REXI is algebraically equivalent to Runge-Kutta methods with a diagonal Butcher matrix A, starting with a decomposition inspired by the solution method developed independently by [START_REF] Butcher | On the implementation of implicit Runge-Kutta methods[END_REF][START_REF] Bickart | An Efficient Solution Process for Implicit Runge-Kutta Methods[END_REF]. Given an eigendecomposition A = EDE -1 (which exists for the collocation methods we will consider [START_REF] Hairer | Solving ordinary differential equations I: Nonstiff problems[END_REF]), we can rewrite (4.3) as (4.4)

u n+1 = 1 + ∆tλb T E(I -∆tλD) -1 E -1 1 u n .
This manuscript is for review purposes only. With W = diag(E -1 1) -1 , we may transform to (4.5)

u n+1 = 1 + ∆tλ b T EW -1 bT (I -∆tλD) -1 W E -1 1 1 u n ,
which is a diagonal Runge-Kutta method with A replaced by D and the original completion vector b replaced by b. Rewriting this to a REXI form leads to

u n+1 = u n + ∆t bT - (∆tD)
-1

I + (∆tλD -I) -1 1u n = 1 - bT D -1 1 γ u n + - bT D -2 β T   ∆tλ -D -1 diag(α)   -1 1u n . (4.6)
Finally, we can write this in the unified REXI formulation (1.3) with

γ = 1 - bT D -1 1 β T = - bT D -2 α = diag D -1 . (4.7)
We have derived a transformation from implicit RK method with nonzero eigenvalues to REXI form with the same stability function. Given a REXI method, one can construct an equivalent diagonal RK method (with complex coefficients) via

D = diag(α) -1 and bT = -β T D 2 . Note that a conventional Butcher table A, b T
is not uniquely determined by this procedure. We remark that standard techniques for analyzing Runge-Kutta methods can readily be applied to REXI methods. This includes barriers such as Theorem 4.3 of [START_REF] Jackson | The potential for parallelism in RK methods. Part 1: RK formulas in standard form[END_REF], which establishes that diagonal (parallel)

RK methods can be no more than second order accurate for nonlinear problems.

Error studies.

We choose the Gauss-Legendre and Chebyshev quadrature points for the error studies, with results given in Figure 4.1. We can observe that increasing the number of stages in the non-diagonalized version (using a dense Butcher table) always improves accuracy per stage. In contrast, B-REXI accuracy degrades when too many stages are used, becoming apparent beyond 8 stages. This effect is related to ill-conditioning that can be interpreted via the condition number of the eigenbasis E that effects diagonalization (4.5) or via the 1-norm of the completion vector b, as shown in Figure 4.2. Note that completion vectors must sum to 1 so b 1 = 1 is optimal (and indeed holds for the original completion vector b);

This manuscript is for review purposes only.

. Condition number of the eigenbasis E for the B-REXI method on Gauss-Legendre collocation points (left) and 1-norm of completion vector b for the diagonalized method (right). The rounding errors incurred by the exponential growth precludes use of this approach for many stages.

a large 1-norm indicates the existence of large positive and negative entries, leading to cancellation errors. Despite this downside, the numerical experiments of §6 will show that these B-REXI (diagonalized Gauss Runge-Kutta) methods with lower stage counts are remarkably efficient compared to the other (better-conditioned) families.

Relation to Crank-Nicolson

. Since this will be relevant for the results section, we would like to show the relation between the B-REXI approximation with a single pole using the Gauss-Legendre quadrature using just a simple pole (centered at the interval). This will lead to the terms γ = -1, α = 2 and β = -4 which yields the REXI approximation exp(x) ≈ -1

+ -4 x-2 = 1+ 1 2 x 1- 1 2 x
with the equation on the right hand side matching the Crank-Nicolson formulation. This REXI approximation with a single term resembles the Crank-Nicolson formulation with a midpoint rule (forward Euler on nominator and backward Euler on denominator for x = ∆tL and a half-time step size). This will explain that later numerical results with B-REXI match the Crank-Nicolson method. Using more REXI poles will result in even higher-order approximations.

T-REXI: Terry's Rational Approximation of the Exponential In-

tegrator. The approach which we will refer to as T-REXI was introduced in [START_REF] Haut | A high-order time-parallel scheme for solving wave propagation problems via the direct construction of an approximate time-evolution operator[END_REF].

Several steps are required to gain the α and β coefficients. Since these steps account for the computational workload and the properties, we briefly describe the derivation, including a discussion on the advantages and limitations of this method.

Derivation. The first step consists of an approximation of a Gaussian

basis function ψ c (x) as follows:

(4.8) ψ h (x) = (4π) -1 2 e -x 2 /(4h 2 ) ≈ Re W k=-W ω k i x h + (µ + ik)
Using W = 11, hence L = 2W + 1 = 23 terms in total, is sufficient for an accurate approximation up to numerical double precision (see [START_REF] Haut | A high-order time-parallel scheme for solving wave propagation problems via the direct construction of an approximate time-evolution operator[END_REF]). The advantage of this representation is an efficient representation of the Gaussian basis function in Fourier space. The proxy with the Gaussian basis function allows for computing the coefficients ν k for an approximation of an oscillatory function within an approximate range This manuscript is for review purposes only.

Absolute error for u( t)

T-REXI method, M=4, N=62

Fig. 4.3. Error studies for the T-REXI method for different M and h values. M relates to the number of Gaussian functions via 2M +1 to approximate an oscillation with N = 2(2M +L) number of REXI terms. Left image: We can observe a very high error for a low number of Gaussian basis bumps, which cannot be improved by changing h. Right image:

Using significantly more Gaussian functions leads to significant improvements. In particular, we observe that optimal values for h influence the quality of the approximation. An optimum can be observed for h ≈ 0.8.

x ∈ [-M h; M h] in Fourier space, yielding (4.9)

exp(ix) ≈ M k=-M ν k ψ h (x + kh).
Both steps are then combined, resulting in the approximation

Re (exp(ix)) ≈ M +W n=-M -W Re β Re n (ix + α n ) -1 (4.10)
where we only showed the Re one. We combine Re and Im to the form

exp(ix) ≈ M +W n=-M -W β n (ix -α n ) -1 (4.11)
eventually leading to the REXI formulation with γ = 0, but x related directly to the imaginary value on the complex plane. So far, we only targeted the ϕ 0 function, and we like to point out that this method can also be used to approximate other ϕ i terms (see [START_REF] Haut | A high-order time-parallel scheme for solving wave propagation problems via the direct construction of an approximate time-evolution operator[END_REF]) or directly with the REXI coefficients (see §3.2). We want to emphasize that this approximation was derived only for purely oscillatory functions and, hence, does not include approximations with non-zero real eigenvalues components. that increasing M leads to a linear increase of the size of the region of high accuracy (see [START_REF] Haut | A high-order time-parallel scheme for solving wave propagation problems via the direct construction of an approximate time-evolution operator[END_REF]) with an optimum value of h ≈ 0.8. For the remainder of this work, we will use h ≈ 1.0 as a compromise between accuracy and total workload.

Error studies. We investigate the errors of the T-REXI method in

CI-REXI: Cauchy Contour Integral method.

Cauchy Contour Integral (CI) methods offer an alternative way to infer the REXI coefficients (see e.g. [START_REF] Trefethen | Talbot quadratures and rational approximations[END_REF][START_REF] Buvoli | A Class of Exponential Integrators Based on Spectral Deferred Correction[END_REF][START_REF] Schreiber | Exponential integrators with parallel-in-time rat. approx. for the shallow-water equations on the rotating sphere[END_REF]). We start with the general CI equation given by

g(x) = 1 2πi Γ g(z) z -x dz (4.12)
This manuscript is for review purposes only.

where g(x) is one of the analytic ϕ i functions (2.4), Γ the contour enclosing the eigenvalue λ for ODEs and all eigenvalues on the diagonal of Λ for PDEs.

g(x) = 1 2πi 1 0 g(σ(w)) σ(w) -x σ (w)dw = 1 0 i(2π) -1 g(σ(w))σ (w) x -σ(w) dw. (4.13)
Using the trapezoidal rule, which is exponentially fast converging on periodic boundaries (see [START_REF] Trefethen | The Exponentially Convergent Trapezoidal Rule[END_REF]) with N trapezoidal points in total, we obtain

g(x) ≈ 1 N N n=1 i(2π) -1 g(σ(w n ))σ (w n ) x -σ(w n ) with w n = n N . (4.14)
Again, we can infer a unified REXI formulation (1.3) by setting (4.15)

α n = σ(w n ) β n = ig(σ(w n ))σ (w n ) N 2π γ = 0.
An ellipse contour is given by σ(w) = R x cos(iw2π) + iR y sin(iw2π) -µ with µ related to the center of the ellipse. This leads to the coefficients [START_REF] Hairer | Solving ordinary differential equations I: Nonstiff problems[END_REF])

α n = R x cos(iw2π) + iR y sin(iw2π) -µ (4.
β n = i N exp(σ(w)) (-R x sin(iw2π) + iR y cos(iw2π)) (4.17)
A study of all kinds of contour shapes (rectangle, bean, polygonal shapes, etc.) is beyond the scope of this work. In the next section, we will mainly focus on the circle to show interesting characterizations and use the ellipse for numerical studies to show its superiority to another REXI method. In the following, we will refer to the special case of a circle as CI-REXI and to the ellipse case as CI-EL-REXI.

Characterization and numerical issues.

Next, we characterize the REXI terms, referred to as the β characterization, with an overview in Figure 4.4.

We remind the reader that REXI approximates functions with a linear combination of rational basis functions. Depending on the placement of these functions (related to α n ) and the weighting of each basis (related to β n ), we have three different cases: very large positive and negative numbers, resulting in severe cancellation errors in this region. Consequently, this region should be avoided.

Examples of different contours are provided in Figure 4.5. Each contour targets a particular problem. The circle can be used for the approximation of a small spectral radius λ∆t < 10. Once requiring a larger approximation along the imaginary axis, the radius cannot be enlarged without sacrificing accuracy due to cancellation errors in β n , see (c) above. This can be avoided by enlarging the radius and choosing the value µ, hence shifting the circle, to exclude a contour across areas with Re(x) > 10, which leads to the shifted circle. Other contours are, e.g., given by the ellipse or rectangle targeting the approximation of a spectrum on or close to the imaginary axis and the bean contour targeting an approximation of diffusive and oscillatory problems.

Studies about these contours are beyond the scope of this work and we will focus on the (shifted) circle and ellipse throughout the remainder of this paper.

Error studies.

We conduct error studies using the shifted circle CI-REXI method. The first study is based on a circle centered at the origin, with studies for different radii. The second is for a circle which is shifted to overcome problems related to the cancellation effects (see (c) above). This manuscript is for review purposes only. The centered circle suffers from cancellation effects for large radii, whereas the shifted circle limits these effects. In particular, for a larger imaginary spectrum to be approximated, adding more REXI poles leads to improved accuracy, which is not the case for the centered circle CI-REXI method. REXI poles. We can observe that the errors significantly increase for the centered circle once the radius exceeds a certain threshold. In particular, errors for a larger radius -including a larger spectrum on the imaginary axis -are outside the plotting range. The results for using a higher number of REXI poles do not significantly improve the results. Using a shifted and enlarged circle, we can gain improved results that overcome cancellation errors.

Results are given in

So far, we only investigated the error itself but neglected the total workload. Pruning β with β (exploiting characterization (a)), we can reduce some workload significantly as depicted in Figure 4.7 for larger radii. For a moderate number of REXI poles N = 128 (left image), we observe a pruning close to the accuracy of REXI itself, hardly impacting the results. In contrast, larger radii already suffer from inaccuracies of the used quadrature, with errors outside the plotting range. For a larger number of REXI poles N = 2048 (right image), we observe very robust pruning, hardly affecting the accuracy of the REXI approximation quality but leading to a significant reduction of the workload.

Stability, normalization & filtering.

So far, we have only studied errors in approximating the ϕ 0 function with REXI methods. However, once we use REXI methods for time integrating differential equations, additional properties such as stability and convergence are assumed to be relevant. We will investigate these properties in this section for the ODE du(t) dt = λu(t).

This manuscript is for review purposes only. Im( ) on the imaginary axis once we reach the boundaries of the approximation range. This 396 can be avoided by an additional T-REXI filter, which could be applied to obtain 397 stability also outside the approximation range (see [START_REF] Haut | A high-order time-parallel scheme for solving wave propagation problems via the direct construction of an approximate time-evolution operator[END_REF]).

398

CI-REXI: Finally, we look at the CI-REXI method based on Cauchy contour 399 integral methods in Figure 5.2. The left image shows an unstable region along the 400 imaginary axis. This is caused by an α pole directly placed on the imaginary axis.

401

We can avoid this by choosing the support points of the trapezoidal rule differently.

402

The right image shows a solution to this by shifting them by a half interval, effectively 403 avoiding this instability, and CI-REXI becomes unconditionally stable for oscillatory 404 systems. To summarize, if using the CI-REXI method, one should avoid placing poles 2 6 2 5 2 4 2 3 2 2 2 1 0 2 1 2 2 2 3 2 4 2 5 2 6

Re( )

2 6 2 5 2 4 2 3 2 2 2 1 0 2 1 2 2 2 3 2 4 2 5 2 6
Im( ) near the eigenvalues of the linear operator.

A contour plot comparing CI methods with an increasing number of poles and approximation range is provided in Figure 5.3. We adopted the contour of the circle to pass through the points ±10i, ±20i, and ±40i on the imaginary axis for increasing the number of poles while keeping the contour never exceeding 10 on the real axis.

Normalization. This section concerns particular problems for stationary

or nearby modes requiring special treatment with the T-REXI method. So far, we only assessed errors for a single time step, and this section will investigate the accuracy and conservation properties of stationary modes concerning REXI methods. We will use Dahlquist's equation ( 5) with λ = 10 -3 i, which is time-integrated until t = 100 using different REXI methods. The particular choice of this low frequency is related to almost stationary modes of PDEs. Such modes play an important role, e.g., for geostrophic balance in atmospheric simulations, and not preserving them might lead to spurious/parasitic modes.

An investigation of the results at the absolute ODE errors at t = 100 is given in The results for this are given in the right column, where we observe relatively small improvements for the CI-REXI method (right top image). However, for the T-REXI method, the errors significantly drop from 10 -8 to about 10 -13 once applying this normalization. We also do not see any accuracy degradation for very large time step sizes. Hence, this normalization can be used without impacting the accuracy of other choices of λ, and we will use it throughout the remainder of this work.

We close with two side comments: First, studies for purely stationary modes (λ = 0, not shown here) showed that the errors are increasing for smaller time step sizes, but only due to round-off errors. Overall, these results still lie within the range of numerical precision; hence, we skipped them here. Second, we skipped the B-REXI method since it is not prone to this problem for a number of REXI terms usable as solvers.

This manuscript is for review purposes only.

. Error studies for different REXI methods of time step size ∆t vs. absolute error at u(t = 100). The left column shows errors without normalization, and the right column shows errors with normalization for near-stationary modes. As we can observe in the first row, the normalization for near-stationary modes with the CI-REXI method does not lead to any significant improvements.

In contrast, significant improvements can be observed for the T-REXI method. See the text for a detailed explanation of the results.

Filtering.

This brief section points out the filtering capabilities of the different REXI methods. We define a filter to apply a reduction of the amplitude of ϕ i (x) for a particular set of eigenvalues. It is, in particular, desirable to filter out (setting them close to zero) the so-called "fast modes" for x starting at a threshold and to have a smooth transition of the change in amplitude towards filtering out modes. Since diffusive problems already have a reduction of amplitude given naturally by their mathematical properties, we will again solely focus on purely oscillatory problems, with results also applicable to a mix of oscillatory/diffusive problems.

Using the B-REXI method, we can observe that the stability contour follows exactly the imaginary axis. Hence, there is no filtering at all. For the T-REXI method, we skip a discussion of filtering due to the inherent instability at the boundaries of the approximation range and point out to an additional filter proposed in [START_REF] Haut | A high-order time-parallel scheme for solving wave propagation problems via the direct construction of an approximate time-evolution operator[END_REF]. The CI-REXI method has a natural filtering. This is due to the property that points outside the contour are rapidly approaching 0 as a property of the Cauchy contour integral.

Comparison of REXI methods.

This section aims to provide guidance about which REXI method is best, and we will explore this in different ways. A This manuscript is for review purposes only.

. Error studies for different REXI methods of time step size ∆t vs. absolute error at u(t = 100). B-REXI (left top) is suitable only for smaller timestep sizes. CI-REXI can be tuned to allow also very large time step sizes. T-REXI requires many poles for small time step sizes and allows also very large time step sizes. CI-EL-REXI allows also very large time step sizes and in addition requires the least number of poles for similar accuracy.

full exploration of all parameter combinations is obviously not possible. Hence, we 456 focused on the ones that were most rational to us based on far more experiments than 457 shown here. We first continue with concrete examples using a linear oscillatory ODE 458 based on the Dahlquist equation followed by a PDE with the nonlinear shallow-water 459 equations on the rotating sphere to gain insight into numerical properties once we 460 apply this to more realistic test cases.

461

Based on the eigendecomposition, we classify linear operators as oscillatory, dif-462 fusive, or both. A purely oscillatory system [START_REF] Cox | Exponential time differencing for stiff systems[END_REF] requires imaginary-only eigenvalues 463 (iλ ∈ R) whereas a diffusive behavior is based on negative real eigenvalues (λ ∈ R and 464 λ < 0). Since oscillatory/hyperbolic systems belong to the most challenging problems 465 for REXI methods, we will solely focus on them.

473

The CI-REXI method (right top image) is tuned with a contour never exceeding a real value of 10 and to include the points on the imaginary axis given by Im(λ max ).

The CI-REXI method clearly outperforms the B-REXI method for medium-sized time step sizes and also allows taking very large time step sizes.

The T-REXI method (left bottom image) requires a significant number of REXI poles if only small time step sizes should be taken. This improves once larger step sizes are taken since the initial overheads of the large number of poles (due to the rational approximation of the Gaussian) have less relative impact on the number of REXI coefficients.

In addition, we also investigated the CI-EL-REXI (right bottom image) method, which is a natural choice for purely oscillatory problems. We chose the semi-major axis of the ellipse along the real axis in an empirical way and never exceeding 10 to avoid numerical issues. This method outperforms both CI-REXI and T-REXI almost everywhere regarding accuracy and number of terms required to solve it.

PDE example.

In this final section, we will investigate different REXI methods with the shallow-water equations (SWE) on the rotating sphere. We decided not to investigate many different PDEs, but to go into depth of exponential integration for a single one which is of purely hyperbolic nature. We chose the SWE since they are frequently used to assess the quality and performance of discretizations in time and space concerning horizontal aspects of the full Euler equations solving the fluid dynamics equations related to the atmosphere. In velocity form, the nonlinear SWE are given by

∂ ∂t Φ V = -Φ∇ • V -∇Φ Lg: linear gravity + 0 -f k × V Lc: linear Coriolis + -∇ • Φ V -V • ∇ V N : nonlinear term (6.1)
with the horizontal velocity V on the longitude/latitude field, geopotential Φ = g • h with height h, average geopotential Φ = g • h with average height h, a linearization around a state h = 10 5 m, Coriolis effect f = 2Ω sin(φ) with latitude φ and angular rate of rotation Ω. We like to emphasize that no (hyper)viscosity is used in this PDE to avoid a simplification of the problem due to diffusive effects.

We use this PDE due to its particularly interesting features: The linear gravity term L g is the stiffest one and can be solved with exponential integrators either analytically or with REXI. We want to point out that a comparison of some methods has already been under investigation in former work [START_REF] Schreiber | Exponential integrators with parallel-in-time rat. approx. for the shallow-water equations on the rotating sphere[END_REF] but solely with the CI-REXI method and the geopotential field, which has also been identified to be the best Strang-split method. Anyhow, this study also lacked comparisons with other variables, particularly other REXI methods, which will lead to new revelations, as presented in the following sections. Since including the T-REXI method would not provide any beneficial insight, since the CI-REXI method is computationally much cheaper and provides additional benefits, we skip this method in the following studies.

Spatial discretization.

We solve these equations using the SWEET software 1 which utilizes spherical harmonics (SH) to solve these equations. Such a global spectral basis leads to a substantial reduction of spatial errors (besides a lack of nonlinear interactions at the limit of resolution), hence allowing us to put the focus on time integration methods. We like to refer to [START_REF] Schreiber | A parallel time integrator for solving the linearized shallow water equations on the rotating sphere[END_REF][START_REF] Hack | Description of a global shallow water model based on the spectral transform method[END_REF] for a detailed description of the spherical harmonics. In particular, we work with the vorticity-divergence formulation in spectral space to avoid spurious modes if one would convert the velocity to spectral space. The standard 2 3 rule [START_REF] Orszag | On the Elimination of Aliasing in Finite-Difference Schemes by Filtering High-Wavenumber Components[END_REF] is used for anti-aliasing to evaluate bi-non-linearities.

Time stepping solvers.

We can find highly efficient solvers in spherical harmonics space for direct exponential integration (without numerical approximations), REXI, and implicit Euler time integrators.

Regarding the direct exponential integration, we can straightforwardly find a direct solution using the vorticity divergence form, see also [START_REF] Thuburn | Numerical representation of geostrophic modes on arbitrarily structured C-grids[END_REF]. Due to orthogonality, each mode can be separately written as

∂ t Φ ∂ t δ = -Φ -∇ 2 Φ δ = G D Φ δ
with the famous identity ∇ 2 = -n(n + 1) for this harmonic. Using D = -∇ 2 and G = -Φ for convenience, we find the eigenvectors Q and eigenvalues diag(Λ)

Q = -G D + G D 1 1
Q -1 =   1 2 D G 1 2 -1 2 D G 1 2   Λ = - √ DG √ DG .
We can then use U (t + ∆t) = Q exp(∆tΛ)Q -1 . From an algebraic perspective, this method matches the method in [START_REF] Harney | Laplace transform integration of a baroclinic model[END_REF], which uses a rather cumbersome derivation using Laplace transforms, whereas our derivation is more elegant and short. This method is also used for investigating errors

For the exponential integration of the full linear terms L = L g + L c , this would relate to the Hough modes [START_REF] Wang | On computation of Hough functions[END_REF] and no direct exponential solution has been derived yet. Hence, it requires evaluations of the form (∆tL -α) -1 with complex-valued α.

The first time this was solved for REXI using spherical harmonics was based on a method requiring transformations to grid space [START_REF] Schreiber | A parallel time integrator for solving the linearized shallow water equations on the rotating sphere[END_REF]. The present work is based on an implicit time stepper [START_REF] Temperton | Treatment of the Coriolis Terms in Semi-Lagrangian Spectral Models[END_REF] of the form (I -∆tL) -1 which has been transformed to solve a REXI term by simply using a complex-valued time step size. We also used it for the implicit time integration of L as it has been originally suggested.

Benchmark.

Our benchmark is based on the barotropic instability test case (see [START_REF] Galewsky | An initial-value problem for testing numerical models of the global shallow-water equations[END_REF]). This benchmark is initialized with a geostropically balanced initial

This manuscript is for review purposes only.

Short notation Description

ERK(X, o = N ) Explicit Runge-Kutta with order N IRK(X) Backward Euler using 2nd order Crank-Nicolson SS(X, Y ) 2nd order Strang-splitting as explained in the text with

F 1 = X and F 2 = Y EXP (X)
Direct exponential integration on X REXI(X)

A particular REXI method on X ET DRK(X, Y ) 2nd order ETDRK method with X being exponentially integrated and Y treated as the nonlinearity X + Y Time tendencies of terms X and Y are added Table 6.1 Overview of time integration methods. Note that they can be composed together. condition, which is perturbed by a small Gaussian bump (see reference for detailed initial conditions). We time integrate this system for 8 days with results in Figure 6.2.

Time integration.

Regarding the particular Runge-Kutta (RK) based time integrators, we used 2nd order midpoint, 3rd order Heun, and classical 4th order RK. The reference solution to compute the errors is based on the 4th order RK with a time step size of ∆t = 5.

Besides the methods already introduced, our investigation also includes the 2nd order Strang splitting (SS) method [START_REF] Strang | On the Construction and Comparison of Difference Schemes[END_REF]. With SS, a PDE given by two terms d dt U = F 1 (U ) + F 2 (U ) can be integrated with 2nd order accuracy if a 2nd order accurate time integrator R ∆t Fi is provided for time step size ∆t by

U (t + ∆t) = R 1 2 ∆t F1 • R ∆t F2 • R 1 2 ∆t F1 .
We use a function-like notation to refer to the particular time integration methods.

An overview of this is given in Table 6.1 where we use X and Y as representatives for either term in the PDE such as L g , L c , and N or to refer to another time integrator.

In the latter case, e.g., ERK, EXP , REXI, and IRK can both be used in the Strang-Splitting SS as arguments. by utilizing the 2nd NUMA domain, then other compute nodes. We gratefully acknowledge the usage of the SHTNS library [START_REF] Schaeffer | Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations[END_REF] which is based on FFTW [START_REF] Frigo | A Fast Fourier Transform Compiler[END_REF]. We precomputed transformation plans and reused them for all studies to ensure the utilization of the same ones over all studies. Each batch job is set to timeout after 1 hour, which follows the idea that the simulations should be finished within a specific time frame.

6.2.6. Performance comparison for splitting L g and L c + N . We start with a comparison of standard methods in Figure 6.3 which we will use as a baseline for further comparisons with REXI-based methods. Plots are given for the three prognostic variables, which we define here as the variables required as input to one time step, since results differ for all of them.

First, the higher-order 3rd-and 4th-order RK method can outperform other lowerorder methods for smaller time step sizes depending on the variable under study. This

This manuscript is for review purposes only. Fig. 6.3. Studies without REXI methods (but using direct exponentiation) on all prognostic variables with error vs. time step size for the barotropic instability benchmark. We also include 2nd, 3rd, and 4th order Runge-Kutta based methods with gray lines. is a known phenomenon for higher-order time integration methods, and we wanted to include it to also see its max. stable time step size. We are primarily interested in very large time step sizes while still having a moderately small error.

The best method concerning the geopotential and the divergence variable is the Strang-split SS(ERK(L c +N ), EXP (L g )), which we account for by the more accurate treatment with the exponential treatment of both variables. Since the vorticity field is not treated exponentially (time tendency for this in L g is null), there's also no benefit visible in the comparison of the vorticity field.

The ETDRK method itself -although assumed to be an excellent off-the-shelf method -does not provide the overall best results compared to the rather straightforward Strang splitting. We can observe it to be the 2nd best for the geopotential and even lower ranked for the other variables. We account for that by the way a 2nd order accurate Strang-splitting is performed. This can be interpreted as a subcycling of time steps by executing two half-time steps for one of the terms (the time step size limiting one).

Next, we will continue with REXI studies by comparing them with the best Strang-split exponential and implicit methods from the previous results in Figure 6.4.

Overall, we can observe a 2nd order convergence even if using only a single pole for the B-REXI methods.

Matching results for SS(REXI,ERK) B-REXI N=1 and SS(IRK,ERK) are observed which is explained in §4.1.3: This particular B-REXI method resembles exactly the Crank-Nicolson method but uses one complex-valued pole to solve the system of equations.

This manuscript is for review purposes only. The other B-REXI methods outperform all alternatives except for the direct exponential integration EXP. We see one particularly interesting and highly important effect: The B-REXI method does not provide any further advantages using more than N = 2 poles. Even using N = 4 poles, the results are not further improved. A particularly important point is the comparison of the CI-REXI method with B-REXI, where absolutely no benefits are visible for N = 128 poles using CI-REXI compared to N = 2 poles using B-REXI. This clearly indicates that significant computational savings of a factor of 64 can be accomplished in this case compared to the former work.

We close this section by HPC studies in Figure 6.5. For sake of better overview, we only plotted the most promising candidates (ETDRK is worse than B-REXI methods, the explicit RK order 3 and 4 methods are better for larger wallclock times (smaller time steps), but unstable otherwise).

We start by comparing the performance of the direct exponential method EXP with the REXI method, where we would expect that the direct method is faster, which is not the case. We account for that by the direct method to be computationally more intensive (square root, exponential, etc., see §6.2.2) in order to solve for this term, whereas the B-REXI methods only require to evaluate two or 4 rational approximations. For the CI-REXI method, which requires N = 128 terms this is again different due to the higher MPI overheads resulting in a lower performance than the others.

Although the Strang-splitting method with the implicit term is computationally quite efficient to evaluate, its overall wallclock time performance is not optimal. 6.2.7. Performance comparison for splitting into L and N . Next, we investigate the performance of REXI methods using a splitting into the linear term This manuscript is for review purposes only. For the geopotential Φ , we can observe significant improvements in terms of accuracy. In particular, we can take very large time step sizes and still observe a convergence, whereas the 2nd order IRK-like methods already stagnate. With respect to the ETDRK scheme, its performance is worse compared to the best (straightforward)

Strang-split methods.

For the vorticity η we can observe that the ETDRK method does not lead to any improvement. The best methods are the Strang-split IRK-based ones and some REXI-based methods. Hence, we do not see any improvement in the accuracy of the vorticity field by using exponential integration methods. This is kind of surprising at first glimpse since we expected a better treatment of the vorticity due to the exponential integration of the Coriolis effect. However, the errors in the nonlinear parts dominate the overall errors. Hence, this does not provide any better results.

The divergence δ study shows REXI methods to be the best ones. Again, the accuracy cannot be improved by using more than N = 2 poles. Everything beyond that would be an additional computational burden. The ETDRK methods again show a poorer performance than the more straightforward approach.

Finally, we investigate the wallclock time vs. errors with results given in Figure 6.7.

We can observe that fully explicit ERK methods actually provide excellent results due

This manuscript is for review purposes only. to their computationally efficient way. In particular, the classical 4th-order accurate ERK method provides excellent results across all prognostic variables.

A closer look at the geopotential Φ errors shows that the B-REXI-based methods with N=2 poles are to be preferred compared to all other methods. Again, the ETDRK method shows no real benefits.

Investigating the vorticity η leads to a different interpretation: Now, the implicit Strang-split method provides the best results which can be easily explained by the situation that the exponential treatment of the L c term did not lead to any beneficial results already in the error vs. time step size plots and additional computational time is required here. Finally, ETDRK are literally the worst in here, not paying off at all. The errors on the divergence δ show similar results compared to the geopotential, which is why we skip a detailed discussion here.

Summary of PDE results.

The CI-REXI method with N = 128 poles is not beneficial at all compared to B-REXI with N = 2 poles. Using only N = 2 poles with the B-REXI method already provides the best results, and no improvement can be gained by adding more poles. This is actually quite surprising, with expectations of exponential integration methods to always provide significantly better results. However, using such a higher-order approximation seems to provide sufficient accuracy so that the errors from the splitting approach dominate the overall errors.

We would like to emphasize that all the statements are specific to the SWE on the rotating sphere PDE and should not be generalized.

This manuscript is for review purposes only. 7. Summary and Conclusions. Exponential integration methods are considered to be a way to integrate with high efficiency. As part of that, ϕ functions need to be solved, which turn out to be computationally rather challenging. This paper investigated different ways to approximate ϕ functions with rational approximations of exponential integration (REXI). The coefficients of REXI methods can be derived in many ways and we introduced a generalized REXI approach, finally allowing to express many different methods in this way. We showed this for the Butcher/Bickard-based REXI, Cauchy Contour integration based REXI and T(erry)-REXI method. All methods have been introduced in a way making its capabilities and limitations easily graspable. With respect to physical properties, the T-REXI method requires special treatment for (quasi-)stationary modes and became obsolete with CI-REXI. In addition, we derived an elegant way to compute higherorder ϕ functions based on REXI coefficients for lower-order ϕ.

An in-depth investigation of the approximation quality of each REXI method has been conducted including an explanation of numerical issues for all of the methods.

Next, we put it into the context of time integration methods. We first used linear ODEs where we studies and discussed properties of stability, convergence and also the filtering capabilities. Second, we performed in-depth studies using the nonlinear shallow-water equations on the rotating sphere. Surprisingly, the best REXI method turned out B-REXI with only N = 2 poles, leading to a significant reduction of computational effort compared to former REXI methods in this context using N = 128 poles. Consequently, regarding demands on computational resources, B-REXI showed a reduction of a factor of 64 compared to previous work. This also means that a This manuscript is for review purposes only.

higher-order implicit Runge-Kutta method is competitive to traditional exponential integration methods for this PDE.

Fig. 4 . 1 .

 41 Fig. 4.1. Error studies for the B-REXI method with (a) Gauss-Legendre and (b) Chebyshev quadrature points for the error given in Eq. (4.1). Each color refers to the same number of stages, marked to B-REXI or RK form. The non-diagonalized version provides significantly better resultscompared to the diagonalized version. In particular, results with B-REXI using 32 or more stages suffer from significant defects in the solution.

1 b 1 :

 11 Fig. 4.2.Condition number of the eigenbasis E for the B-REXI method on Gauss-Legendre collocation points (left) and 1-norm of completion vector b for the diagonalized method (right). The rounding errors incurred by the exponential growth precludes use of this approach for many stages.

  Figure 4.3. On the left image, we can observe that we need a minimum number of Gaussian basis functions to approximate the oscillations. The right image shows exceptionally accurate results for h ≈ 0.8 in the range x ∈ [0; 10] and a rather large region of accuracy of about e(x = 128) ≤ 10 -11 . Other figures (not included) show

4. 3 . 1 .

 31 REXI Derivation. Regarding the contour, we can use different approaches. In what follows, we used parametrized contours Γ = {σ(w)|w ∈ [0; 1]} with the contour function σ(w) : R → C. Using integration by substitution and the contour function, we obtain

  a) Obsolete REXI terms: Contours Re(σ(x)) → -∞ relating to areas of the contour in the distant negative real axis on the complex plane have exponentially fast decaying β coefficients, hence lim σ→-∞ β n = lim σ→-∞ i exp(Re(σ(wn)))σ (wn) N 2π = 0. Once a particular β n coefficient undershoots a threshold β , the corresponding REXI term can be removed if β n < β and β = β /N . The last equation incorporates that a higher numerical resolution results in smaller values of the β weights. b) Regular REXI terms: This characterization refers to those REXI terms that can be incorporated in the approximation in a useful way. c) Cancellation-prone REXI terms: These terms are related to the contour Re(σ(x)) → +∞. Approximating the exp function in the far distance to the right of the origin leads to exponentially increasing the β values. An oscillatory function is also parallel to the imaginary axis, which is approximated. Both effects lead toThis manuscript is for review purposes only.

Fig. 4 . 4 .

 44 Fig. 4.4. Complex plane for the real (left image) and imaginary (right image) value of exp(x).We highlight the different areas related to the different β characterizations.

Fig. 4 . 5 .

 45 Fig. 4.5. Selection of contours used with the CI-REXI method.

  Figure 4.6 with plots based on a fixed number of N = 256

Fig. 4 . 6 .

 46 Fig. 4.6. Error studies for the centered circle (left) and shifted circle (right) CI-REXI method.The centered circle suffers from cancellation effects for large radii, whereas the shifted circle limits these effects. In particular, for a larger imaginary spectrum to be approximated, adding more REXI poles leads to improved accuracy, which is not the case for the centered circle CI-REXI method.

14 10Fig. 4 . 7 .

 1447 Fig. 4.7. Error studies for the shifted circle with N = 128 REXI poles (left) and N = 2048 REXI poles (right) with different β pruning values. We can observe significant reductions in the number of required REXI poles.

Fig. 5 . 1 .

 51 Fig. 5.1. Stability plots. For B-REXI (left): we observe an excellent stability behavior known for collocation methods. T-REXI (right):We observe instabilities at the imaginary axis for the boundaries of the approximation range.

Fig. 5 . 2 .

 52 Fig. 5.2. Stability plots for CI-REXI. The left image depicts the discrete contour points chosen so that one α pole lies on the imaginary axis. This leads to instabilities. The right image depicts a half-shifted variant of it.

393T-

  REXI (right image in Figure 5.1): We can observe that T-REXI provides 394 excellent stability for purely imaginary values. However, we can observe instabilities 395

Fig. 5 . 3 .

 53 Fig. 5.3. Contour lines for stability plots for CI-REXI with half-shifted intervals. We can observe an excellent stability region over this entire range.

Figure 5 . 4 .

 54 Figure 5.4. The left column shows REXI methods as they have been computed with the methods from before. We can observe that the CI-REXI method (top left image) has REXI coefficients preserving the stationary modes. However, the T-REXI suffers from significant defects in it. A normalization can be used to overcome this problem where stationary modes require n βn x-αn = s = 1 and we can ensure this by simply

  Fig. 5.4.Error studies for different REXI methods of time step size ∆t vs. absolute error at u(t = 100). The left column shows errors without normalization, and the right column shows errors with normalization for near-stationary modes. As we can observe in the first row, the normalization for near-stationary modes with the CI-REXI method does not lead to any significant improvements. In contrast, significant improvements can be observed for the T-REXI method. See the text for a detailed explanation of the results.

  Fig. 6.1. Error studies for different REXI methods of time step size ∆t vs. absolute error at u(t = 100). B-REXI (left top) is suitable only for smaller timestep sizes. CI-REXI can be tuned to allow also very large time step sizes. T-REXI requires many poles for small time step sizes and allows also very large time step sizes. CI-EL-REXI allows also very large time step sizes and in addition requires the least number of poles for similar accuracy.

466 6 . 1 .( 5 ) 2 468

 6152 ODE. We investigate the ODE systems again with Dahlquist's equation 467 using λ = 1 and the simulation results at t = 100. We use u(0) = (1 + i)/ √ as an initial condition. We compare various REXI methods in Figure6.1. The total 469 numbers of REXI coefficients are given by N .470The B-REXI method (left upper image) performs extremely well for small step 471 sizes where only a few poles are required. For larger time step sizes of ∆t ≈ 10, using 472 16 poles is sufficient to gain single precision accuracy.

1Fig. 6 . 2 .

 62 Fig. 6.2. Vorticity field of barotropic instability benchmark after 8 days of inviscid shallow-water equations. We see the development of various large and small-scale vortices.

6. 2 . 5 .

 25 Hardware, parallelization & batch configuration. All results have been computed on the Thin Nodes of SUPERMUC-NG. Each node is equipped with two Intel SkylakeXeon Platinum 8174, resulting in two NUMA domains. For the spatial parallelization, we use solely OpenMP on one NUMA domain, resulting in a spatial scaling of up to 24 cores. Scalability for REXI is then based on MPI first

LLLFig. 6 . 4 .

 64 Fig. 6.4. Studies with REXI methods on all prognostic variables with error vs. time step size for the barotropic instability benchmark.

LLLFig. 6 . 5 .

 65 Fig. 6.5. Studies including REXI methods with wallclock time vs. time step size for the barotropic instability benchmark.

LLLFig. 6 . 6 .

 66 Fig. 6.6. Studies using non-REXI methods (using direct exponentiation) on all prognostic variables with error vs. time step size for the barotropic instability benchmark. (ETDRK data points are missing due to 1h timeouts of the job.)

LLLFig. 6 . 7 .

 67 Fig. 6.7. Studies using non-REXI methods (using direct exponentiation) on all prognostic variables with error vs. time step size for the barotropic instability benchmark. (ETDRK data points are missing due to 1h timeouts of the job.)

3.3. Linear solvers for REXI terms. Efficient

	solvers are required for each
	REXI term. Over the last decades, this efficiency aspect turned out to be a very
	challenging task. E.g., in the context of shallow-water equations, this results in the
	original Helmholtz problem (rather than a backward Euler time step) where it is

This manuscript is for review purposes only.

Acknowledgements. Both authors like to thank Pedro S. Peixoto for pointing out the potential relation of exponential integration methods to Laplace transforms and Peter Lynch's work in this context. Martin Schreiber is grateful to NCAR for providing financial support and a very inspiring office space with a splendid view to the flatirons, which strongly supported this work. Both authors thank Matthew Normile for preliminary work as well as Finn Capelle and Raphael Schilling who indirectly contributed to this work with the REXInsight software.

The authors gratefully acknowledge the Gauss Centre for SC e.V. (www.gausscentre.eu) for funding this project by providing computing time on the GCS Supercomputer SUPERMUC-NG at Leibniz Supercomputing Centre (www.lrz.de).

This project has received funding from the Federal Ministry of Education and Research and the European High-Performance Computing Joint Undertaking (JU) under grant agreement No 955701. The JU receives support from the European Union's Horizon 2020 research and innovation programme and Belgium, France, Germany, Switzerland. Martin Schreiber gratefully acknowledges KONWIHR funding as part of the project "Parallel in Time Integration with Rational Approximations targeting Weather and Climate Simulations". Jed Brown acknowledges support from the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, applied mathematics program.