Approximation of Wave Packets on the Real Line

In this paper we compare three different orthogonal systems in L2(R)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{L}_2({\mathbb {R}})$$\end{document} which can be used in the construction of a spectral method for solving the semi-classically scaled time dependent Schrödinger equation on the real line, specifically, stretched Fourier functions, Hermite functions and Malmquist–Takenaka functions. All three have banded skew-Hermitian differentiation matrices, which greatly simplifies their implementation in a spectral method, while ensuring that the numerical solution is unitary—this is essential in order to respect the Born interpretation in quantum mechanics and, as a byproduct, ensures numerical stability with respect to the L2(R)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{L}_2({\mathbb {R}})$$\end{document} norm. We derive asymptotic approximations of the coefficients for a wave packet in each of these bases, which are extremely accurate in the high frequency regime. We show that the Malmquist–Takenaka basis is superior, in a practical sense, to the more commonly used Hermite functions and stretched Fourier expansions for approximating wave packets.


Introduction
Highly oscillatory wave packets play an important role in quantum mechanics.A case in point is the semiclassically scaled time dependent Schrödinger equation (TDSE), which describes the quantum dynamics of the nuclei in a molecule [17].It has the form where ψ(x, t) ∈ C, x ∈ R d and t ∈ R. The Laplacian is ∆ = This equation is important in a wide range of fields such as physical chemistry.However, solving TDSE comes with a number of computational challenges.A major reason is the small parameter ε 1, where ε 2 describes the mass ratio of an electron and the nuclei.Consider equation (1.1) without V (x) (the so-called free Schrödinger equation), then it can be shown that the solution to this problem generates oscillations of frequency O ε −1 in space and time [16].Moreover, the current techniques used to solve equation (1.1) are mainly defined on a finite interval.In that instance Fourier spectral and pseudo-spectral methods work well in practice [6].However, these approaches require the imposition of periodic boundary conditions, which may lead to non-physical behaviour over longer time frames, because the solution inevitably reaches the non-physical periodic boundary of the domain.Modern applications like quantum control also require much longer time frames.One solution to this problem, is to solve TDSE on the whole real line using spectral methods and this underlies the setting of this paper.
In designing a spectral method, we first need to pick basis functions that approximate the solution to equation (1.1).In this paper, we restrict our focus to systems that are orthonormal in L2(R) and have banded skew-Hermitian differentiation matrices.When the basis has a skew-Hermitian differentiation matrix, the solution operator is naturally unitary, so the Born interpretation of quantum mechanics is respected.As a result, these methods are also stable by design [9,11].Heller in 1976 proposed that the solutions to the semi-classical TDSE can be approximated by Gaussian type functions i.e. wave packets [10].Hagedorn proved error bounds for time dependent wave packets in 1980 [8]. Lubich in 2008, proved that for travelling wave packets with fixed width ∼ ε 1/2 that a representation of the solution by wave packets has an L2 error of O ε 1/2 for time intervals of size t ∼ 1 [18].Thus, wave packets are a good first approximation to the solution to semi-classical TDSE.In fact, for a quadratic potential with Gaussian initial data, they yield an exact solution to equation (1.1) [18].Thus, a major yardstick for the suitability of basis functions for equation (1.1) is how well they approximate wave packets.
A univariate wave packet function has the form f (z) = e −α(z−x 0 ) 2 cos(ωz), where α > 0, x0 ∈ R and ω ∈ R. The width of the wave packet is controlled by α, x0 is the position where the wave packet is centred and ω is the frequency of oscillations.For simplicity much of our analysis is restricted to x0, ω > 0 and this should be assumed unless stated otherwise.The expansion coefficients of a wave packet function, f (x), in an orthonormal set of L2(R) basis functions, {ϕn}n∈I where I is a countable index set such as Z, is given by the following inner product, Much of our analysis concerns the asymptotics of bn(ω) in the large-ω regime, from which results about an(ω) follow.Note we use the notation f (x) ∼ g(x) for x → ∞ to denote that functions f (x) and g(x) are asymptotically equivalent as x → ∞.This can be rewritten as f (x) = g(x) [1 + o(1)] [4].
• Malmquist-Takenaka (MT) functions ϕn(x) := i n 2 π In this paper, we derive an asymptotic estimate of the coefficients in each basis, which will allow us to compare and evaluate how well they approximate wave packets on the whole real line.For simplicity and without loss of generality we restrict our analysis to α, x0, ω, n > 0.
For stretched Fourier functions, the asymptotic estimate of coefficients is described in Theorem 2.1.In its more streamlined version, for any given ε > 0 it is possible to choose the stretching constant λ such that A specific choice that satisfies this bound is From this estimate we expect the logarithm of the magnitude of the expansion coefficients to resemble a hump, very similar to a sombrero, and we can observe these characteristics in Fig. 1.1.The sombrero-like hump tells us that the coefficients decay spectrally until the small constant ε eventually dominates the estimate.This drop off in decay rate is real, and for α = 1, x0 = 0 and ω = 50 in Fig. 1.1 occurs at about ε = 10 −20 , in line with equation (1.3).We have a pessimistic certainty that for all n satisfying, (π|n| − λω) 2 ≤ 4αλ 2 log ε −1 2 (2αλ/π) 1/2 − 1 −1 , we have |a n,λ (ω)| ≥ 2ε.The value of the stretching constant λ is discussed further in section 2. One clear advantage of using the stretched Fourier basis is that we can take advantage of the structure seen in Fig. 1.1, where we can confine our calculation to approximating N coefficients, the complement of which is guaranteed to satisfy |an(ω)| ≤ ε, by a single Fast Fourier Transform (FFT) in O(N log 2 N ) operations.2) for n > 0. We can also see the sombrero-like hump that is characteristic to stretched Fourier functions, which shows that our coefficients decay spectrally.The slow decay rate of coefficients (the brim of the sombrero) was purposefully chosen to start after ε = 10 −20 using (1.3).
In summary, stretched Fourier functions hold a great deal of promise; we have proved that the coefficients exhibit spectral decay and can be approximated using FFT.However, in a more practical setting, our method of determining the optimal cut-off of coefficients requires information about our wave packet which is typically unavailable.Moreover, a typical solution resembles an unknown linear combination of wave packets, whereas the nice performance e.g. in Fig. 1.1 is optimized for a single wave packet.This is discussed in more detail in section 2.
The next contender is the basis of Hermite functions [22,Table 18.3.1].Deferring the full asymptotic estimate to Theorem 3.1, a simplified version for n > 0 and α > 1/2 is , where Cosc is a bounded oscillatory component.It is difficult to deduce the 'shape' of these coefficients from this expression.However, to help convince the reader, let us set cn,ω = n/ω 2 to be constant -we see in section 3.2.3 that this is a reasonable assumption and note that this relationship keeps n and ω large.
Then we have • a factor of 1/ √ ω which gently drives the size of an(ω) down.
• A factor, greater than one, to the power n drives the size up geometrically.
• The argument of the exponential is essentially dominated by the competing ω and n term.If we replace n = cn,ωω 2 and note that for α > 1/2 this term is negative, we obtain rapid decay to zero.
Thus, overall we get geometric decay.The above expression matches the envelope of the an's exceedingly well.In Fig. 1.2, we can see that for large ω and small α and x0 we can clearly observe spectral decay.To further make the point that our asymptotic estimate captures the characteristics of the exact coefficients, the plot on the right is zoomed on the coefficients near the peak.
To understand how the number of coefficients grows with respect to ω, consider the following crude inequality, the number of coefficients which are greater than ε in absolute value are all the n's that such that where cn,ω is a constant, and α and x0 are fixed.Note that number of coefficients increases linearly with n but quadratically with ω (as n = cn,ωω 2 ).Fig. 1.3 displays the number of coefficients required for varying ε and fixed ω.
Although it is natural to think of Hermite functions in a quantum mechanical setting, since they are the eigenfunctions of the quantum harmonic oscillator, we show in Section 3, that the convergence properties leave much to be desired.The decay of Hermite function coefficients for highly oscillatory wave packets with small α and x0 is fast indeed, but as α and x0 increase, this soon ceases to be true in practical terms.
Out of our three bases, Malmquist-Takenaka functions are the strongest contender when approximating wave packets on the real line for two reasons: firstly, the asymptotic estimate in Theorem 4.2 demonstrates spectral decay and secondly, the coefficients can be approximated using FFT.To give a brief insight into our findings, the asymptotic estimate of coefficients is for ω ≥ 0. This holds for large |ω| regime and uniformly for n satisfying 1/4 ≤ |n/ω| ≤ c for any given constant c (this is made more precise in Section 4).We can immediately deduce exponential decay with respect to n within this regime relative to ω.This is be confirmed experimentally in Fig. 1.4.This bound highlights how the presence of high oscillations can be an advantage.In the case of no oscillations, ω = 0, Weideman found in [26] that the coefficients decay sub-exponentially.However, counterintuitively, as ω > 0 grows, the estimate in (1.4) exhibits exponential convergence for n ≤ c • ω for some c.In other words, the presence of oscillations means that we see exponential convergence within a domain and the range of the domain increases for larger values of ω.In the context of practical computation, for very highly oscillatory wave packets the coefficients will essentially exhibit exponential decay within the practical range of accuracy.
Our asymptotic estimate shows that we should expect our coefficients to decay spectrally.From (1.4), we can readily deduce the location of the peak which is given by Interestingly, the location of the peak does not depend on how much we stretch the wave packet.The number of coefficients n > 0 that are greater than ε in absolute value is asymptotically equal to That is the number of coefficients required to resolve a wave packet in MT basis grows linearly with ω, which beats Hermite functions where the number of coefficients grow quadratically in ω! Overall, MT functions seem to provide the most practical approach in approximating wave packets.The coefficients can be approximated rapidly using the FFT and we have proved that they exhibit asymptotic exponential decay.A more in-depth discussion about these bases can be found in section 4.
Detailed derivation of the above asymptotic estimates is discussed in the sequel.

Stretched Fourier basis
Stretched Fourier functions are the familiar Fourier functions scaled to the interval [−λ, λ], where λ > 0 can be chosen based on the function f to be approximated and a desired error tolerance ε.Fourier functions are the natural basis in the presence of periodic boundary conditions in a compact interval, but in our setting we are truncating a function which is defined on the whole real line to the interval [−λ, λ], and must contend with an error introduced by the Gibbs effect at the endpoints.The basis functions, for a subspace of L2(R), have the form Figure 1.4: MT functions: Equation (1.4) suggests that for x 0 = 0 (left plot) we expect the rate of decay of coefficients to resemble a straight line and this agrees with the exact coefficients.Note that to aid the reader, we have removed the oscillations in our estimate in the left plot.As we increase x 0 (right plot), we observe that the coefficients form a hump.Our estimate is represented as a dark dashed line and the exact coefficients are denoted by a light solid line.
For a wave packet, f (x) = e −α(x−x 0 ) 2 cos(ωx), where α > 0, x0, ω ∈ R, an integral expression for the expansion coefficient is found by substituting (2.1) into (1.2) to get It is clear that a set of N of these coefficients can be approximated by a single FFT in O(N log 2 N ) operations.
Dietert and Iserles in their unpublished report [5] showed that stretched Fourier expansions for functions of the form f (x) = g(x) exp(−αx 2 ) (where g(x) is an entire function which grows at most exponentially) converge at spectral speed down to a specified error tolerance, provided that care is taken when choosing λ.An appealing result shown therein is the 'sombrero phenomenon', whereby the logarithm of the absolute value of the expansion coefficients is shaped like a sombrero.For the special case of the wave packets considered in this paper, we apply their method to arrive at the result in Theorem 2.1, and give an example to demonstrate the sombrero effect in Figure 2.1.
Before proving the theorem, let us discuss this sombrero phenomenon in more detail.We do not expect the reader to be familiar with millinery terminology, so allow us to explain that the main part of a hat in which the head sits is called the crown and the projecting edge is called the brim.Drawing the analogy with a sombrero seen in Figure 2.1, the crown of the sombrero exhibits super-exponential decay from the central peak and the brim of the sombrero is where the decay rate is very slow -specifically, O(n −1 ), as predicted by standard Fourier theory for a function of bounded variation.Note however, that if λ is sufficiently large, as it turns out to be the case in Fig. 2.1, it is also a region where the magnitude of the coefficients is smaller than the accuracy we need for practical computation.Thus we can truncate the series for an approximation scheme and attain effective spectral convergence down to our desired accuracy.
Proof of Theorem 2.1.Essentially, we evaluate the integrals in Lemma 3 and 4 from [5] exactly.Consider the modified coefficient given by an integral over the whole real line, which is much more amenable to direct calculation: From this we have the bound, A bound for the original coefficient can then be found as follows.
In the last step we have used the known inequality, which holds for t ≥ 0. Substituting the bound for |ǎ n,λ (ω)| obtained above yields the result of the theorem.

The choice of stretching parameter λ
It was found in [5] that the position of the sombrero brim depends on our choice of λ.In order to choose an appropriate λ, we need to understand the characteristics of our bound in (2.3).It requires us to choose λ ≥ |x0|, if we instead choose λ less than |x0| then we obtain a poor approximation to the original wave packet.This can be seen in Fig. 2.2, where the wave packet has much of its 'mass' outside the interval of approximation.The bound also suggests that we are limited by the magnitude of exp −α(λ − |x0|) 2 .If this quantity is not smaller than some desired threshold ε then the coefficients are not guaranteed to decay spectrally below O(ε), but instead we expect them to decay like O n −1 , until they reach such a threshold.This can be seen in There are many effective ways to define λ; for a given > 0, a simple choice is, in order to guarantee Fig. 2.3 demonstrates the dependence of the wave packet coefficients on the choice of λ.If we choose λ to be smaller than λ optimal we would need to choose n = O(ε −1 ) = O(10 20 ) in order to reach our desired accuracy, which is unrealistically large.For λ larger than our optimal value the brim of the sombrero lies lower but the thickness of the crown is greater and we need to choose large (at least, only linearly) n to reach 10 −20 accuracy.Thus even though we get spectral decay with stretched Fourier bases, this approach is extremely sensitive to how well we can define the optimal λ.This effect can also be seen theoretically in Theorem 2.1.
In a more practical setting, a linear combination of wave packets is a more realistic approximate solution to the semi-classical TDSE and due to our current method of determining λ optimal we need more knowledge about the exact solution before we can approximate it.Our analysis indicates that it is safer to overestimate λ than to underestimate it, providing a clue to a good strategy for a linear combination of wave packets.The current practical limitations of stretched Fourier basis notwithstanding, the functions still hold significant potential due to their simplicity, fast approximation of coefficients using FFT and spectral rate of decay.

Number of coefficients required for a given tolerance
Given a tolerance ε > 0, once we choose λ as in equation (2.4) we can guarantee that the coefficients satisfy (2.5).To gain a rough idea on the number of coefficients required for a given ε, we do the following Neglecting the complicated term inside the logarithm (i.e.bounding this inequality from above), we have and, substituting λ = λ optimal from (2.4), we conclude with which emphasises the role different parameters play in the bound.The important factor to note here is the linear dependence on both ω and |x0|, which is extremely favourable compared to the other two bases considered in this paper.The downside, though, is that we require prior knowledge of x0, α and ω to choose λ and that the approximation of a linear combination of wave packets is likely to result in less favourable behaviour.

Hermite functions
It is natural to think of Hermite functions in the context of quantum mechanics since they are solutions of the quantum harmonic oscillator.Hermite functions have a skew-symmetric and tridiagonal differentiation matrix, which implies that their associated spectral methods are stable and preserve unitarity [9].However, computation of coefficients in an expansion in Hermite functions is associated with a number of practical challenges.An obvious method of computing expansion coefficients in Hermite functions with quadrature is by explicit integration of each coefficient at the overall cost of O(N 2 ) operations.Nonetheless they are still a contender when approximating wave packets and, less intuitively, it is possible to compute the expansion coefficients in O(N (log N ) 2 ) operations using the fast multipole algorithm [2].
The behaviour of the coefficients is hard to deduce by a direct examination of this expansion.It is helpful to break down the result into three regimes.Note that in a more practical setting we would expect α > 1/2, so the following analysis considers the asymptotics for this case: , where Cosc is a bounded oscillatory term.The asymptotic expression only depends on ω and we expect exponential decay.
2. When n is large, .
The expression is still fairly opaque but in Fig. 3.1 we can see that the envelope of the logarithm of the absolute value of coefficients is approaching a straight line.3. When both n and ω are large subject to a power law, n = ω λ for some λ > 0, then

Fourier-type integral of Hermite functions
The main idea behind the proof is to consider a Fourier type integral of Hermite functions,

Lemma 3.2. It is true that
The branch choice in eqn , is unimportant because only the even powers of √ 2c − 1 occur on the right-hand side due to symmetry of the Hermite polynomials.
Proof.Consider Hermite functions in a Cauchy integral form, where the contour γ encircles the origin.Substituting this into (3.3) and changing the order of integration, For c = 1/2, the proof of (3.4) is identical (and easier!) to the steps above: we deduce that We proceed by comparing (3.3) with (3.2) to deduce that and After some algebraic manipulation (see A.1) and setting where for simplicity we restrict our analysis to α, x0, ω > 0 so that X and Y are positive real numbers, we note that our expressions are different for different values of α, × Re e i|2α−1|XY +inπ Hn(Y − iX) . (3.8) Thereafter we focus on the asymptotics for α = 1 2 , the remaining case can be also obtained by a limiting process.Note that (3.7) are the exact coefficients for α = 1/2 and the plot of (3.7) can be seen in Fig.

The method of steepest descent
The method of steepest descent (MSD), also known as the saddle point approach, is a key tool used in our Hermite and Malmquist-Takenaka estimates.We give a brief overview here; for a more in-depth review we refer to [4,7,21,29].It is an approximation method that can be applied to integrals of the form where t is a large parameter.Let us assume that F (z) is a complex-valued meromorphic function and G(z) is an entire function.
The main idea is to deform the path of integration over a complex contour Γ to pass through (or near) certain saddle points of F (z), that is, points z * such that F (z * ) = 0.As well as passing through a subset of saddle points, the path of integration is also designed so that Im F (z) is constant, i.e.Im F (z) = Im F (z * ), which ensures that exp(−tF (z)) is non-oscillatory along the contour.If F (z * ) > 0 then exp(−tF (z)) decays rapidly away from z * along this contour -the path of steepest descent.We assume at the outset that G(z) can be expanded as convergent power series about the neighbourhood of z * .
Skipping many technical details, suppose that there exists an ordered set of contours Γ1, Γ2, . . ., Γs such that each contour contains a single saddle point (z * k ∈ Γ k ), within each contour F satisfies Im F (z) = Im F (z * k ), each saddle point in question satisfies F (z * k ) > 0, G(z) can be expanded as a convergent power series about each z * k , and such that b a e −tF (z) dz = z) dz.Then as t → ∞ and G(z * k ) = 0, to first order approximation our integral can be expressed as ([4, eq.5.7.2]).
Remark 3.3.The power of this approach is that the first approximation depends only on the values of F , F and G at a certain subset of the saddle points, and nowhere else.Note that the terms hidden in the O error term in (3.9) is a polynomial in G /G, G /G, F , F (IV ) and 1/F , see [23, eq.2.4.16],evaluated at a subset of saddle points.These terms depend on z and κ, where κ is a carefully chosen parameter with some dependence on t -more importantly they must remain controlled as t → ∞ thus imposing constraints on κ.This ensures that (3.9) holds uniformly as t → ∞.

Steepest descent: Hermite functions
To find the asymptotic behaviour of (3.5), we consider Hn(z) for z given in (3.6) and (3.8),where n is a large parameter and Y is large in the complex argument, z.We start by considering Hermite functions in a Cauchy integral form, where γ is a contour encircling the origin.Let ζ = z/ν, so that where φ(w) = w 2 − 2ζw + 1 2 log w, ν = √ 2n + 1.The saddle points are the zeros of φ , namely w± = Deforming our path of integration to pass through the dominant saddle point, w * , then applying (3.9), gives us . This asymptotic approximation is valid for large values of n and ω.For small values of ω a similar contribution from the remaining saddle point should also be included.
By remark 3.3, in order for the estimate in (3.10) to hold uniformly with respect to ζ for ω → ∞, we require the existence of c > 0 such that c < |w * | < c −1 , which is guaranteed as long as ζ is bounded.We also require |ζ − 1| > δ, δ being a small positive number, to avoid the coalescing of saddle points, in which case we would need to use uniform asymptotic methods, see [25, §23.4].This leads to the combined conditions δ < |ζ − 1| < δ −1 .
We have the key ingredient to find the asymptotic expansion of the Hermite coefficients.For α > 1/2, we want to apply (3.10) to Hn(z) in (3.6), the dominant root is w− and w−φ Recall that for this case, ζ = (X + iY )/ν, where ω is a large parameter, while α and x0 are fixed.For α < 1/2, we want to do the same, apply (3.10) to Hn(z) in (3.8), w+ root accounts for the main contribution to the asymptotic estimate and w+φ (w+ As highlighted in remark 3.3, care must be taken to constrain ζ so that our approximation in (3.10) holds uniformly as ω → ∞.For α = 1/2 and large ω, .
We have two competing large parameters n and ω, and we want to find a substitution which preserves the relationship between the two terms, This satisfies the condition in remark 3.3 for ζ to remain bounded for large n and ω.The additional condition in the remark for non-coalescing saddle points, δ < |ζ − 1| is equivalent to

α > 1/2: Deriving an asymptotic estimate
In a more practical setting, we would expect α > 1/2 and we commence by finding the asymptotics for this case.Equation (3.10) is instrumental in deriving first-order asymptotics when large degree coexists with large argument: higher order estimates can be obtained from [20], but we do not pursue this theme further.
We observe that for reasonable values of x0 we can simplify things by requiring Y X so ζ ∼ iY /ν.However, for more accurate asymptotics we want to retain x0 explicitly in our expansions.As discussed in the last section, to satisfy remark 3.3, (3.12) must hold for ζ to remain bounded for large n and ω.Moreover, as α > 1/2 we have 2(1−4α 2 ) < 0, so condition (3.13) is automatically satisfied.Thus the MSD approximation in (3.10) holds uniformly for n > cn,ωω 2 for any given constant cn,ω.The relationship is demonstrated in our numerical experiments, shown in Fig. 3.4.
Applying (3.10) to (3.6) (with identical notation), we have where  Setting n = cn,ωω 2 where cn,ω is some constant, then expanding about ω we arrive at Substituting A and B back into an and restoring cn,ω = n/ω 2 , we obtain the asymptotic estimate of coefficients, 1+O n −1 .

α < 1/2 asymptotics
We have used the same method as for α > 1/2 so we have omitted some technical details and the bulk of the algebra has been moved to the appendix.
We proceed by applying, (3.10).For α < 1/2 e −ν 2 φ(w + ) 2π w+φ (w+) Note that for this case we apply the MSD on Hn(ζ) for ζ = (Y − iX)/ν where our contour of integration passes from the saddle point w+.Like before, we can rewrite an in terms of A and B, the slow and fast components respectively, Here, we have substituted w+φ (w+) and φ(w+) in terms of X and Y , into A and B (see section A.3).In (3.18) we have used Stirling's formula.
From equations (3.17), (3.19) and (3.20), the key term we need to approximate is Here we have substituted n = cn,ωω 2 , like before, and we want to expand around ω. First thing to note is that the large ω term is given by Observe that there is a change in behaviour at .
This tells us that for a small neighbourhood around c * n,ω , we would expect this region to have different asymptotics to the the region away from it.This agrees nicely with our constraints discussed at the start of this section, where we have identified that the non-coalescing saddle point regime in remark 3.3 is equivalent to considering cn,ω away from c * n,ω .Going through the same steps as in the previous section and sparing the reader more algebra (see section A.3.1 for full derivation), for cn,ω = c * n,ω we obtain The asymptotics are displayed in Fig. 3.7, which confirms the remarkable fit of our formulae with computed values of the coefficients.

Malmquist-Takenaka functions
MT functions have a long history.Discovered simultaneously by Malmquist [19] and Takenaka [24] in 1925, they have since been studied in fields ranging from signal processing (Wiener in 1949 [28]) to spectral methods (Christov in 1982 [3]).In [12], Iserles and Webb proved that this is essentially the only orthonormal and complete system in L2(R) with a banded tridiagonal skew-Hermitian differentiation matrix, and the coefficients can be computed using the Fourier transform.This is part of Iserles and Webb's ongoing work classifying bases with banded skew-Hermitian differentiation matrices [13,14,15].It was the work of Boyd in [1] and Weideman in [26,27] that highlighted the interesting approximation properties of MT functions.The approximating characteristics of MT functions do not obey simple rules.As an example, for f (x) = 1/(1 + x 4 ) the coefficients decay at a spectral rate of an ∼ O (1 + √ 2) −|n| , however if we include an oscillatory factor so that our function is in the form f (x) = sin(x)/(1 + x 4 ), our approximation collapses from spectral to algebraic convergence, an ∼ O |n| −9/4 , see Fig.In [1] and [27] (respectively), Boyd and Weideman proved that MT expansion coefficients decay exponentially if and only if f (z) decays at least linearly as |z| → ∞ in the complex plane, and f is analytic in a region which may exclude ±i/2 but which must include the point at infinity.Analyticity at infinity is a strong condition to impose, since sin z, for example, blows up as z → ∞ along the imaginary axis.This is not unexpected because the proof of exponential convergence in a compact interval (e.g. with orthogonal polynomials) depends on functions being analytic within a Bernstein ellipse surrounding this interval.
Interestingly, Weideman (in [26]) found that the coefficients, an, of a Gaussian function asymptotically decay close to exponential convergence, an ∼ O e − 3 2 |n| 2/3 .Note, we do not have exact spectral decay as the Gaussian function has an essential singularity at ∞.What does this mean for wave packet coefficients?We proved that for wave packet functions the coefficients exhibit asymptotic exponential convergence.In other words as ω gets larger, we get closer to exponential convergence up to some threshold.For example, if f (x) = e −x 2 cos(ωx) the coefficients decay like an ∼ O e −n/ω .This can be seen in Figure 4.2, where even though the cos(50x) factor decreases the rate of convergence compared to just the Gaussian function, we still see asymptotic spectral decay, at least in the regime |an(ω)| ≥ 10 −10 .
To see how we have proved this, let us start by defining our system.From [12] the Malmquist-Takenaka system is defined as Substituting (4.1) into (1.2) yields the expansion coefficients where where α = βω, n = cn,ωω and we take the principal value of tan −1 .We initially restrict ourselves to the case n ≥ 0, ω ≥ 0. The remaining cases will be considered in section 4.3.

4.1
The case ω ≥ 0, n ≥ 0: the method of steepest descent This section draws inspiration from Weideman in [26], who found the rate of decay of MT coefficients for the Gaussian function, e −x 2 by applying MSD to the coefficient integral.Our analysis will be considerably more involved than that of Weideman's for the Gaussian function, because we have two large parameters, n and ω.A brief explanation of the MSD can be found in section 3.2.1.Since β = α/ω, it is bounded for large ω.We stipulate that n and ω vary so that cn,ω is bounded: details will be discussed further in the next section.The main step of the MSD is to deform our integral in (4.3) over a complex contour that passes through (or near) a saddle point of (4.4).This represents the path of steepest descent.We can find the saddle points by solving with some manipulation i.e. multiplying through by 1 + 4z 2 then collecting z terms, this reduces to the cubic Our formulation in (4.3) has a nice symmetry aspect that comes from the tan −1 .For example, when β = 0 there are two saddle points As α is not an asymptotic parameter, β will be small for large ω, and the two shown values z± in (4.6) will be close to two exact saddle points, while the third one approaches infinity when β → 0.
We have for two saddle points z2 and z3 expansions where cn,ω should be bounded away from 1  4 and We obtain this from substituting (4.7) into (4.5),then collecting β terms to find the c k s, the coefficients c k s at saddle point z * ∈ {z1, z2, z3}.For the third saddle point, z1, we use a property of the cubic equation.Writing the equation as a3z 3 + a2z 2 + a1z + a0 = 0, and z3 has the expansion The two saddle points z2 and z3 become nearly equal when cn,ω → 1 4 .That is, there is a special case, when ω ∼ 4n.In this case, we cannot use the MSD and will have to use uniform asymptotic methods -this method can be found in [21,25,29].We do not consider this case in the paper and assume that cn,ω ≥ 1 4 + δ where δ > 0 is small.This is a necessary constraint to satisfy remark 3.3 and we discuss this further in the next section.

Large ω asymptotics
Our main goal is to find out how many MT functions are required to approximate a wave packet with a frequency of ω to prescribed accuracy.We want to estimate the number of {an : |an| = |bn(ω) + bn(−ω)| ≥ ε} for any given ε > 0. Our asymptotic expansion requires a subtle approach as both n and ω are large parameters.
We have learnt from numerical experiments that it is a poor idea to fix one parameter and let the other become large i.e. fix ω and let n 1, as the asymptotics from this regime miss on the spectral decay behaviour that we wish to capture.A simple explanation for this alludes to the relationship between n and ω in the square-root term in equation (4.6).
We found that there was an intermediate range where n and ω followed a relationship of the form for any given constant c.This range satisfies remark 3.3 and ensures that the MSD approximation holds uniformly for ω → ∞.Writing (4.7) and (4.8) in orders of ω, i.e. substituting β = α/ω, the roots have the form We require cn,ω > 1 4 so cn,ω − 1 4 is always positive.For α, x0, n, ω ≥ 0, z3 stays in the third quadrant, z2 can be located in the first or fourth quadrant and z1 is located in first quadrant of the complex plane.Moreover, if x0 = 0, then z2 would only be in the fourth quadrant.From z1, we see that the root is shifted from the pure imaginary axis by x0.By requiring Re z2 > x0, we get cn,ω > 1  4 + x 2 0 and past this point we see that Im z2 < 0, so in this regime z2 is confined to the fourth quadrant.
The plot below suggests that z2 and z3 may be appropriate saddle points through which to deform the contour., z 2 is the reflection of z 3 on the y axis.Right plot: Plot of saddle points when x 0 = 2. z 1 is shifted x 0 = 2 away from the imaginary axis and z 2 travels from the first quadrant to the fourth quadrant as n increases.

Paths of steepest descent
The path of steepest descent is chosen so that Im g(z) is constant and passes through the saddle point z * of g(z).Recall from equation (4.4) that Substituting z = x + iy, and using the identity the imaginary part of g(z) has the form (4.12) We choose the path of steepest descent to be the path where Im g(z) = constant, the constant being chosen so that the path passes through the saddle point z * i.e. path = {z : Im g(z) = Im g(z * )}. Figure 4.4, shows a plot of steepest descent paths, equation (4.12), for certain parameter values.To construct our steepest descent contour Γ, we split Γ = Γ (3) ∪ Γ (2) , where Γ (2) is the contour that passes through z2, i.e. the integral that goes from 0 to ∞.Note that for n > 0 (4.2) has a pole at i/2 and a zero at −i/2 and it is sensible to avoid the pole in the upper half plane.From our discussion on the location of our roots in the last section, z2 can be in the first or fourth quadrant.For cn,ω > 1  4 + x 2 0 , z2 is in the fourth quadrant.For simplicity we will consider the root to be located in the fourth quadrant, but the same principle can be applied when the root is in the first quadrant.We start by defining contour Γ (2) at −i/2: • This was chosen simply because we want to avoid the pole at i/2.We can show that −i/2 is always on our steepest descent path by substituting −i/2 into (4.12)which gives Im g(−i/2) = βx0 = αx0/ω where α, x0, ω > 0 with no dependence on n.For a visual representation see the bottom plot of Fig. 4.4, where the steepest descent contour through z2 passes through both the pole and the zero.We can simply choose to construct our contour from −i/2 as there is nothing to look for along the path from 0 to −i/2.Moreover, the point z = −i/2 is the central point of a valley which is why the path of steepest descent may visit this point.
• By starting our contour from −i/2 means our curve looks nice and smooth.The imaginary part converges to ω/2α whilst the real part tends to infinity.As the real part of our contour tends to infinity, the integrand tends to zero very fast.Although this step is not necessary, we can bring the tail back to the real line by travelling straight down, only varying the imaginary part of our path, to reach the real line and then carrying on towards infinity.MT functions: Plot of (4.12) for given α, ω, x 0 and n.The ×s on the imaginary axis represent the pole at i/2 and zero at −i/2 for n ≥ 0. It can be seen in the top plots that the contour approaches an asymptote that is independent of x 0 .In the bottom plot, we see a zoomed in plot of the top left plot, the steepest descent contours pass through the pole at i/2 and zero at −i/2 for n ≥ 0.
The same idea can be applied to the remaining integral, Γ (3) , from −∞ to 0. For a schematic representation see Fig. 4.5.
The key argument used in constructing Γ (2) is that the imaginary part of the path of steepest descent converges to ω/2α as Re z → ∞.Note, this also eliminates the choice of z1 as Im(z1) > ω/2α which can be seen from eqn (4.9).The same idea holds for Re z → −∞.We can observe this in Fig. 4.4 for a given α, ω x0 and n.Lemma 4.1.For sufficiently large x there exists a function y = γ(x) such that x + iy ∈ Γ.Furthermore, limx→±∞ γ(x) = ω/2α.Proof.Our steepest descent contour satisfies where for some fixed n and ω the RHS is constant.Let z * = c + id, where c and d are constants.We have 0 We retrieve the γ(x) by the implicit function theorem.Dividing through by x and letting x → ±∞, Taking the limit we are left with y = 1/(2β) = ω/(2α), (where β = α/ω) as required.Once the desired path is chosen, we can apply the MSD as discussed in section 3.2.1.The integral in (4.3) can be rewritten as e −ωg(z) dz 1 + 2iz + Γ (2)   e −ωg(z) dz 1 + 2iz , where the path of integration Γ (3) starts at −∞, passes through the saddle point z3 and terminates at −1 2 i, and Γ (2) starts at that point, runs through the saddle point z2 and terminates at +∞.The transformations . (4.14)In the transformations in (4.13) we assume that z on the infinite part of the steepest descent path from −∞ to z3 corresponds to a nonpositive t-value, and that z on the finite part from z3 to − 1 2 i with a nonnegative t-value.Similarly for the integral through z2: the points z on the part from − 1 2 i to z2 correspond to t ≤ 0, the points z on the part from z2 to +∞ correspond to t ≥ 0. In this way there is a one-to-one relation of points z ∈ Γ ( * ) and t ∈ R. Observe that on the paths of steepest decent the quantities g(z) − g(z * ) are real and non-negative and their graphs are convex, like the graph of 1  2 t2 on R.

Substituting the Taylor expansion
where the odd k terms disappear as we are integrating an odd integrand over the real line.The integrals now look like Gamma integrals of the form Let y = ωt 2 /2, y −1/2 / √ 2ω dy = dt and t 2k = (2y/ω) k , then where 1 2 k is the Pochhammer symbol, (x) n = Γ(x + n)/Γ(x).By substituting (4.16) into (4.15)we obtain as ω → ∞, and large n = cn,ωω such that cn,ω ≥ 1 4 + δ, with δ a small positive number.The next step is to determine the a 2k 's, the full derivation can be found in section B.1.The first coefficients have the form , a (2) .

Asymptotic expansion
To complete the asymptotic expansion, we substitute (4.10), (4.11) and β = α/ω into the integrand in (4.18) and expand for large ω.This yields It follows from the first term in (4.19) that, as x0 increases, e −ωg(z3 ) gets smaller than e −ωg(z 2 ) .Substituting the terms inside (4.18), the following asymptotic expansion holds for cn,ω The final line of the above equation vividly demonstrates that our asymptotics exhibit exponential decay.This can also be seen in appendix B.2 where we find a bound for |bn(ω)| which leads to the result in (4.21).All that remains is to consider the n < 0 and ω < 0 cases.

Reduction to the case n ≥ 0 and ω > 0
At the outset we have assumed n, ω > 0, in this subsection we demonstrate that this can be extended to n ∈ R, ω ∈ R.

A symmetry argument
Using this symmetry we deduce the results for n ≤ −1, ω ≤ 0 from n ≥ 0, ω ≥ 0 for all n ∈ Z.

Laguerre coefficients argument
We can further reduce our range by considering the case where n ≤ −1, ω > 0 by using the relationship of MT functions to Laguerre polynomials, considered in [12].Taking an inverse Fourier transform of our wave packet e It is proved in [12] that for n ≤ −1 the MT functions are Therefore, substituting into the expression for the coefficients in terms of e iωx−α(z−x 0 ) 2 and ϕn, we get bn(ω) = −i n 2 As both Ln(ξ)e −ξ/2 and e ix 0 (ξ−ω) have L2(R) unit norm, by Cauchy-Schwartz, Thus by the symmetry argument in (4.20), this result holds for ω ≤ 0, n ≥ 0.
As we have shown in our crude bound above, the coefficients for bn(ω) are small for ω ≥ 0 and n ≤ −1.Moreover, we can use the above symmetry argument to evaluate bn(−ω).
From this point we can see that we obtain the bound in (4.21).

Figure 1 . 1 :
Figure 1.1: Stretched Fourier functions: Wave packet coefficients in a stretched Fourier basis.The solid line represents the exact coefficients and the (light) dashed line is a plot of the spectrally decaying part of our estimate in (2.2) for n > 0. We can also see the sombrero-like hump that is characteristic to stretched Fourier functions, which shows that our coefficients decay spectrally.The slow decay rate of coefficients (the brim of the sombrero) was purposefully chosen to start after ε = 10 −20 using (1.3).

Figure 1 . 2 :
Figure 1.2: Hermite functions: The left figure shows the exact coefficients in the lighter colour and the asymptotic estimate with the bounded oscillatory term removed in the black dashed line.To further make the point that our asymptotic estimate captures the characteristics of the exact coefficients, the plot on the right is zoomed on the coefficients near the peak.

Figure 1 . 3 :
Figure 1.3: Hermite functions: This is the same plot of the decay rate of coefficients as Fig. 1.2 but with a larger range.The plot shows that for varying values of ε = 10 −20 , 10 −30 , 10 −40 the number of coefficients which are greater than ε are 3691, 4584, and 5316 respectively.

Figure 2 . 2 :
Figure 2.2: Stretched Fourier functions: To understand the characteristics of our bound in (2.3), which requires us to choose λ ≥ |x 0 |, this plot shows the effect of choosing λ less than and greater than |x 0 |.The figure shows the plot of the logarithm of the absolute value of coefficients.For λ = 1, green dashed line, we see slow rate of decay as the majority of the 'mass' of the wave packet is outside the interval [-1,1].For λ = 9, we contain the majority of the wave packet in the interval [-9,9] and we see that the coefficients decay spectrally down beyond 10 −20 .

Figure 2 . 3 :
Figure 2.3: Stretched Fourier functions: For |a n (ω)| ≥ 10 −20, λ optimal is roughly 6.7861.The plot shows the effect of under-and overshooting.If we use a smaller λ, no matter how many coefficients we use, we would not arrive at our desired accuracy in reasonable time.Once λ > λ optimal , though, we reach ultimately excessive accuracy -but require more coefficients to reach the 10 −20 tolerance.

Figure 2 . 4 :
Figure 2.4: Stretched Fourier functions: Plot of (2.5) where the optimal λ is used, for fixed accuracy of ε = 10 −20 , and varying α and and x 0 .This demonstrates experimentally that the bound in Theorem 2.1 is representative of the behaviour of the coefficients in the crown of the sombrero, although the estimate of the beginning of the brim is universally pessimistic.Changes in the α parameter have minimal effect on the coefficients, whereas changes in x 0 have a significant effect.

Figure 3 . 1 :
Figure 3.1: Hermite functions: Increasing α makes our coefficients decay much slower.Just as in Fig. 1.2 we can see that the our asymptotic estimate approximates the exact coefficients very well.The figure shows how well the asymptotic estimate envelopes the exact coefficients.

Figure 3 . 2 :
Figure 3.2: Hermite functions: For a large shift in x 0 , the estimate does not do too well for small n.As n increases the estimate gets better. 3.3.

Figure 3 . 3 :
Figure 3.3: Hermite functions: Plot of (3.7) for different x 0 .We see that the decay rate of the exact hermite coefficients for α = 1/2 resembles a hump.
.14b) A is slowly varying for large n, Y and ν = √ 2n + 1, while B is the dominant contribution, the fast component.At this point, the reader can refer to section A.2 for the full derivation of the asymptotics.To summarise, we substitute terms w−φ (w−), φ(w−) and (3.11) into (3.14a) and (3.14b).

Figure 3 . 4 :
Figure 3.4: Hermite functions: Numerical experiments to determine the relationship between n and ω.The first figure is a plot of the position of the peak, n, as ω varies.The relationship follows n = ω 2 /2.The second figure is the experiment where we found the largest of coefficients such that |a n | > 10 −20 .This relationship follows n = 1 2 (ω + 18.5) 2 .The two experiments are consistent with √ n ∝ ω relationship.The last plot shows the peak position in black dotted line and largest n in dashed line for α = 2, x 0 = 2 and ω = 100.

Figure 3 . 5 :
Figure 3.5: Hermite functions: The dashed line represents our estimate without oscillations.This line envelopes the oscillations in the lighter colour, which further confirm that our estimate provides an accurate representation of the exact coefficients.The left column of figures are plots of log 10 |a n | against n, where the vertical scale goes from −20 to 0. The figures in the right column are close-up plots to show how well we are able to bound from above the exact coefficients.

. 15 )Fig. 3 .
Fig. 3.5 demonstrates how remarkably well our estimate models the size of the coefficients.The dominant term in the Hermite estimate with the bounded oscillations (the terms in the curly brackets in (3.15)) removed has the form an(ω) ∼ Cosc ω(2α + 1)

Figure 3 . 6 :
Figure 3.6: Hermite functions: For increasing α, we see that the rate of decay becomes slower.The dashed line indicate our asymptotic estimate as described by (3.16).Note how well our result envelopes the exact coefficients in the lighter plots.

Figure 3 . 7 :
Figure 3.7: Hermite functions: The figures on the left show the decay rate of the coefficients where log 10 |a n | > 10 −20 .In the top row we have α = 1/4, x 0 = 2 and ω = 30.In the first figure (on the left) we can see that our estimate starts to blow up around n = 600, this is expected as c * n,ω for these parameters occurs at n = 600.The plots in the right column are a closeup of the figures to the left.The plots in the second and third row, further show that our asymptotic estimates are accurate for varying α, x 0 and ω.Note, for these figures the c * n,ω occurs outside of the plotted region of log 10 |a n | > 10 −20 .
4.1.Why does multiplying by sin(x), an entire function, have such a large effect on the rate of convergence?How do we get exponential convergence?

Figure 4 . 2 :
Figure 4.2: MT functions: The plot compares the decay rate of the MT coefficients of the Gaussian (dotted black line) and the wave packet (solid green line).Note the slower decay rate of the wave packets relative to the decay rate of the Gaussian coefficients.

Figure 4 . 3 :
Figure 4.3: MT functions: Plot of the roots (4.9), (4.10) and (4.11) for varying n.Left plot: Plot of saddle points when x 0 = 0.All the saddle points are on the imaginary line when c n,ω < 1 4 , when c n,ω > 1 4, z 2 is the reflection of z 3 on the y axis.Right plot: Plot of saddle points when x 0 = 2. z 1 is shifted x 0 = 2 away from the imaginary axis and z 2 travels from the first quadrant to the fourth quadrant as n increases.

Figure 4 . 4 :
Figure 4.4:MT functions: Plot of (4.12) for given α, ω, x 0 and n.The ×s on the imaginary axis represent the pole at i/2 and zero at −i/2 for n ≥ 0. It can be seen in the top plots that the contour approaches an asymptote that is independent of x 0 .In the bottom plot, we see a zoomed in plot of the top left plot, the steepest descent contours pass through the pole at i/2 and zero at −i/2 for n ≥ 0.

Figure 4 . 5 :
Figure 4.5: MT functions: Schematic representation of the contour of integration Γ.The steepest descent path looks like a Gaussian centred around our root.

Figure A. 8 :
Figure A.8: Numerical experiments to determine the relationship between n and ω.The first figure is a plot of the position of the peak, n, as ω varies.The relationship follows n = 2ω 2 .The second figure is the experiment where we found the largest of coefficients such that |a n | > 10 −20 .This relationship follows n = 2.2ω 2 + 200.The two experiments are consistent with a √ n ∝ ω relationship.