A family of orthogonal rational functions and other orthogonal systems with a skew-Hermitian differentiation matrix

In this paper we explore orthogonal systems in $\mathrm{L}_2(\mathbb{R})$ which give rise to a skew-Hermitian, tridiagonal differentiation matrix. Surprisingly, allowing the differentiation matrix to be complex leads to a particular family of rational orthogonal functions with favourable properties: they form an orthonormal basis for $\mathrm{L}_2(\mathbb{R})$, have a simple explicit formulae as rational functions, can be manipulated easily and the expansion coefficients are equal to classical Fourier coefficients of a modified function, hence can be calculated rapidly. We show that this family of functions is essentially the only orthonormal basis possessing a differentiation matrix of the above form and whose coefficients are equal to classical Fourier coefficients of a modified function though a monotone, differentiable change of variables. Examples of other orthogonal bases with skew-Hermitian, tridiagonal differentiation matrices are discussed as well.


Introduction
The motivation for this paper is the numerical solution of time-dependent partial differential equations on the real line.It continues an ongoing project of the present authors, begun in (Iserles & Webb 2019b), which studied orthonormal systems Φ = {ϕ n } n∈Z in L 2 (R) which satisfy the differential-difference relation, for some real, nonzero numbers {b n } n∈Z where b −1 = 0.In other words, the differentiation matrix of Φ is skew-symmetric, tridiagonal and irreducible.The virtues of skew symmetry in this context are elaborated in (Hairer & Iserles 2016, Iserles 2016) and (Iserles & Webb 2019b) -essentially, once Φ has this feature, spectral methods based upon it typically allow for a simple proof of numerical stability and for the conservation of energy whenever the latter is warranted by the underlying PDE.The importance of tridiagonality is clear, since tridiagonal matrices lend themselves to simpler and cheaper numerical algebra.
In this paper we generalise (1), allowing for a skew-Hermitian differentiation matrix.In other words, we consider systems Φ of complex-valued functions such that where {b n } n∈Z+ ⊂ C and {c n } n∈Z+ ⊂ R.
While the substantive theory underlying the characterisation of orthonormal systems in L 2 (R) with skew-Hermitian, tridiagonal, irreducible differentiation matrices is a fairly straightforward extension of (Iserles & Webb 2019b), its ramifications are new and, we believe, important.In Section 2 we establish this theory, characterising Φ as Fourier transforms of weighted orthogonal polynomials with respect to some absolutely-continuous Borel measure dµ.This connection is reminiscent of (Iserles & Webb 2019b) but an important difference is that dµ need not be symmetric with respect to the origin: this affords us an opportunity to consider substantially greater set of candidate measures.
An important issue is that, while the correspondence with Borel measures guarantees orthogonality and the satisfaction of (2), it does not guarantee completeness.In general, once dµ is determinate and supported by the interval (a, b), completeness is assured in the Paley-Wiener space PW (a,b) (R).
So far, the material of this paper represents a fairly obvious generalisation of (Iserles & Webb 2019b).Furthermore, the operation of differentiation for functions on the real line is defined without venturing into the complex plane.Indeed, it is legitimate to challenge why we should allow our differentiation matrices to contain complex numbers.After all, if skew-Hermitian framework is so similar to the (simpler!)skew-symmetric one, why bother?The only possible justification is were (2) to confer an advantage (in particular, from the standpoint of computational mathematics) in comparison with (1).This challenge is answered in Section 3 , where we consider sets Φ associated with generalised Laguerre polynomials, where (a, b) = (0, ∞).We show that a simple tweak to our setting assures the completeness of these Fourier-Laguerre functions, which need be indexed over Z, rather than Z + .
A remarkable property of the MT system (3) is that the computation of the expansion coefficients can be reduced, by an easy change of variables, to a standard Fourier integral.Therefore the evaluation of f−N , . . ., fN−1 can be accomplished with the Fast Fourier Transform (FFT) in O (N log 2 N ) operations: this has been already recognised, e.g. in (Weideman 1994).In Section 4 we characterise all systems Φ, indexed over Z, which tick all of the following boxes: • They are orthonormal and complete in L 2 (R), • They have a skew-Hermitian, tridiagonal differentiation matrix, and • Their expansion coefficients f−N , . . ., fN−1 can be approximated with a discrete Fourier transform by a single change of variables, and hence computed in O (N log 2 N ) operations with fast Fourier transform.
Adding rigorous but reasonable assumptions to these requirements, we prove that, modulo a simple generalisation, the MT system is the only system which bears all three.
We wish to draw attention to (Iserles & Webb 2019a), a companion paper to this one.While operating there within the original framework of (Iserles & Webb 2019b) -skew-symmetry rather than skew-Hermicity -we seek therein to characterise orthonormal systems in L 2 (R) whose first N coefficients can be computed in O (N log 2 N ) operations by fast expansion in orthogonal polynomials.We identify there a number of such systems, all of which can be computed by a mixture of fast cosine and fast sine transforms.Such systems are direct competitors to the Malmquist-Takenaka system, discussed in this paper.
2 Orthogonal systems with a skew-Hermitian differentiation matrix 2.1 Skew-Hermite differentiation matrices and Fourier transforms The subject matter of this section is the determination of verifiable conditions equivalent to the existence of a skew-Hermitian, tridiagonal, irreducible differentiation matrix (2) for a system Φ = {ϕ n } n∈Z+ which is orthonormal in L 2 (R).
Theorem 2.1 (Fourier characterisation for Φ).The set has a skew-Hermitian, tridiagonal, irreducible differentiation matrix (2) if and only if where P = {p n } n∈Z+ is an orthonormal polynomial system on the real line with respect to a non-atomic probability measure dµ with all finite moments1 , g is a square-integrable function which decays superalgebraically fast as |ξ| → ∞, and {θ n } n∈Z+ is a sequence of numbers in [0, 2π).Furthermore, P , g, and {θ n } n∈Z+ are uniquely determined by ϕ 0 , {c n } n∈Z+ , and Remark 2.2.This theorem is a straightforward generalisation of (Iserles & Webb 2019b, Thm. 6), which shows the same result but for real, irreducible skewsymmetric differentiation matrices.The difference is that (2) is replaced by (1), dµ must be even, g must have even real part and odd imaginary part, and θ n is chosen so that e iθn = (−i) n .We will prove sufficiency because it is elementary but enlightening, and leave necessity and uniqueness for the reader to prove by modifying the proof in (Iserles & Webb 2019b).That part of the proof depends on Favard's theorem and properties of the Fourier transform, and we wish to avoid it for the sake of brevity.
Proof.Suppose that ϕ n are given by the equation (1).Then by (Gautschi 2004, Thm. 1.29) there exist real numbers {δ n } n∈Z+ and positive numbers {β n } n∈Z+ such that where β −1 = 0 by convention.3Differentiating under the integral sign and using the above three-term recurrence, we obtain Theorems 2.3 and 2.4 are proved in (Iserles & Webb 2019b) for the real case, as in equation ( 1).The proofs require minimal modification for them to apply to the complex case, as in equation (2).
Theorem 2.3 (Orthogonal systems).Let Φ = {ϕ n } n∈Z+ satisfy the requirements of Theorem 2.1.Then Φ is orthogonal in L 2 (R) if and only if P is orthogonal with respect to the measure |g(ξ)| 2 dξ.Furthermore, whenever Φ is orthogonal, the functions ϕ n / g 2 are orthonormal.
Theorem 2.4 (Orthogonal bases for a Paley-Wiener space).Let Φ = {ϕ n } n∈Z+ satisfy the requirements of Theorem 2.3 with a measure dµ such that polynomials are dense in L 2 (R; dµ).Then Φ forms an orthogonal basis for the Paley-Wiener space PW Ω (R), where Ω is the support of dµ.
The key corollary of Theorem 2.4 is that for a basis Φ satisfying the requirements of Theorem 2.3 to be complete in L 2 (R), it is necessary that the polynomial basis P is orthogonal with respect to a measure which is supported on the whole real line.

Symmetries and the canonical form
Let Φ = {ϕ n } n∈Z+ have a tridiagonal skew-Hermitian differentiation matrix as in equation ( 2).Then the system Φ = { φn } n∈Z+ defined by where ω, A, B, C, κ n ∈ R and A, B = 0, also satisfies equation (2).We can show this directly as follows.
If the differentiation matrix is irreducible then these symmetries permit a unique choice of κ 0 , κ 1 , . . .ensuring that b n is a positive real number for each n ∈ Z + .This corresponds to modifying the choice of θ n in Theorem 2.1 so that e iθn = i n .It is therefore possible for any given g and P to have a canonical choice of Φ, which satisfies b n > 0, by taking We can also produce a unique canonical orthonormal system from an absolutely continuous measure dµ(ξ) = w(ξ)dξ on the real line, where w(ξ) decays superalgebraically fast as |ξ| → ∞.Specifically, the functions form an orthonormal system in L 2 (R) with a tridiagonal, irreducible skew-Hermitian differentiation matrix with a positive superdiagonal.The system

Computing Φ
We proved in ( functions that obey (1) obeys the relation where Our setting lends itself to similar representation, which follows from (2) by induction.
Lemma 2.5.The functions Φ consistent with (2) satisfy the relation where Like (5), the formula ( 6) is often helpful in the calculation of ϕ 1 , ϕ 2 , . . .once ϕ 0 is known.The obvious idea is to compute explicitly the derivatives of ϕ 0 and form their linear combination (6), but equally useful is a generalisation of an approach originating in (Iserles & Webb 2019b).Thus, Fourier-transforming (6), On the other hand, Fourier transforming (4), we have Our first conclusion is that φ0 (ξ) = |w(ξ)| 1/2 /p 0 .Moreover, comparing the two displayed equations, The polynomials p n are often known explicitly.In that case it is helpful to rewrite (6) in a more explicit form.

An example
The next section is concerned with the substantive example of a system Φ with a skew-Hermitian differentiation matrix that originates in the Fourier setting once we use a Laguerre measure.What, though, about other examples?Once we turn our head to generating explicit examples of orthonormal systems in the spirit of this paper and of (Iserles & Webb 2019b), we are faced with a problem: all steps in subsections 2.1-3 must be generated explicitly.Thus, we must choose an absolutely continuous measure for which the recurrence coefficients in (2) are known explicitly, compute explicitly {p n } n∈Z+ and either e ixξ dξ and its derivatives, subsequently forming (8) and manipulating it further into a userfriendly form, or Either course of action is restricted by the limitations on our knowledge of explicit fomulae of orthogonal polynomials for absolutely continuous measures (thereby excluding, for example, Charlier and Lommel polynomials, as well as the Askey-Wilson hierarchy).Thus Hermite polynomials and their generalisations (Iserles & Webb 2019b), Jacobi and Konoplev polynomials (Iserles & Webb 2019b), Carlitz polynomials (Iserles & Webb 2019a) and, in the next section, Laguerre polynomials.Herewith we present another example which, albeit of no apparent practical use, by its very simplicity helps to illustrate our narrative.Let α ∈ R and consider dµ(ξ) = e −(ξ−α) 2 dξ, a shifted Hermite measure.The underlying orthonormal set consists of -we deduce that b n = (n + 1)/2 and c n ≡ α in (2).Moreover, 'twisted' Hermite functions.It is trivial to confirm that they satisfy (2) or derive them directly from (8).

Connections to chromatic expansions
Theorem 1 characterises all systems Φ in L 2 (R) satisfying equation ( 2).These systems depend on a family of orthonormal polynomials on the real line with associated measure dµ and a function g on the real line.Theorems 2 and 3 focussed on the special case in which dµ(ξ) = |g(ξ)| 2 dξ.This special case yields all systems which are orthonormal in the inner standard product, which turn out to be complete in a Paley-Wiener space.
An anonymous referee has made the authors aware of a considerable amount of work devoted to the special case in which dµ(ξ) = g(ξ)dξ, whose systems generate so-called chromatic expansions (Ignjatovic 2007, Zayed 2009, Zayed 2011, Zayed 2014).These systems have some remarkable properties for application to signal processing which we summarise here whilst making connections to the present work.
Given a sequence of orthonormal polynomials {P n } n∈Z+ with respect to a finite Borel measure dµ on the real line, define the function and the operators for all n ∈ Z + acting on C ∞ (R).The chromatic expansion of a smooth function f is the formal series which converges uniformly on the real line if, for example, µ is such that ψ is analytic in a strip around the real axis, f is analytic in this strip too, and the sequence of coefficients {K n [f ](0)} n∈Z+ is in 2 (Zayed 2014).
The connection to the present work is as follows.For all n ∈ Z + , let ϕ n = K n [ψ] be the elements of the system Φ.Then Φ is of the form described in Theorem 1 with g(ξ)dξ = dµ(ξ).
This advantages for signal processing are twofold: • The expansion coefficients are given explicitly and depend locally on the function f centred around the point 0. The expansions can be made local to points other than 0 as in (Ignjatovic 2007).
• The functions Φ are bandlimited, with Fourier transforms supported precisely on the support of µ.
While the basis Φ with ϕ n = K n [ψ] is not orthonormal in the standard inner product on L 2 (R), under some mild assumptions on µ, it is possible to show that Φ is orthonormal with respect to the inner product and is complete in a space of analytic functions on the real line for which the induced norm is finite (Ignjatovic 2007).
3 The Fourier-Laguerre basis

A general expression
A skew-Hermite setting allows an important generalisation of the narrative of (Iserles & Webb 2019b), namely to Borel measures in the Fourier space which are not symmetric.The most obvious instance is the Laguerre measure dµ(ξ) = χ (0,∞) (ξ)ξ α e −ξ dξ, where α > −1.The corresponding orthogonal polynomials are the (generalised) Laguerre polynomials where (z is the Pochhammer symbol and 1 F 1 is a confluent hypergeometric function (Rainville 1960, p. 200).The Laguerre polynomials obey the recurrence relation n−1 (ξ).First, however, we need to recast them in a form suitable to the analysis of Section 2 -specifically, we need to renormalise them so that they are orthonormal and so that the coefficient of ξ n in p n is positive.Since 1960, p. 206) and the sign of ξ n in ( 1) is (−1) n , we set We deduce after simple algebra that 2) and ( 2).(b n = β n because the latter is real and positive.) To compute Φ we note that, letting τ = ( 1 2 − ix)ξ, (4) yields Moreover, and substitution in (8) gives The identity, Olver, Lozier, Boisvert & Clark 2010, 15.8.7), implies that we have where Π n is a polynomial of degree n.Using the substitution x = 1 2 tan θ 2 for θ ∈ (−π, π), which implies (1 + 2ix)/(1 − 2ix) = e iθ , the orthonormality of the basis Φ can be seen to imply that {Π (α) n } n∈Z+ are in fact orthogonal polynomials on the unit circle (OPUC) with respect to the weight To be clear, this means that for all n, m ∈ Z + , n (e iθ )Π (α)  m (e iθ ) cos α θ 2 dθ = δ n,m .
These polynomials are related to the Szegő-Askey polynomials (Olver et al. 2010, 18.33.13),{φ n due to the symmetry of the weight of orthogonality (Szegő 1975, p. 295), (Olver et al. 2010, 18.33.14).Specifically, for some real constants {A n , B n , C n , D n } n∈Z+ .It is therefore possible to express the functions Φ in terms of Jacobi polynomials; this is something we will not pursue here, but could be of interest for further research.In what follows we will restrict ourselves to the case α = 0, which is extremely simple.We are not aware if this connection between the general Laguerre polynomials and Szegő-Askey polynomials (and hence Jacobi polynomials) via the Fourier transform has been acknowledged before in the literature.

The Malmquist-Takenaka system
The expression (2) comes into its own once we let α = 0, namely consider the 'simple' Laguerre polynomials L n .Now W (θ) ≡ 1 and so Π (α) n (z) = z n .We have b n = n + 1, c n = 2n + 1 and The factor of (−1) n which might appear to have been added here comes from the identity 2 F 1 −n, 1; 1; z = (1 − z) n with z = 2/(1 − 2ix).Alternatively, we may apply a formula for the Laplace transform of Laguerre polynomials at an appropriate point in the complex plane (Olver et al. 2010, 18.17.34), to obtain By Theorem 2.4, these functions are dense in the Paley-Wiener space PW [0,∞) (R).To obtain a basis for the whole of L 2 (R), we must add to this a basis for PW (−∞,0] (R).The obvious way to do so is to consider the same functions as above, but for the orthogonal polynomials with respect to χ (−∞,0] (ξ)e ξ dξ, which are precisely L n (−ξ), n ∈ Z + .Using the Laplace transform again, this leads to the functions Letting ϕ n = φ−n−1 , n ≤ −1, we obtain the Malmquist-Takenaka system (3).
As a matter of historical record, Malmquist (1926) and Takenaka (1926) considered a more general system of the form The nature of the questions they have asked was different -essentially, they proved that the above system is a basis (which need not be orthogonal) of H 2 , the Hardy space of complex analytic functions in the open unit disc.In our case the θ k ≡ 2i are all the same and outside the unit circle, yet it seems fair (and consistent with, say, (Pap & Schipp 2015)) to call (3) a Malmquist-Takenaka system.
Fig. 1 displays the real and imaginary parts of few Malmquist-Takenaka functions.
Let us dwell briefly on the properties of (3).
Figure 1: Real and imaginary parts (blue lines with black '+'s and green lines with black '×'s, respectively) of the MT functions ϕ n for n = −3, . . ., 2. The envelope of ± 2 π(1+4x 2 ) is also plotted as a dashed line.
While (3) follows from our construction, it can be also proved directly from (3): • The MT system obeys a host of identities that make it amenable for implementation in spectral methods.The following were identified by Christov, 1982) and the following is apparently new, In particular, (4) implies that allowing for an easy multiplication of expansions in the MT basis.

Expansion coefficients
Arguably the most remarkable feature of the MT system is that expansion coefficients can be computed very rapidly indeed.Thus, let f ∈ L 2 (R).Then We do not dwell here on speed of convergence except for brief comments in subsection 3.4 -this is a nontrivial issue and, while general answer is not available, there is wealth of relevant material in (Weideman 1994).Our concern is with efficient algorithms for the evaluation of fn for −N ≤ n ≤ − 1.
The key observation is that and the term on the right is of unit modulus.We thus change variables in other words x = 1 2 tan θ 2 and, in the new variable We deduce that fn a Fourier integral.Two immediate consequences follow.Firstly, the convergence of the coefficients as |n| → ∞ is dictated by the smoothness of the modified function Secondly, provided F is analytic, ( 6) can be approximated to exponential accuracy by a Discrete Fourier Transform5 and this, in turn, can be computed rapidly with Fast Fourier Transform (FFT): the first N coefficients require O (N log 2 N ) operations.
Proposition 3.1 (Fast approximation of coefficients).The truncated MT system {ϕ n } N −1 n=−N is orthonormal with respect to the discrete inner product, where and can be computed simultaneously in O (N log 2 N ) operations using the FFT.
Proof.Let k, be integers satisfying .
If k = then this is clearly equal to 1. Otherwise, using equation ( 5), we see that, Summing the geometric series, since θ j − θ j−1 = π/N we have This proves that {ϕ n } N −1 n=−N forms an orthonormal basis with respect to the inner product • , • N .It follows that f, h = f, h N for all f and h in the span of {ϕ n } N −1 n=−N .Inserting h(x) = ϕ n (x) into the expression for the discrete inner product and then using equation ( 5) yields (7).

Speed of convergence
for some ρ > 1 if and only if the function t → (1 − 2it)f (t) can be analytically continued to the set where C is the Riemann sphere consisting of the complex plane and the point at infinity, and D r (a) is the disc with centre a ∈ C and radius r > 0, with Proof.See (Weideman 1995) and (Boyd 1987).
As was noted by Weideman, for exponential convergence we require f to be analytic at infinity, which is of meagre practical use.An example for a function f of this kind is Since f is a meromorphic function with singularities at ±i, we obtain exponential decay with ρ = 3 -this is evident from the explicit expansion whose proof we leave to the reader.This is demonstrated in Fig. 2, where we display the errors f (x) − N n=0 fn ϕ n (x) for N = 10, 20, 30 and 40.Compare this with Fig. 3, where we have displayed identical information for an expansion in Hermite functions.Evidently, MT errors decay at an exponential speed, while the error for Hermite functions decreases excruciatingly slowly as N increases.
Meromorphic functions, however, are hardly at the top of the agenda when it comes to spectral methods.In particular, in the case of dispersive hyperbolic equations we are interested in wave packets -strongly localised functions, exhibiting double-exponential decay away from an envelope within which they oscillate rapidly.An example (with fairly mild oscillation) is the function Since f has an essential singularity at infinity, there is no ρ > 1 so that (8) holds -in other words, we cannot expect exponentially-fast convergence.We report errors for MT and Hermite functions in Figs 4 and 5 respectively for N = 10, 40, 70 and 100: definitely, the convergence of MT slows down (part of the reason is also the oscillation) but it still is superior to Hermite functions.The general rate of decay of the error (equivalently, the rate of decay of | f|n| | for n 1 for analytic functions and the MT system) is unknown, although (Weideman 1995) reports interesting partial information, which we display in Table 1 (taken from (Weideman 1995)).The rate of decay does not seem to follow simple rules.For some functions the rate of decay is spectral (faster than a reciprocal of any polynomial), yet sub-exponential.For other functions it is polynomial (and fairly slow).Fig. 6 exhibits MT errors for f (x) = sin x/(1+x 2 ) and N = 20, 40, 60, 80: evidently it is in line with Table 1.It is fascinating that such a seemingly minor change to (9) has such far-reaching impact on the rate of convergence.This definitely calls for further insight.
A future paper will address the rate of approximation of wave packets by both the MT basis and other approximation schemes.

Characterisation of mapped and weighted Fourier bases
The most pleasing feature of the MT basis is that the coefficients can be expressed as Fourier coefficients of a modified function.They can then be approximated using the Fast Fourier Transform.Are there other orthogonal systems like this?Let us consider all orthonormal systems Φ = {ϕ n } n∈Z in L 2 (R) with a tridiagonal skew-Hermitian differentiation matrix such that for all f ∈ L 2 (R), the coefficients are equal to the classical Fourier coefficients of k(θ)f (h(θ)), −π < θ < π, for some functions k and h (with a possible diagonal scaling by {γ n } n∈Z ).Specifically, we consider the ansatz We assume that h : (−π, π) → R is a differentiable function which is strictly Table 1: The rate of decay of the coefficients fn in an MT approximation of different functions.
increasing and onto, whose derivative is a measurable function.This implies the existence of a differentiable, strictly increasing inverse function H : R → (−π, π).
The chain rule implies h (θ)H (h(θ)) ≡ 1 (so that H is also a measurable function).The function k is assumed to be a complex-valued L 2 (−π, π) function (which makes the integral in (1) well defined).The constants γ n are complex numbers.We assume nothing more about k, h and γ n in this section (but deduce considerably more).
Making the change of variables x = h(θ) yields, γ n e −inH(x) k(H(x))H (x)f (x) dx. (2) For this to hold for all f ∈ L 2 (R), we must have where How does this fit in with the MT basis?In the special case of Malmquist-Takenaka we have We prove the following theorem which characterises the Malmquist-Takenaka system as (essentially) the only one of the kind described by equation (3).
Theorem 4.1.All systems Φ = {ϕ n } n∈Z of the form (3), such that Proof.Let us derive some necessary consequences of orthonormality of Φ by applying the change of variables x = h(θ) to the inner product.
Orthogonality implies that the function θ → h (θ)|K(h(θ))| 2 is orthogonal to θ → e ikθ for all k ∈ Z \ {0}.It is therefore a constant function.This constant is positive since h is strictly increasing and K is not identically zero.Normality of the basis implies that |γ n | 2 = 2πh (θ)|K(h(θ))| 2 −1 , which is a constant independent of n.We can absorb this constant into K and assume that Since ϕ 0 (x) = γ 0 K(x) and ϕ 0 is infinitely differentiable (because it is proportional to the inverse Fourier transform of a superalgebraically decaying function g), we deduce that K must be infinitely differentiable.The relationship H (x) = 2π|K(x)| 2 therefore implies that H is infinitely differentiable; in particular H (x) = 4π K (x)K(x) .Furthermore, there exists an infinitely differentiable function κ : R → R such that Let us derive more necessary consequences by taking into account the tridiagonal skew-Hermitian differentiation matrix.For all n ∈ Z, 2H (x) K(x), so dividing through by K(x)γ n ie inH(x) leads to where β n = −ib n γ n+1 /γ n (here we use the fact that γ −1 n = γ n ).Without loss of generality, we can assume that β n ∈ R for all n because the symmetries discussed in subsection 2.2 allow us to choose {γ n } n∈Z (because they are all of the form e iκn for real numbers {κ n } n∈Z ).
Since κ and H are real-valued functions and c n and β n are real for all n ∈ Z, equating real and imaginary parts yields It follows that β n − β n−1 is a constant which is independent of n, so we can write a = 2(β n − β n−1 ) for some real constant a and equation ( 9) becomes which, after integrating with respect to x, becomes for some real constant b.Since H maps R onto (−π, π) in a strictly increasing manner, by the monotone convergence theorem we must have H (±∞) = 0. Since H (±∞) = a cos(±π) + b, we must have a = b.Therefore, using the formula cos θ + 1 = 2 cos 2 (θ/2), we obtain, Integrating with respect to x, we get tan for some real constant c.Hence there exist real constants A and B such that Note that necessarily A = 0.All that remains is to determine K(x), which can be done by determining κ(x).Taking n = 0 in equation ( 8) gives us The antiderivative of cos(2 arctan(Ax + B)) is 2A −1 arctan(2x) − x, so there exist real constants ω, a and b such that κ(x) = ωx + 2a arctan(Ax + B) + b.
Whence we deduce that where γ n = e ib (−1) n+a .Letting δ = a − 1 2 shows that the system Φ must necessarily be of the form in equation (4).To complete the proof we must turn to the question of sufficiency.A derivation exactly as in subsection 3.2 but with n replaced by n + δ verifies the explicit form of the coefficients (5) for the case γ n = (−i) n , λ = i/2 and ω = 0.The symmetry considerations in subsection 2.2 show that the other values of γ n , λ and ω yield orthonormal systems with a tridiagonal skew-Hermitian differentiation matrix too.

Concluding remarks
The subject matter of this paper is the theory of complex-valued orthonormal systems in L 2 (R) with a tridiagonal, skew-Hermitian differentiation matrix.On the face of it, this is a fairly straightforward generalisation of the work of (Iserles & Webb 2019b).Yet, the more general setting confers important advantages.In particular, it leads in a natural manner to the Malmquist-Takenaka system.The latter is an orthonormal system of rational functions, which we have obtained from Laguerre polynomials through the agency of the Fourier transform.The MT system has a number of advantages over, say, Hermite functions, which render it into a natural candidate for spectral methods for the discretization of differential equations on the real line.It allows for an easy calculus, because MT expansions can be straightforwardly multiplied.Most importantly, the calculation of the first N expansion coefficients can be accomplished, using FFT, in O (N log 2 N ) operations.Moreover, the MT system is essentially unique in having the latter feature.
The FFT, however, is not the only route toward 'fast' computation of coefficients in the context of orthonormal systems on L 2 (R) with skew-Hermitian or skew-symmetric differentiation matrices.In (Iserles & Webb 2019a) we characterised all such real systems (thus, with a skew-symmetric differentiation matrix) whose coefficients can be computed with either Fast Cosine Transform, Fast Sine Transform or a combination of the two, again incurring an O (N log 2 N ) cost.We prove there that there exist exactly four systems of this kind.
The connections laid out in Section 3 between the Fourier-Laguerre functions and the Szegő-Askey polynomials (and hence Jacobi polynomials via the Delsarte-Genin transformation), are suggestive of a possible generalisation of Theorem 4.1 on the characterisation of the MT basis.It may be possible to characterise all systems which are orthonormal, have a tridiagonal skew-Hermitian differentiation matrix, and which are of the form ϕ n (x) = Θ(x)Π n e iH(x) , where Θ ∈ L 2 (R), H maps the real line onto (−π, π), and {Π n } n∈Z+ is a system of orthogonal polynomials on the unit circle.The expansion coefficients for a function in such a basis are equal to expansion coefficients of a mapped and weighted function in the orthogonal polynomial basis {Π n } n∈Z+ .The Fourier-Laguerre bases, in particular the MT basis, are certainly within this class of functions, but one can ask if there are more.
From a practical point of view, it is worth noticing that while the MT basis elements decay like |x| −1 as x → ±∞, the Fourier-Laguerre functions decay like |x| −1−α/2 where α > −1 is the parameter in the generalised Laguerre polynomial.For the approximation of functions with a known asymptotic decay rate it may be advantageous to use a basis with the same decay rate.
The jury is out on which is the 'best' orthonormal L 2 (R) system with a skew-Hermitian (or skew-symmetric) tridiagonal differentiation matrix and whose first N coefficients can be computed in O (N log 2 N ) operations.While some considerations have been highlighted in (Iserles & Webb 2019a), probably the most important factor is the speed of convergence.Approximation theory in L 2 (R) is poorly understood and much remains to be done to single out optimal orthonormal systems for different types of functions.Partial results, e.g. in (Ganzburg 2018, Weideman 1994), indicate that the speed of convergence of such systems is a fairly delicate issue.