Exponential Runge-Kutta Methods

• Exponential Runge-Kutta methods are time-stepping schemes that blend exact exponential propagation of stiff linear parts with explicit Runge-Kutta evaluations for the nonlinear term.
• Stiff order theory imposes 16 to 36 operator-valued conditions, leading to high-order schemes that require 8 to 16 stages for fifth- and sixth-order accuracy.
• Efficient implementation leverages Krylov subspace methods, scaling-and-squaring, and parallel stage grouping to compute exponential and ϕ-functions in large-scale applications.

Exponential Runge–Kutta (expRK) methods are time-stepping schemes for semilinear initial-value problems that combine Runge–Kutta stage structure with exponential integration of a stiff linear part. In their standard form, they target problems of the type $u'(t)=Au(t)+g(t,u(t))$, where $A$ is typically large, sparse, and stiff, and $g$ is a milder nonlinear perturbation. Their defining feature is that the contribution of $A$ is propagated through $e^{hA}$ and related $\varphi_k$-functions, while the nonlinear term is handled by explicit or implicit RK-like stage evaluations. This construction places expRK methods between classical explicit RK schemes, which suffer from bounded stability regions, and fully implicit RK schemes, which require nonlinear solves at every step (andò et al., 5 Jul 2025).

1. Core formulation

The basic analytical starting point is the variation-of-constants formula for the semilinear problem

$u'(t)=Au(t)+g(t,u(t)), \qquad u(0)=u_0,$

which over one step $h$ reads

$$u(t_0+h)=e^{hA}u(t_0)+\int_0^h e^{(h-\tau)A}g(t_0+\tau,u(t_0+\tau))\,d\tau.$$

Exponential integrators approximate the integral term while keeping the action of the stiff linear part through $e^{hA}$ exact or highly accurate (andò et al., 5 Jul 2025).

Central to the construction are the $A$0-functions,

$A$1

together with the recurrence

$A$2

These functions appear naturally when the semigroup contribution is integrated against polynomial or RK-type interpolants of the nonlinearity (Luan et al., 2013).

A generic $A$3-stage expRK method has stage values $A$4 and update

$A$5

$A$6

$A$7

where the coefficients $A$8 and $A$9 are operator-valued combinations of $g$0-functions (andò et al., 5 Jul 2025). In the simplest case, exponential Euler,

$g$1

the linear homogeneous part is propagated exactly, and if $g$2 the method reduces to explicit Euler (andò et al., 5 Jul 2025).

For parabolic problems, the standard analytical setting assumes that $g$3 is the infinitesimal generator of an analytic semigroup $g$4 on a Banach space and that $g$5 is sufficiently Fréchet differentiable along the exact solution. In that setting, expRK methods are designed to remove severe stability restrictions caused by the linear stiffness while avoiding the nonlinear solves required by implicit RK schemes (Luan et al., 2013).

2. Stiff order theory and high-order construction

The decisive difference between expRK analysis and classical RK analysis is the distinction between classical order and stiff order. In stiff semilinear problems, order conditions must control the interaction between the semigroup $g$6 and the nonlinear remainder rather than just polynomial moments of the stage abscissae. For explicit expRK methods, local error expansions are expressed in terms of operator-valued functions such as $g$7 and $g$8, and the resulting stiff order conditions generalize classical rooted-tree conditions (Luan et al., 2013).

For fifth-order explicit expRK methods for semilinear parabolic problems, a complete set of stiff order conditions up to order $g$9 was derived, comprising 16 conditions in operator form. These include moment conditions such as

$A$0

together with higher-order coupling conditions involving arbitrary matrices $A$1 and bilinear maps $A$2 (Luan et al., 2013). A structural barrier was proved: there does not exist an explicit exponential Runge–Kutta method of order $A$3 with less than or equal to $A$4 stages. The corresponding constructive result was an 8-stage fifth-order method, expRK5s8, based on weakened order-5 conditions tailored to semilinear parabolic problems (Luan et al., 2013).

The sixth-order theory was subsequently reformulated through exponential B-series on an essential tree set $A$5. In that framework, 36 independent stiff order conditions were identified for order $A$6, with conditions 1–16 matching the order-5 theory and conditions 17–36 providing the genuinely new sixth-order constraints. Two families of sixth-order methods were then constructed: $A$7, which weakens condition No. 17 to evaluation at $A$8, and $A$9, which satisfies all 36 conditions in strong form (Luan et al., 2023).

The sixth-order constructions are notable for their stage organization. Although they use 15 and 16 stages, respectively, they are arranged in four groups of parallel stages, and the authors emphasize that they can be implemented at a cost similar to a 6-stage method. Under the analytic-semigroup assumptions and the corresponding order conditions, the resulting global error satisfies

$e^{hA}$0

uniformly on finite time intervals (Luan et al., 2023).

3. Extensions of the analytical framework

The semigroup-based parabolic setting is not the only environment in which expRK methods have been developed. For semilinear integro-differential equations with memory,

$e^{hA}$1

the linear dynamics are governed not by a $e^{hA}$2-semigroup but by a resolvent family $e^{hA}$3, which in general does not satisfy a semigroup property. In that setting, explicit exponential Runge–Kutta methods were reformulated in terms of operator-valued quadrature weights built from $e^{hA}$4 and generalized resolvent-based $e^{hA}$5-functions

$e^{hA}$6

For the linear case, order conditions were derived for general order $e^{hA}$7 and convergence of order $e^{hA}$8 was proved. For the semilinear case, first- and second-order methods were constructed and analyzed together with spectral Galerkin discretization in space (Ostermann et al., 2022).

Delay equations provide another extension. In the sun-star abstract framework, both delay differential equations and renewal equations can be recast as abstract semilinear evolution equations on Banach spaces, allowing explicit expRK methods to be defined through the corresponding semigroup and $e^{hA}$9-operators. This yields a unified treatment of DDEs, renewal equations, and coupled DDE/RE systems, together with convergence results that mirror the classical expRK theory (Ando' et al., 2024). For parabolic problems with arbitrary state-dependent delay, first- and second-order schemes were constructed by combining explicit expRK discretizations with continuous extensions of the time-discrete solution to approximate delayed values. The same work also gives collocation-type schemes of arbitrary order, as well as unique solvability and convergence results (Huang et al., 10 Sep 2025).

A further extension concerns parabolic problems with non-homogeneous boundary conditions. Classical expRK convergence theory is formulated for homogeneous or periodic boundaries, and direct application in the non-homogeneous case typically leads to order reduction. A correction strategy based on smooth extensions $\varphi_k$0 of the boundary data rewrites the original problem as a homogeneous problem for $\varphi_k$1 with modified source term. For linear problems, the corrected formulation recovers the expected order and, with suitable quadrature rules, allows orders up to $\varphi_k$2 for $\varphi_k$3-stage Gauss collocation methods. For semilinear problems, the same boundary correction preserves the convergence orders guaranteed by stiff order conditions (Arranz-Simón et al., 24 Oct 2025).

4. Specialized variants and structure-preserving formulations

Several research directions adapt expRK ideas to specific structural demands. For multiphysics and partitioned systems, stiffly accurate exponential Runge–Kutta methods have been reformulated so that only matrix functions of the Jacobians of individual components or partitions are required. In the additive and component-partitioned settings, the transformed formulations retain the stiff order properties of the underlying expRK scheme while avoiding matrix functions of the full coupled Jacobian. Numerical tests on partitioned Gray–Scott reaction–diffusion problems showed full design order for partitionings by species, by space, by physical process, and in exponential–explicit configurations (Narayanamurthi et al., 2019).

For highly oscillatory systems, two families of cost-reduction implicit exponential Runge–Kutta methods were introduced: the simplified version ERK (SVERK) and the modified version ERK (MVERK). Their defining feature is that the coefficients $\varphi_k$4 and $\varphi_k$5 are real constants rather than matrix-valued $\varphi_k$6-combinations, with the exponential action confined to factors such as $\varphi_k$7 and $\varphi_k$8. The order conditions of these schemes were shown to be identical to the order conditions of classical Runge–Kutta methods. The same work also derived symplecticity conditions and presented a first-order symplectic SVERK method, together with second- and fourth-order constructions and linear stability studies (Hu et al., 2023).

Structure preservation has also motivated expRK developments for the nonlinear Schrödinger equation. By combining the scalar auxiliary variable approach with a Lawson-type exponential transformation and symplectic RK methods, a family of arbitrarily high-order exponential Runge–Kutta schemes was obtained that preserves discrete mass and a modified energy exactly. The Lawson transformation absorbs the linear Schrödinger operator into an exponential, while the SAV reformulation renders the nonlinear energy quadratic in an augmented variable. With symplectic RK coefficients satisfying

$\varphi_k$9

the resulting ESAV-RK methods preserve both discrete mass and the modified energy (Cui et al., 2020).

In kinetic theory, explicit expRK methods were designed for stiff Boltzmann- and BGK-type equations by decomposing the collision operator into an equilibrium and a non-equilibrium part. These schemes are exact for relaxation operators of BGK type, unconditionally asymptotically stable and asymptotic preserving under stated conditions, and they also admit nonnegativity and entropy-inequality properties. The same formulation is suitable for both deterministic and probabilistic numerical techniques (Dimarco et al., 2010).

5. Implementation and computational patterns

In large-scale applications, the dominant cost of expRK methods lies in computing actions of $u'(t)=Au(t)+g(t,u(t)), \qquad u(0)=u_0,$0 and $u'(t)=Au(t)+g(t,u(t)), \qquad u(0)=u_0,$1 on vectors. Several implementation paradigms recur across the literature. Krylov subspace methods approximate $u'(t)=Au(t)+g(t,u(t)), \qquad u(0)=u_0,$2 or $u'(t)=Au(t)+g(t,u(t)), \qquad u(0)=u_0,$3 by projection onto $u'(t)=Au(t)+g(t,u(t)), \qquad u(0)=u_0,$4. Scaling-and-squaring combined with Taylor or Padé approximants is another standard route, while augmented-matrix formulations make it possible to compute linear combinations $u'(t)=Au(t)+g(t,u(t)), \qquad u(0)=u_0,$5 through a single exponential of a larger matrix. Polynomial interpolation on Leja points and contour-integral approaches are also used, especially when robust evaluation of $u'(t)=Au(t)+g(t,u(t)), \qquad u(0)=u_0,$6-functions is required (andò et al., 5 Jul 2025).

Parallelism and coefficient reuse have become increasingly important in high-order designs. The sixth-order methods $u'(t)=Au(t)+g(t,u(t)), \qquad u(0)=u_0,$7 and $u'(t)=Au(t)+g(t,u(t)), \qquad u(0)=u_0,$8 were implemented with phipm_simul_iom, using incomplete orthogonalization Krylov techniques for simultaneous $u'(t)=Au(t)+g(t,u(t)), \qquad u(0)=u_0,$9-function actions, and their stage structure was explicitly organized into four parallel groups (Luan et al., 2023). In partitioned formulations, the savings come instead from replacing matrix functions of the full Jacobian by matrix functions of the Jacobians of individual components or processes (Narayanamurthi et al., 2019).

For resolvent-based memory equations, the principal objects are $h$0 and the corresponding generalized $h$1, which are computed through the spectral decomposition of the elliptic operator and scalar resolvent functions $h$2 satisfying Volterra equations. For Riesz kernels these scalar resolvent functions reduce to Mittag–Leffler functions, which makes exact coefficient evaluation possible in the numerical experiments (Ostermann et al., 2022).

A distinct implementation line appears in the fast nonlinear Fourier transform. There, exponential RK and exponential linear multistep discretizations of the Zakharov–Shabat scattering problem are analyzed through transfer matrices. When the ERK stage nodes are equispaced, the transfer matrices become polynomial in $h$3, which enables FFT-based fast polynomial arithmetic. In that setting, the scattering coefficients can be computed with complexity $h$4 and convergence rate $h$5, where $h$6 is the order of the underlying discretization (Vaibhav, 2018).

6. Limitations, misconceptions, and current directions

A recurrent misconception is that expRK methods solve stiffness in general. The literature is more specific: they are most effective when the dominant stiffness resides in a linear part that can be treated through exponentials or resolvent operators. When stiffness is mainly nonlinear, explicit treatment of the nonlinear term can still impose practical restrictions, and the survey’s Duffing example shows that expRK advantages can become less decisive in that regime (andò et al., 5 Jul 2025).

Another limitation is the rapid growth of the order theory. The number of stiff order conditions rises sharply with order: 16 at order $h$7 and 36 at order $h$8. This growth explains both the stage-count barrier at order $h$9 and the large stage numbers of the first sixth-order constructions (Luan et al., 2013, Luan et al., 2023). High-order methods are therefore available, but their construction is algebraically intricate and often tailored to specific analytical settings.

Order reduction remains a central issue. Non-homogeneous boundary conditions typically cause order reduction unless a boundary correction is introduced through smooth extensions of the boundary data (Arranz-Simón et al., 24 Oct 2025). In sixth-order parabolic experiments, marked order reduction was reported for certain inhomogeneous boundary conditions, while smoother boundary data did not trigger the same effect (Luan et al., 2023). State-dependent delays create a related difficulty: derivative discontinuities can propagate through the delay, and higher-order convergence may require detecting or aligning the time mesh with those discontinuity points (Huang et al., 10 Sep 2025).

Finally, several extensions remain only partially developed. For semilinear integro-differential equations with memory, the linear theory reaches general order $$u(t_0+h)=e^{hA}u(t_0)+\int_0^h e^{(h-\tau)A}g(t_0+\tau,u(t_0+\tau))\,d\tau.$$0, but the detailed semilinear analysis in the cited work is restricted to first and second order (Ostermann et al., 2022). The survey literature also identifies variable stepsizes, non-parabolic operators, nonlinear stiffness, and efficient matrix-function evaluation on modern architectures as continuing research directions (andò et al., 5 Jul 2025). In that sense, expRK methods form a mature but still actively expanding class: their central principle is stable, yet their analytical and computational envelopes continue to widen.