# Time and band limiting for exceptional polynomials

Castro M. M., Grünbaum F. A. and Zurrián I  Department of Mathematics,
University of California, Berkeley CA 94705<EMAIL_ADDRESS>Departamento de Matemática Aplicada II, Universidad de Sevilla, Seville 41011
<EMAIL_ADDRESS>Departamento de Matemática Aplicada II, Universidad de Sevilla,
Seville 41011<EMAIL_ADDRESS>
###### Abstract.

The “time-and-band limiting” commutative property was found and exploited by
D. Slepian, H. Landau and H. Pollak at Bell Labs in the 1960’s, and
independently by M. Mehta and later by C. Tracy and H. Widom in Random matrix
theory. The property in question is the existence of local operators with
simple spectrum that commute with naturally appearing global ones.

Here we give a general result that insures the existence of a commuting
differential operator for a given family of exceptional orthogonal polynomials
satisfying the “bispectral property”. As a main tool we go beyond
bispectrality and make use of the notion of Fourier Algebras associated to the
given sequence of exceptional polynomials. We illustrate this result with two
examples, of Hermite and Laguerre type, exhibiting also a nice Perline’s form
for the commuting differential operator.

###### Key words and phrases:

Time-band limiting, exceptional polynomials

###### 2010 Mathematics Subject Classification:

33C45, 22E45, 33C47

##  1\. Introduction

The time-and-band limiting problem has its origins in signal processing, and
was motivated by a question posed by Claude Shannon in [35]: what is the best
use one can make of the values of the Fourier transform $\mathcal{F}f(k)$, for
values of $k$ in the band $\left[-\mathcal{W},\mathcal{W}\right]$, to recover
the time-limited signal $f(x)$ with finite support
$\left[-\mathcal{T},\mathcal{T}\right]$. The solution was given by three
workers at Bell Labs during the 1960s: David Slepian, Henry Landau, and Henry
Pollak. See [40, 26, 27, 36, 38, 39, 37].

These papers show the central role of an integral operator $S$\- with a kernel
built from the exponentials that feature in the Fourier transform - which
exhibits a “spectral gap”. Only about $2\mathcal{W}\times 2\mathcal{T}$ of its
eigenvalues are truly nonzero and one should find the projection of $f(x)$ on
the span of the corresponding eigenfunctions. Any effort to compute the
components of $f(x)$ along the remaining eigenfunctions will result in
numerical instability. In practice the eigenfunctions of the integral operator
in question can only be computed because they happen to be the eigenfunctions
of a commuting differential operator $T$. This unexpected algebraic accident
has played a role in other areas, even as far afield as the study of the
Riemann zeta function, see [11]. For a further discussion of this and other
historical points the reader can consult [10].

The “bispectral problem” introduced in [12] is an effort to understand the
reason for the existence of this commuting operator $T$ and to extend this
commutativity miracle. The are numerous papers linking bispectrality with the
building of a commuting operator, the most recent ones may be [3, 4, 5, 6,
10].

One has a “bispectral situation” when a function $f(x,n)=q_{n}(x)$ is a
simultaneous eigenfunction of an operator $D$, acting in the variable $x$, and
an operator $L$, acting in the variable $n$. Namely,

$Dq_{n}(x)=\lambda_{n}q_{n}(x),\qquad Lq_{n}(x)=\mu_{x}q_{n}(x),$

where the eigenvalue $\lambda_{n}$ (respectively $\mu_{x}$) depends only on
$n$ (respectively $x$). Both variables may be discrete or continuous. In this
work we shall consider a continuous-discrete setup, given by a family of
orthogonal polynomials that are eigenfunctions of a common differential
operator. Previous works in the same direction include [19], where the
situation of the classical orthogonal polynomials is considered, and for
matrix valued orthogonal polynomials we point out [22, 7, 8, 9, 20, 21]. As a
continuation of the work started in [10], in this paper we deal with sequences
of exceptional orthogonal polynomials in connection with time-and-band
limiting. We manage to get results that were not reached in [10].

The relationship between exceptional polynomials (to be defined in Section 3),
bispectrality and the use of the Darboux process has been considered by
different authors, [17, 32, 34, 16, 18, 15, 24, 6, 10]. The study of
exceptional orthogonal polynomials has grown into a very active area in recent
years. The literature is large and we just mention a few papers, [13, 14, 15,
16, 17, 18, 24, 32, 34]. A more detailed account is given in [10].

##  2\. The contents of the paper

The case of exceptional Jacobi Polynomials was already covered in [10, Section
5], where an effort is made to write down the commuting operator in “Perline
form” as explained now.

Quite a while back the case of polynomials defined on finite sets was
considered in [30, 31, 29]. In [29] R. Perline obtained a simple form for the
commuting local operator. It was subsequently seen, see [20, 23] that this
simple form of Perline applies unchanged to other situations. In [23] one
shows that a more complicated form of the commuting operator, still built with
an appropriate extension of Perline’s construction, yields one with simple
spectrum.

We will see in this paper that to obtain a Perline form of the commuting
operator one has to go beyond bispectrality and make use of the more general
notion of Fourier Algebras (see Section 4 below), first considered in [1] in
connection with bispectral meromorphic functions and the commutativity
property and in [2] in the framework of matrix-valued orthogonal polynomials.
This concept was also exploited in [3, 4, 5, 6].

Specifically, in Theorem 4.2 we give a general result that guarantees the
existence of a commuting differential operator $T$ for an integral operator
$S$ built from a family of exceptional orthogonal polynomials satisfying the
bispectral conditions (1) below. The main tool turns out to be the computation
of the dimensions of certain subspaces of the Fourier Algebras associated with
the given sequence of exceptional polynomials (see Section 4). The method of
considering the Fourier algebras to find a commuting operator is not new and
was used already in [6]. Here, we dive into a more detailed study and we
proceed with the considerations needed in the particular case given by
exceptional orthogonal polynomials, which do not fit neatly with the general
results.

##  3\. Preliminaries

In this paper we deal with sequences of real valued polynomials
{$q_{n}(x)\\}_{n\in\mathbb{N}_{0}}$ for which there exist a nontrivial
differential operator $D$ and a nontrivial difference operator $L$ such that:

(1) $Dq_{n}(x)=\lambda_{n}q_{n}(x),\qquad Lq_{n}(x)=\mu_{x}q_{n}(x).$

We allow $q_{n}\equiv 0$ only for $n$ in a finite set $X$, these are the so
called exceptional degrees. On the other hand, we require $\deg({q_{n}})=n$
for $n\in Z:=\mathbb{N}_{0}\setminus X$ and {$q_{n}(x)\\}_{n\in Z}$ to be a
complete orthonormal system with respect to a weight $w(x)$ defined on a
(possibly infinite) open interval $(a,b)$. In the particular situation when
$\operatorname{ord}(D)=2$, we are dealing with what are nowadays called
“exceptional orthogonal polynomials”.

We consider the Hilbert spaces $\ell^{2}(Z)$ and
$L^{2}(w)=L^{2}((a,b),w(x)dx)$ given by the real valued sequences
$\\{\gamma_{n}\\}_{n\in Z}$ such that $\sum_{n\in Z}\gamma_{n}^{2}<\infty$ and
all measurable functions $f(x)$, $x\in(a,b)$, satisfying
$\int_{a}^{b}|f(x)|^{2}w(x)\,dx<\infty$, respectively. A natural analog of the
Fourier transform is $F:\ell^{2}(Z)\longrightarrow L^{2}(w)$ given by

(2) $\\{\gamma_{n}\\}_{n\in Z}\overset{F}{\longrightarrow}\sum_{n\in
Z}\gamma_{n}q_{n}(x).$

Its adjoint $F^{*}:L^{2}(w)\longrightarrow\ell^{2}(Z)$ is given by

$f\overset{F^{*}}{\longrightarrow}\gamma_{n}=\int_{a}^{b}f(x)\,w(x)\,q_{n}(x)dx.$

In the spirit of Shannon we consider the problem of determining a function $f$
from the following data: $f$ is supported on the compact set $\\{n\in
Z:\,n\leq N\\}$, for a value $N\in Z$, and its transform $Ff$ is known on a
compact set $[a,\Omega]$, for $\Omega\in(a,b)$. This leads us to compute the
singular vectors (and values) of the operator $E:\ell^{2}(Z)\longrightarrow
L^{2}(w)$ defined by

$Ef=\chi_{\Omega}F\chi_{N}f,$

where $\chi_{N}$ represents the time-limiting operator on $\ell^{2}(Z)$ and
$\chi_{\Omega}$ represents the band-limiting operator on $L^{2}(w)$. Namely,
$\chi_{N}$ acts on $\ell^{2}(Z)$ by zeroing out all components with indices
greater than $N$. On the other hand, $\chi_{\Omega}$ acts on $L^{2}(w)$ by
multiplying by the characteristic function of the interval $(a,\Omega)$.

We are thus required to study the eigenvectors of the operators

$E^{*}E=\chi_{N}F^{*}\chi_{\Omega}F\chi_{N}\qquad\text{ and }\qquad
EE^{*}=\chi_{\Omega}F\chi_{N}F^{*}\chi_{\Omega}.$

The operator $E^{*}E$, acting on $\ell^{2}(Z)$, is just a finite dimensional
matrix with each entry given by

(3) $(E^{*}E)_{m,n}=\int_{a}^{\Omega}q_{m}(x)w(x)q_{n}(x)dx,\qquad m,n\in
Z,\quad m,n\leq N.$

The operator $EE^{*}$ acts on $L^{2}((a,\Omega),w(x)dx)$ by means of the
integral kernel

(4) $k(x,y)=\sum_{n\in Z,\,n\leq N}q_{n}(x)q_{n}(y),$

namely, the integral operator $S=EE^{*}$ is given by

(5) $(Sf)(x)=\int_{a}^{\Omega}f(y)w(y)k(x,y)dy.$

This is the analog of the integral operator considered by Slepian, Landau and
Pollak.

For general $N$ and $\Omega$, finding the eigenfunctions of $EE^{*}$ and
$E^{*}E$ analytically is an impossible task. However, a miracle comes to our
rescue: we identify a local operator with simple spectrum that shares the same
eigenfunctions as the operators $EE^{*}$ or $E^{*}E$. This approach was
successfully employed by Slepian, Landau, and Pollak in the scalar case when
dealing with traditional Fourier analysis. They made the following
discoveries:

  * •

For each $N$ and $\Omega$, a symmetric tridiagonal matrix $\hat{T}$ exists
with a simple spectrum that commutes with $E^{*}E$.

  * •

For each $N$ and $\Omega$, there exists a self-adjoint second-order
differential operator $T$ with a simple spectrum that commutes with the
integral operator $S=EE^{*}$.

In the case of traditional Fourier analysis, the operators $D$ and $L$
involved in (1) may be taken of very low order. In this work, we need to
consider arbitrary orders, which leads us to the construction of a higher
order differential operator $T$.

For a full discussion of important numerical aspects in the traditional
Fourier case, see [33].

##  4\. The Fourier Algebras

A difference operator of order $2k\geq 0$ is denoted by
${\tilde{L}}=\sum_{j=-k}^{k}{\tilde{L}}_{j}\,\delta^{j}$, with
${\tilde{L}}_{k}$ or ${\tilde{L}}_{-k}$ nonzero, where $\delta$ is the usual
shift. The operator ${\tilde{L}}$ acts on functions
$h:\mathbb{N}_{0}\longrightarrow\mathbb{R}$ by

(6) $({\tilde{L}}h)_{n}=\sum_{j=-k}^{k}{\tilde{L}}_{n,n+j}\,h_{n+j},$

with the understanding that ${\tilde{L}}_{n,j}=0$ for $j<0$.

Let us consider the Fourier algebras associated with $q_{n}(x)$ defined by

$\mathcal{F}_{x}=\\{\text{differential operators
}\tilde{D}:\tilde{D}q_{n}(x)=\tilde{L}q_{n}(x)\text{ for some difference
operator }\tilde{L},\text{ for all }n,\,x\\},$
$\mathcal{F}_{n}=\\{\text{difference operators
}\tilde{L}:\tilde{D}q_{n}(x)=\tilde{L}q_{n}(x)\text{ for some differential
operator }\tilde{D},\text{ for all }n,\,x\\}.$

By a slight abuse of notation, given any $\tilde{L}\in\mathcal{F}_{n}$ we will
also denote by $\tilde{L}$ the semi-infinite $(2k+1)$-diagonal matrix whose
entries are given by (6). We say that the operator $\tilde{L}$ is symmetric if
the corresponding matrix is symmetric. For any $x\in\mathbb{(}a,b)$, let us
consider the semi-infinite vector
$\mathbf{q}(x):=\left(q_{0}(x),q_{1}(x),\dots\right)^{T}$. Then, if we have
$\tilde{D}q_{n}(x)=\tilde{L}q_{n}(x)$ for some $\tilde{D}\in\mathcal{F}_{x}$,
in vector notation we will have

$\tilde{D}\mathbf{q}(x)=\tilde{L}\mathbf{q}(x).$

It is worth noticing that we have made an important notational choice: the
vector $\mathbf{q}(x)$ has zeros in every entry corresponding to an
exceptional degree. The reader will realize that we could have skipped all
these entries and that the results and statements will hold with the due
technical adaptations. Nevertheless, we have made this choice resulting in a
more succinct and shorter presentation, not only in the proofs in this
section, but also in the examples in Sections 5 and 6. The matrix
representation of $\tilde{L}$ in the basis obtained by including only the
$q_{n}$ with $n$ in $Z$ is obtained by “compressing” the matrices going with
the choice above, i.e. we need to throw away the rows and columns
corresponding to the exceptional degrees.

From (1) we know that $D\in\mathcal{F}_{x}$ and, since $\\{q_{n}(x)\\}_{n\in
Z}$ is a orthogonal complete system, that $D$ is symmetric in $L^{2}(w)$. This
leads us to consider, for given $d$ and $k$ in $\mathbb{N}_{0}$, the subspaces

$\mathcal{F}_{x,sym}^{d,k}=\left\\{\text{Differential operators
}\tilde{D}:\begin{aligned} \tilde{D}q_{n}(x)=\tilde{L}q_{n}(x)\text{ for a
symmetric difference operator }\tilde{L},\\\ \text{ with
}\operatorname{ord}(\tilde{D})\leq 2d\text{ with
}\operatorname{ord}(\tilde{L})\leq 2k.\end{aligned}\right\\}$

Notice that if $\tilde{D}q_{n}(x)=\tilde{L}q_{n}(x)$ with $\tilde{L}$
symmetric, since $\\{q_{n}(x)\\}_{n\in Z}$ is a orthogonal complete system,
$\tilde{D}$ is a symmetric differential operator in $L^{2}(w)$. Hence, all the
operators in $\mathcal{F}_{x,sym}^{d,k}$ are symmetric in $L^{2}(w)$. Recall
that an operator $\tilde{D}$ is symmetric in $L^{2}(w)=L^{2}((a,b),w(x)dx)$ if
for every $f$ and $g$ in an appropriate dense domain we have

$\int_{a}^{b}\tilde{D}f(x)w(x)g(x)dx=\int_{a}^{b}f(x)w(x)\tilde{D}g(x)dx.$

One defines in an analogous manner the notion of symmetry for a differential
operator $\tilde{D}$ in $L^{2}((a,\Omega),w(x)dx)$.

The identity above requires the vanishing of several terms resulting from
integration by parts. This can come about either with the help of certain
coefficients in the operator $\tilde{D}$ or because the weight $w(x)$ and some
of its derivatives vanish at the end points. In the absence of this, one needs
to restrict the domain of the operator to functions $f(x),g(x)$ which, maybe
along with their derivatives, vanish at these end points. One should keep this
remark in mind when reading the comment after the proof of the lemma below.

Now we give the general results that allow us to find commuting operators in
our time-and-band limiting setup.

######  Lemma 4.1.

If $T\in\mathcal{F}_{x,sym}^{d,k}$ is also symmetric in
$L^{2}((a,\Omega),w(x)dx)$ and $\deg(Tq_{n}(x))\leq N$ for any $n\leq N$, then
$T$ commutes with the integral operator $S$ in (5).

###### Proof.

If we denote by $I_{N}$ the semi-infinite diagonal matrix whose entries
$(I_{N})_{j,j}$ are $1$ for $0\leq j\leq N$ and zero elsewhere, we have

$k(x,y)=\sum_{n\in Z,\,n\leq
N}q_{n}(x)q_{n}(y)=\mathbf{q}(x)^{T}\,I_{N}\,\mathbf{q}(y).$

By hypothesis, there is a symmetric semi-infinite $(2k+1)$-diagonal matrix
$\hat{T}$ such that $T\mathbf{q}(x)=\hat{T}\mathbf{q}(x)$. This, combined with
the fact $\deg(Tq_{n}(x))\leq N$ for any $n\leq N$, gives us
$\hat{T}\,I_{N}=I_{N}\,\hat{T}$. Hence,

(7)
$T(x)(k(x,y))=({\hat{T}}\mathbf{q}(x))^{T}\,I_{N}\,\mathbf{q}(y)=\mathbf{q}(x)^{T}\,I_{N}\,{\hat{T}}\mathbf{q}(y)=T(y)(k(x,y)).$

Notice that to avoid ambiguity in the previous identities, we have used $T(x)$
and $T(y)$ to denote the action of the differential operator in the variable
$x$ and $y$, respectively.

Now, let $f:(a,b)\longrightarrow\mathbb{R}$ be a sufficiently differentiable
function in the domain of $T$, then

$\displaystyle T(S(f)(x))$
$\displaystyle=\int_{a}^{\Omega}f(y)w(y)T(x)k(x,y)dy\stackrel{{\scriptstyle
i}}{{=}}\int_{a}^{\Omega}f(y)w(y)T(y)k(x,y)dy$
$\displaystyle\stackrel{{\scriptstyle
ii}}{{=}}\int_{a}^{\Omega}Tf(y)w(y)k(x,y)dy{=}(S(Tf)(x)),$

with $(i)$ given by (7) and $(ii)$ a consequence of $T$ being symmetric with
respect to $w$ in $(a,\Omega)$. This completes the proof.∎

One more comment before our main result. Since the integral operator $S$ has
as its domain all of $L^{2}((a,\Omega),w(x)dx)$ it is important that the
commuting differential operator that we produce should have a large domain
that includes the eigenfunctions of $S$. These eigenfunctions, or their
derivatives, do not vanish at either end point. The weight $w(x)$ will not
vanish at $\Omega$ and thus the symmetry of $T$ will require conditions such
as (9) below.

######  Theorem 4.2.

If for some $d,k$ in $\mathbb{N}_{0}$, we have that

(8) $\dim\left(\mathcal{F}_{x,sym}^{d,k}\right)>\frac{k(k+1)+d(d+1)}{2}+1,$

then there exists a nonscalar operator $T\in\mathcal{F}_{x,sym}^{d,k}$ that
commutes with the integral operator $S$ in (5).

###### Proof.

Given an arbitrary $T\in\mathcal{F}_{x,sym}^{d,k}$ (with real coefficients),
we can express it as

$T=w(x)^{-1/2}\left(\sum_{j=0}^{d}\,\partial_{x}^{j}\,t_{j}(x)\,\partial_{x}^{j}\right)w(x)^{1/2},$

for some functions $t_{j}$. Given that it is symmetric with respect to $w$ in
$(a,b)$, if we also want it to be symmetric in $(a,\Omega)$ it is enough to
impose the following $\frac{d(d+1)}{2}$ linear conditions:

(9) $t_{j}^{(i)}(\Omega)=0,\quad\text{ for }j=1,\dots,d,\quad i=0,\dots,j-1.$

Since $T\in\mathcal{F}_{x,sym}^{d,k}$, we have a symmetric semi-infinite
$(2k+1)$-diagonal matrix $\hat{T}$ such that
$T\mathbf{q}(x)=\hat{T}\mathbf{q}(x)$,

(10) $\hat{T}=\left[\begin{array}[]{ccccc
;{2pt/2pt}ccc}{\hat{T}}_{00}&{\hat{T}}_{01}&\dots&{\hat{T}}_{0k}&&&&\\\
{\hat{T}}_{10}&\ddots&\ddots&\ddots&\ddots&&&\\\
\vdots&\ddots&\ddots&\ddots&\ddots&{\hat{T}}_{N+1-k,N+1}&&\\\
{\hat{T}}_{k0}&\ddots&\ddots&\ddots&\ddots&\vdots&\ddots&\\\
&\ddots&\ddots&\ddots&\ddots&{\hat{T}}_{N,N+1}&\ldots&{\hat{T}}_{N,N+k}\\\
\hdashline[2pt/2pt]&&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots\end{array}\right].$

Hence, if we want $\deg(Tq_{n}(x))\leq N$ for any $n\leq N$, it is enough to
satisfy the following $\frac{k(k+1)}{2}$ linear conditions:

(11) $\hat{T}_{i,j}=0,\quad\text{ for }j=N+1,\dots,N+k,\quad i=j-k,\dots,N,$

i.e., the elements of the top right-hand side of the matrix (10) should
vanish. We understand that $\hat{T}_{i,j}=0$ if $i$ is negative.

Then, for $\mathcal{F}_{x,sym}^{d,k}$ to contain a nonzero operator $T$
satisfying (9) and (11), it is enough that its dimension should be at least
$\frac{k(k+1)+d(d+1)}{2}+1$. If we require that $T$ should not be a scalar
multiplication, it is sufficient that the dimension should be larger than
$\frac{k(k+1)+d(d+1)}{2}+1$. By Lemma 4.1, any operator $T$ satisfying (9) and
(11) commutes with $S$. This completes the proof which is in part an extension
of the one in [29]. ∎

In the following sections, we consider two different families of exceptional
orthogonal polynomials. We apply Theorem 4.2 to show the existence of a
commuting differential operator, and we use the concepts in Lemma 4.1 to build
it explicitly.

A careful look at the appropriate selfadjoint extension of the symmetric
commuting differential operator $T$ has, to the best of our knowledge, only
been addressed in [25] and only for the case of traditional Fourier analysis.

##  5\. Hermite Exceptional Polynomials

Let us consider the sequence of Hermite exceptional polynomials given by

$q_{n}(x)=\frac{\sqrt{(n-1)(n-2)}}{\sqrt[4]{\pi}\sqrt{2^{n}n!}}\left(H_{n}(x)+4nH_{n-2}(x)+4n(n-3)H_{n-4}(x)\right),\quad
n\in\mathbb{N}_{0},$

where $H_{n}$ are the classical Hermite polynomials given by the Rodrigues
formula

$H_{n}=(-1)^{n}e^{x^{2}}D_{x}^{n}e^{-x^{2}}.$

We have that $q_{n}\equiv 0$ only for $n$ in $X=\\{1,2\\}$ and
{$q_{n}(x)\\}_{n\notin X}$ is a complete orthonormal system with respect to
the weight

$w(x)=\frac{e^{-x^{2}}}{(1+2x^{2})^{2}},\quad x\in(-\infty,\infty).$

Now we use Theorem 4.2 to find a differential operator $T$ that commutes with
the operator $S$ given by

$(Sf)(y)=\int_{-\infty}^{\Omega}\left(\sum_{k=0}^{N}q_{k}(x)q_{k}(y)\,w(x)\right)\,f(x)\,dx.$

Namely, we look for $d$ and $k$ satisfying (8). Solving the system of
equations $\tilde{D}q_{n}(x)=\tilde{L}q_{n}(x)$, with $\deg(\tilde{D})\leq 2$
and $\deg(\tilde{T})\leq 6$, one finds that $\mathcal{F}_{x,sym}^{3,3}$
contains the differential operators

$\displaystyle D_{0}$
$\displaystyle=\dfrac{d^{2}}{dx^{2}}-\frac{2\,{\left(2\,x^{3}+5\,x\right)}}{2\,x^{2}+1}\dfrac{d}{dx},$
$\displaystyle D_{1}$
$\displaystyle=x\dfrac{d^{2}}{dx^{2}}-\frac{4\,x^{4}+8\,x^{2}-1}{2\,x^{2}+1}\dfrac{d}{dx},$
$\displaystyle D_{2}$
$\displaystyle=\left(-\frac{1}{56}\,x^{2}+\frac{17}{112}\right)\dfrac{d^{2}}{dx^{2}}+\left(\frac{4\,x^{5}-28\,x^{3}-87\,x}{56\,{\left(2\,x^{2}+1\right)}}\right)\dfrac{d}{dx}+1,$
$\displaystyle D_{3}$ $\displaystyle=\frac{4x^{3}}{3}+2x,$

with $D_{j}q_{n}(x)=L_{j}q_{n}(x)$ for

$\displaystyle L_{0}=$ $\displaystyle-2n,$ $\displaystyle L_{1}=$
$\displaystyle\sqrt{2}\sqrt{n-1}\sqrt{n-3}\sqrt{n}\delta^{-1}+\sqrt{2}\sqrt{n+1}\sqrt{n-2}\sqrt{n}\delta,$
$\displaystyle L_{2}=$
$\displaystyle\frac{\sqrt{{\left(n-2\right)}{\left(n-3\right)}{\left(n-4\right)}n}}{56}\delta^{-2}+\frac{{{\left(n-4\right)}{\left(n-7\right)}}}{28}+\frac{\sqrt{{\left(n+2\right)}{\left(n-1\right)}{\left(n-2\right)}n}}{56}\delta^{2},$
$\displaystyle L_{3}=$
$\displaystyle\frac{\sqrt{2}}{3}\sqrt{n(n-4)(n-5)}\delta^{-3}+\sqrt{2n(n-1)(n-3)}\,\delta^{-1}+\sqrt{2(n+1)n(n-2)}\,\delta$
$\displaystyle+\frac{\sqrt{2}}{3}\sqrt{(n+3)(n-1)(n-2)}\delta^{3}.$

Furthermore, one can verify that the following $14$ operators are linearly
independent:

$\displaystyle Y$
$\displaystyle=\\{Id,D_{0},D_{1},D_{2},D_{3},D_{0}^{2},D_{0}^{3},D_{1}^{2},D_{3}D_{0}+D_{0}D_{3},D_{3}D_{0}^{2}+D_{0}^{2}D_{3},D_{0}D_{1}+D_{1}D_{0},D_{3}D_{0}^{3}+D_{0}^{3}D_{3},$
$\displaystyle
D_{0}^{2}D_{1}+D_{1}D_{0}^{2},D_{0}D_{1}^{2}+D_{1}^{2}D_{0}\\}.$

Since $Y\subset\mathcal{F}_{x,sym}^{3,3}$ , and
$14>\frac{3(3+1)+3(3+1)}{2}+1$, by Theorem 4.2 we know that there exists an
operator $T\in\mathcal{F}_{x,sym}^{3,3}$ that commutes with $S$. To find such
a $T$, we consider a linear combination of the operators in $Y$ that satisfies
both hypothesis of Lemma 4.1. After solving the corresponding system of linear
equations one obtains that

$\displaystyle T=$
$\displaystyle-\,\Omega^{3}D_{0}^{3}+\tfrac{3}{4}\left(2\,\Omega^{2}-1\right)\left({D_{1}}D_{0}^{2}+D_{0}^{2}{D_{1}}\right)-\tfrac{3}{2}\,\Omega({D_{1}}^{2}D_{0}+D_{0}{D_{1}}^{2})+\tfrac{3}{8}\left({D_{3}}D_{0}^{3}+D_{0}^{3}{D_{3}}\right)$
$\displaystyle+3(\,N\Omega^{2}+2\,\Omega^{2}-\tfrac{3}{2}\,N-4)({D_{1}}D_{0}+D_{0}{D_{1}})-3\,{\left(2\,\Omega^{2}-1\right)}\Omega\
D_{0}^{2}$
$\displaystyle-2\Omega\,{\left(54\,N^{2}+4\,\Omega^{2}+48\,N-3\right)}D_{0}-6\,{\left(2\,N+1\right)}\Omega\
{D_{1}}^{2}+\tfrac{9}{8}(2\,N+1)({D_{3}}D_{0}^{2}+D_{0}^{2}{D_{3}})$
$\displaystyle+\tfrac{3}{2}(3\,N^{2}+3\,N-7)(D_{0}{D_{3}}+{D_{3}}D_{0})+3(6\,N\Omega^{2}-4\,N^{2}+2\,\Omega^{2}-17\,N-7){D_{1}}$
$\displaystyle+672{\left(N+1\right)}N\Omega\
{D_{2}}+\tfrac{3}{2}{\left(2\,N^{2}+2\,N-15\right)}{\left(2\,N+1\right)}{D_{3}}$

satisfies the desired conditions. Then, $T$ commutes with $S$.

This extension of the Perline representation of a commuting $T$ is not unique,
and a different choice of enough linearly independent operators in
$\mathcal{F}_{x,sym}^{3,3}$ would produce a different representation.

Returning to the operator $T$ given above, if one in interested in an explicit
expression for the coefficients one can use

$T=w^{-1/2}\left(\sum_{k=0}^{3}\partial_{x}^{k}f_{k}\partial_{x}^{k}\right)w^{1/2},$

with

$f_{3}=(x-\Omega)^{3}$
$f_{2}={3}{{\left(2\,x^{2}+1\right)}^{-2}}{\left(\Omega-x\right)}^{2}\,\left(4\,\Omega
x^{6}-4\,x^{7}+8\,Nx^{5}+8\,\Omega x^{4}-12\,x^{5}+8\,Nx^{3}+21\,\Omega
x^{2}-17\,x^{3}+2\,Nx-7\,\Omega+10\,x\right),$
$f_{1}=-{3}{{\left(2\,x^{2}+1\right)}^{-4}}{\left(\Omega-x\right)}\,\left(16\,\Omega^{2}x^{12}-32\,\Omega
x^{13}+16\,x^{14}+64\,N\Omega
x^{11}-64\,Nx^{12}+64\,N^{2}x^{10}+64\,\Omega^{2}x^{10}-160\,\Omega
x^{11}+96\,x^{12}+192\,N\Omega
x^{9}-256\,Nx^{10}+160\,N^{2}x^{8}+168\,\Omega^{2}x^{8}-272\,\Omega
x^{9}+72\,x^{10}+416\,N\Omega
x^{7}-448\,Nx^{8}+160\,N^{2}x^{6}+152\,\Omega^{2}x^{6}-192\,\Omega
x^{7}-120\,x^{8}+160\,N\Omega
x^{5}-128\,Nx^{6}+80\,N^{2}x^{4}-279\,\Omega^{2}x^{4}+334\,\Omega
x^{5}-359\,x^{6}-76\,N\Omega
x^{3}+92\,Nx^{4}+20\,N^{2}x^{2}+1210\,\Omega^{2}x^{2}-1662\,\Omega
x^{3}+324\,x^{4}-36\,N\Omega x+32\,Nx^{2}+2\,N^{2}-83\,\Omega^{2}+610\,\Omega
x-497\,x^{2}-2\,N+18\right),$
$f_{0}={{\left(2\,x^{2}+1\right)}^{-6}}(64\,\Omega^{3}x^{18}-192\,\Omega^{2}x^{19}+192\,\Omega
x^{20}-64\,x^{21}+384\,N\Omega^{2}x^{17}-768\,N\Omega
x^{18}+384\,Nx^{19}+768\,N^{2}\Omega
x^{16}+384\,\Omega^{3}x^{16}-768\,N^{2}x^{17}-1344\,\Omega^{2}x^{17}+1536\,\Omega
x^{18}-576\,x^{19}+512\,N^{3}x^{15}+1920\,N\Omega^{2}x^{15}-4608\,N\Omega
x^{16}+2688\,Nx^{17}+3456\,N^{2}\Omega
x^{14}+624\,\Omega^{3}x^{14}-4224\,N^{2}x^{15}-1296\,\Omega^{2}x^{15}+336\,\Omega
x^{16}+336\,x^{17}+2304\,N^{3}x^{13}+4128\,N\Omega^{2}x^{13}-8640\,N\Omega
x^{14}+4000\,Nx^{15}+50496\,N^{2}\Omega
x^{12}+112\,\Omega^{3}x^{12}-10176\,N^{2}x^{13}+3552\,\Omega^{2}x^{13}-10896\,\Omega
x^{14}+7616\,x^{15}+4224\,N^{3}x^{11}+3840\,N\Omega^{2}x^{11}+37824\,N\Omega
x^{12}-1344\,Nx^{13}+137184\,N^{2}\Omega
x^{10}-1332\,\Omega^{3}x^{10}-12192\,N^{2}x^{11}+12156\,\Omega^{2}x^{11}-29052\,\Omega
x^{12}+18996\,x^{13}+4160\,N^{3}x^{9}-1464\,N\Omega^{2}x^{9}+134928\,N\Omega
x^{10}-10584\,Nx^{11}+165360\,N^{2}\Omega
x^{8}+1200\,\Omega^{3}x^{8}-7440\,N^{2}x^{9}+9804\,\Omega^{2}x^{9}-31272\,\Omega
x^{10}+20940\,x^{11}+2400\,N^{3}x^{7}+11160\,N\Omega^{2}x^{7}+153312\,N\Omega
x^{8}-9176\,Nx^{9}+107784\,N^{2}\Omega
x^{6}+17521\,\Omega^{3}x^{6}-1944\,N^{2}x^{7}-24975\,\Omega^{2}x^{7}-1917\,\Omega
x^{8}+9851\,x^{9}+816\,N^{3}x^{5}+8670\,N\Omega^{2}x^{5}+109500\,N\Omega
x^{6}-12858\,Nx^{7}+39708\,N^{2}\Omega
x^{4}-139845\,\Omega^{3}x^{4}+108\,N^{2}x^{5}+266580\,\Omega^{2}x^{5}-160785\,\Omega
x^{6}+31338\,x^{7}+152\,N^{3}x^{3}-1044\,N\Omega^{2}x^{3}+51324\,N\Omega
x^{4}-9792\,Nx^{5}+7818\,N^{2}\Omega
x^{2}+69189\,\Omega^{3}x^{2}+162\,N^{2}x^{3}-247383\,\Omega^{2}x^{3}+242145\,\Omega
x^{4}-68367\,x^{5}+12\,N^{3}x-1134\,N\Omega^{2}x+11466\,N\Omega
x^{2}-1964\,Nx^{3}+642\,N^{2}\Omega-2097\,\Omega^{3}+24\,N^{2}x+27432\,\Omega^{2}x-56673\,\Omega
x^{2}+29088\,x^{3}+582\,N\Omega+144\,Nx+1134\,\Omega-1512\,x).$

The construction above follows the general arguments given in the previous
section. This is not the only way to produce such a commuting operator $T$.
One can avail oneself of enough many linearly independent difference operators
in the Fourier algebra and, given $N$ and $\Omega$, look for a linear
combination of these that would commute with the matrix of size $N$ of inner
products of the exceptional polynomials in the interval $(-\infty,\Omega)$.
When we replace in this linear combination the difference operators by their
corresponding differential ones in the Fourier algebra we obtain a commuting
$T$.

##  6\. Laguerre Exceptional Polynomials

Let $\alpha>0$ and $L^{(\alpha)}_{n}(x)$, $n\geq 0$, denote the classical
Laguerre polynomials orthogonal with respect to the weight $e^{-x}x^{\alpha}$
in the interval $(0,+\infty)$. We consider the sequence of exceptional
Laguerre polynomials

$q_{n}(x)=\sqrt{\frac{(\alpha+n-1)(n-1)!}{(\alpha+n)\Gamma(n+\alpha)}}\left(-(x+\alpha+1)L_{n-1}^{(\alpha)}+L_{n-2}^{(\alpha)}\right),\quad
n\in\mathbb{N}_{0}.$

We have that $q_{n}\equiv 0$ only for $n=0$ and {$q_{n}(x)\\}_{n>0}$ is a
complete orthonormal system with respect to the weight

$w(x)=\frac{e^{-x}x^{\alpha}}{(x+\alpha)^{2}}.$

In this section, just as we did in the Hermite case, we find a differential
operator $T$ that commutes with the operator $S$ given by

$(Sf)(y)=\int_{0}^{\Omega}\left(\sum_{k=0}^{N}q_{k}(x)q_{k}(y)\,w(x)\right)\,f(x)\,dx.$

In this situation, and proceeding as in the previous example, we find that
$\mathcal{F}_{x,sym}^{2,2}$ contains the differential operators

$\displaystyle D_{0}=$
$\displaystyle\frac{1}{4}\,x\dfrac{d^{2}}{dx^{2}}+\frac{{\left(\alpha+x+1\right)}{\left(\alpha-x\right)}}{4\,{\left(\alpha+x\right)}}\dfrac{d}{dx}+\frac{\alpha+3\,x}{4\,{\left(\alpha+x\right)}},$
$\displaystyle D_{1}=$
$\displaystyle{\left(3\,\alpha-x+10\right)}\frac{x}{42}\dfrac{d^{2}}{dx^{2}}$
$\displaystyle+\frac{3\,\alpha^{3}-\alpha^{2}x-3\,\alpha
x^{2}+x^{3}+13\,\alpha^{2}-5\,\alpha
x-10\,x^{2}+10\,\alpha-10\,x}{42\,{\left(\alpha+x\right)}}\dfrac{d}{dx}$
$\displaystyle-\frac{\alpha^{3}-\alpha x^{2}+5\,\alpha^{2}-3\,\alpha
x+10\,\alpha-10\,x}{42\,{\left(\alpha+x\right)}},$ $\displaystyle D_{2}=$
$\displaystyle(x+\alpha)^{2},$

with $D_{j}q_{n}(x)=L_{j}q_{n}(x)$ for

$\displaystyle L_{0}=$ $\displaystyle(3-n)/4,$ $\displaystyle L_{1}=$
$\displaystyle-\frac{\sqrt{\alpha+n}\sqrt{\alpha+n-1}\sqrt{\alpha+n-2}\sqrt{n-1}}{42}\delta^{-1}+\frac{{\left(n-1\right)}{\left(n-5\right)}}{21}$
$\displaystyle-\frac{\sqrt{\alpha+n+1}\sqrt{\alpha+n}\sqrt{\alpha+n-1}\sqrt{n}}{42}\delta,$
$\displaystyle L_{2}=$
$\displaystyle\sqrt{\alpha+n}\sqrt{\alpha+n-3}\sqrt{n-1}\sqrt{n-2}\,\delta^{-2}-4\,\sqrt{\alpha+n}\sqrt{\alpha+n-1}\sqrt{\alpha+n-2}\sqrt{n-1}\,\delta^{-1}$
$\displaystyle+(4\alpha^{2}+10\,\alpha
n+6\,n^{2}-5\,\alpha-6\,n)-4\sqrt{\alpha+n+1}\sqrt{\alpha+n}\sqrt{\alpha+n-1}\sqrt{n}\,\delta^{1}$
$\displaystyle+\sqrt{\alpha+n+2}\sqrt{\alpha+n-1}\sqrt{n+1}\sqrt{n}\,\delta^{2},$

In this case, we see that $\mathcal{F}_{x,sym}^{2,2}$ already contains at
least $8$ operators :

$Y=\\{Id,D_{0},D_{1},D_{2},D_{0}^{2},D_{1}^{2},D_{0}D_{1}+D_{1}D_{0},D_{2}D_{0}+D_{0}D_{2}\\}.$

One can verify that the operators in $Y$ are linearly independent. Hence,
since $8>\frac{2(2+1)+2(2+1)}{2}+1$, by Theorem 4.2 we know that there exists
an operator $T\in\mathcal{F}_{x,sym}^{2,2}$ that commutes with $S$. As before,
to find such a $T$, we consider a linear combination of the operators in $Y$
that satisfies the hypothesis of Lemma 4.1, obtaining

$\displaystyle T=$ $\displaystyle
1764D_{1}^{2}+168(D_{0}D_{1}+D_{1}D_{0})+16\left(\Omega-3\,\alpha-10\right)^{2}D_{0}^{2}+4\left(N+\alpha\right)\left(D_{0}D_{2}+D_{2}D_{0}\right)$
$\displaystyle+42\,{\left(2\,N\Omega+2\,N\alpha+4\,\alpha^{2}-5\,\Omega+15\,\alpha+40\right)}D_{1}$
$\displaystyle+4\Big{(}N\Omega^{2}-6\,N\Omega\alpha-\Omega^{2}\alpha-7\,N\alpha^{2}+2\,\Omega\alpha^{2}-13\,\alpha^{3}-20\,N\Omega-5\,\Omega^{2}-20\,N\alpha+30\,\Omega\alpha-85\,\alpha^{2}$
$\displaystyle\quad\quad+90\,\Omega-270\alpha-400\Big{)}D_{0}+{\left(N+\alpha\right)}{\left(N-\alpha-5\right)}D_{2}.$

An explicit expression is given by

$T=w^{-1/2}\left(\sum_{k=0}^{2}\partial_{x}^{k}f_{k}\partial_{x}^{k}\right)w^{1/2},$

with

$f_{2}=\left(\Omega-x\right)^{2}x^{2},$
$f_{1}=\left(\Omega-x\right){\left(\alpha+x\right)}^{-2}\Big{(}\Omega\alpha^{4}+2\,N\Omega\alpha^{2}x-2\,\Omega\alpha^{3}x+\alpha^{4}x+4\,N\Omega\alpha
x^{2}-2\,N\alpha^{2}x^{2}-6\,\Omega\alpha^{2}x^{2}+4\,\alpha^{3}x^{2}+2\,N\Omega
x^{3}-4\,N\alpha x^{3}-2\,\Omega\alpha
x^{3}+4\,\alpha^{2}x^{3}-2\,Nx^{4}+\Omega
x^{4}-x^{5}+3\,\Omega\alpha^{3}-11\,\Omega\alpha^{2}x+6\,\alpha^{3}x-19\,\Omega\alpha
x^{2}+30\,\alpha^{2}x^{2}-5\,\Omega x^{3}+30\,\alpha
x^{3}+6\,x^{4}-10\,\Omega\alpha x+38\,\alpha^{2}x-2\,\Omega x^{2}+96\,\alpha
x^{2}+50\,x^{3}\Big{)},$
$f_{0}=\left(\alpha+x\right)^{-3}\Big{(}2\,N^{2}\alpha^{5}-4\,N\Omega\alpha^{5}-4\,N\alpha^{6}-9\,\alpha^{7}+10\,N^{2}\alpha^{4}x-8\,N\Omega\alpha^{4}x-8\,N\alpha^{5}x-25\,\alpha^{6}x+20\,N^{2}\alpha^{3}x^{2}-22\,\alpha^{5}x^{2}+20\,N^{2}\alpha^{2}x^{3}+8\,N\Omega\alpha^{2}x^{3}+8\,N\alpha^{3}x^{3}-6\,\alpha^{4}x^{3}+10\,N^{2}\alpha
x^{4}+4\,N\Omega\alpha
x^{4}+4\,N\alpha^{2}x^{4}-\alpha^{3}x^{4}+2\,N^{2}x^{5}-\alpha^{2}x^{5}+2\,N\Omega^{2}\alpha^{3}-32\,N\Omega\alpha^{4}-2\,\Omega^{2}\alpha^{4}-40\,N\alpha^{5}+13\,\Omega\alpha^{5}-86\,\alpha^{6}+10\,N\Omega^{2}\alpha^{2}x-88\,N\Omega\alpha^{3}x-10\,\Omega^{2}\alpha^{3}x-128\,N\alpha^{4}x+34\,\Omega\alpha^{4}x-249\,\alpha^{5}x+14\,N\Omega^{2}\alpha
x^{2}-80\,N\Omega\alpha^{2}x^{2}-14\,\Omega^{2}\alpha^{2}x^{2}-154\,N\alpha^{3}x^{2}+24\,\Omega\alpha^{3}x^{2}-240\,\alpha^{4}x^{2}+6\,N\Omega^{2}x^{3}-24\,N\Omega\alpha
x^{3}-6\,\Omega^{2}\alpha
x^{3}-90\,N\alpha^{2}x^{3}-2\,\Omega\alpha^{2}x^{3}-82\,\alpha^{3}x^{3}-30\,N\alpha
x^{4}-5\,\Omega\alpha x^{4}-10\,\alpha^{2}x^{4}-6\,Nx^{5}-5\,\alpha
x^{5}-80\,N\Omega\alpha^{3}-7\,\Omega^{2}\alpha^{3}-78\,N\alpha^{4}+97\,\Omega\alpha^{4}-374\,\alpha^{5}-240\,N\Omega\alpha^{2}x-43\,\Omega^{2}\alpha^{2}x-234\,N\alpha^{3}x+293\,\Omega\alpha^{3}x-1058\,\alpha^{4}x-240\,N\Omega\alpha
x^{2}-57\,\Omega^{2}\alpha
x^{2}-234\,N\alpha^{2}x^{2}+311\,\Omega\alpha^{2}x^{2}-946\,\alpha^{3}x^{2}-80\,N\Omega
x^{3}-21\,\Omega^{2}x^{3}-78\,N\alpha x^{3}+115\,\Omega\alpha
x^{3}-222\,\alpha^{2}x^{3}+40\,\alpha
x^{4}+2\,\Omega^{2}\alpha^{2}+242\,\Omega\alpha^{3}-920\,\alpha^{4}-6\,\Omega^{2}\alpha
x+790\,\Omega\alpha^{2}x-2684\,\alpha^{3}x+884\,\Omega\alpha
x^{2}-2596\,\alpha^{2}x^{2}+320\,\Omega x^{3}-840\,\alpha
x^{3}-1200\,\alpha^{3}-3600\,\alpha^{2}x-3600\,\alpha
x^{2}-1200\,x^{3}\Big{)}.$

## Acknowledgements

The research of M. M. Castro was partially supported by PID2021-124332NB-C21
(FEDER(EU) / Ministerio de Ciencia e Innovación-Agencia Estatal de
Investigación) and FQM-262 (Junta de Andalucía). The research of I. Zurrián
was supported by CONICET PIP grant 112-200801-01533, SECyT grant
30720150100255CB, Ministerio de Ciencia e Innovación PID2021-124332NB-C21 and
Universidad de Sevilla VI PPIT - US.

## References

  * [1] Casper, W. R., Yakimov, M.,  Integral operators, bispectrality and growth of Fourier algebras, J. Reine Angew. Math. 766 (2020), 151–194.
  * [2] Casper, W. R., Yakimov, M., The matrix Bochner problem, Amer. J. Math. 144 (2022), no. 4, 1009–1065.
  * [3] Casper, W. R., Grünbaum F. A., Yakimov, M., Zurrian, I., Reflective prolate-spheroidal operators and the KP/KdV equations. Proc. Nat. Academy of Sciences, USA 116 (2019) , no. 37, 1831–18315.
  * [4] Casper, W., Grünbaum, F. A., Yakimov, M., Zurrian, I., Reflective prolate-spheroidal operators and the adelic Grassmannian. Comm. Pure and Applied Math (2023),1–42. DOI: 10.1002/cpa.22118
  * [5] Casper, W., Grünbaum, F. A., Yakimov, M., Zurrian, I.,  Algebras of commuting differential operators for kernels of Airy type. Toeplitz operators and Random Matrices, in memory of Harold Widom. Birkhouser series: Operator Theory, Advances and Applications. E. Basor, A. Böttcher, T. Ehrhardt, C. Tracy editors, 2023.
  * [6] Casper, W., Grünbaum, F. A., Yakimov, M., Zurrian, I., Matrix valued discrete-continuous functions with the prolate spheroidal property and bispectrality (2023), arXiv:2302.05750.
  * [7] Castro, M., Grünbaum, F. A., The Darboux process and time-and-band limiting for matrix orthogonal polynomials. Linear Algebra Appl. 487 (2015), 328-341.
  * [8] Castro, M., Grünbaum, F. A., Time and band limiting for matrix orthogonal polynomials of Jacobi type . Random Matrices: Theory and Appl. 6, no. 4 (2017): 1740001–12.
  * [9] Castro, M. ,Grünbaum, F.A.,Pacharoni, I., Zurrián, I., A further look at time-and-band limiting for matrix orthogonal polynomials, “Frontiers in Orthogonal Polynomials and q-Series”, 139–153. Contemp. Math. Appl. Monogr. Expo. Lect. Notes, 1, World Scientific 2018.
  * [10] Castro, M. M., Grünbaum, F. A.,  A new property of exceptional orthogonal polynomials (2023), arXiv:2210.13928.
  * [11] Connes, A., Moscovici, H., The UV prolate spectrum matches the zeros of zeta. PNAS (2022) 119 (22).
  * [12] Duistermaat J. J., Grünbaum F. A., Differential equations in the spectral parameter. Comm. Math. Phys. 103 (1986), 177–240.
  * [13] Durán, A., Higher order recurrence relation for Exceptional Charlier, Meixner, Hermite and Laguerre orthogonal polynomials. Integral Transforms Spec. Funct. 26 (2015), no. 5, 357–376.
  * [14] Durán, A. Exceptional orthogonal polynomials via Krall discrete polynomials. Lectures on orthogonal polynomials and special functions, 1-75, London Math. Soc. Lecture Note Ser., 464, Cambridge U. Press, Cambridge 2021.
  * [15] Garcia Ferrero, M.A., Gomez Ullate, D., Milson, R. A Bochner type classification theorem for exceptional orthogonal polynomials. J. Math. Anal. Appl. 472 (2019) 584-626.
  * [16] Gomez Ullate, D., Grandati, Y., N., Milson, R., Rational extensions of the quantum harmonic oscillator and exceptional Hermite polynomials. J. Phys. A 47 (2014), no. 1, 015203, 27 pp.
  * [17] Gomez Ullate, D., Kamran, N., Milson, R., An extended class of orthogonal polynomials defined by a Sturm-Liouville problem. J. of Math. Analysis and Appl. 359 (2009) 352-367.
  * [18] Gomez-Ullate, D., Kasman, A., Kuijlaars, A.B.J., Milson, R., Recurrence relations for exceptional Hermite polynomials. J. of Approx. Theory 204 2016, 1-16.
  * [19] Grünbaum, F. A., A new property of reproducing kernels of classical orthogonal polynomials. J. Math. Anal. Applic. 95 (1983), 491–500.
  * [20] Grünbaum, F. A., Pacharoni I., Zurrian, I., Time and band limiting for matrix valued functions, an example. SIGMA 11 (2015), 044, 14 pages.
  * [21] Grünbaum, F. A., Pacharoni I., Zurrian, I., Time and band limiting for matrix valued functions: an integral and a commuting differential operator. Inverse Problems 33, No. 2 (2017), 025005.
  * [22] Grünbaum, F. A., Pacharoni I., Zurrian, I., Bispectrality and time and band limiting: matrix valued polynomials. Int. Math. Res. Not. (2020), no. 13, 4016–4036.
  * [23] Grünbaum, F.A. , Vinet, L., Zhedanov, A., Algebraic Heun operator and Band and Time limiting. Com. Math. Physics 364 (2018), 1041–1068.
  * [24] Kasman, A., Milson R., The adelic Grassmannian and exceptional Hermite polynomials, Math Physics Anal. Geometry 23 40 (2020), 51 pp.
  * [25] Katsnelson V., Selfadjoint boundary conditions for the prolate spheroidal differential operator, arXiv:1603.07542, (2016).
  * [26] Landau, H. J., Pollak, H. O., Prolate spheroidal wave functions, Fourier Analysis and Uncertainty, II. Bell System Tech. Journal, Vol. 40, No. 1 (1961), 65–84.
  * [27] Landau, H. J., Pollak, H. O., Prolate spheroical wave functions, Fourier Analysis and Uncertainty, III. Bell System Tech. Journal, Vol. 41, No. 4 (1962), 1295–1336.
  * [28] Osipov, A. Rokhlin, V. Xiao H.,  Prolate spheroidal wave functions of order zero.  Springer Ser. Appl. Math. Sci, 187, 2013.
  * [29] Perline, R.K., Discrete Time-Band Limiting Operators and Commuting Tridiagonal Matrices. SIAM. J. on Algebraic and Discrete Methods, 8(2) (1987), 192–195.
  * [30] Perlstadt, M., Chopped orthogonal polynomial expansions – some discrete cases. SIAM J. Disc. Meth. 4 (1983), 94–100.
  * [31] Perlstadt, M., A property of orthogonal polynomial families with polynomial duals. SIAM J. Math. Anal. 15 (1984), 1043–1054.
  * [32] Quesne, C., Solvable rational potentials and exceptional orthogonal polynomials in supersymmetric quantum mechanics. SIGMA 5 (2009), 24 pp.
  * [33] Osipov, A. Rokhlin, V. Xiao H.,  Prolate spheroidal wave functions of order zero.  Springer Ser. Appl. Math. Sci, 187, 2013.
  * [34] Sasaki,R., Tsujimoto, S., Zhedanov, A., Exceptional Laguerre and Jacobi polynomials and the corresponding potentials through Darboux-Crum transformations. J. Phys. A. Math. 43 (2010), no. 31, 315204, 20 pp.
  * [35] Shannon, C., A mathematical theory of communication. Bell Tech. J., vol 27, 1948, 379–423 (July) and 623–656 (Oct).
  * [36] Slepian, D., Prolate spheroidal wave functions, Fourier Analysis and Uncertainty, IV. Bell System Tech. Journal, Vol. 43, No. 6 (1964), 3009–3058.
  * [37] Slepian, D., On bandwidth. Proc. of IEEE, Vol. 64, no. 3, March 1976.
  * [38] Slepian, D., Prolate spheroidal wave functions, Fourier Analysis and Uncertainty, V. Bell System Tech. Journal, Vol. 57, No. 5 (1978), 1371–1430.
  * [39] Slepian, D., Some comments on Fourier analysis, Uncertainty and Modeling SIAM Review, Vol. 25, 3, July 1983, 379–393.
  * [40] Slepian, D., Pollak H. O., Prolate spheroidal wave functions, Fourier Analysis and Uncertainty, I. Bell System Tech. Journal, Vol. 40, No. 1 (1961), 43–64.


# Ultracool Subdwarfs: The Halo Population Down to the Substellar Limit

Adam J. Burgasser  Sebastién Lépine  Nicolas Lodieu  Ralf-Dieter Scholz
Phillippe Delorme  Wei-Chun Jao  Brandon J. Swift  Michael C. Cushing

###### Abstract

Ultracool subdwarfs are low luminosity, late-type M and L dwarfs that exhibit
spectroscopic indications of subsolar metallicity and halo kinematics. Their
recent discovery and ongoing investigation have led to new insights into the
role of metallicity in the opacity structure, chemistry (e.g. dust formation)
and evolution of low-temperature atmospheres; the long-term evolution of
magnetic activity and angular momentum amongst the lowest-mass stars; the form
of the halo luminosity and mass functions down to the hydrogen-burning mass
limit; and even fundamental issues such as spectral classification and
absolute brightness scales. This Splinter Session was devoted to bringing
advances in observational and theoretical ultracool subdwarf research to the
attention of the low-mass stellar and brown dwarf communities, as well as to
share results among ultracool subdwarf enthusiasts.

###### Keywords:

stars:subdwarfs – stars: abundances – techniques: spectroscopic

###### :

97.10.Ex,97.10.Ri,
97.10.Tk,97.10.Vm,97.10.Wn,97.20.Tr,97.20.Vs,97.80.Di,98.35.Gi

##  1 Introduction

The late-type M and L dwarf classes of ultracool dwarf stars and brown dwarfs
are now well-sampled, due in large part to wide-field red and infrared imaging
surveys such as 2MASS, DENIS and SDSS. We are now beginning to uncover the
metal-poor, halo counterparts to these disk dwarfs, the so-called ultracool
subdwarfs 2005ESASP.560..237B . These sources have been found primarily in
wide-field imaging surveys and the first generations of red/near-infrared
proper motion surveys. Discoveries include the first L-type—and possibly
T-type—subdwarfs, extending our knowledge of the Galactic halo down to the
substellar (brown dwarf) regime. As the low-mass ultracool subdwarfs, with
their extremely long nuclear burning lifetimes, were presumably formed early
in the Galaxy’s history, they are important tracers of Galactic structure and
chemical enrichment history. In addition, detailed studies of their complex
spectral energy distributions are facilitating new insights on the role of
metallicity in the opacity structure, chemistry and evolution of cool
atmospheres; the long-term evolution of magnetic activity and angular
momentum; and fundamental issues of spectral classification and
temperature/luminosity scales.

This Splinter was devoted to highlighting advances in ultracool subdwarf
research, both observational and theoretical. It was organized into three
subtopics: Discoveries, Classification and Fundamental Parameters. Individual
presentations will be made available at http://www.browndwarfs.org/cs15.

##  2 Discoveries

###  2.1 An Extremely Wide and Very Low-Mass Common Proper Motion Pair -
Representatives of a Nearby Halo Stream?
(R.-D. Scholz et al.)

A pair of faint stars sharing exactly the same very large proper motion of
about 860 mas/yr and separated by about six degrees has been discovered in a
high proper motion survey of the southern sky using multi-epoch positions and
photometry from the SuperCOSMOS Sky Surveys. The two stars, SSSPM J2003$-$4433
and SSSPM J1930$-$4311, have been classified as a late-type (M7) dwarf and an
ultracool subdwarf (sdM7) (Figure 1) with individually estimated spectroscopic
distances of 38 pc and 72 pc, respectively. In view of the accurate agreement
in their large proper motions a common distance of about 50 pc and a projected
physical separation of about 5 pc has been assumed, ruling out a physical
binary. The mean heliocentric space velocity of the pair ($U,V,W$) =
($-232,-170,+74$) km s-1 is typical of the Galactic halo population. These
values rely on a preliminary radial velocity measurement and on the assumption
of a common distance and velocity vector. The large separation and the
different metallicities of dwarf and subdwarf make a common formation scenario
as a wide binary (later disrupted) improbable. It seems more likely that this
wide pair is part of an old halo stream 2008A&A…487..595S .

Figure 1: Low-resolution classification spectra of SSSPM J2003$-$4433 (M7) and
SSSPM J1930$-$4311 (sdM7) and some comparison sources (from 2008A&A…487..595S
).

###  2.2 A Metal-poor Mid-T dwarf from the CFBDS Survey
(P. Delorme et al.)

We report the discovery of CFBDS J150000-182407, a T subdwarf candidate which
we identified during the Canada France Brown Dwarf Survey (CFBDS). CFDBS
2008A&A…484..469D  is an $i^{\prime}-z^{\prime}$ wide-field search for
ultracool brown dwarfs which uses the MegaCam wide-field camera on the Canada-
France-Hawaii Telescope (CFHT).

CFBDS1500 is a peculiar T4.5 dwarf, with strong spectroscopic evidence for a
subsolar metallicity. Comparison of the overall shape of its spectrum with
synthetic spectra from 2006ApJ…640.1063B  suggests [M/H]$\sim$-0.3, while the
complete absence of the $\sim$1.25$\mu$m K I doublet (Figure 2) rather
suggests that [M/H]$\leq$-0.5 when comparing with the same models. The
kinematics of CFBDS1500 imply an 80% probability that it belongs to the thick
disk, but leaves $\sim$10% probability that it instead is an older member of
the thin disk or a member of the halo. [M/H]$\sim$-0.5 is consistent with
either thick disc membership or an older thin disk population, while
[M/H]$\sim$-0.3 would lean slightly towards the older thin disk but remains
easily compatible with thick disk membership.

Figure 2: Spectra of T subdwarf candidate CFBDS1500 (black line) and 2M2254
(red line, from 2008ApJ…678.1372C ), the latter a solar-metallicity T dwarf
template that best matches the spectrum of CFBDS1500. Left panel shows the
whole NIR spectrum while the right panel zooms in on $J$-band region. The
vertical lines mark the locations of the K I doublet.

Spectroscopy and kinematics together make CFBDS1500 a strong thick disk
candidate, and as such the first incontrovertibly substellar object which
would not belong to the thin disk. As one of the few low-metallicity brown
dwarfs known, it is a major benchmark for synthetic spectra. The inconsistency
between the metallicity determined from the strength of the K I doublet and
from the general shape of the spectrum will need to be explained. We tend to
trust the metallicity constrain from the K I doublet more, since it depends
(besides, admittedly, the structure of the atmosphere) on a single line pair
with very well understood atomic physics, rather than on millions of
incompletely characterized molecular lines.

If the [M/H]$\sim$-0.5 metallicity holds, the observed spectrum would suggest
that the models overestimate the sensitivity of the general shape of the
spectrum to metallicity, especially in the $J$-band. That would hardly be a
surprise, given the many poorly constrained physical inputs that need to be
taken into account. Since the discovery of brown dwarfs belonging to old
Galactic populations is now possible, it becomes crucial to interpret their
spectra with state-of-the-art evolutionary models. The only one for low
metallicity brown dwarfs (to our knoweldge) was calculated without taking into
account clouds in the atmosphere 2001RvMP…73..719B . New theoretical work on
this topic (e.g., 2008arXiv0808.2611S ) is now a critical point of T sub-
dwarfs characterization.

##  3 Classification

###  3.1 A Spectral Sequence of K- and M-type Subdwarfs
(W.-C. Jao et al.)

Using new spectra of 88 K and M-type subdwarfs, we have considered novel
methods for assigning spectral types and take steps toward developing a
comprehensive spectral sequence for subdwarf types K3.0 to M6.0. The types are
assigned based on the overall morphology of spectra covering 6000Å to 9000Å
through the understanding of GAIA model grids. The types and sequence
presented link the spectral types of cool subdwarfs to their main sequence
counterparts, with emphasis on the relatively opacity-free region from
8200–9000Å. When available, supporting abundance, kinematic, and trigonometric
parallax information is used to provide more complete portraits of the
observed subdwarfs. We find that the CaHn (n= 1–3) and TiO5 indices often used
for subdwarf spectral typing are affected in complicated ways by combinations
of subdwarfs’ temperatures, metallicities, and gravities, and we use model
grids to evaluate the trends in all three parameters. Because of the complex
interplay of these three characteristics, it is not possible to identify a
star as an “extreme” subdwarf simply based on very low metallicity, and we
suggest that the modifiers “extreme” or “ultra” only outline locations on
spectroscopic indices plots, and should not be used to imply low or very low
metallicity stars. In addition, we propose that “VI” be used to identify a
star as a subdwarf, rather than the confusing “sd” prefix, which is also used
for hot O and B subdwarfs that are unrelated to the cool subdwarfs. These
results have been published in 2008AJ….136..840J .

###  3.2 Spectroscopic Sequences of Cool and Ultra-cool Subdwarfs
(sdM/esdM/usdM) from the Sloan Digital Sky Survey
(S. Lépine et al.)

We present new spectra of cool and ultra-cool subdwarfs, obtained from the
Sloan Digital Sky Survey (SDSS). Thousands of these sources have now been
observed in SDSS (see 2008ApJ…681L..33L  and poster by Lépine, these
proceedings). The stars were categorized in a new, recently expanded
classification system 2007ApJ…669.1235L  as subdwarfs (sdM), extreme subdwarfs
(esdM), and ultrasubdwarfs (usdM), based on the ratio of the TiO to CaH
molecular bands. It is argued that the M/sdM/esdM/usdM subclasses represent a
metallicity sequence (cf. 1997PASP..109.1233G ; 2006PASP..118..218W ), with
the usdM being the most metal-poor of the subdwarfs.

While it has been suggested by others that the morphology of the subdwarfs
depends on both their metallicity and gravity, observational evidence supports
the idea that metallicity is the dominant parameter, and that the current
classification system effectively ranks the stars according to mass (subtype)
and metallicity (subclass). As demonstrated in Figures 2 and 3 of
2007ApJ…669.1235L , stars kinematically selected from the Galactic disk
display only small variations in their TiO/CaH ratio, while stars
kinematically selected from the Galactic halo display a very broad range of
TiO/CaH values. The disk population has a broader range of stellar ages, and
should also have a broader range of surface gravities than the halo stars,
which are uniformly old. The fact that the disk stars have a much smaller
range of TiO/CaH than the halo stars indicates than metallicity effects must
dominate in the subdwarfs. Furthermore, atmospheric models indicate that
gravity variations not only change the TiO/CaH ratio, but they also notably
affect the strength of atomic spectral lines such as KI and NaI. However SDSS
spectra show no clear variation in the equivalent widths of the KI, NaI, and
CaII lines for stars of a given subclass and subtype, which again suggest that
stars of similar subtypes have similar surface gravities. Gravity, in the
nearby halo subdwarfs, is not a free parameter but is essentially set by the
initial mass and metallicity of the subdwarf.

Ultra-cool subdwarfs are particularly useful for refining the current
classification scheme, because their molecular bands and atomic lines are more
prominent. Spectra of many cool and ultra-cool subdwarfs are now available in
the SDSS database, and new ones should become available in the future as halo
subdwarfs are now specifically targeted for follow-up spectroscopy. The new
SDSS spectra also have a broad spectral coverage (4000Å-9000Å), and spectra
from many late-type subdwarfs can be combined to obtain high signal-to-noise
spectral templates. A fit to these templates yield a more reliable
classification than the narrowly defined TiO5, CaH2, and CaH3 spectral
indices. The use of the spectral indices in the classification is useful, but
should be phased out in favor of the new SDSS classification templates
(Lépine, in preparation).

###  3.3 L Subdwarfs: Classification, Distance Scale and Low-Metallicity
Condensate Formation
(A. Burgasser et al.)

L subdwarfs are the metal-poor counterparts to the L dwarf class of low mass
stars and brown dwarfs, effectively extending our sampling of halo stars to
the (metallicity-dependent) hydrogen-burning mass limit. With strong metal
hydride bands, weak metal oxides and red optical spectral energy
distributions, the L subdwarfs share many of the same spectral characteristics
as L dwarfs, but differ in their considerably bluer near-infrared colors
($J-K_{s}$ $\approx$ 0 vs 1.5–2.5) and much stronger metal hydride, metal
oxide and alkali line absorption in the red optical (Figure 3). The first L
subdwarfs were identified in 2003; today there are at least three reported in
the literature 2003ApJ…592.1186B ; 2004ApJ…614L..73B ; siv04 , not including
the unusual source LSR 1610-0040 2003ApJ…591L..49L  discussed in detail below.

One of the recent advances in L subdwarf research is a formalized spectral
classification scheme for these objects. While the classification of M
subdwarfs is actively debated due to the large samples now available (see
contributions from S. Lépine and W.-C. Jao), there are too few L subdwarfs to
accurately define the class. 2007ApJ…657..494B  have proposed that L subdwarfs
be classified according to the closest match to the L dwarf standards of
1999ApJ…519..802K  in the 7300–9000 Å range (see Figure 3), a region in which
peculiarities are minimized. This had provided preliminary types for 2MASS
0532+82 (sdL7 2003ApJ…592.1186B ), 2MASS 1626+39 (sdL4 2004ApJ…614L..73B ) and
SDSS 1256-02 (sdL3.5 siv04 ).

Figure 3: Left: Comparison of the red optical spectra of the sdL4 2MASS
1626+39 to the L4 spectral standard 2MASS 1155+23, illustrating their
similarities in the 7300–9000 Å range (from 2007ApJ…657..494B ). Right:
$M_{K_{s}}$/spectral type relation for late-type M and L dwarfs (black points,
polynomial fit indicated by the solid line with dashed lines indicating
$\pm$1$\sigma$) compared to measurements for the sdM7 LHS 377, sdM8 LSR
1425+71 and sdL7 2MASS 0532+82 (red points; 1992AJ….103..638M ;
2008arXiv0806.2336D ; 2008ApJ…672.1159B ). The predicted $M_{K_{s}}$ for the
sdL3.5 SDSS 1256-02 is indicated by the hatched region (from bur1256 ).

As the first L subdwarf identified, 2MASS 0532+82 has been the best studied.
It is the first L subdwarf to have its astrometric parallax measured ($\pi$ =
37.5$\pm$1.7 mas 2008ApJ…672.1159B ), and combined with photometry from 2MASS
and Spitzer 2006ApJ…651..502P  we have measured a luminosity of
$\log_{10}{L_{bol}/L_{\odot}}$ = $-4.24{\pm}0.06$ and estimate Teff =
1730$\pm$90 K for this source. These parameters are comparable to those of
mid-type L field dwarfs, suggesting that L subdwarf classifications are
“later” than corresponding L dwarf types. Interesting, these values make the
hydrogen burning status of 2MASS 0532+82 somewhat ambiguous, since depending
on its metallicity (which is largely unconstrained) this source may be a star,
a brown dwarf, or a brown dwarf which will ultimately reach the main sequence
and become a star! Examining absolute magnitude relations, it appears that the
$M_{K}$/spectral type relation of L subdwarfs is roughly coincident with those
of L dwarfs (Figure 3), although more astrometric work is needed (see also
2008arXiv0806.2336D ). 2MASS 0532+82 is also the only L subdwarf with high
resolution spectral observations 2006AJ….131.1806R , which indicate rapid
rotation ($v\sin{i}=65{\pm}15$ km s-1) suggesting the absence of angular
momentum evolution for very low mass stars. The $UVW$ velocities of this
source unambiguously confirm it as a halo object, with a retrograde Galactic
orbit relative to the Galactic disk.

We are currently examining theoretical spectral models for the sdL3.5 SDSS
1256-02 (bur1256 ; also see next contribution by B. Swift). Using the Phoenix-
Drift models 2008ApJ…675L.105H , we have attempted to reproduce the observed
optical and near-infrared colors and spectra of this source. Our results so
far indicate a Teff $\approx$ 2100–2500 K and [M/H] $\approx$ -1.5 – -1.0 dex,
although there remain strong discrepancies particularly in the red optical.
Interestingly, we find that even with the most advanced treatment of
(metallicity-dependent) cloud formation, clouds in the models are too thick
and condensation too efficient to match the observations of SDSS 1256-02,
consistent with the “cloud suppression” in ultracool subdwarf atmospheres
suggested in prior studies 2003ApJ…592.1186B ; 2006AJ….131.1806R ;
2006AJ….132.2372G ; 2007ApJ…657..494B

##  4 Fundamental Parameters

###  4.1 Theoretical Modeling of L Subdwarf Spectra
(B. Swift et al.)

We present fits of atmospheric models of varying metallicity to the published
optical and near-infrared spectrum and IRAC photometry of the sdL7 2MASS
J05325346+8246465 2003ApJ…592.1186B . This source is the best-observed member
of a growing population of what are suspected to be metal-poor ultracool
dwarfs, and the only one with a parallax distance. The model selection was
made using a goodness-of-fit statistic in conjunction with Monte Carlo
simulated data to account for uncertainties in the absolute flux; this
statistic was then examined in various bands to study systematic issues in the
models and to produce a concordance fit for the properties of this object in
Teff, $\log{g}$ and [M/H] space.

###  4.2 Astrometric Observations of LSR 1610-0040
(M. Cushing et al.)

Since its discovery 2003ApJ…591L..49L , LSR 1610-0040 has defied explanation.
Its red optical spectrum and high proper motion suggest that LSR 1610-0040 is
an early-type L subdwarf, yet its near-infrared spectrum indicates a mid-type
dM or sdM, albeit with numerous peculiar spectral features 2006AJ….131.1797C ;
2006AJ….131.1806R . Based on a comparison of its near-infrared spectrum to
that of field M dwarfs, 2006AJ….131.1797C  assign LSR 1610-0040 a spectral
type of sd/dM6.

We present new astrometric observations of LSR 1610-0040 2008arXiv0806.2336D .
At a measured distance of $\sim$32 pc, its position in color/ absolute
magnitude diagrams ($M_{V}$ vs $V-I$, $M_{K_{s}}$ vs $I-K_{s}$, $M_{K_{s}}$ vs
$J-K_{s}$) is consistent with mid-type M dwarfs. However, its $B-V$ color is
1.2 mags redder than both M subdwarfs and dwarfs of similar spectral types. We
speculate that the $B$-band magnitude of LSR 1610-0040 is suppressed due to
enhanced AlH absorption since Al appears over-abundant relative to solar based
on the strength of the 1.313 $\mu$m Al I doublet.

Perhaps most interesting is that the astrometric observations indicate that
LSR 1610-0040 is an unresolved binary with a period of 1.66 yr. The
photocentric orbit has a moderate eccentricity of 0.44$\pm$0.02, a semi-major
axis of 0.28$\pm$0.01 AU, and an inclination of 83$\pm$1o. Under the
assumption that the secondary contributes little to no light to the system,
and using a near-infrared mass luminosity relation, we estimate the masses of
the components to be MA = 0.095 M⊙ and MB = 0.059 to 0.082 M⊙. We speculate
that LSR 1610-0040A was originally a 0.05 Msun star with [Fe/H] $\sim$ $-2$
which later accreted $<$0.05 M⊙ of material from a massive AGB star that has
undergone hot bottom burning. Pollution by such material, enhanced in Al and
Na and depleted in O, would then explain the peculiar spectrum of LSR
1610-0040A. LSR 1610-0040B has too little mass to be the remnant white dwarf
of such a hypothetical AGB star so the AGB star in question must have been a
more distant companion that has since been lost from the system.

We thank the conference organizers for providing a forum to discuss this
topic.

## References

  * (1) A. J. Burgasser, J. D. Kirkpatrick, and S. Lépine, “Ultracool subdwarfs: metal-poor stars and brown dwarfs extending into the late-type M, L and T dwarf regimes,” in _13th Cambridge Workshop on Cool Stars, Stellar Systems and the Sun_ , edited by F. Favata, G. A. J. Hussain, and B. Battrick, 2005, vol. 560 of _ESA Special Publication_ , pp. 237.
  * (2) R.-D. Scholz, N. V. Kharchenko, N. Lodieu, and M. J. McCaughrean, _A &A_ 487, 595–599 (2008).
  * (3) P. Delorme, C. J. Willott, T. Forveille, X. Delfosse, C. Reylé, E. Bertin, L. Albert, E. Artigau, A. C. Robin, F. Allard, R. Doyon, and G. J. Hill, _A &A_ 484, 469–478 (2008).
  * (4) A. Burrows, D. Sudarsky, and I. Hubeny, _ApJ_ 640, 1063–1077 (2006).
  * (5) M. C. Cushing, M. S. Marley, D. Saumon, B. C. Kelly, W. D. Vacca, J. T. Rayner, R. S. Freedman, K. Lodders, and T. L. Roellig, _ApJ_ 678, 1372–1395 (2008).
  * (6) A. Burrows, W. B. Hubbard, J. I. Lunine, and J. Liebert, _Reviews of Modern Physics_ 73, 719–765 (2001).
  * (7) D. Saumon, and M. S. Marley, _ArXiv e-prints_ 808 (2008).
  * (8) W.-C. Jao, T. J. Henry, T. D. Beaulieu, and J. P. Subasavage, _AJ_ 136, 840–880 (2008).
  * (9) S. Lépine, and R.-D. Scholz, _ApJL_ 681, L33–L36 (2008).
  * (10) S. Lépine, R. M. Rich, and M. M. Shara, _ApJ_ 669, 1235–1247 (2007).
  * (11) J. Gizis, and I. Reid, _PASP_ 109, 1232-1236 (1997).
  * (12) V. M. Woolf, and G. Wallerstein, _PASP_ 118, 218–226 (2006).
  * (13) A. J. Burgasser, J. D. Kirkpatrick, A. Burrows, J. Liebert, I. N. Reid, J. E. Gizis, M. R. McGovern, L. Prato, and I. S. McLean, _ApJ_ 592, 1186–1192 (2003).
  * (14) A. J. Burgasser, _ApJL_ 614, L73–L76 (2004).
  * (15) T. Sivarani, A. K. Kembhavi, and J. Gupchup, _ArXiv e-prints_ (2004).
  * (16) S. Lépine, R. M. Rich, and M. M. Shara, _ApJL_ 591, L49–L52 (2003).
  * (17) A. J. Burgasser, K. L. Cruz, and J. D. Kirkpatrick, _ApJ_ 657, 494–510 (2007).
  * (18) J. D. Kirkpatrick, I. N. Reid, J. Liebert, R. M. Cutri, B. Nelson, C. A. Beichman, C. C. Dahn, D. G. Monet, J. E. Gizis, and M. F. Skrutskie, _ApJ_ 519, 802–833 (1999).
  * (19) D. G. Monet, C. C. Dahn, F. J. Vrba, H. C. Harris, J. R. Pier, C. B. Luginbuhl, and H. D. Ables, _AJ_ 103, 638–665 (1992).
  * (20) C. C. Dahn, H. C. Harris, S. E. Levine, T. Tilleman, A. K. B. Monet, R. C. Stone, H. H. Guetter, B. Canzian, J. R. Pier, W. I. Hartkopf, J. Liebert, and M. Cushing, _ArXiv e-prints_ 806 (2008).
  * (21) A. J. Burgasser, F. J. Vrba, S. Lépine, J. A. Munn, C. B. Luginbuhl, A. A. Henden, H. H. Guetter, and B. C. Canzian, _ApJ_ 672, 1159–1166 (2008).
  * (22) A. J. Burgasser, S. Witte, C. Helling, and P. H. Hauschildt, _ApJ_ (in preparation).
  * (23) B. M. Patten, J. R. Stauffer, A. Burrows, M. Marengo, J. L. Hora, K. L. Luhman, S. M. Sonnett, T. J. Henry, D. Raghavan, S. T. Megeath, J. Liebert, and G. G. Fazio, _ApJ_ 651, 502–516 (2006).
  * (24) A. Reiners, and G. Basri, _AJ_ 131, 1806–1815 (2006).
  * (25) C. Helling, M. Dehn, P. Woitke, and P. H. Hauschildt, _ApJL_ 675, L105–L108 (2008).
  * (26) J. E. Gizis, and J. Harvin, _AJ_ 132, 2372–2375 (2006).
  * (27) M. C. Cushing, and W. D. Vacca, _AJ_ 131, 1797–1805 (2006).


# Turbulent and adaptive optics corrected point-spread functions as
convolutive orders of the phase power spectral density

Romain JL. Fétick,1,2 Benoit Neichel,1 Laurent M Mugnier,2 Aurélie Montmerle-
Bonnefois,2 and Thierry Fusco1,2
1Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France
2ONERA, The French Aerospace Lab BP72, 29 avenue de la Division Leclerc, 92322
Châtillon Cedex, France
E-mail<EMAIL_ADDRESS>
(Accepted 2018 September 13. Received 2018 September 13; in original form 2018
July 20)

###### Abstract

Ground-based astronomy is severely limited by the atmospheric turbulence,
resulting in a large Point-Spread-Function (PSF) and poor imaging resolution.
Even imaging with Adaptive Optics (AO) cannot completely correct the aberrated
wavefront, and a residual turbulence still corrupts the observation. Thus the
consequences of the turbulence on the PSF is of first interest when building
any ground-based telescope. The Power Spectral Density (PSD) of a spatially
stationary turbulent phase carries all the information needed for describing
the long-exposure PSF. We then develop an analytical description of the long-
exposure PSF as a series expansion of the aberrated phase PSD. Our description
of the PSF given the PSD of the phase is a simple theoretical way to describe
the impact of turbulence on the PSF. We also show accordance with previous
papers when restricting our model to its first expansion order. Finally we
derive applications of our formula to some particular cases, such as
Kolmogorov or von-Kármán models, or the AO correction impact on the PSF.

###### keywords:

instrumentation: adaptive optics – high angular resolution – methods:
analytical

††pubyear: 2018††pagerange: Turbulent and adaptive optics corrected point-
spread functions as convolutive orders of the phase power spectral density–A

##  1 Introduction

The Point Spread Function is a critical information for any optical system
since it is directly related to the imaging resolution. Diffraction limited
optics with a circular aperture results in the famous Airy pattern. However,
in practice, the wavefront of the observed scene suffers from phase
aberrations during its propagation through an inhomogeneous medium, reducing
imaging performance. In case of ground-based telescopes, the atmospheric
turbulence distorts the wavefront and severely limits the resolution of
telescopes (Roddier, 1981). Adaptive Optics (AO) is a technique using
deformable mirrors to reduce phase aberrations by flattening the wavefront
back (Roddier, 1999). Nevertheless the AO correction is not perfect and some
phase aberrations remain, especially the high spatial frequencies which are
not seen by the wavefront sensor and thus not corrected.

On the one hand, Roddier (1981) developed a relation to find the long-exposure
PSF from the phase Power Spectral Density (PSD). The formula is exact but the
link PSD/PSF is not direct since it requires a double Fourier transform. On
the other hand Racine et al. (1999); Bloemhof et al. (2001); Jolissaint &
Veran (2002) derived a first order approximation of the Roddier method. The
advantage is to explain the direct link between PSD and PSF, by canceling the
double Fourier transformation, however the drawback is an inaccuracy due to
the first order approximation.

The two methods above illustrate the very close relation between the phase PSD
and the PSF, each one with its advantages and drawbacks. In this paper we use
the Roddier method to derive an exact expression of the PSF, with a direct
analytical link PSD/PSF without requiring the multiple Fourier transforms.
Indeed, our method relies on a convolutional formalism, where the PSF is
described as an infinite sum of convolutive orders of the phase PSD. Finally,
when restricting our expression to its first order, we find the approximations
of Racine et al.; Bloemhof et al.; Jolissaint & Veran.

Straightforward applications involve practical cases where the turbulent, or
partially corrected AO, phase is not known for each time-step, but its PSD
follows well-known statistics (e.g. Kolmogorov, von-Kármán law, AO residuals
PSD). As a consequence our model can be easily fed with available information
on the turbulence. Our direct expression between PSD and PSF also allows for
simple estimation of atmospheric parameters directly from the PSF.

Section 2 develops our analytical method that writes the PSF as orders of the
residual aberrated phase PSD. This method is general and can be used for any
spatially stationary phase. Section 3 shows some applications related to
Kolmogorov or von-Kármán statistics. We give an explicit dependence of the PSF
shape with the Fried parameter $r_{0}$ (Fried, 1966). We also derive a simple
description of the halo for long-exposure PSF using AO correction. Section 4
concludes our work. We discuss about future applications of our model.

##  2 Method

###  2.1 Factorization of the OTF

Let’s call $P$ the entrance pupil transmission function (that may include
static aberrations), $h$ the long-exposure PSF and $\psi$ the complex phasor
(Goodman, 1968) that is the complex scalar field when observing a point
source. The scintillation is neglected so the complex phasor has unit
amplitude. The long-exposure Optical Transfer Function (OTF) is given by the
autocorrelation of the phasor on the pupil (Goodman)

$\tilde{\mathrm{h}}(\vec{\rho}/\lambda)=\frac{1}{S}\iint
P(\vec{r})P(\vec{r}+\vec{\rho})\langle\psi(\vec{r},t)\psi(\vec{r}+\vec{\rho},t)\rangle_{t}\mathrm{d}\vec{r}$
(1)

where $S$ ensures normalization on the surface of the pupil,
$\langle\cdot\rangle_{t}$ denotes the time average and $\tilde{h}$ stands for
the Fourier transform of the PSF $h$. We then write the phasor $\psi$ as the
complex exponential of the aberrated phase $\phi$

$\psi(\vec{r},t)=e^{i\phi(\vec{r},t)}$ (2)

Roddier (1981) has shown that assuming Gaussian statistics of the phase, the
long-exposure OTF writes

$\tilde{\mathrm{h}}(\vec{\rho}/\lambda)=\frac{1}{S}\iint
P(\vec{r})P(\vec{r}+\vec{\rho})e^{-\frac{1}{2}D_{\phi}(\vec{r},\vec{\rho})}\mathrm{d}\vec{r}$
(3)

where $D_{\phi}(\vec{r},\vec{\rho})$ is the structure function of the residual
phase, defined as

$D_{\phi}(\vec{r},\vec{\rho})=\langle|\phi(\vec{r},t)-\phi(\vec{r}+\vec{\rho},t)|^{2}\rangle_{t}$
(4)

Assuming spatially stationary residual phase, the dependency in the position
$\vec{r}$ vanishes and only remains the separation $\vec{\rho}$. This
stationarity hypothesis is verified for a purely tubulent phase, and is a
reasonable approximation in case of partial AO correction (Conan, 1994). The
approximation is better for larger aperture telescopes, since the non-
stationarity of the residual phase essentially affects the edges of the pupil.
The structure function then simplifies

$D_{\phi}(\vec{r},\vec{\rho})=D_{\phi}(\vec{\rho})=\langle\langle|\phi(\vec{r},t)-\phi(\vec{r}+\vec{\rho},t)|^{2}\rangle_{t}\rangle_{r}$
(5)

The OTF can then be separated into two OTFs, the telescope one and the
atmospheric one

$\tilde{\mathrm{h}}(\vec{\rho}/\lambda)=\tilde{\mathrm{h}}_{T}(\vec{\rho}/\lambda)\cdot\tilde{\mathrm{h}}_{A}(\vec{\rho}/\lambda)$
(6)

with the telescope OTF

$\tilde{\mathrm{h}}_{T}(\vec{\rho}/\lambda)=\frac{1}{S}\iint
P(\vec{r})P(\vec{r}+\vec{\rho})\mathrm{d}\vec{r}$ (7)

and the atmospheric OTF

$\tilde{\mathrm{h}}_{A}(\vec{\rho}/\lambda)=e^{-\frac{1}{2}D_{\phi}(\vec{\rho})}$
(8)

Even though this expression may give the PSF by a numerical Fourier transform,
it doesn’t exhibit an explicit dependence in the PSD of the phase. It is then
necessary to modify the equation to find the direct link we are looking for.

###  2.2 Expression of the PSF

Assuming the phase variance is finite, let’s define the autocorrelation of the
phase

$B_{\phi}(\vec{\rho})=\langle\langle\phi(\vec{r},t)\phi(\vec{r}+\vec{\rho},t)\rangle_{t}\rangle_{r}$
(9)

Directly follows the equality

$D_{\phi}(\vec{\rho})=2[B_{\phi}(\vec{0})-B_{\phi}(\vec{\rho})]$ (10)

Using this identity and the series expansion of the exponential, one gets the
atmospheric OTF

$\tilde{\mathrm{h}}_{A}(\vec{\rho}/\lambda)=e^{-B_{\phi}(\vec{0})}e^{B_{\phi}(\vec{\rho})}=e^{-B_{\phi}(\vec{0})}\sum_{n=0}^{+\infty}\frac{B_{\phi}(\vec{\rho})^{n}}{n!}$
(11)

We denote $\mathcal{F}\\{f\\}(\vec{u})$ the Fourier transform operator applied
on the function $f$, evaluated in $\vec{u}$. Using the linearity of the
Fourier operator on the infinite sum, the PSF writes

$\mathrm{h}_{A}(-\vec{u})=e^{-B_{\phi}(\vec{0})}\sum_{n=0}^{+\infty}\frac{\mathcal{F}\\{B_{\phi}(\vec{\rho})^{n}\\}(\vec{u})}{n!}$
(12)

where $\vec{\rho}/\lambda$ and $\vec{u}$ are conjugated variable by the
Fourier transform. Since the Fourier transform links multiplication and
convolution, we define the convolutive orders of a function $W$ as

$\left\\{\star^{n}~{}W\right\\}(u)=\left\\{\begin{array}[]{cl}\delta(u)&\text{
if }n=0\\\ W(u)&\text{ if }n=1\\\ (W\star W\star...\star W)(\vec{u})&\text{ if
}n\geq 2\\\ \end{array}\right.$ (13)

where $\delta(\vec{u})$ denotes the Dirac distribution. Then the PSF writes

$\mathrm{h}_{A}(-\vec{u})=e^{-B_{\phi}(\vec{0})}\sum_{n=0}^{+\infty}\frac{\star^{n}~{}\mathcal{F}\\{B_{\phi}(\vec{\rho})\\}(\vec{u})}{n!}$
(14)

By the Wiener-Khintchine theorem, the Fourier transform of the autocorrelation
is equal to the PSD of the phase. We write $W_{\phi}$ the PSD of the phase,
and taking care of coordinates dilation for the Fourier transform

$\mathrm{h}_{A}(-\vec{u})=e^{-B_{\phi}(\vec{0})}\sum_{n=0}^{+\infty}\frac{\left\\{\star^{n}~{}W_{\phi}/\lambda^{2}\right\\}(\vec{u}/\lambda)}{n!}$
(15)

This equation states that the atmospheric long-exposure PSF can be written as
an infinite sum of convolutive orders of the phase PSD, under assumption of
spatially stationary phase. In other words, the long-exposure PSF can be
retrieved when the turbulent phase for each time-step is unknown but its PSD
is well-known. This is the case for long-exposure observation through the
turbulent atmosphere, without or with adaptive optics. Then, using Eq. (6) to
retrieve the total PSF, atmospheric and diffraction, directly writes

$h(\vec{u})=h_{T}(\vec{u})\star h_{A}(\vec{u})$ (16)

The convolutive formalism of Eq. 15 can be found in the literature
(Sivaramakrishnan et al., 2002; Perrin et al., 2003), but for short-exposure
PSF only. These authors described the short-exposure PSF as a double infinite
sum on the instantaneous phase. As we show in this paper, based on Roddier’s
method, the long-exposure PSF allows simplifications to a unique infinite sum
on the phase PSD.

##  3 Applications

###  3.1 Factorization of the $r_{0}$ dependency

For usual turbulence spectrum such as Kolmogorov or von-Kármán, there is a
direct dependency of the PSF in the Fried parameter $r_{0}$. Let’s consider
the von-Kármán spectrum

$W_{\phi}(f)=0.023r_{0}^{-5/3}\left(\frac{1}{L_{0}^{2}}+f^{2}\right)^{-11/6}$
(17)

where $L_{0}$ is the external scale of the turbulence. The von-Kármán spectrum
for $L_{0}=8$ m is plotted on Fig. 1. Letting $L_{0}\to\infty$ one would find
the Kolmogorov spectrum.

Figure 1: Three models of phase PSD considered in the following applications:
von-Kármán (red), perfect AO as truncated Kolmogorov (green) and more
realistic AO correction (blue).

Let’s now extract the $r_{0}$ dependency from the PSD. We define $Y_{\phi}$
the PSD of the phase with unit $r_{0}$. Then, by definition of $Y_{\phi}$

$\frac{W_{\phi}(f)}{\lambda^{2}}=r_{0}^{-5/3}Y_{\phi}(f)$ (18)

The convolutive orders of $Y_{\phi}$ can be numerically computed only once for
a given $L_{0}$. The PSF for different values of $r_{0}$ is then a weighted
sum of these components. Plugging $Y_{\phi}$ into Eq. (15) gives the PSF

$\mathrm{h}_{A}(-\vec{u})=e^{-B_{\phi}(\vec{0})}\sum_{n=0}^{+\infty}\frac{r_{0}^{-5n/3}}{n!}\left\\{\star^{n}~{}Y_{\phi}\right\\}(\vec{u}/\lambda)$
(19)

The impact of the strength of the turbulence on the PSF is then more explicit:
for high $r_{0}$ (weak turbulence) the high order convolutive terms have
little importance. Consequently the shape of the PSF is dominated by the low
convolutive orders. For low $r_{0}$ (strong turbulence) the higher convolutive
orders get more importance for describing the PSF shape.

For a practical application, the infinite sum is truncated up to a number
$N_{\mathrm{max}}$ of orders. One would like to get an idea on this number of
terms to consider for having an accurate PSF model. The number
$N_{\mathrm{max}}$ depends on $r_{0}$ but also on the evolution of the
$Y_{\phi}$ convolutive orders with $n$. Let’s call $E$ the integral of
$Y_{\phi}$, or equivalently the phase variance for a unit $r_{0}$. We find a
semi-empirical approximation (see explanations in appendix A) for values of
$\vec{u}$ in the halo

$\left\\{\star^{n}~{}Y_{\phi}\right\\}(\vec{u}/\lambda)\approx
E^{n-1}Y_{\phi}(\vec{u}/\lambda)$ (20)

Dropping the term $n=0$ (particular behavior of the Dirac in
$\vec{u}=\vec{0}$), we get a rough estimate of the PSF behavior towards the
orders

$\forall\vec{u}\in\\{\text{halo}\\}~{},~{}~{}\mathrm{h}_{A}(-\vec{u})\approx
Y_{\phi}(\vec{u}/\lambda)e^{-B_{\phi}(\vec{0})}\sum_{n=1}^{+\infty}\frac{E^{n-1}r_{0}^{-5n/3}}{n!}$
(21)

The terms of the series are decreasing in magnitude when the following
constraint on $n$ is satisfied

$\frac{E^{n}r_{0}^{-5(n+1)/3}}{(n+1)!}<\frac{E^{n-1}r_{0}^{-5n/3}}{n!}\Longleftrightarrow
n>Er_{0}^{-5/3}-1$ (22)

Let’s define $N_{\mathrm{crit}}$ the smaller integer satisfying the
inequality. One cannot truncate the sum before $N_{\mathrm{crit}}$ or would
miss terms of higher magnitude. This basic constraint gives a lower bound on
the number of terms to be considered in the series for accurate description of
the PSF, mathematically we should ensure
$N_{\mathrm{max}}>N_{\mathrm{crit}}$.The criterion can be visualized on Fig. 2
for different values of $r_{0}$ and a chosen $E=0.18$ (corresponding to the
perfect AO correction up to the cutoff frequency $f_{\mathrm{AO}}=0.625$ m-1,
see Sec. 3.3). Once again we see that more terms are needed for high phase
variance $\sigma^{2}=B_{\phi}(0)=Er_{0}^{-5/3}$. In case of AO correction, for
a given $r_{0}$, $E$ depends mainly on the AO cutoff frequency.

Figure 2: $E^{n-1}r_{0}^{-5n/3}/n!$ versus $n$, given for $E=0.18$. The
different black curves are plotted for different values of the Fried parameter
$r_{0}$. All curves pass by a maximum before decreasing. The colored vertical
lines help for visualizing the $N_{\mathrm{crit}}$ criterion for which the
maximum of the curves is reached.

###  3.2 Validation on a von-Kármán turbulence

In order to test the validity of the formula, we consider a von-Kármán
spectrum for the turbulent phase, given Eq. (17). On one hand the PSF is
generated using the atmospheric OTF as exponential of the phase structure
function Eq. (8). On the other hand, we use our expression giving directly the
PSF with the PSD of the phase Eq. (15). Figure 3 shows the PSF for
$N_{\mathrm{max}}=5$ and $N_{\mathrm{max}}=30$ convolutive orders. For a
wavelength of $\lambda=1600$ nm, the convergence requires approximately 7
orders. This is in accordance with our minimal criterion Eq. (22) that gives
$N_{\mathrm{crit}}=4$ at this wavelength. For a shorter wavelength
($\lambda=800$ nm) the PSD is higher and the convergence is much slower and
requires approximately 30 orders. The criterion gives $N_{\mathrm{crit}}=22$,
it conforms the fact that the criterion is only a lower bound, but in practice
many more orders can be required to reach convergence.

Figure 3: PSF computed from a von-Kármán spectrum $L_{0}=8$ m and $r_{0}=12$
cm at $\lambda=400$ nm. Colored lines show PSF for two wavelengths
$\lambda=800$ nm and $\lambda=1600$ nm. Plain: PSF computed by the Roddier OTF
method. Dashed: PSF computed by the convolutive orders method. Left panel: 7
orders. Right panel: 30 orders.

This test shows the validity of our formula with respect to the Roddier method
on the OTF. This example illustrates the difficulty to approximate the PSF,
especially for shorter wavelengths, by the first order of the PSD for a full
von-Kármán spectrum. The method is more suited for low energetic PSD obtained
with AO correction, as developed in the following paragraphs.

###  3.3 Truncated Kolmogorov & Adaptive Optics

Let’s consider a simplified adaptive optics system with a linear number of
actuators $N_{AO}$ in the pupil plane of aperture diameter $D$. All spatial
frequencies of the turbulent phase below the AO cutoff frequency
$f_{AO}=N_{AO}/(2D)$ are assumed to be perfectly corrected by the AO system.
All frequencies above $f_{AO}$ are not seen by the AO system and follow
Kolmogorov statistics. The PSD of the residual phase writes

$W_{\phi}(f)=\left\\{\begin{array}[]{cl}0&\text{ if }f<f_{\mathrm{AO}}\\\
0.023r_{0}^{-5/3}f^{-11/3}&\text{ if }f\geq f_{\mathrm{AO}}\\\
\end{array}\right.$ (23)

This shape of the AO corrected PSD is plotted in green in Fig. 1, for a chosen
AO cutoff frequency $f_{\mathrm{AO}}=0.625$ m-1, corresponding to $N_{AO}=10$
actuators for a $D=8$ m aperture. In practice, the center of the PSF is
dominated by highly energetic low frequency modes, such as residual tip-tilt.
Our simple model of a perfect AO high-pass filter doesn’t take into account
these residues, inducing poor fidelity in the center of the computed PSF.
However we can expect a correct description of the halo. Indeed Fig. 4 shows
good agreements between the two PSF computation methods, using a numerical
Fourier transform Eq. (8) or using the convolutive orders Eq. (15). Let’s
notice that the $r_{0}$ factorization Sec. 3.1 is still applicable for
numerical computation efficiency with our method on convolutive orders of the
PSD. The atmospheric PSF is then convolved by $h_{T}$ to take into account the
telescope contribution Eq. (16). At first order, we retrieve the shape
developed in Racine et al. (1999) with a coherent core made of the Airy
pattern and a turbulent halo depending on the residual PSD. The difference is
that our method includes all the convolutive orders of the PSD in the
description of the PSF. Higher orders can be of high magnitude when the PSD is
of high energy $E$ or for strong turbulence (Fig. 5). In particular the
convolution folds back energy from the PSF halo into the corrected area,
reducing contrast in this zone. Magnitude of higher orders can be reduced (and
then the contrast in the corrected area increased) by decreasing the integral
$E$ of $Y_{\phi}$, which is strongly dependent on the cutoff frequency
$f_{\mathrm{AO}}$ and the wavelength $\lambda$ (Fig. 4). Few orders,
$N_{\mathrm{max}}\simeq 3-8$ in our example, are required to describe the AO
corrected PSF since there is few energy in the PSD, to be compared with the
$N_{\mathrm{max}}>30$ orders required to describe the pure von-Kármán
turbulent phase. Our criterion gives respectively $N_{\mathrm{crit}}=6$ for
$\lambda=400$ nm and $N_{\mathrm{crit}}=1$ for $\lambda=800$ nm. This is in
agreement with the PSFs plotted on the figure since the convergence is reached
at $N_{\mathrm{max}}=8$ for $\lambda=400$ nm, and $N_{\mathrm{max}}=3$ for
$\lambda=800$ nm.

Figure 4: PSF computed for the given AO spectrum. Plain curve: PSF computed by
the OTF method. Dashed: PSF computed by the convolutive orders method. Left
panel: 3 orders. Right panel: 8 orders. Figure 5: Impact of the $W_{\phi}$
convolutive orders in the PSF. Each order is also convolved with the pupil
PSF, called $h_{T}$ in the text. Left: $r_{0}=20$ cm. Right: $r_{0}=10$ cm.
The three orders $n\in\\{0,1,2\\}$ are respectively plotted in colored dashed
lines. The PSF based on the Roddier method, from the OTF in $e^{-D_{\phi}/2}$,
is plotted in plain black.

###  3.4 Modeling errors in the AO system

We previously considered a perfect AO system up to the cutoff frequency, with
a phase PSD identically null for $f<f_{\mathrm{AO}}$. However in practice a
residue of phase PSD remains at low frequencies, due to errors of the AO
system such as WFS noise, corrective-loop delay, or deformable mirror fitting
error. These error terms can be measured or modeled (Fusco et al., 2006) for
finer description of the AO partial correction. Plugging in these errors in
our PSD, we can obtain a better approximation of the PSF, especially in the
central peak. Methods directly based on Roddier’s OTF formula are already
implemented (Rigaut et al., 1998; Jolissaint, 2010; Correia et al., 2017). Our
convolutive formalism then offers a direct understanding of the link between
phase PSD and the PSF on the results provided by Roddier’s OTF based
algorithms. For extreme AO correction we find that the PSF behaves as the
phase PSD (first order approximation) whereas the following PSD convolutive
orders are required for lower AO correction quality.

We choose here to model the residual phase PSD due to AO errors by a Moffat
function (Moffat, 1969). The shape of the residual PSD writes

$W_{\phi}(f)=\left\\{\begin{array}[]{cl}A(1+f^{2}/\alpha^{2})^{-\beta}+B&\text{
if }f<f_{\mathrm{AO}}\\\ 0.023r_{0}^{-5/3}f^{-11/3}&\text{ if
}f>f_{\mathrm{AO}}\\\ \end{array}\right.$ (24)

where $A$, $B$, $\alpha$ and $\beta$ are the Moffat parameters associated to
the error budget on the corrected phase. The shape of the PSD for arbitrary
values of Moffat parameters is plotted in blue on Fig. 1. The resulting PSD is
in-between the von-Kármán uncorrected turbulence, and the perfect AO
correction.

Figure 6: PSF computed for the given AO spectrum. Plain curve: PSF computed by
the OTF method. Dashed: PSF computed by the convolutive orders method. Left
panel: 3 orders. Right panel: 15 orders.

Figure 6 shows the convolutive orders associated to the discussed PSD, for
$N_{\mathrm{max}}=3$ and $N_{\mathrm{max}}=15$. The number of terms required
is higher than in the case of a perfect AO system due to a higher energy under
the PSD, but lower than in the case of uncorrected turbulence.

##  4 Conclusions

For long-exposure observations, the highly time-variable phase aberrations
average, and only remains the statistical behavior of the turbulence. We
consequently derived a general method to compute the long-exposure PSF, given
only the PSD of the turbulent phase. When restricting our model to its first
order approximation, we find agreement with previous works on the subject
(Racine et al., 1999; Bloemhof et al., 2001). We have shown that the higher
convolutive orders have more impact on the PSF for stronger turbulence and
must be taken into account for low Strehl systems.

A hypothesis of spatial stationarity of the phase over the pupil is required,
that is in agreement with Kolmogorov and von-Kármán statistics . Our method is
also applicable to partially corrected phase with adaptive optics. Finally, we
extracted the $r_{0}$ dependency for better understanding of the PSF shape
with the turbulence strength. It is then possible to get a direct estimation
of the $r_{0}$ given the PSF halo, or reciprocally an estimation of the PSF
halo given the $r_{0}$. We plan to use this model for an accurate description
of the halo for AO corrected PSF. The link we make between the PSD and the PSF
can be used as a direct measurement of the Fried parameter $r_{0}$ from a PSF
observed on a star. Finally the strong correlation between the AO corrected
PSD and the shape of the PSF can provide a tool for AO system diagnostics
(Jolissaint & Veran, 2002).

## Acknowledgements

This work was supported by the French Direction Générale de l’Armement (DGA)
and Aix-Marseille Université (AMU).

This project has received funding from the European Union’s Horizon 2020
research and innovation program under grant agreement No 730890. This material
reflects only the authors views and the Commission is not liable for any use
that may be made of the information contained therein.

This study has been partly funded by the French Aerospace Lab (ONERA) in the
frame of the VASCO Research Project.

## References

  * Bloemhof et al. (2001) Bloemhof E., Dekany R., Troy M., Oppenheimer B., 2001, The Astrophysical Journal Letters, 558, L71
  * Conan (1994) Conan J.-M., 1994, PhD thesis, Université Paris XI Orsay
  * Correia et al. (2017) Correia C. M., Bond C. Z., Sauvage J.-F., Fusco T., Conan R., Wizinowich P. L., 2017, JOSA A, 34, 1877
  * Fried (1966) Fried D. L., 1966, JOSA, 56, 1372
  * Fusco et al. (2006) Fusco T., et al., 2006, Optics Express, 14, 7515
  * Goodman (1968) Goodman J. W., 1968, Appendix B Sec. B, 3
  * Jolissaint (2010) Jolissaint L., 2010, Journal of the European Optical Society – Rapid Publications, 5, 10055
  * Jolissaint & Veran (2002) Jolissaint L., Veran J.-P., 2002, ESO Conference and Workshop Proceedings. Vol. 58
  * Moffat (1969) Moffat A., 1969, Astronomy and Astrophysics, 3, 455
  * Perrin et al. (2003) Perrin M. D., Sivaramakrishnan A., Makidon R. B., Oppenheimer B. R., Graham J. R., 2003, The Astrophysical Journal, 596, 702
  * Racine et al. (1999) Racine R., Walker G. A., Nadeau D., Doyon R., Marois C., 1999, Publications of the Astronomical Society of the Pacific, 111, 587
  * Rigaut et al. (1998) Rigaut F. J., Véran J.-P., Lai O., 1998, Proc. SPIE. Vol. 3353, pp 1038–1048
  * Roddier (1981) Roddier F., 1981, in , Vol. 19, Progress in optics. Elsevier, pp 281–376
  * Roddier (1999) Roddier F., 1999, Adaptive optics in astronomy. Cambridge university press
  * Sivaramakrishnan et al. (2002) Sivaramakrishnan A., Lloyd J. P., Hodge P. E., Macintosh B. A., 2002, The Astrophysical Journal Letters, 581, L59

##  Appendix A On the convolutive approximation

Let’s consider a PSD shape verifying the two following hypothesis

  1. 1.

The majority of the PSD energy is in the center
($|\vec{u}|\ll 1$) rather than in the halo ($|\vec{u}|\gg 1$).

  2. 2.

The halo is slowly varying.

These hypothesis are verified in particular for the $-11/3$ power law of the
Kolmogorov spectrum.

Let’s now consider $\vec{u}$ a coordinate in the halo. The PSD auto-
convolution then writes

$Y_{\phi}\star
Y_{\phi}(\vec{u})=\iint_{\mathbb{R}^{2}}Y_{\phi}(\vec{t})Y_{\phi}(\vec{u}-\vec{t})d\vec{t}$
(25)

The $\mathbb{R}^{2}$ domain is then separated into three subdomains

  * •

$|\vec{t}|\ll 1$ where $Y_{\phi}(\vec{t})$ is large due to the first
hypothesis, and $Y_{\phi}(\vec{u}-\vec{t})$ is slowly varying due to the
second one (when choosing $\vec{u}\gg 1$)

  * •

$|\vec{u}-\vec{t}|\ll 1$ where the opposite behavior comes out

  * •

The remaining domain of $\mathbb{R}^{2}$, called $\mathcal{A}$, where both
$Y_{\phi}(\vec{t})$ and $Y_{\phi}(\vec{u}-\vec{t})$ are small due to the first
hypothesis

Figure 7 shows an example of these domains for a given $\vec{u}$ and a PSD
shape as described in Sec. 3.4. The convolutive integral of $Y_{\phi}$ is
separated on the three domains described above. The integral on the first
domain is

$\iint_{|\vec{t}|\ll
1}Y_{\phi}(\vec{t})Y_{\phi}(\vec{u}-\vec{t})d\vec{t}\simeq
Y_{\phi}(\vec{u})\iint_{|\vec{t}|\ll 1}Y_{\phi}(\vec{t})d\vec{t}$ (26)

The approximation we just made rely on the fact that $|\vec{u}|\gg 1$ is in
the halo and $|\vec{t}|\ll 1$ so $\vec{u}-\vec{t}$ is still in the halo. The
slow variations in the halo (second hypothesis) allows to approximate
$Y_{\phi}(\vec{u}-\vec{t})\simeq Y_{\phi}(\vec{u})$. We make the same
reasoning on the second domain, with a coordinate change
$\vec{x}=\vec{u}-\vec{t}$, so we find

$\iint_{|\vec{u}-\vec{t}|\ll
1}Y_{\phi}(\vec{t})Y_{\phi}(\vec{u}-\vec{t})d\vec{t}\simeq
Y_{\phi}(\vec{u})\iint_{|\vec{x}|\ll 1}Y_{\phi}(\vec{x})d\vec{x}$ (27)

The two first integrals are then identical. Regarding the third domain, due to
the first hypothesis, it is a second order terms as the convolution of two
small terms. We simply write it

$\iint_{\mathcal{A}}Y_{\phi}(\vec{t})Y_{\phi}(\vec{u}-\vec{t})d\vec{t}=\eta(\vec{u})E$
(28)

where $\eta(\vec{u})$ is a small valued function, and $E$ denotes the integral
of $Y_{\phi}$ on $\mathbb{R}^{2}$. Summing the partial convolutions on the
three domains, one gets the full convolution

$Y_{\phi}\star Y_{\phi}(\vec{u})\simeq
2Y_{\phi}(\vec{u})\left[\iint_{|\vec{t}|\ll
1}Y_{\phi}(\vec{t})d\vec{t}\right]+\eta(\vec{u})E$ (29)

Finally, using the first hypothesis, the integral above contains the majority
of the energy. Then we write

$\iint_{|\vec{t}|\ll 1}Y_{\phi}(\vec{t})d\vec{t}=\kappa E$ (30)

where $\kappa\in]0,1]$ represents the ratio of energy in the area
$|\vec{t}|\ll 1$. It depends on the shape of $Y_{\phi}$, in practice we find
$\kappa\simeq 0.85$ for our AO truncated Kolmogorov. Using $\kappa$, the
convolution rewrites

$Y_{\phi}\star Y_{\phi}(\vec{u})\simeq 2\kappa
Y_{\phi}(\vec{u})E+\eta(\vec{u})E$ (31)

If we neglect the $\eta(\vec{u})E$ term, one can easily iterate by recurrence
over the convolutive orders to find

$\\{\star^{n}~{}Y_{\phi}\\}(\vec{u})\simeq(2\kappa E)^{n-1}Y_{\phi}(\vec{u})$
(32)

Errors in this equation come from the different approximations we made and the
modification of shape for $\\{\star^{n}~{}Y_{\phi}\\}$, especially in the
center ($\vec{u}\ll 1$). One should remind that this expression is only true
in order of magnitude since the approximations propagate through the
convolutive orders.

Figure 7: Visualization of the three different domains for a given convolutive
shift $\vec{u}$. Colormap show the product
$Y_{\phi}(\vec{t})Y_{\phi}(\vec{u}-\vec{t})$. Intensity plot is in logarithmic
scale.