Complex Time Integrators

• Complex Time Integrators are algorithms that use complex-valued time steps to solve ODEs/PDEs, offering higher-order accuracy and expanded stability regions.
• They employ modified Runge–Kutta and Euler schemes with complex coefficients to break traditional order barriers and efficiently handle stiff, dispersive, and quantum systems.
• These integrators enhance computational performance by reducing storage requirements and enabling effective parallelization in large-scale simulations.

Complex time integrators are a class of algorithms for the numerical solution of ordinary and partial differential equations (ODEs/PDEs) that employ time steps with both real and imaginary components. This approach extends the standard real-valued time-stepping paradigm, offering enhanced stability, higher-order accuracy, storage efficiency, and tailored adaptation to stiff, dispersive, and quantum-evolution problems. By leveraging the degrees of freedom available in the complex time plane, complex time integrators can outperform traditional explicit and implicit schemes, particularly when the spectrum of a problem is non-real or strongly clustered, and for computationally demanding contexts such as quantum mechanics, stiff multiphysics, and nonlinear dispersive wave propagation (George et al., 12 Jan 2026, George et al., 2021, Buvoli, 2020, Casas et al., 2021).

1. Mathematical Formulation and Order Conditions

Complex time integrators generalize explicit and implicit schemes by permitting complex-valued increments and coefficients in one-step or multistep algorithms. For an initial value problem $dy/dt = f(y)$, a generic explicit one-step update with $\Delta t \in \mathbb{C}$ is

$$y^{n+1} = y^n + \Delta t \cdot \Phi(f, y^n, \Delta t)$$

where, for Forward Euler, $\Phi(f, y^n, \Delta t) = f(y^n)$. In the Runge–Kutta context, s-stage schemes allow complex Butcher tableaux $(a_{ij}, b_i)$ and stage nodes $c_i$, with updates of the form

$$k_i = f\Big(y^n + \Delta t \sum_{j<i} a_{ij} k_j \Big), \quad y^{n+1} = y^n + \Delta t \sum_{i=1}^s b_i k_i$$

and $a_{ij}, b_i, c_i \in \mathbb{C}$.

Complex time stepping can achieve higher order via contour selection. For instance, explicit Euler along n distinct complex substeps $w_i \Delta t$ (with $\sum w_i = 1$), matches Taylor coefficients of the exponential by solving algebraic constraints:

• $n=2$: $w_1 + w_2 = 1$, $w_1 w_2 = 1/2$ (giving $w_{1,2} = 1/2 \pm i/2$).
• $n=3$: $w_1 + w_2 + w_3 = 1$, $w_1 w_2 + w_2 w_3 + w_3 w_1 = 1/2$, $w_1 w_2 w_3 = 1/6$.

This enables arbitrary high order, including breaking the classical explicit Runge–Kutta order barrier (e.g., 5-stage, 5th-order explicit RK with complex coefficients) (George et al., 2021).

2. Stability Analysis in the Complex Time Plane

The stability of a complex time integrator is typically analyzed via the scalar test problem $y' = \lambda y$, with $z = \lambda \Delta t$. The amplification factor $R(z)$ describes the propagation of the solution:

• Forward Euler: $R(z) = 1 + z$
• RK2: $R(z) = 1 + z + \frac{1}{2} z^2$
• General complex Euler: $R(z) = \prod_{i=1}^n (1 + w_i z)$

The absolute stability region is $$\mathcal{S} = \{ z \in \mathbb{C} : |R(z)| \leq 1 \}$$. By judicious selection of $w_i \in \mathbb{C}$, complex time-integrators reproduce classic regions (e.g., RK2, RK3) and, crucially, can generate asymmetric or purely imaginary-axis-elongated regions. For linear Schrödinger equations, where $\lambda \in i \mathbb{R}$, the optimally-chosen complex coefficients yield stability regions double those of standard real methods, allowing for time steps up to $2/|\lambda_{\max}|$ (George et al., 12 Jan 2026, George et al., 2021).

3. Complex Time-Stepping for Quantum and Dispersive Evolution

Complex time integrators demonstrate unique advantages for evolutions whose spectral content lies off the real axis. For the linear Schrödinger equation, after spatial discretization ($u_t = -i \Lambda u$), the optimal forward Euler path uses complex weights $(w_1, w_2) = (1/2 \pm i/2)$, forming a stability polynomial with maximal $\|z\|$ along the imaginary axis. These methods maintain approximate unitarity and dispersive accuracy, showing no secular drift in $\|u\|^2$ even at large steps.

Block polynomial methods constructed with imaginary nodes (e.g., BDF/GBDF/Adams architectures with complex stencils) yield full imaginary-axis (A(90°)) stability, essential for dispersive PDEs such as nonlinear Schrödinger or wave equations. Geometric contour selection replaces classical order conditions, leading to high order (often $q+1$ or $q+2$ for block width $q$), while preserving strong stability (Buvoli, 2020).

4. Applications to Stiff Problems and Projective Integration

Complex time integrators are beneficial for ODEs/PDEs with stiff or multiscale structure. For systems with spectral clusters (fast $\lambda_f$, slow $\lambda_s$), projective integration (PI) traditionally requires small explicit steps to resolve the fast modes. By choosing complex inner-step sizes $\delta t = -1/\lambda_f$ and aligning the substep in the complex plane, one can move the fast eigenvalue into the damped region of the amplification function, allowing the outer step $\Delta t$ to be chosen based on slow modes only. This eliminates the small-step restriction, providing cost reduction by a factor $|\lambda_f / \lambda_s|$ and preserving first order global accuracy (George et al., 12 Jan 2026).

Multi-adaptive Galerkin methods (mcG(q), mdG(q)), while not strictly complex time-step schemes, employ componentwise step adaptation and explicit stabilization techniques resembling the effect of complex time pathways, particularly for systems with localized fast scales or reaction fronts (Logg, 2012).

5. Implementation, Storage Efficiency, and Parallelization

Complex integrators can confer storage and parallel advantages. For $n$-step complex Euler or complex-coefficient RK, only one solution vector (with real and imaginary parts) need be retained, in contrast to $s$ stage vectors for real $s$-stage RK. This reduction is significant for large-scale problems (George et al., 2021).

Higher-order complex integrators built as linear combinations of symmetric-conjugate compositions (the "T-methods") parallelize efficiently: $2^{k-1}$ palindromic compositions for order $2(n+k)$ can be executed concurrently, minimizing wall-clock time relative to real recursive triple-jump or composition techniques (Casas et al., 2021).

6. Benchmark Results and Performance Comparisons

Representative benchmarks illustrate the efficacy of complex time integrators:

• For nonlinear Schrödinger solitons, a complex time-stepping method achieves the same accuracy as optimized real methods with %%%%45$$y^{n+1} = y^n + \Delta t \cdot \Phi(f, y^n, \Delta t)$$46%%%% larger $\Delta t$, reducing computation time by nearly half (George et al., 12 Jan 2026).
• For stiff real systems (Prothero–Robinson, oscillators with large complex eigenvalues), complex inner steps confer stability where real-projective FE fails.
• Polynomial block methods with imaginary/symmetric nodes outperform BDF and ESDIRK schemes for both dispersive and parabolic equations (Buvoli, 2020).
• For multi-scale Klein-Gordon equations, tailored Magnus-based complex splitting maintains error constants uniformly in the oscillation frequency, avoiding the severe CFL-type restrictions typical of explicit real methods (Kropielnicka et al., 2021).

Example Problem	Method	Max $\Delta t$	Accuracy	CPU time
NLSE soliton	FE (real)	0.002	$1.1 \times 10^{-3}$	15 ms
NLSE soliton	Optimized complex FE	0.014	$7.1 \times 10^{-4}$	9 ms
Prothero–Robinson ($\lambda$ large)	cPFE	0.05	Stable	—
Viscous Burgers	BDF-SMFC $q=8$	—	$10^{-10}$ error in $1/2$ time of BDF6	—

7. Design Principles, Limitations, and Outlook

The primary advantage of complex time integrators lies in the additional geometric degrees of freedom for contour, weight, and node selection, which enable:

• Expanded or tailored stability regions aligned with non-real spectra
• Arbitrary high-order accuracy for single-stage explicit methods
• Relaxation of classical order barriers for Runge-Kutta-type integrators
• Storage efficiency via reduction in intermediate stage variables
• Parallelization owing to block and linear-combination formulations

A plausible implication is that further development of complex time-stepping, especially when coupled with adaptive error control and multi-component integration, will continue to expand the explicit integrator toolbox for quantum, stiff, and dispersive systems, reducing the reliance on fully implicit methods in many regimes (George et al., 12 Jan 2026, George et al., 2021, Buvoli, 2020, Casas et al., 2021). Limitations include increased complexity in algorithmic implementation, the necessity of robust complex arithmetic support, and potential round-off sensitivity when contour nodes stray far from the real axis.

Complex time integrators represent a conceptually and practically significant extension of traditional time-marching paradigms, offering versatile design and substantial performance benefits for a wide range of applied mathematics and computational physics problems.