text:
\begin{document}

\title{An extended Filon--Clenshaw--Curtis method for high-frequency wave scattering problems in two dimensions}

\author{
G. Maierhofer\thanks{Email: g.maierhofer@maths.cam.ac.uk}, A. Iserles and N. Peake}
\affil{DAMTP, Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge, CB3 0WA, UK}

\maketitle

\begin{abstract}

{We study the efficient approximation of integrals involving Hankel functions of the first kind which arise in wave scattering problems on straight or convex polygonal boundaries. Filon methods have proved to be an effective way to approximate many types of highly oscillatory integrals, however finding such methods for integrals that involve non-linear oscillators and frequency-dependent singularities is subject to a significant amount of ongoing research. In this work, we demonstrate how Filon methods can be constructed for a class of integrals involving a Hankel function of the first kind. These methods allow the numerical approximation of the integral at uniform cost even when the frequency $\omega$ is large. In constructing these Filon methods we also provide a stable algorithm for computing the Chebyshev moments of the integral based on duality to spectral methods applied to a version of Bessel's equation. Our design for this algorithm has significant potential for further generalisations that would allow Filon methods to be constructed for a wide range of integrals involving special functions.
These new extended Filon methods combine many favourable properties, including robustness in regard to the regularity of the integrand and fast approximation for large frequencies. As a consequence, they are of specific relevance to applications in wave scattering, and we show how they may be used in practice to assemble collocation matrices for wavelet-based collocation methods and for hybrid oscillatory approximation spaces in high-frequency wave scattering problems on convex polygonal shapes.}
\end{abstract}
\keywords{highly oscillatory integrals; numerical integration; wave scattering.}

\section{Introduction}
\label{sec:introduction}
The efficient approximation of highly oscillatory integrals is a crucial step in many numerical simulations of physical systems, including in high-frequency wave scattering. Although efficient methods for the computation of highly oscillatory integrals have been discovered as early as the first half of the twentieth century by Louis Napoleon George \cite{filon1930iii}, and thorough research over the past decades has lead to an increase in efficiency and applicability of such methods \citep{deano2017computing}, integrals arising from numerical wave scattering still pose significant difficulty, due to the nature of the Green's function of the Helmholtz operator: The Green's function involves special functions that are not only oscillatory but also singular, with frequency-dependent singularities and, in some cases, stationary points. In this work we will consider two types of integrals arising from wave scattering and in particular from the use of hybrid approximation spaces in this context \citep{chandler2012numerical}:
\begin{align}\label{eqn:form_of_integrals_alpha_zero}
I_{\omega,\beta}^{(1)}[f]&=\int_{0}^{1}f(x)H_0^{(1)}(\omega x)e^{i\omega\beta x}dx,
\end{align}
where $\beta\in\mathbb{R}$, and
\begin{align}
\label{eqn:form_of_integrals_alpha_non_zero}
I_{\omega,\alpha,\beta}^{(2)}[f]&=\int_{-1}^{1}f(x)H_0^{(1)}\left(\omega \sqrt{x^2-2\alpha\beta+\alpha^2}\right)e^{i\omega\beta x}dx,
\end{align}
where $\alpha>  0, |\beta|<1,\omega>0$. In both cases $f$ is a sufficiently differentiable amplitude function and $H_0^{(1)}$ is the Hankel function of the first kind of order zero. These two integrals arise from a collocation method on screens and convex polygonal wave scatterers when the approximation spaces proposed by \cite{chandler-wilde2007} and \cite{hewett2014} are used.

One class of highly oscillatory quadrature methods are so-called Filon methods, named after Louis Napoleon George Filon. These have been developed further and analysed by \cite{iserles2004,iserles2005} and \cite{iserlesnorsett2004} who demonstrated their favourable asymptotic properties in the high-frequency regime. With renewed interest in the study of these methods, significant progress was made over the last two decades, amongst many by Sheehan \cite{olver_2007} who described a moment-free version of the Filon method that is applicable to algebraic singularities and stationary points, and in the non-singular setting by \cite{Dominguez2011} who provided a stable and efficient way for the moment computation for Filon--Clenshaw--Curtis methods and error estimates which are explicit in both the number of interior approximation points and in the asymptotic order of the method. Improvements to the asymptotic order in Filon methods were introduced by \cite{iserlesnorsett2005efficient} by including information about the derivative values of the amplitude function $f$ and resulted in the development of the extended Filon method by \cite{Gao2017a,Gao2017b} which can achieve any desired asymptotic order.

In the context of singular oscillatory integrals \cite{piessens1983modified} introduced a Filon-type method that permits the efficient computation of integrals involving Bessel functions of linear arguments. This work was further extended by \cite{Zhenhua2016} to allow for Hankel kernels, which resulted in the treatment of integrals that effectively correspond to \eqref{eqn:form_of_integrals_alpha_zero} in the special case when $\beta=0$. A related Filon method that can be applied to logarithmic singularities of the form $\log|x-\alpha|$ was developed by \cite{dominguez2014}. In parallel, \cite{keller1999,keller2007method} described a fairly general framework for computing indefinite integrals with singular and oscillatory kernels arising from a range of special functions using Clenshaw--Curtis points. He also provided a way of constructing finite difference equations for the Chebyshev moments based on a differential equation satisfied by the kernel function. However, while for certain cases Keller was able to provide conditions under which stable algorithms for the moment computation can be found from this recurrence, in general this stable computation of the moments has to be studied on a case-by-case basis. More recently \cite{dominguez2013filon} constructed a composite (graded) version of the Filon method that can be applied to arbitrary algebraic and logarithmic singularities. While this method is already significantly better than traditional quadrature and, in fact, the application of non-singular Filon quadrature, we will see in section \ref{sec:numerical_examples_comparison_to_previous_work} that our direct Filon-based work further improves on this approach through a direct exploitation of the asymptotic properties of the integrals under consideration.

The work on applying these Filon methods in the context of efficient hybrid numerical methods for wave scattering is more limited. \citep{dominguez2013filon} has been successfully applied to hybrid methods in wave scattering by \cite{chandler2012numerical}, \cite{kim2012} and \cite{parolin2015} using the advantages of robustness and flexibility of the Filon approach which requires only very weak assumptions on the type of singularity and the regularity of $f$. Notable in the context of the efficient solution of integral equations based on Filon methods is also work by \cite{Shuhuang2011} where a Filon-type method for Bessel transforms is applied in the solution of Volterra integral equations.

An alternative approach to computing highly oscillatory integrals is numerical steepest descent which was introduced by \cite{huybrechs2006evaluation}. Numerical steepest descent has recently been applied by \cite{gibbs2019fast} to wave scattering problems on multiple screens and can also be applied to convex polygonal scatterers (see \cite{software_hnabemlab,software_pathfinder}). These results serve as a comparison to our methods in section \ref{sec:numerical_examples_scattering_on_screen}. While numerical steepest descent in general outperforms Filon methods in terms of asymptotic order and computational expense, Filon methods have the advantage that they are more robust towards changes in the behaviour of $f$, and in particular require no analyticity assumptions on the integrand as opposed to the requirement of at least local analyticity for numerical steepest descent. Furthermore, as long as moments can be computed accurately, Filon methods permit uniform error estimates that are explicit both in the number of interior interpolation points and in the frequency of the integral (see \citep{Dominguez2011} and \citep{Gao2017b}).

In the present work we analyse the asymptotic behaviour of \eqref{eqn:form_of_integrals_alpha_zero} and \eqref{eqn:form_of_integrals_alpha_non_zero} carefully to construct a direct Filon-type method that can achieve any desired asymptotic order while remaining convergent also for small values of the frequency $\omega$. The advantages of this method are significantly smaller errors for large frequencies in comparison to previous applications of Filon methods to this setting and removing the requirement for grading. In addition, this method requires very weak assumptions on $f$, even continuously differentiable $f$ will yield an asymptotic decay in the error, and in contrast to numerical steepest descent no region of analyticity needs to be assumed. We extend results by \cite{Gao2017a} to show how the interpolation problem at Clenshaw--Curtis points can be solved efficiently in this setting and provide a stable and efficient way of computing the Chebyshev moments using a duality relation to a spectral method. In this context our theorem \ref{thm:recurrence_by_duality} provides an extension of \cite[Lemma 2.4]{keller2007method} to arbitrary basis functions. Finally, we provide an application of our new Filon methods to a collocation method for high-frequency wave scattering problems on screens and convex polygons, and show how with this quadrature the collocation matrix can be assembled at frequency independent cost, resulting in a very efficient method for these scattering problems.

We begin this paper in section \ref{sec:extended_filon_method} by recalling the main features of the extended Filon method and by motivating the necessity of constructing a new Filon method in the context of frequency-dependent singularities such as \eqref{eqn:form_of_integrals_alpha_zero}. This is followed by our extension of the Filon paradigm to the current setting and finally an outline of our new Filon methods for the integrals \eqref{eqn:form_of_integrals_alpha_zero} and \eqref{eqn:form_of_integrals_alpha_non_zero} in section \ref{sec:extended_filon_method_for_Hankel_kernels}. This method requires a way of computing Chebyshev moments, and we provide a robust and efficient algorithm that achieves this for both types of integrals in section \ref{sec:moment_computation_by_duality}. The construction is based on a duality argument and we provide asymptotic and numerical justification for the stability in section \ref{sec:algorithm_for_stable_moment_computation}. We follow this with extensive numerical examples in section \ref{sec:numerical_examples_filon_methods} providing comparison of our new methods to previous work and two applications of the quadrature to high-frequency wave scattering on a screen and convex polygon. Our results are summarised and an outlook towards future research directions is provided in the concluding remarks in section \ref{sec:conclusions}.
\section{The extended Filon method}\label{sec:extended_filon_method}
In order to accurately approximate integrals of the form \eqref{eqn:form_of_integrals_alpha_zero} and \eqref{eqn:form_of_integrals_alpha_non_zero} at uniform cost for all frequencies $\omega\geq0$ we seek to construct a version of the extended Filon method (which for non-singular integrals was introduced by \cite{Gao2017a,Gao2017b}) that captures the important features of the asymptotic behaviour of the integrals. Gao and Iserles achieved this in the non-singular case by observing that for integrals of the form
\begin{align*}
I[f;g]=\int_{-1}^{1}f(x)e^{i\omega g(x)}dx,\quad \text{where\ }g'(x)\neq 0, -1<x<1,
\end{align*}
the asymptotic expansion of $I[f;g]$ for large $\omega$ depends only on the values
\begin{align*}
\mathcal{S}=\{f^{(j)}(\pm 1):j=0,1,\dots\}.
\end{align*}
\subsection{What goes wrong with classical singular Filon?}\label{sec:why_singular_Filon_cannot_be_applied}
As mentioned in the introduction there are Filon methods that deal with singularities in $f(x)$ and stationary points in $g(x)$, for instance \citep{olver_2007}. The Hankel function $H_0^{(1)}(z)=h_0(z)\exp(iz)$ where $h_0(z)$ is non-oscillatory in the sense of lemma \ref{lem:phase_extraction_in_H_0}, and $h_0(z)$ has a logarithmic singularity at the origin. This essentially suggests to use an approximation basis of the form
\begin{align*}
\mathcal{B}=\left\{x^j\right\}_{j=1}^n\cup \left\{x^j\log x\right\}_{j=1}^n
\end{align*}
to interpolate $f$ and to approximate
\begin{align*}
I^{(1)}_{\omega,\beta}[f]\approx \int_{0}^{1}(p(x)+q(x)\log x)e^{i\omega (\beta+1)x}dx,
\end{align*}
where $p,q$ are some appropriate choice of polynomials. The pitfall with this approach is that it is not possible to find an approximation for $h_0(\omega x)$ in the space $\mathcal{B}$ that remains uniformly accurate as $\omega \rightarrow \infty$. This motivates our search for a more direct Filon method that allows us to incorporate the singular and oscillatory behaviour of $H^{(1)}_0$ at the same time.
\subsection{Filon paradigm for the Hankel kernel}
\label{sec:filon_paradigm_for_hankel_kernel}
Indeed we are able to determine the dominant contributions to $I^{(1)}_{\omega,\beta}$ and $I^{(2)}_{\omega,\alpha,\beta}$ when $\omega\rightarrow\infty$ and exploit this knowledge to construct a well-behaved quadrature method. In order to do so we briefly recall the following property of the Hankel function of the first kind of order 0, $H_0^{(1)}$:
\begin{lemma}[Phase extraction of $H_0^{(1)}$, see Lemma 4.6 in \cite{chandler2012numerical}]\label{lem:phase_extraction_in_H_0}
	Let $h_0(z):=\exp(-iz)H_0^{(1)}(z)$, then for each $n\geq 0$ there is a constant $C_n$ such that
	\begin{align*}
	\left|\frac{\D^n}{\D z^n}h_0(z)\right|\leq C_n\begin{cases}
	\max\left\{1+\log(1/z),z^{-n}\right\},&z\in(0,1]\\
	z^{-(n+1/2)},&z\in[1,\infty).
	\end{cases}
	\end{align*}
\end{lemma}
The above lemma \ref{lem:phase_extraction_in_H_0} essentially tells us that $h_0(z)$ changes slowly in its argument, i.e. it does not oscillate.
\begin{proposition}[Filon paradigm for \eqref{eqn:form_of_integrals_alpha_zero}]\label{prop:filon_paradigm_for_alpha_zero}
	Let $\beta\neq -1$, and suppose $f\in C^{k+2}[0,1]$ has $k$ derivatives that vanish at the endpoints $x=0,1$, then
	\begin{align*}
	I_{\omega,\beta}^{(1)}[f]=\Ob_\delta\left(\omega^{-(k+2)+\delta}\right),\text{\ for\ any\ }\delta>0,
	\end{align*}
	where, by $\Ob_\delta$ we mean that the asymptotic constant depends on $\delta$.
\end{proposition}
\begin{proof}
	Choose $\epsilon=\epsilon(\omega)>0$ such that $\epsilon\rightarrow 0\ \wedge\  \omega\epsilon\rightarrow\infty$ as $\omega\rightarrow\infty$. Then
	\begin{align*}
	I_{\omega,\beta}^{(1)}[f]=\int_{0}^\epsilon H^{(1)}_0(\omega x)f(x)e^{i\omega \beta x}dx+\int_{\epsilon}^1H^{(1)}_0(\omega x)f(x)e^{i\omega\beta x}dx
	\end{align*}
	To bound the first integral note that by our assumptions on $f$ there is a constant $\tilde{C}>0$ such that $|f(x)|\leq \tilde{C} x^{k+1}$ in a neighbourhood of $0$, therefore for $\epsilon$ sufficiently small:
	\begin{align*}
	\left|\int_{0}^\epsilon H^{(1)}_0f(x)e^{i\omega\beta x}dx\right|\lesssim \int_{0}^{\epsilon}\left| H^{(1)}_0(\omega x)\right|x^{k+1}dx.
	\end{align*}
	By Lemma \ref{lem:phase_extraction_in_H_0} for $n=0$ there is a constant $C_0$ such that for all $x>0$:
	\begin{align*}
	\left|H_0^{(1)}(\omega x)\right|&\leq C_0\begin{cases}
	\max\left\{1+\log(1/(\omega x)),1\right\},&\omega x\in(0,1]\\
	(\omega x)^{-1/2},&\omega x\in[1,\infty)
	\end{cases}\\
	&\leq 2C_0(\omega x)^{-1/2}
	\end{align*}
	Since this bound applies uniformly in $\omega x$ we thus have:
	\begin{align*}
	\left|\int_{0}^\epsilon H^{(1)}_0(\omega x)f(x)e^{i\omega g(x)}dx\right|\lesssim\omega^{-\frac{1}{2}} \int_{0}^{\epsilon}x^{k+\frac{1}{2}}dx\lesssim\omega^{-\frac{1}{2}}\epsilon^{k+\frac{3}{2}}.
	\end{align*}
	Moreover by integration by parts we have (noting that $h_0$ is non-singular on $(0,1]$ and that $f^{(j)}(1)=0,j=0,\dots,k$):
	\begin{align*}
	\int_{\epsilon}^1H^{(1)}_0(\omega x)f(x)e^{i\omega \beta x}dx&=\int_{\epsilon}^1h_0(\omega x)f(x)e^{i\omega(\beta+1) x}dx\\
	&=\sum_{j=0}^k\left(\frac{-1}{i\omega (\beta+1)}\right)^{j+1}\left[\frac{\D^j}{\D x^j}\left(h_0(\omega x)f(x)\right)\right]_{x=\epsilon}\\
	&\quad+\left(\frac{-1}{i\omega (\beta+1)}\right)^{k+2}\left[\frac{\D^{k+1}}{\D x^{k+1}}\left(h_0(\omega x)f(x)\right)\right]_{x=\epsilon}^{1}\\
	&\quad+\left(\frac{-1}{i\omega (\beta+1)}\right)^{k+2}\int_{\epsilon}^1\left[\frac{\D^{k+2}}{\D x^{k+2}}\left(h_0(\omega x)f(x)\right)\right]e^{i\omega(\beta+1) x}dx
	\end{align*}
	Now observe that since the first $k$ derivatives of $f$ vanish at $0$ we have
	\begin{align*}
	\left|f^{(l)}(\epsilon)\right|\lesssim \epsilon^{k+1-l}, l=0,\dots,k
	\end{align*}
	as $\epsilon\rightarrow 0$. Further observe that since $\epsilon\omega\rightarrow\infty$ as $\omega\rightarrow\infty$, we have by lemma \ref{lem:phase_extraction_in_H_0} for $j=0,\dots, k+1$:
	\begin{align*}
	\left|\left[\frac{\D^j}{\D x^j}\left(h_0(\omega x)f(x)\right)\right]_{x=\epsilon}\right|&\leq \sum_{l=0}^j \binom{j}{l}\omega^l \underbrace{\left|\left[\frac{\D^lh_0}{\D x^l}\right]_{x=\omega\epsilon}\right|}_{\lesssim \omega^{-l-1/2}\epsilon^{-l-1/2}}\quad \underbrace{\left|f^{(j-l)}(\epsilon)\right|}_{\lesssim \epsilon^{k+1-(j-l)}}\\
	&\lesssim \omega^{-1/2}\epsilon^{k-j+1/2},\quad\text{as\ }\epsilon\rightarrow 0, \omega\epsilon\rightarrow\infty.
	\end{align*}
	Similarly we find:
	\begin{align*}
	\left|\left[\frac{\D^{k+1}}{\D x^{k+1}}\left(h_0(\omega x)f(x)\right)\right]_{x=1}\right|&=\left|\left(h_0(\omega )f^{k+1}(1)\right)\right|\lesssim \omega^{-1/2},
	\end{align*}
	and
	\begin{align*}
	\left|\frac{\D^{k+2}}{\D x^{k+2}}\left(h_0(\omega x)f(x)\right)\right|\lesssim \omega^{-1/2}x^{-1/2}\left|f^{(k+2)}(x)\right|+\sum_{l=1}^{k+2}\omega^{-1/2}x^{-l-1/2}\left|f^{(k+2-l)}(x)\right|\lesssim x^{-1},
	\end{align*}
	uniformly in $\epsilon\leq x\leq 1$ as $\epsilon\rightarrow 0, \omega\epsilon\rightarrow\infty$. Therefore
	\begin{align*}
	\int_{\epsilon}^1\left[\frac{\D^{k+2}}{\D x^{k+2}}\left(h_0(\omega x)f(x)\right)\right]e^{i\omega(\beta+1) x}dx\lesssim \int_{\epsilon}^{1}x^{-1}dx\lesssim |\log\epsilon|.
	\end{align*}
	Thus in total we find
	\begin{align*}
	\left|I_{\omega,\beta}^{(1)}[f]\right|\lesssim \sum_{j=0}^{k+2}\omega^{-\frac{1}{2}-j}\epsilon^{k+\frac{3}{2}-j}+\omega^{-k-2-1/2}+\omega^{-k-2}\log|\epsilon|\lesssim_{\delta}\omega^{-(k+2)+\delta},\text{\ for\ any\ }\delta>0.
	\end{align*}
	Where the final line follows by taking $\epsilon=\omega^{-1+\tilde{\delta}}$, for a suitable choice of $\tilde{\delta}>0$.
\end{proof}

Similarly we can establish the Filon paradigm for $I^{(2)}_{\omega,\alpha,\beta}$, which is non-singular but has a stationary point at the origin. We follow the motivation and derivation in \cite[pp. 8-9]{deano2017computing} with a little extra care that treats the benign frequency dependence of $h_0$:
\begin{proposition}[Filon paradigm for \eqref{eqn:form_of_integrals_alpha_non_zero}]\label{prop:filon_paradigm_for_alpha_non_zero}
	Let $\beta\neq 1$, and suppose $f\in C^{k+2}[0,1]$ has $k$ derivatives that vanish at each of the points $x=-1,0,1$, then
	\begin{align*}
	I_{\omega,\alpha,\beta}^{(2)}[f]=\Ob\left(\omega^{-(k+2)}\right),\text{\ as\ }\omega \rightarrow \infty.
	\end{align*}
\end{proposition}
\begin{proof}
	Integrating by parts $n+1$ times we find
	\begin{align*}
	I_{\omega,\alpha,\beta}^{(2)}[f]&=\int_{-1}^{1}f(x)h_0\left(\omega \sqrt{x^2-2\alpha\beta+\alpha^2}\right)\exp\left[i\omega \underbrace{\left(\sqrt{x^2-2\alpha\beta+\alpha^2}+\beta x\right)}_{=:g(x)}\right]dx\\
	&=f_0(0;\omega)\int_{-1}^{1}\exp\left[i\omega g(x)\right]dx+\int_{-1}^{1}(f_0(x;\omega)-f_0(0;\omega))\exp\left[i\omega g(x)\right]dx\\
	&=f_0(0;\omega)\int_{-1}^{1}\exp\left[i\omega g(x)\right]dx+\frac{1}{i\omega}\left[\frac{f_0(1;\omega)-f_0(0;\omega)}{g'(1)}e^{i\omega g(1)}-\frac{f_0(-1;\omega)-f_0(0;\omega)}{g'(-1)}e^{i\omega g(-1)}\right]\\
	&\quad-\frac{1}{i\omega}\int_{-1}^{1}f_1(x;\omega)\exp\left[i\omega g(x)\right]dx
	\end{align*}
	\begin{align*}	 &\hspace{-1cm}=\cdots=(-1)\left(\frac{-1}{i\omega}\right)^{n+1}\left(f_{n+1}(0;\omega)\int_{-1}^{1}\exp\left[i\omega g(x)\right]dx\right.\\
	&\hspace{-1cm}\quad+\frac{1}{i\omega}\left[\frac{f_{n+1}(1;\omega)-f_{n+1}(0;\omega)}{g'(1)}e^{i\omega g(1)}-\frac{f_{n+1}(-1;\omega)-f_{n+1}(0;\omega)}{g'(-1)}e^{i\omega g(-1)}\right]\\
	&\hspace{-1cm}\quad\left.-\frac{1}{i\omega}\int_{-1}^{1}f_{n+2}(x;\omega)\exp\left[i\omega g(x)\right]dx\right)
	\end{align*}
	where we use the notation that is also introduced by \cite{deano2017computing}:
	\begin{align*}
	f_0(x;\omega)=f(x)h_0\left(\omega \sqrt{x^2-2\alpha\beta+\alpha^2}\right), \quad\quad f_n(x;\omega)=\frac{d}{dx}\frac{f_{k-1}(x;\omega)-f_{k-1}(0;\omega)}{g'(x)},\quad k\geq 1.
	\end{align*}
	Now we by lemma \ref{lem:phase_extraction_in_H_0} we have the estimate
	\begin{align*}
	\left|f_0(x;\omega)\right|\leq C_{f} \omega^{-1/2}
	\end{align*}
	for some constant depending on (the uniform norm of) $f$. Consequently, by the quotient rule and lemma \ref{lem:phase_extraction_in_H_0}, there is a constant $D_{f}>0$ such that
	\begin{align*}
	\left|f_{j}(x;\omega)\right|\leq D_f \omega^{-1/2},\quad 0\leq j\leq n+2.
	\end{align*}
	Thus, since $\int_{-1}^{1}\exp\left[i\omega g(x)\right]dx=\Ob(\omega^{-1/2})$ we have overall
	\begin{align*}
	I_{\omega,\alpha,\beta}^{(2)}[f]=\Ob\left(\omega^{-(k+2)}\right),\text{\ as\ }\omega \rightarrow \infty.
	\end{align*}
\end{proof}

\subsection{The extended Filon method for Hankel kernels}\label{sec:extended_filon_method_for_Hankel_kernels}
This motivates the following construction of the Filon method. For the first type of integral $I_{\omega,\beta}^{(1)}[f]$ we begin by finding $p_1$, the unique polynomial of degree $2s+\nu+1$ such that
\begin{align}\label{eqn:interpolation_problem_p1}
p_1^{(j)}(0)=f^{(j)}(0), p_1^{(j)}(1)=f^{(j)}(1), j=0,\dots,s\,\,\,\,\text{\ and\ }\,\,\,\,p_1((c_n+1)/2)=f((c_n+1)/2), n=1,\dots,\nu,
\end{align}
where $c_n=\cos(n\pi/(\nu+1)),n=0,\dots,\nu+1$ are the Clenshaw--Curtis points. Then we approximate the integral by the quadrature rule $\Q^{(1)}_{[s,\nu]}[f]$, where
\begin{align*}
\Q ^{(1)}_{[s,\nu]}[f]=I_{\omega,\beta}^{(1)}[p_1].
\end{align*}
For the second type of integrals we compute $p_2$, the unique polynomial of degree $3s+\nu+1$ such that
\begin{align}\label{eqn:interpolation_problem_p2}
p_2^{(j)}(0)=f^{(j)}(0), p_2^{(j)}(\pm1)=f^{(j)}(\pm1), j=0,\dots,s\,\,\,\,\text{\ and\ }\,\,\,\,p_2(c_n)=f(c_n), n=1,\dots,\nu,
\end{align}
where $c_n$ are again the Clenshaw--Curtis points, \textit{and this time we require $\nu$ to be odd}. Then we approximate the integral by the quadrature rule $\Q^{(2)}_{[s,\nu]}[f]$, where
\begin{align*}
\Q ^{(2)}_{[s,\nu]}[f]=I_{\omega,\alpha,\beta}^{(2)}[p_2].
\end{align*}
\begin{remark}
	Of course to compute $I_{\omega,\beta}^{(1)}[p_1], I_{\omega,\alpha,\beta}^{(2)}[p_2]$ we are required to have a way of computing the moments corresponding to our basis of interpolation polynomials. We will choose a basis of \textit{shifted} Chebyshev polynomials, which requires us to find a way of computing
	\begin{align*}
	\sigma_n^{(1)}&=\int_{0}^{1}T_n\left(2x-1\right)H_0^{(1)}(\omega x)e^{i\omega\beta x}dx,\\
	\sigma_n^{(2)}&=\int_{-1}^{1}T_n(x)H_0^{(1)}\left(\omega \sqrt{(x-\alpha\beta)^2+\alpha^2(1-\beta^2)}\right)e^{i\omega\beta x}dx.
	\end{align*}
	We provide a stable way for computing these moments from recurrences in section \ref{sec:moment_computation_by_duality}.
\end{remark}
By the Filon paradigm described in Propositions \ref{prop:filon_paradigm_for_alpha_zero} \& \ref{prop:filon_paradigm_for_alpha_non_zero} we find immediately the asymptotic error of these quadrature methods to be
\begin{align*}
\Q ^{(1)}_{[s,\nu]}[f]-I_{\omega,\beta}^{(1)}[f]&=\Ob_\delta\left(\omega^{-(s+2)+\delta}\right)\\
\Q ^{(2)}_{[s,\nu]}[f]-I_{\omega,\alpha,\beta}^{(2)}[f]&=\Ob\left(\omega^{-(s+2)}\right).
\end{align*}
for any $\delta>0$ as $\omega\rightarrow\infty$.
\subsubsection{Fast interpolation at Clenshaw--Curtis points and endpoint derivative values}\label{sec:fast_interpolation_CC_points_sigma1}
It was shown by \cite{Gao2017a} that the interpolation problem \eqref{eqn:interpolation_problem_p1} for $p_1$ in a basis of (shifted) Chebyshev polynomials can be very efficiently solved by a single application of a \textit{Discrete Cosine Transform}, DCT-I, and solving a small auxiliary linear system of size $2s\times 2s$. We summarize their results here very briefly and extend those to the case of interpolating also at midpoints in the next section. We write

\begin{align*}
p_1(x)=\sum_{n=0}^{\nu+2s+1}p_nT_n(2x-1).
\end{align*}

Let us denote by $\mathcal{C}_{\nu+1}\mathbf{u}$ the DCT-I of the vector $\mathbf{u}=(u_0,\dots,u_{\nu+1})$. The inverse DCT-I then takes the form:
\begin{align}\label{eqn:def_DCT-I}
\mathcal{C}_{\nu+1}^{-1}\mathbf{u} &= \frac{2}{\nu+1}\sideset{}{''}\sum_{k=0}^{\nu+1} u_k \cos \left[\frac{nk\pi}{\nu+1} \right] \quad \quad \text{for\ }n = 0, \dots, \nu+1,
\end{align}
where $\displaystyle\sideset{}{''}\sum_{k=0}^{\nu+1}a_k$ simply means that for $k=0$ and $k=\nu+1$ the terms are halved. Then let $\mathbf{\check{p}}=(\check{p}_0,\dots,\check{p}_{\nu+1})$ be the inverse DCT-I of $\mathbf{f}=(f((c_0+1)/2),\dots,f((c_{\nu+1}+1)/2))$, i.e.
\begin{align*}
\check{p}_n=\mathcal{C}_{\nu+1}^{-1}\mathbf{f},
\end{align*}
then
\begin{align}
\begin{split}\label{eqn:form_of_pn_without_derivatives}
p_n&=\frac{1}{2}\check{p}_n\quad\quad n=0,\nu+1\\
p_n&=\check{p}_n\quad\quad n=1,\dots,\nu-2s\\
p_n&=\check{p}_n-p_{2\nu-n+2}\quad\quad n=\nu-2s+1,\dots,\nu
\end{split}
\end{align}
and the coefficients $p_{\nu+2},\dots,p_{\nu+2s+1}$ are given as the solution to the following well-posed linear \mbox{system}:
\begin{align}
\begin{split}\label{eqn:linear_system_for_chebyshev_coefficients}
\sum_{n=1}^{2s}&\left[\frac{(\nu+n+1)(\nu+n+j)!}{(\nu+n-j+1)!}-\frac{(\nu-n+1)(\nu-n+j)!}{(\nu-n-j+1)!}\right]p_{\nu+n+1}\\
&\hspace{4cm}=\frac{(2j)!}{2^jj!}\left[\frac{f^{(j)}(1)}{2^j}-\sideset{}{''}\sum_{n=0}^{\nu+1}T^{(j)}_n(1)\check{p}_n\right]\\
\sum_{n=1}^{2s}(-1)^{n}&\left[\frac{(\nu+n+1)(\nu+n+j)!}{(\nu+n-j+1)!}-\frac{(\nu-n+1)(\nu-n+j)!}{(\nu-n-j+1)!}\right]p_{\nu+n+1}\\
&\hspace{4cm}=(-1)^{\nu+j+1}\frac{(2j)!}{2^jj!}\left[\frac{f^{(j)}(0)}{2^j}-\sideset{}{''}\sum_{n=0}^{\nu+1}T^{(j)}_n(-1)\check{p}_n\right]
\end{split}
\end{align}
for $j=1,\dots,s$.

\subsubsection{Fast interpolation at Clenshaw--Curtis points, mid- and endpoint derivative values}\label{sec:fast_interpolation_sigma2}
It is indeed possible to extend the above results to the interpolation problem \eqref{eqn:interpolation_problem_p2} for $p_2$. In particular the Chebyshev coefficients of $p_2$ can be found through a single application of DCT-I and the solution of a small auxiliary system of size $3s\times 3s$ (a full proof is given in appendix A). To begin with, let us highlight that in this case we take $\nu$ to be odd. This is such that $c_{(\nu+1)/2}=0$ is among the Clenshaw--Curtis points. Note that under these circumstances this is a Hermite interpolation problem and consequently well-posed with a unique solution (see for instance chapter 5.5 in \cite{powell1981}).

Let us denote the coefficients again by $p_n$, such that (in this case no shift is required)
\begin{align*}
p_2(x)=\sum_{n=0}^{\nu+3s+1}p_nT_n(x).
\end{align*}
Let further be $\mathbf{\check{p}}=(\check{p}_0,\dots,\check{p}_{\nu+1})$ the DCT-I of $\mathbf{f}=(f(c_0),\dots,f(c_{\nu+1}))$, i.e.
\begin{align*}
\mathbf{\check{p}}=\mathcal{C}_{\nu+1}^{-1}\mathbf{f},
\end{align*}
then
\begin{align*}
p_n&=\frac{1}{2}\check{p}_n\quad\quad n=0,\nu+1\\
p_n&=\check{p}_n\quad\quad n=1,\dots,\nu-3s\\
p_n&=\check{p}_n-p_{2\nu-n+2}\quad\quad n=\nu-3s+1,\dots,\nu
\end{align*}
and the coefficients $p_{\nu+2},\dots,p_{\nu+3s+1}$ are given as the solution to the following well-posed linear system:
\begin{align*}
\sum_{n=1}^{3s}p_{\nu+1+n}\left[T_{\nu+1+n}^{(j)}(-1)-T_{\nu+1-n}^{(j)}(-1)\right]&=f^{(j)}(-1)-\sideset{}{''}\sum_{n=0}^{\nu+1}\check{p}_nT_n^{(j)}(-1)\\
\sum_{n=1}^{3s}p_{\nu+1+n}\left[T_{\nu+1+n}^{(j)}(0)-T_{\nu+1-n}^{(j)}(0)\right]&=f^{(j)}(0)-\sideset{}{''}\sum_{n=0}^{\nu+1}\check{p}_nT_n^{(j)}(0)\\
\sum_{n=1}^{3s}p_{\nu+1+n}\left[T_{\nu+1+n}^{(j)}(1)-T_{\nu+1-n}^{(j)}(1)\right]&=f^{(j)}(1)-\sideset{}{''}\sum_{n=0}^{\nu+1}\check{p}_nT_n^{(j)}(1).
\end{align*}
Note that the coefficients in this linear system can be easily determined to be
\begin{align*}
T_n^{(j)}(\pm1)&=(\pm 1)^{n-j}\frac{2^j j! n(n+j-1)!}{(2j)!(n-j)!}, \quad 0\leq j\leq n,n+j\geq 1
\end{align*}
and the values of $T_{n}^{(j)}(0), j=1,\dots,s$ can be efficiently computed from the recurrence
\begin{align*}
T_{n+1}^{(j)}(0)&=2T_n^{(j-1)}(0)-T_{n-1}^{(j)}(0), \quad T_{0}(x)=1,T_{1}(x)=x.
\end{align*}

\section{Moment computation via duality}\label{sec:moment_computation_by_duality}
A crucial step in the Filon method described in section \ref{sec:extended_filon_method_for_Hankel_kernels} is the computation of the moments
\begin{align*}
\sigma_n^{(1)}&=\int_{0}^{1}T_n\left(2x-1\right)H_0^{(1)}(\omega x)e^{i\omega\beta x}dx,\\
\sigma_n^{(2)}&=\int_{-1}^{1}T_n(x)H_0^{(1)}\left(\omega \sqrt{(x-\alpha\beta)^2+\alpha^2(1-\beta^2)}\right)e^{i\omega\beta x}dx.
\end{align*}
In order to do so we make a few general observations. We note that a similar way of deriving a recurrence for the moments can be found in previous work by \cite{keller2007method}, but in Theorem \ref{thm:recurrence_by_duality} we extend Keller's results to more general Hilbert bases. This means that ultimately our method for computing moments can be applied to a more general setting, although in the particular instance of this paper we only require the result for Chebyshev polynomials.
\begin{theorem}\label{thm:recurrence_by_duality}
	Let $\{\phi_n\}_{n=0}^{\infty}$ be a Hilbert basis of $L^2([a,b],W)$ with the non-negative weight function $W(x)$. Consider the moments
	\begin{align*}
	\sigma_n=\int_a^b \phi_n(x)h(x)W(x)dx
	\end{align*}
	where $\mathcal{L}h=0$, for some linear differential operator $\mathcal{L}$. Suppose further that the action of $\mathcal{L}$  on the basis functions is given by a banded (infinite) matrix $B_{mn}$ with bandwidth $k$, such that
	\begin{align*}
	\mathcal{L}\phi_n=B_{mn}\phi_n,
	\end{align*}
	then the moments satisfy a $k+1$-term recurrence relation,
	\begin{align*}
	B^T_{mn}\sigma_n=0,
	\end{align*}
	which together with $k$ initial or boundary conditions uniquely determines all moments.
\end{theorem}
\begin{proof}
	By orthonormality the moments $\sigma_n$ are the coefficients of $h$ in this basis, thus we can solve $\mathcal{L}h=0$ by expanding in the basis $\phi_n$:
	\begin{align*}
	\sum_n \sigma_n\mathcal{L}\phi_n=\mathcal{L}h&=0,\\
	\text{\ thus\ \ }\sum_n\sigma_n \left(\sum_m B_{nm}\phi_m\right)&=0.
	\end{align*}
	Therefore $\sum_m\sum_n\left(\sigma_n  B_{nm}\right)\phi_m=0$ and the result follows since the functions $\phi_m$ are orthogonal.
\end{proof}
Theorem \ref{thm:recurrence_by_duality} provides a generalisation of Lemma 2.4 from \cite{keller2007method} and of results from \cite{Lewanowicz1991}.
\subsection{Recurrence relations for the Chebyshev moments}
The Chebyshev moments
\begin{align*}
\sqrt{s_n}\sigma_n^{(j)}=\sqrt{s_n}\int_{-1}^{1}T_n(y)h^{(j)}_\omega(y)dy
\end{align*}
correspond to the Chebyshev coefficients of $\sqrt{1-y^2}h^{(j)}_\omega(y)$, $j=1,2$. In this definition we introduced the scaling constants which ensure that our basis functions $\left\{\sqrt{s_n}T_n\right\}_{n=0}^\infty$, with $s_0=1/\pi,s_n=2/\pi,n\geq1$, are normalised with respect to the weight $1/\sqrt{1-y^2}$. Furthermore for our two types of integrals we have
\begin{align*}
h^{(1)}_\omega(y)&=\frac{1}{2}H_0^{(1)}\left(\omega \frac{y+1}{2}\right)e^{i\omega\beta (y+1)/2},\\
h^{(2)}_\omega(y)&=H_0^{(1)}\left(\omega \sqrt{(y-\alpha\beta)^2+\alpha^2(1-\beta^2)}\right)e^{i\omega\beta y}.
\end{align*}
Suppose that
\begin{align}\label{eqn:operator_de_for_weight_functions}
L^{(j)}_\omega\left[\sqrt{1-y^2}h^{(j)}_\omega(y)\right]=0,\quad j=1,2,
\end{align}
then computing the Chebyshev coefficients is equivalent to solving the above equation with a spectral method. From standard trigonometric identities one can show:
\begin{lemma}The Chebyshev polynomials satisfy the following relations for all $n\geq 1$:
	\begin{align*}
	yT_0(y)&= T_1(y),\quad\quad\quad  yT_n(y)=\frac{1}{2}T_{n-1}(y)+\frac{1}{2}T_{n+1}(y)\\
	(1-y^2)T_0'(y)&=0,\quad\quad (1-y^2)T_n'(y)=\frac{n}{2}T_{n-1}(y)-\frac{n}{2}T_{n+1}(y)
	\end{align*}
	In particular the actions of $y,(1-y^2)d/dy$ on $\left\{\sqrt{s_n}T_n\right\}_{n=0}^\infty$ are both banded, with bandwidth 3.
\end{lemma}
We now observe that from Bessel's equation it can be seen that our weight functions $h_\omega^{(j)}$ satisfy the differential equation in \eqref{eqn:operator_de_for_weight_functions} with the following differential operators:
\begin{align}\label{eqn:L_omega_1}
\begin{split}
L_\omega^{(1)}&=\left((1-y^2)\frac{d}{dy}\right)^2+i  \left(\beta  \omega \left(y^2-1\right)-i (3
y+1)\right)\left(\left(1-y^2\right) \frac{d}{dy}\right)\\
&\quad-\frac{1}{4} \left((\beta ^2-1) \omega^2 \left(y^2-1\right)^2-2 i \beta  \omega
(y-1) (y+1)^2-4 \left(y^2+y+1\right)\right)
\end{split}
\end{align}
and

\begin{align}\label{eqn:L_omega_2}
\begin{split}
L_\omega^{(2)}&=
(y-\alpha  \beta ) \left(\alpha ^2-2 \alpha  \beta
y+y^2\right)\left((1-y^2)\frac{d}{dy}\right)^2\\
&\quad-\left[2 \alpha ^3 \beta  \left(2 y+i \beta  \omega
\left(y^2-1\right)\right)\right.\\
&\quad\quad+\alpha ^2 \left(-2 i \left(2 \beta ^3+\beta \right)
\omega  y^3+2 i \left(2 \beta ^3+\beta \right) \omega  y-\left(6 \beta
^2+5\right) y^2-2 \beta ^2+1\right)\\
&\quad\quad\left.+2 \alpha  \beta  y \left(3 i \beta  \omega
\left(y^2-1\right) y+5 y^2+1\right)+y^2 \left(-2 i \beta  \omega
\left(y^2-1\right) y-3 y^2-1\right)\right]\left(\left(1-y^2\right) \frac{d}{dy}\right)\\
&\quad+\left[-2 \alpha ^3 \beta  y^2-\alpha ^3 \beta +2 \alpha ^2 \beta ^2 y^3+4 \alpha
^2 \beta ^2 y+3 \alpha ^2 y^3-4 \alpha  \beta  y^4-5 \alpha  \beta  y^2+y^5+2
y^3\right]\\
&\quad+\omega i\beta\left(1-y^2\right) \left[2 \alpha ^3 \beta y -2 \alpha ^2 \beta ^2
\left(1+y^2\right)+\alpha ^2 \left(-3 y^2+1\right)+2  \alpha
\beta  y \left(2y^2+1\right)-   y^2 \left(1+y^2\right)\right]\\
&\quad+\omega^2(1-\beta^2)y\left(1-y^2\right)^2\left[ 2 \alpha ^2 \beta ^2
-3 \alpha
\beta  y +y^2
\right]
\end{split}
\end{align}

Therefore their actions on the basis $\left\{\sqrt{s_n}T_n\right\}_{n=0}^\infty$ are represented by banded matrices $B^{(1)}_{mn}, B^{(2)}_{mn}$, of bandwidth eight and fourteen respectively. Let us write $\sigma^{(j)}_{-n}:=\sigma^{(j)}_{n}, n\geq1$, then
the condition $	B^T_{mn}\sigma_n=0$ can be conveniently expressed in terms of the following recurrences (valid for all $n\in\mathbb{Z},n\neq 0$):
\begin{align}\label{eqn:recurrence_for_sigma_1}
\begin{split}
&\frac{\omega ^2 (1-\beta^2) }{4 n^2}\sigma^{(1)}_n+\omega i\beta\left[\frac{2 i \beta }{n}-\frac{i \beta }{n^2}\right]\sigma^{(1)}_{n-1}+\left[4-\frac{8}{n}+\frac{\omega^2(\beta^2-1)+\omega 2 i\beta+4}{n^2}\right]\sigma^{(1)}_{n-2}\\&
+\left[\frac{\omega (-6i\beta)+8}{n}+\frac{\omega 17 i \beta-16}{n^2}\right]\sigma^{(1)}_{n-3}+\left[-8+\frac{64}{n}+\frac{\omega^2 3(1-\beta^2)+\omega(-8i\beta)-208}{2n^2}\right]\sigma^{(1)}_{n-4}\\&
+\left[\frac{\omega 6i\beta -8}{n}+\frac{\omega(-31i\beta)+48}{n^2}\right]\sigma^{(1)}_{n-5}+\left[4-\frac{56}{n}+\frac{\omega^2(\beta^2-1)+\omega 2i\beta+196}{n^2}\right]\sigma^{(1)}_{n-6}\\&
+\omega i \beta \left[-\frac{2  }{n}+\frac{15  }{n^2}\right]\sigma^{(1)}_{n-7}+\frac{\omega ^2 (1-\beta^2) }{4
	n^2}\sigma^{(1)}_{n-8}=0,
\end{split}
\end{align}
and
\begin{align*}
\begin{split}
&\frac{\omega ^2 (1-\beta^2)}{n^2}\sigma^{(2)}_n+\left[\frac{4 i \beta  \omega
}{n}+\frac{6\omega ^2  \alpha  \beta (\beta^2-1) -2 i \beta  \omega }{n^2}\right]\sigma^{(2)}_{n-1}\\&
+\left[4+\frac{ -24i\omega  \alpha  \beta ^2-8}{n}+\frac{\omega ^2(8 \alpha ^2 \beta ^2-1)(1-\beta^2)+32i\omega   \alpha  \beta ^2 +4}{n^2}\right]\sigma^{(2)}_{n-2}\\&
+\left[-24 \alpha  \beta+\frac{16i\omega \alpha^2\beta (2\beta ^2+1) +88 \alpha  \beta }{n}\right.\\
&\left.+\frac{12\omega ^2 \alpha  \beta (1-\beta^2) +i\omega \beta(-80  \alpha ^2 \beta ^2  -24  \alpha ^2  +12    )-80 \alpha
	\beta}{n^2}\right]\sigma^{(2)}_{n-3}\\&
+\left[4+32 \alpha ^2 \beta ^2+16 \alpha ^2+\frac{ 8i\omega \alpha  \beta ^2(3 -4  \alpha ^2)+8-64 \alpha^2-192\alpha ^2 \beta ^2}{n}\right.\quad\quad\quad\quad\quad\quad\\&
+\left.\frac{\omega^2(1-\beta^2)(-3)(8\alpha ^2 \beta ^2+1)+\omega(96 i \alpha ^3 \beta ^2 -144 i \alpha  \beta ^2)+288\alpha^2\beta^2+48\alpha^2-44 }{n^2}\right]\sigma^{(2)}_{n-4}\\
\end{split}
\end{align*}
\begin{align}
\begin{split}\label{eqn:recurrence_for_sigma_2}&
+\left[-32 \alpha ^3 \beta +\frac{-2i\omega \beta( 32 \alpha ^2 \beta ^2 +16\alpha ^2+6) + 224 \alpha ^3 \beta-176\alpha  \beta}{n}\right.\quad\quad\quad\quad\quad\quad\\&
\left.  +\frac{6\omega ^2 \alpha\beta (1-\beta^2)+i\omega \beta(128 \alpha ^2   +58+384 \alpha ^2 \beta ^2)-384 \alpha ^3 \beta+592 \alpha  \beta}{n^2}\right]\sigma^{(2)}_{n-5}\\
&+\left[-8 -16 \alpha ^2-32 \alpha ^2 \beta ^2 +\frac{48i\omega   \alpha  \beta ^2( 1+2\alpha ^2)+576 \alpha ^2 \beta ^2 +192\alpha ^2+128}{n}\right.\quad\quad\quad\quad\quad\quad\\&
+\left.\frac{\omega ^2(1-\beta^2) (3+16\alpha ^2 \beta ^2)+i\omega \alpha\beta^2(-544 \alpha ^2-224 )-1952 \alpha ^2 \beta ^2-432 \alpha ^2-344}{n^2}\right]\sigma^{(2)}_{n-6}\\&
+\left[64 \alpha ^3 \beta+48 \alpha  \beta +\alpha  \beta \frac{-672  -896 \alpha ^2 }{n}\right.\quad\quad\quad\quad\quad\quad\\&
+\left.\frac{-24\omega^2\alpha  \beta (1-\beta^2) +i\omega \beta(16\alpha^2 -160\alpha^2 \beta ^2 -24)+2880 \alpha ^3 \beta+1840 \alpha  \beta}{n^2}\right]\sigma^{(2)}_{n-7}\\&
+\left[-8-16 \alpha ^2-32 \alpha ^2 \beta ^2+\frac{-48i\omega\alpha  \beta ^2 ( 1+2 \alpha ^2)+320 \alpha ^2 \beta^2+256\alpha^2+96}{n}\right.\quad\quad\quad\quad\quad\quad\\&
+\left.\frac{\omega^2(1-\beta^2)(16\alpha ^2 \beta ^2+3)+i\omega \alpha  \beta ^2 (800  \alpha ^2+448 )-160 \alpha ^2 \beta ^2-880 \alpha ^2-120}{n^2}\right]\sigma^{(2)}_{n-8}\\&
+\left[-32 \alpha ^3 \beta +\frac{i\omega
	\beta (12 +64 \alpha ^2 \beta ^2 +32  \alpha ^2)+672\alpha ^3 \beta +176 \alpha  \beta}{n}\right.\quad\quad\quad\quad\quad\quad\\&
\left.+\frac{6\omega ^2 \alpha \beta  (1-\beta^2)+i\omega\beta(-512\alpha ^2 \beta ^2 -110 -320 \alpha ^2 )-3520 \alpha ^3 \beta-1872 \alpha  \beta}{n^2}\right]\sigma^{(2)}_{n-9}\\&
+\left[4+32 \alpha ^2 \beta ^2 +16 \alpha ^2 +\frac{8i\omega \alpha\beta^2(4  \alpha ^2 -3) -120  -704 \alpha
	^2 \beta ^2 -384 \alpha ^2  }{n}\right.\quad\quad\quad\quad\quad\quad\\&
+\left.\frac{\omega ^2(1-\beta^2)(-24 \alpha ^2 \beta ^2-3)+\omega(192 i \alpha  \beta ^2 -352 i \alpha ^3 \beta ^2 )+852   +3872 \alpha ^2 \beta ^2 +2288 \alpha ^2 }{n^2}\right]\sigma^{(2)}_{n-10}\\&
+\left[-24 \alpha  \beta +\frac{-16i\omega \alpha ^2 \beta  (1+2\beta^2) +584 \alpha  \beta}{n}\right.\\
&+\left.\frac{12 \omega ^2  \alpha  \beta (1-\beta^2) +i\omega  \beta(12  +368 \alpha ^2 \beta ^2 +200 \alpha ^2  )  -3552 \alpha
	\beta   }{n^2}\right]\sigma^{(2)}_{n-11}\\&
+\left[4+\frac{24i \omega \alpha  \beta ^2   -104 }{n}+\frac{\omega ^2(1-\beta^2)(8 \alpha ^2 \beta ^2-1)-304 i \alpha  \beta ^2 \omega+676}{n^2}\right]\sigma^{(2)}_{n-12}\\&
+ \left[-\frac{4 i \beta  \omega
}{n}+\frac{-6\omega^2 \alpha  \beta(1-\beta^2)+54 i \beta  \omega }{n^2}\right]\sigma^{(2)}_{n-13}+\frac{\omega ^2 (1-\beta^2)}{n^2}\sigma^{(2)}_{n-14}=0,
\end{split}
\end{align}
where we divided by $n^{2}$ for easier access to the dominant terms for our later analysis of the stability of these recurrences.
\begin{remark}
	Although the recurrences \eqref{eqn:recurrence_for_sigma_1} \& \eqref{eqn:recurrence_for_sigma_2} involve eight and fourteen terms each, the fact that they are also valid when $n=4,5,6,7$ (for the first case) and $n=7,\dots,13$ (for the second case) allows us to start the moment computation from just four and seven initial values respectively. This is a significant improvement over having to start with eight and fourteen initial values, and efficient expressions for computing those are provided in section \ref{sec:initial_conditions}.
\end{remark}

\subsection{Initial conditions}\label{sec:initial_conditions}
In order to be able to compute the moments we need to find simple expressions for the initial moments.
\subsubsection{Initial conditions for $I_{\omega,\beta}^{(1)}[f]$}
Therefore let us begin by proving an expression for $\sigma_0^{(1)}$ that allows for efficient and accurate computation.
\begin{lemma} For $\omega>0$:
	\begin{align}\label{eqn:closed_form_expression_rho_0}
	\sigma_0^{(1)}&=-\frac{2i}{\pi}\omega^{-1}\int_{0}^{\infty} \frac{1}{\sqrt{t}}\frac{e^{i(\beta+1)\omega}e^{-t\omega}-1}{\sqrt{2i-t}(1+\beta+it)}dt.
	\end{align}
\end{lemma}
\begin{proof}Denoting the Fourier transform by $\mathcal{F}$, and noting due to the oscillations and $\mathcal{O}(x^{-\frac{1}{2}})$ decay of $H^{(1)}_0$ the Fourier transform is well-defined as an integral (where we take in the following WLOG $H^{(1)}_0(-x):=H^{(1)}_0(x)$) we have
	\begin{align}\label{eqn:convolution_expression_rho_0}
	\sigma_0^{(1)}(\omega)&=\omega^{-1}\int_{0}^{\omega}H^{(1)}_0( x)e^{i \beta x}dx=\omega^{-1}\left(\mathcal{F}H^{(1)}_0\star\mathcal{F}\mathbbm{1}_{[0,\omega]} \right)(\beta).
	\end{align}
	Now we have
	\begin{align*}
	\mathcal{F}\mathbbm{1}_{[0,\omega]} (k)=\frac{(-i)(-1+e^{ik\omega})}{k\sqrt{2\pi}}.
	\end{align*}
	And to find the Fourier transform of our (even version of the) Hankel function we note that for $\alpha\notin\mathbb{Z}$ we have
	\begin{align*}
	H^{(1)}_\alpha(x)=\frac{J_{-\alpha}(x)-e^{i\alpha\pi}J_\alpha(x)}{i\sin(\alpha\pi)},
	\end{align*}
	and for $\alpha$ any real number:
	\begin{align*}
	\mathcal{F}J_\alpha(k)=\begin{cases}
	\frac{\textstyle(-1)^{-\alpha} \left(\sqrt{k^2-1}-k\right)^\alpha \left(\left((-1)^\alpha+1\right) \sin
		\left(\frac{\pi  \alpha}{2}\right)-i \left((-1)^\alpha-1\right) \cos \left(\frac{\pi
			\alpha}{2}\right)\right)}{\textstyle\sqrt{2 \pi } \sqrt{k^2-1}} & k<-1 \\
	\frac{\textstyle(-1)^{-\alpha} \left(\left((-1)^\alpha+1\right) \cos \left(\alpha \sin ^{-1}(k)\right)-i
		\left((-1)^\alpha-1\right) \sin \left(\alpha \sin ^{-1}(k)\right)\right)}{\textstyle\sqrt{2 \pi }
		\sqrt{1-k^2}} & -1<k<0,\\
	&\text{or } 0<k<1\\
	0&\text{otherwise.}
	\end{cases}
	\end{align*}
	Thus for $\alpha\notin\mathbb{Z}$:
	\begin{align}\label{eqn:Fourier_transform_of_H_alpha}
	\displaystyle\mathcal{F}H^{(1)}_\alpha(k)=\begin{cases}
	-\frac{\displaystyle i e^{-\frac{3}{2} i \pi  \alpha} \left(-1+e^{2 i \pi  \alpha}\right)
		\left(\sqrt{k^2-1}-k\right)^{-\alpha} \left(\left(\sqrt{k^2-1}-k\right)^{2
			\alpha}+1\right)}{\displaystyle i\sqrt{2 \pi }\sin(\alpha\pi) \sqrt{k^2-1}} & k<-1 \\
	0 & \text{otherwise.}
	\end{cases}
	\end{align}
	In this expression we are able to take the limit $\alpha\rightarrow 0^+$ to find
	\begin{align} \label{eqn:Fourier_transform_of_H_0}
	\mathcal{F}H^{(1)}_\alpha(k)=\begin{cases}
	-\frac{\displaystyle 2 i \sqrt{2/\pi}}{\displaystyle\sqrt{k^2-1}} & k<-1 \\
	0 & \text{otherwise.}
	\end{cases}
	\end{align}
	Taking this limit is justified by regarding the functions $H^{(1)}_\alpha$ as temperate distributions. $\mathcal{F}$ is a continuous function on the space of temperate distributions $\mathscr{S}'$ and $H^{(1)}_\alpha\stackrel{\alpha\rightarrow0}{\rightarrow}H^{(1)}_0$ in the topology of $\mathscr{S}'$. Similarly the RHS of \eqref{eqn:Fourier_transform_of_H_alpha} converges to the RHS of \eqref{eqn:Fourier_transform_of_H_0} in $\mathscr{S}'$, so \eqref{eqn:Fourier_transform_of_H_0} must indeed hold. Thus by \eqref{eqn:convolution_expression_rho_0} we have
	\begin{align*}
	\sigma_0^{(1)}&=-\frac{2}{\pi}\omega^{-1}\int_{-\infty}^{-1} \frac{1}{\sqrt{y^2-1}}\frac{e^{i(\beta-y)\omega}-1}{(\beta-y)}dy.
	\end{align*}
	By the decay of the integrand (which has only a square root singularity at $y=-1$ which can be dealt with in the usual way, by excluding a small neighbourhood of the origin during the change of contour), we can deform the contour of integration to the following, where we take the branch of the square root to be such that $\Re\sqrt{z}\geq0$:
	\begin{align*}
	\sigma_0^{(1)}(\omega)&=-\frac{2i}{\pi}\omega^{-1}\int_{0}^{\infty} \frac{1}{\sqrt{t}}\frac{e^{i(\beta+1)\omega}e^{-t\omega}-1}{\sqrt{2i-t}(1+\beta+it)}dt
	\end{align*}
	as required.
\end{proof}

We can further simplify the expression found in \eqref{eqn:closed_form_expression_rho_0}, by changing variables and performing part of the integral exactly in the case when $\beta>-1,\beta\neq 1$:
\begin{align*}
\sigma_0^{(1)}(\omega)&=-\frac{2i}{\pi}\omega^{-1}\int_{0}^{\infty} \frac{1}{\sqrt{t}}\frac{e^{i(\beta+1)\omega}e^{-t\omega}}{\sqrt{2i-t}(1+\beta+it)}dt+\frac{2}{\pi}\omega^{-1}\frac{2 \tanh ^{-1}\left(\sqrt{\frac{\beta-1}{\beta+1}}\right)}{\sqrt{\beta^2-1}}\\
&=-\frac{2i}{\pi}\omega^{-\frac{3}{2}}e^{i(\beta+1)\omega}\int_{0}^{\infty}f_\beta\left(\frac{t}{\omega}\right) \frac{1}{\sqrt{t}}e^{-t}dt+\frac{2}{\pi}\omega^{-1}\frac{2 \tanh ^{-1}\left(\sqrt{\frac{\beta-1}{\beta+1}}\right)}{\sqrt{\beta^2-1}},
\end{align*}
where $f_\beta(t)=\frac{1}{\sqrt{2i-t}(1+\beta+it)}$. Let us now define the standard moments by
\begin{align*}
\rho^{(1)}_n:=I^{(1)}_{\omega,\beta}[x^n]
\end{align*}
then we find similarly through differentiating with respect to $\beta$:
\begin{corollary}\label{cor:initial_moments_sigma_1}
	Letting again $f_\beta(t)=\frac{1}{\sqrt{2i-t}(1+\beta+it)}$, we have
	\begin{align*}
	\rho^{(1)}_0&=\sigma^{(1)}_0=-\frac{2}{\pi}\omega^{-1}\begin{cases}
	i\omega^{-\frac{1}{2}}e^{i(\beta+1)\omega}\int_{0}^{\infty}f_\beta\left(\frac{t}{\omega}\right) \frac{1}{\sqrt{t}}e^{-t}dt+\begin{cases}-
	\frac{2 \tanh ^{-1}\left(\sqrt{\frac{\beta-1}{\beta+1}}\right)}{\sqrt{\beta^2-1}},&\beta>-1,\beta\neq1\\
	-1,&\beta=1\\
	\frac{2 \tanh ^{-1}\left(\sqrt{\frac{\beta-1}{\beta+1}}\right)}{\sqrt{\beta^2-1}},&\beta<-1
	\end{cases}\\
	i\int_{0}^\infty\sqrt{2i-t/\omega}\frac{1}{\sqrt{t}}e^{-t}dt+1,&\hspace{-2.55cm}\beta=-1
	\end{cases}
	\end{align*}
	\begin{align*}
	i\omega\rho_1&=-\frac{2}{\pi}\omega^{-1}\begin{cases}
	i\omega^{-\frac{1}{2}}e^{i(\beta+1)\omega}\int_{0}^{\infty}\left[i\omega f_\beta\left(\frac{t}{\omega}\right)+\left(\partial_{\beta}f_\beta\left(\frac{t}{\omega}\right)\right)\right] \frac{1}{\sqrt{t}}e^{-t}dt\\
	\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+\begin{cases}\frac{1}{1-\beta ^2}+\frac{2 \beta  \tanh ^{-1}\left(\sqrt{\frac{\beta -1}{\beta
				+1}}\right)}{(\beta^2 -1)^{3/2}},&\beta>-1,\beta\neq 1\\
	1/3,&\beta=1\\
	\frac{1}{1-\beta ^2}-\frac{2 \beta  \tanh ^{-1}\left(\sqrt{\frac{\beta -1}{\beta+1}}\right)}{(\beta^2 -1)^{3/2}},&\beta<-1
	\end{cases}\\
	-i/3\omega^{3/2}\int_{0}^{\infty}(-i-t/\omega)\sqrt{2i-t/\omega}\frac{1}{\sqrt{t}}e^{-t}dt+1/3,&\hspace{-2.55cm}\beta=-1
	\end{cases}
	\end{align*}
	\begin{align*}
	-\omega^2\rho_2&=-\frac{2}{\pi}\omega^{-1}\begin{cases}
	i\omega^{-\frac{1}{2}}\int_{0}^{\infty}\left[-\omega^2 f_\beta\left(\frac{t}{\omega}\right)+2i\omega\left(\partial_{\beta}f_\beta\left(\frac{t}{\omega}\right)\right)+\left(\partial_{\beta}^2f_\beta\left(\frac{t}{\omega}\right)\right)\right] \frac{1}{\sqrt{t}}e^{-t}dt\\
	\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad +\begin{cases}
	\frac{3 \beta }{\left(\beta ^2-1\right)^2}-\frac{\left(4 \beta ^2+2\right) \tanh
		^{-1}\left(\sqrt{\frac{\beta -1}{\beta +1}}\right)}{(\beta^2 -1)^{5/2}},&\beta>-1,\beta\neq1\\
	-4/15,&\beta=1\\
	\frac{3 \beta }{\left(\beta ^2-1\right)^2}+\frac{\left(4 \beta ^2+2\right) \tanh
		^{-1}\left(\sqrt{\frac{\beta -1}{\beta +1}}\right)}{(\beta^2 -1)^{5/2}},&\beta<-1
	\end{cases}\\
	\omega^{5/2}i/15\int_{0}^\infty\sqrt{2i-t/\omega}(2(t/\omega)^2+2it/\omega-3)\frac{1}{\sqrt{t}}e^{-t}dt+4/15,&\hspace{-2.55cm}\beta=-1
	\end{cases}\\
	\end{align*}
	and when $\beta\neq\pm1$:
	\begin{align*}
	-i\omega^3\rho_3=& -\frac{2i}{\pi}\omega^{-\frac{3}{2}}e^{i(\beta+1)\omega}\int_{0}^{\infty}\left[-i\omega^3 f_\beta\left(\frac{t}{\omega}\right)-3\omega^2\left(\partial_{\beta}f_\beta\left(\frac{t}{\omega}\right)\right)\right.\\
	&\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad \left.+3i\omega\left(\partial_{\beta}^2f_\beta\left(\frac{t}{\omega}\right)\right)+\left(\partial_{\beta}^3f_\beta\left(\frac{t}{\omega}\right)\right)\right] \frac{1}{\sqrt{t}}e^{-t}dt\\
	&\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad +\frac{2}{\pi}\omega^{-1}\begin{cases}
	\frac{11 \beta ^2+4}{\left(\beta ^2-1\right)^3}-\frac{6 \beta  \left(2 \beta
		^2+3\right) \tanh ^{-1}\left(\sqrt{\frac{\beta -1}{\beta +1}}\right)}{(\beta
		-1)^{7/2} (\beta +1)^{7/2}},&\beta>-1,\beta\neq1\\
	\frac{11 \beta ^2+4}{\left(\beta ^2-1\right)^3}+\frac{6 \beta  \left(2 \beta
		^2+3\right) \tanh ^{-1}\left(\sqrt{\frac{\beta -1}{\beta +1}}\right)}{(\beta
		-1)^{7/2} (\beta +1)^{7/2}},&\beta<-1
	\end{cases}
	\end{align*}
\end{corollary}

These expressions of the first moments as given in corollary \ref{cor:initial_moments_sigma_1} are useful because they involve only simple functions and exponentially decaying integrals. In particular integrals of the form
\begin{align*}
\int_{0}^{\infty}f\left(\frac{t}{\omega}\right) \frac{1}{\sqrt{t}}e^{-t}dt=\int_{-\infty}^{\infty}f\left(\frac{z^2}{\omega}\right)e^{-z^2}dz,
\end{align*}
which can be efficiently evaluated using Gauss--Hermite quadrature. The error $\mathcal{E}_m$ of Gaussian quadrature using $m$ points satisfies a bound of the form
\begin{align}\label{eqn:error_of_general_gaussian_quadrature}
|\mathcal{E}_m|&\leq C_m \sup_{[0,\infty)}\left|\left[\frac{d}{dx}\right]^{2m}f\left(\frac{x^2}{\omega}\right)\right|\\\nonumber
&\leq \omega^{-m}\tilde{C}_m\max_{m\leq j\leq 2m}\sup_{[0,\infty)}\left|f^{(j)}(x)\right|, \quad \omega\geq 1.
\end{align}
Thus we see that while Gaussian quadrature converges for all $\omega>0$ the convergence rate becomes particularly favourable for larger frequencies and evaluating
\begin{align*}
\sigma^{(1)}_0&=\rho^{(1)}_0\\
\sigma^{(1)}_1&=2\rho^{(1)}_1-\rho^{(1)}_0\\
\sigma^{(1)}_2&=8\rho^{(1)}_2-8\rho^{(1)}_1+\rho^{(1)}_0\\
\sigma^{(1)}_3&=32\rho^{(1)}_3-48\rho^{(1)}_2+18\rho^{(1)}_1-\rho^{(1)}_0
\end{align*}
in this way allows for to approximation of the initial conditions to high accuracy at frequency-independent cost. Note that Gauss--Laguerre quadrature performs well even for moderate values of $\omega$, but when $\omega\ll1$ Gauss--Laguerre requires too many quadrature nodes to capture the behaviour of $f$ along the contour of integration. However in this case the original integral $I^{(1)}_{\omega,\beta}$ is no longer oscillatory and thus alternative methods are available to evaluate this integral.
\subsubsection{Initial conditions for $I_{\omega,\alpha,\beta}^{(2)}[f]$}
In order to compute the moments for the second type of integral we transform the integral expressions to steepest descent contours.
\begin{remark}
	This essentially amounts to numerical steepest descent as described by \cite[Chapter 5]{deano2017computing}, and it is appropriate to apply this to our moments because we know the analyticity behaviour of the entire functions $T_n(x)$, in contrast to the analyticity of $f$ in $I_{\omega,\alpha,\beta}^{(2)}[f]$ which only needs to be differentiable up to a given order, and may not have an analytic extension. Thus applying steepest descent to $I_{\omega,\alpha,\beta}^{(2)}[f]$ directly would be more restrictive than calculating the moments in this way.
\end{remark}
In this setting we aim to compute
\begin{align*}
\sigma_n^{(2)}&=\int_{-1}^{1}\underbrace{T_n(x)h_0\left(\omega \sqrt{(x-\alpha\beta)^2+\alpha^2(1-\beta^2)}\right)}_{f_n(x)}\exp\Big[\omega\underbrace{i\left(\sqrt{(x-\alpha\beta)^2+\alpha^2(1-\beta^2)}+\beta x\right)}_{=:g(x)}\Big]dx
\end{align*}
where we choose the branch of the root
\begin{align*}
\sqrt{x^2-2\alpha\beta+\alpha^2}=\sqrt{x-\alpha\beta+\sqrt{\alpha^2(\beta^2-1)}}\sqrt{x-\alpha\beta-\sqrt{\alpha^2(\beta^2-1)}}
\end{align*}
such that it takes non-negative real part, i.e. we have two branch cuts extending along $\{z|\Re z=\alpha\beta,|\Im z|>\sqrt{\alpha^2(1-\beta^2)}\}$ as sketched in figure \ref{fig:steepest_descent_contour_sigma_2}.

The behaviour is also such that $\Re g(z)$ decreases in the first and third quadrant, so we can find a contour of steepest descent by connecting $-1$ and $1$ through a contour of three parts $\mathcal{C}_{-1}+\mathcal{C}_{0}+\mathcal{C}_{1}$, which are a combination of steepest descent contours (i.e. $\Im g(z)$ remains constant along each of the contours) of the following behaviour as seen in figure \ref{fig:steepest_descent_contour_sigma_2}.
\begin{figure}
	\centering
	\begin{tikzpicture}

	\draw[thin,->] (-5,0) -- (5,0) node[right] {$\Re z$};
	\draw[thin,->] (0,-4) -- (0,4) node[above] {$\Im z$};
	\draw[line width=0.75mm,decorate, decoration={snake, segment length=3mm, amplitude=0.5mm},-] (1,1) -- (1,4);
	\draw[line width=0.75mm,decorate, decoration={snake, segment length=3mm, amplitude=0.5mm},-] (1,-1) -- (1,-4);
		\draw [line width=1.5pt] plot [smooth, tension=1] coordinates { (-3,0) (-3.0,-1.5)(-3.2,-4)  };
			\draw [-angle 60, line width=1.5pt] (-3.05,-2.0)--(-3.06,-2.15)node[left] {$\mathcal{C}_{-1}$};
		\draw [line width=1.5pt] plot [smooth, tension=1] coordinates { (3,0) (3.0,1.5)(3.2,4)  };
			\draw [-angle 60, line width=1.5pt] (3.055,2.15)--(3.04,2.0) node[right] {$\mathcal{C}_{1}$};
	\draw [line width=1.5pt] plot [smooth, tension=0.5] coordinates { (-2,-4) (-1.5,-0.75)  (0,0) (1.5,0.75)  (2,4)};
	\draw [-angle 60, line width=1.5pt] (-1.0,-0.325)--(-.85,-0.24) node[below right] {$\mathcal{C}_{0}$};

\draw[thin] (1,0.12) -- (1,-0.12) node[below] {$\alpha\beta$};
\draw[thin] (0.12,1) -- (-0.12,1) node[left] {$\alpha\sqrt{1-\beta^2}$};
\draw[thin] (-3,-0.12) -- (-3,0.12) node[above] {$-1$};
\draw[thin] (3,0.12) -- (3,-0.12) node[below] {$1$};
	\end{tikzpicture}
	\caption{Steepest descent contour for $\sigma_n^{(2)}$. The curled lines represent that branch-cuts of $g(z)$.}
	\label{fig:steepest_descent_contour_sigma_2}
\end{figure}
\newpage
Then we change coordinates on each of them as follows:
\begin{itemize}
	\item $\mathcal{C}_{-1}$ starts at $-1$ and ends at complex $\infty$ with $\arg z\in [\pi,3\pi/2)$. Along this contour we apply the change-of-variable $g(-1)-it=g(x)$. Then $x=g^{-1}(g(-1)-it)$. Thus the contribution from this integral can be written as
	\begin{align*}
	I_1=\int_{-i\infty}^{0}\frac{f_n(g^{-1}(g(-1)-it))}{g'(g^{-1}(g(-1)-it))}e^{\omega g(-1)-i\omega t}idt=-\omega^{-1}\int_0^\infty\frac{f_n(g^{-1}(g(-1)-t\omega^{-1}))}{g'(g^{-1}(g(-1)-t\omega^{-1}))}e^{\omega g(-1)}e^{-t}dt.
	\end{align*}
	\item $\mathcal{C}_0$ connects $\infty$ in $\arg z\in[\pi,3\pi/2)$ with $\infty$ in $\arg z\in [0,\pi/2)$ and passes through the saddle point $0$. Along this contour we change variables to $g(0)+it^2=g(x)$, which allows us to write the contribution as
	\begin{align*}
	I_2=\int_{-\infty e^{i\pi/4}}^{\infty e^{i\pi/4}}\frac{f_n(g^{-1}(g(0)+it^2))}{g'(g^{-1}(g(0)+it^2))}e^{\omega g(0)+i\omega t^2}2itdt=-\omega^{-1}\int_{-\infty}^{\infty }\frac{f_n(g^{-1}(g(0)-t^2\omega^{-1}))}{g'(g^{-1}(g(0)-t^2\omega^{-1}))}e^{\omega g(0)}e^{-t^2}2tdt
	\end{align*}
	where we have chosen the branches of $g^{-1}$ appropriately such that the path $\mathcal{C}_{0}$ maps to the straight line $\arg z=\pi/4$.
	\item $\mathcal{C}_{1}$ starts from complex $\infty$ at $\arg z\in [0,\pi/2)$ and ends in $1$. Along this contour we apply the change-of-variable $g(-1)+it=g(x)$. Then $x=g^{-1}(g(1)+it)$. Thus the contribution from this integral can be written as
	\begin{align*}
	I_3=-\int_{0}^{i\infty}\frac{f_n(g^{-1}(g(1)+it))}{g'(g^{-1}(g(1)+it))}e^{\omega g(1)+i\omega t}idt=\omega^{-1}\int_0^\infty\frac{f_n(g^{-1}(g(1)-t\omega^{-1}))}{g'(g^{-1}(g(1)-t\omega^{-1}))}e^{\omega g(1)}e^{-t}dt.
	\end{align*}
\end{itemize}
Then in summary the moments can be written as the following sum of exponentially decaying integrals
\begin{align}\label{eqn:NSD_expressions_for_sigma2}
\sigma_n^{(2)}=I_1+I_2+I_3.
\end{align}
We note that $I_1,I_3$ can be computed using standard Gauss--Laguerre quadrature, whereas $I_2$ can be computed using Gauss--Hermite quadrature. The error satisfies a similar bound as \eqref{eqn:error_of_general_gaussian_quadrature} and hence this approximation improves in quality as $\omega$ increases. As before, Gaussian quadrature performs well even for moderate values of $\omega$, however, when $\omega\ll 1$, too many quadrature nodes are required for numerical steepest descent to be a suitable way of computing the initial conditions. We note that in this case the original integral is non-oscillatory and so could be evaluated with alternative methods. This is a commonly known feature of numerical steepest descent.

\subsection{Asymptotic behaviour of homogeneous solutions to the recurrences}
In theory the aforementioned initial conditions, together with the recurrences \eqref{eqn:recurrence_for_sigma_1} and \eqref{eqn:recurrence_for_sigma_2}, determine all moments uniquely. However, in practice the process of computing the moments using forward propagation on these recurrences is exponentially unstable after a certain cut-off value. To determine the underlying reason, the cut-off value and a way of overcoming this instability we consider the asymptotic behaviour of homogeneous solutions to the recurrences.
\subsubsection{Stability of forward propagation when $\omega\gg n$}\label{sec:stability_of_forward_propagation}
The dominant terms in this case arise from the $\omega^2$ terms in the recurrences, thus from \eqref{eqn:L_omega_1} we see that the action of $L^{(1)}_\omega$ is dominated by
\begin{align}\label{eqn:recurrence_1st_moments_omega_larger_n}
-\frac{1}{4}(\beta^2-1)\left(\mathbf{M}^2-1\right)^2\sigma^{(1)}_n\approx 0
\end{align}
and from \eqref{eqn:L_omega_2} the action of $L^{(2)}_\omega$ is dominated by
\begin{align}\label{eqn:recurrence_2nd_moments_omega_larger_n}
(1-\beta ^2)\mathbf{M}\left(\mathbf{M}^2-1\right)^2\left[2 \alpha ^2 \beta
^2  -3 \alpha
\beta  \mathbf{M}  +\mathbf{M}^2\right]\sigma_n^{(2)}\approx 0,
\end{align}
where $\mathbf{M}_{m,n}=1/2\delta_{m,n-1}+1/2\delta_{m,n+1}$ is the matrix representing the action of multiplying our basis functions by $y$.
To understand the stability let us now look at the characteristic equations for both recurrences \eqref{eqn:recurrence_1st_moments_omega_larger_n} \& \eqref{eqn:recurrence_2nd_moments_omega_larger_n}: By the definition of $\mathbf{M}T_n=\frac{1}{2}T_{n-1}+\frac{1}{2}T_{n+1}$ the characteristic equations are
\begin{align*}
\left(\frac{1}{4}\left(\lambda+\lambda^{-1}\right)^2-1\right)^2=0
\end{align*}
and
\begin{align*}
\frac{1}{2}\left(\lambda+\lambda^{-1}\right)\left(\frac{1}{4}\left(\lambda+\lambda^{-1}\right)^2-1\right)^2\left[2 \alpha ^2 \beta
^2  -3 \alpha
\beta  \frac{1}{2}\left(\lambda+\lambda^{-1}\right)  +\frac{1}{4}\left(\lambda+\lambda^{-1}\right)^2\right]=0
\end{align*}
respectively. These have solutions
\begin{align*}
\lambda^{(1)}_{1,2,3,4}=-1,\lambda^{(1)}_{5,6,7,8}=1
\end{align*}
and
\begin{align*}
\lambda^{(2)}_{1,2,3,4}&=-1,\lambda^{(2)}_{5,6,7,8}=1,\lambda^{(2)}_{9}=-i,\lambda^{(2)}_{10}=i,\\
\lambda^{(2)}_{11}&= \alpha \beta-\sqrt{\alpha^2 \beta^2-1},\quad \lambda^{(2)}_{12}=\alpha\beta+\sqrt{\alpha^2 \beta^2-1},\\
\lambda^{(2)}_{13}&= 2\alpha \beta-\sqrt{4\alpha^2 \beta^2-1},\quad
\lambda^{(2)}_{14}=2\alpha\beta+\sqrt{4\alpha^2 \beta^2-1}.
\end{align*}
Polynomial stability\footnote{Note, here we mean that the error grows at most polynomially. Although this is not as strong as boundedness in the error (for which we would require the characteristic polynomial to satisfy the root condition) this guarantees sufficient accuracy in the moment computation as we discuss in section \ref{sec:justification_of_stability}.} of the recurrence in forward propagation is equivalent to $|\lambda^{(1)}_i|\leq1$ and $|\lambda^{(2)}_i|\leq1$. For the first case (i.e. to compute the moments $\sigma_n^{(1)}$) this is always satisfied. For the second case (i.e. to compute the moments $\sigma_n^{(2)}$) this is satisfied precisely when $2|ab|\leq 1$.
\subsubsection{Asymptotic behaviour of homogeneous solutions for $n\gg\omega$}\label{sec:asymptotic_behaviour_of_homogeneous solutions_for_n_gg_omega}
Having understood the behaviour of homogeneous solutions to \eqref{eqn:recurrence_for_sigma_1} and \eqref{eqn:recurrence_for_sigma_2} for small n (i.e. when $n\ll\omega$) we now proceed to determine their behaviour when $n\gg\omega$. Note that the coefficients in the recurrences are analytic at $n=\infty$, and in particular have the form $a_0+a_1n^{-1}+a_2n^{-2}$ for some constants $a_0,a_1,a_2$. According to the results by  \cite{birkhoff1933analytic} for such recurrence equations it is possible to find a complete set of homogeneous solutions with known asymptotic behaviour as $n\rightarrow+\infty$. A practical way of determining these is provided by \cite{denef1974asymptotic} and \cite{branders1974asymptotic}. Using those methods we find that there are eight linearly independent homogeneous solutions to \eqref{eqn:recurrence_for_sigma_1} which have the following asymptotic behaviour as $n\rightarrow\infty$:
\begin{align}
\begin{split}\label{eqn:asymptotic_behaviour_homogeneous_solutions_sigma1}
y_n^{(1)}&\sim\left(\frac{i\omega(\beta-1)}{4}\right)^nn^{-n}e^{n}\sim \left(\frac{i\omega(\beta-1)}{4}\right)^n\frac{1}{n!n^{-\frac{1}{2}}},\\
y_n^{(2)}&\sim\left(\frac{i\omega(\beta+1)}{4}\right)^nn^{-n}e^{n}\sim\left(\frac{i\omega(\beta+1)}{4}\right)^n\frac{1}{n!n^{-\frac{1}{2}}},\\
y_n^{(3)}&\sim n^{-4},\quad y_n^{(4)}\sim n^{-2},\quad y_n^{(5)}\sim(-1)^nn^{-2},\quad y_n^{(6)}\sim(-1)^nn^{-2}\log n,\\
y_n^{(7)}&\sim\left(\frac{4i}{\omega(1+\beta)}\right)^nn^{n}e^{-n}\sim\left(\frac{4i}{\omega(1+\beta)}\right)^nn!n^{-\frac{1}{2}},\\
y_n^{(8)}&\sim\left(-\frac{4i}{\omega(1-\beta)}\right)^nn^{n}e^{-n}\sim \left(-\frac{4i}{\omega(1-\beta)}\right)^nn!n^{-\frac{1}{2}}.
\end{split}
\end{align}
\begin{remark}
	In the special case $\beta=0$ these reduce exactly to the results presented by \cite{piessens1983modified} for Hankel transforms.
\end{remark}
Using the same method we find there are fourteen linearly independent homogeneous solutions to \eqref{eqn:recurrence_for_sigma_2} which have the following asymptotic behaviour as $n\rightarrow\infty$:
\begin{align}
\begin{split}\label{eqn:asymptotic_behaviour_homogeneous_solutions_sigma2}
y_n^{(1)}&\sim\left(\frac{i\omega(\beta-1)}{2}\right)^nn^{-n}e^{n}\sim \left(\frac{i\omega(\beta-1)}{2}\right)^n\frac{1}{n!n^{-\frac{1}{2}}},\\
y_n^{(2)}&\sim\left(\frac{i\omega(\beta+1)}{2}\right)^nn^{-n}e^{n}\sim\left(\frac{i\omega(\beta+1)}{2}\right)^n\frac{1}{n!n^{-\frac{1}{2}}},\\
y_n^{(3)}&\sim n^{-3}\lambda_3^n,\quad
y_n^{(4)}\sim n^{-1}\lambda_5^n,\quad
y_n^{(5)}\sim  n^{-1}\lambda_7^n,\quad
y_n^{(6)}\sim n^{-4},\quad
y_n^{(7)}\sim n^{-4}(-1)^n,\\
y_n^{(8)}&\sim n^{-2},\quad
y_n^{(9)}\sim n^{-2}(-1)^n,\quad
y_n^{(10)}\sim n^{-3}\lambda_4^n,\quad
y_n^{(11)}\sim n^{-1}\lambda_6^n,\quad y_n^{(12)}\sim n^{-1}\lambda_8^n ,\\
y_n^{(13)}&\sim\left(\frac{2i}{\omega(1+\beta)}\right)^nn^{n}e^{-n}\sim\left(\frac{2i}{\omega(1+\beta)}\right)^nn!n^{-\frac{1}{2}},\\
y_n^{(14)}&\sim\left(-\frac{2i}{\omega(1-\beta)}\right)^nn^{n}e^{-n}\sim \left(-\frac{2i}{\omega(1-\beta)}\right)^nn!n^{-\frac{1}{2}},
\end{split}
\end{align}
where in the above we denoted $\lambda_3=\alpha  \beta -\sqrt{\alpha ^2 \beta^2-1}, \lambda_4=\alpha  \beta +\sqrt{\alpha ^2 \beta^2-1}$ and
\begin{align*}
\lambda_5&= \alpha\beta-\sqrt{\alpha ^2 \beta ^2-\alpha ^2}-\sqrt{2 \alpha
	^2 \beta ^2-\alpha ^2+\frac{2 \alpha ^3 \beta }{\sqrt{\alpha ^2 \beta ^2-\alpha
			^2}}-\frac{2 \alpha ^3 \beta ^3}{\sqrt{\alpha ^2 \beta ^2-\alpha ^2}}-1},\\
\lambda_6&=\alpha\beta-\sqrt{\alpha ^2 \beta ^2-\alpha ^2}+\sqrt{2
	\alpha ^2 \beta ^2-\alpha ^2+\frac{2 \alpha ^3 \beta }{\sqrt{\alpha ^2 \beta
			^2-\alpha ^2}}-\frac{2 \alpha ^3 \beta ^3}{\sqrt{\alpha ^2 \beta ^2-\alpha
			^2}}-1} ,
		\end{align*}
\begin{align*}
\lambda_7&=\alpha\beta+\sqrt{\alpha ^2 \beta
	^2-\alpha ^2}-\sqrt{2 \alpha ^2 \beta ^2-\alpha ^2-\frac{2 \alpha ^3 \beta
	}{\sqrt{\alpha ^2 \beta ^2-\alpha ^2}}+\frac{2 \alpha ^3 \beta ^3}{\sqrt{\alpha
			^2 \beta ^2-\alpha ^2}}-1} ,\\
\lambda_8&=\alpha\beta+\sqrt{\alpha ^2 \beta ^2-\alpha ^2}+\sqrt{2 \alpha ^2 \beta ^2-\alpha
	^2-\frac{2 \alpha ^3 \beta }{\sqrt{\alpha ^2 \beta ^2-\alpha ^2}}+\frac{2
		\alpha ^3 \beta ^3}{\sqrt{\alpha ^2 \beta ^2-\alpha ^2}}-1}.
\end{align*}
Note there is a choice of the parameters $\alpha,\beta$ for which all of $|\lambda_3|,\dots|\lambda_8|\neq 1$. This tells us that for large $n$ there are at most 5 exponentially growing and 5 exponentially decaying linearly independent homogeneous solutions (and potentially less for certain choices of $\alpha,\beta$). This information allows us to construct a robust way for computing the moments in this setting as well (see section \ref{sec:algorithm_for_stable_moment_computation}).
\subsubsection{Asymptotic behaviour for homogeneous solutions for $n\propto\omega\gg1$}\label{sec:asymptotic_behaviour_solutions_n_propto_omega}
We observed above that for small $n\ll \omega$ the homogeneous solutions are all polynomially behaved, i.e. $\left|y_{n+1}/{y_n}\right|\approx 1$, whilst for $n\gg \omega$ there are some homogeneous solutions of exponential growth leading to exponential instabilities. In this section we aim to understand the cut-off value of $n$ at which this phenomenon occurs (i.e. after which forward propagation is unstable) as a function of $\omega,\alpha,\beta$.
\paragraph{The first recurrence \eqref{eqn:recurrence_for_sigma_1}:} In order to find this cut-off point we begin by letting $C_{\omega,n}=\omega/n$, then the recurrence \eqref{eqn:recurrence_for_sigma_1} is equivalent to
\begin{align}\begin{split}\label{eqn:full_recurrence_with_C_sigma_1}
&\frac{C_{\omega,n}^2(1 - \beta^2)}{4}\sigma_n+\left(2-\frac{1}{n}\right)i\beta C_{\omega,n}\sigma_{n-1}+\left(4-C_{\omega,n}^2(1-\beta^2)+\frac{-8+2i\beta C_{\omega,n}}{n}+\frac{4}{n^2}\right)\sigma_{n-2}\\
&+\left(-6 i \beta C_{\omega,n}+\frac{8+17i\beta C_{\omega,n}}{n}-\frac{16}{n^2}\right)\sigma_{n-3}+\left(-8+\frac{3}{2}C_{\omega,n}^2(1-\beta^2)+\frac{64-4i\beta C_{\omega,n}}{n}-\frac{104}{n^2}\right)\sigma_{n-4}\\
&+\left(6i\beta C_{\omega,n}+\frac{-8-31i\beta C_{\omega,n}}{n}+\frac{48}{n^2}\right)\sigma_{n-5}+\left(4-C_{\omega,n}^2(1-\beta^2)+\frac{-56+2i\beta C_{\omega,n}}{n}+\frac{196}{n^2}\right)\sigma_{n-6}\\
&+\left(-2+\frac{15}{n}\right)i\beta C_{\omega,n}\sigma_{n-7}+\frac{C_{\omega,n}^2(1 - \beta^2)}{4}\sigma_{n-8}=0.
\end{split}
\end{align}

We assume that initially we are in the regime $\omega\gg n\gg 1$ (i.e. initially $C_{\omega,n}\gg 1$), and wish to understand when the initial polynomial behaviour transitions into exponential behaviour as $\omega$ is fixed but $n$ grows. We expect this to occur when $C_{\omega,n}=\Ob(1)$ (if it did not we would simply have to include higher order terms in the following analysis). Thus we expand the ratio of consecutive elements as
\begin{align*}
\frac{y_{n+1}}{y_n}=\lambda_{0}+\Ob\left(n^{-1}\right)
\end{align*}
Now $\lambda_0$ must satisfy
\begin{align}\label{eqn:leading_order_recurrence_with C_sigma_1}
\begin{split}
\frac{C_{\omega,n}^2(1 - \beta^2)}{4}\lambda_0^8+2i\beta C_{\omega,n}\lambda_0^7+\left(4-C_{\omega,n}^2(1-\beta^2)\right)\lambda_0^6+\left(-6 i \beta C_{\omega,n}\right)\lambda_0^5+\left(-8+\frac{3}{2}C_{\omega,n}^2(1-\beta^2)\right)\lambda_0^4&\\
+\left(6i\beta C_{\omega,n}\right)\lambda_0^3+\left(4-C_{\omega,n}^2(1-\beta^2)\right)\lambda_0^2-2i\beta C_{\omega,n}\lambda_0+\frac{C_{\omega,n}^2(1 - \beta^2)}{4}&=0.
\end{split}
\end{align}
\begin{remark}
	The polynomial on the LHS of \eqref{eqn:leading_order_recurrence_with C_sigma_1} is conjugate reciprocal, i.e. such that
	\begin{align*}
	\text{Coefficient\  of\ } \lambda_{0}^{8-j}= \text{Complex\ conjugate\ to\ coefficient\  of\ }\lambda_{0}^{j},\quad j=0,\dots 8.
	\end{align*}
	This means that if $\lambda_0$ solves this equation then so does $1/\lambda_{0}^*$, i.e. the reciprocal of its complex conjugate. Therefore the solutions to the leading order characteristic equation \eqref{eqn:leading_order_recurrence_with C_sigma_1} come in pairs: For any solution with modulus greater than 1 there is precisely one of modulus less than 1. This observation will provide a justification for the stability of our algorithm for computing the moments in section \ref{sec:algorithm_for_stable_moment_computation}.
\end{remark}

We can solve \eqref{eqn:leading_order_recurrence_with C_sigma_1} exactly and find (to leading order in $n$) the following eight solutions, including multiplicities:

\begin{align*}
\left\{-1,-1, 1,1,\frac{2i\pm\sqrt{C_{\omega,n}^2(1-\beta)^2-4}}{ C_{\omega,n}(\beta-1)},\frac{2i\pm\sqrt{C_{\omega,n}^2(1+\beta)^2-4}}{ C_{\omega,n}(\beta+1)}\right\}.
\end{align*}

This tells us the following about the behaviour of $y_{n+1}/y_n$: As long as $C_{\omega,n}$ is large the modulus of each $\lambda_0$ equals 1, so all the homogeneous solutions $y_n$ grow at worst polynomially. However once $C_{\omega,n}=\Ob(1)$ the algebraic behaviour is disturbed and the solutions split into some that are exponentially growing and some that are exponentially decaying. The earliest place this occurs is when
\begin{align*}
C_{\omega,n}\approx \max\left\{\frac{2}{|1-\beta|},\frac{2}{|1+\beta|}\right\}.
\end{align*}
Additionally, there are four solutions, $\lambda_0=-1,-1,1,1$ which have algebraic behaviour even when $C_{\omega,n}=\Ob(1)$.

In summary we have the following behaviour of solutions as visualised in figure \ref{fig:behaviour_of_homogeneous_solutions_sigma_1}.

\begin{figure}[h!]
	\centering
	\begin{tikzpicture}[ scale=0.9]
	\draw[|-|] (0,0) -- (5,0) node[midway, fill=white] {8 algebraic};

	\draw[|-|] (5,0) -- (10,0) node[midway, fill=white] {6 algebraic};
	\draw[|-|] (5,1) -- (10,1) node[midway, fill=white] {1 exponentially growing};
	\draw[|-|] (5,-1) -- (10,-1) node[midway, fill=white] {1 exponentially decaying};

	\draw[|->] (10,0) -- (15,0) node[midway, fill=white] {4 algebraic};
	\draw[|->] (10,1) -- (15,1) node[midway, fill=white] {2 exponentially growing};
	\draw[|->] (10,-1) -- (15,-1) node[midway, fill=white] {2 exponentially decaying};

	\draw[-, dashed] (5,3) -- (5,-2) node[below, fill=white] {$n\approx \omega\min\left\{\frac{|1-\beta|}{2},\frac{|1+\beta|}{2}\right\}$};
	\draw[-, dashed] (10,3) -- (10,-2) node[below, fill=white] {$n\approx \omega\max\left\{\frac{|1-\beta|}{2},\frac{|1+\beta|}{2}\right\}$};

	\draw[color=black] (2.5,2.5) node {$n\ll\omega$};
	\draw[color=black] (7.5,2.5) node {$n\approx\omega$};
	\draw[color=black] (12.5,2.5) node {$n\gg\omega$};

	\draw[-, dashed] (0,2.0) -- (15,2);
	\end{tikzpicture}
	\caption{The growth behaviour of eight linearly independent homogeneous solutions in three\\ regimes for $n,\omega$.}
	\label{fig:behaviour_of_homogeneous_solutions_sigma_1}
\end{figure}

Note in this case it is actually possible to write down the behaviour of the solutions in terms of $C_{\omega, n}$ using the methods given by \cite{denef1974asymptotic} and \cite{branders1974asymptotic}, the corresponding results verify above analysis and are given in appendix B.

\paragraph{The second recurrence \eqref{eqn:recurrence_for_sigma_2}:} We are again interested in the first point when this splitting of the leading order growth ratio occurs. Thus let again $C_{\omega,n}=\omega/n$ and write
\begin{align*}
\frac{y_{n+1}}{y_n}=\lambda_{0}+\Ob\left(n^{-1}\right).
\end{align*}
Then from \eqref{eqn:recurrence_for_sigma_2} we find after a few steps of simplification that $\lambda_{0}$ must satisfy:

\begin{align}
\begin{split}\label{eqn:leading_order_recurrence_with C_sigma_2}
&\left(\lambda_{0} ^2-1\right)^2 \left(\lambda_{0} ^2-2 \alpha  \beta  \lambda_{0} +1\right)
\left(\left(1-\beta ^2\right)
C_{\omega,n}^2 \lambda_{0} ^8+\left(4 \alpha  \beta  \left(\beta ^2-1\right)C_{\omega,n}^2 +4 i \beta C_{\omega,n}\right)\lambda_{0}^7\right.\\
&\left. +\left(4-16 i \alpha  \beta^2 C_{\omega,n}\right) \lambda_{0} ^6+\left(16 i \alpha ^2 \beta  C_{\omega,n}-16 \alpha\beta +4 i
\beta  C_{\omega,n}-4 \alpha  \beta  \left(\beta ^2-1\right) C_{\omega,n}^2 \right) \lambda_{0} ^5\right.\\
&\left. +\left(16 \alpha ^2 +8+2 \left(\beta ^2-1\right) C_{\omega,n}^2 \right)\lambda_{0}
^4+\left(-16 i \alpha ^2 \beta  C_{\omega,n}-4 \alpha  \beta  \left(\beta ^2-1\right) C_{\omega,n}^2 -4 i \beta  C_{\omega,n}-16 \alpha  \beta \right) \lambda_{0} ^3\right.\\
&\left.+\left(4+16 i \alpha  \beta ^2 C_{\omega,n} \right)\lambda_{0} ^2+\left(4\alpha  \beta  \left(\beta ^2-1\right) C_{\omega,n}^2-4 i \beta
C_{\omega,n}\right)\lambda_{0}+(1-\beta^2)C_{\omega,n}^2 \right)=0
\end{split}
\end{align}
As was the case for \eqref{eqn:leading_order_recurrence_with C_sigma_1} we observe that the LHS of \eqref{eqn:leading_order_recurrence_with C_sigma_2} defines a conjugate reciprocal polynomial. This means that if $\lambda_0$ solves this equation then so does $1/\lambda_{0}^*$, i.e. the reciprocal of its complex conjugate. Therefore the solutions to the leading order characteristic equation \eqref{eqn:leading_order_recurrence_with C_sigma_2} come in pairs: For any solution with modulus greater than 1 there is precisely one of modulus less than 1. In particular for any exponentially increasing homogeneous solution in this regime there must be an exponentially decreasing one.

While from \eqref{eqn:leading_order_recurrence_with C_sigma_2} we can spot six solutions (with multiplicities)
\begin{align*}
-1,-1,1,1,\alpha\beta-\sqrt{\alpha^2\beta^2-1},\alpha\beta+\sqrt{\alpha^2\beta^2-1}
\end{align*}
we are still left with a polynomial of degree eight. Instead of finding its exact factors we can use this as a method of checking \textit{approximately} whether for a given $n$ we have reached the regime of exponential growth, and how many solutions have started to exhibit this behaviour.

Indeed for fixed $\omega\gg 1$ there are two cut-off values $N_{min}(\omega)\leq N_{max}(\omega)$ as follows:

\begin{itemize}
	\item In the regime $n\leq N_{min}(\omega)$:
	\begin{itemize}
		\item If $2|\alpha\beta|\leq 1$: All solutions to \eqref{eqn:leading_order_recurrence_with C_sigma_2} are on the complex unit circle and consequently all fourteen homogeneous solutions are algebraically behaved.
		\item If $1<|\alpha\beta|\leq 2$: 12 solutions to \eqref{eqn:leading_order_recurrence_with C_sigma_2} are on the complex unit circle, 1 is inside and 1 outside. Thus there are 12 algebraically behaved, 1 exponentially decaying and 1 exponentially growing solutions.
		\item If $2<|\alpha\beta|$: 10 solutions to \eqref{eqn:leading_order_recurrence_with C_sigma_2} are on the complex unit circle, 2 is inside and 2 outside. Hence there are 10 algebraically behaved and 2 pairs of exponentially growing and decaying solutions.
	\end{itemize}
	\item In the regime $n>N_{max}(\omega)$: For any given choice of $\alpha,\beta$ there are 4/6/8/10/12 algebraically behaved homogeneous solutions and, correspondingly, 5/4/3/2/1 pairs of exponentially growing and decaying solutions.
\end{itemize}

Thus for each fixed choice of $n,\omega$, \eqref{eqn:leading_order_recurrence_with C_sigma_2} is a convenient way of checking approximately in which asymptotic regime $n$ is located.

\subsection{Asymptotic behaviour of the Chebyshev moments}\label{sec:asymptotic_behaviour_chebyshev_moments}

In order to understand how we can stably compute the moments we also need to know their asymptotic behaviour. After a few changes of variable it is possible to write both moments in the form
\begin{align*}
\sigma^{(j)}_n=\frac{(-1)^{n}}{2}\underbrace{\int_{0}^\pi e^{in\theta}\tilde{g}_j(\theta) d\theta}_{=:I_n}+\frac{(-1)^{n}}{2}\underbrace{\int_{0}^\pi e^{-in\theta}\tilde{g}_j(\theta) d\theta}_{=:I_{-n}},\quad j=1,2,
\end{align*}
where
\begin{align*}
\tilde{g}_j(\theta)&=g_j(-\cos\theta)\sin\theta,\quad j=1,2\\
g_1(y)&=\frac{1}{2}H^{(1)}_0(\omega(y+1)/2)e^{i\omega\beta (y+1)/2},\\
g_2(y)&=H_0^{(1)}\left(\omega \sqrt{(y-\alpha\beta)^2+\alpha^2(1-\beta^2)}\right)e^{i\omega\beta y}.
\end{align*}
Note that $g_1$ has a weak singularity at $y=-1$, nevertheless we can change our contours of integration (at least locally) to paths of steepest descent, which shows that
\begin{align*}
I_n&\sim i\int_{0}^{\infty}\tilde{g}_j(it)e^{-nt}dt-i(-1)^n\int_{0}^{\infty}\tilde{g}_j(\pi+it)e^{-nt}dt\\
I_{-n}&\sim -i\int_{0}^{\infty}\tilde{g}_j(-it)e^{-nt}dt+i(-1)^n\int_{0}^{\infty}\tilde{g}_j(\pi-it)e^{-nt}dt
\end{align*}
as $n\rightarrow+\infty$. Now we note the following asymptotic expansions of $\tilde{g}_j$ around $\theta=0,\pi$:
\begin{align*}
\tilde{g}_1(\pm it)&\sim\mp \frac{t \left(2 \log \left(\omega/8\right)+4 \log (t)+ i \pi/2\pm i\pi/4 +2 \gamma
	\right)}{2 \pi }+\Ob(t^3),\quad &t\rightarrow 0^+,\\
\tilde{g}_1(\pi\pm it)&\sim \mp ti \frac{1}{2}H^{(1)}_0(\omega)e^{i\omega \beta}+\Ob(t^3),\quad &t\rightarrow 0^+,\\
\tilde{g}_2(\pm it)&\sim \pm t i H^{(1)}_0(\omega\sqrt{1+2\alpha\beta+\alpha^2})e^{-i\omega \beta} +\Ob(t^3),\quad& t\rightarrow 0^+,\\
\tilde{g}_2(\pi\pm it)&\sim \mp t i H^{(1)}_0(\omega\sqrt{1-2\alpha\beta+\alpha^2})e^{i\omega \beta} +\Ob(t^3),\quad &t\rightarrow 0^+.
\end{align*}
Thus we can apply Laplace's method to $I_{\pm n}$ to find the following asymptotic expansions as $n\rightarrow \infty$:
\begin{align}\label{eqn:asymptotic_behaviour_sigma_1}
\sigma^{(1)}_n&\sim n^{-2}\log n \frac{2i(-1)^{n}}{\pi}+n^{-2}\left(\frac{(-1)^{n+1}}{2}-i(-1)^{n}\frac{\log(\omega/8)+2-\gamma}{\pi}-\frac{1}{2}H^{(1)}_0(\omega)e^{i\omega\beta} \right)+\Ob(n^{-4}),\\\label{eqn:asymptotic_behaviour_sigma_2}
\sigma^{(2)}_n&\sim n^{-2}\left((-1)^{n+1} H_{0}^{(1)}(\omega\sqrt{1+2\alpha\beta+\alpha^2})e^{-i\omega\beta}-H_{0}^{(1)}(\omega\sqrt{1-2\alpha\beta+\alpha^2})e^{i\omega\beta}\right)+\Ob(n^{-4}).
\end{align}

\subsection{Algorithm for the stable computation of the moments $\sigma_n^{(j)}, j=1,2$}\label{sec:algorithm_for_stable_moment_computation}

To begin with let us briefly summarize an observation made by \cite{oliver1968numerical} (see also \cite[chapter 8.3]{wimp1984computation}), concerning the stable solution of recurrence relations. Suppose we are given a generic $m+1$-term recurrence of the form
\begin{align}\label{eqn:generic_recurrence}
a_m(n)y_{n+m}+a_{m-1}(n)y_{n+m-1}+\cdots +a_0(n)y_{n}=0, \quad n\geq 0,
\end{align}
and suppose we are given the values of $y_0, \dots, y_{j-1}$ and $y_{N-m+j+1},\dots, y_{N}$. Solving this discrete boundary value problem is equivalent to the \textit{sparse} matrix system
\begin{align}\label{eqn:BVP_as_matrix_system}
A\mathbf{y}=\mathbf{b}
\end{align}
\begin{align*}
A&=\left[\begin{matrix}
a_{j}(0) & a_{j+1}(0) & \cdots & a_{m}(0) &  & & \\
a_{j-1}(1) & a_{j} (1) & \cdots & a_{m-1}(1) & a_{m}(1) & &\\
\\
& & & \ddots& \ddots &  \ddots & \\
\\
& & & a_{0}(N-m-1) & a_{1}(N-m-1) & \cdots & a_{j+1}(N-m-1)\\
& & & & a_{0}(N-m) & \cdots & a_{j}(N-m)
\end{matrix}\right]\\
\mathbf{y}&=\left[\begin{matrix}
y_{j+1}\\
y_{j+2}\\
\vdots\\
y_{N-m+j}
\end{matrix}\right],\quad \mathbf{b}=\left[\begin{matrix}
-\sum_{l=0}^{j-1}a_l(0)y_{l}\\

\vdots\\
-a_{0}(j-1)y_{j-1}\\
0\\
\vdots\\
0\\
-a_m(N-2m+j+1)y_{N-m+j+1}\\
\vdots\\
-\sum_{l=0}^{m-j-1}a_{m-l}(N-m)y_{N-l}
\end{matrix}\right].
\end{align*}

This matrix system allows us, based on known values of the $(j,m-j)$ boundary conditions, to solve the recurrence \eqref{eqn:generic_recurrence} by LU-factorising \eqref{eqn:BVP_as_matrix_system} using partial pivoting. It is suggested by \cite{oliver1968numerical} that under certain conditions on the recurrence \eqref{eqn:generic_recurrence} this process is numerically stable and we will comment in more detail on theses conditions in section \ref{sec:justification_of_stability}. Let us first however describe our algorithms for computing the moments $\sigma^{(1)}_n,\sigma^{(2)}_n, n=0,\dots, N$ for a fixed $N$ specified in advance. We proceed as follows:

\begin{algorithm}[Computation of the moments $\sigma_{n}^{(1)}$ when $\beta\neq \pm1$]\label{alg:computing_sigma1}\
\begin{enumerate}
	\item Compute the initial conditions $\sigma_0^{(1)},\sigma_1^{(1)},\sigma_2^{(1)},\sigma_3^{(1)}$ using Gauss--Hermite quadrature on the expressions given in corollary \ref{cor:initial_moments_sigma_1}.
	\item If $N\leq \omega\min\{|1-\beta|/2,|1+\beta|/2\}$ then we compute $\sigma_4^{(1)},\dots,\sigma_N^{(1)}$ directly by forward propagation from the recurrence \eqref{eqn:recurrence_for_sigma_1}.
	\item If $N> \omega\min\{|1-\beta|/2,|1+\beta|/2\}$:
	\begin{enumerate}
		\item  We firstly compute $\sigma_4^{(1)},\sigma_5^{(1)}$ directly from \eqref{eqn:recurrence_for_sigma_1}.
		\item Then we pick $M\geq N, M\gg \omega \max\{|1-\beta|/2,|1+\beta|/2\}$ and set $\sigma_{M-1}^{(1)},\sigma_M^{(1)}$ equal to the asymptotic expression given in \eqref{eqn:asymptotic_behaviour_sigma_1}.
		\item Finally, we solve the difference equation \eqref{eqn:recurrence_for_sigma_1} together with 6 initial values $\sigma_0^{(1)},\dots, \sigma_5^{(1)}$ and 2 endpoint values $\sigma_{M-1}^{(1)},\sigma_M^{(1)}$, by solving the sparse matrix system \eqref{eqn:BVP_as_matrix_system} using LU factorisation to yield $\sigma_6^{(1)},\dots,\sigma_{M-2}^{(1)}$.
	\end{enumerate}
	\item From this we output $\sigma_0^{(1)},\dots, \sigma_N^{(1)}$.
\end{enumerate}
\end{algorithm}
	When $\beta=\pm1$ the recurrence \eqref{eqn:recurrence_for_sigma_1} simplifies to a 7-term recurrence. In these cases we can simplify our algorithm slightly:
	\begin{algorithm}[Computation of the moments $\sigma_{n}^{(1)}$ when $\beta=\pm1$]\label{alg:computing_sigma1_beta_pm1}\
		\begin{enumerate}
			\item Compute the initial conditions $\sigma_0^{(1)},\sigma_1^{(1)},\sigma_2^{(1)},\sigma_3^{(1)}$ using Gauss--Hermite quadrature on the expressions given in corollary \ref{cor:initial_moments_sigma_1}.
			\item If $N\leq \omega$ then we compute $\sigma_4^{(1)},\dots,\sigma_N^{(1)}$ directly by forward propagation from the recurrence \eqref{eqn:recurrence_for_sigma_1}.
			\item If $N> \omega$:
			\begin{enumerate}
				\item  We firstly compute $\sigma_4^{(1)},\sigma_5^{(1)}$ directly from \eqref{eqn:recurrence_for_sigma_1}.
				\item Then we pick $M\geq N, M\gg \omega $ and set $\sigma_M^{(1)}$ equal to the asymptotic expression given in \eqref{eqn:asymptotic_behaviour_sigma_1}.
				\item Finally, we solve the difference equation \eqref{eqn:recurrence_for_sigma_1} together with 5 initial values $\sigma_0^{(1)},\dots, \sigma_4^{(1)}$ and 1 endpoint value $\sigma_M^{(1)}$, by solving the sparse matrix system \eqref{eqn:BVP_as_matrix_system} using LU factorisation to yield $\sigma_5^{(1)},\dots,\sigma_{M-1}^{(1)}$.
			\end{enumerate}
			\item From this we output $\sigma_0^{(1)},\dots, \sigma_N^{(1)}$.
		\end{enumerate}
	\end{algorithm}

For the computation of the moments $\sigma_n^{(2)}$ we proceed in a very similar manner:
\begin{algorithm}[Computation of the moments $\sigma_{n}^{(2)}$]\label{alg:computing_sigma2}\
\begin{enumerate}
	\item Compute the initial conditions $\sigma_0^{(2)},\sigma_1^{(2)},\sigma_2^{(2)},\sigma_3^{(2)},\sigma_4^{(2)},\sigma_5^{(2)},\sigma_6^{(2)}$ using Gauss--Laguerre and Gauss--Hermite quadrature on \eqref{eqn:NSD_expressions_for_sigma2}.
	\item Numerically compute the eight non-trivial solutions $\lambda_0^{(1)},\dots,\lambda_0^{(8)}$ of the leading order characteristic equation \eqref{eqn:leading_order_recurrence_with C_sigma_2} by solving
	\begin{align*}
	&\left(1-\beta ^2\right)
	\frac{\omega^2}{N^2} \lambda_{0} ^8+\left(4 \alpha  \beta  \left(\beta ^2-1\right)\frac{\omega^2}{N^2}+4 i \beta \frac{\omega}{N}\right)\lambda_{0}^7\\
	& +\left(4-16 i \alpha  \beta^2 \frac{\omega}{N}\right) \lambda_{0} ^6+\left(16 i \alpha ^2 \beta  \frac{\omega}{N}-16 \alpha\beta +4 i
	\beta  \frac{\omega}{N}-4 \alpha  \beta  \left(\beta ^2-1\right) \frac{\omega^2}{N^2}\right) \lambda_{0} ^5\\
	& +\left(16 \alpha ^2 +8+2 \left(\beta ^2-1\right) C_{\omega,n}^2 \right)\lambda_{0}
	^4+\left(-16 i \alpha ^2 \beta  \frac{\omega}{N}-4 \alpha  \beta  \left(\beta ^2-1\right) \frac{\omega^2}{N^2} -4 i \beta  \frac{\omega}{N}-16 \alpha  \beta \right) \lambda_{0} ^3\\
	&+\left(4+16 i \alpha  \beta ^2 \frac{\omega}{N} \right)\lambda_{0} ^2+\left(4\alpha  \beta  \left(\beta ^2-1\right) \frac{\omega^2}{N^2}-4 i \beta
	\frac{\omega}{N}\right)\lambda_{0}+(1-\beta^2)\frac{\omega^2}{N^2}=0.
	\end{align*}
	\item If $\min\left\{\left|\alpha\beta\pm\sqrt{\alpha^2\beta^2-1}\right|,\left|\lambda_0^{(1)}\right|,\dots,\left|\lambda_0^{(8)}\right|\right\}>1-\epsilon$, for some fixed small\footnote{In practice we might choose $\epsilon=10^{-8}$.} $\epsilon$ then we assume that all homogeneous solutions are of polynomial growth/decay, and so we compute $\sigma_7^{(2)},\dots,\sigma_N^{(2)}$ directly from the recurrence \eqref{eqn:recurrence_for_sigma_2}.
	\item If $\min\left\{\left|\alpha\beta\pm\sqrt{\alpha^2\beta^2-1}\right|,\left|\lambda_0^{(1)}\right|,\dots,\left|\lambda_0^{(8)}\right|\right\}\leq 1-\epsilon$:
	\begin{enumerate}
		\item Firstly we compute $\sigma_7^{(2)},\sigma_8^{(2)}$ directly from \eqref{eqn:recurrence_for_sigma_2}.
		\item Then pick $M\gg N$ and set $\sigma_{M-4}^{(2)},\dots,\sigma_{M}^{(2)}$ equal to the asymptotic expression \eqref{eqn:asymptotic_behaviour_sigma_2}.
		\item Finally, we solve the difference equation \eqref{eqn:recurrence_for_sigma_2} together with 9 initial values $\sigma_0^{(2)},\dots, \sigma_8^{(2)}$ and 5 endpoint values $\sigma_{M-4}^{(2)},\dots,\sigma_{M}^{(2)}$, by solving the sparse matrix system \eqref{eqn:BVP_as_matrix_system} using LU factorisation to yield $\sigma_7^{(2)},\dots,\sigma_{M-5}^{(2)}$.
	\end{enumerate}
\end{enumerate}
\end{algorithm}

\subsubsection{Justification of stability}\label{sec:justification_of_stability}

We have already established in section \ref{sec:asymptotic_behaviour_solutions_n_propto_omega} that forward propagation is (at least polynomially) stable when $n$ is less than the cut-off given in our algorithms. We note that in practice polynomial growth in error is acceptable because Chebyshev coefficients of the interpolated amplitude $f$ decay rapidly. In particular, because $f$ is at least $C^2$ for the Filon paradigm and correspondingly for our quadrature methods to make sense, its Chebyshev coefficients must decay at least like $\check{p}_n=\Ob(n^{-2})$. Thus even if the error in the moment computation grows polynomially, in practice it will not influence the accuracy of the overall quadrature scheme. More worrying, however, would be the exponentially unstable behaviour of forward propagation when $n$ is greater than the cut-off. Thus we need to justify why this is avoided by solving the recurrences using $(6,2),(5,1),(9,5)$ boundary conditions respectively.

Let us focus on the case $|\beta|\neq 1$ in \eqref{eqn:recurrence_for_sigma_1}, since the case $|\beta|=1$ for this recurrence and the justification for \eqref{eqn:recurrence_for_sigma_2} work analogously. Suppose we are given the initial and endpoint values for the moments with a small error $\tilde{\sigma}_n^{(1)}=\sigma_n^{(1)}+\epsilon_n, n=0,\dots,5,M-1,M$. Then, assuming the solution of the $(6,2)$ boundary value problem can be found accurately using LU-factorisation, the sequence $\left(\tilde{\sigma}_n^{(1)}-\sigma_n^{(1)}\right)_{n= 1}^M$ is a homogeneous solution to the recurrence and as such can be expanded in the basis of homogeneous solutions whose asymptotic behaviour we found in section \ref{sec:asymptotic_behaviour_of_homogeneous solutions_for_n_gg_omega}:
\begin{align*}
\tilde{\sigma}_n^{(1)}-\sigma_n^{(1)}=\sum_{j=1}^{8}\eta_j y^{(j)}_n, \quad 0\leq n\leq M.
\end{align*}
The coefficients $\eta_j$ can therefore be found from the following linear system:
\begin{align*}
Y\eta=\left(\begin{matrix}
y^{(1)}_0&y^{(2)}_0&\cdots&y^{(8)}_0\\
\vdots&\vdots&\ddots&\vdots\\
y^{(1)}_5&y^{(2)}_5&\cdots&y^{(8)}_5\\
y^{(1)}_{M-1}&y^{(2)}_{M-1}&\cdots&y^{(8)}_{M-1}\\
y^{(1)}_{M}&y^{(2)}_{M}&\cdots&y^{(8)}_{M}\\
\end{matrix}\right)\left(\begin{matrix}
\eta_1\\
\vdots\\
\eta_8
\end{matrix}\right)=\left(\begin{matrix}
\epsilon_0\\
\vdots\\
\epsilon_5\\
\epsilon_{M-1}\\
\epsilon_M
\end{matrix}\right)
\end{align*}
For large $M$ the asymptotic behaviour of $y_M^{(j)}$ is known, and we thus infer that the matrix $Y$ behaves asymptotically like
\begin{align*}
Y\sim\left(\begin{matrix}
A&\star\\
0&B_M
\end{matrix}\right)+\Ob(M^{-2}),\quad M\rightarrow \infty,
\end{align*}
where $A$ is a $6\times6$ non-singular matrix (independent of $M$), and $B_M$ is given by
\begin{align*}
B_M=\left(
\begin{matrix}
\left(\frac{4i(M-1)}{\omega e(1+\beta)}\right)^{M-1}&\left(-\frac{4i(M-1)}{\omega e(1-\beta)}\right)^{M-1}\\
\left(\frac{4iM}{\omega e(1+\beta)}\right)^M&\left(-\frac{4iM}{\omega e(1-\beta)}\right)^M
\end{matrix}\right)
\end{align*}
\begin{align*}
B_M^{-1}=\left(
\begin{matrix}
\left(\frac{\omega e (1+\beta)}{4 i (M-1)}\right)^{M-1}\frac{1+\beta}{2}&\left(\frac{\omega e (1+\beta)}{4 i M}\right)^{M}\frac{1-\beta}{2}\\
\left(\frac{\omega e (1-\beta)}{4 i (M-1)}\right)^{M-1}\frac{1-\beta}{2}&\left(\frac{\omega e (\beta-1)}{4 i M}\right)^{M}\frac{1-\beta}{2}
\end{matrix}\right).
\end{align*}
Therefore to leading order
\begin{align*}
\left(\begin{matrix}
\eta_7\\
\eta_8
\end{matrix}\right)=B_M^{-1}\left(\begin{matrix}\epsilon_{M-1}\\\epsilon_M\end{matrix}\right)\rightarrow 0,\quad M\rightarrow\infty,
\end{align*}
and indeed this decay is spectrally fast in $M$. Thus if we choose $M$ sufficiently large the solution of the $(6,2)$ boundary value problem acts as a projection onto the span of $y^{(1)},\dots,y^{(6)}$, i.e.
\begin{align*}
\left(\begin{matrix}
\eta_1\\
\vdots\\
\eta_6
\end{matrix}\right)\approx A^{-1}\left(\begin{matrix}
\epsilon_0\\
\vdots\\
\epsilon_5
\end{matrix}\right), \quad \left(\begin{matrix}
\eta_7\\
\eta_8
\end{matrix}\right)\approx 0.
\end{align*}
By the asymptotic analysis from section \ref{sec:asymptotic_behaviour_chebyshev_moments} the moments $\sigma_n^{(1)}$ are a linear combination of those six homogeneous solutions, and so our algorithm ensures that (as long as $M\gg N$) the absolute error is bounded by
\begin{align*}
|\tilde{\sigma}^{(1)}_n-\sigma^{(1)}_n|\leq 6\max_{1\leq j\leq 6} |\eta_j y^{(j)}_n|, \quad 0\leq n\leq N\ll M.
\end{align*}
In other words the absolute error is uniformly bounded in $n$. This is also true for the relative error since to leading order $\sigma^{(1)}_n$ behaves similar to the algebraically decaying homogeneous solutions:
\begin{align*}
\max_{1\leq j\leq 6}|y^{(j)}_n/\sigma^{(1)}_n|\lesssim \frac{n^{-2}\log n}{n^{-2}\log n}\lesssim 1,\quad n\rightarrow \infty.
\end{align*}
\begin{remark}
In this justification we implicitly assumed the homogeneous solutions to remain well-behaved relative to the true moments even for moderate values of $n$. This is not immediately obvious from the asymptotic analysis and we therefore provide numerical evidence in the next section to demonstrate that this is indeed the case in practice. In a similar setting  \cite{piessens1983modified} justified the stability of their algorithm for moment computation numerically and the problem that local behaviour cannot be directly inferred from asymptotic behaviour and hence stability often requires numerical justification was mentioned in \cite[p. 233]{keller2007method}.
\end{remark}
\subsubsection{Moment computation in practice}\label{sec:numerical_demonstration_of_stability}
In figure \ref{fig:absolute_errors_moment_computations} we demonstrate how the absolute error of our moment computations performs in practice if we compute the initial conditions to machine accuracy using Gaussian quadrature as described in section \ref{sec:initial_conditions} and determine the endpoint conditions (where applicable) from the asymptotic behaviour of the moments. Of course, forward propagation (blue dotted curves) is computationally cheaper, which is why, if only a small number of moments are required, our algorithms resolve to use this process instead. To ensure accuracy we have determined the cut-off when solutions start to split into exponentially increasing and decreasing ones, which is highlighted using the vertical red dashed lines. This is the point at which our algorithms switch from forward propagation to LU-factorisation of the relevant matrix system.

\begin{figure}[h!]
	\centering
	\begin{subfigure}[h]{0.49\linewidth}
		\centering
		\includegraphics[width=0.95\textwidth]{figure_4a.eps}
		\caption{Absolute error in $\sigma_n^{(1)}$, $\beta=1$, with \\$(8,0)$ and $(6,2)$ boundary conditions.}
		\label{fig:absolute_errors_moment_computations_linear}
	\end{subfigure}
	\hfill
	\begin{subfigure}[h]{0.49\linewidth}
		\centering
		\includegraphics[width=0.95\textwidth]{figure_4b.eps}
		\caption{Absolute error in $\sigma_n^{(2)}$, $\alpha=0.5,\beta=0.5$, with \\$(14,0)$ and $(9,5)$ boundary conditions.}
		\label{fig:absolute_errors_moment_computations_nonlinear}
	\end{subfigure}
	\caption{The absolute errors $|\sigma_{n,approx}^{(j)}-\sigma_{n,true}^{(j)}|, j=1,2,$ for our moment computations without initial perturbations for $0\leq n\leq N=1000$. The vertical red dashed line marks our predicted change of the stability regime.}
	\label{fig:absolute_errors_moment_computations}
\end{figure}

In figure \ref{fig:absolute_errors_moment_computations_linear} we predict the exponential instability to occur at $n\approx \omega=500$, which is indeed the case - forward propagation is stable up until this point, but incurs an exponential error thereafter. In figure \ref{fig:absolute_errors_moment_computations_nonlinear} we predict the exponential instability to already occur at $n\approx 85$ and we notice that indeed beyond this point the forward propagation is exponentially unstable. However,  the $(9,5)$ boundary distribution succeeds in computing all moments to good accuracy.

Now that we have described and analysed stable algorithms for the moment computation, we can combine them with our Filon methods described in section \ref{sec:extended_filon_method_for_Hankel_kernels} and look at their performance in practice.

\section{Numerical examples of the Filon methods}
\label{sec:numerical_examples_filon_methods}
\subsection{Asymptotic order of the method}
\label{sec:numerical_examples_asymptotic_order}
To begin with, let us confirm that our bounds on the asymptotic error in proposition \ref{prop:filon_paradigm_for_alpha_zero} and \ref{prop:filon_paradigm_for_alpha_non_zero} were tight. To see this we look at the absolute error
\begin{align*}
\mathcal{E}^{(1)}_{[s,\nu]}[f]&=\left|I^{(1)}_{\omega,\beta}-\mathcal{Q}^{(1)}_{[s,\nu]}\right|, \\
\mathcal{E}^{(2)}_{[s,\nu]}[f]&=\left|I^{(1)}_{\omega,\alpha,\beta}-\mathcal{Q}^{(2)}_{[s,\nu]}\right|,
\end{align*}
as we increase the frequency $\omega$ and keep $s,\nu$ fixed. In figure \ref{fig:asymptotic_order_our_filon} we show the resulting error for both methods when $\omega\in[1,300]$. Here we chose $f(x)=x\cos x/(1+x^4)$, and we computed the moments using the algorithms in section \ref{sec:algorithm_for_stable_moment_computation}. The black dash-dotted curves correspond to the asymptotic orders $\Ob(\omega^{-k-2})$ and confirm the results from section \ref{sec:filon_paradigm_for_hankel_kernel}. The error of the Filon methods remain small both for large as well as for moderate frequencies $\omega$.
\begin{figure}[h!]
	\centering
	\begin{subfigure}[h]{0.49\linewidth}
		\centering
		\includegraphics[width=0.95\textwidth]{figure_5a.eps}
		\caption{$\mathcal{E}^{(j)}_{[s,\nu]}[f]$, with $\nu=8,\beta=1.5$.}
		\label{fig:asymptotic_order_our_filon_integral1}
	\end{subfigure}
	\begin{subfigure}[h]{0.49\linewidth}
		\centering
		\includegraphics[width=0.95\textwidth]{figure_5b.eps}
		\caption{$\mathcal{E}^{(j)}_{[s,\nu]}[f]$, with $\nu=9,\alpha=0.2,\beta=0.5$.}
	\end{subfigure}
	\caption{The absolute errors $\mathcal{E}^{(j)}_{[s,\nu]}[f], j=1,2$ for our methods for large $\omega$. The black dash-dotted lines correspond to the asymptotic orders $\Ob(\omega^{-s-2})$, predicted in propositions \ref{prop:filon_paradigm_for_alpha_zero} \& \ref{prop:filon_paradigm_for_alpha_non_zero}.}
	\label{fig:asymptotic_order_our_filon}
\end{figure}

\subsection{Direct comparison to previous work}
\label{sec:numerical_examples_comparison_to_previous_work}
Let us now evaluate how our methods perform in practice when compared to previous work. In particular we compare against the composite Filon method described by \cite{dominguez2013filon} and a graded Gauss--Clenshaw--Curtis quadrature. The relevant errors are shown in figures \ref{fig:comparison_against_dominguez_integral1} and \ref{fig:comparison_against_dominguez_integral2}, where in each case we use the same number of 27 function evaluations of $f=x\cos x/(1+x^4)$. We only compare against our method when $s=0$, i.e. when no derivative values are used, since this already suffices to achieve machine accuracy for the first type of integral $I^{(1)}_{\omega,\beta}[f]$. Our method with $s\geq 1$ would be able to achieve the same accuracy with an even smaller number of function evaluations of $f$, and indeed exhibit a rapid decay in error as $\omega$ increases. We observe that the composite Filon method performs significantly better than the graded Gauss--Clenshaw--Curtis, nevertheless both are outperformed by our direct Filon method. In particular we notice that the asymptotic error decays faster for our method than the composite Filon quadrature, whose relative error appears to remain constant or even increase slightly in $\omega$.

\begin{figure}[h!]
	\centering
	\begin{subfigure}[h]{0.49\linewidth}
		\centering
		\includegraphics[width=0.95\textwidth]{figure_6a.eps}
		\caption{Absolute error, $\mathcal{E}^{(1)}_{[0,27]}[f]$ .}
	\end{subfigure}
	\begin{subfigure}[h]{0.49\linewidth}
		\centering
		\includegraphics[width=0.95\textwidth]{figure_6b.eps}
		\caption{Relative error, $\mathcal{E}^{(1)}_{[0,27]}[f]/I^{(1)}_{\omega,\beta}[f]|$.}
	\end{subfigure}
	\caption{The error of our current method (blue) in comparison to the composite Filon method (red) and a graded Gauss--Clenshaw--Curtis quadrature (green) for $\omega\in[1,300]$.}
	\label{fig:comparison_against_dominguez_integral1}
\end{figure}

\begin{figure}[h!]
	\centering
	\begin{subfigure}[h]{0.49\linewidth}
		\centering
		\includegraphics[width=0.95\textwidth]{figure_7a.eps}
		\caption{Absolute error, $\mathcal{E}^{(2)}_{[0,27]}[f]$ .}
	\end{subfigure}
	\begin{subfigure}[h]{0.49\linewidth}
		\centering
		\includegraphics[width=0.95\textwidth]{figure_7b.eps}
		\caption{Relative error, $\mathcal{E}^{(2)}_{[0,27]}[f]/I^{(2)}_{\omega,\beta}[f]|$.}
	\end{subfigure}
	\caption{The error of our current method (blue) in comparison to the composite Filon method (red) and a graded Gauss--Clenshaw--Curtis quadrature (green) for $\omega\in[1,300]$.}
	\label{fig:comparison_against_dominguez_integral2}
\end{figure}
\subsection{Application to wavelet collocation methods}

To give an example of the robustness of our method with respect to the amplitude function $f$ we consider an example from wavelet collocation methods, as considered for instance by \cite{hsaio2002}. The idea of applying a wavelet basis to acoustic scattering problems is used in order to allow for significant compression of the collocation stiffness matrix, and hence significantly reduces the complexity of the method over standard basis functions. According to \cite[p. 301]{hsaio2002} the most difficult part of numerical implementation of such a scheme on a computer is the quadrature of the non-zero entries in the stiffness matrix. In the aforementioned work the computation of these entries is performed using standard singular quadrature and Gauss--Legendre methods (based on whether the integral in question is (near-)singular or not). The entries are given as the application of the single layer potential $\mathcal{S}$ to basis functions on the boundary $\Gamma$, based on a piecewise smooth parametrisation $p:[a,b]\rightarrow\Gamma$, in the form
\begin{align}\label{eqn:collocation_example_single_layer_wavelet}
\mathcal{S}\chi(t)=\frac{i}{4}\int_a^b H_0^{(1)}(\omega |\mathbf{p}(t)-\mathbf{p}(s)|)\chi(\mathbf{p}(s))|\mathbf{p}'(s)|ds, \quad \text{where }t\in [a,b].
\end{align}
In wavelet-based methods, $\chi\circ \mathbf{p}:[a,b]\rightarrow \mathbb{R}$ is chosen to be a wavelet (supported on some partition of $[a,b]$) which through its vanishing moment property facilitates compression of the collocation matrix. Thus we look the evaluate \eqref{eqn:collocation_example_single_layer_wavelet} when $\chi$ is a wavelet with a large number of vanishing moments. This leads us to consider the family of Daubechies' wavelets - the family of wavelets with the highest number of vanishing moments for a given support. The first of those wavelets is the Haar wavelet, which was used in \citep{hsaio2002}. Here we briefly wish to demonstrate how our method can be used to compute these integrals in frequency independent time if $\Gamma$ is a polygonal boundary. We begin by splitting $\mathcal{S}\chi(t)$ into integrals over the edges of the polygon, and have the following cases:
\begin{itemize}
	\item If $\mathbf{p}(t)$ lies on the edge over which we are integrating, we need to compute an integral of the form
	\begin{align*}
	\int_{c}^{d}H_0^{(1)}(\omega|t-s|)\psi(s)ds,
	\end{align*}
	\item If $\mathbf{p}(s)$ lies on another edge than the one over which we are integrating, we need to compute an integral of the form
	\begin{align*}
	\int_{c}^{d}H_0^{(1)}\left(\omega \sqrt{(p_1(t)-p_1(s))^2+(p_2(t)-p_2(s))^2}\right)\psi(s)ds,
	\end{align*}
\end{itemize}
where $\psi$ is (part of) a wavelet, and $c,d$ are the limits of the parametrisation of the current edge. We notice that these integrals correspond (after rescaling and potentially subdividing into multiple integrals of similar kind) to the forms of the integrals $I^{(1)}_{\omega,\beta},I^{(2)}_{\omega,\alpha,\beta}$ plus some non-singular oscillatory integrals. Thus in effect evaluating the collocation entries amounts to computing the integrals considered in the present work.
\begin{figure}[h!]
	\centering
	\begin{subfigure}[t]{0.49\linewidth}
		\centering
		\includegraphics[width=0.935\textwidth]{figure_8a.eps}
		\caption{Daubechies wavelet $f(x)=\psi^{Daub}_{10}(x)$ produced using \texttt{WaveLab850}.}
	\end{subfigure}\hfill
	\begin{subfigure}[t]{0.49\linewidth}
		\centering
		\includegraphics[width=0.99\textwidth]{figure_8b.eps}
		\caption{Absolute error, $\mathcal{E}^{(1)}_{[0,\nu]}[f]$, the black dash-dotted curve corresponds to $\Ob(\omega^{-2})$ .}
		\label{fig:wavelet_approx_error}
	\end{subfigure}
	\caption{The error of our current method when applied to a wavelet amplitude. The amplitude function is twice continuously differentiable, but has no analytic extension.}
	\label{fig:wavelet_approximation}
\end{figure}

To give one simple example in practice let us consider the Daubechies' wavelet $\psi=\psi^{Daub}_{10}(x)$. The wavelet is twice continuously differentiable, however nowhere smooth. As such it does not have an analytic extension, which means numerical steepest descent cannot be applied. Nevertheless according to Proposition \ref{prop:filon_paradigm_for_alpha_non_zero} we should find that our method is convergent with asymptotic order $\Ob(\omega^{-2})$. This is indeed confirmed in figure \ref{fig:wavelet_approximation}. We note the cost of increasing $\nu$ as shown in figure \ref{fig:wavelet_approx_error} grows like $\Ob(\nu\log\nu)$, since we can compute the Chebyshev coefficients of the amplitude function in $\Ob(\nu\log\nu)$ operations using the method described in section \ref{sec:fast_interpolation_CC_points_sigma1} and the computation of Chebyshev moments is performed using algorithm \ref{alg:computing_sigma1} in $\Ob(\nu)$ operations.

Note that our method allows for the inclusion of plane wave oscillations in the basis functions $\chi$, i.e. we can similarly compute the integrals when the basis is of the form $\chi(\mathbf{p})\exp(\omega\hat{\mathbf{a}}\cdot\mathbf{p})$, for some direction $\hat{\mathbf{a}}$. This would facilitate a tailored approach for a basis in high-frequency scattering.

\subsection{Wave scattering on a screen}
\label{sec:numerical_examples_scattering_on_screen}
Finally, let us consider an example from hybrid methods for high-frequency wave scattering. We look at the scattered field from a finite plate $\Gamma$ extending from $(-1,0)$ to $(1,0)$ in $\mathbb{R}^2$. The scattering problem can (as in the previous example) be formulated as the following boundary integral equation:
\begin{align*}
\psi_i(\mathbf{p})=\frac{i}{4}\int_{\Gamma} H_0^{(1)}\left(\omega|\mathbf{p}-\mathbf{q}|\right)\frac{\partial \psi}{\partial n}(\mathbf{q})d\mathbf{q}.
\end{align*}
Similar to the approach used by \cite{parolin2015} and by \cite{gibbs2019fast} we consider a collocation method for solving this integral equation. In the recent work \citep{gibbs2019fast} the entries of the collocation matrix were computed using numerical steepest descent, which lead to a frequency-independent and very efficient scheme. In this example we demonstrate that our Filon-based approach can also be used to assemble the collocation matrix at frequency independent cost. We follow the construction of a hybrid-basis from \cite{hewett2014} and expand $\frac{\partial \psi}{\partial n}(\mathbf{q})$ along the plate as
\begin{align*}
V_0(s;\omega)+\sum_{n=1}^N V_n^+(s)e^{i\omega s}+V_n^-(s)e^{-i\omega s}
\end{align*}
where $s$ is arclength along the plate, $V_0(s;\omega)=2\partial\psi_i/\partial n$ is the geometrical optics approximation and $V_m^\pm$ are splines graded towards the endpoints of the plate, and $s$ is arclength on the plate.

In this setting we follow a graded mesh approach for the non-oscillatory basis components, $V_m^\pm$. In particular we let $V_m^\pm$ be cubic splines, defined on a mesh $t_{n,j}^\pm$ of spline points in the form
\begin{align*}
-1=t_{1,1}^\pm<\cdots<t_{1,J_1}^\pm<t^\pm_{2,1}<\cdots<t_{N+1,1}^\pm=1
\end{align*}
where $t_{n,1}^+=-1+2\sigma^{N_g+1-n}, n=1,\dots, N_g+1$, i.e. these points are exponentially graded towards $s=-1$ (with parameter $\sigma=0.3$ in our numerical examples). Having defined this rough mesh, the points $t_{n,j},j=1,\dots J_n$ are chosen equidistantly in $[t_{n,1},t_{n+1,1})$, where
\begin{align*}
J_n=p-\left\lfloor \frac{N_g+2-n}{N_g+1} p\right\rfloor+1,\quad  n=1,\dots, N_g
\end{align*} for a given parameter $p$. The mesh $\left\{t_{n,j}^-\right\}$ is constructed in the same way, but flipped such that it is graded towards $+1$. This leads to a basis of $N=1+\sum_{j=0}^{N_g} J_n$ degrees of freedom. In previous work \cite{hewett2014} and also \cite{gibbs2019fast} considered a piecewise-polynomial basis on each subgrid $[t_{n,1},t_{n+1,1})$, each of degree $J_n$. They showed both theoretically (theorem 5.1 in \cite{hewett2014}) and in practice that this basis leads to exponentially fast convergence in the maximal polynomial degree $p$ and that, for fixed $p$, the error increases only at an algebraic rate in $\omega$. We will see in the below graphs that in practice the same is observed when cubic splines are used on the mesh $t_{n,j}^\pm$. This means that in order to solve the scattering problem to a fixed error we need only increase $p$ like $\log \omega$ as $\omega\rightarrow\infty$. The collocation points we denote by $\left\{y_l\right\}_{l=1}^M$ which are chosen such that $\left\{t_{n,j}^+\right\}_{n,j}\cup \left\{t_{n,j}^-\right\}_{n,j}\subset \left\{y_l\right\}_{l=1}^M$. In practice the collocation matrix is ill-conditioned, but this can be resolved (as demonstrated by \cite{gibbs2019fast}) by oversampling $M\gg N$. Thus we introduce additional equispaced sampling points between the spline points
\begin{align*}
\left\{t_{n,j}^+\right\}_{n,j}\cup \left\{t_{n,j}^-\right\}_{n,j}.
\end{align*}
In our examples we do this at a rate three, i.e. $M=3N$, so we introduce two new collocation points between every pair of neighbouring spline points. We then have to compute the entries of the collocation matrix as the following integrals
\begin{align}\label{eqn:col_matrix_entry}
\int_{\mathrm{supp}\,V_n^{\pm}} V_{n}^\pm(s)e^{\pm i\omega s}H_0^{(1)}\left(\omega |s-y_l|\right)ds.
\end{align}
In effect this leads to three types of integrals, and corresponding application of suitable quadrature rules:
\begin{itemize}
	\item Non-oscillatory: If $\omega |\mathrm{supp}\,V_{n}^{\pm}|\leq \omega_0$, for some cut-off $\omega_0$, the integral is non-oscillatory, and we evaluate it using standard adaptive quadrature methods. The same is true when considering $V_n^+$ and $y_l\geq \max(\mathrm{supp}\,V_{n}^{+})+\epsilon$ for some tolerance $\epsilon>0$ and when considering $V_n^-$ and $y_l\leq\min(\mathrm{supp}\,V_{n}^{-})$. In both theses cases the relevant integral takes the non-oscillatory form
	\begin{align*}
	\int_\epsilon^b H^{(1)}_0(\omega x)e^{-i\omega x}f(x)dx.
	\end{align*}
	\item Singular and oscillatory: If $y_l\in \mathrm{supp}\,V_n^{\pm}+(-\epsilon,\epsilon)$, for some tolerance $\epsilon>0$, the integral can be split into a singular oscillatory integral in the form of $I^{(1)}_{\tilde{\omega},1}$ and non-oscillatory but singular integral of the form
	\begin{align*}
	\int_{0}^{b}H^{(1)}_0(\omega x)e^{-i\omega x}f(x)dx.
	\end{align*}
	We apply our new method to the former integral and a standard adaptive method to the latter.
	\item Non-singular and oscillatory:  If we consider $V_n^-$ and $y_l\geq \max(\mathrm{supp}\,V_{n}^{+})+\epsilon$ for some tolerance $\epsilon>0$ and when $V_n^+$ and $y_l\leq\min(\mathrm{supp}\,V_{n}^{-})$, we have an oscillatory, but non-singular integral of the form:
	\begin{align*}
	\int_\epsilon^b H^{(1)}_0(\omega x)e^{i\omega x}f(x)dx.
	\end{align*}
	This we evaluate using the non-singular FCC method as described by \cite{dominguez2013filon}.
\end{itemize}

By evaluating our collocation matrix entries in this way we ensure that as $\omega$ increases the integrals can be computed at uniform cost. That this is indeed the case can be observed in figure \ref{fig:wave_scattering_on_screen}. In figure \ref{fig:wave_scattering_on_screen_time} we show the CPU time it took to compute all entries of the matrix to an absolute tolerance of no more than $10^{-6}$. We notice that the time is uniformly bounded and indeed it becomes cheaper to assemble the matrix for larger frequencies. This is partly due to the asymptotic decay in error in the Filon method, and partly because as $\omega$ increases more of the entries in the collocation matrix become oscillatory ($\omega |\mathrm{supp}V_{n}^{\pm}|\geq \omega_0$), so that we can apply our more efficient Filon method instead of the relatively expensive adaptive quadrature methods for larger frequencies. In this example we used $\omega_0=2.0, \epsilon=10^{-4}, N_{g}=2p, p=3,4,5$ and the computations were performed on an Intel Core i5-8500 CPU processor, on a single core (all computations were performed in series without parallelisation).

In figure \ref{fig:wave_scattering_on_screen_L1} we show the relative $L_1$ error of the scattered part of the field on the plate, i.e.
\begin{align*}
\text{Rel. $L_1$ error}=\frac{\|u-u_{true}\|_{L^1(\Gamma)}}{\|u_{true}\|_{L^1(\Gamma)}},
\end{align*}
where $u(s)=\sum_{n=1}^N V_n^+(s)e^{i\omega s}+V_n^-(s)e^{-i\omega s}$ is our approximation to $-\partial\psi/\partial n$ on $\Gamma$ and the corresponding reference value was computed using the collocation method provided by \cite{gibbs2019fast} using $576$ degrees of freedom. We observe that for fixed $p$ the error grows at an algebraic rate in $k$, whilst as we increase $p$ the error decays exponentially, similar to the observations made by \cite{hewett2014} and \cite{gibbs2019fast}. Hence, together with the hybrid basis provided by \cite{hewett2014}, our direct Filon method can indeed be used to solve a 2D scattering problem at near uniform cost even for large frequencies.

\begin{figure}[h!]
	\centering
	\begin{subfigure}[h]{0.49\linewidth}
		\centering
		\includegraphics[width=0.9\textwidth]{cpu_time_fcc_plate.eps}
		\caption{CPU time of collocation matrix assembly.}
		\label{fig:wave_scattering_on_screen_time}
	\end{subfigure}\hfill
	\begin{subfigure}[h]{0.49\linewidth}
		\centering
		\includegraphics[width=0.9\textwidth]{relative_error_col_field_plate.eps}
		\caption{Relative $L_1$-error on the screen.}
		\label{fig:wave_scattering_on_screen_L1}
	\end{subfigure}
	\caption{The application of our quadrature routine to a collocation method for wave scattering on a screen.}
	\label{fig:wave_scattering_on_screen}
	\vspace{-0.4cm}
\end{figure}
\newpage To conclude this example, we highlight the conceptual difference between our current Filon-based method, and the numerical steepest descent method used in \cite{gibbs2019fast}. The latter method is highly efficient but requires some knowledge and care in order to split the domain of integration into subintervals on which the basis functions have analytic extensions. Our Filon-based method provides additional robustness by avoiding this requirement of analyticity, instead replacing it by simple twice continuous differentiability. For cubic splines this means no splitting of the domain of integration is required. This, of course, comes at a cost: in particular our method (when $s=0$) achieves an asymptotic error of size no larger than $\Ob_\delta(\omega^{-2+\delta})$ for any $\delta>0$, whilst numerical steepest descent incurs significantly smaller asymptotic errors of $\Ob(\omega^{-n})$, where $n$ is proportional to the number of quadrature points used. Nevertheless, our method provides a frequency-independent quadrature method for computing the entries of the collocation matrix.

\section{Concluding remarks}
\label{sec:conclusions}
In this work we present a direct version of a Filon--Clenshaw--Curtis method that can be applied to integrals arising in 2D wave-scattering problems. These integrals are generally difficult to compute since they involve the 2D Green's function of the Helmholtz equation, which has both high-frequency oscillations and (frequency-dependent) singularities. The singularities prevent an application of classical Filon methods, and in particular require the development of quadrature techniques that are tailored to the kernel functions at hand. In the setting of Filon methods a particular challenge arises from the need to compute moments accurately and efficiently. In this work we achieve this stable moment computation using a simple recurrence that is found via duality to a spectral method applied to the kernel function. Exploiting information about the asymptotic behaviour of solutions to this recurrence allows us to avoid exponential instabilities in this process, essentially by projecting onto the space of algebraically behaved solutions. The advantages of our approach include a rapid asymptotic decay in error and stability with respect to the regularity of the amplitude function. These properties are demonstrated in a number of practical examples which include an application to wavelet-based compression for wave scattering problems, and to hybrid numerical-asymptotic methods on screens and polygons in 2D.

The moment computation by duality as described in theorem \ref{thm:recurrence_by_duality} generalises a corresponding result for Chebyshev basis functions provided by \cite{keller2007method} and could provide a basis for future work on the moment computation of Filon methods in higher dimensions. Additionally, the same principles as in our current work can be applied to any kernel function satisfying a differential equation with polynomial coefficients, which means they could also be exploited to compute integrals over a wide range of special functions. A particular problem of interest, that is amongst these kernel functions, is the moment computation of the corresponding Filon method for wave scattering problems on more complicated shapes and curves.

\section*{Acknowledgements}
The authors would like to thank  Alfredo Dea\~no (University of Kent), Andrew Gibbs (University College London), Daan Huybrechs (KU Leuven), Anastasia Kisil (University of Manchester) and Sheehan Olver (Imperial College London) for several interesting discussions about highly oscillatory quadrature and special functions. We also thank Victor Domínguez (Universidad P\'ublica de Navarra) for sharing the Matlab implementation for his Filon method with us. Finally, the authors gratefully acknowledge support from the UK Engineering and Physical Sciences Research Council (EPSRC) grant EP/L016516/1 for the University of Cambridge Centre for Doctoral Training, the Cambridge Centre for Analysis.

\bibliography{IMANUM-refs}

\appendix

\section*{Appendix A.\ Fast interpolation at Clenshaw--Curtis points, mid- and endpoint derivatives}\label{app:proof_of_fast_interpolation_at_CC_and_midpoints}
Here we provide a brief derivation of the results which are given in section \ref{sec:fast_interpolation_sigma2}. Recall that we wish to solve the following interpolation problem:
\begin{align}\label{eqn:interpolation_problem_p2_appendix}
p_2^{(j)}(0)=f^{(j)}(0), p_2^{(j)}(\pm1)=f^{(j)}(\pm1), j=0,\dots,s\,\,\,\,\text{\ and\ }\,\,\,\,p_2(c_n)=f(c_n), n=1,\dots,\nu,
\end{align}
using an expansion in Chebyshev polynomials:
\begin{align*}
p_2(x)=\sum_{n=0}^{\nu+3s+1}p_nT_n(x).
\end{align*}
Let us adopt the notation used by \cite{Gao2017a} and define
\begin{align*}
\hat{p}_0=2p_0, \quad \hat{p}_k=p_k, \quad k=1,\dots,\nu,\quad \hat{p}_{\nu+1}=2p_{\nu+1}\\
h_j=f_j-\sum_{m=\nu+2}^{\nu+3s}p_m\cos\left(\frac{jm\pi}{\nu+1}\right),\quad  f_j=f\left(\cos\frac{j\pi}{\nu+1}\right),\quad j=0,\dots,\nu+1
\end{align*}
Then the interpolation conditions $p_2(c_n)=f(c_n), n=0,\dots,\nu+1$ are equivalent to saying that $\mathcal{C}_{\nu+1}\hat{\mathbf{p}}=\mathbf{h}$, where $\mathcal{C}_{\nu+1}$ is the DCT-I. The inverse is
\begin{align}\label{eqn:inverse_DCT-I_for_pmhat}
\hat{p}_m=\left(\mathcal{C}_{\nu+1}^{-1}\mathbf{h}\right)_m=\frac{2}{\nu+1}\sideset{}{''}\sum_{j=0}^{\nu+1}h_j \cos \left(\frac{jm\pi}{\nu+1} \right) \quad \text{for\ }m = 0, \dots, \nu+1.
\end{align}
Note now
\begin{align}\nonumber
h_j&=f_j-\sum_{m=\nu+2}^{\nu+3s+1}\cos\left(\frac{jm\pi}{\nu+1}\right)p_m\\\label{eqn:simplification_hj}
&=f_j-\sum_{m=1}^{3s}(-1)^j\cos\left(\frac{jm\pi}{\nu+1}\right)p_{\nu+1+m},\quad j=0,\dots,\nu+1
\end{align}
Using \eqref{eqn:simplification_hj} in \eqref{eqn:inverse_DCT-I_for_pmhat} we find for $m=0,\dots, \nu+1$ (and with $\check{\mathbf{p}}=\mathcal{C}_{\nu+1}^{-1}\mathbf{f}$):
\begin{align*}
\hat{p}_m&=\check{p}_m-\frac{2}{\nu+1}\sum_{n=1}^{3s}p_{\nu+1+n}\left[\sideset{}{''}\sum_{j=0}^{\nu+1}(-1)^j\cos\left(\frac{jn\pi}{\nu+1}\right) \cos \left(\frac{jm\pi}{\nu+1} \right) \right]\\
&=\check{p}_m-\frac{1}{\nu+1}\sum_{n=1}^{3s}p_{\nu+1+n}\left[\sideset{}{''}\sum_{j=0}^{\nu+1}(-1)^j\cos\left(\frac{j(n+m)\pi}{\nu+1}\right)+ \sideset{}{''}\sum_{j=0}^{\nu+1}(-1)^j\cos\left(\frac{j(n-m)\pi}{\nu+1}\right) \right].
\end{align*}
Now because $\nu$ is odd, one can quickly check using standard trigonometric identities:
\begin{align*}
\sideset{}{''}\sum_{j=0}^{\nu+1}(-1)^j\cos\left(\frac{j(n+m)\pi}{\nu+1}\right)=\begin{cases}
0,&n+m\neq \nu+1\\
\nu+1,&n+m=\nu+1
\end{cases}.
\end{align*}
Thus we find
\begin{align*}
\hat{p}_m&=\check{p}_m-\sum_{n=1}^{3s}p_{\nu+1+n}(\delta_{n+m,\nu+1}+\delta_{n-m,\nu+1})
\end{align*}
which implies
\begin{align*}
p_n&=\frac{1}{2}\check{p}_n\quad\quad n=0,\nu+1\\
p_n&=\check{p}_n\quad\quad n=1,\dots,\nu-3s\\
p_n&=\check{p}_n-p_{2\nu-n+2}\quad\quad n=\nu-3s+1,\dots,\nu
\end{align*}
The remaining interpolation conditions $p_2^{(j)}(0)=f^{(j)}(0), p_2^{(j)}(\pm1)=f^{(j)}(\pm1),j=1,\dots,s$ are equivalent to
\begin{align*}
\sum_{n=1}^{3s}p_{\nu+1+n}\left[T_{\nu+1+n}^{(j)}(-1)-T_{\nu+1-n}^{(j)}(-1)\right]&=f^{(j)}(-1)-\sideset{}{''}\sum_{n=0}^{\nu+1}\check{p}_nT_n^{(j)}(-1)\\
\sum_{n=1}^{3s}p_{\nu+1+n}\left[T_{\nu+1+n}^{(j)}(0)-T_{\nu+1-n}^{(j)}(0)\right]&=f^{(j)}(0)-\sideset{}{''}\sum_{n=0}^{\nu+1}\check{p}_nT_n^{(j)}(0)\\
\sum_{n=1}^{3s}p_{\nu+1+n}\left[T_{\nu+1+n}^{(j)}(1)-T_{\nu+1-n}^{(j)}(1)\right]&=f^{(j)}(1)-\sideset{}{''}\sum_{n=0}^{\nu+1}\check{p}_nT_n^{(j)}(1)
\end{align*}
where the values of $T_n^{(j)}(\pm1)$ are given by
\begin{align*}
T_n^{(j)}(\pm1)&=(\pm 1)^{n-j}\frac{2^j j! n(n+j-1)!}{(2j)!(n-j)!}, \quad 0\leq j\leq n,n+j\geq 1
\end{align*}
and the values of $T_{n}^{(j)}(0), j=1,\dots,s$ can be efficiently computed from the recurrence
\begin{align*}
T_{n+1}^{(j)}(0)&=2T_n^{(j-1)}(0)-T_{n-1}^{(j)}(0), \quad T_{0}(x)=1,T_{1}(x)=x.
\end{align*}

\section*{Appendix B.\ Detailed asymptotic behaviour of the homogeneous solutions to \eqref{eqn:recurrence_for_sigma_1} when $\omega/n=\Ob(1)$}\label{app:detailed_asymptotic_behaviour_intermediate_regime_sigma1}
If we consider \eqref{eqn:full_recurrence_with_C_sigma_1} we can again look at the asymptotic behaviour of solutions and we notice that if $C_{\omega,n}=\omega/n=\Ob(1)$, we have homogeneous solutions with the following asymptotic behaviour as $n\rightarrow\infty$:

\begin{align}
\begin{split}\label{eqn:homogeneous_solutions_intermediate_regime}
y_n^{(1)}&\sim n^{\gamma_1}\left(\frac{2i-i\sqrt{4-(\beta-1)^2C_{\omega,n}^2}}{(\beta-1)C_{\omega,n}}\right)^n,\quad
y_n^{(2)}\sim n^{\gamma_2}\left(\frac{2i-i\sqrt{4-(\beta+1)^2C_{\omega,n} ^2}}{(\beta+1)C_{\omega,n}}\right)^n,\\
y_n^{(3)}&\sim n^{-4}, \quad y_n^{(4)}\sim n^{-2},\quad y_n^{(5)}\sim(-1)^nn^{-2},\quad y_n^{(6)}\sim(-1)^nn^{-2}\log n,\\
y_n^{(7)}&\sim n^{\gamma_7}\left(\frac{2i+i\sqrt{4-(\beta+1)^2C_{\omega,n} ^2}}{(\beta+1)C_{\omega,n}}\right)^n,\quad
y_n^{(8)}\sim n^{\gamma_8}\left(\frac{2i+i\sqrt{4-(\beta-1)^2C_{\omega,n}^2}}{(\beta-1)C_{\omega,n}}\right)^n,
\end{split}
\end{align}
where
\begin{align*}
\gamma_1&=\frac{1}{2} \left(-\frac{(\beta -1) C_{\omega,n}}{\sqrt{(\beta -1)^2
		C_{\omega,n}^2-4}}+\frac{16 i}{\sqrt{(\beta -1)^2 C_{\omega,n}^2-4}}+\beta
+6\right)\\
\gamma_2&=\frac{1}{4} \left(-\frac{(1+8 i) (\beta  C_{\omega,n}+C_{\omega,n}-2)}{\sqrt{(\beta
		+1)^2 C_{\omega,n}^2-4}}-\frac{(1-8 i) (\beta
	C_{\omega,n}+C_{\omega,n}+2)}{\sqrt{(\beta +1)^2 C_{\omega,n}^2-4}}-2 \beta
+12\right)\\
\gamma_7&=\frac{1}{4} \left(\frac{(1+8 i) (\beta  C_{\omega,n}+C_{\omega,n}-2)}{\sqrt{(\beta
		+1)^2 C_{\omega,n}^2-4}}+\frac{(1-8 i) (\beta
	C_{\omega,n}+C_{\omega,n}+2)}{\sqrt{(\beta +1)^2 C_{\omega,n}^2-4}}-2 \beta
+12\right)\\
\gamma_8&=\frac{1}{2} \left(\frac{(\beta -1) C_{\omega,n}}{\sqrt{(\beta -1)^2
		C_{\omega,n}^2-4}}-\frac{16 i}{\sqrt{(\beta -1)^2 C_{\omega,n}^2-4}}+\beta
+6\right).
\end{align*}
Note that as we have fixed the ratio $C_{\omega,n}$ for this derivation we cannot in general expect to recover the precise behaviour of solutions of the original difference equation \eqref{eqn:recurrence_for_sigma_1} (since we have absorbed some $1/n$ terms into $C_{\omega,n}$). Nevertheless this provides a way of looking at the properties of the homogeneous solutions for the case when $\omega,n$ are of the same order of magnitude. In particular looking at $y_n^{(j)}, j=1,2,7,8$ we observe the following split in the growth rates of the homogeneous solutions
\begin{itemize}
	\item If $C_{\omega,n}>\max\left\{|2/(\beta+1)|,|2/(\beta-1)|\right\}$ all growth ratios are of modulus 1, i.e.
	\begin{align*}
	\left|\frac{y_{n+1}^{(j)}}{y_{n}^{(j)}}\right|\approx 1,\quad j=1,\dots, 8.
	\end{align*}
	\item If $\min\left\{|2/(\beta+1)|,|2/(\beta-1)|\right\}< C_{\omega,n} < \max\left\{|2/(\beta+1)|,|2/(\beta-1)|\right\}$ there is one exponentially decreasing, and one exponentially decreasing one, depending on the sign of $\beta$: If $\beta<0$:
	\begin{align*}
	\left|\frac{y_{n+1}^{(8)}}{y_{n}^{(8)}}\right|&>1,\quad\left|\frac{y_{n+1}^{(1)}}{y_{n}^{(1)}}\right|<1,\quad \left|\frac{y_{n+1}^{(j)}}{y_{n}^{(j)}}\right|\approx 1,\quad j=2,\dots, 7,
	\end{align*}
	and if $\beta>0$ similar applies to $y_n^{(7)},y_n^{(2)}$.
	\item If $C_{\omega,n}<\min\left\{|2/(\beta+1)|,|2/(\beta-1)|\right\}$ then there are two exponentially increasing and two exponentially decreasing solutions:
	\begin{align*}
	\left|\frac{y_{n+1}^{(8)}}{y_{n}^{(8)}}\right|,\left|\frac{y_{n+1}^{(7)}}{y_{n}^{(7)}}\right|&>1,\quad\left|\frac{y_{n+1}^{(1)}}{y_{n}^{(1)}}\right|,\left|\frac{y_{n+1}^{(2)}}{y_{n}^{(2)}}\right|<1,\quad \left|\frac{y_{n+1}^{(j)}}{y_{n}^{(j)}}\right|\approx 1,\quad j=3,\dots, 6.
	\end{align*}
\end{itemize}
This confirms the result on the change of the growth rates occurs at $\min\left\{|2/(\beta+1)|,|2/(\beta-1)|\right\}$ that we derived in section \ref{sec:asymptotic_behaviour_solutions_n_propto_omega}.

\newpage

\end{document}

config: arxiv