Fast Computation of Orthogonal Systems with a Skew-symmetric Differentiation Matrix

Orthogonal systems in $\mathrm{L}_2(\mathbb{R})$, once implemented in spectral methods, enjoy a number of important advantages if their differentiation matrix is skew-symmetric and highly structured. Such systems, where the differentiation matrix is skew-symmetric, tridiagonal and irreducible, have been recently fully characterised. In this paper we go a step further, imposing the extra requirement of fast computation: specifically, that the first $N$ coefficients {of the expansion} can be computed to high accuracy in $\mathcal{O}(N\log_2N)$ operations. We consider two settings, one approximating a function $f$ directly in $(-\infty,\infty)$ and the other approximating $[f(x)+f(-x)]/2$ and $[f(x)-f(-x)]/2$ separately in $[0,\infty)$. In each setting we prove that there is a single family, parametrised by $\alpha,\beta>-1$, of orthogonal systems with a skew-symmetric, tridiagonal, irreducible differentiation matrix and whose coefficients can be computed as Jacobi polynomial coefficients of a modified function. The four special cases where $\alpha, \beta= \pm 1/2$ are of particular interest, since coefficients can be computed using fast sine and cosine transforms. Banded, Toeplitz-plus-Hankel multiplication operators are also possible for representing variable coefficients in a spectral method. In Fourier space these orthogonal systems are related to an apparently new generalisation of the Carlitz polynomials.


Introduction
This paper continues a project which we have commenced in [21], to investigate complete systems in L 2 (R) with a skew-symmetric differentiation matrix. Let us begin by explaining briefly the underlying concepts, motivating our work and revisiting in a very abbreviated manner the contents of [21].
The concept of a differentiation matrix originates in the theory of finite differences -essentially, this is the matrix taking a vector of function values to those of an approximation of the first derivative. An elementary example, is obtained from a second-order central difference approximation on a uniform grid with spacing ∆x > 0 -note that it is skew-symmetric. Skew-symmetry mimics the self-adjointness of the first derivative operator in the standard L 2 Hilbert space with either zero Dirichlet or periodic boundary conditions, and it confers a wide range of advantages on a numerical method. We refer to [15,17,18] for a wealth of specific examples: in essence, once a differentiation matrix is skew symmetric, it is often easy to prove stability for linear PDEs, as well as conservation of energy whenever it is mandated by the underlying equation. The matrix in (1.1) is skew symmetric, yet it is a sobering thought that this second-order approximation of the derivative is as good as it gets: no skew-symmetric finite-difference differentiation matrix on a uniform grid may exceed order 2 [17]. It is possible (but not easy) to obtain higher-order differentiation matrices of this kind carefully choosing specific non-uniform grids [15,16], but this is far from easy for high orders and, at any rate, finite differences are not the approach of choice in this paper.
An example that motivates much of our work is the linear Schrödinger equation in a semiclassical regime, (for simplicity, we focus here on the univariate case), typically given with periodic boundary conditions. The parameter ε > 0 is small and V is the potential energy [22]. Typically, (1.2) is solved by a spectral method (or spectral collocation) in space, followed by a combination of splittings (e.g. Strang splitting or Zassenhaus splitting) and, in the case of time-dependent potential, the Magnus expansion or a version thereof in time [3,4,19]. The definition of a differentiation matrix extends naturally into the realm of spectral methods. In any spectral method we expand the unknown in a (truncated) orthonormal basis {ϕ m } m∈Z + of the underlying L 2 Hilbert space, where the ϕ m s are suitably regular. The (infinite-dimensional) linear map D taking {ϕ m } m∈Z + into {ϕ m } m∈Z + is called a differentiation matrix. In other words, ϕ m = ∑ ∞ k=0 D m,k ϕ k , m ∈ Z + . The virtues of skew-symmetry remain undimmed in this setting and they immediately imply stability and energy conservation for (1.2).
Because of periodic boundary conditions, it is natural to use a Fourier basis for (1.2), leading to a skew-Hermitian (and diagonal) differentiation matrix. 1 To all intents and purposes, a skew-Hermitian matrix shares the benefits of a skewsymmetric one. Moreover, a Fourier basis has another critical virtue: we can approximate the first N expansion coefficients using FFT in just O (N log 2 N) operations. Therefore, there is little incentive to seek alternative bases.
However, the traditional setting -finite interval, periodic boundary conditionsof (1.2) and other dispersive equations of quantum mechanics is being increasingly challenged, since it is inappropriate for long-term integration (or even short-term integration in the presence of quantum scattering) because sooner or later, wave packets reach a boundary and periodicity leads to non-physical solutions. Modern applications of quantum mechanics, not least quantum control [30], are an important example of this behaviour. Arguably, the most natural setting for (1.2) is that of a Cauchy problem on the entire real line, with solutions confined to L 2 (R). Zero Dirichlet boundary conditions on an interval correspond to a closed quantum system which can be interpreted as a quantum system on the real line with u(x, 0) = 0 for x / ∈ (a, b) or where the potential V is sufficiently large at the endpoints that no tunnelling is possible [27].
In [21] we have developed a comprehensive theory of orthonormal bases of L 2 (R) with a skew-symmetric, tridiagonal, irreducible differentiation matrix where without loss of generality b m > 0, m ∈ Z. There exists a one-to-one relationship between absolutely continuous real Borel measures dµ with all their moments bounded, symmetric with respect to the origin and supported by the entire real line on the one hand, orthonormal bases Φ = {ϕ m } m∈Z + , complete in L 2 (R) and with a differentiation matrix of the form (1.3) on the other. The emphasis on a tridiagonal differentiation matrix is motivated by our wish to minimise the cost of eventual computations. It is standard to solve equations like (1.2) by splittings [2,4]. Thus, eventually we need to compute terms of the form e hD µu or e ihD 2 µu, and this becomes considerably cheaper once D is tridiagonal -a forthcoming paper of Celledoni and Iserles, for example, will describe an O (N log 2 N) algorithm to this end.
Specifically, let dµ(ξ ) = w(ξ )dξ be such Borel measure and {p m } m∈Z + the underlying monic orthogonal polynomials. Because of symmetry, they obey a threeterm recurrence relation of the form 1 We typically index a Fourier basis by integers, rather than non-negative integers, but this causes no difficulty in our framework.
where λ 0 = 0 and λ m > 0, m ∈ N. Set b m = λ m+1 , m ∈ Z + , and and g is a complex-valued function with even real part and odd imaginary part such that |g(ξ )| 2 = w(ξ ). It has been proved in [21] that {ϕ m } m∈Z + is indeed orthonormal and complete in L 2 (R) (essentially due to the Plancharel theorem) and has the differentiation matrix (1.3). On the other hand, given such a basis, using the Favard Theorem we can always associate it with a symmetric Borel measure supported by the real line: such measure is unique if the corresponding Hamburger moment problem is determinate [8, p. 73].
Realistic implementation of the aforementioned bases, however, requires that two computations can be performed explicitly, namely explicit formation of the Orthogonal Polynomial System (OPS) {p m } m∈Z + given dµ and the integrals (1.5), since otherwise we will not have the basis Φ in an explicit form. Here we are rapidly running against the limitations of current theory of orthogonal polynomials. An obvious and familiar example of Hermite functions where the H m are Hermite polynomials, aside, very few symmetric measures supported by the entire real line are known for which the underlying OPS (or even the recurrence coefficients λ m ) can be written explicitly: [21] list just the generalised Hermite polynomials and Carlitz polynomials. Indeed, the resolution of the Freud problem, precise asymptotic determination of the λ m s for w(ξ ) = e −ξ n +q(ξ ) , where n ≥ 2 is an even integer and q an even polynomial of degree ≤ n − 2, using the Riemann-Hilbert transform [10,13], is considered a recent triumph of the theory of orthogonal polynomials -and in our setting we need explicit values, not just asymptotics! Thus, our basic challenge in this paper is to search for systems Φ by alternative means. We seek systems that confer a tangible advantage when compared to the default of Hermite functions, and in this paper we find systems that allow for faster computation of expansion coefficients via Fast Cosine Transforms (FCT) or Fast Sine Transforms (FST). Given f ∈ L 2 (R) ∩ C ∞ (R), we wish to computê to sufficiently high accuracy in the least possible number of operations. We can identify three options. Firstly, it is possible to use the fast multipole method to compute expansion coefficients in O N log N + N log ε −1 operations to accuracy ε, provided that a table of nodes and weights of a suitable Gauss-Christoffel quadrature is available [12]. In the Hermite case, Gauss-Hermite quadrature is required and we note that its nodes and weights can be computed in O (N) operations [31]. For other bases of the form (1.5) quadrature rules may take up to O N 2 operations to compute.
Secondly, by the Fourier-transform on L 2 (R), the evaluation of thef m s reduces to an expansion in the polynomial basis {p m } m∈Z + . This, however, confers an advantage only if we can compute the polynomial coefficients in a fast manner. Of all fast and stable methodologies known to the present authors, this restricts us to dµ with bounded support, which has been ruled out earlier. The reason why we require supp µ = R is the completeness of the expansion: if, without loss of generality, supp µ = (−1, 1) then Φ is complete not in L 2 (R) but in the Paley-Wiener space PW [−1,1] (R) [21]. This is of no clear interest in the design of spectral methods for PDEs but might be relevant, for example, to signal processing.
In this paper we follow a third option. We seek a basis of ϕ m s so that the expansion coefficientsf m are equal to standard orthogonal polynomial coefficients of a modified function. If the standard polynomials are, for example, the Chebyshev polynomials (which a priori is not necessarily even possible) then the first N coefficients can be approximated in O (N log 2 N) operations using the FCT. In Section 2, this line of reasoning leads us to consider the model, where Θ ∈ L 2 (R), {q m } m∈Z + is an orthonormal OPS with respect to a measure W (t)dt supported in (−1, 1), and H is a differentiable, strictly monotonically increasing function mapping (−∞, ∞) onto (−1, 1). We prove in Section 2 that enforcing orthonormality on such a model means that the coefficients are equal to the coefficients of f (h(t))/Θ(h(t)) in the {q m } m∈Z + basis, where h is the inverse function of H. When we further require that our system has the differentiation matrix in (1.3), all options, up to an affine change of variables, are whittled down to is the mth Jacobi polynomial and (1.8) The first N expansion coefficients of a function with respect to Φ as in equation (1.7) can always be approximated in O N(log N) 2 operations by fast polynomial transform techniques [32]. The cases α, β = ±1 correspond to the Chebyshev polynomials of four different kinds (see [24, 18.3]), all of whose coefficients can be approximated in O (N log 2 N) operations via various flavours of FCT and FST. Furthermore, a given expansion f (x) = ∑ N−1 k=0 c k ϕ (α,β ) k (x) can be evaluated at a single point x ∈ R using Clenshaw's algorithm [9].
In Section 3 we turn our gaze to a more complicated model, Here Θ E and Θ O are given functions, the first even and the second odd, {r m } m∈Z + and {s m } m∈Z + are orthonormal OPSs with respect to the measures w E dt and w O dt, both supported by (−1, 1), and H maps strictly monotonically (0, ∞) to (−1, 1).
Neither w E nor w O need to be symmetric with respect to the origin. Advancing along similar lines to our analysis of (1.6), computation of coefficients corresponding to the system (1.7) can be reduced to the expansion in the two OPS, {r m } m∈Z + and {s m } m∈Z + . The exploration of (1.9) is a lengthier process.
Step after step we winnow the possibilities for H, Θ E , Θ O , w E and w O -at the end we are left with just a oneparameter family of such systems. We prove that (1.9) is consistent with orthonormality and with a skew-symmetric, tridiagonal, irreducible differentiation matrix if and only if w 1 2 for some α > −1 -in other words, the r m s and s m s are (normalised) Jacobi polynomials P (α,− 1 2 ) m and P (α, 1 2 ) m respectively. However, looking deeper, we demonstrate that these systems are identical to the full-range system (1.7) with α = β , implemented in a half-range mode.
Why do we need to discuss both the 'full-length' systems of Section 2 and the 'two half-length' systems of Section 3, in particular as the latter are nothing but an alternative implementation of the former? As things stand, they exhibit the same range of benefits: they are orthonormal and complete in L 2 (R), have skewsymmetric, tridiagonal differentiation matrices, while their expansion coefficients can be approximated in a fast and stable manner. We touch upon the advantages of either approach in Section 5. However, these are early days in the investigation of orthogonal systems with skew-symmetric differentiation matrices and their applications to spectral methods. There is great merit, we believe, in exploring their theory, potential and limitations in the broadest possible sense.
In Section 4 we give the functions g and measures dµ(ξ ) which appear in the formula (1.5) for these orthonormal bases, for all α, β > −1. They turn out to be the weights for a generalised version of the Carlitz polynomials (now with two parameters) discussed in [21], which have the explicit form, The case β = −α corresponds to the one-parameter family of Carlitz polynomials.
The original interest in Carlitz polynomials stems from the fact that their moments can be written in terms of Bernoulli polynomials. We are not aware of work on this two-parameter generalisation, nor can we anticipate its relevance to the theory of orthogonal polynomials on the real line.
In Section 5 we list some conclusions of our work, and flesh out the basis for applications of these functions to computing the Fourier transform and to a fast Olver-Townsend-type spectral method on the real line [26]. This is the place to mention [20], a companion paper to the current one, where we explore complex-valued orthonormal systems with skew-Hermitian tridiagonal differentiation matrices. Inter alia we demonstrate there the existence of another system with coefficients that can be computed in O (N log 2 N) operations, this time using the Fast Fourier Transform.
An alternative to the approach of this paper is to map a differential equation from L 2 (R), say, to a finite domain and solve the new equation there with a spectral method. This is often accomplished with conformal maps, providing a convenient means to study the speed of convergence [5,14,23,35]. This very natural and useful approach is genuinely different from the theme of this paper: the procedures of approximation and 'shrinkage' in general do not commute. As an example, consider a Cauchy problem for (1.2). A major feature of the linear Schrödinger equation is that the Euclidean norm of the solution remains constant as the time evolves, reflecting the quantum-mechanical interpretation of the solution: |u( ·,t)| 2 is a probability density. Mapping the equation into a finite interval, the L 2 norm corresponds to a fairly nonstandard Sobolev norm. While good enough for the proof of convergence by means of the Lax Equivalence Theorem, this makes matters quite difficult insofar as energy conservation is concerned.

The full-range setting
In this section we characterise all orthonormal systems Φ = {ϕ m } m∈Z + with a skew-symmetric, irreducible, tridiagonal differentiation matrix, with the additional constraint that the expansion coefficients of a function f ∈ L 2 (R) are equal to expansion coefficients for a weighted and mapped version of f , in an orthonormal polynomial basis on (−1, 1). Explicitly, for the differentiation matrix, we require where b −1 = 0 and, without loss of generality, b m > 0 for m ∈ Z + . For the expansion coefficients, we require where {q m } m∈Z are orthonormal polynomials with respect to the absolutely continuous measure W (t)dt on (−1, 1), and h, θ : (−1, 1) → R. We will make the following assumptions about θ , h and W . Later we deduce much more indeed about these functions as necessary consequences of our model.
• h maps onto the whole of R, is differentiable with h being a measurable, positive function.
• This implies the existence of an inverse function H : R → (−1, 1) which is differentiable with h (t)H (h(t)) = 1, therefore H is also positive and measurable.
of the weighted Cauchy-Schwarz inequality with weight h shows that this implies thatf m is finite for all m ∈ Z + . Changing variables to x = h(t) yields . For this to hold for all f ∈ L 2 (R), ϕ m must be of the form, The rest of the section is devoted to proving the following surprisingly simple result.
Theorem 2.1. All orthonormal systems of the form in equation (2.4) with a skewsymmetric, tridiagonal, irreducible differentiation matrix are, up to an affine change of variables, of the form The coefficients b m are given by (2.8).
We will discuss properties of these orthonormal systems in Section 4. In particular, we will show that they are in fact complete orthonormal bases for L 2 (R) by deriving the measure dµ(ξ ) and function g(ξ ) in the Fourier transform expression which these systems must satisfy as in equation (1.5) and [21].

Necessary conditions
Let us work out the necessary consequences of Φ being an orthonormal system in L 2 (R). For all n, m ∈ Z + , It follows at once that ). We will return to this later.
It is considerably more complicated to ensure the existence of a skew-symmetric, tridiagonal differentiation matrix. First note that by setting m = 0 in equation (2.4), we see that Θ is infinitely differentiable and in L 2 (R), because it is proportional to ϕ 0 (which is an infinitely differentiable function in L 2 (R), being the Fourier transform of a superalgebraically decaying L 2 (R) function by (1.5)).
Inserting equation (2.4) into equation (2.1), we obtain for all m ∈ Z + , , which can be substituted back into the equation for general m ∈ Z + to find Dividing through by H (x)Θ(x) and changing variables to t = H(x) gives us Here we used the fact discussed earlier that H (h(t))h (t) = 1. Now, the left-hand-side of equation (2.7) is a polynomial of degree exactly m−1, while the right-hand-side is h (t) times a polynomial of degree at most m+1. It follows from the case m = 1 that h (t) = 1/(at 2 + bt + c) for some real numbers a, b and c. Since h : (−1, 1) → R and h > 0, necessarily h (±1) = ∞ and it follows that, up to an affine change of variables, h(t) = arctanht and H(x) = tanh x.
Back to equation (2.6). Substituting in the necessary form of h we have, Using the fact that Θ (x) = b 0 q 1 (H(x))Θ(x)/q 0 (H(x)), we conclude after simple algebra that Writing 2b 0 q 1 (t)/q 0 + 2t = (1 + t)α + (1 − t)β for real numbers α and β , the solution to this ODE with W (0) = 1 without loss of generality, is It is none other than a Jacobi weight. We deduce that, if at all possible, the polynomials in the system (2.4) must be the (orthonormal) Jacobi polynomials, where g (α,β ) m has been defined in (1.8) and k m is an integer -for the time being, to allow for simpler algebra, we assume that k m = 0 but will change this later.
Substituting W (t) into equation (2.6) shows that

Sufficient conditions
All that remains is to check whether these systems actually satisfy the requirements set out at the start of the section. Firstly, the functions are in L 2 (R) if and only if α, β > −1 (which also ensures that W (t)dt is a finite measure). A quick change of variables will show that these functions are indeed orthonormal for all α, β > −1. For the final requirement on the differentiation matrix, we have, According to [24, 18.9.17-18] (1 − t 2 ) dP while the three-term recurrence relation for Jacobi polynomials is [28, p. 263]. In other words, We require that η m = −b m−1 and θ m = 0. However, before we look further into the above quantities, we note that necessarily b m < 0. Fortunately, this is an artefact of our choice of k m = 0: once we replace this with k m = m, we obtain b m > 0, m ∈ Z + , as required. This requires long, yet simple algebra, best performed using a symbolic algebra package. It follows from (1.8) that hence all the conditions are satisfied for

A half-range setting
In this section we take the model from the previous section a step further. Rather than requiring that our coefficients are equal to coefficients of a modified function in a single orthogonal polynomial basis, we consider two polynomial bases which are linked to the odd and even basis elements. In principle this provides a much richer range of possibilities over the more restrictive setting of the previous section, but we will discuss in subsection 3.3 that we have not actually gained anything in the end.
Specifically, suppose we have an orthonormal system Φ in L 2 (R) such that ϕ m is an even function if m is even and an odd function if m is odd; let us call this an even-odd system. Then we can re-express the expansion coefficients in the form In other words, we can 'translate' the setting from the real line to (0, ∞) and this forms the theme of this section. We suppose that the odd and even coefficients can each be expressed in a similar fashion to the previous section, but with separate expressions for each case, as follows. Let {r m } m∈Z + and {s m } m∈Z + be orthonormal polynomial systems with respect to the absolutely continuous measures w E (t)dt and w O (t)dt, respectively, on (−1, 1). We wish to find all orthonormal, even-odd systems Φ with a tridiagonal, skew-symmetric differentiation matrix, i.e.
and such that the coefficients are equal tô We will make the following assumptions about h, θ E , θ O , w E and w O . Just as in the previous section, we will deduce more from these basic assumptions and our model.
• h maps onto the whole of (0, ∞), is differentiable with h a measurable function, and h (t) > 0 for all t ∈ (−1, 1). • This implies the existence of an inverse function H : (0, ∞) → (−1, 1) which is differentiable with h (t)H (h(t)) = 1, which implies that H is also positive and measurable. • θ E is such that t → θ E (t) w E (t)/h (t) ∈ L ∞ (R), and θ O satisfies the analogous property. The motivation is exactly as in the previous section.
Changing variables to x = h(t) yieldŝ For this to hold for all f ∈ L 2 (R), we must necessarily have the 'half-range model', We extend H to the whole of R by setting H(−x) = H(x). Since ϕ 0 is an even function, we must have that Θ E is even function and likewise Θ O is an odd function. Note that ϕ m is an infinitely differentiable function for all m ∈ Z + because it is the Fourier transform of a superalgebraically decaying function by equation (1.5). Therefore, ϕ 1 (0) = 0 in order to have oddness. It follows that Θ O (0) = 0. Allow us to place one more assumption into the mix: Θ E (0) = 0. Otherwise all basis functions vanish at the origin, rendering it clearly unsuitable for the approximation of functions which are in general nonzero at the origin. The rest of the section is devoted to proving the following.
We show that these bases are equal to the bases in Theorem 2.1 with β = α, something which is far from obvious at first sight. All discussions of mathematical properties of these functions such as completeness in L 2 (R) can therefore be derived from properties of the functions in Theorem 2.1.

Necessary conditions
The first condition to explore is orthonormality. Since the ϕ m s have the parity of m on the real line, we have m, n ∈ Z + and need to check orthogonality only within each set. Changing variables, it follows from (3.4) that Since {r m } m∈Z + is an orthonormal set with respect to w E (t)dt, we must have and, by the same token, Note that, by the monotonicity of h, both weight functions are nonnegative, as required. These expressions can be inverted: changing back to x, (3.8) Observe that we need to be very careful in our choice of sign. As things stand, we allow for both options. Our next, considerably more challenging task is to ensure the existence of a differentiation matrix of the correct form. Recall that both Θ E and Θ O are smooth functions on R. Substituting (3.4) into (3.1), we have Setting m = 0 yields Substituting these back into (3.9) and changing variables back to t, we have, The way forward rests upon the observation that the left-hand sides of both (3.11) and (3.12) are (m − 1)-degree polynomials, and this places important constraints upon their right-hand sides. This is similar to the analysis of Section 2 yet considerably more complicated. Setting m = 1, taking products of equations (3.11) and (3.12) and performing some simple algebra, we obtain, Since r 1 , s 1 , and A 1 are polynomials of degree 1 and B 1 is a polynomial of degree 2, there exist constants a, b, c, γ 1 , γ 2 (which are real except that b and c may be complex conjugates) such that and Substituting the ratio of w O (t) and w E (t) for t ∈ (−1, 1) into equations (3.6) and (3.7), we find Since h(−1) = 0, Θ O (0) = 0 and Θ E (0) = 0 we deduce (without loss of generality) that c = −1. Since lim t→1 h(t) = +∞, we must have, integrating h (t), This is only possible if a = b = 1. Therefore, The formula for h (t) is readily integrated using the substitution tanh u = (1 + s)/2, giving Inverting, and ignoring a linear change of variables in x, we have whittled our way down to a single option, Now let us find Θ E (x). We know that Θ E (x) = Θ O (x)b 0 s 0 r 0 by equation (3.10) and Θ O (x) = γ 2 (1 + H(x))Θ E (x) by equation (3.14). Hence, for some α > −1, and, converting to the x variable, we have Equations (3.6) and (3.7) give us the weights, where without loss of generality w E (0) = w O (0) = 1.

Sufficient conditions
As things stand, we have identified one -and just one -one-parameter family of weights {w E , w O } for which we might be able to obtain an orthonormal system (3.4) with a tridiagonal skew-symmetric differentiation matrix: where we have used (3.8) (with a minus sign) to determine Θ E and Θ O . The question is, do the systems here actually satisfy our requirements for all α > −1? The following subsection answers this question in the affirmative by relating these functions to the full-range systems from Section 2.

The connection to full-range systems
Further investigation of these half-range systems leads one to the conclusion the half-range systems of Theorem 3.1 are a special case of full-range systems of Theorem 2.1 with β = α. The formulae look completely different, but as we will now show, they are identical.
For half-range functions while for full range (with β set to α), Now, by [24, 18.7.13], P (α,α) 2m (X) ∝ P (α,− 1 2 ) m (2X 2 − 1). Using this with X = tanh x, along with the identity tanh 2 x = 1−sech 2 x it follows readily that P (α,α) (1 − 2/ cosh 2 x). In addition, we know that (1 − tanh 2 x) α+1 2 ∝ sech 2 x. Combining the identities discussed in this paragraph, we arrive at the proportionality statement, The fact that the constants which depend on α but not x are uniquely determined to ensure that ϕ 2m has L 2 (R) norm equal to 1 proves that the two expressions for ϕ 2m given above are identical.
Likewise, in a half-range formalism, while the full-range expression is This time, we have by [24, 18.7.14] that P (α,α) 2m+1 (X) ∝ XP (α, 1 2 ) m (2X 2 − 1). Using this with X = tanh x, along with the identity tanh 2 x = 1 − sech 2 x it follows readily that P (α,α) (1 − 2/ cosh 2 x). In addition, we know that (1 − Combining the identities discussed in this paragraph, we arrive at the proportionality statement, The fact that the constants which depend on α but not x are uniquely determined to ensure that ϕ 2m+1 has L 2 (R) norm equal to 1 proves that the two expressions for ϕ 2m+1 given above are identical.
Theorem 2 therefore provides us with no new systems over the systems given by Theorem 1. This leads us to ask the question, could we have deduced Theorem 2 from Theorem 1, instead of conducting the full derivation of subsection 3.2?
We believe the answer is no, but in principle if one could show that the model in equation (3.4) necessarily reduces to the model in equation (2.4), then Theorem 2 would indeed follow directly from Theorem 1.

The tanh-Jacobi functions
In Section 2 we identified the following orthonormal bases Φ (α,β ) = {ϕ We call these function the tanh-Jacobi functions, for obvious reasons. These orthonormal bases have a tridiagonal, irreducible, skew-symmetric differentiation matrix as in equation In Section 3 we showed that in the special case where β = α, there is an alternative, equivalent expression which separates the functions into odd and even parts. For all α > −1, let for all m ∈ Z + , at least mathematically speaking. For computation of coefficients, this basis is different, as is discussed below.

Expansion coefficients
Using the change of variables t = tanh x, the expansion coefficients of a function f ∈ L 2 (R) in the orthonormal basis Φ (α,β ) can be expressed as which are the Jacobi coefficients of the modified function with a diagonal scaling. The regularity of F determines the convergence of the coefficients, since Jacobi polynomials bases have spectral convergence properties (this is an elementary consequence of integration by parts and derivative identity [24, 18.19.6]). For the half-range model, while the coefficients are equal to those given above with β = α, they take the following form. Let be the even and odd portions of f , respectively. Then, using the transformation The convergence of the coefficients is determined by the functions Exactly how properties of the function f itself determine the rate of convergence is a topic for another paper. It is clear from the outset that it is some combination of regularity and decay at infinity for the function f which will determine the regularity of F -if for example f only decays algebraically at infinity, then F will be unbounded in (−1, 1)! With regards to computation, the first N Jacobi coefficients of a given function can be approximated in O N(log N) 2 operations using fast polynomial transform techniques described in [32]. An efficient and straightforward Julia implementation exists in the software package, APPROXFUN [25].
The cases where α, β = ± 1 2 correspond to Chebyshev polynomials of various kinds. There are four kinds of Chebyshev polynomials: (1) Chebyshev polynomials of the first kind T m , orthogonal in (−1, 1) with the weight function (1 − t 2 ) − 1 2 ; (2) Chebyshev polynomials of the second kind U m , orthogonal in (−1, 1) with the weight function (1 − t 2 ) 1 2 ; (3) Chebyshev polynomials of the third kind V m , orthogonal in (−1, 1) with the weight function (1 − t) Because of this algorithmic advantage, and odd-even symmetry, the bases Φ (α,β ) with (α, β ) = − 1 2 , − 1 2 and (α, β ) = 1 2 , 1 2 , are our preferred choice of parameters. The Tanh-Chebyshev-T functions are given by The Tanh-Chebyshev-U functions are given by For the half-range model, equations (4.2) and (4.3) can also be converted into Sine and Cosine transforms. Specifically, with α = − 1 2 we obtain the combination of Chebyshev polynomials of the first and the fourth kind, T m and W m , which can be computed with FCT-I and FCT-II respectively -cf. [11] for different flavours of FCT and FST. Likewise, α = 1 2 results in a combination of U m and V m , computable with FST-I and FST-II respectively.

The Fourier transform representation and completeness
We will now prove that these bases are complete in L 2 (R). By [21,Thms 6,8] (see also equation (1.5)) there exists a complex-valued function g α,β (ξ ) with even real part and odd imaginary part such that where P = {p m } m∈Z + are the orthonormal polynomials with respect to dµ(ξ ) = w(ξ ) dξ = |g α,β (ξ )| 2 dξ . By [21,Thm. 9], the functions Φ (α,β ) are complete in L 2 (R) if w(ξ ) > 0 for all ξ ∈ R and polynomials are dense in the space L 2 (R, dµ(ξ )). Now, applying the Fourier transform to both sides of equation (4.4) and setting m = 0, we have We can normalise this function after we find an expression. Let us use the change of variables τ = e x to manipulate this integral into a more reasonable form.
The change of variables σ = τ 2 transforms this into dσ This is one of the standard integral formulae for Euler's Beta function [24, 5.12.3], meaning that Using the identity B(a, b) = Γ(a)Γ(b)/Γ(a + b) for all a, b ∈ C, we deduce the existence of a real constant C α,β such that Note that g α,β is complex-valued in general, but has an even real part and odd imaginary part, which implies that ϕ α,β m (x) is real-valued for all x ∈ R. Barnes's Beta Integral [24, 5.13.3] can be used to compute the constant C α,β which makes the measure dµ α,β (ξ ) = |g α,β (ξ )| 2 dξ have mass equal to one.
We have thus proved the following theorem.
Theorem 4.1. The Fourier-space measure associated with {ϕ There is little hope of simplifying the expression for dµ α,β from Theorem 4.1 for general α, β > −1 except in some special cases.
The polynomials associated with this measure are known as the Carlitz polynomials (after mapping them from a line in the complex plane to the real line) [7]. They were mentioned in [21], but now we have the full picture of their relationship to orthogonal systems with a tridiagonal, skew-symmetric differentiation matrix. The case α = β yields Lemma 4.2. In the case α = β ≥ 0 being integers we have Proof. In the case n = 0 (4.6) is confirmed by direct computation. Otherwise we use the standard recurrence formula Γ(z + 1) = zΓ(z). (4.5) implies that g n,n (ξ ) C n,n = (n − 1) 2 + ξ 2 4 g n−2,n−2 (ξ ) C n−2,n−2 , n ≥ 2, and (4.6) follows by simple induction.
Regarding completeness of Φ (α,β ) in L 2 (R), as mentioned above, g α,β being nonzero for all ξ ∈ R (because the Γ function has no roots in the complex plane), all that remains is the question of density of polynomials in the space L 2 (R, dµ). It would be sufficient that w(ξ ) have exponential decay as ξ → ±∞ [1, pp. 45,86,86].
To show that these measures have exponential decay, we use the asymptotic formula [24, 5.11.9], This implies, as required. An intriguing fact is that (4.5) can be alternatively derived from a little-known formula due to Ramanujan [29], namely Taking a = (α + 1)/2, a trivial change of variable takes g α,α to the correct ϕ 0 and the proof follows by inverting the Fourier transform. Going in the reverse direction, we can obtain a generalisation of Ramanujan's formula, to the Fourier transform of Γ(a + iξ )Γ(b − iξ ) for a, b > 0. Except for Carlitz polynomials and their immediate generalisations, orthogonal polynomials associated with measures dµ α,β do not appear to have been studied in the literature. They might be an interesting object for further study, being examples of measures supported by the entire real line, yet distinct from the more familiar Freud-type measures.

Conclusions
We set ourselves a goal in this paper: to identify and characterise orthonormal systems, complete in L 2 (R) and with a skew-symmetric, tridiagonal, irreducible differentiation matrix whose expansion coefficients can be computed rapidly. In particular we are interested in the computation of the first N expansion coefficients in O (N log 2 N) operations, utilising familiar transforms, e.g. FFT, FCT or FST.
In Section 2 we introduced the Tanh-Jacobi functions, of which the four cases where α, β = ± 1 2 have expansion coefficients which can be computed using FCTs or FSTs. Subsequently, in Section 3 we described two kinds of half-range expansions (i.e., treating the even and the odd part of a function separately) which can be computed rapidly with either FCT or FST: the first based on a combination of Chebyshev polynomials of the first and the fourth kind, the second on such polynomials of the second and third kind.
Mathematically, the two approaches are identical when α = β , although their computation is somewhat different. Which is preferable? As things stand, there is no clear answer (and things are complicated by the availability of yet another approach of this kind, using skew-Hermitian differentiation matrices, which is described in [20]). The full-range approximation has the virtue of simplicity, hence of easier implementation. A possible advantage of a half-range approximation is more subtle. Once D approximates the first derivative, D 2 approximates the second one and, provided D is skew symmetric, D 2 is negative semi-definite -this is in line with the Laplace operator being negative semi-definite and is vital for stability.
As a square of a skew-symmetric, tridiagonal matrix, D 2 neatly separates even and odd functions. Specifically, let E ⊕ O = L 2 (R) ∩ C 2 (R) be a representation of square-integrable, twice differentiable functions on the line as a direct sum of even functions E and odd functions O. A derivative takes E to O and vice versa, hence a second derivative is invariant in both E and O. There is thus a virtue, at least once both first and second derivatives are present, to work separately in E and O, as is the case with half-range approximations.
Is simplicity preferable to even-odd separation? Are there additional considerations at play? By this stage it is impossible to provide a definitive answer. The purpose of the paper is to present a range of new results that improve our knowledge of approximation on the real line in the context of spectral methods.
The outlook for spectral methods on the real line using Tanh-Chebyshev-T functions appears promising. Consider the basic first order differential operator, where a is a bounded function such that a rapidly convergent expansion of the form a(x) = ∑ ∞ m=0 a mTm (tanh x) is possible 2 . In coefficient space for the Tanh-Chebyshev-T functions, the operator L becomes the infinite-dimensional matrix D + A , where D is the skew-symmetric tridiagonal differentiation matrix (1.3) with b m = 1 2 m + 1 2 , and a 0 a 1 a 2 a 3 · · · a 1 a 0 a 1 a 2 . . . a 2 a 1 a 0 a 1 . . .
The matrix is Toeplitz-plus-Hankel and, if a is sufficiently regular then the matrix is effectively banded, because if a m = 0 for m > M for some integer M then the resulting operator has bandwidth M. This implies that an Olver-Townsend type approach for infinite-dimensional QR solution is possible in principle. An rth order differential operator with variable coefficients whose expansions have a maximum of M terms yields a matrix with bandwidth of at most r + M. Furthermore, the resulting matrices respect certain symmetries that the operator L may have, something which is not true of the original Olver-Townsend ultraspherical spectral method. For example, if L is self-adjoint on L 2 (R) then L is a self-adjoint operator on 2 .
One final, credible, but as of yet unexplored, application of this work is for the computation of the Fourier transform of certain functions on the real line. Let f ∈ L 2 (R) have a rapidly convergent expansion in one of the Tanh-Jacobi bases, f (x) = ∑ m=0 c m ϕ (α,β ) m (x). An approximation f N (x) = ∑ N m=0 c N m ϕ (α,β ) m (x) computed using either fast polynomial transforms or the FCT/FST as discussed in Section 4. By the identity in equation (1.5) and [21], the Fourier transform of f N is an expansion in the generalised Carlitz polynomials weighted by g α,β (see equation (4.2)), By Theorem 6 of [21], these generalised Carlitz polynomials, which are orthonormal with respect to dµ α,β (ξ ) in Theorem 3, satisfy where b m is given in equation (2.8). Clenshaw's algorithm can be used to evaluate this expansion at a single point ξ ∈ R, using purely the coefficients {b m } m∈Z + [9]. This is similar to the work of Weber in which one computes a series expansion with orthogonal rational functions whose Fourier transforms are Laguerre functions, utilising the Fast Fourier Transform [33]. A major issue that we have not pursued in this paper is the speed of convergence of different orthonormal bases in L 2 (R), not least Tanh-Jacobi bases. While approximation theory of analytic functions in a finite interval by means of an analytic orthonormal basis is well known, this is not the case on the real line. The issue (which persists under a map from R to a finite interval) is that this theory requires both the underlying function and the orthonormal basis to be analytic in a larger ellipse surrounding a finite interval: it is the size and shape of that ellipse that determines the exponential speed of convergence. This becomes problematic on the real line. Some orthonormal bases (e.g. Hermite functions) have an essential singularity at infinity (viewed as the North Pole of the Riemann sphere), as does, for example, the tanh map. Even when, like the Malmquist-Takenaka system [20], a basis is analytic in a strip surrounding R, our problems are not over because most analytic functions of interest are likely to have an essential singularity at infinity. As a striking example, while the Malmquist-Takenaka expansion coefficients of 1/(1 + x 2 ) (which is analytic ina strip about R decay exponentially, as 3 −|n| , the speed of decay for sin x/(1 + x 2 ) is O |n| −5/4 , barely better than linear [34]! While important results have been published by Boyd [6], much more needs be done to understand approximation theory on the real line. Appendix: