Approximate Exponential Integrators for Time-Dependent Equation-of-Motion Coupled Cluster Theory

With a growing demand for time-domain simulations of correlated many-body systems, the development of efficient and stable integration schemes for the time-dependent Schrödinger equation is of keen interest in modern electronic structure theory. In this work, we present two approaches for the formation of the quantum propagator for time-dependent equation-of-motion coupled cluster theory based on the Chebyshev and Arnoldi expansions of the complex, nonhermitian matrix exponential, respectively. The proposed algorithms are compared with the short-iterative Lanczos method of Cooper et al. [J. Phys. Chem. A 2021 125, 5438–5447], the fourth-order Runge–Kutta method, and exact dynamics for a set of small but challenging test problems. For each of the cases studied, both of the proposed integration schemes demonstrate superior accuracy and efficiency relative to the reference simulations.


Introduction
In recent years, there has been renewed interest in the development of efficient numerical methods to study the quantum dynamics of correlated electrons in molecular and materials systems (see, e.g., Refs.1,2 and references therein).Under particular approximations, it is possible to circumvent the direct solution of the time-dependent Schrödinger equation (TDSE) in favor of time-dependent perturbation theory (or "frequency-domain" methods) which aims to implicitly access quantum dynamics through probing the spectral structure of the Hamiltonian operator.4][15] While these methods can often be a powerful tool for the simulation and prediction of observable phenomena such as spectroscopies, their veracity depends on the applicability of their various approximations to accurately characterize queried physical conditions.Further, the vast majority of these perturbative methods serve to access the equilibrium behaviour of electronic dynamics, leaving non-equilibrium phenomena, such as charge migration, 16 inaccessible.From a theoretical perspective, time-domain simulations do not suffer from these deficiencies and may be straightforwardly extended to non-perturbative and non-equilibrium regimes. 1,2ven the ability to faithfully represent physical conditions by a chosen Hamiltonian, wave-function ansatz, and initial condition, the primary challenges of time-domain electronic structure methods are practical rather than theoretical.In contrast to frequency-domain methods which trade the problem of temporal dynamics for the tools of numerical linear algebra, [17][18][19][20][21][22] time-domain methods require explicit integration of the TDSE, which is generally more resource intensive.For hermitian discretizations of molecular Hamiltonians, such as Hartree-Fock (real-time time-dependent HF, RT-TDHF 23,24 ), density functional theory (RT-TDDFT), 25 and configuration interaction (TD-CI), [26][27][28][29] significant research effort has been afforded to the development of efficient numerical methods to integrate the TDSE. 30,31 particular, approximate exponential integrators based on polynomial (Chebyshev 30,[32][33][34][35] and Krylov subspace (short-iterative Lanczos, 36 SIL) expansions of the quantum propagator are among the most widely used integration techniques for hermitian quantum dynamics.
Exponential integrators are powerful geometric techniques for the solution of linear ordinary differential equations (ODE), such as the TDSE, as they preserve their exact flow, 37 thereby allowing for much larger time-steps than simpler, non-geometric integrators such as the fourth-order Runge-Kutta method (RK4).In addition, these methods may also be formulated in such a way as to only require knowledge of the action of a matrix-vector product, 30,[38][39][40] thereby avoiding explicit materialization of the Hamiltonian matrix which is generally large for correlated many-body wave-functions.
The situation is significantly more complex for non-hermitian Hamiltonian discretizations such as those arising from coupled-cluster (CC) theory (see Ref. 41 for a recent review).Due to its simplicity and low memory requirement, RK4 is generally the integrator of choice for time-domain CC methods in the recent past. 41Symplectic, [42][43][44] multistep, 45 and adaptive 46 integrators for time-domain CC methods have been developed, and have yielded significant efficiency improvements over their non-symplectic counterparts.Exponential Runge-Kutta integrators have been explored in the context of nonlinear time-dependent CC theory (TD-CC), 47 but have yet to see wider adoption.Recently, Cooper, et al. 48suggested an approximate exponential integration scheme for time-dependent equation-of-motion CC theory (TD-EOM-CC) 29,41,43,[49][50][51][52] based on the hermitian SIL method to efficiently generate linear absorption spectra for molecular systems.Despite being only valid for hermitian matrices, the proposed SIL approach was demonstrated to produce sufficiently accurate spectra with relatively low subspace dimensions.However, the ability of this scheme to produce faithful, long-time dynamics within TD-EOM-CC has not been assessed, and is unlikely due to its hermitian ill-formation.In this work, we pursue the development and assessment of polynomial and non-hermitian Krylov subspace (short-iterative Arnoldi, SIA) methods for the complex matrix exponential to enable the efficient and accurate simulation of TD-EOM-CC.
The remainder of thie work is organized as follows.In Sec.2.1, we review the salient aspects of TD-EOM-CC theory relevant to the development of efficient exponential integrators.In Secs.2.2 and 2.3 we examine the properties of exact and approximate dynamics for the TD-EOM-CC ODE and present the developmed integration schemes based on the Chebyshev (Sec.2.3.1) and SIA (Sec.2.3.2) expansions of the complex matrix exponential.In Sec. 3, we apply the developed integration schemes to a set of small test problems and compare their verasity with exact dynamics as well as previously employed SIL and RK4 methods.We conclude this work in Sec. 4 and offer outlook on future directions for approximate exponential integrator development in TD-EOM-CC in the years to come.
2 Theory and Methods 2.1 Time-Dependent Equation-of-Motion Coupled-Cluster Theory Time-dependent equation-of-motion coupled-cluster (TD-EOM-CC) theory is a general timedomain reformulation of many-body quantum mechanics capable of simulating the dynamics of both time-dependent 29,41,49,50 and time-independent 51,52 Hamiltonians.In this work, we consider the moment-based formulation 51 of TD-EOM-CC to compute the spectral function, where ) is (the dual of) the time-dependent moment function which describes the propagation of weak perturbations throughout the many-body system.We note for clarity that, due to the nonhermiticity of the CC formalism, ⟨ M (t)| is not the complex conjugate of |M (t)⟩.Additionally, throughout this paper, we chose S(t) to be ⟨ M (0)|M (−t)⟩, although ⟨ M (t)|M (0)⟩ is also valid.|M (t)⟩ (⟨ M (t)|) may generally be described via a linear expansion of (de-)excitations from a reference state |0⟩ (typically taken to be HF), where m 0 ( m0 ), m a i ( mi a ) and m ab ij ( mij ab ) are time-dependent (de-)excitation amplitudes, c p (c † p ) is the fermionic annihilation (creation) operator associated with the spin-orbital p, and the indices i, j, . . .and a, b. . . .denote occupied and virtual spin-orbitals relative to |0⟩.In this work, we truncate Eq. ( 2) to only include up to double excitations from the reference, resulting in the TD-EOM-CCSD approach.
Within the TD-EOM-CC formalism, the moment excitation and de-excitation amplitudes obey the following set of coupled, linear-time-invariant (LTI) ODEs 51 and their left-hand counterparts where H N ∈ C n×n is the non-hermitian, normal-ordered, similarity-transformed Hamiltonian represented in the basis of Slater determinants. 10,11From the moment state-vectors, m(t) and m(t), S(t) of Eq. (1) may be evaluated as where we have taken m ≡ m(0).It is worth mentioning that the TD-EOM-CCSD formalism used here requires propagating only the right-or left-hand moment amplitudes [in this case, the right-hand amplitudes following Eq.( 4)].While Eq. ( 1) is perturbatively derived from Fermi's Golden Rule, 51 time evolution of |M (t)⟩ via Eq.( 4) also serves as a useful model for the development of both LTI and non-LTI integration techniques for TD-EOM-CC methods as it formally consists of the same algorithmic components that are required for the simulation of time-dependent Hamiltonians. 29,41,49,50en specified as an initial value problem, Eq. ( 4) admits an analytic solution where exp −iH N t is the quantum propagator and exp is the matrix exponential defined in the canonical way. 53We refer the reader to Refs.51,52 for discussions pertaining to the choices of initial conditions for Eq. ( 7) to simulate various spectroscopic properties.In this work, we consider the dipole initial conditions 51 induced by where T and Λ are the ground-state CC excitation and de-excitation operators (again truncated at double excitation/de-excitations in this work), and μ is a particular component of the electronic dipole operator.

Exact Matrix Exponential
When H N is small enough to be formed explicitly in memory, Eq. ( 7) may be directly evaluated as where Ω ∈ C n×n is the diagonal matrix of EOM-CC eigenvalues, Ω = {ω I ∈ C} n I=1 , and L, R ∈ C n×n are the full, biorthogonal set of corresponding left and right eigenvectors safisfying the equations 10,11 where I ∈ C n×n is the identity-matrix.As Ω is a diagonal matrix, exp(−iΩt) is simply the diagonal matrix with entries e −iω I t .Insertion of Eq. ( 9) into Eq.( 6) yields the following simple expression for the exact autocorrelation function As a non-hermitian matrix, H N is not guaranteed to have real eigenvalues if the manyelectron basis is truncated, and as such, Eq. ( 9) (and by extension Eq. ( 7)) is not guaranteed to be unitary (norm-preserving) and will generally yield dissipative or divergent dynamics along EOM-CC modes with ℑω I ̸ = 0 (see, e.g., a recent study in Ref. 54).However, it has been shown that, 55,56 barring suboptimal ground-state CC solutions or the presence of conical intersections, H N typically admits a real spectrum representing physical excited states and thus, Eq. ( 9) is unitary in exact arithmetic.

Approximate Exponential Integrators
While Eq. ( 9) is an exact solution to the LTI TD-EOM-CC dynamics considered in this work, it requires the full diagonalization of H N .As the memory requirement associated with the EOM-CCSD H N grows O(N 8 ) with system size, full diagonalization is impractical for all but the smallest problems.For some systems, it is possible to integrate the TD-EOM-CC equations in a subspace spanned by a small number of states such that full diagonalization is not required. 29,41,49,50However, if a large number of states are required or spectral regions of interest are densely populated or spectrally interior, this approach also becomes impractical.
Matrix exponentiation is a challenging numerical linear algebra problem, and the past half century has yielded a wealth of research into the development of efficient implicit 30,[38][39][40] and direct 53 methods both for hermitian and non-hermitian matrices.In this work, we will consider subspace approaches for evaluation of the complex, non-hermitian matrix exponential generally taking the form where V ∈ C n×k is a k-dimensional subspace (with k ≪ n) generated by the action of −iH N onto the current state vector, m(t), and c(δt) ∈ C k is a time-varying coefficient vector.Given the ability to implicitly form σ ← H N v (i.e. a "σ build"), which is a standard algorithmic component of any EOM-CC implementation, 10,11 the implementations of Eq. ( 12) considered in this work will not require materialization of H N in memory.Within the subspace ansatz, Eq. ( 6) becomes where w is time-independent for fixed V .
For a particular expansion order k and state vector m(t), Eq. ( 12) will generally be valid for |δt| ≤ |∆t|, where ∆t will be referred to as a macro time-step in the following.
Within this prescription, the total simulation length, T , will be partitioned into subintervals ]} where t 0 = 0, t i = t i−1 + ∆t i , and ∆t i is the macro time step for the i-th interval.The relationship between k and ∆t is method-dependent, and will be discussed for both the Chebyshev and Arnoldi integrators below.Due to the factorization of the time-dependence into c(t), a general property of truncated expansions such as Eq. ( 12) is in their ability to interpolate within each T i without requiring additional σ builds. 30is property is particularly advantageous for methods such as EOM-CCSD in which the computational complexity of σ formation scales O(N 6 ) with system size. 10,11For each T i , a single V is computed and the propagator may be interpolated to arbitrary temporal resolution by varying the corresponding coefficients.For each of the intermediate time intervals (i > 0), the approximation of m(t i+1 ) generated from the endpoint of T i is used as the starting vector to generate V for T i+1 .

Chebyshev Time Integration
The use of the Chebyshev expansion to evaluate the quantum propagator for hermitian 3][34][35] In this work, we demonstrate that this approach is also applicable to non-hermitian Hamiltonians with real spectra.Chebyshev polynomials of the first kind, {Φ p }, given by the recurrence are a powerful tool in the approximation of scalar and matrix functions on the real-line as they form the unique approximation basis which minimizes the uniform (infinity) norm on [−1, 1] at a particular order. 57In the Chebyshev basis, the TD-EOM-CC propagator acting on a general vector v may be exactly expanded as 30,32 exp where γ ± = 1 2 (ω max ± ω min ), ω min/max are the minimum/maximum eigenvalues of H N , δ k0 is a Kronecker delta, J p is the p-th Bessel function of the first kind, and HN is an auxilary matrix that scales the spectrum of the image of Φ p remains on the unit disk.Practically, HN need not be formed explicitly (see Alg. 1) and γ ± need not be computed from exact eigenvalues and can be approximated using standard techniques [58][59][60][61][62] as long as the mapped spectral bounds are contained in [−1, 1].
In practice, the sum in Eq. ( 15) is truncated to a finite order k, yielding a compact representation of the propagator in the Chebyshev basis, given by The truncation error at the interval endpoint (t + ∆t) of the Chebyshev expansion can be shown 63,64 to be bounded by For fixed argument, J p (z) is highly oscillatory for p < z but decays exponentially for p > z, as depicted in Fig. 1.We note that for even (odd) p, J p is an even (odd) function about zero.Therefore, for p sufficiently larger than |γ − ∆t|, we may approximate C(∆t) ≈ 2∥v∥|J p (γ − ∆t)|.Given a desired step size, ∆t cheb , and error threshold ε cheb , we may use this As ∆t cheb is fixed, T may be evenly partitioned into ⌈ T |∆t cheb | ⌉ intervals.The Chebyshev subspace vectors may be efficiently evaluated using only k σ-builds (Alg.1), thus the total σbuild cost for this method is ⌈ T |∆t cheb | ⌉•k.Another important aspect of the Chebyshev method is that, due to fact that the expressions in Eq. ( 16) are analytic, one need not materialize V cheb in memory.Instead, one may evaluate wcheb = V T cheb m (Eq.( 13)) directly as the Algorithm 1: Evaluation of Eqs. ( 13) and ( 15) via the Chebyshev Expansion subspace is generated, as is shown in Alg. 1, thus changing the memory requirement from O(kn) to O(3n).As it is often the case that one requires high-order Chebyshev polynomials (≫ 3) to accurately approximate the matrix exponential, this realization leads to a drastic reduction in memory consumption for large systems.

Short Iterative Arnoldi Time Integration
Considering the spectral decomposition of the exact propagator given in Sec.2.2, it is expected that the Chebyshev method discussed in Sec.2.3.1 will be most effective when Ω is nearly uniformly distributed within [ω min , ω max ], due to the fact that the Chebyshev basis minimizes the uniform function norm.If Ω is clustered, Krylov subspace techniques for the formation of the exponential propagator are often more effective. 38The basic principle behind Krylov approximation techniques for matrix-functions is rooted in the generation of a k-dimensional, orthonormal basis, where v 0 ∈ C n is an arbitrary vector with ∥v 0 ∥ = 1.Given V krlv , one may form a subspaceprojected Hamiltonian, and approximate the action of the matrix exponential as where e 1 is the first column of a k × k identity matrix and V krlv is the Krylov subspace generated from v 0 = v/∥v∥.Given that k ≪ n, the exponential in Eq. ( 20) may be efficiently evaluated via Eq.( 9).
For hermitian matrices, V krlv can be efficiently generated by the Lanczos iteration, 65 H krlv is a tridiagonal matrix, and both H krlv and V krlv may be formed implicitly via a simple three-term recursion.For the approximation of the propagator, this approach has come be known as the short-iterative Lanczos (SIL) method. 36Here, we present an analogous scheme for the exponential propagator based on the Arnoldi iteration, 65,66 which is a general Krylov subspace technique which extends to both hermitian and non-hermitian matrices.We will refer to this approach as the short-iterative Arnoldi (SIA) method in the following.Instead of a tridiagonal matrix, the Arnoldi method produces an upper Hessenburg matrix via the recursion where e k is the k-th column of the k × k identity matrix and β k+1 v krlv k+1 is the residual Algorithm 2: The Arnoldi Iteration Input: If H N were an hermitian matrix, H krlv would be tridiagonal and V krlv would span the same subspace as the one produced by the Lanczos iteration in exact arithmetic.
Much like the Lanczos iteration, H krlv may also be formed incrementally via the Arnoldi iteration as shown in Alg. 2. However, unlike the 3-term recurrence used in the Lanczos method, the Arnoldi iteration requires explicit orthogonalization of newly produced subspace vectors as opposed to the implicit orthgonalization generated by Lanczos.8][69] In this work, we have utilized the classical Gram-Schmidt method with reorthogonaliztation to perform the explicit basis orthogonalization. 70There exist non-hermitian extensions of the Lanczos method 71 which produce simultaneous, biorthogonal approximations for the left-and right-hand eigenspaces of non-hermitian matrices and have seen successful applications in both frequency domain CC applications 19 as well as in state selection for TD-EOM-CC. 50However, the biorthogonalization requirements of these methods can often be numerically unstable, [72][73][74] and as such, we expect the Arnoldi method to yield superior numerical stability in finite precision. 75 has been shown 38 that the error produced by Eq. ( 20) can be bounded by the right hand side of the following inequality where µ(H N ) is the largest eigenvalue of (H N + H † N )/2 and ρ = ∥H N ∥ 2 .Although tighter bounds can be found, 40 the bound given in ( 23) is more instructive.It shows that the approximation error made in an Arnold time integrator depends on the departure of V krlv from an invariant subspace of H N , which is measured by β k+1 , the step size or time window δt as well as the spectral radius of H N , measured by ρ and µ(H N ).
Unlike the Chebyshev method, where the expansion coefficients are known ahead of time, the coefficients for SIA are related to the spectrum of H krlv , which itself is dependent on v (the current state vector, m(t), in the context of Eq. ( 12)).As such, it is canonical to adopt a dynamic time-stepping approach where the Krylov subspace dimenion (k) is fixed before the simulation and each ∆t i corresponding to T i is determined dynamically throughout the time propagation.As Eq. ( 23) is only a loose bound, its practical ability to determine ∆t is limited.Given that the Arnoldi method produces successively more accurate Krylov subspaces with increasing k, a more practical error bound is given by c krlv k (∆t), which measures the potential for projections of the exact matrix-exponential onto vectors outside the Krylov subspace.Therefore, as has been successfully applied to the SIL method, 48 a reasonable choice for the step size is the largest ∆t such that |c krlv k (∆t)| < ε krylov , where ε krylov ∈ R + is a chosen error threshold.
Another side effect of the non-analytic nature of the SIA coefficients is that, unlike V cheb , V krlv must be materialized in memory and Eqs. ( 13) and ( 15) must be evaluated explicitly.As such, the memory requirement assococaited with SIA will grow O(kn) with basis dimension.
However, as will be demonstrated in Sec. 3, the SIA method will generally require fewer σ builds than the Chebyshev method to achieve commensurate integration accuracy.

Results
To assess the efficacy of the Chebyshev and SIA TD-EOM-CC integrators developed in this work, we compare the accuracy and efficiency of these methods for two test systems, N 2 (1.1 Å) and MgF (1.6 Å), relative to exact dynamics (Eq.( 9)) as well as RK4 and the TD-EOM-CC SIL method of Ref. 48.Each of these systems were treated at the EOM-CCSD level of theory with the minimum STO-3G basis set 76,77 to allow for practical comparisons with exact dynamics.All ground-state CC calculations were performed using a prototype Python implementation interfaced with the HF and integral transformation routines in the Psi4 software package 78 and geometries were aligned along the z-cartesian axis without the use of point-group symmetry.At their respective geometries, both of these systems exhibit realvalued EOM-CC spectra.All simulations in this work were performed using ε cheb = 10 −16 and ε krylov = 10 −6 (for both SIL and SIA) for a duration of T = 1350 E −1 h (≈ 32 fs).First, we examine the temporal error accumulation in the autocorrelation function (Eq.( 1)) using the normalized root-mean-square-deviation (RMSD) metric where S ex is given in Eq. ( 11) and δt is the temporal resolution of the integrated time series.
For the Chebyshev, SIA, and SIL integrators, δt = 0.05 E −1 h .As the temporal resolution and step-size coincide for RK4, we have compared our methods with 3 different RK4 stepsizes to illustrate convergence: RK4-1 (δt = 0.05 E −1 h ), RK4-2 (δt = 0.01 E −1 h ), and RK4-3 (δt = 0.001 E −1 h ).In the following, we will use E(T ) (i.e. the total accumulated autocorrelation error) as a global error metric to assess each integrators' relative accuracy.Figure 2 illustrates the accumulated autocorrelation error for each of the integrators considered.Parameters for Chebyshev (∆t cheb ), SIA (k), and SIL (k) simulations in Fig. 2 were selected to    minimize E(T ) for each method.For N 2 , the Chebyshev, SIA and RK4 integrators exhibit near constant error accumulation over the full simulation.SIL exhibits a sharp error increase between 1-10 E −1 h which is of the same order as ε krylov .For k = 36, SIA yields an invariant subspace up to an error of O(ε krylov ), and as such, the entire simulation (t < T ) can be performed using a single Krylov subspace.For MgF, SIL and RK4-1 diverge, while Chebyshev, SIA, RK4-2 and RK4-3 exhibit similar error accumulation characteristics as were observed for N 2 .However, unlike N 2 , SIA does not yield an invariant subspace even with largest subspace of k = 400, and thus multiple Krylov subspaces must be generated over the course of the simulation.As such, error O(ε krylov ) is compounded at each macro-time step, which explains the overtaking of SIA by Chebyshev in the long-t limit.
Figure 3 presents the cost-to-accuracy ratio, characterized by E(T ) as a function of σ builds emitted by each integrator, for a range of parameter choices.For N 2 (MgF), Chebyshev results were obtained for ∆t cheb ∈ {1, 5} (∆t cheb ∈ {1, 5, 10, 30, 50}).As discussed in Sec.2.3.1, the number of required σ builds for the Chebyshev is fixed at m cheb T /∆t cheb and m cheb generally increases as a function of ∆t cheb .This behaviour is shown explicitly for MgF in Fig. 4a.For both systems studied, neither E(T ) nor to the total number of  As is shown in Fig. 4b, the achievable time step (σ build count) subject to ε krylov is (inversely) proportional to k and thus the SIA and SIL data points in Fig. 3 are plotted in order of decreasing k.Unlike the Chebyshev method, the accuracy of SIA consistently improves with increased k, and thus k should be maximized subject to available memory resources to improve both accuracy and efficiency of the SIA method.
For N 2 , SIL results were also obtained with k ∈ {5, 10, 20, 36, 50} for a direct order-byorder comparison with SIA.At each order, SIA achieves between 2-3 orders of magnitude better accuracy over SIL, and requires >50% fewer σ builds in cases where SIA is able to take time-steps lager than δt (k ≥ 10).This is due to the fact that the Arnoldi method generates a faithful Krylov subspace representation H N while the Lanczos method, being only valid for Hermitian matrices, does not.This fact is particularly apparent in SIA's generation of an invariant subspace for k = 36 while SIL fails to demonstrate similar convergence.
For all problems considered, the proposed SIA and Chebyshev integrators exhibit superior accuracy and efficiency over analogous SIL and RK4 simulations.While it is possible for RK4  to yield reasonable accuracy at small time-steps (RK4-3), these simulations require excessive number of σ builds and would not be practical for the simulation of realistic TD-EOM-CC problems.

Conclusions
In this work, we have presented two approximate exponential time-integrators for TD-EOM-CC theory based on Chebyshev and Arnoldi (SIA) expansions of the quantum propagator.
The efficacies of these integrators were demonstrated via comparison with exact exponential dynamics for two small test problems.The Chebyshev and SIA integrators were demonstrated to yield superior accuracy and efficiency when compared to RK4 and the recently developed SIL method for TD-EOM-CC. 48As both of the presented methods are built from standard algorithmic components required for any implementation of (TD-)EOM-CC, the implementation of these methods has a low barrier for entry and holds the potential to yield significant performance and accurate improvements for these simulations in the future.
The practical application of the presented schemes requires consideration of the balance between desired integration accuracy and available computational resources.If memory capacity allows, the SIA method would be preferred for most chemistry applications due to its systematic improvability with respect to truncation order.However, the memory requirement of SIA quickly becomes prohibitive for large problems and the explicit orthogonalization requirement complicated efficient distributed memory implementations.In these instances, the Chebyshev method would be preferred due to its low memory requirement and the simplicity of its implementation.
While the results presented in this work have focused on the moment-based formalism of TD-EOM-CC, the presented efficacy experiments serve as an important proof-of-concept to demonstrate the the proposed methods for general TD-EOM-CC simulations.Future work to extend these methods to large scale TD-EOM-CC simulations is currently being pursued by the authors.Further, extension of these methods for use with time-dependent Hamiltonians, such as those required to study field-driven dynamics of molecular systems, are currently under development.

Figure 1 :
Figure 1: Graphical depiction of the order decay behaviour of Bessel functions of the first kind for fixed argument.The function is highly oscillatory for p < z but decays exponentially for p > z.

Figure 4 :
Figure 4: Assessment of the variance of cost and accuracy of (a) Chebyshev and (b) SIA integrators as a function of parameter selection.SIA results are presented as the average time-step ∆t as a function of k.