TITLE: Colon property of Gorenstein rings
QUESTION [4 upvotes]: Let $(R, \mathfrak{m})$ be a Gorenstein local ring of characteistic $p>0$. Let $x_1,...,x_d$ be a system of parameters of $R$. Let $I$ be an $\mathfrak{m}$-primary ideal containg $(x_1,...,x_d)$. Then there exists an ideal $J$ containng $(x_1,...,x_d)$ such that
$(x_1,...,x_d):J = I$.
Question. Is it true that for all power of $p$, say $q$, we have
$$(x_1^q,...,x_d^q):J^{[q]} = I^{[q]}?$$

REPLY [2 votes]: We prove an characterization of regularity based on the result of Wenliang Zhang (arxiv.org/abs/0709.0943). It should be metion his result
Theorem 1 (W. Zhang). Let $(R, \mathfrak{m})$ be a local ring of characteristic $p>0$. The following are equivalent.

$R$ is regular
$I^{[q]}:J^{[q]} = (I:J)^{[q]}$ for all pair of ideals $I, J$ and for all powers $q$ of $p$.
$I^{[q]}:x^{[q]} = (I:x)^{[q]}$ for all $\mathfrak{m}$-primary ideals $I$, all elements $x\in R$, and for all powers $q$ of $p$.

Definition 2. Let $R$ be a Notherian ring of characteristic $p$, $I$ an ideal. Then

The tight closure of $I$ is $I^{*} = \{ x :\exists c \in R^o, cx^q \in I^{[q]} \,\, \text{for all} \,\, q =p^e \gg 0\}$, where $R^o = R \setminus \cup_{\mathfrak{p}\in min(R)} \mathfrak{p}$.
The Frobenius closure of $I$ is $I^F = \{ x: x^q \in I^{[q]}\,\, \text{for all}\,\,  q =p^e \gg 0\}$

Remark 3. It is easy to see that $I \subset I^F \subset I^*$. If $I^{[q]}:x^{[q]} = (I:x)^{[q]}$ for all elements $x\in R$, and for all powers $q$ of $p$, then $I = I^*$ (see Zhang's paper, Lemma 3.1).
Lemma 4. Let $I$ and $J$ is Frobenius closed i.e. $I = I^F$ and $J=J^F$. Then $(I \cap J)^{[q]} = I^{[q]} \cap J^{[q]}$ for all $q$.
proof. easy.
Result 5. Let $(R, \mathfrak{m})$ be a local ring of characteristic $p>0$. The following are equivalent.

$R$ is regular
(a) $R$ is Gorenstein;
(b) $\mathfrak{a}^{[q]}:x^{[q]} = (\mathfrak{a}:x)^{[q]}$ for all parameter ideals $\mathfrak{a}$, all elements $x\in R$, and for all powers $q$ of $p$.

Proof. $(1) \to (2)$ is clear.
$(2) \to (1)$, By Remark 3 we have every parameter ideal of $R$ is tight closed. Moreover $R$ is Gorenstein so every ideal of $R$ is tight closed and hence Frobenius closed. It is enough to prove (3) of Zhang's theorem. Let $I$ be an $\mathfrak{m}$-primary ideal and $x$ an element. Choose a parameter ideal $\mathfrak{a}$ contained in $I$. Since $R$ is Gorenstein we have an ideal $\mathfrak{b} = (b_1,...,b_t)$ such that $I = \mathfrak{a}:\mathfrak{b} = \cap_{i=1}^t(\mathfrak{a}:b_i)$. By Lemma 4 we have
$$I^{[q]} = \cap_{i=1}^t(\mathfrak{a}:b_i)^{[q]} = \cap_{i=1}^t(\mathfrak{a}^{[q]}:b_i^q).$$
Therefore
$$I^{[q]}:x^q = \big(\cap_{i=1}^t(\mathfrak{a}^{[q]}:b_i^q) \big):x^q = \cap_{i=1}^t(\mathfrak{a}^{[q]}:(xb_i)^q) = \cap_{i=1}^t(\mathfrak{a}:(xb_i))^{[q]}.$$
The last equation follows from our asumption. So by Lemma 4 we have
$$I^{[q]}:x^q = \big(\cap_{i=1}^t(\mathfrak{a}:(xb_i)) \big)^{[q]} = (I:x)^{[q]}.$$
The proof is complete.
Remark 6. I do not know the charateriztion of rings satisfying the property 2(b) of Result 5.


text:
\begin{document}

\title{On fractional GJMS operators}
\author{Jeffrey S. Case}
\thanks{JSC was partially supported by NSF Grant No.\ DMS-1004394}
\address{Department of Mathematics \\ Princeton University \\ Fine Hall, Washington Road \\ Princeton, NJ 08544}
\email{jscase@math.princeton.edu}
\author{Sun-Yung Alice Chang}
\thanks{SYAC was partially supported by NSF Grant No.\ DMS-1104536}
\address{Department of Mathematics \\ Princeton University \\ Fine Hall, Washington Road \\ Princeton, NJ 08544}
\email{chang@math.princeton.edu}
\date{\today}
\keywords{fractional Laplacian; Poincar\'e--Einstein manifold; conformally covariant operator; smooth metric measure space}
\subjclass[2010]{Primary 58J70; Secondary 35J70, 53A30}
\begin{abstract}
We describe a new interpretation of the fractional GJMS operators as generalized Dirichlet-to-Neumann operators associated to weighted GJMS operators on naturally associated smooth metric measure spaces.  This gives a geometric interpretation of the Caffarelli--Silvestre extension for $(-\Delta)^\gamma$ when $\gamma\in(0,1)$, and both a geometric interpretation and a curved analogue of the higher order extension found by R.\ Yang for $(-\Delta)^\gamma$ when $\gamma>1$.  We give three applications of this correspondence.  First, we exhibit some energy identities for the fractional GJMS operators in terms of energies in the compactified Poincar\'e--Einstein manifold, including an interpretation as a renormalized energy.  Second, for $\gamma\in(1,2)$, we show that if the scalar curvature and the fractional $Q$-curvature $Q_{2\gamma}$ of the boundary are nonnegative, then the fractional GJMS operator $P_{2\gamma}$ is nonnegative.  Third, by assuming additionally that $Q_{2\gamma}$ is not identically zero, we show that $P_{2\gamma}$ satisfies a strong maximum principle.
\end{abstract}
\maketitle

\input{intro}

\subsection*{Acknowledgments} The authors would like to thank Mar\'ia del Mar Gonz\'alez, Robin Graham, Jie Qing, Paul Yang, and Ray Yang for helpful conversations while this paper was being put together.

\input{scattering}
\input{smms}
\input{extension}
\input{energy}
\input{adapted}
\input{positivity}

\appendix
\input{product}

\bibliographystyle{abbrv}
\bibliography{../bib}
\end{document}

config: arxiv


TITLE: Let $S^2$ denote the unit sphere in $\mathbb{R}^3$ and $e^3=(0,0,1)$. Show that $X:= S^2\setminus \{e_3, -e_3\}$ is not compact.
QUESTION [0 upvotes]: Let $S^2$ denote the unit sphere in $\mathbb{R}^3$ and $e^3=(0,0,1)$. Show that $X:= S^2\setminus \{e_3, -e_3\}$ is not compact by exhibiting an open cover which does not admit a finite subcover.
My attempt: One can project down the sphere on the plane via sterographic projection.
But I cannot procced further. Please help. Thanks.

REPLY [1 votes]: More generally if $X$ is a Hausdorff space and $x\in X$ is not isolated then $X\backslash \{x\}$ is not compact. This is simply because there is a net $(x_\alpha)$ convergent to $x$ totally contained in $X\backslash \{x\}$. And thus it has no convergent subnet in $X\backslash \{x\}$, by uniqueness of limits in Hausdorff spaces.
The above should be applied to for example $X=S^2\backslash\{e_3\}$ and $x=-e_3$.


text:
\section{More examples}\label{sec:misc}
We discuss a few more families of cyclic branched covers of knots.

\subsection{Knot epimorphisms}
Other examples of knots $K$ with $\pi_1(\Sigma_n(K))$ left-orderable arise from the following easy lemma.

\begin{lemma}\label{lem:knotepi}
Let $K, K'$ be prime knots in $S^3$.  Suppose there exists an epimorphism $\varphi: \pi_1(S^3 - K) \to \pi_1(S^3 - K')$ such that $\varphi(\mu) = \mu'$ for meridians $\mu, \mu'$ of $K, K'$ respectively.  If $\pi_1(\Sigma_n(K'))$ is left-orderable, then so is $\pi_1(\Sigma_n(K))$.
\end{lemma}
\begin{proof}
Since there is an exact sequence
\[
1 \to \pi_1(\Sigma_n(K)) \to \pi_1(S^3 - K)/\langle \langle \mu^n \rangle \rangle \to \Z/n \to 1,
\]
and similarly for $K'$, we have that $\varphi$ induces an epimorphism from $\pi_1(\Sigma_n(K))$ to $\pi_1(\Sigma_n(K'))$.  The result now follows from Theorem~\ref{thm:brw}.
\end{proof}

Examples to which Lemma~\ref{lem:knotepi} applies are given in \cite{KS1,KS2}, which list all pairs of distinct prime knots $K, K'$ with at most 10 crossings such that there exists an epimorphism $\pi_1(S^3 - K) \to \pi_1(S^3 - K')$ (written $K \geq K'$ in the notation of \cite{KS1,KS2}).  The authors note in \cite{KS2} that in all cases, there is actually a meridian-preserving epimorphism.  Also, the knot $K'$ is either $3_1$, $4_1$, or $5_2$.  It follows from Theorem~\ref{thm:torus} that for the knots with $K \geq 3_1$, we have $\pi_1(\Sigma_n(K))$ is left-orderable for $n \geq 6$.  It also follows from \cite{Hu2013} (see also \cite{Tran2013}) that for the knots $K$ with $K \geq 5_2$, we have $\pi_1(\Sigma_n(K))$ is left-orderable for $n \geq 9$.  Since all cyclic branched covers of $4_1$ have non-left-orderable fundamental group \cite{DPT2005}, we are not able to say anything about $\pi_1(\Sigma_n(K))$ for $K \geq 4_1$.

\subsection{Unknotting number one, determinant one knots}

\begin{proposition}
Let $K$ be an unknotting number one knot with $det(K) = 1$.  If the Li-Roberts Conjecture is true, then $\Sigma_2(K)$ is an excellent manifold.
\end{proposition}
\begin{proof}
Suppose that $K$ is an unknotting number one knot with $det(K) = 1$.  Since $K$ has unknotting number 1, it follows that $\Sigma_2(K)$ is given by a half-integral surgery on a knot $K'$ in $S^3$.  Since $det(K) = 1$, $H_1(\Sigma_2(K)) = 0$ and thus $\Sigma_2(K) = S^3_{\pm \frac{1}{2}}(K')$.  The result now follows from Conjecture~\ref{conj:lr} and Lemma~\ref{lem:toroidalzs3}.
\end{proof}


text:
\section{Equivalence of Definitions of Ceiling Function}
Tags: Ceiling Function

\begin{theorem}
Let $x$ be a [[Definition:Real Number|real number]].
{{TFAE|def = Ceiling Function}}
\end{theorem}

\begin{proof}
=== Definition 1 equals Definition 2 ===
Follows from [[Infimum of Set of Integers equals Smallest Element]].
{{qed|lemma}}
=== Definition 1 equals Definition 3 ===
Let $S$ be the [[Definition:Set|set]]:
:$S = \left\{ {m \in \Z: m \ge x}\right\}$
Let $n = \inf S$.
By [[Infimum of Set of Integers is Integer]], $n \in \Z$.
By [[Infimum of Set of Integers equals Smallest Element]], $n\in S$.
Because $n\in S$, we have $n \geq x$.
Because $n-1 < n$, we have by definition of [[Definition:Infimum of Subset of Real Numbers|supremum]]:
:$n-1 \notin S$
Thus $n-1 < x$.
Thus $n$ is an [[Definition:Integer|integer]] such that:
:$n-1 < x \leq n$
So $n$ is the [[Definition:Ceiling Function/Definition 3|ceiling function by definition 3]].
{{qed|lemma}}
=== Definition 3 equals Definition 2 ===
Let $n$ be an [[Definition:Integer|integer]] such that:
:$n-1 < x \leq n$
We show that $n$ is the [[Definition:Smallest Element|smallest element]] of the [[Definition:Set|set]]:
:$S = \left\{ {m \in \Z: m \ge x}\right\}$
Let $m \in \Z$ such that $m \ge x$.
We show that $n\leq m$.
{{AimForCont}} $m < n$.
By [[Weak Inequality of Integers iff Strict Inequality with Integer minus One]]:
:$m \le n - 1$
and so from the definition of $g$ it follows that $m < x$.
By [[Proof by Contradiction]] it follows that $m \ge n$.
Because $m \in S$ was arbitrary, $n$ is the [[Definition:Smallest Element|smallest element]] of $S$.
Thus $n$ is the [[Definition:Ceiling Function/Definition 2|ceiling function by definition 2]].
{{qed}}
\end{proof}

config: wiki


TITLE: Proof of overall momentum conservation
QUESTION [0 upvotes]: In $\phi^4$ momentum conservation at each vertex comes from integration in the vertex.
How can from momentum conservation at each vertex proof overall momentum conservation that is if $k_1,k_2 ,k_3,k_4$ are external lines how can we proof that  the deltas function  at each vertex turn into one delta function, namely $\delta(k_1+k_2 -k_3-k_4)$

REPLY [1 votes]: The number of 4-momentum conserving delta functions is one more than the number of (independent) internal momenta. After performing integrations over the internal momenta, you will be left with one remnant delta function. This is the total 4-momentum conserving delta function every amplitude must have. Can you see this for $\phi^4$ theory?


text:
\begin{document}
\maketitle
\begin{abstract}
 The value of a highly oscillatory integral is typically determined asymptotically by the behaviour of the integrand near a small number of critical points. These include the endpoints of the integration domain and the so-called stationary points or saddle points -- roots of the derivative of the phase of the integrand -- where the integrand is locally non-oscillatory. Modern methods for highly oscillatory quadrature exhibit numerical issues when two such saddle points coalesce. On the other hand, integrals with coalescing saddle points are a classical topic in asymptotic analysis, where they give rise to uniform asymptotic expansions in terms of the Airy function. In this paper we construct Gaussian quadrature rules that remain uniformly accurate when two saddle points coalesce. These rules are based on orthogonal polynomials in the complex
plane. We analyze these polynomials, prove their
existence for even degrees, and describe an accurate and efficient numerical scheme for the evaluation of oscillatory integrals with coalescing saddle points.
\end{abstract}

\section{Introduction and statement of results}
\label{sect_intro}

\subsection{Introduction}

Highly oscillatory integrals are a challenge for numerical integration methods, as the oscillatory nature of an integrand typically necessitates a large number of quadrature points. However, efficient numerical methods have been described for oscillatory integrals of the form
\begin{equation}\label{eq:I}
 I[f]=\int_a^b f(x)e^{i\omega g(x)} {\rm d}x,
\end{equation}
where both $f$ and $g$ are smooth functions and the oscillations can be attributed to a large value of the frequency parameter $\omega$. Examples include Filon-type quadrature, Levin quadrature, the numerical method of steepest descent and others \cite{iserles2005quad,olver2006levin,huybrechs2006osc1,huybrechs2009hoq,deano2018book}. Integrals of this form, and variations thereof, arise in a variety of applications typically involving wave phenomena.

An advantageous property of integral \eqref{eq:I} is that it is amenable to asymptotic analysis, and its value can be approximated using the method of stationary phase or the method of steepest descent \cite{bleistein1975asymptotic,wong2001asymptotic}. This approximation, being asymptotic, improves with increasing $\omega$. This is unlike classical quadrature schemes, which deteriorate with increasing $\omega$. The main goal of highly oscillatory quadrature methods, such as the methods mentioned above, is to combine improved accuracy for increasing $\omega$ with numerical convergence for any value of $\omega$, at a computational cost that is independent of $\omega$.

Asymptotic analysis of oscillatory integrals becomes more involved in the presence of \emph{stationary points} or \emph{saddle points}. These are roots of the derivative of the phase function, $g'(\xi)=0$, around which the integrand is locally non-oscillatory. Highly oscillatory quadrature techniques have to explicitly take such stationary points into account. The situation worsens when two stationary points are close to each other or when two stationary points coalesce for a particular value of a parameter. A canonical example is given by the cubic oscillator with a linear perturbation,
\begin{equation}\label{eq:osc}
 g(x,c)=\frac{x^3}{3}-c x.
\end{equation}
This oscillator features stationary points at $\pm \sqrt{c}$, coalescing at $c=0$.

Classical Poincar{\'e}-type asymptotic expansions (i.e., using just integer powers of $\omega^{-1})$ break down in the presence of two coalescing saddle points. This problem is well-known in asymptotic analysis and the solution is to consider uniform asymptotic expansions, see for example \cite{olver1974specialfunctions,bleistein1975asymptotic,wong2001asymptotic}. They are uniform in the sense that they are valid for a range of the parameter $c$, including its critical value $c=0$. Uniform asymptotic expansions typically involve a special function that captures the special transitional behaviour around a critical value of a parameter. In the case of two coalescing saddle points, the special function is the classical Airy function, which itself has an integral representation that involves a cubic oscillator. An example of such an expansion in this paper is \eqref{eq:airy_expansion} further on.

Numerical methods based on the existence of asymptotic expansions do not necessarily break down completely, but they certainly deteriorate in the presence of coalescing saddle points. In this paper we explore the analogue of uniform asymptotic expansions for one scheme, the numerical method of steepest descent. Though a great variety of uniform asymptotic expansions have been described in the literature for variations of integral \eqref{eq:I}, the Airy case is a canonical example, which has received the most study. For that reason, we pursue this case in detail in this paper. Experiments indicate that other cases that require uniform asymptotics are amenable to similar numerical techniques.

The goal of this paper is the construction and analysis of a uniformly applicable quadrature rule, uniform in the parameter $c$, for the canonical integral
\begin{equation}\label{eq:quadrature_rule}
 \int_{-1}^1 f(x) e^{i \omega \left( \frac{x^3}{3}-cx \right)} {\rm d}x \approx \sum_{k=1}^n w_k f(x_k),
\end{equation}
where $n$ is small and independent of $\omega$. We focus on a quadrature rule with optimal asymptotic order, in the sense that the error decays at the fastest algebraic rate in $\omega^{-1}$ among all quadrature rules with $n$ points. A consequence is that the points and weights depend on $\omega$ and $c$ and, furthermore, that the points $x_k$ typically lie in the complex plane. For that reason, we assume that $f$ is an analytic function at least in an open neighbourhood of $[-1,1]$ in the complex plane. Optimal quadrature rules involving real quadrature points in the interval $[-1,1]$ only are the subject of ongoing research, but in any case their convergence for large $\omega$ is slower than that of rules with complex points.

\subsection{Main results and outline of the paper}

The analysis in this paper centers around a family of orthogonal polynomials $p_{n,\delta}(z)$. They are monic and orthogonal with respect to a complex-valued oscillatory weight function,
\begin{equation}\label{eq:orthogonality}
 \int_\Gamma p_{n,\delta}(z) z^k e^{i \left( \frac{z^3}{3} - \delta z \right)} {\rm d}z=0,\qquad k=0,\ldots, n-1.
\end{equation}
Here, $\Gamma$ is any contour in the complex plane that connects the points $\infty \times e^{\frac{5i\pi}{6}}$ and $\infty \times e^{\frac{i\pi}{6}}$ at infinity. After a suitable rescaling, that will be detailed further on, the roots of these polynomials in combination with the roots of some other (known) polynomials give rise to the sought Gaussian quadrature rule of the form \eqref{eq:quadrature_rule}.

This is a non-classical setting of orthogonal polynomials, because the weight function is oscillatory. Hence, unique existence of the polynomials is not guaranteed for each value of $\delta$. We proceed by analyzing the corresponding Hankel determinants
\[
h_n = \det H_n,
\]
where $H_n$ is the Hankel matrix given by
\begin{equation}\label{eq:hankel}
  H_n=\left[
  \begin{array}{cccc}
         \mu_0 & \mu_1 & \cdots & \mu_{n-1}\\
         \mu_1 & \mu_2 & \cdots & \mu_{n}\\
         \vdots & \vdots & & \vdots\\
         \mu_{n-1} & \mu_{n} & \cdots & \mu_{2n-2}
  \end{array}
  \right]
\end{equation}
in terms of the \emph{moments} $\mu_k$ of the weight function,
\begin{equation} \label{eq:moments}
 \mu_k = \int_\Gamma z^k e^{i \left( \frac{z^3}{3} - \delta z \right)} {\rm d}z.
\end{equation}
These quantities are all functions of $\delta$, which we omit in our notation. Existence of the polynomial $p_{n,\delta}(z)$ for a particular value of $\delta$ is equivalent to the non-vanishing of $h_n$, since the latter appears in the denominator of the well-known determinant formula for orthogonal polynomials that remains valid in our setting,
\begin{equation}\label{eq:pn_determinant}
  p_{n,\delta}(x)=\frac{1}{h_n}\det\!
  \left[
  \begin{array}{ccccc}
         \mu_0 & \mu_1 & \cdots & \mu_{n-1} & 1\\
         \mu_1 & \mu_2 & \cdots & \mu_n & x\\
         \vdots & \vdots & & \vdots & \vdots\\
         \mu_n & \mu_{n+1} & \cdots & \mu_{2n-1} & x^n
  \end{array}
  \right].
\end{equation}

In the theoretical part of this paper, we show the following results.

\begin{theorem}\label{thm:existence_even}
 For any $n \in \mathbb{N}$ and $\delta \in \mathbb{R}$, $h_{2n} \neq 0$.
\end{theorem}
This implies that all monic even-degree polynomials $p_{2n,\delta}(x)$ are free of singularities as a function of $\delta$. Furthermore:

\begin{theorem}\label{thm:smalldelta}
For any $n \in \mathbb{N}$ and $\delta \in (-\infty,\delta_0)$, where $\delta_0\approx 2.338$ is
the smallest root of $\Ai(-\delta)$, $h_n \neq 0$.
\end{theorem}
Here, $\Ai$ is the classical Airy function \cite{NIST:DLMF,Olver:2010:NHMF}. The result shows that all polynomials of any degree exist for all values of $\delta$ on an interval that includes the negative halfline. Of particular interest in practice is that the interior of the interval includes the origin $\delta=0$, a critical value that corresponds to the case where two stationary points coalesce. Hence, as we will see, when two stationary points are sufficiently close, we can find quadrature rules with any desired number of quadrature points.

We denote the roots of $p_{n,\delta}(z)$, if the polynomial exists, by $t_{k,\delta}$, for $k=1,\ldots,n$. Thus, we have $p_{n,\delta}(z) = \prod_{k=1}^n (z-t_{k,\delta})$. As in the real-valued case, the Gaussian quadrature rule corresponds to the exact integral of the interpolating polynomial in the points $t_{k,\delta}$. Using Lagrange interpolation, this leads to a standard expression for Gaussian quadrature weights
\begin{equation}\label{eq:weights_lagrange}
 w_{k,\delta} = \int_\Gamma \frac{p_{n,\delta}(z)}{(z-t_{k,\delta})p_{n,\delta}'(t_{k,\delta})}e^{i \left( \frac{z^3}{3} - \delta z \right)}  {\rm d}z.
\end{equation}
An alternative expression is given further on in \eqref{eq:weights}. The weights are finite if all $t_{k,\delta}$ are distinct, such that the $p_{n,\delta}'(t_{k,\delta})$ does not vanish. This is always the case for real-valued polynomials with strictly positive weight function -- in which case all points are real and distinct and the weights are all positive --  but it is not necessarily so in the complex and oscillatory case.

In order to connect the oscillator of the original integral in \eqref{eq:I} with that of the oscillatory weight in \eqref{eq:orthogonality}, we introduce the scaling $\delta = c \omega^{2/3}$. Finally, we consider the quadrature formula
\begin{equation}\label{eq:rule}
 Q[f]= \frac{1}{\omega^{\frac{1}{3}}} \sum_{k=1}^n w_{k,\delta} f\left( \frac{t_{k,\delta}}{\omega^{1/3}} \right).
\end{equation}
This quadrature rule approximates the contribution of the two stationary points to the original integral \eqref{eq:I}, at least for large $\omega$. This contribution can be singled out by considering the integrand along a path in the complex plane, that is restricted to a neighbourhood of the origin. This leads to the following asymptotic error estimate.

\begin{theorem}\label{thm:error}
Let $f$ be analytic in a disk $D_r$ around the origin of radius $r > 0$. Furthermore, let $\Gamma$ be the concatenation of a straight line from $\infty \times e^{\frac{5i\pi}{6}}$ to $0$, and a straight line from $0$ to $\infty \times e^{\frac{i\pi}{6}}$. Finally, let $\tilde{\Gamma} = \Gamma \cap D_r$ and consider the integral
\begin{equation}\label{eq:Igamma}
 I_{\tilde{\Gamma}}[f] = \int_{\tilde{\Gamma}} f(z)e^{i\omega \left(\frac{z^3}{3}-cz\right)}dz.
\end{equation}

Assume that the polynomial $p_{n,\delta}$ of degree $n$ satisfies \eqref{eq:orthogonality} with $\delta=c \omega^{2/3}$ and that it has $n$ distinct finite roots $t_{n,k}$. For fixed $\delta$ and increasing $\omega$, which implies $c = {\mathcal O}(\omega^{-2/3})$, the error for the quadrature rule $Q[f]$ given by \eqref{eq:rule} is
 \[
  I_{\tilde{\Gamma}}[f] - Q[f] = {\mathcal O}\left(\omega^{-\frac{2n+1}{3}}\right), \qquad \omega \to \infty.
 \]
\end{theorem}

We describe the numerical method in \S\ref{sect_numerical_method}. Numerical examples are included in \S\ref{sect_experiments} and the  experiments in that section highlight several interesting features of the quadrature approach. These features are analyzed and explained theoretically in \S\ref{sect_analysis}, which includes also proofs of Theorems \ref{thm:existence_even} and \ref{thm:smalldelta} above. A brief asymptotic error analysis is carried out in \S\ref{sect_error_analyis}, leading to the proof of Theorem \ref{thm:error}. Finally, the construction of the quadrature rule is detailed in \S\ref{sect_construction}. Here, the major complication is the fact that the rule depends on two parameters, $\omega$ and $c$, and has to be computable efficiently on the fly in applications. We end the paper in \S\ref{sect_application} with an example application to a problem that attracts current interest in the literature on numerical methods for oscillatory integrals \cite{dominguez2011filonclenshawcurtis,ledoux2012interpolatory,huybrechs2012superinterpolation}.

\section{The numerical method}\label{sect_numerical_method}

The method of steepest descent is one classical way to derive asymptotic expansions for oscillatory integrals, in which the integration path is explicitly deformed into the complex plane. It goes back to Riemann and Cauchy \cite{bleistein1975asymptotic,wong2001asymptotic}. The purpose of the so-called \emph{numerical} method of steepest descent is to evaluate the resulting line integrals numerically, rather than asymptotically, using Gaussian quadrature rules. Our description is based on \cite{huybrechs2006osc1}, but earlier methods similar in spirit have been described in literature before for specific applications (e.g. \cite{todd54eei,gautschi70ecc,chandlerwilde95ecg}). In this section we develop a generalization of this approach to integrals with coalescing saddle points. Path deformation is only implicit in this case, since the contour is determined by connecting the roots of orthogonal polynomials in the complex plane.

\subsection{The numerical method of steepest descent (NSD)}\label{ss:NSD}

In the method of steepest descent for integral \eqref{eq:I}, the path of integration is deformed onto the steepest descent paths for the oscillator $g(x)$. For integral \eqref{eq:I} this results in:
\begin{itemize}
\item a half infinite path integral $\Gamma_{\{a,b\}}$ through each endpoint of the interval $[a,b]$,
\item and a double infinite path $\Gamma_\xi$ through each stationary point $\xi$ of the osillator $g(x,c)$.
\end{itemize}
For the details of the numerical scheme, we refer to the references \cite{huybrechs2006osc1,deano2009complexgauss}. The paths are such that the weight function $e^{i \omega g(x)}$ is non-oscillatory along these paths. This is achieved by following a level curve of the real part of $g$. For the case of $g(x,c)$, given by \eqref{eq:osc}, there are two stationary points: $\pm \sqrt{c}$. We denote the associated paths by $\Gamma_+$ and $\Gamma_-$. A typical illustration of the deformed paths is shown in Figure~\ref{fig:paths}.

\begin{figure}[t]
\begin{center}
\begin{overpic}[scale=0.4,unit=1mm]{Deform_SD_arrow.pdf}
\put(20,38){$\Gamma_a$}
\put(77,38){$\Gamma_b$}
\put(30,45){$\Gamma_-$}
\put(63,45){$\Gamma_+$}
\end{overpic}
\caption{Steepest descent paths for an oscillatory integral on $[-1,1]$ with two stationary points: two half infinite paths $\Gamma_a$ and $\Gamma_b$ originate from the endpoints $\pm 1$, and two doubly-infinite paths $\Gamma_{+/-}$ pass through the stationary points.}\label{fig:paths}
\end{center}
\end{figure}

The steepest descent integrals can be parameterized in a way that makes them suitable for Gaussian quadrature. In particular, the half-infinite paths can be written in the form
\[
 \int_0^\infty u(t) e^{-\omega t} {\rm d}t.
\]
Up to a scaling by $\omega$, this integral can be evaluated with Gauss-Laguerre quadrature~\cite{huybrechs2006osc1}. Interestingly, it can be shown that an $n$-point Gauss-Laguerre rule carries an error of the order ${\mathcal O}\left(\omega^{-2n-1}\right)$ for this integral. Truncating the asymptotic expansion of the same integral after $n$ terms leads to an error of size ${\mathcal O}\left(\omega^{-n-1}\right)$. The difference by a factor of nearly two in the exponents is due to the Gaussian nature of the quadrature and is the reason for the (asymptotic) optimality of this approach. For large values of $\omega$, i.e. for very highly oscillatory integrals, the approximation error is likely to be very small even when $n$ is a small number.

Similarly, the doubly-infinite paths can be written in the form (see \cite{deano2009complexgauss})
\begin{equation}\label{eq:hermite}
 \int_{-\infty}^\infty u(t) e^{-\omega t^2} {\rm d}t,
\end{equation}
well suited for Gauss-Hermite quadrature. The error in this case behaves as ${\mathcal O}\left(\omega^{-\frac{2n+1}{2}}\right)$. Here, too, the exponent is twice as large as in the error term of an $n$-term truncated asymptotic expansion. In both cases, the quadrature points correspond to function evaluations of $f$ at points that lie exactly on one of the steepest descent paths shown in Figure \ref{fig:paths}.

However, the integrand $u(t)$ in the integral corresponding to the stationary point at $-\sqrt{c}$ has a singularity in the complex plane, that arises from the other stationary point at $\sqrt{c}$ (and vice-versa). As the parameter $c$ decreases, this singularity of $u(t)$ moves closer towards the real axis in \eqref{eq:hermite}. As such, though the method in principle applies for any $c > 0$, accuracy deteriorates for small $c$ as the convergence rate of Gauss-Hermite quadrature decreases. We illustrate this with a numerical experiment further on. At $c=0$, the saddle points coincide and $g'(x)$ has a double root. In this case, Gauss-Hermite quadrature no longer applies. Like the asymptotic expansions themselves, the numerical method for $c>0$ does not apply to the case $c=0$.

\subsection{Uniform asymptotic expansion of oscillatory integrals}\label{ss:uniform}

Uniform asymptotic expansions for oscillatory integrals with coalescing saddle points are usually formulated in terms of the Airy function and its derivative \cite{bleistein1975asymptotic,olver1974specialfunctions}. The uniform expansion is not unique and variations are possible, but one statement is:
\begin{equation}\label{eq:airy_expansion}
\int_\Gamma f(x)e^{i\omega\left(\frac{x^3}{3}-cx\right)}dx\sim \frac{1}{\omega^{\frac{1}{3}}}\sum_j\frac{a_{2j}}{\omega^{j}}\Ai(-c\omega^{\frac{2}{3}})+\frac{1}{\omega^{\frac{2}{3}}}\sum_j\frac{a_{2j+1}}{\omega^{j}}\Ai'(-c\omega^{\frac{2}{3}}).
\end{equation}
Here, the coefficients $a_j$ are determined by $f,c$ and $\omega$ in a complicated way as follows. First, the function $f$ is written in the form $f(x) = (x^2-c) f_1(x) + f_2(x)$. Next, integration by parts is performed for the first term in this sum (note that $(x^2-c) = g'(x,c)$). This results in an Airy function and a new integral along $\Gamma$, after which the process is repeated recursively.

There are several numerical issues with expansion \eqref{eq:airy_expansion}. First, the asymptotic expansion does not necessarily converge. This is true in general: asymptotic expansions rarely converge \cite{boyd1999devil}. One exception is for polynomials $f$, in which case the expansion terminates and the above expression becomes exact. The second issue is less often reported, but equally profound: the numerical computation of the coefficients $a_j$ is exceedingly difficult and the evaluation of the truncated expansion is numerically unstable \cite{temme1997airy}. This is due essentially to the fact that the coefficients $a_j$ are obtained after repeated application of L'H{\^o}spital's rule, leading in practice for small $c$ to extensive cancellation errors.

In spite of its numerical issues, expansion \eqref{eq:airy_expansion} is valid uniformly for $c \in [-C_1,C_2]$ for any $C_1,C_2 > 0$. For large $c$, the Airy functions can be expanded asymptotically and the Poincar{\'e}-type asymptotic expansion of the integral is recovered. For small $c$, the Airy function precisely captures the behaviour of the coalescing saddle points.

\subsection{A uniformly valid quadrature rule (UNSD)}\label{sect_quad_uniform}

The numerical method of steepest descent starts with an explicit path deformation, followed by the application of a Gaussian quadrature rule. The points of the quadrature rule lie exactly on the steepest descent paths. Here, we reverse the order of the steps. We formulate a Gaussian quadrature rule, one that evaluates the path integrals $\Gamma_+$ and $\Gamma_-$ simultaneously. The roots of the orthogonal polynomials implicitly correspond to a contour in the complex plane. The resulting path deformation is illustrated in Figure~\ref{fig:combined_paths}.

\begin{figure}[t]
\begin{center}
\begin{overpic}[scale=0.5,unit=1mm]{deformed_contour1cut.pdf}
\put(16,16){$\Gamma_a$}
\put(55,16){$\Gamma$}
\put(94,25){$\Gamma_b$}
\end{overpic}
\caption{The figure shows the paths through the endpoints $\Gamma_a$ and $\Gamma_b$ as in Figure \ref{fig:paths}, but the paths $\Gamma_+$ and $\Gamma_-$ are combined into a single contour $\Gamma$.}\label{fig:combined_paths}

\end{center}
\end{figure}

Let us be more precise. For the time being, we assume that $f$ is sufficiently analytic and does not grower faster than exponentially at infinity, such that all path deformations are justified by Cauchy's integral theorem. We assume that the endpoint integrals are treated using Gauss-Laguerre quadrature as before and focus on the doubly-infinite paths from now on. Define a contour $\Gamma$ that connects the endpoints of the half-infinite paths, for example $\Gamma = \Gamma_+ \cup \Gamma_-$ where $\Gamma_{+}$ and $\Gamma_{-}$ are the steepest descent paths at both stationary points. Equivalently, we can consider any contour $\Gamma$ that connects the points $\infty \times e^{\frac{5i\pi}{6}}$ and $\infty \times e^{\frac{i\pi}{6}}$ at infinity, such as the one defined in Theorem~\ref{thm:error}.

Upon the change of variables
\begin{equation}\label{eq:change_of_variables}
 t=\omega^{\frac{1}{3}}x \qquad \mbox{and}\qquad \delta=c\omega^{\frac{2}{3}},
\end{equation}
we have (recall the definition \eqref{eq:osc} of $g$)
\[
 \int_\Gamma f(x) e^{i \omega g(x,c)} {\rm d}x = \omega^{-1/3} \int_\Gamma f\left( \frac{t}{\omega^{1/3}}\right) e^{i g(t,\delta)} {\rm d}t.
\]
Note that the path of integration does not change under this scaling, assuming $f$ is analytic in a sufficiently large region. The advantage of the latter reformulation is that the oscillator depends on only one parameter, $\delta$, rather than two.

Consider a family of monic orthogonal polynomials $p_{n,\delta}(x)$ satisfying the orthogonality conditions \eqref{eq:orthogonality}. These polynomials should for the time being be considered only formally orthogonal, as their existence is not guaranteed for this oscillatory weight function in the complex plane. Assuming the monic orthogonal polynomial of degree $n$ exists uniquely for some particular value of $\delta$, i.e. $h_{n-1} \neq 0$, there is a corresponding Gaussian quadrature rule with $n$ points and weights. This rule is suitable for weighted integration along $\Gamma$,
\[
 \int_\Gamma u(t) e^{i g(t,\delta)} {\rm d}t \approx \sum_{k=1}^n w_{k,\delta} \, u(t_{k,\delta}).
\]
It is emphasized in this notation that the points and weights depend on the parameter $\delta$. Note that this rule plays the same role as the Laguerre and Hermite rules before.

For the original integral we arrive, after undoing the transformation \eqref{eq:change_of_variables}, at the quadrature rule \eqref{eq:rule}.

This is the quadrature rule we propose and investigate in this paper. It is denoted in the following by UNSD, for uniform numerical method of steepest descent. Note that any phase function with two stationary points can be mapped to the canonical case $g(x,c)$ by a smooth change of variables.

This quadrature rule can be seen as the numerical equivalent of the uniform asymptotic expansion. It is clear that the rule applies to any value of $c$, including $c=0$. That case corresponds to $\delta=0$ in our notation, and the stationary points coalesce. A special case treatment for such a degenerate stationary point was described before \cite{deano2009complexgauss}, and the current quadrature method simply reduces to that case exactly. As it turns out, for large $c$ the rule above also reduces (numerically) to two individual applications of Gauss-Hermite along the steepest descent paths $\Gamma_+$ and $\Gamma_-$ separately, shown in Fig. \ref{fig:paths}. For large $c$, and thus for large $\delta$, the orthogonal polynomial $p_{2n,\delta}(x)$ of degree $2n$ is close to the product of two Hermite polynomials of degree $n$. Like the uniform asymptotic expansion reduces to the regular Poincar{\'e}-type expansions for large $c$, the uniform numerical scheme reduces to the regular NSD scheme in the same regime.

A disadvantage of the quadrature rule that should be noted is its dependence on the parameter $\delta$. Since this parameter may have any value in applications, the quadrature rule has to be computable on the fly. Thus, unless a priori computations are feasible, the construction of the quadrature rule is an integral part of the cost of the application of the scheme. For phase functions different from \eqref{eq:osc}, the cost of the change of variables has to be taken into account as well.

\section{Numerical experiments}\label{sect_experiments}

Before analyzing the polynomials and the associated quadrature rule, we perform a number of experiments to illustrate the accuracy that can be achieved for varying values of $c$ and $\omega$. We show results for the simple analytic function
\[
 f(x) = \sin(4x).
\]
This function is simple, yet not entirely innocent as it grows exponentially in the complex plane. This growth, as we shall see, offsets the benefit of steepest descent deformation for small $\omega$. A more interesting example is given in \S\ref{sect_application} later on.

\setlength{\figurewidth}{9cm}
\setlength{\figureheight}{7cm}

\begin{figure}[t]
\centering

\begin{small}
\input{unsd_abs_error.tikz}
\end{small}
\caption{Approximation error (absolute error) of the UNSD method as a function of $\omega$ for $f(x)=\sin(4x)$ with a $6$-point rule ($n=6$) and for various values of $c$.}\label{fig:error_omega}
\end{figure}

\subsection{Asymptotic order of convergence}

Our first observation is that the error of the UNSD scheme decreases algebraically with increasing $\omega$. This is illustrated in Figure~\ref{fig:error_omega}. The error decreases as a function of $\omega$ at a similar rate for all $c$, though it appears to be smaller for smaller $c$. Though the number of quadrature points is very modest, only $n=6$, absolute errors on the order of $1e-10$ are reached for $\omega$ as small as $100$.

In contrast, for small values of $\omega$ the error is fairly large. As mentioned above, this is exacerbated by the exponential growth of $f$ in the upper half of the complex plane. However, the lack of convergence for small $\omega$ is typical. This cause is most obvious from the scaling of the roots in \eqref{eq:rule}: the term  $\omega^{-1/3} t_{k,\delta}$ implies that the quadrature points $t_{k,\delta}$ are mapped closer to the real line for increasing $\omega$, but away from the real line for decreasing $\omega$. The rule as formulated here does not have a proper limit $\omega \to 0$. Yet, recall that for small $\omega$ the original integral is non-oscillatory, and a straightforward approach is to evaluate it by other means in this regime.

\subsection{Accuracy for small $c$ and comparison to NSD}

In Figure~\ref{fig:unsd_vs_hermite} we compare the error of the proposed UNSD-method with the usual NSD scheme based on Hermite polynomials described in \S\ref{ss:NSD}, and analyzed in \cite{deano2009complexgauss}. For NSD we evaluate the two line integrals through $\pm \sqrt{c}$ using Gauss-Hermite quadrature with $n$ points each. For UNSD we evaluate their sum using a single quadrature rule with $2n$ points in total. It is clear that the NSD method fails for small $c$: the error blows up as $c$ tends to $0$. In contrast, the error of the UNSD scheme is uniformly small in $c$, and that is the motivation of this work.

\setlength{\figurewidth}{9cm}
\setlength{\figureheight}{7cm}

\begin{figure}[t]
\centering

\begin{small}
\input{nsd_vs_unsd.tikz}
\end{small}
\caption{A comparison of the accuracy of two applications of the Gauss-Hermite rule versus the UNSD rule, for $n=6$ and $f(x)=\sin (4x)$, as a function of $c$.}\label{fig:unsd_vs_hermite}
\end{figure}

It is also interesting to observe in this figure that the UNSD and NSD rules appear to exhibit very similar errors for large $c$. It turns out that in this regime the quadrature points are nearly the same. We illustrate the location of the quadrature points in Figure~\ref{fig:roots_split}.

The left panel shows the location in the complex plane of the quadrature points for various $\delta$. When $\delta$ is small, the roots lie on a single curve. For larger values of $\delta$, the roots seemingly cluster in two separate curves. These two clusters are illustrated again in the right panel and compared to the quadrature points of the NSD scheme: they are indeed very close to each other.

\setlength{\figurewidth}{7cm}
\setlength{\figureheight}{6cm}

\begin{figure}[!t]
\centering
\subfigure{

\input{delta_points.tikz}
}
\subfigure{

\input{unsd_hermite_points.tikz}
}

\caption{Illustration of the roots of $p_{n,\delta}$. The left panel shows that the roots seem to lie on a single curve for small $\delta$, but split in two groups for larger $\delta$ ($n=14$, $\delta$ increases from $0$ to $12$ from top to bottom). The right panel shows that the roots of the UNSD-rule and the NSD rule for a large value of $\delta$ are close to each other ($\delta=10$).}\label{fig:roots_split}
\end{figure}

\subsection{Issues with quadrature rules with an odd number of points}

The previous experiments have shown that the proposed UNSD method works well for all values of $c$. However, thus far we have used only quadrature rules with an even number of points. One quickly observes that there are serious issues with the quadrature rules with an odd number of points. This is illustrated in Figure \ref{fig:odd_vs_even}: the quadrature rules with odd values of $n$ lead to large errors, at least for some isolated values of $c$. It turns out that the polynomials introduced in the previous section do not exist for all values of $c$, when $n$ is odd. Experiments seem to indicate that all even-degree polynomials do exist, regardless of the value of $c$. There seem to be no problems of existence for small values of $c$, which is the intended regime for the quadrature rule. Yet, it seems advisable to use only an even number of points in computations in general.

\setlength{\figurewidth}{7cm}
\setlength{\figureheight}{6cm}

\begin{figure}[t]
\begin{center}
\subfigure[Even $n$]{

\begin{small}
\input{error_even_n.tikz}
\end{small}
}
\subfigure[Odd $n$]{

\begin{small}
\input{error_odd_n.tikz}
\end{small}
}
\caption{Illustration of the accuracy of the UNSD rule applied to $f(x)=\sin x + \cos x$, for various numbers of points $n$ and as a function of $c$ (using $\omega=1$). For odd values of $n$, the quadrature error exhibits large spikes at particular values of $c$. These spikes correspond to non-existence of the underlying orthogonal polynomials.}\label{fig:odd_vs_even}
\end{center}

\end{figure}

\section{Proofs of Theorems \ref{thm:existence_even} and \ref{thm:smalldelta} }
\label{sect_analysis}

In this section we establish analytical results that completely describe the features we have observed with the experiments. The proofs rely on the relation of the moments \eqref{eq:moments} with the Airy function and its derivative.

In order to simplify the expressions, as well as to make connections to existing literature on orthogonal polynomials, in this section we adopt the integration contour $C$ of the standard Airy function (defined below), rather than contour $\Gamma$ in \eqref{eq:orthogonality}. We note that $C$ is simply a rotation of $\Gamma$, and the polynomials $P_n(s)$ in this section are related to $p_{n,\delta}$ by
\[
 P_n(s) = i^n p_n^{-x}(-i s).
\]

\subsection{Airy function and determinants}
The Airy function $\Ai$ is the solution of the Airy differential equation $y'' = xy$
with the asymptotic behavior
\begin{equation} \label{eq:Airyasymp}
	\Ai(x) = \frac{1}{2\sqrt{\pi}x^{1/4}} e^{-\frac{2}{3} x^{3/2}} \left(1 + O(x^{-3/2})\right)
\end{equation}
as $x \to +\infty$. It has an integral representation
\begin{equation} \label{eq:Airyintegral}
	\Ai(x) = \frac{1}{2\pi i} \int_C e^{- \frac{1}{3} s^3 + xs} ds
	\end{equation}
where $C$ is an infinite contour in the complex $s$-plane from $ \infty e^{-2\pi i/3}$
to $\infty e^{2\pi i/3}$. Then
\[ \Ai^{(k)}(x) = \frac{1}{2\pi i} \int_C s^k e^{- \frac{1}{3} s^3 + xs} ds \]
and a comparison with the definition \eqref{eq:moments} of $\mu_k$ shows that
\begin{equation} \label{eq:Airymuk}
	\mu_k =  2 \pi (-i)^k \Ai^{(k)}(-\delta).
	\end{equation}

We define $D_0(x) = 1$ and for $n \geq 1$,
\begin{equation} \label{eq:Airydet}
	D_n(x) = \det \left[ \frac{d^{j+k-2}}{dx^{j+k-2}} \Ai(x) \right]_{j,k=1, \ldots, n}.
	\end{equation}
Then $D_n$ is an entire function, which is related to the determinant $h_n = \det H_n$
as follows.
\begin{lemma}
We have
\begin{equation} \label{eq:Airydethn}
	h_n = \det H_n =  (2\pi)^n (-1)^{n(n-1)/2}  D_n(-\delta)
	\end{equation}
\end{lemma}
\begin{proof}
Because of \eqref{eq:hankel} and \eqref{eq:Airymuk} we have
\[ h_n = (2 \pi)^n \det \left[ (-i)^{j+k-2} \Ai^{(j+k-2)}(-\delta) \right]_{j,k=1, \ldots, n}. \]
We take out a factor $(-i)^{j-1}$ out of row $j$ and a factor $(-i)^{k-1}$ out of column $k$
of the determinant, for each $j,k = 1, \ldots, n$.
This results in the factor
\[ \prod_{j=1}^n (-i)^{j-1} \cdot \prod_{k=1}^n (-i)^{k-1} = (-1)^{n(n-1)/2}. \]
and the formula \eqref{eq:Airydethn} follows.
\end{proof}

In view of \eqref{eq:Airydethn} we need to prove that $D_n$ has no real zeros if $n$ is even,
and has no zeros in $(\iota_1, \infty)$, where $\iota_1$ is the largest zero of the Airy function
if $n$ is odd.
Recall that $\Ai$ has an infinite number of zeros, all  negative and simple, that are
usually denoted as $0 > \iota_1 > \iota_2 > \iota_3 > \cdots$.

This claim is easy to verify for small values of $n$. Clearly $D_1 = \Ai$
has no zeros in $(\iota_1, \infty)$. For $n=2$ we have by \eqref{eq:Airydet}
\[ D_2(x) = \Ai(x) \Ai''(x) - \Ai'(x)^2 = x \Ai(x)^2 - \Ai'(x)^2, \]
where we used the Airy differential equation $\Ai'' = x \Ai$. Then by
a simple calculation
\begin{equation} \label{eq:D2formula}
	D_2'(x) = \Ai(x)^2,
	\end{equation}
and, since $D_2(x) \to 0$ as $x \to +\infty$, we get
\[  D_2(x) = - \int_x^{\infty} \Ai(s)^2 ds < 0, \qquad \text{for } x \in \mathbb R. \]
Thus $D_2$ has no real zeros indeed.

The proof for general $n$ will follow from certain identities for the Airy
determinants.

\begin{remark}
\normalfont
The Airy determinant \eqref{eq:Airydet} appears in the work of Forrester and Witte \cite{FW} within the context of random matrix theory, see formula (1.25). It arises in the soft edge scaling limit of the expectation value of the $n$-th power of the characteristic polynomial of a GUE matrix. From this interpretation it is natural that $D_n$ does not vanish if $n$ is even. The quantity $D_2$ is (up to a sign) the density of the Airy point process.
\end{remark}

\begin{remark}
\normalfont
The Airy determinant \eqref{eq:Airydet} is also a building block for the construction of Airy like solutions of the Painlev{\'e} II equation, see for example Theorem 4 in \cite{Clarkson} and references therein.
\end{remark}

\subsection{Identities for Airy determinants}

\begin{lemma} \label{lem:Dnidentity}
We have the differential identity
\begin{equation} \label{eq:Airydetdiffid}
	D_{n-1} D_{n+1} = D_n D_n'' - (D_n')^2.
	\end{equation}
\end{lemma}
\begin{proof}
This can be proved as in [9, Lemma 6] using the connection with
orthogonal polynomials, and their recurrence coefficients (see also
the proof of Proposition \ref{prop:Airydet}).

However, it can also be proved from Jacobi's identity for determinants,
see \cite[formula (0.8.4.1)]{horn1985}. A special case of this formula is the following. For a square matrix
$A$ we use $A[i;j]$ to denote the matrix obtained from $A$ by deleting
row $i$ and column $j$, and we use $A[ij;kl]$ for the matrix obtained
from $A$ by deleting rows $i$ and $j$ and columns $k$ and $l$,
where $i<j$ and $k< l$. Then
\begin{equation} \label{eq:JacobiId}
	\det A  \cdot= \det A[ij;kl] = \det A[i;k] \cdot \det A[j;l]
	- \det A[i;l] \cdot \det A[j;k].
	\end{equation}
Now let $A$ be the $(n+1) \times (n+1)$ matrix
\[ A = \left[ \frac{d^{j+k-2}}{dx^{j+k-2}} \Ai(x) \right]_{j,k=1, \ldots, n+1}. \]
Then $\det A = D_{n+1}$, $\det A[n+1; n+1] = D_n$,
$\det A[n n+1; n n+1] = D_{n-1}$, $\det A[n; n+1] = \det A[n+1; n] = D_n'$
and $\det A[n;n] = D_n''$ and using this in \eqref{eq:JacobiId}
we obtain \eqref{eq:Airydetdiffid}.
\end{proof}

The differential identity \eqref{eq:Airydetdiffid} is well-known in integrable
systems where it is related to the Toda lattice equations, see e.g. \cite{hirota1987}. However the following identity seems to be new.
It is specific for the Airy determinants \eqref{eq:Airydet} and it
reduces to the identity \eqref{eq:D2formula} for $n=1$.

\begin{proposition} \label{prop:Airydet}
For each $n \geq 1$ we have
\begin{equation} \label{eq:Airydetid2}
	\left( \frac{D_{n+1}}{D_{n-1}} \right)' = n \left( \frac{D_{n}}{D_{n-1}}\right)^2.
	\end{equation}
\end{proposition}
\begin{proof}
Let $P_n$ be the monic orthogonal polynomial of degree $n$ for the complex weight
$e^{-\frac{1}{3}s^3 + xs}$ (with $x$ as a parameter) on the contour $C$.
That is
\begin{equation} \label{eq:Pnortho}
	\int_C P_n(s) s^k e^{-\frac{1}{3} s^3 + xs} ds = 0, \qquad k=0, \ldots, n-1. \end{equation}
The polynomial $P_n$ exists if and only if $D_n(x) \neq 0$.
Note that $P_n(s) = i^n p_{n,\delta}^{-x}(-is)$, where $p_{n,\delta}$ is the orthogonal polynomial from \eqref{eq:orthogonality}.

Three consecutive polynomials (if they exist) are related by a three term recurrence
\begin{equation}
\label{eq:recurrence_arno}
sP_n(s) = P_{n+1}(s) + b_n P_n(s) + a_n P_{n-1}(s)
\end{equation}
where
\begin{equation} \label{eq:anbn}
	\begin{aligned}
	a_n & = \frac{D_{n-1} D_{n+1}}{D_n^2},   \\
	b_n & = \frac{D_{n+1}'}{D_{n+1}} - \frac{D_n'}{D_n}.
	\end{aligned}
	\end{equation}
As a function of $x$ they satisfy the Toda differential equations
\begin{equation} \label{eq:anbnToda}
	\begin{aligned}
	\frac{da_n}{dx} & = a_n(b_n-b_{n-1}), \\
	\frac{db_n}{dx} & = a_{n+1}-a_n.
	\end{aligned}
	\end{equation}
The equations \eqref{eq:anbn} and \eqref{eq:anbnToda} are well-known
and  can be found in various forms in the literature, see e.g.
equations (2.9)  and (2.10) in \cite{sogo1993}, where  $a_n^2$ is used instead of $a_n$ and $\tau_n$
is used for  $D_{n+1}$. The equations \eqref{eq:anbnToda} hold
for recurrence coefficients whenever $x$-dependence of
the orthogonality weight comes in the form $w(s) e^{xs}$.

For the cubic weight $w(s) = -\frac{1}{3} s^3$ as in \eqref{eq:Pnortho}
the recurrence coefficients satisfy nonlinear difference equations
(sometimes called string equations)
\begin{equation} \label{eq:string}
	\begin{aligned}
	a_n + a_{n+1} & = x - b_n^2, \\
	a_n(b_n+b_{n-1}) & = n.
	\end{aligned}
	\end{equation}
These equations are also given in the recent paper \cite{clarkson2016}, where they were identified as an alternative discrete Painlev\'e I equation.
Note that the recurrence coefficients in our paper differ by a sign from the ones used by Clarkson et al.
See also the appendix of \cite{vanassche2015}. These are recent papers, but the idea that recurrence coefficients for semi-classical orthogonal polynomials satisfy nonlinear equations
of Painlev\'e type dates back to at least \cite{magnus1995}.

Inserting \eqref{eq:anbn} into the second string equation \eqref{eq:string}, we find
\begin{equation}
 n = \frac{D_{n-1} D_{n+1}}{D_n^2}  \left( \frac{D_{n+1}'}{D_{n+1}} - \frac{D_n'}{D_n} \right).
\end{equation}
which can be rearranged to give \eqref{eq:Airydetid2}.
\end{proof}

\subsection{Behavior for $x \to +\infty$}

As the last preparation for the proofs of the theorems, we have to investigate
the behavior of $D_n(x)$ as $x \to +\infty$.

The Airy function $\Ai$ and all its derivatives tend to zero as the argument tends
to infinity along the positive real line. Thus it is clear from \eqref{eq:Airydet}
that $D_n(x) \to 0$ as $x \to +\infty$. We have the more precise behavior.

\begin{lemma}
We have
\begin{equation} \label{eq:Airydetasymp}
	D_n(x) =  c_n  x^{-n^2/4}  	e^{-\frac{2}{3}n^{3/2}} \left(1+ O(x^{-3/2}) \right)
\end{equation}
as $x \to +\infty$, with the constant
\[
 c_n = (-1)^{n(n-1)/2}  \frac{\prod_{k=0}^{n-1} k!}{2^{n(n+1)/2} \pi^{n/2}}.
\]
\end{lemma}
\begin{proof}
The lemma is true for $n=0$ and $n=1$, see \eqref{eq:Airyasymp}.

For general $n$, we use $D_{n+1} = \frac{D_n^2}{D_{n-1}}  a_n$,  see \eqref{eq:anbn}, and the fact that
\[
a_n(x) =  - \frac{n}{2} x^{-1/2} \left(1 + O(x^{-3/2})\right)   \quad \text{ as } x \to +\infty
\]
see \cite[formula (5.4.a)]{clarkson2016}.
Then \eqref{eq:Airydetasymp} follows by an induction argument.
\end{proof}

\subsection{Proofs of Theorems \ref{thm:existence_even} and \ref{thm:smalldelta} }

\begin{proof}
From \eqref{eq:Airydetasymp} it follows that
$ \frac{D_{n+2}(x)}{D_{n}(x)}  \to 0$ as $ x \to +\infty$.
Thus, after changing $n \to n+1$ in \eqref{eq:Airydetid2} and integrating, we obtain
\begin{equation} \label{eq:Dn2integral}
	D_{n+2}(x) = - (n+1) D_n(x) \int_x^{+\infty} \left( \frac{D_{n+1}(s)}{D_n(s)} \right)^2 ds
	\end{equation}
provided that $D_n$ has no zeros on the interval $[x,\infty)$.

Now we can easily prove the theorems by an induction argument.
It is true that $D_0 =1$ has no zeros on the real line.
Suppose $n \geq 0$ is even, and assume that $D_n$ has no zeros on the real line.
Then the integral \eqref{eq:Dn2integral} is valid for every $x \in \mathbb R$.
Since $D_{n+1}$ clearly does not vanish identically (this follows from \eqref{eq:Airydetasymp}),
we immediately find that $D_{n+2}$ has no zeros on the real line as well.

For $n=1$, we have $D_1 = \Ai$, and the Airy function has no zeros on $(\iota_1, \infty)$.
Suppose $n \geq 1$ is odd, and assume that $D_n$ has no zeros on $(\iota_1, \infty)$.
Then again it follows from \eqref{eq:Dn2integral} that $D_{n+2}(x)$ has no zeros
 for $x > \iota_1$ as well.

This proves Theorems \ref{thm:existence_even} and \ref{thm:smalldelta}, where
we recall \eqref{eq:Airydethn}.
\end{proof}

The proof in this paper is based on the relation between the Hankel determinant and Airy determinants. An alternative proof, including more detailed analysis of the recurrence coefficients based on the string equations, can be found in \cite{lejon2016phd}.

\section{Asymptotic error analysis of the quadrature method}
\label{sect_error_analyis}

After showing the existence of the orthogonal polynomials, hence of the proposed quadrature rules, we set out to analyze its convergence characteristics. Since, like the NSD scheme itself, the current quadrature scheme is asymptotic, we focus on the behaviour for $\omega \gg 1$. We consider nearly coalescing saddle points, i.e. $0 < c \ll 1$, since that is the intended regime of the quadrature rule. We want to show that the error decays with $\omega$ at an algebraic rate, as observed in the experiments.

To that end, we investigate the $\omega$-dependence of the quadrature error $I_\Gamma[f] - Q[f]$. Recall the definitions \eqref{eq:Igamma} and \eqref{eq:rule}.

\begin{proof}[Proof of Theorem \ref{thm:error}]
In order to avoid stringent analyticity requirements on $f$, the contour $\Gamma$ is in Theorem \ref{thm:error} truncated to a neighbourhood of the origin in which $f$ is analytic. This does not alter the asymptotic behaviour of the integral, as can be confirmed using integration by parts and noting that the endpoints of the finite contour $\tilde{\Gamma} = \Gamma \cup D_r$ are in the sectors of the complex plane where the oscillator $e^{i \omega g(z,c)}$ becomes exponentially small. In the large $\omega$ limit with fixed $\delta$, the quadrature points $\frac{t_{k,\delta}}{\omega^{1/3}}$ in \eqref{eq:rule} are also in $D_r$. Since the quadrature rule is assumed to be well-defined for given $n$ and $\delta$, it can be applied for any $\omega$ and we investigate the large $\omega$ limiting behaviour.

For simplicity, in what follows we assume that $f$ is extended to a smooth function $\tilde{f}$ along $\Gamma$ in such a way that the extended integral
\[
 I_\Gamma[\tilde{f}] = \int_\Gamma \tilde{f}(z) e^{i \omega g(z,c)} dz
\]
exists and is finite. Regardless of the particular extension, we have that $I_{\tilde{\Gamma}}[f] - I_\Gamma[\tilde{f}] = {\mathcal O}(\omega^{-p})$ for all $p > 0$. It follows that $I_{\tilde{\Gamma}}[f] - Q[f] = I_\Gamma[\tilde{f}] - Q[f] + {\mathcal O}(\omega^{-p})$, for all $p > 0$.

Using the change of variables \eqref{eq:change_of_variables} and formula \eqref{eq:rule} the approximation error becomes, up to the above asymptotically small error,
\begin{align*}
 I_\Gamma[\tilde{f}]-Q[f] = \frac{1}{\omega^{\frac{1}{3}}} \int_\Gamma \tilde{f}\left(\frac{t}{\omega^{\frac{1}{3}}}\right)e^{i\left(\frac{t^3}{3}-\delta t\right)} {\rm d}t- \frac{1}{\omega^{\frac{1}{3}}} \sum_{k=1}^n  f\left(\frac{t_{k,\delta}}{\omega^{\frac{1}{3}}}\right)w_{k,\delta}.
\end{align*}
Recall that the quadrature points $t_{k,\delta}$ are the roots of the orthogonal polynomials $p_{n,\delta}$ and the coefficients $w_{k,\delta}$ are the associated weights.

Since the quadrature points are zeros of orthogonal polynomials, polynomials up to degree $2n-1$ are integrated exactly. Because $f$ is analytic in the disk $D_r$, it can be expanded in a Taylor series,
\[
f(x)=\sum_{j=0}^\infty a_j x^j, \qquad |x| < r.
\]
Since the first $2n-1$ terms are integrated exactly, we have in an asymptotic sense that
\begin{small}\begin{align*}
I_\Gamma[\tilde{f}]-Q[f]& \sim \frac{1}{\omega^{\frac{1}{3}}}\sum_{j=0}^\infty a_j\int_\Gamma \frac{t^j}{\omega^{\frac{j}{3}}} e^{i\omega g(t,\delta)} {\rm d}t - \frac{1}{\omega^{\frac{1}{3}}} \sum_{j=0}^\infty a_j \sum_{k=1}^n \frac{t_{k,\delta}^j}{\omega^{\frac{j}{3}}}w_{k,\delta} \\
&=\frac{1}{\omega^{\frac{1}{3}}} \sum_{j=2n}^\infty \frac{a_j}{\omega^{\frac{j}{3}}}\left[\int_\Gamma t^j e^{ig(t,\delta)} {\rm d}t - \sum_{k=1}^n t_{k,\delta}^j w_{k,\delta} \right]\\
&=\frac{1}{\omega^{\frac{1}{3}}} \sum_{j=2n}^\infty \frac{a_j}{\omega^{\frac{j}{3}}}\left[\mu_j(\delta)-\sum_{k=1}^n t_{k,\delta}^j w_{k,\delta} \right].
\end{align*}\end{small}
Note that interchanging summation and integration is not justified here, as the Taylor series of $f$ does not converge outside the circle with radius $R$, yet the integration contour is infinite. The result only holds as an asymptotic series, in exactly the same way as in the well-known Watson's Lemma \cite{bleistein1975asymptotic}.

For fixed $\delta$, the moments $\mu_j(\delta)$ and the sum $\sum_{k=1}^n t_{k,\delta}^j w_{k,\delta}$ are independent of $\omega$. Thus, again for \emph{fixed $\delta$} and for increasing $\omega$, we conclude from the size of the first term in the summation that
\[
I_\Gamma[f]-Q[f] = \mathcal{O}\left(\omega^{-\frac{2n+1}{3}}\right).
\]
This concludes the proof.
\end{proof}
Note that we have taken a rather simplified approach here. There are two parameters in the original problem, $c$ and $\omega$, and we have chosen to fix $\delta=c \omega^{2/3}$. The case of fixed $\delta$ and increasing $\omega$ implies that $c$ shrinks like $\omega^{-2/3}$. The rate of decay as a function of $\omega$ we arrive at in this regime is exactly the same as in the case where $\delta=0$ \cite{deano2009complexgauss}.

An analysis for the case of fixed $c$ and increasing $\omega$ would be considerably more involved, as this regime implies that $\delta$ grows like $\omega^{2/3}$. Hence, this regime includes the transition from a single cluster of quadrature points to a double cluster. Still, we know the outcome in the regime of larger $c$ as well, from the existing analysis of the NSD scheme. In NSD we can use $n$ Gauss-Hermite points for a single stationary point integral and achieve $\mathcal{O}\left(\omega^{-\frac{2n+1}{2}}\right)$ error. Using $n/2$ points for the first stationary point and $n/2$ points for the second, we achieve $\mathcal{O}\left(\omega^{-\frac{n+1}{2}}\right)$ error in total. Thus, the UNSD scheme transitions from $\mathcal{O}\left(\omega^{-\frac{2n+1}{3}}\right)$ behaviour for small $c$ to
$\mathcal{O}\left(\omega^{-\frac{n+1}{2}}\right)$ behaviour for large $c$. The error for large $c$ is larger, at least asymptotically, and this effect was clearly visible in Fig.~\ref{fig:unsd_vs_hermite}.

\section{Construction of the quadrature rule}
\label{sect_construction}

The Gaussian quadrature rule in this paper is unconventional because it has points in the complex plane. Unfortunately, this means that standard methods for the computation of Gaussian quadrature rules have to be amended. Two issues we address here are the stable numerical computation of recurrence coefficients, and the accurate computation of quadrature weights from the Jacobi matrix eigenvalue problem. We start with the latter, in some detail in order to properly explain the modification.

\subsection{The complex-symmetric Jacobi matrix}

Usually, for real-valued problems, a Gaussian quadrature rule is found from an eigenvalue problem involving the symmetric Jacobi matrix. This is a tridiagonal matrix that is defined in terms of the three-term recurrence coefficients of the orthogonal polynomials. In this context, the conventional form in which the three-term recurrence relation is expressed (for the \emph{monic} polynomials $p_k$) is
\begin{equation}
\label{eq:recurrence_monic}
 p_{k+1}(x) = (x - \alpha_k) p_k(x) - \beta_k p_{k-1}(x).
\end{equation}
Here, the recurrence coefficients are given by the standard formulas
\begin{equation}\label{eq:alpha_beta}
 \alpha_k = \frac{ (x p_k, p_k) }{(p_k,p_k) } \qquad \mbox{and} \qquad \beta_k = \frac{(p_k,p_k)}{(p_{k-1},p_{k-1})},
\end{equation}
where the bilinear form $(f,g)$ is given by
\begin{equation}\label{eq:bilinear_form}
 (f,g) = \int_\Gamma f(z) g(z) e^{i \left(\frac{z^3}{3}-\delta z \right) } {\rm d}z.
\end{equation}
The complex-valued orthogonal polynomials satisfy a three-term recurrence relation because the orthogonality condition \eqref{eq:orthogonality} is non-hermitian, i.e., we have $(x p, q) = (p, x q)$, rather than $(x p, q) = \overline{(p, x q)}$. Hence, the standard formulas remain valid in the complex case.

Note that the bilinear form \eqref{eq:bilinear_form} is not an inner product, due to the oscillatory weight function, hence the problem of existence of $p_n$.

The corresponding recurrence relation for the normalized polynomials $\pi_k(x)$, normalized in the sense that $(\pi_k,\pi_k) = 1$, is \cite{gautschi2004opq},
\begin{equation}
\label{eq:recurrence_normalized}
 \sqrt{\beta_{k+1}} \pi_{k+1}(x) = (x - \alpha_k) \pi_k(x) - \sqrt{\beta_k} p_{k-1}(x),
\end{equation}
where we set $\beta_0 = \mu_0$. This leads to the complex-symmetric Jacobi matrix,
\begin{align*}
J_n=\begin{pmatrix}
\alpha_0 &\sqrt{\beta_1} &0&0&\ldots &0\\
\sqrt{\beta_1} & \alpha_1 &\sqrt{\beta_2}  &0&\ldots &0\\
0& \sqrt{\beta_2} &\alpha_2 & \sqrt{\beta_3}  &\ldots &0\\
\vdots &\vdots &\vdots & \ddots &\ddots &\vdots\\
0& 0&0&\ldots &\alpha_{n-1}&\sqrt{\beta_n} \\
0& 0&0&\ldots &\sqrt{\beta_n} &\alpha_n
\end{pmatrix}.
\end{align*}
The roots of the polynomial are the eigenvalues of this matrix. The matrix is complex-symmetric, i.e. it is symmetric ($J_n = J_n^T$) but not hermitian ($J_n = J_n^*$). One can show from the string equations (see \eqref{eq:string_beta}-\eqref{eq:string_alpha} further below) that $\alpha_k$ is purely imaginary, and $\beta_k$ is real and positive, hence the square root is well defined \cite{lejon2016phd}.

\subsection{Eigenvalue decomposition and quadrature weights}

The eigenvector $\mathbf{v}_j$ corresponding to eigenvalue $x_j$ is a normalized vector, where each entry is proportional to $\pi_i(x_j)$:
\begin{equation}\label{eq:eigenvector}
 \mathbf{v}_j = c \left[\begin{array}{c}
                   \pi_0(x_j) \\
                   \pi_1(x_j) \\
                   \ldots \\
                   \pi_{n-1}(x_j)
                  \end{array}\right].
\end{equation}
The constant $c$ is such that $\Vert \mathbf{v}_j \Vert = 1$. In the complex-valued case, this means that
\begin{equation}\label{eq:c}
 c = \left( \frac{1}{\sum_{k=0}^{n-1} \pi_k(x_j) \overline{\pi_k(x_j)}} \right)^{1/2}.
\end{equation}
An explicit formula for the corresponding quadrature weight $w_j$ is \cite{gautschi2004opq}
\begin{equation}\label{eq:weights}
 w_j = \frac{1}{\sum_{k=0}^{n-1} \pi_k(x_j)^2},
\end{equation}
involving as above the orthonormal polynomials given by \eqref{eq:recurrence_normalized}.\footnote{Note that one can not in general rule out division by zero in this expression for the weights: the sum of squares in the denominator is not necessarily strictly positive in the complex case, as it is in the real case. An alternative and generally applicable expression for the weights is the weighted integral of the Lagrangian basis polynomials $l_{n,j}(t)$ for polynomial interpolation: $w_j = \int_\Gamma l_{n,j}(t) e^{i g(t,\delta)} {\rm d}t$. Hence, the weights are well defined as long as the Lagrangian polynomials are well defined, and this is the case if all roots $t_{k,\delta}$ are distinct. An issue may arise if the orthogonal polynomial $p_{n,\delta}$ has roots of higher multiplicity. In that case, however, a Gaussian-like quadrarure rule would still exist based on Hermite interpolation using derivatives, rather than Lagrangian interpolation. Strictly speaking, we have not excluded that case for the polynomials of this paper. Yet, coalescence of roots was never observed in experiments.}

In the real-valued case, we find the simple relation $w_j = c^2$. Moreover, from $\pi_0(x) = \sqrt{\mu_0}$, we can deduce $c$ and, hence, compute the quadrature weight from the first entry of the eigenvector $\mathbf{v}_j$:
\begin{equation}\label{eq:first_entry}
 w_j = c^2 = \mu_0 v_{j1}^2.
\end{equation}
This is a standard step in the well-known Golub--Welsch algorithm, which has $O(n^2)$ computational complexity overall by exploiting the tridiagonal structure of the Jacobi matrix \cite{golub1969gauss}. The final step \eqref{eq:first_entry} needs modification in the complex case, because the expressions for $w_j$ and $c^2$ are no longer the same due to the appearance of complex conjugates in \eqref{eq:c}.

In the complex case, we can still determine $c = \sqrt{\mu_0} v_{j1}$ from the first entry of the eigenvector $\mathbf{v}_j$. Hence, from \eqref{eq:eigenvector}, we can find the values $\pi_k(x_j)$ and we can explicitly evaluate the sum in \eqref{eq:weights} to find the weight $w_j$. This does require the full eigenvalue decomposition of the Jacobi matrix and leads to $O(n^3)$ computational complexity.

Alternatively, once the points $x_j$ have been determined (in $O(n^2)$ operations as with the Golub--Welsch algorithm), one can evaluate the orthonormal polynomials $\pi_k(x_j)$ using their recurrence relation \eqref{eq:recurrence_normalized}. This requires only $O(n^2)$ operations, and so does the explicit summation of \eqref{eq:weights} for the $n$ quadrature weights. With this modification, the computational cost of the quadrature rule remains $O(n^2)$.

\subsection{Computation of the recurrence coefficients for varying $\delta$}

It remains to compute the recurrence coefficients, which is a complication because they depend on the parameter $\delta$. A popular strategy is the Stieltjes algorithm \cite{gautschi2004opq}, which is based on a discrete approximation to the bilinear form \eqref{eq:bilinear_form}. For example, we may use the trapezoidal rule on a parabolic contour $P$ that connects the right points at infinity. However, we found this to be numerically stable only up to moderate values of $n$. The main complication is the growth of the polynomials $p_{n,\delta}$ in the complex plane, which is similar to the growth of real-valued orthogonal polynomials when evaluated away from the real line.

The recurrence coefficients satisfy the so-called string equations, recall \eqref{eq:string}. In the notation of \eqref{eq:recurrence_monic}, these become
\begin{eqnarray}\label{eq:string2}
\beta_{n+1} &=\delta-\beta_n-\alpha_n^2, \label{eq:string_beta} \\
\alpha_{n+1} &=\frac{i(n+1)}{\beta_{n+1}}-\alpha_n. \label{eq:string_alpha}
\end{eqnarray}
These relations are very simple to evaluate. Unfortunately, these relations also quickly become numerically unstable for increasing $n$, with notable loss of accuracy beyond $n=20$.

On the other hand, assuming an accurate evaluation of the recurrence coefficients, the eigenvalue computation of the symmetric Jacobi matrix does appear to be accurate for larger values of $n$. Thus, one can evaluate the simple string equations in high-precision arithmetic, followed by an eigenvalue computation with standard floating point precision.

In practice, we found it convenient to approximate the recurrence coefficients a priori as a function of $\delta$, using Chebyshev expansions on an interval $[0,\delta_{\textrm{max}}]$. These expansions are computed once using the string equations \eqref{eq:string_beta}--\eqref{eq:string_alpha} in high-precision arithmetic, but the resulting expansions (for $\alpha_i(\delta)$ and $\beta_i(\delta)$) can be stored in standard precision. The quadrature rule for a given value of $\delta$ can then be computed on the fly using an efficient eigenvalue computation as described above. Assuming $N$ recurrence coefficients are precomputed, one can compute quadrature rules with any number of points between $1$ and $N$.

\section{Example: the computation of Chebyshev moments}
\label{sect_application}

\setlength{\figurewidth}{6.8cm}
\setlength{\figureheight}{6cm}

\begin{figure}[t]
\centering

\subfigure[$n=4$]{
\begin{small}
\input{relerror_chebyshev_n4.tikz}
\end{small}
}
\subfigure[$n=8$]{
\begin{small}
\input{relerror_chebyshev_n8.tikz}
\end{small}
}
\caption{Oscillatory moments of Chebyshev polynomials are computed using the UNSD quadrature rule, with polynomial degree $k$ and frequency $\omega$. Relative error is shown as a function of $k$, for several values of $m = \frac{k}{\omega}$. The known recurrence for the moments is unstable when $m \approx 1$, but UNSD is very accurate using only a very small number of points.}\label{fig:chebyshev}
\end{figure}

We conclude the paper with an example of recent interest in the literature on oscillatory integrals. We consider the evaluation of oscillatory moments of Chebyshev polynomials,
\[
 \int_a^b T_k(x)e^{i\omega x}dx
\]
where the polynomials have moderate to high degree $k$ and the frequency $\omega$ is large. These Chebyshev moments regularly appear in numerical methods for Fourier-type integrals of the form \eqref{eq:I} with $g(x)=x$, in which $f$ is approximated by a Chebyshev expansion \cite{dominguez2011filonclenshawcurtis,ledoux2012interpolatory,huybrechs2012superinterpolation}.

These moments satisfy a recursion, based on the recurrence relation of the Chebyshev polynomials themselves \cite{piessens1971oscillatory}. This recursion is numerically stable for small $k$ \cite{dominguez2011filonclenshawcurtis}. However, it becomes unstable once $k \sim \omega$. The remaining moments can then be computed by solving a linear system of equations \cite{dominguez2011filonclenshawcurtis}, which can become expensive.

The resonance in the regime $k \sim \omega$ can be reformulated as a problem with coalescing saddle points. This is seen by rewriting the Chebyshev moment as follows, using well-known properties of the Chebyshev polynomials of the first kind:
\begin{align*}
\int_{-1}^1T_k(x) e^{i\omega x}dx &= \int_0^\pi \sin t \, T_k(\cos t)e^{i\omega \cos t}dt
\\&=\int_{-1}^1\sin t \cos (kt)e^{i\omega \cos t}dt
\\&=\frac{1}{2}\int_{-1}^1 \sin t \, e^{i kt}e^{i\omega \cos t}dt+\frac{1}{2}\int_{-1}^1\sin t \, e^{-ikt}e^{i\omega \cos t}dt \\
&= I_1 + I_2.
\end{align*}
The integral is written as the sum of two integrals of the form \eqref{eq:I}, with oscillatory factors
\[
e^{\pm ikt}e^{i\omega \cos t}=e^{i\omega(\pm mt+\cos t)},
\]
where $m=\frac{k}{\omega}$.
It can be verified that the second oscillator does not have stationary points in the interval $[0,\pi]$. The first oscillator, $mt+\cos t$, has a stationary point whenever
\[
\sin t= m.
\]
This equation has two solutions close to $\frac{\pi}{2}$ when $m\approx 1$. This implies that the second integral can be approximated by the simpler NSD-method, and the first by the UNSD-algorithm.

In order to apply UNSD, the oscillator of the second integral is mapped to a cubic oscillator of the form $\frac{t^3}{3}-Kt$ by a change of variables. Since this is a non-linear problem, we have used a scheme based on Newton-Raphson iteration to compute this change of variables numerically.

The case of particular practical interest here arises when $m = \frac{k}{\omega} \approx 1$. Figure \ref{fig:chebyshev} illustrates the accuracy of applying the UNSD scheme to integral $I_1$. It turns out that the method is very accurate precisely in the regime $m\approx 1$, for large $k$, using just $4$ quadrature points. Increasing the quadature points to $n=8$ leads to high relative accuracy even for moderately small values of $k$ (and correspondingly small values of $\omega = \frac{k}{m}$).

\section*{Acknowledgements}
All authors were supported by FWO Flanders project G.0641.11. In addition, Arno Kuijlaars is supported by a long term structural funding-Methusalem grant of the Flemish Government, and via support of FWO Flanders through projects G.0864.16 and EOS G0G9118N. Daan Huybrechs had additional support from FWO Flanders project G.A004.14N.

\bibliographystyle{abbrv}
\bibliography{cubicosc}

\end{document}


TITLE: Back-and-forth harmonic progression -- do we get all of $\mathbb{Q}?$
QUESTION [4 upvotes]: I was attempting a problem earlier and this side question came about.
Set up:
We create a set of rationals $A$ by hopping along $\mathbb{Q}$. It's awkward to write out, but the idea is very simple, I attached a picture of the first couple steps to illustrate.

Start at $0\in A$
Jump to the right by $1$, add that value to set (in other words, $1\in A$)
Jump to left by $\frac12$, add value to set (i.e. $\frac{1}{2}\in A$), jump again to the left by $\frac13$, add value to set (i.e. $\frac16\in A$), jump to the left by $\frac14$ and so on.
... till we reach just before $-1$ (this happens after left-jumping $\frac{1}{10}$). Now continue to the right, following the progression i.e. to the right by $\frac{1}{11}$, then $\frac{1}{12}$ etc.
... till we reach just before $+2$. Now continue to the left, following the progression etc. till we reach just over $-2$, then back to right till just under $+3$ etc.

In other words, we are going back and forth over progressively larger intervals where the jumps are getting smaller, following a harmonic progression. Since $\sum 1/n$ diverges this never converges.
Now the resulting set $A$ does not contain all rationals, for instance it misses all integers other than $0, 1$. It seems messy and probably not very interesting to look at $A$ alone, but it led me to think about the following:
Question
If we label the above $A$ as $A_0$, and consider sets $A_k$ obtained in the same way but starting at the integer $k$, rather than $0$, does:
$$\bigcup_{k\in\mathbb{Z}} A_k = \mathbb{Q}?$$
And moreover, is $\left|A_i\cap A_j\right|\leqslant 1$ for $i\neq j$?
This appears to get messy rather quickly if one tries a hands-on approach, but I am having trouble getting a feel of what is happening. I wrote a program to get an idea but the numerators/denominators get large very quickly

REPLY [3 votes]: Whatever is meant by "in the same way but starting at the integer k" (I can think of at least 3 different ways to do "in the same way" and they don't give the same $A_k$s)
I'll give a proof of why I said it doesn't hit $\frac13$ or any integral translates of it.
Claim: The procedure does not hit any rational number whose lowest form has $3$ in the denominator.
Proof: Whatever "in the same way" may be, the procedure only hits rationals that are of the form
$$
k\pm\frac11\pm\frac12\pm\frac13\pm\dots\pm\frac1n
$$
We check that none of the signs, for $n\leq 4$, hits a rational with 3 in the denominator unless it is an integer.  For $n>4$, note that Bertrand's postulate gives a prime $p$ between $\lceil(n+1)/2\rceil$ and $n$ inclusive (in particular, $>3$ starting from $p=n=5$) and so its appearance in the denominator is guaranteed.  QED.
As to whether $\lvert A_i\cap A_j\rvert\leq 1$, that really depends on how you define $A_k$ "in the same way".  If you just mean translation by $k$ then it is true.  If you mean some other ways, e.g., starting at $k$ move $\pm 1$ to near $0$ (+ if at 0), move $\pm\frac12$ ... reversing only when you travelled more than say $(\lvert k\rvert+2)$ times the previous distance then I just defeated this. by having $1,0,\frac12\in A_0\cap A_1$.


TITLE: Function and graph
QUESTION [1 upvotes]: A relation is a subset of a carthesian product.
A function from $X$ to $Y$ is a relation with particular properties. It follows that a function is a subset of the cartesian product $X\times Y$.
The graph of f is
$$
\{(x,y) \in  X\times Y\ |\  y=f(x)\}.
$$
Then what is the difference between the graph and the function if both are the same subset of $X\times Y$?

REPLY [1 votes]: I'm with Terence Tao when he writes

functions are not sets, and sets are not functions; it does not make
sense to ask whether an object $x$ is an element of a function $f$,
and it does not make sense to apply a set $A$ to an input $x$ to
create an output $A(x)$.

Taking it slowly, note first that even if, for certain purposes, we implement a function as a set of ordered pairs, it is arbitrary (i) to  fix on e.g. Kuratowski's particular implementation of pairs as sets (there other perfectly workable options).  And even relative to a choice of set-theoretic pairing scheme, (ii) we could equally well implement a function using the set of pairs $\langle y, x\rangle$ where $f(x) = y$, rather than by the set of pairs $\langle x, y\rangle$ -- some textbooks do just this. (Again, other choices are also possible.)
However, if various permutations of choices at stages (i) and (ii) are pretty much as workable as each other, then we surely can't suppose that -- when we choose to implement a function as its graph as conventionally defined  -- we have made the uniquely right choice, i.e. the choice that correctly identifies which set that function really is.
But if there is no determinate fact of the matter about which sets functions are, then functions aren't sets.
It is worth digging a bit deeper. The key underlying point is that functions just aren't the right logical type of thing to be sets. As  the great logician Alonzo Church puts it:

it lies in the nature of any given [one-place] function to be
applicable to certain things and, when applied to one of them as
argument, to yield a certain value.

For example, a function such as the factorial defined over the natural numbers is, of its nature, the type of thing which yields a numerical value when given a number as argument. By contrast a set doesn't, of its nature, take an argument or yield a value, as Tao is pointing out; and what applies to sets in general applies to e.g. sets of ordered pairs of numbers in particular. (A well-known Fregean metaphor might help. Functions of their nature are 'unsaturated', having a number of empty slots waiting to be filled appropriately when the function is applied to the right number of arguments. By contrast, a set is already 'saturated', it is self-standing, with no empty slots waiting to be filled.)
Which isn't for a moment to deny that we can make use of the graph $\hat{f}$ (a glorified input-output look-up table) for the function $f$ in mapping an input object to an output value. But to do this, we need to deploy another function, namely a two-place evaluation function which takes an object $x$ and the graph $\hat{f}$, and outputs $y$ if and only if the pair $\langle x, y\rangle$ is in the graph. And unless we are planning to set off on an infinite regress, we had better not seek to  again trade in this evaluation function for another set.
So yes, for many purposes, appropriate sets (graphs) can be used as proxies for functions ('implement' functions, if you prefer). But that doesn't make functions sets. 'Everything is what it is and not another thing'.


TITLE: Given P(A) intersecting itself, when does the resulting probability equal P(A)^2 versus 1 - P(A')^2?
QUESTION [0 upvotes]: P(A) is the probability that something occurred. If it occurred again would the probability of both occurring (intersection) be P(A)^2 or would it be 1 - P(A')^2? What are the conditions that would distinguish one from the other?

REPLY [1 votes]: Suppose we have two independent trials, with an identical probability for an event happening in each particular trial.   Let $\mathsf P(A)$ be that probability.
Then the probability that the event occurs in both trials will be $\mathsf P(A)^2$, by the independence.
Whereas $(1-\mathsf P(A))^2$ would be the probability that the event does not occur in both trials.
Also $2\,\mathsf P(A)\,(1-\mathsf P(A))$ would be the probability for the event occuring exactly one trial.
And $1-(1-\mathsf P(A))^2$ the probability for at least one occurance among the two trials.


text:
\begin{document}

\title[Zeta Functions of Cyclic covers in Large Characteristic]{Computing Zeta Functions of Cyclic covers in Large Characteristic}

\begin{abstract}
  We describe an algorithm to compute the zeta function of a cyclic cover of the projective line over a finite field of characteristic $p$ that runs in time $p^{1/2 + o(1)}$.
  We confirm its practicality and effectiveness by reporting on the performance of our \textsc{SageMath} implementation on a range of examples.
  The algorithm relies on Gon\c{c}alves's generalization of Kedlaya's algorithm for cyclic covers, and Harvey's work on Kedlaya's algorithm for large characteristic.
\end{abstract}

\maketitle

\blfootnote{2010 {\normalfont\itshape
		Mathematics Subject Classification}: 11G20 (primary)
	11Y16, 11M38, 14G10 (secondary)}

\section{Introduction}
\label{sec:introduction.tex}
\input{introduction.tex}

\section{Setup and notation}
\label{sec:setup}
\input{setup.tex}

\section{The Frobenius action on differentials}
\label{sec:frob_action}
\input{frob_action.tex}

\section{Reducing differentials}
\label{sec:reduction}
\input{reduction.tex}

\section{Main algorithm} \label{sec:algorithm}
\input{complexity.tex}

\section{Sample Computations} \label{sec:sample_computations}
\input{examples.tex}

\bibliographystyle{alpha}
\bibliography{biblio}

\end{document}


TITLE: Are parenthesis required for subsitution?
QUESTION [0 upvotes]: Do we need to use parenthesis to substitute into an expression, for example if I have an expression like $2\div4(3)$ I might wish to replace the $4(3)$ with $4\times3$ it is doable without parameters, so I simply replace one set of symbols with another, giving $2\div4\times3$, however if we simply replace one symbol with another it is misleading as there is no actual operation on the value $4\times3$ (there is a multiplication between ($2\div4$) and $3$) it is just that there appears a $\times$ between '$4$' and the $'3'$, is the parenthesis required to avoid this?

REPLY [1 votes]: Addition and multiplication have the nice property that $(a+b)+c=a+(b+c)$ and $(ab)c=a(bc)$ called associativity. However both subtraction and division fail to have this property since in general $(a-b)-c \neq a-(b-c)$ and $(a/b)/c \neq a/(b/c)$. There is another property called commutativity which addition and multiplication have which means $a+b=b+a$ and $ab=ba$ but again, subtraction and division fail to have this property.
To avoid this problem in higher math we just get rid of subtraction by adding negative numbers and division by multiplying reciprocals. So this means when I write $a-b$ I really mean $a+-b$ but it seems excessive to add the extra symbol all the time. The same happens for multiplication too where $a/b$ means $a(1/b)$ instead but again, we don't want all those extra symbols so we use the shorter version. In practice this is a very good way to avoid errors so it's worth thinking about things from this point of view.