# Algebras of commuting differential operators
for integral kernels of Airy type

W. Riley Casper  Dept of Mathematics
California State University Fullerton
Fullerton, CA 92831
U.S.A<EMAIL_ADDRESS>,  F. Alberto Grünbaum  Dept of Mathematics
University of California, Berkeley
Berkeley, CA 94720
U.S.A<EMAIL_ADDRESS>,  Milen Yakimov  Dept of Mathematics
Northeastern University
Boston, MA 02115
U.S.A<EMAIL_ADDRESS>and  Ignacio Zurrián  Dept of Mathematics
Universidad Nacional de Córdoba, Córdoba, X5016HUA, Argentina
<EMAIL_ADDRESS>
###### Abstract.

Differential operators commuting with integral operators were discovered in
the work of C. Tracy and H. Widom [37, 38] and used to derive asymptotic
expansions of the Fredholm determinants of integral operators arising in
random matrix theory. Very recently, it has been proved that all rational,
symmetric Darboux transformations of the Bessel, Airy, and exponential
bispectral functions give rise to commuting integral and differential
operators [6, 7, 8], vastly generalizing the known examples in the literature.
In this paper, we give a classification of the the rational symmetric Darboux
transformations of the Airy function in terms of the fixed point submanifold
of a differential Galois group acting on the Lagrangian locus of the (infinite
dimensional) Airy Adelic Grassmannian and initiate the study of the full
algebra of differential operators commuting with each of the integral
operators in question. We leverage the general theory of [8] to obtain
explicit formulas for the two differential operators of lowest orders that
commute with each of the level one and two integral operators obtained in the
Darboux process. Moreover, we prove that each pair of differential operators
commute with each other. The commuting operators in the level one case are
shown to satisfy an algebraic relation defining an elliptic curve.

##  1\. Introduction

Our contribution to this volume bears a connection with a phenomenon uncovered
by Craig Tracy and Harold Widom [38] in their work on level spacing in Random
Matrix Theory. For a double scaling limit at the “edge of the spectrum” they
observed that the resulting integral operator with the Airy kernel acting on
an appropriate interval admits a commuting second order differential operator.
This highly exceptional fact is put to good use in section IV of their paper
where a number of asymptotic results for several quantities of interest are
given.

In the context of Random Matrix Theory the existence of such a commuting pair
of operators had been exploited earlier, for instance in work by M. Mehta [29]
and W. Fuchs [14]. In this case one is interested in the “bulk of the
spectrum” and the role of the Airy kernel is taken up by the more familiar
sinc kernel. Both of these situations deal with the Gaussian Unitary Ensemble.

The consideration of either the Laguerre or the Jacobi ensembles at the “edge
of the spectrum” gives rise to the Bessel kernel. This case, as well as the
corresponding commuting pair of integral-differential operators is considered
by C. Tracy and H. Widom in [37]. There, once again, this exceptional fact is
exploited in section III to derive a number of important asymptotic results.

In this paper we concentrate on the “exceptional fact” mentioned above in
three different situations relevant to Random Matrix theory. This fact had
emerged in other areas of mathematics. In a ground-breaking collection of
papers by D. Slepian, H. Landau and H. Pollak done at Bell labs in the 1960’s
[27, 28, 32, 33, 34, 35, 36] instances of this phenomenon were discovered and
used in a key way in communication-signal processing theory. In fact, some
precedents can be traced further back, see [4, 22]. For an up-to-date
treatment of the numerical issues involving the prolate spheroidal function,
one can see [30].

Incidentally in the case of the Bessel kernel the existence of a commuting
operator was already proved by D. Slepian, while the situation of the Airy
kernel appears for the first time in C. Tracy and H. Widom’s paper mentioned
above. The so called “prolate spheroidal wave functions,” which arise in the
case of the sinc kernel and their corresponding integral-differential pair of
operators, have played an important role in areas far removed from signal
processing that motivated the research of Slepian and collaborators. We give
only two instances of this, but we are sure that other people can provide
other examples: the paper by J. Kiukas and R. Werner [24] in connection with
Bell’s inequalities, and the program by A. Connes in connection with the
Riemann hypothesis with C. Consani, M. Marcolli and H. Moscovici [10, 11, 12].

One should also mention that the Airy function itself and variants of it have
played an important role in other very active areas of current research, such
as quantum gravity and intersection theory on moduli space of curves, see [25,
41].

In all the three instances discussed above (namely the sinc, Bessel and Airy
kernels), the commuting differential operator has been found by a direct
computation that relies heavily on integration by parts. The interest in
understanding and extending this exceptional phenomenon in a variety of other
situations has produced some few more examples, see [5, 9, 15, 16, 17, 19,
20].

The bispectral problem formulated 1986 in [13] aimed at a conceptual
understanding of the phenomenon of integral operators admitting a commuting
differential operator. The idea is that all known kernels with this property
are built from bispectral functions, that is functions in two complex
variables that are eigenfunctions of differential operators in each of them.
There has been a substantial amount of research on this problem [18, 21],
which started with the classification of all bispectral differential operators
of second order [13] and culminated in the classification of bipectral
functions of rank 1 in [39] and the construction of bispectral functions of
arbitrary rank via Darboux transformations [2, 23] and automorphisms of the
first Weyl algebra [1, 3].

Since the mid 80s, the belief that the two problems, bispectrality and the
existence of a commuting pair made up of a differential and an integral
operator were closely connected has been driving research on both fronts.
However, for a long time there no general argument proving that bispectral
functions give rise to kernels of integral operators with the commutativity
property. This was finally settled in [8] where it was proved to be the case
for self-adjoint bispectral functions of rank 1 and 2.

More recently we proved in [6, 7] that all bispectral functions of rank 1 give
rise to integral operators that reflect a differential operator rather than
plain commute with it.

All of the previous results on integral operators address the construction of
a single differential operator commuting with it. The purpose of this paper is
to initiate the systematic study of the _algebras_ of differential operators
that commute with a given integral operator. We start with the Airy example
considered by C. Tracy and H. Widom and consider all self-adjoint bispectral
Darboux transformations. This is an infinite dimensional manifold which sits
canonically in the infinite dimensional Grassmannian of all Darboux
transformation from the Airy function, obtained from factorizations of
polynomials of the Airy operator

(1) $L(x,\partial_{x})=\partial^{2}-x.$

We give a conceptual classification of the former manifold as the fixed point
set of a Lagrangian Grassmannian with respect to the canonical action of the
associated differential Galois group. The Lagrangian Grassmannian in question
is the sub-Grassmannian with respect to a canonical symplectic form. We
consider the first two instances of self-adjoint bispectral Darboux
transformations coming from factorizations of

$(L-t_{1})^{2}\quad\mbox{and}\quad(L-t_{2})^{4}$

of the form $P^{*}P$ for a differential operator $P(x,\partial_{x})$ with
rational coefficients. The corresponding bispectral functions, referred to
here as level one and level two bispectral functions, are significantly more
complicated than the bispectral Airy function $\mathrm{Ai}(x+z)$. The integral
operators that they give rise to depend on parameters classifying different
factorizations. For each integral operator, we compute explicitly the
differential operators of the lowest two orders and prove that they are
algebraically dependent. In the level one situation, the commuting operators
have order $4$ and $6$. They generate the algebra of all differential
operators commuting with the integral operator and satisfy an algebraic
relation which happens to be an elliptic curve. In the level two situation,
the lowest two commuting operators have order $10$ and $12$. However, we are
also able to find commuting operators of order $14$, $16$, and $18$ and to
prove that these differential operators commute with each other. In a future
publication, we will return to the problem of studying algebras of
differential operators commuting with a fixed integral operator and will
present general structural results for the algebra of differential operator
commuting with all integral operators which are built from bispectral
functions, and which are motivated by the examples in this paper.

This paper is written as a small token of admiration and gratitude to the
amazing mathematical work of Harold Widom. Widom started mathematical life as
an algebraist working with Irving Kaplansky at Chicago, before becoming mainly
an analyst through the influence of Mark Kac at Cornell. This paper uses tools
from both analysis and algebra, uniting Widom’s dual mathematical history. His
influence will be a lasting one, and we will miss him badly.

##  2\. Bispectral functions, Fourier algebras and prolate spheroidal type
commutativity

###  2.1. Bispectrality and Fourier Algebras

For an open subset $U\subseteq\mathbb{C}$, denote by $\mathfrak{D}(U)$ the
algebra of differential operators on $U$ with meromorphic coefficients.

######  Definition 2.1.

[13] Let $U$ and $V$ be two domains in $\mathbb{C}$. A nonconstant meromorphic
function $\Psi(x,z)$ defined on $U\times V\subseteq\mathbb{C}^{2}$ is called
bispectral if there exist differential operators
$B(x,\partial_{x})\in\mathfrak{D}(U)$ and
$D(z,\partial_{z})\in\mathfrak{D}(V)$ such that

$\displaystyle B(x,\partial_{x})\Psi(x,z)$ $\displaystyle=g(z)\Psi(x,z)$
$\displaystyle D(z,\partial_{z})\Psi(x,z)$ $\displaystyle=f(x)\Psi(x,z)$

for some nonconstant functions $f(x)$ and $g(z)$ meromorphic on $U$ and $V$,
respectively.

Denote by $\operatorname{{\mathrm{Ai}}}(x)$ the classical Airy function. The
function

$\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z):=\operatorname{{\mathrm{Ai}}}(x+z)$

is bispectral because

(2)
$L(x,\partial_{x})\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)=z\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)\quad\mbox{and}\quad
L(z,\partial_{z})\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)=x\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z),$

where $L(x,\partial_{x})$ is the Airy operator (1). The differential equations
satisfied by a bispectral function are captured by the following definition.

######  Definition 2.2.

[1] Let $\Psi(x,z)$ be a bispectral meromorphic function defined on $U\times
V\subseteq\mathbb{C}^{2}$. Define the left and right Fourier algebras of
differential operators for $\Psi$ by

$\displaystyle{\mathfrak{F}}_{x}(\Psi)=\\{R(x,\partial_{x})\in\mathfrak{D}(U):\,$
there exists a differential operator $S(z,\partial_{z})\in\mathfrak{D}(V)$
$\displaystyle\text{satisfying
$R(x,\partial_{x})\Psi(x,z)=S(z,\partial_{z})\Psi(x,z)$}\\}$

and

$\displaystyle{\mathfrak{F}}_{z}(\Psi)=\\{S(z,\partial_{z})\in\mathfrak{D}(V):\,$
there exists a differential operator $R(x,\partial_{x})\in\mathfrak{D}(U)$
$\displaystyle\text{satisfying
$R(x,\partial_{x})\Psi(x,z)=S(z,\partial_{z})\Psi(x,z)$}\\}.$

By [8, Proposition 2.4], for every bispectral meromorphic function
$\Psi:U\times V\to\mathbb{C}$, there exists a canonical anti-isomorphism

$b_{\Psi}:{\mathfrak{F}}_{x}(\Psi)\to{\mathfrak{F}}_{z}(\Psi),$

given by $b_{\Psi}(R(x,\partial_{x}))=S(z,\partial_{z})$ if

$R(x,\partial_{x})\Psi(x,z)=S(z,\partial_{z})\Psi(x,z).$

We call this the generalized Fourier map associated to $\Psi(x,z)$. Define the
co-order of an element $R(x,\partial_{x})\in{\mathfrak{F}}_{x}(\Psi)$ by

$\operatorname{\mathrm{cord}}R:=\operatorname{\mathrm{ord}}b_{\Psi}(R).$

Analogously, we define the co-order of
$S(z,\partial_{z})\in{\mathfrak{F}}_{z}(\Psi)$ by
$\operatorname{\mathrm{cord}}S:=\operatorname{\mathrm{ord}}b_{\Psi}^{-1}(S)$.
The Fourier algebras of $\Psi(x,z)$ have natural
$\mathbb{N}\times\mathbb{N}$-filtrations:

$\displaystyle{\mathfrak{F}}_{x}(\Psi)^{\ell,m}=\\{R(x,\partial_{x})\in{\mathfrak{F}}_{x}(\Psi):\,\operatorname{\mathrm{ord}}R\leq\ell,\operatorname{\mathrm{cord}}R\leq
m\\},$
$\displaystyle{\mathfrak{F}}_{z}(\Psi)^{m,\ell}=\\{S(z,\partial_{z})\in{\mathfrak{F}}_{z}(\Psi):\,\operatorname{\mathrm{ord}}S\leq
m,\operatorname{\mathrm{cord}}S\leq\ell\\},$

where $\mathbb{N}=\\{0,1,\ldots\\}$ and
$b_{\Psi}({\mathfrak{F}}_{x}(\Psi)^{\ell,m})={\mathfrak{F}}_{z}(\Psi)^{m,\ell}$.
The commutative algebras

${\mathfrak{B}}_{x}(\Psi):=\bigcup_{\ell\geq
0}{\mathfrak{F}}_{x}(\Psi)^{\ell,0}\quad\mbox{and}\quad{\mathfrak{B}}_{z}(\Psi):=\bigcup_{m\geq
0}{\mathfrak{F}}_{z}(\Psi)^{0,m}$

are precisely the algebras of differential operators in $x$ and $z$,
respectively, for which $\Psi(x,z)$ is a eigenfunction.

######  Example 2.3.

The Airy bispectral function $\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)$
satisfies

$\displaystyle L(x,\partial_{x})\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)$
$\displaystyle=z\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z),$
$\displaystyle\partial_{x}\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)$
$\displaystyle=\partial_{z}\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z),$
$\displaystyle x\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)$
$\displaystyle=L(z,\partial_{z})\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z).$

The Fourier algebras ${\mathfrak{F}}_{x}(\Psi_{\operatorname{{\mathrm{Ai}}}})$
and ${\mathfrak{F}}_{z}(\Psi_{\operatorname{{\mathrm{Ai}}}})$ coincide with
the first Weyl algebra in the variables $x$ and $z$, respectively, and the
generalized Fourier map $b_{\Psi_{\operatorname{{\mathrm{Ai}}}}}$ is the anti-
isomorphism from the first Weyl algebra in $x$ to the first Weyl algebra in
$z$ given by

$b_{\Psi_{\operatorname{{\mathrm{Ai}}}}}(x)=\partial_{z}^{2}-z,\;\;b_{\Psi_{\operatorname{{\mathrm{Ai}}}}}(\partial_{x})=\partial_{z}.$

Furthermore,

$\dim{\mathfrak{F}}_{x}(\Psi_{\operatorname{{\mathrm{Ai}}}})^{2\ell,2m}=\ell
m+\ell+m+1,$

see [8, Sect. 3.1 and Lemma 5.5]. On the level of Wilson’s adelic
grassmannian, the anti-isomorphism $b_{\psi}$ is equivalent to Wilson’s
bispectral involution [39]. More generally, every anti-automorphism of the
first Weyl algebra determines a bispectral function as proved in [3].

######  Definition 2.4.

A rational Darboux transformation from the bispectral Airy function
$\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)$ is a function of the form

(3)
$\Psi(x,z):=\frac{P(x,\partial_{x})\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)}{q(z)p(x)}$

such that

(4)
$\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)=Q(x,\partial_{x})\frac{\Psi(x,z)}{\widetilde{q}(z)\widetilde{p}(x)}$

for some differential operators $P$ and $Q$ with polynomial coefficients and
polynomials $p(x)$, $\widetilde{p}(x)$, $q(z)$ and $\widetilde{q}(z)$ with
coefficients in $\mathbb{C}$. We define the bidegree of such a transformation
to be the pair $(\operatorname{\mathrm{ord}}P,\operatorname{\mathrm{cord}}P)$.

In this setting we have
$Q,P\in{\mathfrak{F}}_{x}(\Psi_{\operatorname{{\mathrm{Ai}}}})$,
$\widetilde{p}(x),p(x)\in{\mathfrak{F}}_{x}(\Psi_{\operatorname{{\mathrm{Ai}}}})^{0,m}$
and
$\widetilde{q}(x),q(x)\in{\mathfrak{F}}_{z}(\Psi_{\operatorname{{\mathrm{Ai}}}})^{0,\ell}$
for some $\ell,m\in\mathbb{N}$. Furthermore, eqs. (3)–(4) imply that

$Q(x,\partial_{x})\frac{1}{\widetilde{p}(x)p(x)}P(x,\partial_{x})\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)={\widetilde{q}(z)q(z)}\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z),$

and thus by Example 2.3,

(5)
$Q(x,\partial_{x})\frac{1}{\widetilde{p}(x)p(x)}P(x,\partial_{x})=\widetilde{q}(L(x,\partial_{x}))q(L(x,\partial_{x})).$

######  Theorem 2.5.

[1, 3, 23] All rational Darboux transformations of the bispectral Airy
function $\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)$ are bispectral functions.
More precisely, if $\Psi(x,z)$ is as in Definition 2.4, then it satisfies the
spectral equations

$\displaystyle\frac{1}{p(x)}P(x,\partial_{x})Q(x,\partial_{x})\frac{1}{\widetilde{p}(x)}\Psi(x,z)$
$\displaystyle=q(z)\widetilde{q}(z)\Psi(x,z),$
$\displaystyle\frac{1}{q(z)}b_{\Psi_{\operatorname{{\mathrm{Ai}}}}}(P)(z,\partial_{z})b_{\Psi_{\operatorname{{\mathrm{Ai}}}}}(S)(z,\partial_{z})\frac{1}{\widetilde{q}(z)}\Psi(x,z)$
$\displaystyle=p(x)\widetilde{p}(x)\Psi(x,z).$

###  2.2. Prolate Spheroidal Type Commutativity

A rational Darboux transformation $\Psi(x,z)$ of the bispectral Airy function
of bidegree $(d_{1},d_{2})$ is called self-adjoint if it has a presentation as
in Definition 2.4 such that

$Q(x,\partial_{x})=P^{*}(x,\partial_{x})$

and $\widetilde{p}(x)=p(x)$, $\widetilde{q}(z)=q(z)$. Here $P\mapsto P^{*}$
denotes the formal adjoint. It follows from (5) that $P$ has even order. A
rational Darboux transformation $\Psi(x,z)$ of the Airy bispectral function
$\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)$ is self-adjoint if and only if the
spectral algebras ${\mathfrak{B}}_{x}(\Psi)$ and ${\mathfrak{B}}_{z}(\Psi)$
are preserved under the formal adjoint, and this condition is satisfied if and
only if $\Psi(x,z)$ is an eigenfunction of nonconstant, formally symmetric
differential operators in $x$ and $z$ (i.e., operators that are fixed by the
formal adjoint), see [8, Remark 3.17 and Proposition 3.18].

For self-adjoint rational Darboux transformations $\Psi(x,z)$ of
$\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)$, both Fourier algebras
${\mathfrak{F}}_{x}(\Psi)$ and ${\mathfrak{F}}_{z}(\Psi)$ are preserved under
the formal adjoint and

(6) $(b_{\Psi}(R))^{*}=b_{\Psi}(R^{*})\quad\mbox{for all}\quad
R\in{\mathfrak{F}}_{x}(\Psi),$

see [8, Proposition 3.24 and 3.25]. Define

${\mathfrak{F}}_{x,\operatorname{\mathrm{sym}}}(\Psi):=\\{R\in{\mathfrak{F}}_{x}(\Psi):\,R^{*}=R\\}.$

By (6) for all $R\in{\mathfrak{F}}_{x,\operatorname{\mathrm{sym}}}(\Psi)$,

$(b_{\Psi}(R))^{*}=b_{\Psi}R.$

######  Example 2.6.

[8, Lemma 5.5] For all $\ell,m\in\mathbb{N}$,
${\mathfrak{F}}_{x,\operatorname{\mathrm{sym}}}^{2\ell,2m}(\Psi_{\operatorname{{\mathrm{Ai}}}})$
has a basis given by

$\\{L(x,\partial_{x})^{j}x^{k}+x^{k}L(x,\partial_{x})^{j}:0\leq j\leq l,\
0\leq k\leq m\\},$

and in particular,

${\mathfrak{F}}_{x,\operatorname{\mathrm{sym}}}^{2\ell,2m}(\Psi_{\operatorname{{\mathrm{Ai}}}})=(\ell+1)(m+1).$

For $\epsilon>0$ consider the sector

$\Sigma_{\epsilon}=\\{re^{i\theta}\in\mathbb{C}:r>0,\
|\theta|<\pi/6-\epsilon\\}.$

The Airy function $\operatorname{{\mathrm{Ai}}}(x)$ of the first kind is
holomorphic on this domain and has the asymptotic expansion

$\operatorname{{\mathrm{Ai}}}(x)=e^{-\frac{2}{3}x^{3/2}}\Bigg{(}\sum_{j=1}^{\infty}c_{j}x^{-j/4}\Bigg{)}$

for some $c_{j}\in\mathbb{R}$ where $x^{1/4}$ is interpreted as the principal
$4$th root of $x$. Furthermore, any rational Darboux transformation of
$\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)$ equals
$\Psi(x,z)=\frac{1}{p(x)q(z)}P(x,\partial_{x})\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)$
for some polynomials $p(x),q(z)$ and a differential operator
$P(x,\partial_{x})$ with polynomial coefficients. Thus, for any bispectral
Darboux transformation of $\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)$ we have
the asymptotic estimate

$\|\partial_{x}^{j}\partial_{z}^{k}\Psi(x,z)\|=e^{-\frac{2}{3}(x+z)^{3/2}}\mathcal{O}((|x|+|z|)^{(j+k)/2+m})$

on $\Sigma_{\epsilon}$ for some integer $m$. The transformation
$z\mapsto(2/3)z^{3/2}$ sends $\Sigma_{\epsilon}$ to the sector
$\\{re^{i\theta}\in\mathbb{C}:r>0,\ |\theta|<\pi/4-3\epsilon/2\\}$. Therefore
if $\Gamma_{1},\Gamma_{2}\subseteq\Sigma_{\epsilon}$ are smooth, semi-infinite
curves inside this domain with parametrizations
$\gamma_{i}(t):[0,\infty)\rightarrow\mathbb{C}$ then the real part of
$-2(\gamma_{1}(t)+\gamma_{2}(s))^{3/2}/3$ will go to $-\infty$ as
$t\rightarrow\infty$ or $s\rightarrow\infty$. The above asymptotic estimate
now shows that $\Psi(x,z)$ satisfies

$\int_{\Gamma_{1}}|x^{m}z^{n}\partial_{x}^{j}\partial_{z}^{k}\Psi(x,z)|dx\in
L^{\infty}(\Gamma_{2})\quad\mbox{and}\quad\int_{\Gamma_{2}}|x^{m}z^{n}\partial_{x}^{j}\partial_{z}^{k}\Psi(x,z)|dz\in
L^{\infty}(\Gamma_{1}),$

for every pair of smooth, semi-infinite curves
$\Gamma_{1},\Gamma_{2}\subseteq\Sigma_{\epsilon}$.

Recall that the bilinear concomitant of a differential operator

$R(x,\partial_{x})=\sum_{j=0}^{m}d_{j}(x)\partial_{x}^{j}.$

is the bilinear form $\mathcal{C}_{R}(-,-;p)$ defined on pairs of functions
$f(x),g(x)$, which are analytic at $p\in\mathbb{C}$ by

$\displaystyle\mathcal{C}_{R}(f,g;p)$
$\displaystyle=\sum_{j=1}^{m}\sum_{k=0}^{j-1}(-1)^{k}f^{(j-1-k)}(x)(d_{j}(x)g(x))^{(k)}|_{x=p}$
$\displaystyle=\sum_{j=1}^{m}\sum_{k=0}^{j-1}\sum_{\ell=0}^{k}\binom{k}{\ell}(-1)^{k}f^{(j-1-k)}(x)d_{j}(x)^{(k-\ell)}g(x)^{(\ell)}|_{x=p}.$

See for example [22, Chapter 5, Section 3].

######  Theorem 2.7.

[8] Let $\Psi(x,z)$ be a self-adjoint bispectral Darboux transformation of the
Airy bispectral function $\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)$ of
bidegree $(d_{1},d_{2})$ and let $\Gamma_{1}$ and $\Gamma_{2}$ be two semi-
infinite, smooth curves in $\Sigma_{\epsilon}$ for some $\epsilon>0$, whose
finite endpoints are $t_{1}$ and $t_{2}$, respectively. Assume moreover that
$\Psi(x,z)$ is holomorphic in a neighborhood of $\Gamma_{1}\times\Gamma_{2}$
and that the operators in ${\mathcal{F}}_{x}(\Psi)$ and
${\mathcal{F}}_{z}(\Psi)$ have holomorphic coefficients in a neighborhood of
$\Gamma_{1}$ and $\Gamma_{2}$, respectively. Then the following hold:

  1. (1)

$\dim{\mathfrak{F}}_{x,\operatorname{\mathrm{sym}}}^{2\ell,2m}(\Psi)\geq(\ell+1)(m+1)+1-d_{1}d_{2}$.

  2. (2)

If a differential operator
$S(z,\partial_{z})\in{\mathfrak{F}}_{z,\operatorname{\mathrm{sym}}}(\Psi)$
satisfies

$\mathcal{C}_{S}(-,-;t_{1})\equiv
0\quad\mbox{and}\quad\mathcal{C}_{b_{\Psi}^{-1}(S)}(-,-;t_{2})\equiv 0,$

then it commutes with the integral operator

$T:f(z)\mapsto\int_{\Gamma_{1}}K(z,w)f(w)dw\quad\text{with kernel}\quad
K(z,w)=\int_{\Gamma_{2}}\Psi(x,z)\Psi(x,w)dx.$

  3. (3)

If
$\dim\mathfrak{F}_{z,\operatorname{\mathrm{sym}}}^{2\ell,2\ell}\geq\ell(\ell+1)+2$,
in particular if $\ell=d_{1}d_{2}$, then there exists a differential operator
$S(z,\partial_{z})\in\mathfrak{F}_{z,\operatorname{\mathrm{sym}}}^{\ell,\ell}(\Psi)$
of positive order satisfying the assmption and conclusion in part (2).

As a special case of this theorem, we are able to recover the commuting
integral and differential operators studied by Tracy and Widom in [38]. In
particular, if we take $\Psi=\Psi_{\operatorname{{\mathrm{Ai}}}}$ in the
theorem, then it guarantees the existence of a differential operator of order
$2$ commuting with the integral operator

$T_{\operatorname{{\mathrm{Ai}}}}f(z)\mapsto\int_{t_{1}}^{\infty}K_{\operatorname{{\mathrm{Ai}}}}(z,w)f(w)dw,$

with kernel

$K_{\operatorname{{\mathrm{Ai}}}}(z,w)=\int_{t_{2}}^{\infty}\operatorname{{\mathrm{Ai}}}(x+z)\operatorname{{\mathrm{Ai}}}(x+w)dx=\frac{\operatorname{{\mathrm{Ai}}}^{\prime}(t_{2}+z)\operatorname{{\mathrm{Ai}}}(t_{2}+w)-\operatorname{{\mathrm{Ai}}}(t_{2}+z)\operatorname{{\mathrm{Ai}}}^{\prime}(t_{2}+w)}{z-w}.$

Solving the associated system of linear equations for the vanishing
concomitant, we discover that the differential operator

$S_{\operatorname{{\mathrm{Ai}}}}(z,\partial_{z}):=\partial_{z}(t_{1}-z)\partial_{z}+(t_{2}-t_{1})z+z^{2}$

satisfies the condition that
$\mathcal{C}_{S_{\operatorname{{\mathrm{Ai}}}}}(f,g;t_{1})=0$ for all
functions $f,g$ analytic at $t_{1}$. Its preimage under the generalized
Fourier map

$b_{\Psi}^{-1}(S_{\operatorname{{\mathrm{Ai}}}}(z,\partial_{z}))=\partial_{x}(t_{2}-x)\partial_{x}+(t_{1}-t_{2})x+x^{2}$

also satisfies the condition
$\mathcal{C}_{b_{\Psi}^{-1}S_{\operatorname{{\mathrm{Ai}}}}}(f,g;t_{2})=0$ for
all functions $f,g$ analytic at $t_{2}$. Therefore the differential operator
$S_{\operatorname{{\mathrm{Ai}}}}(z,\partial_{z})$ commutes with
$T_{\operatorname{{\mathrm{Ai}}}}$. This is precisely the differential
operator discovered by Tracy and Widom in [38].

##  3\. Classification of self-adjoint rational Darboux transformations of the
bispectral Airy function

In this section, we will classify the self-adjoint rational Darboux
transformations of the bispectral Airy function by leveraging two tools: (1)
the technology of differential Galois theory, and (2) the classification of
self-adjoint Darboux transformations in terms of Lagrangian subspaces of
symplectic vector spaces found in [8]. A similar classification is performed
in [2] using the entirely different technique of performing an explicit
asymptotic analysis of Wronskians associated to subspaces of the kernel. More
explicity, in this section we wish to classify factorizations of the form

(7)
$P(x,\partial_{x})^{*}\frac{1}{p(x)^{2}}P(x,\partial_{x})=q(L(x,\partial_{x}))^{2}$

where here $p$ and $q$ are polynomials and $P(x,\partial_{x})$ is a
differential operator with polynomial coefficients. Without loss of
generality, we take $q(z)$ to be monic so that $p(x)$ is the leading
coefficient of the operator $P(x,\partial_{x})$. The associated self-adjoint
rational Darboux transformation of the bispectral function
$\operatorname{{\mathrm{Ai}}}(x+z)$ is then defined by

$\Psi(x,z)=\frac{1}{p(x)q(z)}P(x,\partial_{x})\cdot\operatorname{{\mathrm{Ai}}}(x+z).$

###  3.1. Lagrangian Subspaces and Concomitant

We begin by recalling the classification of self-adjoint factorizations of
self-adjoint differential operators found in [8]. To begin, let
$A(x,\partial_{x})$ be a differential operator and recall the standard fact
that the concomitant $\mathcal{C}_{A}(f,g;x)$ is independent of $x$ for all
$f\in\ker(A)$ and $g\in\ker(A^{*})$.

######  Lemma 3.1 ([40], Section 3).

Let $A(x,\partial_{x})$ be a linear differential operator. Then the
concomitant of $A$ defines a canonical nondegenerate pairing

$\ker(A)\times\ker(A^{*})\rightarrow\mathbb{C},\
(f,g)\mapsto\mathcal{C}_{A}(f,g).$

Combining this with the identity
$\mathcal{C}_{A}(f,g)=-\mathcal{C}_{A^{*}}(g,f)$, we see that the concomitant
restricts to a symplectic bilinear form on $\ker(A)$ when $A(x,\partial_{x})$
is formally symmetric.

We will also rely on the following formula for concomitants of differential
operator products.

######  Lemma 3.2 ([40], Lemma 3.6).

Let $A(x,\partial_{x})=A_{1}(x,\partial_{x})A_{2}(x,\partial_{x})$. Then

$\mathcal{C}_{A}(f,g;x)=\mathcal{C}_{A_{1}}(A_{2}f,g;x)+\mathcal{C}_{A_{2}}(f,A_{1}^{*}g;x).$

From this, we see that if $A=A^{*}$, then $\ker(A_{2})\subseteq\ker(A)$ and
$\ker(A_{1}^{*})\subseteq\ker(A)$ are orthogonal under the pairing defined by
the concomitant of $A(x,\partial_{x})$.

As is well-known in the theory of factorizations of linear differential
operators, a factorization of a differential operator

$A(x,\partial_{x})=A_{1}(x,\partial_{x})A_{2}(x,\partial_{x})$

corresponds to a choice of a subspace $V\subseteq\ker(A)$. The subspace $V$
corresponds to the kernel of $A_{2}(x,\partial_{x})$ and determines the value
of the operator $A_{2}(x,\partial_{x})$ up to a left multiple by a function of
$x$. As is readily seen from the previous lemma, the kernel of
$A_{1}(x,\partial_{x})^{*}$ is completely determined by $V$ and given by the
orthogonal complement

$V^{\perp}=\\{g\in\ker(A^{*}):\mathcal{C}_{A}(f,g)=0\ \forall f\in V\\}.$

Thus to obtain factorizations of the form (7), we search in particular for
subspaces $V\subseteq\ker(q(L)^{2})$ satisfying $V^{\perp}=V$. In other words,
we search for Lagrangian subspaces of the symplectic vector space
$\ker(q(L)^{2})$. To summarize, we have the following proposition.

######  Proposition 3.3.

Factorizations of the form (7) with $p(x)$ and the coefficients of
$P(x,\partial_{x})$ not necessarily rational functions, correspond precisely
to Lagrangian subspaces of the symplectic vector space $\ker(q(L)^{2})$ whose
symplectic form is defined by the concomitant of $q(L)^{2}$.

###  3.2. Differential Galois Theory

Our next task is to determine the symmetric factorizations obtained in the
previous section which are rational. For the convenience of the reader, we
briefly outline the requisite basics of Picard-Vessiot extensions and the
Fudamental Theorem of Differential Galois Theory. We direct the interested
reader to [31] for a more thorough treatment.

######  Definition 3.4.

Let $(K,\partial)$ be a differential field and let $A\in K[\partial]$ be a
linear differential operator with coefficients in $K$. The Picard-Vessiot
extension of $K$ associated with $A(x,\partial_{x})$ is a differential field
extension $(F,\partial)$ of $K$ whose constants all belong to $K$ and which is
generated by the solutions of the homogeneous equation $Ag=0$.

Picard-Vessiot extensions of a differential field play precisely the role of
Galois extensions in field theory. Likewise, the usual Galois group is
replaced by a similar object consisting of field automorphisms respecting
differentiation.

######  Definition 3.5.

The differential Galois group $\operatorname{\mathrm{DGal}}(F/K)$ consists of
all $K$-linear field automorphisms $\sigma:F\rightarrow F$ of $F$ satisfying
$\sigma(\partial\cdot a)=\partial\cdot\sigma(a)$ for all $a\in F$.

Analogous to the case of Galois extensions of fields, we have the following
theorem relating differential subextensions and Zariski-closed subgroups of
the differential Galois group (see [31, Proposition 1.34]).

######  Theorem 3.6 (Fundamental Theorem of Differential Galois Theory).

Let $(K,\partial)$ be a differential field whose subfield of constants is
algebraically closed and let $(F,\partial)$ be a Picard-Vessiot extension of
$K$. Then there is a bijective correspondence between differential subfields
$K\subseteq F^{\prime}\subseteq F$ and Zariski-closed subgroups
$G^{\prime}\subseteq\operatorname{\mathrm{DGal}}(F/K)$ given by

$G^{\prime}\subseteq\operatorname{\mathrm{DGal}}(K/F)\mapsto
K^{G^{\prime}}=\\{a\in K:\sigma(a)=a,\ \forall\sigma\in G^{\prime}\\}.$
$K\subseteq F^{\prime}\subseteq
F\mapsto\operatorname{\mathrm{DGal}}(F^{\prime}/K)=\\{\sigma\in\operatorname{\mathrm{DGal}}(F/K):\sigma(a)=a,\
\forall a\in F^{\prime}\\}.$

Furthermore, this correspondence restricts to a correspondence between Picard-
Vessiot subextensions of $F/K$ and normal subgroups of
$\operatorname{\mathrm{DGal}}(F/K)$.

We will not rely on the full force of this correspondence, and therefore will
not have to recall the precise nature of the topological structure of
$\operatorname{\mathrm{DGal}}(F/K)$ as a group subscheme of a general linear
group. Instead, we will use only the immediate fact that

(8) $K=\\{a\in F:\sigma(a)=a,\
\forall\sigma\in\operatorname{\mathrm{DGal}}(F/K)\\}.$

Since differential operators are determined (up to a multiple) by their
kernels, rationality of a differential operator may be characterized by
differential Galois invariance of the associated kernel.

######  Theorem 3.7.

Let $A(x,\partial_{x})$ be a differential operator with rational coefficients
and let $F$ be the Picard-Vessiot extension of $\mathbb{C}(x)$ for $A$.
Consider a factorization
$A(x,\partial_{x})=A_{1}(x,\partial_{x})A_{2}(x,\partial_{x})$ with $A_{2}$
monic. Then $A_{1}(x,\partial_{x})$ and $A_{2}(x,\partial_{x})$ have rational
coefficients if and only if $\ker(A_{2})\subseteq\ker(A)$ is invariant under
the action of $\operatorname{\mathrm{DGal}}(F/\mathbb{C}(x))$.

###### Proof.

For $\sigma\in\operatorname{\mathrm{DGal}}(F/\mathbb{C}(x))$, let
$\sigma(A_{j}):=\sigma(A_{j})(x,\partial_{x})$ denote the operator obtained by
applying the automorphism to the coefficients. Since the automorphism
preserves differentiation, we know that

$\sigma(A_{j})(x,\partial_{x})\cdot\sigma(a)=\sigma(A_{j}(x,\partial_{x})\cdot
a),\ \forall a\in F.$

If $A_{1}(x,\partial_{x})$ and $A_{2}(x,\partial_{x})$ have rational
coefficients, then clearly $\sigma(A_{j})=A_{j}$ and therefore
$\ker(\sigma(A_{j}))=\ker(A_{j})$. Thus the kernel of $A_{j}(x,\partial_{x})$
is invariant under the action of
$\operatorname{\mathrm{DGal}}(F/\mathbb{C}(x))$.

Conversely, suppose that $\ker(A_{2})\subseteq\ker(A)$ is invariant under the
action of the differential Galois group, ie. $\sigma(\ker(A_{2}))=\ker(A_{2})$
Then $\sigma(A_{2})\cdot\sigma(a)=\sigma(A_{2}\cdot a)=\sigma(0)=0$ for all
$a\in\ker(A)$ and therefore $\ker(A_{2})\subseteq\ker(\sigma(A_{2}))$. Since
the order of $A_{2}$ and the order of $\sigma(A_{2})$ are the same, their
kernels will have the same dimension. Therefore
$\ker(\sigma(A_{2}))=\ker(A_{2})$ and consequently $\sigma(A_{2})=bA_{2}$ for
some $b\in F$. Since $A_{2}$ has leading coefficient $1$, it follows that
$b=1$. Hence $\sigma(A_{2})=A_{2}$ and from the Fundamental Theorem of
Differential Galois Theory, the coefficients of $A_{2}$ must all be rational
functions. Lastly, since $A$ and $A_{2}$ have rational coefficients, it
follows that $A_{1}$ has rational coefficients. ∎

######  Corollary 3.8.

Let $A(x,\partial_{x})$ be a self-adjoint differential operator with rational
coefficients and let $F$ be the Picard-Vessiot extension for $A$. Then the
self-adjoint, rational factorization of $A(x,\partial_{x})$ correspond
precisely with the $\operatorname{\mathrm{DGal}}(F/\mathbb{C}(x))$-invariant
Lagrangian subspaces of $\ker(A)$.

###### Proof.

This follows immediately from the theorem and the results of the previous
subsection. ∎

###  3.3. The classification

Now let $a_{1},\dots,a_{r}\in\mathbb{C}$ be the distinct roots of $q(z)$ and
write

$q(z)=(z-a_{1})^{d_{1}}\dots(z-a_{r})^{d_{r}}$

for some positive integers $d_{1},\dots,d_{r}$ and distinct
$a_{1},\ldots,a_{r}\in\mathbb{C}$. The kernel of $q(L)^{2}$ for
$L(x,\partial_{x})=\partial_{x}^{2}-x$ the Airy operator is given by the
following lemma.

######  Lemma 3.9.

The kernel of $q(L)^{2}$ has basis given by

$\\{\operatorname{{\mathrm{Ai}}}^{(j)}(x+a_{i}),\operatorname{\mathrm{Bi}}^{(j)}(x+a_{i}):1\leq
k\leq r,\ 0\leq j\leq 2d_{k}\\}$

where here $\operatorname{{\mathrm{Ai}}}(x)$ and
$\operatorname{\mathrm{Bi}}(x)$ are the Airy functions of the first and second
kind, respectively.

###### Proof.

To prove this, we will rely on the fundamental relation

$L(x,\partial_{x})\partial_{x}=\partial_{x}L(x,\partial_{x})+1,$

which implies that

$(L(x,\partial_{x})-a_{k})^{m}\partial_{x}^{n}=\sum_{j=0}^{m\wedge
n}\binom{m}{j}\frac{n!}{(n-j)!}\partial_{x}^{n-j}(L(x,\partial_{x})-a_{k})^{m-j}.$

Thus for all $0\leq n<2d_{k}-1$, we have

$(L(x,\partial_{x})-a_{k})^{2d_{k}}\operatorname{{\mathrm{Ai}}}^{(n)}(x+a_{k})=\sum_{j=0}^{n}\binom{2d_{k}}{j}\frac{n!}{(n-j)!}\partial_{x}^{n-j}(L(x,\partial_{x})-a_{k})^{2d_{k}-j}\operatorname{{\mathrm{Ai}}}(x+a_{k})=0.$

Hence
$\operatorname{{\mathrm{Ai}}}^{(n)}(x+a_{k})\in\ker((L(x,\partial_{x})-a_{k})^{2d_{k}})\subseteq\ker(q(L)^{2})$
for all $0\leq n<2d_{k}$. The same calculation shows that
$\operatorname{\mathrm{Bi}}^{(n)}(x+a_{k})\in\ker(q(L)^{2})$ for all $0\leq
n<2d_{k}$. ∎

Thus the Picard-Vessiot extension of the differential field
$(\mathbb{C}(x),\partial_{x})$ corresponding to the linear differential
operator $q(L(x,\partial_{x}))^{2}$ is finitely generated by $4r$ elements

$F^{q}=\mathbb{C}(x)(\operatorname{{\mathrm{Ai}}}(x+a_{k}),\operatorname{{\mathrm{Ai}}}^{\prime}(x+a_{k}),\operatorname{\mathrm{Bi}}(x+a_{k}),\operatorname{\mathrm{Bi}}^{\prime}(x+a_{k}):1\leq
k\leq r).$

Using this, we see that the differential Galois group of $F$ is isomorphic to
$r$ copies of $\text{SL}_{2}(\mathbb{C})$.

######  Lemma 3.10.

The differential Galois group consists of all differential
$\mathbb{C}(x)$-linear morphisms

$\sigma:F^{q}\rightarrow
F^{q},\quad\left\\{\begin{array}[]{cc}\operatorname{{\mathrm{Ai}}}(x+a_{k})\mapsto\alpha_{k}\operatorname{{\mathrm{Ai}}}(x+a_{k})+\beta_{k}\operatorname{\mathrm{Bi}}(x+a_{k})\\\
\operatorname{\mathrm{Bi}}(x+a_{k})\mapsto\gamma_{k}\operatorname{{\mathrm{Ai}}}(x+a_{k})+\delta_{k}\operatorname{\mathrm{Bi}}(x+a_{k})\end{array}\right.\
\forall 1\leq k\leq r,$

where here $\alpha_{k},\beta_{k},\gamma_{k},\delta_{k}\in\mathbb{C}$ with
$\alpha_{k}\delta_{k}-\beta_{k}\gamma_{k}=1$.

###### Proof.

The fact that

$\begin{array}[]{c}\operatorname{{\mathrm{Ai}}}(x+a_{k})\mapsto\alpha\operatorname{{\mathrm{Ai}}}(x+a_{k})+\beta\operatorname{\mathrm{Bi}}(x+a_{k})\\\
\operatorname{\mathrm{Bi}}(x+a_{k})\mapsto\gamma\operatorname{{\mathrm{Ai}}}(x+a_{k})+\delta\operatorname{\mathrm{Bi}}(x+a_{k})\end{array},\quad\left[\begin{array}[]{cc}\alpha&\beta\\\
\gamma&\delta\end{array}\right]\in\text{SL}_{2}(\mathbb{C}),$

is a differential automorphism is standard. See for example [31, Example
8.15]. Therefore, we need only show that this accounts for all differential
automorphisms.

If $\sigma:F^{q}\rightarrow F^{q}$ is a differential automorphism fixing
$\mathbb{C}(x)$, then

$\frac{\sigma(\operatorname{{\mathrm{Ai}}}(x+a_{k}))^{\prime\prime}}{\sigma(\operatorname{{\mathrm{Ai}}}(x+a_{k}))}=\sigma\left(\frac{\operatorname{{\mathrm{Ai}}}^{\prime\prime}(x+a_{k})}{\operatorname{{\mathrm{Ai}}}(x+a_{k})}\right)=\sigma(x+a_{k})=x+a_{k}.$

Thus $\sigma(\operatorname{{\mathrm{Ai}}}(x+a_{k}))$ must be a solution of the
differential equation $y^{\prime\prime}=(x+a_{k})y$, and therefore a linear
combination of $\operatorname{{\mathrm{Ai}}}(x+a_{k})$ and
$\operatorname{\mathrm{Bi}}(x+a_{k})$ for all $k$. A similar statement holds
for $\sigma(\operatorname{\mathrm{Bi}}(x+a_{k}))$ so that

$\sigma:\left\\{\begin{array}[]{cc}\operatorname{{\mathrm{Ai}}}(x+a_{k})\mapsto\alpha_{k}\operatorname{{\mathrm{Ai}}}(x+a_{k})+\beta_{k}\operatorname{\mathrm{Bi}}(x+a_{k})\\\
\operatorname{\mathrm{Bi}}(x+a_{k})\mapsto\gamma_{k}\operatorname{{\mathrm{Ai}}}(x+a_{k})+\delta_{k}\operatorname{\mathrm{Bi}}(x+a_{k})\end{array}\right.\
\forall 1\leq k\leq r$

for some $\alpha_{k},\beta_{k},\gamma_{k},\delta_{k}\in\mathbb{C}$. Lastly,
the Wronskian identity implies

$W(\operatorname{{\mathrm{Ai}}}(x+a_{k}),\operatorname{\mathrm{Bi}}(x+a_{k}))=\operatorname{{\mathrm{Ai}}}^{\prime}(x+a_{k})\operatorname{\mathrm{Bi}}(x+a_{k})-\operatorname{{\mathrm{Ai}}}(x+a_{k})\operatorname{\mathrm{Bi}}^{\prime}(x+a_{k})=\frac{1}{\pi}\cdot$

Since the Wronskian is skew-symmetric, we can conclude that

$\displaystyle\frac{1}{\pi}$
$\displaystyle=\sigma(W(\operatorname{{\mathrm{Ai}}}(x+a_{k}),\operatorname{\mathrm{Bi}}(x+a_{k})))$
$\displaystyle=W(\sigma(\operatorname{{\mathrm{Ai}}}(x+a_{k})),\sigma(\operatorname{\mathrm{Bi}}(x+a_{k})))$
$\displaystyle=(\alpha\delta-\beta\gamma)W(\operatorname{{\mathrm{Ai}}}(x+a_{k}),\operatorname{\mathrm{Bi}}(x+a_{k}))=(\alpha\delta-\beta\gamma)/\pi.$

Hence $\alpha\delta-\beta\gamma=1$. ∎

Using this, we can obtain the following characterization of the Galois-
invariant subspaces of $\ker(q(L)^{2})$.

######  Lemma 3.11.

Suppose that $V\subseteq\ker(q(L)^{2})$ is a subspace. Then $V$ is invariant
under the action of the differential Galois group if and only if $V$ is
spanned by _pairs_ of elements of the form

$\sum_{j=0}^{2d_{k}-1}\alpha_{kj}\operatorname{{\mathrm{Ai}}}^{(j)}(x+a_{k}),\quad\sum_{j=0}^{2d_{k}-1}\alpha_{kj}\operatorname{\mathrm{Bi}}^{(j)}(x+a_{k})$

###### Proof.

Clearly any subspace spanned by pairs of elements of this form is invariant
under the action of the Galois group, since the the action restricts to an
action sending each of the functions in the pair to a linear combination of
the functions in the pair. Thus it suffices to show the converse.

Let $f(x)$ be a nonzero element of $V$. Then

$f(x)=\sum_{k=1}^{r}\sum_{j=0}^{2d_{k}-1}\alpha_{kj}\operatorname{{\mathrm{Ai}}}^{(j)}(x+a_{k})+\beta_{kj}\operatorname{\mathrm{Bi}}^{(j)}(x+a_{k})$

Consider the differential automorphisms $\sigma_{k}$ and $\tau_{k}$ which fix
$\operatorname{{\mathrm{Ai}}}(x+a_{j})$ and
$\operatorname{\mathrm{Bi}}(x+a_{j})$ and satisfy

$\sigma_{k}:\operatorname{{\mathrm{Ai}}}(x+a_{k})\mapsto-\operatorname{\mathrm{Bi}}(x+a_{k}),\
\operatorname{\mathrm{Bi}}(x+a_{k})\mapsto\operatorname{{\mathrm{Ai}}}(x+a_{k}),$
$\tau_{k}:\operatorname{{\mathrm{Ai}}}(x+a_{k})\mapsto\operatorname{{\mathrm{Ai}}}(x+a_{k})+\operatorname{\mathrm{Bi}}(x+a_{k}),\
\operatorname{\mathrm{Bi}}(x+a_{k})\mapsto\operatorname{\mathrm{Bi}}(x+a_{k}).$

We see that

$\tau_{k}(f(x))-f(x)=\sum_{j=0}^{2d_{k}-1}\alpha_{kj}\operatorname{\mathrm{Bi}}^{(j)}(x+a_{k})\in
V.$

Following up by applying $\sigma_{k}$, we see that

$\sigma_{k}(\tau_{k}(f(x))-f(x))=\sum_{j=0}^{2d_{k}-1}\alpha_{kj}\operatorname{{\mathrm{Ai}}}^{(j)}(x+a_{k})\in
V.$

Likewise, one may show
$\sum_{j=0}^{2d_{k}-1}\beta_{kj}\operatorname{{\mathrm{Ai}}}^{(j)}(x+a_{k}),\
\sum_{j=0}^{2d_{k}-1}\beta_{kj}\operatorname{\mathrm{Bi}}^{(j)}(x+a_{k})\in V$
and since $k$ was arbitrary, the statement of the Lemma follows immediately. ∎

Our explicit description of the kernel of $q(L)^{2}$ allows us to give a
concrete formula for the symplectic form on $\ker(q(L)^{2})$ defined by the
bilinear concomitant. We start with a combinatorial Lemma.

######  Lemma 3.12.

Let $a,b,m$ be integers. Then

$\sum_{k=0}^{m}(-1)^{k}\binom{k+a}{k}\binom{b}{m-k}=\binom{b-1-a}{m}.$

###### Proof.

We use the binomial series expansion on the identity

$(1-z)^{-a-1}(1-z)^{b}=(1-z)^{b-a-1}$

to find

$\sum_{j,k=0}^{\infty}\binom{-a-1}{k}\binom{b}{j}(-1)^{j+k}z^{j+k}=\sum_{m=0}^{\infty}\binom{b-1-a}{m}(-1)^{m}z^{m}.$

Now comparing coefficients of $z^{m}$:

$\sum_{k=0}^{m}\binom{-a-1}{k}\binom{b}{m-k}(-1)^{j+k}=\binom{b-1-a}{m}.$

Noting that $\binom{-a-1}{k}=(-1)^{k}\binom{k+a}{a}$, the statement of the
lemma follows immediately. ∎

######  Proposition 3.13.

Let
$f(x),g(x)\in\\{\operatorname{{\mathrm{Ai}}}(x),\operatorname{\mathrm{Bi}}(x)\\}$
and choose $0\leq m<2d_{j}$ and $0\leq n<2d_{k}$. Then

$\mathcal{C}_{q(L)^{2}}(f^{(m)}(x+a_{j}),g^{(n)}(x+a_{k}))=\delta_{jk}\frac{m!n!W(f,g)}{(m+n-2d_{k}+1)!}\left.\partial_{z}^{m+n-2d_{k}+1}\right\rvert_{z=a_{k}}\cdot\left(\frac{q(z)^{2}}{(z-a_{k})^{2d_{k}}}\right)$

for all nonnegative integers $m,n$ with $m+n\geq 2d_{k}-1$ and is zero
otherwise.

###### Proof.

For simplicity of notation, we will let $h(z)=q(z)^{2}$ and write $f$ and $g$
in place of $f(x+a_{k})$ and $g(x+b_{k})$, respectively. First note that if
$j\neq k$ then $f^{(m)}\in\ker((L-a_{j})^{2d_{j}})$ and
$g^{(n)}\in\ker(h(L)(L-a_{j})^{-2d_{j}})$, which is the orthogonal complement
of the subspace $\ker((L-a_{j})^{2d_{j}})$ of $\ker(h(L))$. Hence
$\mathcal{C}_{h(L)}(f^{(m)},g^{(n)})=0$. Thus it suffices to consider the case
when $j=k$.

Let $\widetilde{h}(z)=h(z)/(z-a_{k})^{2d_{k}}$. Applying Lemma 3.2 and the
fundamental relation $L\partial_{x}=\partial_{x}L+1$ we see that

$\displaystyle\mathcal{C}_{h(L)}(f^{(m)},g^{(n)})$
$\displaystyle=\mathcal{C}_{h(L)\partial_{x}^{m}}(f,g^{(n)})$
$\displaystyle=\mathcal{C}_{\widetilde{h}(L)(L-a_{k})^{2d_{k}}\partial_{x}^{m}}(f,g^{(n)})$
$\displaystyle=\sum_{s=0}^{m}\binom{m}{s}\frac{(2d_{k})!}{(2d_{k}-s)!}\mathcal{C}_{\widetilde{h}(L)\partial_{x}^{m-s}(L-a_{k})^{2d_{k}-s}}(f,g^{(n)}).$

Now using the fact that the concomitant of $L$ is the Wronskian and again
applying Lemma 3.2 and the more general relation

$\widetilde{h}(L(x,\partial_{x}))\partial_{x}^{m}=\sum_{s=0}^{m}\binom{m}{s}\partial_{x}^{m-s}\widetilde{h}^{(s)}(L(x,\partial_{x}))$

we see that

$\displaystyle\mathcal{C}_{h(L)}$ $\displaystyle(f^{(m)},g^{(n)})$
$\displaystyle=\sum_{s=0}^{m}\binom{m}{s}\frac{(2d_{k})!}{(2d_{k}-s)!}(-1)^{m-s}W(f,(L-a_{k})^{2d_{k}-s-1}\partial_{x}^{m-s}\widetilde{h}(L)\partial_{x}^{n}\cdot
g)$
$\displaystyle=\sum_{t=0}^{n}\binom{n}{t}\widetilde{h}^{(t)}(a_{k})\sum_{s=0}^{m}\binom{m}{s}\frac{(2d_{k})!}{(2d_{k}-s)!}(-1)^{m-s}W(f,(L-a_{k})^{2d_{k}-s-1}\partial_{x}^{n+m-s-t}\cdot
g).$

From this it is clear that if $n+m<2d_{k}-1$ then the concomitant is zero.
Thus without loss of generality we take $m+n\geq 2d_{k}-1$. Then for
$\ell=n+m-2d_{k}+1$

$\displaystyle\mathcal{C}_{h(L)}$ $\displaystyle(f^{(m)},g^{(n)})$
$\displaystyle=\sum_{t=0}^{n}\binom{n}{t}\widetilde{h}^{(t)}(a_{k})\sum_{s=0}^{m}\binom{m}{s}\frac{(2d_{k})!}{(2d_{k}-s)!}(-1)^{m-s}W(f,(L-a_{k})^{2d_{k}-s-1}\partial_{x}^{n+m-s-t}\cdot
g)$
$\displaystyle=\sum_{t=0}^{n\wedge\ell}\binom{n}{t}\widetilde{h}^{(t)}(a_{k})\sum_{s=0}^{m}\binom{m}{s}\frac{(2d_{k})!}{(2d_{k}-s)!}(-1)^{m-s}\frac{(m+n-s-t)!}{(\ell-t)!}W(f,\partial_{x}^{\ell-t}\cdot
g)$
$\displaystyle=\sum_{t=0}^{n\wedge\ell}\binom{n}{t}\widetilde{h}^{(t)}(a_{k})\frac{m!(n-t)!}{(\ell-t)!}\sum_{s=0}^{m}\binom{2d_{k}}{s}\binom{m+n-s-t}{m-s}(-1)^{m-s}W(f,\partial_{x}^{\ell-t}\cdot
g).$

Now reindexing the sum and applying the previous lemma, we obtain

$\displaystyle\mathcal{C}_{h(L)}$ $\displaystyle(f^{(m)},g^{(n)})$
$\displaystyle=\sum_{t=0}^{n\wedge\ell}\binom{n}{t}\widetilde{h}^{(t)}(a_{k})\frac{m!(n-t)!}{(\ell-t)!}\sum_{s=0}^{m}\binom{2d_{k}}{m-s}\binom{s+n-t}{s}(-1)^{s}W(f,\partial_{x}^{\ell-t}\cdot
g)$
$\displaystyle=\sum_{t=0}^{n\wedge\ell}\binom{n}{t}\widetilde{h}^{(t)}(a_{k})\frac{m!(n-t)!}{(\ell-t)!}\binom{2d_{k}-1-n+t}{m}W(f,\partial_{x}^{\ell-t}\cdot
g).$

The binomial coefficient in the last sum is nonzero if and only if $\ell\leq
t$. Since the sum is taken between $t=0$ and $t=\ell$, the only nonzero term
comes from when $t=\ell$. Thus

$\mathcal{C}_{h(L)}(f^{(m)},g^{(n)})=\frac{m!n!}{(n+m-2d_{k}+1)!}\widetilde{h}^{(n+m-2d_{k}+1)}(a_{k})W(f,g).$

∎

The rational Darboux transformations of $\operatorname{{\mathrm{Ai}}}(x+z)$
come directly from factorizations of the form (7) with $P(x,\partial_{x})$
having rational coefficients. As we have outlined above, these correspond
precisely to the Galois-invariant Lagrangian subspaces of $\ker(q(L)^{2})$.
This characterization is made explicit in the next theorem.

######  Theorem 3.14 (Classification Theorem).

Let $f_{m},g_{m}\in\ker(q(L)^{2})$ for $1\leq m\leq d$ be $2d$ linearly
independent functions of the form

$f_{i}(x)=\sum_{m=0}^{2d_{\ell_{i}}-1}\alpha_{im}\operatorname{{\mathrm{Ai}}}^{(m)}(x+a_{\ell_{i}}),\quad
g_{i}(x)=\sum_{n=0}^{2d_{\ell_{i}}-1}\alpha_{in}\operatorname{\mathrm{Bi}}^{(n)}(x+a_{\ell_{i}})$

satisfying the condition that

$\sum_{m+n\geq
2d_{k}-1}^{2d_{k}-1}\alpha_{im}\alpha_{jn}\frac{m!n!}{(m+n-2d_{k}+1)!}\left.\partial_{z}^{m+n-2d_{k}+1}\right\rvert_{z=a_{k}}\cdot\left(\frac{q(z)^{2}}{(z-a_{k})^{2d_{k}}}\right)=0$

for all $k$ and for all $i,j$ with $\ell_{i}=\ell_{j}=k$. Then the
differential operator $P(x,\partial_{x})$ of order $2d$ defined in terms of a
Wronskian by

$P(x,\partial_{x})\cdot
f:=W(f_{1},f_{2},\dots,f_{d},g_{1},g_{2},\dots,g_{d},f)$

has rational coefficients and satisfies

$P(x,\partial_{x})^{*}\frac{1}{p(x)^{2}}P(x,\partial_{x})=q(L(x,\partial_{x}))^{2}$

for some rational function $p(x)$. Furthermore every self-adjoint rational
factorization of $q(L(x,\partial_{x}))^{2}$ is of this form.

###### Proof.

This follows directly from our direct calculation of the concomitant along
with our characterization of the Galois-invariant subspaces of the kernel. ∎

This result is particularly nice in the situation that $q(z)=(z-s_{1})^{d}$,
so that the concomitant has the simple form

$\mathcal{C}_{q(L)^{2}}(f^{(m)}(x+s_{1}),g^{(n)}(x+s_{1}))=\left\\{\begin{array}[]{cc}m!n!/\pi,&f=\operatorname{{\mathrm{Ai}}},\
g=\operatorname{\mathrm{Bi}},\ m+n=2d-1\\\
-m!n!/\pi,&f=\operatorname{\mathrm{Bi}},\ g=\operatorname{{\mathrm{Ai}}},\
m+n=2d-1\\\ 0,&\text{otherwise}.\end{array}\right.$

The payout of our dive through all the differential Galois theory and
symplectic geometry above is that we immediately provide _explicit
factorizations_ of $(L(x,\partial_{x})-s_{1})^{2}$ and
$(L(x,\partial_{x})-s_{1})^{4}$.

######  Corollary 3.15.

Let $s_{1}\in\mathbb{C}$. Then up to a function multiple, the only self-
adjoint rational factorizations of $(L(x,\partial_{x})-s_{1})^{2}$ are the
trivial one and

$P_{1}(x,\partial_{x})^{*}\frac{1}{(x+s_{1})^{2}}P_{1}(x,\partial_{x})=(L(x,\partial_{x})-s_{1})^{2}$

for

$P_{1}(x,\partial_{x})=(x+s_{1})\partial_{x}^{2}-\partial_{x}-(x+s_{1})^{2}.$

###### Proof.

From the previous Theorem, we know we must choose functions

$f_{1}(x)=\alpha_{11}\operatorname{{\mathrm{Ai}}}(x+s_{1})+\alpha_{12}\operatorname{{\mathrm{Ai}}}^{\prime}(x+s_{1}),\quad
g_{1}(x)=\alpha_{11}\operatorname{\mathrm{Bi}}(x+s_{1})+\alpha_{12}\operatorname{\mathrm{Bi}}^{\prime}(x+s_{1})$

satisfying $2\alpha_{11}\alpha_{12}1!1!/\pi=0$. Thus either $\alpha_{11}=0$ or
$\alpha_{12}=0$ and without loss of generality we may take the remaining
coefficient to be $\pi$. In the first case, the operator $P(x,\partial_{x})$
is

$P(x,\partial_{x})\cdot
f=W(\operatorname{{\mathrm{Ai}}}^{\prime}(x+s_{1}),\operatorname{\mathrm{Bi}}^{\prime}(x+s_{1}),f)=(x+s_{1})f^{\prime\prime}(x)-f^{\prime}(x)-(x+s_{1})^{2}f.$

In the second case, the operator $P(x,\partial_{x})$ is

$P(x,\partial_{x})\cdot
f=W(\operatorname{{\mathrm{Ai}}}(x+s_{1}),\operatorname{\mathrm{Bi}}(x+s_{1}),f)=f^{\prime\prime}(x)-(x+s_{1})f.$

Thus in this second case $P(x,\partial_{x})=L(x,\partial_{x})-s_{1}$, giving
us the trivial factorization of $(L(x,\partial_{x})-s_{1})^{2}$. ∎

######  Corollary 3.16.

Let $s_{1}\in\mathbb{C}$. Then up to a function multiple, the self-adjoint
rational factorizations of $(L(x,\partial_{x})-s_{1})^{4}$ are of the form

$P_{2}(x,\partial_{x})^{*}\frac{1}{(x+s_{1})^{2}}P_{2}(x,\partial_{x})=(L(x,\partial_{x})-s_{1})^{4}$

for

$P_{2}(x,\partial_{x})\cdot f=W(f_{1},f_{2},g_{1},g_{2},f),$

where here

$f_{k}(x)=\sum_{j=0}^{3}\alpha_{kj}\operatorname{{\mathrm{Ai}}}^{(j)}(x+s_{1}),\quad
g_{k}(x)=\sum_{j=0}^{3}\alpha_{kj}\operatorname{\mathrm{Bi}}^{(j)}(x+s_{1})$

for some constants $\alpha_{kj}$ satisfying the three relations

$6\alpha_{m3}\alpha_{n0}+2\alpha_{m1}\alpha_{n2}+2\alpha_{m2}\alpha_{n1}+6\alpha_{m3}\alpha_{n0}=0,\
\ 1\leq m\leq n\leq 2.$

###### Proof.

This follows immediately from the Classification Theorem. ∎

The operator $P_{2}(x,\partial_{x})$ in this latter situation is more
complicated. First of all, it features the factorizations from the previous
corollary, as may be obtained from taking $\alpha_{13}=\alpha_{23}=0$. Thus to
get new factorizations, we can without loss of generality take
$\alpha_{13}=\alpha_{23}=1$. Then the three relations simplify to
$\alpha_{m0}=-\alpha_{m1}\alpha_{m2}/3$ for $m=1,2$ plus a choice of either
$\alpha_{11}=\alpha_{21}$ or $\alpha_{12}=\alpha_{22}$. For sake of
concreteness, we choose $\alpha_{11}=\alpha_{21}$ and take $\alpha_{22}=1$,
$\alpha_{12}=0$, and $s_{1}=0$. This determines all parameters, except for
$\alpha_{11}$ and the associated operator $P(x,\partial_{x})$ is explicitly
computed to be

$\displaystyle P_{2}(x,\partial_{x})$
$\displaystyle=\left(x^{4}-4x^{3}\alpha_{11}+\frac{10}{3}x^{2}\alpha_{11}^{2}+\left(\frac{4}{3}\alpha_{11}^{3}+4\right)x+\frac{1}{9}\alpha_{11}^{4}-8\alpha_{11}\right)\partial_{x}^{4}$
$\displaystyle+\left(-4x^{3}+12x^{2}\alpha_{11}-\frac{20}{3}\alpha_{11}^{2}x-\frac{4}{3}\alpha_{11}^{3}-4\right)\partial_{x}^{3}$
(9)
$\displaystyle+\left(-2x^{5}+8x^{4}\alpha_{11}-\frac{20}{3}\alpha_{11}^{2}x^{3}-\left(\frac{8}{3}\alpha_{11}^{3}+2\right)x^{2}-\left(\frac{2}{9}\alpha_{11}^{4}-4\alpha_{11}\right)x+\frac{10}{3}\alpha_{11}^{2}\right)\partial_{x}^{2}$
$\displaystyle+\left(2x^{4}-4x^{3}\alpha_{11}-\frac{4}{3}x\alpha_{11}^{3}-16-\frac{2}{9}\alpha_{11}^{4}+36\alpha_{11}\right)\partial_{x}$
$\displaystyle+x^{6}-4x^{5}\alpha_{11}+\frac{10}{3}x^{4}\alpha_{11}^{2}+\left(\frac{4}{3}\alpha_{11}^{3}+8\right)x^{3}+\left(\frac{1}{9}\alpha_{11}^{4}-22\alpha_{11}\right)x^{2}+\frac{16}{3}x\alpha_{11}^{2}+2\alpha_{11}^{3}+16$

##  4\. Commuting differential operators for the level one kernels

In this section, we explore the commuting differential operators for integral
operators with level one Airy kernels, ie. those defined by bispectral
functions $\Psi$ obtained from self-adjoint rational Darboux transformations
of $(L(x,\partial_{x})-s_{1})^{2}$. There is only one such bispectral
function, determined by the factorization of $(L(x,\partial_{x})-s_{1})^{2}$
in Corollary 3.15. Using the operator $P_{1}(x,\partial_{x})$ described in
this Corollary, the associated bispectral function is

$\Psi_{1}(x,z)=\frac{1}{(x+s_{1})(z-s_{1})}P_{1}(x,\partial_{x})\cdot\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)=\operatorname{{\mathrm{Ai}}}(x+z)-\frac{1}{(x+s_{1})(z-s_{1})}\operatorname{{\mathrm{Ai}}}^{\prime}(x+z).$

Let
$\widetilde{P}_{1}(z,\partial_{z})=b_{\Psi_{\operatorname{{\mathrm{Ai}}}}}(P_{1}(x,\partial_{x}))$,
$p_{1}(x)=x+s_{1}$ and $q_{1}(z)=z-s_{1}$. For every
$R(x,\partial_{x})\in\mathfrak{F}_{x}(\Psi_{\operatorname{{\mathrm{Ai}}}})$
and
$S(z,\partial_{z})=b_{\Psi_{\operatorname{{\mathrm{Ai}}}}}(R(x,\partial_{x}))$,
we have the identities

(10)
$\frac{1}{p_{1}(x)}P_{1}(x,\partial_{x})^{*}R(x,\partial_{x})P_{1}(x,\partial_{x})\frac{1}{p_{1}(x)}\cdot\Psi_{1}(x,z)=(q_{1}(z))S(z,\partial_{z})(q_{1}(z))\cdot\Psi_{1}(x,z),$
(11)
$\frac{1}{q_{1}(z)}\widetilde{P}_{1}(z,\partial_{z})^{*}S(z,\partial_{z})\widetilde{P}_{1}(z,\partial_{z})\frac{1}{q_{1}(z)}\cdot\Psi_{1}(x,z)=(p_{1}(x))R(x,\partial_{x})(p_{1}(x))\cdot\Psi_{1}(x,z),$

and the more complicated identity

(12)
$\displaystyle\left(\frac{1}{p_{1}(x)}P_{1}(x,\partial_{x})R(x,\partial_{x})p(x)+p(x)R(x,\partial_{x})^{*}P_{1}(x,\partial_{x})^{*}\frac{1}{p_{1}(x)}\right)\cdot\Psi_{1}(x,z)$
$\displaystyle=\left(\frac{1}{q_{1}(z)}\widetilde{P}_{1}(z,\partial_{z})S(z,\partial_{z})q(z)+q(z)S(z,\partial_{z})^{*}\widetilde{P}_{1}(z,\partial_{z})^{*}\frac{1}{q_{1}(z)}\right)\cdot\Psi_{1}(x,z)$

Comparing the orders of these operators, we see that
$\mathfrak{F}_{x,\operatorname{\mathrm{sym}}}^{2\ell,2m}(\Psi_{1})$ contains
the direct sum

$\displaystyle\mathfrak{F}_{x,\operatorname{\mathrm{sym}}}^{2\ell,2m}(\Psi_{1})$
$\displaystyle\supseteq\frac{1}{p_{1}(x)}P_{1}(x,\partial_{x})^{*}\mathfrak{F}_{x,\operatorname{\mathrm{sym}}}^{2\ell-4,2m}(\Psi_{\operatorname{{\mathrm{Ai}}}})P_{1}(x,\partial_{x})\frac{1}{p_{1}(x)}$
$\displaystyle\oplus
p_{1}(x)\mathfrak{F}_{x,\operatorname{\mathrm{sym}}}^{2,2m-4}(\Psi_{\operatorname{{\mathrm{Ai}}}})p_{1}(x)\oplus\mathfrak{E}\oplus\mathbb{C}$

for all $\ell,m\geq 2$, where here $\mathfrak{E}$ is a set of additional
operators stemming from the equation (12)

$\mathfrak{E}=\left\\{\frac{1}{p_{1}(x)}P_{1}(x,\partial_{x})R(x,\partial_{x})p(x)+p(x)R(x,\partial_{x})^{*}P_{1}(x,\partial_{x})^{*}\frac{1}{p_{1}(x)}:R(x,\partial_{x})\in\mathfrak{F}^{1,1}_{x}(\Psi_{\operatorname{{\mathrm{Ai}}}})\right\\}.$

Explicit calculation shows that $\mathfrak{E}$ is two dimensional.
Consequently the dimension of
$\mathfrak{F}_{x,\operatorname{\mathrm{sym}}}^{2\ell,2m}(\Psi_{1})$ is at
least $(\ell-1)(m+1)+2(m-1)+2+1=(\ell+1)(m+1)-1$. One can show that this is
precisely the dimension for all $m,n>1$ and that both $\mathfrak{F}_{x}(\Psi)$
and $\mathfrak{F}_{z}(\psi)$ are equal to algebras of differential operators
on a rational curve with a cusipdal singularity of degree $2$ at the origin.

Let $T_{1}$ be the integral operator

$T_{1}:f(z)\mapsto\int_{t_{1}}^{\infty}K_{1}(z,w)f(w)dw,\quad
K_{1}(z,w)=\int_{t_{2}}^{\infty}\Psi_{1}(x,z)\Psi_{1}(x,w)dx.$

The specific value of the kernel $K_{1}(z,w)$ is determined via integration by
parts to be

$K_{1}(z,w)=\frac{q_{1}(w)}{q_{1}(z)}K_{\operatorname{{\mathrm{Ai}}}}(z,w)+\mathcal{C}_{P_{1}}(\psi_{\operatorname{{\mathrm{Ai}}}}(x,z),\psi_{1}(x,w)/p_{1}(x);t_{2}).$

From the previous estimate of the dimension of
$\mathfrak{F}_{x,\operatorname{\mathrm{sym}}}^{2\ell,2m}(\Psi_{1})$, we see
that $T_{1}$ will commute with a differential operator $S_{1}(z,\partial_{z})$
in $\mathfrak{F}_{z,\operatorname{\mathrm{sym}}}^{4,4}(\Psi_{1})$.

The values of the commuting integral and differential operators will in
general depend on $s_{1}$, albeit predictably. If we make the
$s_{1}$-dependence of $\Psi(x,z)=\Psi(x,z;s_{1})$ explicit, we see
$\Psi(x,z;s_{1})=\Psi(x+s_{1},z-s_{1};0)$ and consequently the differential
operator $S_{1}(z,\partial_{z})$ commuting with $T_{1}$ for arbitrary $s_{1}$
is the same as in the case $s_{1}=0$, but with $z$ replaced by $z-s_{1}$ and
$t_{2}$ replaced by $t_{2}+s_{1}$. Thus without loss of generality we will
take $s_{1}=0$.

Explicitly computing the condition of the vanishing of the concomitant and
solving the resulting linear system of equations yields the operator of order
$4$

$S_{1}(z,\partial_{z})=\frac{1}{z}\left(\sum_{k=0}^{2}\partial_{z}^{k}a_{k}(z)(z-t_{1})^{k}\partial_{z}^{k}\right)\frac{1}{z}$

where here

$\displaystyle a_{2}(z)$ $\displaystyle=z^{2},$ $\displaystyle a_{1}(z)$
$\displaystyle=-2(z^{4}+(t_{2}-t_{1})z^{3}-3t_{1}),$ $\displaystyle a_{0}(z)$
$\displaystyle=z^{3}(z^{3}+2(t_{2}-t_{1})z^{2}+(t_{2}-t_{1})^{2}z-8)+(t_{1}+t_{2})z^{2}/3.$

The dimension estimates also imply the existence of a commuting differential
operator of order $6$, which we find to be

$\widetilde{S}_{1}(z,\partial_{z})=\frac{1}{z}\left(\sum_{k=0}^{3}\partial_{z}^{k}\widetilde{a}_{k}(z)(z-t_{1})^{k}\partial_{z}^{k}\right)\frac{1}{z}$

where here

$\displaystyle\widetilde{a}_{3}(z)=$ $\displaystyle\ z^{2},$
$\displaystyle\widetilde{a}_{2}(z)=$ $\displaystyle\
-3(z^{4}+(t_{2}-t_{1})z^{3}-4t_{1}),$ $\displaystyle\widetilde{a}_{1}(z)=$
$\displaystyle\
3(z^{6}+2(t_{2}-t_{1})z^{5}+(t_{2}-t_{1})^{2}z^{4}-10z^{3}+(5t_{1}-4t_{2})z^{2}-3t_{1}(t_{2}-t_{1})z),$
$\displaystyle\widetilde{a}_{0}(z)=$ $\displaystyle\
-z^{8}-3(t_{2}-t_{1})z^{7}-3(t_{2}-t_{1})^{2}z^{6}-((t_{2}-t_{1})^{3}-32)z^{5}$
$\displaystyle+(42t_{2}-63t_{1})z^{4}+(36t_{1}^{2}-48t_{1}t_{2}+12t_{2}^{2})z^{3}+t_{1}t_{2}(t_{1}+t_{2})z^{2}+12t_{1}^{2}-6t_{1}t_{2}.$

The operators $S_{1}(z,\partial_{z})$ and
${\widetilde{S}}_{1}(z,\partial_{z})$ commute and thus satisfy an algebraic
relation. The relation is

$\widetilde{S}_{1}^{2}=S_{1}^{3}-\frac{t_{1}^{2}-t_{1}t_{2}+t_{2}^{2}}{3}S_{1}+\frac{(t_{1}-2t_{2})(2t_{1}-t_{2})(t_{1}+t_{2})}{3^{3}}.$

The discriminant of the polynomial on the right hand side is

$\Delta=-\frac{16}{27}(260t_{1}^{6}-780t_{1}^{5}t_{2}-627t_{1}^{4}t_{2}^{2}+2554t_{1}^{3}t_{2}^{3}-627t_{1}^{2}t_{2}^{4}-780t_{1}t_{2}^{5}+260t_{2}^{6}),$

so for generic values of $t_{1}$ and $t_{2}$, the associated algebraic variety
is an elliptic curve.

##  5\. Commuting differential operators for the level two kernels

In this section, we explore the commuting differential operators for integral
operators with level two Airy kernels, ie. those defined by bispectral
functions $\Psi$ obtained from self-adjoint rational Darboux transformations
of $(L(x,\partial_{x})-s_{1})^{2}(L(x,\partial_{x})-s_{2})^{2}$. We will focus
on the particular case when $s_{1}=s_{2}$, leaving the other situation to a
future publication. Note also that due to the nice translation behavior of
$\psi_{\operatorname{{\mathrm{Ai}}}}(x,z)$, we can easily rederive the formula
for general values of $s_{1}$ from the case when $s_{1}=0$. So for sake of
simplicity, we will take $s_{1}=0$.

There are many bispectral functions in the level two case, all of which are
determined by the factorizations of $L(x,\partial_{x})^{4}$ in Corollary 3.16,
which in turn are determined by a choice of $\alpha_{jk}$ for $j=1,2$ and
$0\leq k\leq 3$ satisfying the constraints of the Corollary. The precise value
$P_{2}(x,\partial_{x})$ and the commuting operator is very complicated in
general. To facilitate our computations, and the inclusion of exact formulas
in our paper, we will take $\alpha_{31}=\alpha_{32}=1$,
$\alpha_{11}=\alpha_{21}$ and take $\alpha_{22}=1$, $\alpha_{12}=0$, so that
$P_{2}(x,\partial_{x})$ is given by (3.3). Additionally we will take
$\alpha_{11}=0$ so that $P_{2}(x,\partial_{x})$ has the simplified formula

$P_{2}(x,\partial_{x})=x(x^{3}+4)\partial_{x}^{4}-4(x^{3}+1)\partial_{x}^{3}-2x^{2}(x^{3}+1)\partial_{x}^{2}+2x(x^{3}-8)\partial_{x}+x^{6}+8x^{3}+16.$

Let $q_{2}(z)=z^{2}$ and $p_{2}(x)=x(x^{3}+4)$. The corresponding bispectral
function is defined by

$\displaystyle\Psi_{2}(x,z)$
$\displaystyle=\frac{1}{p_{2}(x)q_{2}(z)}P_{2}(x,\partial_{x})\cdot\Psi_{\operatorname{{\mathrm{Ai}}}}(x,z)$
$\displaystyle=\operatorname{{\mathrm{Ai}}}(x+z)+\frac{6(x^{3}+x^{2}z+2)}{p_{2}(x)q_{2}(z)}\operatorname{{\mathrm{Ai}}}(x+z)-\frac{4(x^{3}w+3x+w)}{p_{2}(x)q_{2}(z)}\operatorname{{\mathrm{Ai}}}^{\prime}(x+z).$

The Fourier algebras for $\Psi_{2}(x,z)$ are given by algebras of differential
operators on some rank $1$, torsion-free modules over certain rational curves
with cuspidal singularities. Specifically, let
$\mathscr{A}_{x}=\\{f(x)\in\mathbb{C}[x]:p(x)|f^{\prime}(x)\\}$ be the
coordinate ring of a singular rational curve $X$ with cusps of order $2$ at
the roots of $p(x)$. Then

$\mathfrak{F}_{x}(\Psi_{2})=\\{D(x,\partial_{x}):D(x,\partial_{x})\cdot\mathscr{A}_{x}\subseteq\mathscr{A}_{x}\\}$

is the algebra of differential operators on $X$. Likewise, let
$\mathscr{A}_{z}=\mathbb{C}[z^{4},z^{5}]$ be the affine coordinate ring of a
rational curve $Z$ with a higher-order cusp at $0$ and consider the torsion-
free rank $1$ $\mathscr{A}_{z}$-module
$\mathscr{M}_{z}=\operatorname{\mathrm{Span}}_{\mathbb{C}}\\{z^{-2},z^{-1}\\}\oplus
z^{2}\mathbb{C}[z]$. Then

$\mathfrak{F}_{z}(\Psi_{2})=\\{D(z,\partial_{z}):D(z,\partial_{z})\cdot\mathscr{M}_{z}\subseteq\mathscr{M}_{z}\\}$

is the algebra of differential operators on the line bundle $\mathcal{L}$ over
$Z$ associated to $\mathscr{M}_{z}$.

The generalized Fourier map $b_{\Psi}$ may be described in terms of
$b_{\Psi_{\operatorname{{\mathrm{Ai}}}}}$ by

$b_{\Psi}(A(x,\partial_{x}))=\frac{1}{q_{2}(z)}b_{\Psi_{\operatorname{{\mathrm{Ai}}}}}\left[P_{2}(x,\partial_{x})^{*}\frac{1}{p_{2}(x)}A(x,\partial_{x})\frac{1}{p_{2}(x)}P_{2}(x,\partial_{x})\right]\frac{1}{q_{2}(z)}.$

Let $T_{2}$ be the integral operator

$T_{2}:f(z)\mapsto\int_{t_{1}}^{\infty}K_{2}(z,w)f(w)dw,\quad
K_{2}(z,w)=\int_{t_{2}}^{\infty}\Psi_{2}(x,z)\Psi_{2}(x,w)dx.$

The specific value of the kernel $K_{2}(z,w)$ is determined via integration by
parts to be

$K_{2}(z,w)=\frac{q_{2}(w)}{q_{2}(z)}K_{\operatorname{{\mathrm{Ai}}}}(z,w)+\mathcal{C}_{P_{2}}(\psi_{\operatorname{{\mathrm{Ai}}}}(x,z),\psi_{2}(x,w)/p_{2}(x);t_{2}).$

Computer calculation finds
$\dim\mathfrak{F}_{x,\operatorname{\mathrm{sym}}}^{10,10}(\Psi_{2})=32$, and
therefore $T_{2}$ will commute with a differential operator
$S_{2}(z,\partial_{z})$ in
$\mathfrak{F}_{z,\operatorname{\mathrm{sym}}}^{10,10}(\Psi_{2})$.

Taking $t_{1}=t_{2}=1$, and solving the linear system describing the vanishing
of the concomitants, we find differential operators of order $10$, $12$, $14$,
$16$, and $18$ commuting with $T_{2}$. The operators $S_{2}(z,\partial_{z})$
and $\widetilde{S}_{2}(z,\partial_{z})$ of order $10$ and $12$ are given by

$S_{2}(z,\partial_{z})=\frac{1}{z^{2}}\left(\sum_{k=0}^{5}\partial_{z}^{k}(1-z)^{k}a_{k}(z)\partial_{z}^{k}\right)\frac{1}{z^{2}},$
$\displaystyle a_{0}(z)$
$\displaystyle=z^{14}-200z^{11}+170z^{10}+5640z^{8}-7360z^{7}+2160z^{6}-11520z^{5}-2880z+4320,$
$\displaystyle a_{1}(z)$
$\displaystyle=5z^{12}-580z^{9}+380z^{8}+6240z^{6}-3700z^{5}-960z^{2}-9600z+4800,$
$\displaystyle a_{2}(z)$
$\displaystyle=10z^{10}-560z^{7}+180z^{6}+960z^{4}+1800z^{3}+300z^{2},$
$\displaystyle a_{3}(z)$ $\displaystyle=10z^{8}-180z^{5}-100z^{4}-420z+1260,$
$\displaystyle a_{4}(z)$ $\displaystyle=5z^{6}-70z^{2},$ $\displaystyle
a_{5}(z)$ $\displaystyle=z^{4};$
$\widetilde{S}_{2}(z,\partial_{z})=\frac{1}{z^{2}}\left(\sum_{k=0}^{6}\partial_{z}^{k}(1-z)^{k}\widetilde{a}_{k}(z)\partial_{z}^{k}\right)\frac{1}{z^{2}},$
$\displaystyle\widetilde{a}_{0}(z)$
$\displaystyle=z^{16}-340z^{13}+504z^{12}+21040z^{10}-52200z^{9}+28812z^{8}$
$\displaystyle\ -192000z^{7}+490464z^{6}-328320z^{5}-201600z+130464,$
$\displaystyle\widetilde{a}_{1}(z)$
$\displaystyle=6(z^{14}-220z^{11}+300z^{10}+7000z^{8}-14212z^{7}+5148z^{6}$
$\displaystyle\ -16800z^{5}+13568z^{4}+13568z^{3}+2368z^{2}-6240z+12480,$
$\displaystyle\widetilde{a}_{2}(z)$
$\displaystyle=3(5z^{12}-640z^{9}+760z^{8}+7800z^{6}-8792z^{5}-2996z^{4}-3120z^{2}-36000z+50400,$
$\displaystyle\widetilde{a}_{3}(z)$
$\displaystyle=4z^{2}(5z^{8}-310z^{5}+270z^{4}+600z^{2}+1566z-2268),$
$\displaystyle\widetilde{a}_{4}(z)$
$\displaystyle=3(5z^{8}-100z^{5}-224z+784),$
$\displaystyle\widetilde{a}_{5}(z)$ $\displaystyle=6(z-2)z^{2}(z+2)(z^{2}+4),$
$\displaystyle\widetilde{a}_{6}(z)$ $\displaystyle=z^{4}.$

From Burchnall-Chaundy Theory and its extensions (see for example [26]), we
know that each pair of operators must satisfy a polynomial relation. Together,
the algebra they generate is the coordinate ring of an affine curve. However,
the precise relations that are satisfied are sufficiently complicated so as to
be omitted from the paper.

Acknowledgements. The research of W.R.C. was supported by CSU Fullerton RSCA
intramural grant 0359121. The research of M.Y. was supported by NSF grant
DMS-2131243. The research of I.Z. was supported by CONICET PIP grant
112-200801-01533.

## References

  * [1] B. Bakalov, E. Horozov, and M. Yakimov, General methods for constructing bispectral operators, Phys. Lett. A 222 (1996), 59–66.
  * [2] B. Bakalov, E. Horozov, and M. Yakimov, Bispectral algebras of commuting ordinary differential operators, Comm. Math. Phys. 190 (1997), 331–373.
  * [3] B. Bakalov, E. Horozov, and M. Yakimov, Automorphisms of the Weyl algebra and bispectral operators, in: The bispectral problem, CRM Proc. Lecture Notes, pp. 3–10, AMS 1998.
  * [4] H. Bateman, On the inversion of a definite integral, Proc. London Math. Soc. 2 (1907), 461–498.
  * [5] M. Bertero and F. A. Grünbaum, Commuting differential operators for the finite Laplace transform, Inverse Problems 1 (1985), 181–192.
  * [6] W. R. Casper, F. A. Grünbaum, M. Yakimov and I. Zurrián, Reflective prolate-spheroidal operators and the KP/KdV equations, Proc. Natl. Acad. Sci. USA 116 (2019), no. 37, 18310–18315.
  * [7] W. R. Casper, F. A. Grünbaum, M. Yakimov and I. Zurrián, Reflective prolate-spheroidal operators and the Adelic Grassmannian, arXiv:2003.11616, to appear in Comm. Pure Appl. Math.
  * [8] W. R. Casper and M. Yakimov,  Integral operators, bispectrality and growth of Fourier algebras, J. Reine Angew. Math. 766 (2020), 151–194.
  * [9] M. Castro and F. A. Grünbaum, Time-and-band limiting for matrix orthogonal polynomials of Jacobi type, Random Matrices Theory Appl. 6, No. 04, 1740001, (2017), 12 p.
  * [10] A. Connes and C. Consani, Spectral triples and zeta-cycles, arXiv:2106.01715, 2021.
  * [11] A. Connes and M. Marcoli, Noncommutative geometry, quantum fields and motives, AMS Colloq. Pub. vol 55.
  * [12] A. Connes and H. Moscovici. Prolate spheroidal operator and zeta, arXiv:2112.05500, 2021.
  * [13] J. J. Duistermaat and F. A. Grünbaum, Differential equations in the spectral parameter, Comm. Math. Phys. 103 (1986), 177–240.
  * [14] Fuchs J. Math Anal, Appl. 9 (1964) pp 317–330.
  * [15] F. A. Grünbaum, Eigenvectors of a Toeplitz matrix: discrete version of the prolate spheroidal wave functions, SIAM J. Algebraic Discrete Methods 2 (1981), 136–141.
  * [16] F. A. Grünbaum, Band-time-band limiting integral operators and commuting differential operators, Algebra i Analiz 8 (1996), 122–126.
  * [17] F. A. Grünbaum, Commuting Integral and differential operators and the master symmetries of the Korteweg-de Vries equation, Inverse Problems 37 (2021), no. 8, Paper No. 085010, 21 pp.
  * [18] F. A. Grünbaum and L. Haine, The Wilson bispectral involution: some elementary examples, In: Canterbury meeting, P. Clarkson and F. Nijhoff eds, Cambridge Univ Press, pp 353-368,1999.
  * [19] F. A. Grünbaum, L. Longhi and M. Perlstadt, Differential operators commuting with finite convolution integral operators: some nonabelian examples, SIAM J. Appl. Math. 42 (1982), 941–955.
  * [20] F. A. Grünbaum, I. Pacharoni and I. Zurrián, Bispectrality and time-band limiting: matrix-valued polynomials, Int. Math. Res. Not. 2020, No. 13 (2020), 4016–4036.
  * [21] P. Iliev, $q$-KP hierarchy, bispectrality and Calogero–Moser systems, J. Geom. Phys. 35 (2000), 157–182.
  * [22] E. L. Ince, Ordinary Differential Equations, Dover Publ., New York, 1944.
  * [23] A. Kasman and M. Rothstein, Bispectral Darboux transformations: the generalized Airy case, Phys. D 102 (1997), 159–176.
  * [24] J. Kiukas and R. Werner, Maxima violation of Bell inequalities by position measurements, J. Math. Phys. 51 (2010), Art 072105.
  * [25] M. Kontsevich, Intersection theory on the moduli space of curves and the matrix Airy function, Comm. Math. Phys. 147 (1992), 1–23.
  * [26] I. M. Krichever, Integration of nonlinear equations by the methods of algebraic geometry, Funkcional. Anal. i Priložen. 11 (1977), 15–31; transl. Funct. Anal. Appl., 11:1 (1977), 12–26.
  * [27] H. J. Landau and H. O. Pollak, Prolate spheroidal wave functions, Fourier analysis and uncertainty. II, Bell System Tech. J. 40 (1961), 65–84.
  * [28] H. J. Landau and H. O. Pollak, Prolate spheroidal wave functions, Fourier analysis and uncertainty. III, Bell System Tech. J. 41 (1962), 1295–1336.
  * [29] M. L. Mehta, Random matrices, 3rd ed, Elsevier Publ., 2004.
  * [30] A. Osipov, V. Rokhlin, and H. Xiao, Prolate spheroidal wave functions of order zero, Springer Ser. Appl. Math. Sci 187, 2013.
  * [31] M. Van der Put and M.F. Singer, Galois theory of linear differential equations, vol. 328. Springer Science & Business Media, 2012.
  * [32] D. Slepian, Prolate spheroidal wave functions, Fourier analysis and uncertainity. IV. Extensions to many dimensions; generalized prolate spheroidal functions, Bell System Tech. J. 43 (1964), 3009–3057.
  * [33] D. Slepian, Prolate spheroidal wave functions, Fourier analysis, and uncertainty V. The discrete case Bell system technical journal, 52 (1978) pp 1371–1430.
  * [34] D. Slepian, On bandwith, Proc. IEEE 64 (1976), 292–300.
  * [35] D. Slepian, Some comments on Fourier analysis, uncertainty and modelling, SIAM Review 25 (1983), 379–393.
  * [36] D. Slepian and H. O. Pollak, Prolate spheroidal wave functions, Fourier analysis and uncertainty. I, Bell System Tech. J. 40 (1961), 43–63.
  * [37] C. A. Tracy and H. Widom, Fredholm determinants, differential equations and matrix models, Comm. Math. Phys. 163 (1994), 33–72.
  * [38] C. A. Tracy and H. Widom, Level-spacing distributions and the Airy kernel, Comm. Math. Phys. 159 (1994), 151–174.
  * [39] G. Wilson, Bispectral commutative ordinary differential operators, J. Reine Angew. Math. 442 (1993), 177–204.
  * [40] G. Wilson, On the antiplectic pair connected with the Adler-Gelfand-Dikii braket, Nonlinearity 5 (1992), 109–131.
  * [41] E. Witten, Two-dimensional gravity and intersection theory on moduli space, Surveys in diff. geom. (Cambridge, MA, 1990), pp. 243–310, Lehigh Univ. 1991.