Some Open Problems in Random Matrix Theory and the Theory of Integrable
Systems. II

Some Open Problems in Random Matrix Theory
and the Theory of Integrable Systems. II††This paper is a contribution to the
Special Issue on Asymptotics and Universality in Random Matrices, Random
Growth Processes, Integrable Systems and Statistical Physics in honor of Percy
Deift and Craig Tracy. The full collection is available at
http://www.emis.de/journals/SIGMA/Deift-Tracy.html

Percy DEIFT

P. Deift

Department of Mathematics, New York University, 251 Mercer Str., New York, NY
10012, USA<EMAIL_ADDRESS>
Received October 11, 2016, in final form March 10, 2017; Published online
March 14, 2017

We describe a list of open problems in random matrix theory and the theory of
integrable systems that was presented at the conference Asymptotics in
Integrable Systems, Random Matrices and Random Processes and Universality,
Centre de Recherches Mathématiques, Montréal, June 7–11, 2015. We also
describe progress that has been made on problems in an earlier list presented
by the author on the occasion of his 60th birthday in 2005 (see [Deift P.,
Contemp. Math., Vol. 458, Amer. Math. Soc., Providence, RI, 2008, 419–430,
arXiv:0712.0849]).

integrable systems; numerical algorithms; random matrices; random particle
systems; Riemann–Hilbert problems

37K15; 15B52; 34M55; 35Q53; 60F05; 60K35; 62H25; 82B44

At my 60th birthday conference in 2005, I was asked to present a list of
unsolved problems [34]. At my 70th birthday conference in 2015, I was asked to
update the list, and also present any other problems that I found interesting
and important. What follows is, more or less, the list of updates and problems
that I presented. As in [34] the list is

> “written in an informal style. Detailed references are readily available on
> the web. In the text, various authors are mentioned by name: This is to aid
> the reader in researching the problem at hand, and is not meant in any way
> to be a detailed, historical account of the development of the field. I ask
> the reader for his/her indulgence on this score. In addition, the list of
> problems is not meant to be complete or definitive: It is simply an
> (unordered) collection of problems that I think are important and
> interesting, and which I would very much like to see solved”.

I will proceed in the following way. For any given problem listed in [34], I
will first discuss the progress that has been made since 2005, and then I will
discuss some aspects of the problems that remain open. The motivation for
these problems and their significance will be discussed only briefly: For more
information, the reader should consult [34]. I will also discuss some open
problems not on the previous list.

######  Problem 1 (KdV with almost periodic initial data).

Consider the Korteweg–de Vries (KdV) equation

$\displaystyle u_{t}+uu_{x}+u_{xxx}=0$ (1)

with initial data

$\displaystyle u(x,t=0)=u_{0}(x),\qquad x\in\mathbb{R}.$

In the 1970’s, McKean and Trubowitz [76] proved the remarkable result that if
the initial data $u_{0}$ is periodic, $u_{0}(x+p)=u_{0}(x)$ for some $p>0$,
then the solution $u(x,t)$ of (1) is almost periodic in time. This result
leads to the following natural conjecture: The same is true if $u_{0}(x)$ is
almost periodic, i.e., if the initial data is almost periodic in space, the
solution evolves almost periodically in time. This is a very challenging
problem with many new and unique features. For example, for $u_{0}(x)$ almost
periodic, one does not know, by standard PDE methods, whether a solution
$u(x,t)$ exists, even for a small time, never mind globally. In spectacular
and seminal work over the last few years, D. Damanik and M. Goldstein [31,
32], joined later by I. Binder and M. Lukic [7], have partially resolved this
conjecture in the affirmative, namely, for small quasi-periodic analytic
initial data with Diophantine frequencies, unique global solutions exist and
are almost periodic in time. The key insight of Damanik et al. was how to
utilize the formal integrability of KdV in an effective analytical manner. Not
only are the arguments in Damanik et al. applicable to other (formally)
integrable systems (for example [8] proves a Toda lattice analog of the KdV
theorem), but some of their arguments are also applicable to non-integrable
systems. It is a natural, open and very challenging conjecture that the result
of Damanik et al. remains true for (suitable) perturbations of KdV, i.e., for
(suitable) perturbations of KdV, the Cauchy problem with (suitable) almost
periodic initial data $u_{0}(x)$, has a unique global solution $u(x,t)$, which
evolves almost periodically in time. After all, the finite-dimensional analog
of this result is precisely the content of classical KAM theory. We note that
such results for perturbations of KdV and the nonlinear Schrödinger equation
(NLS) with periodic initial data $u_{0}(x)$, have already been proved by S.
Kuksin, T. Kappeler and others.

With the discovery of quasi-crystals, one can anticipate, in particular, that
interest in PDE problems on quasi-periodic backgrounds, will surely grow.

######  Problem 2 (universality for random matrix theory (RMT)).

The question of universality for the distribution of the eigenvalues of a
random matrix on the microscopic scale, goes back to the earliest days of RMT
in the 1950s and the work of E. Wigner, F. Dyson and M. Mehta, in particular.
By 2005, universality for unitary invariant ensembles had been established,
both in the bulk and at the edge, for very general analytic matrix
distributions. Universality had also been established for invariant orthogonal
and symplectic ensembles, both in the bulk and at the edge, for distributions
of the form $e^{-\operatorname{tr}V(M)}dM$, where the matrices $M$ are
$n\times n$ and

$\displaystyle
V(x)=\gamma_{k}x^{2k}+\gamma_{k-1}x^{2k-1}+\dots+\gamma_{0},\qquad\gamma_{k}>0$
(2)

is fixed. For such distributions, the underlying equilibrium measure is
supported on a single interval. However, in the important case with
distribution $e^{-n\operatorname{tr}V(M)}dM$ where $V(x)$ is again fixed as in
(2), and where the underlying equilibrium measure for the problem can now be
supported on more than one interval, a challenging technical issue for
orthogonal and symplectic invariant ensembles involving a certain determinant
(see $D_{n}$ in [34]) remained unresolved. This issue was subsequently
resolved in 2010–2011 by T. Kriecherbauer and M. Shcherbina [70] and M.
Shcherbina [92].

The challenge par excellence that remained for RMT as of 2005, was to prove
universality in the bulk for general Wigner ensembles for which the entries
$M_{ij}$ of the matrices are independent random variables, subject only to the
given symmetry of the matrices. For such ensembles, the algebraic
simplifications and powerful analytical techniques (orthogonal polynomials,
Riemann–Hilbert/steepest-descent methods) developed for invariant ensembles,
are no longer applicable. In 1999, A. Soshnikov [99] proved universality at
the edge for Wigner ensembles. Soshnikov’s proof utilized the so-called moment
method introduced by E. Wigner in the 1950’s, together with certain trace
estimates obtained earlier by Ya. Sinai and A. Soshnikov [95]. Apart from the
technical, combinatorial virtuosity required to carry through his proof,
Soshnikov introduced the strategy that foreshadowed all later work on
universality for Wigner ensembles. Rather than computing explicitly the
distribution of the eigenvalues at the top of the spectrum as $n=\dim
M\to\infty$, as one does for invariant ensembles, Soshnikov just showed that
the limiting distribution, whatever it was for any given Wigner ensemble,
depended only on the first two moments of the distribution of the
(independent) entries $M_{ij}$ of the matrices $M$. So in order to see
explicitly what the limiting eigenvalue distribution was, all one had to do
was to match up the distribution of the $M_{ij}$’s with a Gaussian
distribution with the same two moments. For such Gaussian ensembles the
limiting distributions had famously been computed explicitly by C. Tracy and
H. Widom in 1994 [102] in terms of Painlevé functions. In this way
universality at the edge was proved, by inference rather than by direct
computation. This strategy is of course similar in spirit to standard proofs
of the classical central limit theorem for the sum $S_{n}=x_{1}+\dots+x_{n}$
of iid variables $\\{x_{i}\\}$ with distribution $\chi$, where one shows that
the sum of the variables, suitably scaled, converges to a limiting
distribution which only depends on the first two moments of $\chi$.

The first proof of universality for eigenvalues in the bulk of the spectrum
for a class of matrices broader than the class of invariant matrices, was
given by K. Johansson in 2001 [65]. Johansson considered a Wigner ensemble of
random matrices of the form
$H_{t}=\frac{1}{\sqrt{n}}\big{(}H+\sqrt{t}V\big{)}$, $t>0$. Here Hermitian $H$
is taken from a given Wigner ensemble of $n\times n$ Hermitian matrices with
mean zero and variance $\frac{1}{4}$, and $V$ is an $n\times n$ GUE matrix,
also with mean zero and variance $\frac{1}{4}$. For such ensembles, there are
explicit formulae for the eigenvalue correlation functions closely related to
analogous formulae for invariant Hermitian ensembles, and for fixed $t>0$,
Johansson was able to use these explicit formulae to prove universality as
$n\to\infty$ for the eigenvalues of such matrices $H_{t}$ in the bulk by
direct computation. En route to his proof of universality, Johansson observed
that the matrices $H_{t}$, $t\geq 0$, could be thought of in terms of Dyson’s
Brownian motion (DBM) model, where $H_{t}$ is obtained from $H$ by letting the
matrix elements execute independent Brownian motions, subject only to
Hermitian symmetry, for a time $t$. While Johansson did not utilize this
observation directly in his proof, it turned out that the observation played,
and continues to play, a central role, in the analysis of universality
phenomena in random matrix theory. The realization that Dyson Brownian motion
was related to the issue of universality can already be traced to the 1962
paper [49] where Dyson introduced DBM. Instead of Brownian motion, Dyson
considered the process $\tilde{H}_{t}$ obtained by letting the entries of
$\tilde{H}_{t}$ evolve as independent Ornstein–Uhlenbeck processes. The
process $\tilde{H}_{t}$ has an invariant measure, and most fortunately, that
measure is precisely GUE. As $t\to\infty$, or, as it turns out, for fixed $t$
as $n\to\infty$, the distribution for $\tilde{H}_{t}$ converges to the fixed
point, equilibrium distribution, viz., GUE! Putting aside the differences
between the Brownian process and the Ornstein–Uhlenbeck process, Dyson’s
observations provided compelling motivation for Johansson’s result. The
breakthrough came in 2009 when two independent groups, L. Erdős, S. Peche, J.
Ramirez, B. Schlein and H.-T. Yau [54], on the one hand, and T. Tao and V. Vu
[101] on the other, almost simultaneously announced proofs of universality in
the bulk for broad classes of Hermitian Wigner ensembles.

The basic idea in Yau et al. is to use DBM $\tilde{H}_{t}$ as above. As in the
non-Ornstein–Uhlenbeck case considered by Johansson, for fixed $t>0$, the
distribution of the eigenvalues of $\tilde{H}_{t}$ in the bulk exhibit
universal behavior as $n\to\infty$. However, as emphasized by Dyson in his
1962 paper, the time $t\sim 1$ is needed to reach global equilibrium for all
the eigenvalues of $\tilde{H}_{t}$: In order to reach local equilibrium, times
$t\sim\frac{1}{n}$ should be sufficient. On the other hand, for such small
times $t\sim\frac{1}{n}$, or a little longer, one expects the eigenvalue
distribution of $\tilde{H}_{t}$ to be essentially the same as the distribution
of the eigenvalues of $H$. So one seeks a window given by $0<b<a<1$ such that
for times $\leq t^{-b}$, the eigenvalue distribution of $\tilde{H}_{t}$ is
essentially the same as that of $H$, while for times $\geq t^{-a}$, the
eigenvalues of $\tilde{H}_{t}$ are already in local GUE equilibrium. As
$\frac{1}{t^{a}}<\frac{1}{t^{b}}$, one draws the conclusion that the
eigenvalues of $H$ have the same local distribution as GUE and it is this
window that Yau et al. exploited to prove universality in the bulk for
Hermitian Wigner matrices. Over the years this technique, now called the “3
step” method, has been greatly extended and refined by Yau and his many
collaborators, and has now been used to prove universality in the bulk, and
also at the edge, for orthogonal and symplectic ensembles, and more generally
for certain $\beta$-ensembles for any $\beta>1$, certain Erdős–Renyi random
graphs and certain banded random matrices, and also, most recently [1] for
certain random matrices with short range correlations.

The basic idea of Tao–Vu to prove universality for Hermitian Wigner ensembles
was to use Lindeberg’s replacement trick. In 1922, Lindeberg gave a new proof
of the central limit theorem for sums $S_{n}=x_{1}+\dots+x_{n}$ of iid
variables $\\{x_{i}\\}$ with distribution $\chi$, by replacing the $x_{i}$’s
one at a time by iid Gaussian variables $\\{\hat{x}_{i}\\}$ with the same mean
and variance as $\chi$,

$\displaystyle
S_{n}=x_{1}+x_{2}+\dots+x_{n}\to\hat{x}_{1}+x_{2}+\dots+x_{n}\to\hat{x}_{1}+\hat{x}_{2}+x_{3}+\dots+x_{n}$
$\displaystyle\hphantom{S_{n}}{}\to\dots\to
S_{n}=\hat{x}_{1}+\hat{x}_{2}+\dots+\hat{x}_{n}$

with negligible error as $n\to\infty$. By analogy, starting with a Hermitian
Wigner matrix $H=\\{H_{ij}\\}$, Tao and Vu showed that the $H_{ij}$’s with
distribution $\chi$ could be replaced one at a time by iid Gaussian variables
with the same mean and variance as $\chi$, obtaining in the end, with
negligible error in the bulk eigenvalue statistics as $n\to\infty$, a GUE
matrix $\hat{H}$.

An open problem of great, ongoing interest is the question of the behavior of
random $n\times n$ band matrices $H$ with bandwidth $w\leq n$, $H_{ij}=0$ if
$|i-j|\geq w$. The long-standing conjecture of V. Fyodorov and A. Mirlin [59]
from 1991 is that for iid $H_{ij}$, $|i-j|<w$, with $w\gg\sqrt{n}$, the local
statistics of the eigenvalues of $H$ are GUE, but for $w\ll\sqrt{n}$, the
behavior is Poisson. Indeed for $w=n$, $H$ is a GUE matrix, and for $w=2$, the
tridiagonal (Jacobi) case, results proving Poisson statistics go back to the
90’s (see [79]). This conjecture is of interest both in the real symmetric
case, $H=\overline{H}=H^{T}$, and in the complex Hermitian case, $H=H^{\ast}$.
The primary obstacle in addressing the problem is that the powerful methods of
invariant ensembles and the “3-step”method of Yau et al. using DBM, as well as
the replacement method of Tao–Vu, are no longer applicable/effective. New
methods are needed. The only rigorous result indicating the critical
$\sqrt{n}$ threshold is due to T. Shcherbina [93, 94] concerning mixed moments
of the characteristic polynomial for certain Gaussian ensembles. Shcherbina
used so-called super-symmetric methods first introduced in the physics
community, and then pioneered in the mathematical community by T. Spencer, and
which are, unfortunately, applicable only to Gaussian ensembles.

In all approaches to proving universality in the bulk for general Wigner
ensembles $\\{H\\}$, the behavior of the eigenfunctions of the matrices $H$
plays an important role: Localized eigenfunctions are associated with Poisson
statistics and delocalized eigenfunctions are associated with random matrix
statistics. For banded random matrices, localization of the eigenfunctions has
been proved for $w\ll n^{\frac{1}{8}}$ [91] and delocalization of the
eigenfunctions for $w\gg n^{4/5}$ [53]. Recently [20] have found a way to
overcome some of the obstacles to implementing the “3-step” process, and they
have proved bulk universality for random band matrices with independent
entries in the case $w=\frac{1}{p}n$ for some fixed integer $p$. This is where
the matter stands as of 2016 regarding bulk universality . On the other hand,
Tracy–Widom statistics have been proved for extreme eigenvalues near the edge
of the spectrum for $w\gg n^{5/6}$, an essentially optimal result [97].

Perhaps the most challenging problem currently in the spectral theory of
Schrödinger operators $H=-\Delta+V$, $\Delta=\text{Laplacian}$, is to analyze
the spectrum of $H$ when $V$ is a random function. On the lattice
$\mathbb{Z}^{d}$, a standard model (the Anderson model) is to assume that
$\\{V(x)\\}_{x\in\mathbb{Z}^{d}}$ are iid. In seminal work in 1983, J.
Fröhlich and T. Spencer [58] showed that near the edge of the spectrum of $H$,
the eigenfunctions are localized. It is believed that in the bulk of the
spectrum for $d=3$, the eigenfunctions of $H$ are delocalized and the
eigenvalues exhibit random matrix statistics. There is also believed to be a
so-called mobility edge $\mu$ in the spectrum at which the statistics switch
from a Poisson distribution to a random matrix distribution. So far this
conjecture is wide open. It is of interest to construct and analyze models
where such a transition, Poisson $\to$ random matrix statistics, takes place,
and apart from the intrinsic interest in random band matrices, much of the
interest in random band matrices and the putative $\sqrt{n}$ threshold, lies
in the role of random band matrices as models for localization/delocalization
in random Schrödinger theory. An interesting model that exhibits the Poisson-
GUE transition was introduced by O. Bohigas and M. Pato [10] who considered
the gap probability of a random particle system obtained from the Coulomb gas
by dropping at random a fraction $1-\gamma$ of levels, a procedure known as
“thinning”. As $\gamma\downarrow 0$ the gap probability becomes Poissonian,
and as $\gamma\uparrow 1$ the gap probability follows GUE (see also [18] for a
rigorous analysis). More information on the thinning process in RMT can be
found, in particular, in [23].

Another open problem in RMT is to prove universality, in the bulk and at (all)
the edges, for Wigner ensembles in the case that the equilibrium measure for
the ensemble is supported on more than one interval. Furthermore, one would
like to analyze and control situations where parameters in the underlying
matrix distributions vary, and intervals in the support of the equilibrium
distribution for the ensemble open and close. By contrast with the Wigner
case, these problems are well understood for unitary invariant ensembles, and
to a slightly lesser extent, for invariant orthogonal and symplectic
ensembles. For example, for the unitary invariant ensemble with probability
distribution proportional to $e^{-n\operatorname{tr}V(M)}dM$, where
$V(x)=-ax^{2}+bx^{4}$, with $b>0$ fixed and $a\in\mathbb{R}$, there is a
critical value $a_{c}$ such that for $a<a_{c}$, the equilibrium measure for
the ensemble is supported on one interval, and for $a>a_{c}$ the measure is
supported on two intervals. In both cases, the eigenvalues exhibit universal
behavior in the bulk away from the end points of support of the equilibrium
measure. In the critical case $a=a_{c}$, when the two intervals merge at some
point $x=x_{c}$, the behavior of the eigenvalues in a neighborhood of $x_{c}$
corresponds to a different universality class, which can be analyzed and
described explicitly using, for example, the results in [36]. Similar problems
are completely open for Wigner ensembles.

Double-scaling limits and, more generally, multi-scaling limits, are of great
interest in physics and in mathematics. For example, regarding the classical
Ising model in two dimensions, much work has been done on analyzing the
double-scaling limit for the two-point correlation function as the temperature
$T$ goes to the critical temperature $T_{c}$ and simultaneously the separation
$x$ between the spins goes to infinity. It turns out, quite remarkably, that
the limit is described by a particular solution of the Painlevé III equation
[113]. Apart from the intrinsic interest in this result, it is widely believed
in the physics community that the behavior of the correlation function in the
double-scaling limit is universal, and precisely the same behavior should be
expected for a very wide class of statistical spin models. In other words,
double-scaling limits can provide access to universal behavior that transcends
the problem at hand. Indeed, the random band matrix problem described above as
the band width $w$ and the matrix size $n$ go to infinity, provides a model
from which we hope to learn about localization in the bulk for random
Schrödinger operators. In RMT, various double-scaling problems have been
considered. The first rigorous analysis is due to P. Bleher and A. Its [9] who
considered the unitary invariant ensemble with potential $V(x)=-ax^{2}+bx^{4}$
as above, in the limit as $a$ converges to the critical value $a_{c}$ and the
size $n$ of the matrices goes to infinity. Again Painlevé equations and
associated functions arise in the analysis. In 2006, T. Claeys and A.
Kuijlaars [25] generalized the results of Its and Bleher to a wider class of
potentials $V(x)$, and in the present volume T. Claeys and B. Fahs [24] have
analyzed a very interesting double-scaling problem where the Painlevé I
equation plays a role. There have also been numerous non-rigorous analyses of
double-scaling problems in RMT, both within the mathematical community and the
physics community. In RMT almost all problems involve matrix distributions
which depend on some external parameters, say $p$. Typically there are
critical values $p_{c}$ at which the nature of the equilibrium measure for the
ensemble changes. The double-scaling problems, $p\to p_{c}$ and matrix size
$n\to\infty$, are of basic interest. Most of these problems remain open.

An open problem of long-standing interest concerns explicit formulae for the
Tracy–Widom distribution functions for general $\beta$-ensembles. By the
universality results mentioned above, we know that for each $\beta$, the
distribution is universal, but up till now explicit expressions were known
only for the invariant ensembles, $\beta=1,2$ and $4$. However, recently
Rumanov [88, 89] derived an explicit formula for the Tracy–Widom distribution
in the case $\beta=6$ in terms of Painlevé II and the solution of an auxiliary
ODE. Although parts of Rumanov’s analysis are not rigorous (see [60] where
some of the gaps are filled), Rumanov’s work constitutes a very significant
development. As noted by Rumanov, his analysis extends, in principle, to the
Tracy–Widom distributions for all even $\beta>0$. It is of great interest to
complete Rumanov’s analysis rigorously, not only for $\beta=6$, but for all
even $\beta>0$. Using loop equation methods, G. Borot, B. Eynard, S. Majumdar
and C. Nadal [16] have written down very detailed formulae for the asymptotic
behavior of the Tracy–Widom distributions for general $\beta>0$. It is a
challenging open problem to verify rigorously these formulae of Borot et al.

There are two aspects to random matrix theory. The first aspect concerns
internal questions within RMT, and most of the problems that we have
considered so far, are of this nature. The second aspect concerns applications
of RMT. The situation is analogous to the theory of orthogonal polynomials,
where researchers study orthogonal polynomials per se, on the one hand, and
their applications, on the other. As new applications arise, they can, and
often do, raise new structural questions within the theory of orthogonal
polynomials. The same is true for RMT, with RMT playing the role of a
“stochastic special function theory”, so extending the reach of classical
special function theory to the nonlinear problems of modern mathematical
physics. In Problem 3 we will discuss applications of RMT to interacting
particle systems and in Problem 4 we will discuss applications of RMT to
numerical analysis. As we will see, numerical analysis raises, in particular,
many new, deep and very interesting problems in RMT.

######  Problem 3 (interacting particle systems and KPZ).

In 1999, J. Baik, P. Deift and K. Johansson [3] solved the longest increasing
subsequence problem, establishing en route a connection between combinatorics
and RMT (Tracy–Widom distribution). Shortly thereafter, many statistical
models/interacting particle systems were solved in related ways, first by K.
Johansson [63], and then by A. Borodin, A. Okounkov and G. Olshanski [13],
among many others. All these early examples were “integrable” in the sense
that certain powerful algebraic tools, e.g., determinantal particle systems,
the Robinson–Schensted–Knuth correspondence, Fredholm operator theory,
$\dots$, could be applied to reduce the problems at hand to analytically
viable forms from which their asymptotic behavior could be deduced using
standard asymptotic methods such as the Laplace method/stationary phase, or
sometimes, Riemann–Hilbert/nonlinear steepest-descent methods. In a very
significant and unanticipated development during 2007–2009, C. Tracy and H.
Widom [103, 104, 105] analyzed ASEP (the asymmetric simple exclusion process)
with step initial data, a system to which the above algebraic tools cannot be
applied. Nevertheless, in an algebraic/combinatorial/analytical tour de force,
they were able to show that ASEP exhibited the same asymptotic RMT behavior as
the integrable models.

A striking feature of all these models, is that fluctuations are of order
$N^{1/3}$, where $N\to\infty$ is the size of the system, thereby showing that
the fluctuations for these models are in the KPZ (Kardar–Parisi–Zhang)
universality class. Note that, by contrast, systems in the Gaussian
universality class have fluctuations of order $N^{1/2}$. The notion of KPZ
universality comes from the work of M. Kardar, G. Parisi and Y.-C. Zhang [69]
who in 1986 introduced an equation, the KPZ equation, to describe the growth
of random interfaces with height function $h(x,t)$, $x\in\mathbb{R}$, $t\geq
0$. The KPZ equation is stochastic and non-linear, and using various physical
arguments they showed that as $t\to\infty$, the surface should be correlated
over distances of order $t^{2/3}$ and the fluctuations of $h$ should be of
order $t^{1/3}$. In other words there should be a $3:2:1$ scaling for time :
spatial correlations : fluctuations. It is widely believed that a large class
of physical systems and mathematical models obey KPZ $3:2:1$ universality. The
above integrable models were the first to be shown rigorously to have KPZ
fluctuations $3:1$ (i.e., $N:N^{1/3}$). In addition, in 2000, K. Johansson
[64] showed that for the last passage percolation problem (LPP), spatial
correlations were of order $t^{2/3}$ in accord with $3:2$ KPZ universality,
and in later more detailed work Johansson [66] showed that as $t\to\infty$ the
fixed $t$ marginal for $h(x,t)$ converged to the Airy 2 process. We note that
M. Prähofer and H. Spohn [86] were the first to show that $h(x,t)$ converged
to the Airy 2 process, but their argument was not rigorous (see [26] for more
information on spatial correlations).

A basic question is whether KPZ is in the KPZ universality class. In other
words, can one show rigorously that solutions of KPZ indeed have the KPZ
scaling $3:2:1$? The first difficulty here is to show that the Cauchy problem
for the non-linear stochastic KPZ equation indeed has a solution. Following L.
Bertini and N. Cancrini [5], one approach is to use the Cole–Hopf
transformation which converts the KPZ equation into a heat equation with a
random potential which can then be solved via a generalized Feynmann–Kac
formula. In breakthrough work in 2010, G. Amir, I. Corwin and J. Quastel [2],
used results from the 2008–2009 work of Tracy and Widom [103, 104, 105] on
ASEP to find an explicit formula for the solution of the Cole–Hopf transformed
KPZ equation with so-called “narrow wedge” initial data. From this formula
they were able to extract the $t^{1/3}$ KPZ scaling for the fluctuations of
the solution, but their methods give no information on spatial correlations.
Simultaneous with the work of Amir, Corwin and Quastel on the KPZ equation,
there was also the independent, related, but not rigorous work, of T. Sasamoto
and H. Spohn [90] also relying on Tracy–Widom’s ASEP formula, as well as the
work of V. Dotsenko [48], and P. Calabrese, P. Le Doussal and A. Russo [22]
using replicas. The varying levels of rigor of these papers are discussed in
[26]. Regarding full $3:2:1$ KPZ universality, the strongest result is due to
I. Corwin and A. Hammond [28] (see [30] for the strongest and most precise
statement of the KPZ universality conjecture for the KPZ equation). Again for
the stochastic heat equation with narrow wedge initial data, Corwin and
Hammond were able to show tightness of the suitably $3:2:1$ scaled height
function $h(x,t)$ in a “good” sense which implied convergence for suitable
subsequences to limiting processes with “good” properties, but they were not
able to show that the limits were unique and given by the physically
anticipated process. In particular, they could not identify the anticipated
Airy2 marginal. We note here that M. Hairer [61] was able to show directly
without using the Cole–Hopf transformation that the Cauchy problem for the KPZ
equation has a solution. This was a major achievement that figured prominently
in the citation for Hairer when he won the Fields Medal in 2014. Under the
appropriate conditions, Hairer’s solution agrees with a Cole–Hopf solution.
For more information on KPZ, see, for example [12, 15, 26, 27].

The most challenging open problem regarding the KPZ equation is to prove
$3:2:1$ KPZ universality in the strong sense of Corwin, Quastel and Remenik
[30] alluded to above, for solutions of the Cauchy problem of the KPZ equation
for a wide set of initial data, and to identify the limiting process. Another
open problem is to analyze solutions of KPZ at more than one space-time point.
Regarding these issues, much significant progress has been made for various
integrable particle systems. For example in 2015 K. Johansson [67] computed
correlations for two space-time points for the Brownian LPP problem, from
which it easily follows, in particular, that Brownian LPP exhibits $3:2:1$ KPZ
scaling. We note that this is the only system for which $3:2:1$ KPZ scaling
has been verified rigorously. Also, in 2015, I. Corwin, Z. Liu and D. Wang
[29] demonstrated KPZ universality for the fluctuations for LPP and TASEP with
general initial data. However, many basic problems remain open. In particular
showing that LPP with general weights has $3:2:1$ scaling, is far from being
resolved. Also, regarding the original longest increasing subsequence problem,
KPZ universality for the fluctuations has been proved so far only for the
uniform distribution on the permutations: The problem for the fluctuations
with more general distributions is very much open. In another direction, it
remains an open, and somewhat mysterious, problem to understand how Tracy and
Widom’s solution of ASEP actually “works”, and to extend their method of
solution to more general initial data (however, see [12], and more recently,
[14]).

Finally we note that, with few exceptions, all the rigorous asymptotic scaling
results that have been obtained so far for KPZ/interacting particle systems,
have been obtained only for restricted classes of initial data. In addition,
it is particularly frustrating that there is no known perturbation theory,
analogous to KAM theory for integrable dynamical systems, that one can use to
investigate the persistence of “integrability” for non-solvable models and
demonstrate robustness for asymptotic scaling phenomena. Developing such a
perturbation theory is an open conceptual and technical problem of the first
order.

######  Problem 4 (numerical computation with random data).

Standard algorithms to compute the eigenvalues of a random matrix $H$ are
completely integrable Hamiltonian systems (see [38, 39, 100]). The situation
is as follows. Let $\Sigma_{n}$ denote the set of real $n\times n$ symmetric
matrices. Associated with each algorithm $\mathcal{A}$, there is, in the
discrete case such as the QR algorithm, a map
$\varphi=\varphi_{\mathcal{A}}\colon\Sigma_{n}\to\Sigma_{n}$, with the
properties

$\displaystyle\bullet\text{(isospectral)}\quad$ $\displaystyle\text{spec
}(\varphi_{\mathcal{A}}(H))=\text{spec $(H)$},$
$\displaystyle\bullet\text{(convergence)}\quad$ $\displaystyle\text{the
iterates}\ X_{k+1}=\varphi_{\mathcal{A}}(X_{k}),\ k\geq 0,$ $\displaystyle
X_{0}=H\quad\text{given, converge to a diagonal matrix $X_{\infty}$ as
$k\to\infty$},$

and in the continuum case, such as the Toda algorithm, there is a flow $t\to
X(t)\in\Sigma_{n}$ with the properties

$\displaystyle\bullet\text{(isospectral) }\quad$
$\displaystyle\text{spec}(X(t))\text{ is constant, $t\geq 0$},$
$\displaystyle\bullet\text{(convergence)}\quad$ $\displaystyle\text{the f\/low
$X(t)$, $t\geq 0$, $X(0)=H$ given},$ $\displaystyle\text{converges to a
diagonal matrix $X_{\infty},X(t)\to X_{\infty}$ as $t\to\infty$}.$

In both cases, necessarily the (diagonal) entries of $X_{\infty}$ are the
eigenvalues of the given matrix $H$. The flow $X(t)$ in the continuum case is
completely integrable and in the discrete case $X_{k}$ is the $\text{time}-k$
evaluation of a completely integrable flow (see references above). Given
$\epsilon>0$, it follows, in the discrete case, that for some $m$ the off-
diagonal entries of $X_{m}$ are $O(\epsilon)$ and hence the diagonal entries
of $X_{m}$ give the eigenvalues of $X_{0}=H$ to $O(\epsilon)$. The situation
is similar for continuous algorithms $t\to X(t)$.

On the other hand, matrices carry another natural integrable structure. This
is the structure of RMT in which basic statistics, such as the gap probability
or the nearest-neighbor distribution, converge as the size $n$ of the matrices
becomes large to distributions that are expressed in terms of the solution of
integrable Hamiltonian systems (in particular, Painlevé equations). The
question raised in [34] was: What happens if we try to marry these two
integrabilities? Or more precisely, what can be said about the computation of
the eigenvalues of random matrices using standard algorithms?

This question was taken up by P. Deift, G. Menon and C. Pfrang in 2012 [85].
What they found was that the computations exhibited universality, with a
different universality class for each algorithm $A$. More precisely, for an
$n\times n$ matrix $H$ chosen from an ensemble $\mathcal{E}$ on $\Sigma_{n}$,
and a given accuracy $\epsilon>0$, let
$T(H)=T_{\epsilon,n,\mathcal{A},\mathcal{E}}(H)$ be the smallest $k$, or
smallest $t$ in the continuum case, for which $X_{k}$ with $X_{0}=H$, evolving
under algorithm $\mathcal{A}$, has off diagonal entries (see “deflation time”
in [85] and also $T^{(1)}$ below) less than $\epsilon$. The authors in [85]
then recorded

$\displaystyle\tau_{\epsilon,n,\mathcal{A},\mathcal{E}}(H)\equiv\frac{T_{\epsilon,n,\mathcal{A},\mathcal{E}}(H)-\langle
T_{\epsilon,n,\mathcal{A},\mathcal{E}}\rangle}{\sigma_{\epsilon,n,\mathcal{A},\mathcal{E}}},$

where $\langle T_{\epsilon,n\mathcal{A},\mathcal{E}}\rangle$ and
$\sigma^{2}_{\epsilon,n,\mathcal{A},\mathcal{E}}$ are respectively the sample
average and sample variance for $T_{\epsilon,n,\mathcal{A},\mathcal{E}}(H)$
computed for a very large sample of matrices $H\in\mathcal{E}$. What the
authors found was that, for a given algorithm $\mathcal{A}$, and $\epsilon$,
$n$ in a suitable scaling region, the histogram for
$\tau_{\epsilon,n,\mathcal{A},\mathcal{E}}$ was universal, independent of the
ensemble $\mathcal{E}$. This phenomenon was present, in particular, for such
disparate distributions as GOE and Bernoulli. The histogram for
$\tau_{\epsilon,n,\mathcal{A},\mathcal{E}}$ does, however, depend on
$\mathcal{A}$ and is different, in particular, for the $QR$ and Toda
algorithms.

In [40] the authors P. Deift, G. Menon, S. Olver and T. Trogdon raised the
question of whether the universality results in [85] were limited to
eigenvalue algorithms, or whether they were present more generally in
numerical computations. And indeed in [40] the authors found similar
universality results for a wide variety of numerical algorithms, including

  * •

more general eigenvalue algorithms such as the Jacobi algorithm, and also
algorithms for Hermitian ensembles,

  * •

the conjugate gradient and GMRES algorithms to solve linear $n\times n$
systems $Hx=b$, $H$ and $b$ random,

  * •

an iterative algorithm to solve the Dirichlet problem $\Delta u=0$ on a random
star-shaped region $\Omega\subset\mathbb{R}^{2}$ with random boundary data $f$
on $\partial\Omega$, and

  * •

a genetic algorithm to compute the equilibrium measure for orthogonal
polynomials on the line, and was also for a decision making process in recent
laboratory experiments of Yu. Bakhtin and J. Correll [4].

Until recently, all the evidence for universality in computation was numerical
and experimental. In order to establish the phenomenon rigorously, P. Deift
and T. Trogdon [42] focused on the Toda algorithm, $\mathcal{A}=\text{Toda}$,
applied to $n\times n$ Hermitian matrices $H$ chosen from some ensemble
$\mathcal{E}$, with stopping time
$T^{(1)}(H)=T^{(1)}_{\epsilon,n,\text{Toda},\mathcal{E}}(H)$, which is the
first time that

$\displaystyle\sum^{n}_{j=2}|X_{1j}(t)|^{2}=\sum^{n}_{i=2}|X_{i1}(t)|^{2}<\epsilon^{2},$

where $X(t)$ solves the Toda flow with $X(0)\\!=\\!H$. For
$t\\!=\\!T^{(1)}(H)$ we clearly have
$|X_{11}(t)\\!-\\!\lambda^{\ast}|\\!<\\!\epsilon$ for some eigenvalue
$\lambda^{\ast}$ of $H$. By the ordering properties of the Toda algorithm,
$\lambda^{\ast}=\lambda_{\max}$, the largest eigenvalue of $H$, with high
probability as $n\to\infty$. Thus the Toda algorithm with stopping time
$T^{(1)}(H)$ is an algorithm to compute the largest eigenvalue of a given
matrix $H$.

In [42] the authors proved the following result:

######  Theorem 1 (universality for the Toda algorithm).

Let $0<\tau<1$ be fixed and let $(\epsilon,n)$ be in the scaling region

$\displaystyle\frac{\log\epsilon^{-1}}{\log n}\geq\frac{5}{3}+\frac{\tau}{2}.$

Then if $H$ is distributed according to any real $(\beta=1)$ or complex
$(\beta=2)$ invariant or Wigner ensemble111For a precise description of the
ensembles allowed, see [42].

$\displaystyle\lim_{n\to\infty}\operatorname{Prob}\left(\frac{T^{(1)}}{C^{2/3}_{v}2^{-2/3}N^{2/3}\left(\log\epsilon^{-1}-\frac{2}{3}\log
n\right)}\right)=F^{\text{\rm gap}}_{\beta}(t).$

Here $C_{v}$ is an explicit ensemble dependent constant and for $\beta=1,2$,

$\displaystyle
F^{\text{gap}}_{\beta}(t)=\lim_{n\to\infty}\operatorname{Prob}\left(\frac{1}{C^{2/3}_{v}2^{-2/3}n^{2/3}(\lambda_{n}-\lambda_{n-1})}\leq
t\right),$

where
$\lambda_{\max}=\lambda_{n}>\lambda_{n-1}>\dots>\lambda_{j}>\dots>\lambda_{1}$
are the eigenvalue of a random matrix $H$. Properties of $G_{\beta}(t)\equiv
1-F_{\beta}^{\text{gap}}(t)$, the distribution function for the first gap are
examined, for example, in [112]. The authors also prove as a corollary that
$\epsilon^{-1}\big{|}\lambda_{\max}-X_{11}\big{(}T^{(1)}\big{)}\big{|}$
converges to zero as $n\to\infty$ in the scaling region.

The Toda system has the form

$\displaystyle\frac{dX}{dt}=[X,B(X)],\qquad B(X)=X_{-}-X^{*}_{-},\qquad
X(0)=H,$

where $X_{-}$ is the strictly lower-triangular part of $X$, $X^{*}_{-}$ is the
adjoint of $X_{-}$, and $[A,B]$ is the standard matrix commutator. The proof
of the above theorem relies in a crucial way on the integrability of the Toda
system which makes it possible to reduce the proof to the analysis of explicit
formulae. Indeed, following Moser [81] one finds that

$\displaystyle
E(t)\equiv\sum^{n}_{j=2}|X_{ij}(t)|^{2}=\sum^{n}_{j=1}(\lambda_{j}-X_{11}(t))^{2}|u_{1j}(t)|^{2},$
(3)
$\displaystyle\lambda_{j}-X_{11}(t)=\sum^{n-1}_{i=1}(\lambda_{j}-\lambda_{i})|u_{1i}(t)|^{2},$
(4)

where $u_{j}(t)=(u_{1j}(t),\dots,u_{nj}(t))^{T}$ is the $j^{\rm th}$
normalized eigenvector of $X(t)$, $X(t)u_{j}(t)=\lambda_{j}u_{j}(t)$, which
evolves as

$\displaystyle|u_{1j}(t)|=\frac{\beta_{j}e^{\lambda_{j}t}}{\left(\sum\limits^{n}_{i=1}\beta^{2}_{i}e^{2\lambda_{i}t}\right)^{\frac{1}{2}}},$
(5)

where $\beta_{j}=|u_{ij}(0)|$ are the moduli of the first components of the
eigenvectors of $X(0)=H$. Substituting (4), (5) into (3) we obtain an explicit
formula for $E(t)$ in terms of the eigenvalues $\\{\lambda_{j}\\}$ and the
moduli $\\{\beta_{j}\\}$ of the first components of the normalized
eigenvectors of $H=X(0)$. To obtain $T^{(1)}$, we must solve for the first $t$
such that $E(t)=\epsilon^{2}$. We see, in particular, that the statistics of
$T^{(1)}$ depends on the statistics of $\\{\lambda_{j}\\}$ and
$\\{\beta_{k}\\}$. More precisely, in order to analyze the statistics of
$T^{(1)}$, it turns out that one needs to control the probabilities for the
following events as $n\to\infty$:

######  Condition 1.

For $0<p<\frac{1}{3}$

  1. (i)

$\lambda_{n-1}-\lambda_{n-2}\geq p(\lambda_{n}-\lambda_{n-1})$.

Let $G_{n,p}$ denote the set of matrices that satisfy this condition.

######  Condition 2.

For any fixed $0<s<\min\\{\sigma/44,p/8\\}$

  1. (i)

$\beta_{j}\leq n^{-\frac{1}{2}+s/2}$ for all $j$,

  2. (ii)

$n^{-\frac{1}{2}-s/2}\leq\beta_{j}$ for $j=n,n-1$,

  3. (iii)

$n^{-2/3-s}\leq\lambda_{n}-\lambda_{j}\leq n^{-2/3+s}$ for $j=n-1,n-2$, and

  4. (iv)

$|\lambda_{j}-\gamma_{j}|\leq n^{-2/3+s}(\min\\{j,n-j+1\\})^{-1/3}$ for all
$j$.

Here the $\gamma_{j}$’s represent the quantiles of the equilibrium measure
$d\mu$, i.e., $\gamma_{j}$ is defined to be the smallest value of $b$ such
that $j/n=\int^{b}_{-\infty}d\mu$. Let $R_{n,s}$ denote the set of matrices
that satisfy these conditions (i)–(iv).

In order to prove Theorem 1 and the corollary, one needs to show that for both
Wigner ensembles (WE’s) and invariant ensembles (IE’s), Condition 1 holds with
probability 1 as $p\downarrow 0$, i.e.,

$\displaystyle\lim_{p\downarrow
0}\limsup_{n\to\infty}\operatorname{Prob}\big{(}G^{c}_{np}\big{)}=0$ (6)

and Condition 2 holds with high probability as $n\to\infty$, that is, for any
$s>0$

$\displaystyle\operatorname{Prob}(R_{n,s})=1+o(1)$ (7)

as $n\to\infty$.

In addition one needs the following results: For both WE’s and IE’s, as
$n\to\infty$

$\displaystyle n^{-1/2}\\{\beta_{n},\beta_{n-1},\beta_{n-2}\\}$ (8)

converge jointly in distribution to $\\{|X_{1}|,|X_{2}|,|X_{3}|\\}$ where
$\\{X_{1},X_{2},X_{3}\\}$ are iid real $(\beta=1)$ or complex $(\beta=2)$
standard normal random variables. Additionally

$\displaystyle
2^{-2/3}n^{2/3}(b_{v}-\lambda_{n},b_{v}-\lambda_{n-1},b_{v}-\lambda_{n-2})$
(9)

converge jointly in distribution to random variables
$\\{\Lambda_{1,\beta},\Lambda_{2,\beta},\Lambda_{3,\beta}\\}$ which are the
smallest three eigenvalues of the stochastic Airy operator. Furthermore
$(\Lambda_{1,\beta},\Lambda_{2,\beta},\Lambda_{3,\beta})$ are distinct with
probability one. Here $b_{v}$ is the top of the support of the equilibrium
measure $d\mu$. In order to prove (6)–(9), one needs to draw on results from
the forefront of research in RMT and that were proved only very recently. In
particular, one needs to draw on results from [19, 21, 52, 55, 87]. How these
results are used to prove (6)–(9) is described in detail in [42].

It is clear from the above considerations that the analysis of numerical
algorithms with random data is a rich source of new and challenging problems
in RMT with immediate practical applications. Here are some open problems. In
order to prove universality for the centered and scaled variable
$\frac{T^{(1)}-\langle T^{(1)}\rangle}{\sigma_{T^{(1)}}}$ (cf.
$\tau_{\epsilon,n,\mathcal{A},\mathcal{E}}$ at the beginning of this section)
one needs to control the statistics of the $\lambda_{j}$’s and the
$\beta_{k}$’s as in (8), (9) jointly as $n\to\infty$. For $\beta=2$ invariant
ensembles, this can be done, but for $\beta=1$ invariant ensembles, and for
$\beta=1$ or 2 WE’s, this problem is still very much open. For the $QR$
algorithm, the role of the eigenvalue gaps $\lambda_{k}-\lambda_{j}$ is now
played by the eigenvalue ratios $\lambda_{k}/\lambda_{j}$. So instead of the
behavior of $T^{(1)}=T^{(1)}_{\epsilon,n,QR,\mathcal{E}}$ being controlled by
the inverse gap $(\lambda_{n}-\lambda_{n-1})^{-1}$ as in Theorem 1, $T^{(1)}$
will now be controlled by $\big{|}\frac{\lambda_{j+1}}{\lambda_{j}}\big{|}$
and $\big{|}\frac{\lambda_{j-1}}{\lambda_{j}}\big{|}$ where $\lambda_{j}$ is
the eigenvalue of $H$ closest to zero. It is very much an open problem to
control, as $n\to\infty$, the statistics of
$\big{|}\frac{\lambda_{j+1}}{\lambda_{j}}\big{|}$,
$\big{|}\frac{\lambda_{j-1}}{\lambda_{j}}\big{|}$ for $\lambda_{j}\sim 0$ deep
in the bulk of the spectrum, together with their associated normalized
eigenvectors, as in (8), (9). The $QR$ algorithm is the “go to” algorithm for
eigenvalue computation, and the proof of the analog of Theorem 1 for the $QR$
algorithm is an open problem of great interest. For recent work in proving
universality for the $QR$ algorithm and related algorithms for strictly
positive definite matrices, see [43].

The most challenging problem for eigenvalue computation with random data is to
prove the analog of Theorem 1 for the deflation time
$T_{\epsilon,n,\mathcal{A},\mathcal{E}}$ (see [42] and the references therein)
for some algorithm $\mathcal{A}$. If $\mathcal{A}=\text{ Toda}$, for example,
there is an analog $E_{\text{def}}(t)$ of $E(t)$ in (3). Again one computes
$T_{\epsilon,n,\mathcal{A},\mathcal{E}}$ as the smallest time for which
$E_{\text{def}}(t)=\epsilon^{2}$, but now $E_{\text{def}}(t)$ is expressed in
terms of all the components of the eigenvectors of $X(0)=H$, in addition to
the eigenvalues. We refer the reader to [42] where one can see precisely what
statistics on the eigenvalues and eigenvectors are needed to prove
universality for $T_{\epsilon,n,\mathcal{A},\mathcal{E}}$, and how the problem
can be broken down into smaller sub-problems of interest. Indeed,
$T_{\epsilon,n,\mathcal{A},\mathcal{E}}$ is the smallest time $t$ for which
$X(t)$ has the block form $\left(\begin{smallmatrix}A&B\\\
B^{*}&C\end{smallmatrix}\right)$ where $A$ is $j\times j$, $C$ is
$(n-j)\times(n-j)$ and $\|B\|=\|B^{*}\|\leq\epsilon$ for some $1\leq j\leq
n-1$. For such $t$, the eigenvalues of the deflated matrix
$\left(\begin{smallmatrix}A&0\\\ 0&C\end{smallmatrix}\right)$ agree with the
eigenvalues of $X(0)=H$ to $O(\epsilon)$. The algorithm with deflation then
continues by applying $\mathcal{A}$ to the (smaller) matrices $A$ and $C$,
etc. Note that the time $T^{(1)}$ above controls deflation to the form
$\left(\begin{smallmatrix}A&0\\\ 0&C\end{smallmatrix}\right)$ where $A$ is
$1-1$ and $C=(n-1)\times(n-1)$. The next step in analyzing
$T_{\epsilon,n,\mathcal{A},\mathcal{E}}$ is to control $T^{(2)}$, the time for
deflation to the form $\left(\begin{smallmatrix}A&0\\\
0&C\end{smallmatrix}\right)$ where $A$ is $2\times 2$ and $C\text{ is
}(n-2)\times(n-2)$. To control $T^{(2)}$ would require further advances in
RMT, in particular, statistical information on the second components of the
eigenvectors, and so on.

It is an interesting and challenging open problem to prove universality
rigorously for any of the other numerical algorithms discussed in [40]. The
conjugate gradient algorithm, in particular, to compute the solution of the
linear system $Ax=b$ with $A$ and $b$ random, seems to be “within reach”. In
this connection, we refer the reader to [44].

Finally we consider an open statistical problem. In order to test a given
algorithm for universality, one needs to sample various ensembles. This is a
standard matter if one is dealing with Wigner ensembles and the entries are
independent, but when the entries are dependent, as in the case for general
invariant ensembles, this is a non-trivial matter. In 2013, X. Li and G. Menon
[74] developed an approach for sampling invariant ensembles based on
simulating Dyson Brownian motion. In 2014, for invariant unitary ensembles, S.
Olver, R. Nadakuditi and T. Trogdon [83] developed an alternative approach
based on computing the orthogonal polynomials naturally associated with the
ensembles. While the method of Li and Menon is significantly faster, the
method of Olver et al. samples the distribution to essentially machine
accuracy, a task which is computationally impractical using Dyson Brownian
motion. It is an open problem to extend the method of S. Olver et al. to
orthogonal and symplectic invariant ensembles.

######  Problem 5 (initial boundary value problems for integrable systems
(IBVP)).

Initial boundary value problems for integrable systems in $1+1$ dimensions are
of great interest. In seminal work in the late 1990s A. Fokas [56] introduced
a new, more flexible approach to inverse scattering, and over the years many
researchers (A. Fokas, A. Its, A. Boutet de Monvel, D. Shepelsky, …) have
applied Fokas’ approach to IBVPs for various integrable systems. It turns out,
however, that in Fokas’ approach the general initial value problem is
overdetermined. For example, in order to solve NLS

$\displaystyle iu_{t}=u_{xx}\pm 2|u|^{2}u$

in $(x\geq 0$, $t\geq 0)$ one needs to know $u_{x}(0,t)$, $t\geq 0$, in
addition to the standard determining data

$\displaystyle u(x,0)=f(x),\qquad$ $\displaystyle x\geq 0,$ $\displaystyle
u(0,t)=g(t),\qquad$ $\displaystyle t\geq 0.$

In certain “integrable” cases (see below) $u_{x}(0,t)$ can be determined a
priori, but this is not the case in general. As demonstrated by Fokas, his
approach leads to a Riemann–Hilbert reformulation of the IBVP, but one can
only apply the powerful non-linear steepest descent method to determine the
asymptotic behavior of $u(x,t)$ as $t\to\infty$ in the case that $u_{x}(0,t)$
is known a priori. A particular stand out open problem is to compute the long-
time behavior of the solution $u(x,t)$ of the focusing NLS equation,
$iu_{t}=u_{xx}+2|u|^{2}u$, with $u(x,0)=f(x)$, $x\geq 0$, given and
$u(0,t)=e^{iwt}$, $t\geq 0$ with $w\in\mathbb{R}$. Interesting partial results
for this problem as $w\to\infty$ are given in A. Fokas and J. Lenells [57].

What kind of behavior do we anticipate for $u(x,t)$ as $t\to\infty$,
$u(0,t)=e^{iwt}$? From a physical point of view one can think of a long rope
which is shaken a one end with frequency $w/2\pi$. Experience tells us that
for $w$ very large, only a very small disturbance will be generated in the
rope. On the other hand, as we shake the rope slower and slower, more
complicated and sizable wave forms propagate along the rope. In an analogous
situation, P. Deift, T. Kriecherbauer and S. Venakides [37] considered a
lattice of particles $x_{k}$, $k\geq 0$, interacting with Toda forces

$\displaystyle\ddot{x}_{k}=e^{x_{k-1}-x_{k}}-e^{x_{k}-x_{k+1}},\qquad k\geq 1$

and initial/boundary value conditions

$\displaystyle x_{k}(0)=\dot{x}_{k}(0)=0,\qquad k\geq 1,$ $\displaystyle
x_{0}(t)=2at+h(\gamma t).$

Here $h(\cdot)$ is $2\pi$-periodic and $a$ and $\gamma>0$ are given constants,
with $a$ below a certain critical value $a_{\text{crit}}$. The authors in [37]
observed the following phenomena: In a frame moving with the average velocity
of the drive $x_{0}(t)$, as $t\to\infty$ the asymptotic motion of the
particles behind the shock front is $2\pi/\gamma$-periodic in time,

$\displaystyle
x_{k}\left(t+\frac{2\pi}{\gamma}\right)=x_{k}(t),\qquad\text{for f\/ixed
$k\geq 1$.}$

Moreover, there is a sequence of thresholds

$\displaystyle\gamma_{1}>\gamma_{2}>\dots>\gamma_{j}>\dots>0,\qquad\gamma_{j}\to
0\quad\text{as}\quad j\to\infty$

with the following properties:

  * •

If $\gamma>\gamma_{1}$, there exist constants $c,d$ such that $x_{k}-ck-d$
converges exponentially to zero for $1\ll k\ll t$. In other words, the effect
of the oscillatory component of the driver does not propagate into the lattice
and away from the boundary at $k=0$.

  * •

If $\gamma_{1}>\gamma>\gamma_{2}$, then the asymptotic motion is described by
a travelling wave

$\displaystyle x_{k}(t)=c_{1}k+Y_{1}(\beta_{1}k+\gamma t),\qquad 1\ll k\ll t$

transporting energy away from the driver. Here $c_{1}$ and $\beta_{1}$ are
certain determined constants and $Y_{1}(1)$ is $2\pi$-periodic 1-gap solution
of the periodic Toda lattice.

  * •

More generally, if $\gamma_{j}>\gamma>\gamma_{j+1}$, a multi-phase wave
emerges of the form

$\displaystyle x_{k}=c_{j}k+Y_{j}(\beta_{1}k+\gamma t,\beta_{2}k+\gamma
t,\dots,\beta_{j}k+\gamma t),\qquad 1\ll k\ll t$

again transporting energy away from the driver, for certain constants
$c_{j},\beta_{1},\dots,\beta_{j}$. Now $Y_{j}(\cdot,\dots,\cdot)$ is $j$-gap
solution of the periodic Toda lattice, $2\pi$-periodic in each of its
variables.

Thus, at the phenomenological level, we see that indeed the periodically
driven lattice behaves like a long rope which one shakes up and down at the
one end. We expect similar behavior for focusing NLS, with $k$-gap solutions
of periodic NLS emerging from the driver $e^{iwt}$ as $w\downarrow 0$. This is
an open problem of the first order, illustrating how the natural modes of an
extended physical system ($k$-gap solutions in the case of NLS) can be excited
by driving the system externally at one end.

In 2010, P. Deift and J. Park [41] considered the IPVP for focusing NLS with
the boundary condition $u(0,t)=g(t)$ replaced by a Robin boundary condition,

$\displaystyle u_{x}(0,t)+qu(0,t)=0$

with $q$ fixed. This problem arises, in particular, in analyzing the
Gross–Pitaevskii equation with a delta function external potential at $x=0$
and even initial data. It turns out that the NLS equation with Robin boundary
conditions is an example of an integrable system which is not overdetermined
from the point of view of Fokas’ method, and indeed for this problem the
method can be used to obtain a well-defined and effective Riemann–Hilbert
representation for the solution $u(x,t)$ of the IBVP (see [62]). In their 2010
article, however, Deift and Park [41] used an earlier method introduced by R.
Bikbaev and V. Tarasov in 1991 [6]. The method of Bikbaev and Tarasov is a
non-linear version of the method of images in linear theory, where the
solution in $x>0$ is extended to $x<0$, not by reflection as in the linear
case, but by a suitable Bäcklund transformation. In this way the IBVP is
converted smoothly into a standard problem on the whole line, to which
Riemann–Hilbert/steepest-descent methods apply directly. It is an interesting
open problem to use the Bikbaev–Tarasov method to solve the IBVP for
defocusing NLS, $iu_{t}=u_{xx}-2|u|^{2}u$, with Robin boundary conditions at
$x=0$. In this case the Bäcklund transformation extending the solution from
$x>0$ to $x<0$ introduces certain singularities which must be controlled and
the long-time behavior of the solution of the IBVP will be very different from
the case of focusing NLS. Another open problem, which has broad implications,
concerns the smoothness of the solutions $u(x,t)$ of the IBVP for NLS. One
expects that if one solves the IBVP for NLS with $f(x)=u(x,0)$ and
$g(t)=u(0,t)$ sufficiently smooth, and with the right matching conditions at
the corner $x=0$, $t=0$, then the solution $u(x,t)$ should be classical in
$x>0$, $t>0$, in the sense that $u(x,t)$ is $C^{2}$ in $x$ and $C^{1}$ in $t$.
However in the case of focusing NLS with smooth $f(x)=u(x,0)$ and Robin
boundary conditions with $q\neq 0$ (the case $q=0$ is trivial – use simple
reflection), the only way in which Deift and Park are able to show that
$u(x,t)$ is a classical solution is by using the Bäcklund transform method of
Bikbaev and Tarasov. It is easy to show that the solution has one spatial
derivative in $L^{2}$, but any attempt to iterate the solution in the usual
way in some suitable Sobolev space, $\mathcal{H}$ say, fails, as the mismatch
between the non-linearity of the equation and the Robin boundary condition
immediately kicks one out of the putative iterating space $\mathcal{H}$. This
means, in particular, if one replaces the specific non-linearity $2|u|^{2}u$
in the equation by any other reasonable and smooth non-linearity, for which
there is no known Bäcklund method, for example $2|u|^{4}u$, one does not know
that the solution of the IBVP is classical. So the open problem here is a PDE
problem, viz., to prove that the IPVP for NLS with a general class of non-
linearities, with Robin boundary conditions and matching, smooth initial data,
is classical in $x>0$, $t>0$. Much more generally, the open problem is to
understand what can be said about smoothness of solutions of IBVP’s for
dispersive equations in any dimension. This problem seems to be largely
unexplored.

######  Problem 6 (numerical solution of integrable systems).

Solutions of the Cauchy problem for linear dispersive equations can be
expressed in terms of Fourier integral operators. For example, the solution
$u(x,t)$ of the free Schrödinger problem

$\displaystyle iu_{t}=u_{xx},\qquad$ $\displaystyle-\infty<x<\infty,\quad
t>0,$ $\displaystyle u(x,0)=f(x),\qquad$ $\displaystyle-\infty<x<\infty$

has the form

$\displaystyle u(x,t)=\frac{1}{\sqrt{2\pi}}\int\hat{f}(z)e^{i(xz+tz^{2})}dz,$
(10)

where $\hat{f}(z)=\frac{1}{\sqrt{2\pi}}\int f(x)e^{-ixz}dz$ is the Fourier
transform of $f$. As $t\to\infty$, we can evaluate $u(x,t)$ using the
classical steepest descent method, to obtain

$\displaystyle u(x,t)\sim
e^{i\pi/4}\hat{f}(z_{0})\frac{e^{ix^{2}/4t}}{(2t)^{1/2}},$ (11)

where $z_{0}=-\frac{x}{2t}$ is the stationary phase point for the phase
function $\theta=xz+tz^{2}$. However, if one wants to know the behavior of
$u(x,t)$ for moderate values $t$ for which (11) is no longer a good
approximation, one must resort to the numerical evaluation of (10). The
problem of determining $u(x,t)$ for moderate values of $t$ is in two parts.
The first part is the “forward” problem and involves the evaluation of the
Fourier transform $f\to\hat{f}$. The second part is the “inverse” problem and
involves the evaluation of (10). The second part is considerably harder than
the first part because of the high oscillations coming from the multiplier
$e^{i(xz+tz^{2})}$. For non-linear integrable systems such as the NLS, KdV,
sine-Gordon, …, equations, the solution $u(x,t)$ of the equation is expressed
in terms of the solution of an associated Riemann–Hilbert problem RHP (see,
for example, P. Deift and X. Zhou [47], and the references therein. S. Manakov
[75] was the first to observe that the solution of integrable systems (in
particular, KdV) could be expressed in terms of a RHP). The associated RHP
depends on a product of functions $r(z)e^{\theta_{x,t}(z)}$ where $r(z)$ (the
“reflection” coefficient) encodes the initial data $u_{0}(x)=u(x,t=0)$ for the
equation, and the multiplier $e^{i\theta_{x,t}(z)}$ depends on the equation
itself. For example, for NLS $\theta_{x,t}(z)=\theta(z)=xz+tz^{2}$ as in the
case of the free Schrödinger equation, and for the KdV equation

$\displaystyle u_{t}+6uu_{x}+u_{xxx}=0,\qquad-\infty<x<\infty,$

$\theta_{x,t}(z)=xz+tz^{3}$, which is the same multiplier as for the linear
KdV equation $u_{t}+u_{xxx}=0$. We may think of the mapping
$r(z)e^{i\theta_{x,t}(z)}\to u(x,t)$ as a generalized Fourier integral
transformation, but as opposed to (10), the transformation is highly
nonlinear. As $t\to\infty$, it is possible to use the Deift–Zhou steepest-
descent method for RHP’s in order obtain an expression for the solution of the
equation at hand that is as precise as (11) in the linear case. For moderate
values of $t$, however, one must again resort to numerical methods. This is a
very challenging problem in numerical analysis. Again the problem of the
solution of the Cauchy problem is in two parts. The first part, the forward
problem, involves evaluating the reflection map222In general, one must also
append $L^{2}$ eigenvalues $\\{\lambda_{j}\\}$ and norming constants $(c_{j})$
and consider the extended map $u_{0}\to(r(z),\\{\lambda_{j}\\},\\{c_{j}\\})$.,
$u_{0}(x)\to r(z)$. This is a spectral/scattering theoretic problem. For
example, in the case of KdV, one must solve the spectral/scattering problem
for the one-dimensional Schrödinger operator
$H_{0}=-\frac{d^{2}}{dx^{2}}+u_{0}(x)$. The second part, the inverse problem,
now involves the numerical solution of a RHP in a situation where the
difficulty is again compounded by the high oscillations in the data
$r(z)e^{it\theta_{x,t}(z)}$. The second part is considerably harder than the
first part, and also considerably harder than the evaluation of (10) in the
case of the linear problem. Techniques to solve the forward problem
numerically are readily available (see, for example, the recent text [107] and
the references therein. In situations in which the forward problem itself
contains a small parameter $\epsilon$, $\epsilon\downarrow 0$, there are still
many unresolved issues – see below). Concerning the second problem, in a very
significant development starting in 2011, S. Olver [82] introduced a
methodology to solve RHP’s numerically and applied this methodology to compute
solutions of the Painlevé II equation, which has its own Riemann–Hilbert
formulation. This methodology turned out to be sufficiently robust to capture
the high oscillations arising from factors such as $e^{i\theta_{x,t}(z)}$, and
in 2012, B. Deconinck, S. Olver and T. Trogdon [108] used the methodology to
solve the Cauchy problem for KdV with remarkable accuracy for remarkably long
periods of time, indeed until such time as the steepest-descent method kicks
in. In the text [107] of Olver and Trogdon mentioned above, the authors
present a framework in which to understand the accuracy of their RHP
methodology, and they also discuss applications to other systems.

Rather than a single, or several, open problems, what we have here is an “open
program”, viz., to apply the Olver/Deconinck–Olver–Trogdon methodology to
solve the many different numerical problems for integrable systems as they
arise. One of the most challenging problems concerns situations mentioned
above where the forward problem contains a small parameter $\epsilon$,
$\epsilon\downarrow 0$ (for KdV, $\epsilon$ is the small dispersion parameter
and for NLS, $\epsilon=\hbar$, Planck’s constant). Such problems have been
studied with great interest for more than 30 years, first by P. Lax and D.
Levermore [71, 72, 73], then by S. Venakides [109, 110], then P. Deift, S.
Venakides and X. Zhou [45], then by P. Miller and his collaborators [68, 78],
and many others. Particularly difficult are situations where the forward
problem is not self-adjoint, as in the case of focusing NLS. It remains an
open problem (see [34]) to solve the forward problem for focusing NLS with
general smooth initial data, and then implement the RHP methodology, to obtain
the solution of the Cauchy problem in the small dispersion limit
$\epsilon=\hbar\downarrow 0$. At the technical level, the difficulty in
solving the forward problem is that the relevant physical quantities only
appear, in the language of Martin Kruskal, “beyond all orders”. For example,
every point $\lambda$ in the numerical range of the scattering operator
$T(\hbar)$ is an eigenvalue of $T(\hbar)$ to all orders, i.e.,
$\|(T(\hbar)-\lambda)u\|=\mathcal{O}(|h|^{k})$ for any $k\geq 1$, for some
$u=u(k)$, $\|u\|=1$ (see [34] for more information). The forward problem for
focusing NLS should, of course, be seen as one more example of the general
open problem of computing the spectrum of non-self-adjoint problems (see,
e.g., N. Trefethen and E. Embree [106], …, and particularly E.B. Davies [33]).

Many problems in modern mathematical physics can be expressed in terms of a
Fredholm determinant. For example (see M. Mehta’s classic text [77]),
$F_{s}=\det(1-K_{s})$ is the probability that there are no eigenvalues in a
gap of size $2s/\pi$ in the bulk scaling limit for GUE. Here $K_{s}$ has
kernel $\frac{\sin s(\lambda-\mu)}{\pi(\lambda-\mu)}$, acting on
$L^{2}(-1,1)$. Many authors, starting with F. Dyson [50] in the 1960’s,
followed by [51] in the 70’s, have contributed to the analysis of $F_{s}$ as
$s\to\infty$ (see P. Deift, A. Its, I. Krasovsky and X. Zhou [35] for a
history of this problem), but in order to evaluate $F_{s}$ at finite values of
$s$, we must again turn to numerical evaluation. In a striking development in
2008, F. Bornemann [11] devised a very general and flexible method to evaluate
Fredholm determinants, and he then applied the method to compute a variety of
determinants arising in modern mathematical physics, including $F_{s}$.
Bornemann’s method has emerged as the state of the art for the computation of
determinants in RMT and integrable systems. Although Fredholm determinants
were introduced more than 100 years, it is surprising that, prior to
Bornemann’s work, no such general method to evaluate Fredholm determinants
numerically and accurately was developed. The operator $K_{s}$ above on
$L^{2}(-1,1)$ is also of great interest in many other fields of pure and
applied mathematics, in particular in signal processing where one needs to
know detailed information about the asymptotic behavior as $s\to\infty$ of the
eigenvalues $\\{\lambda_{j}(s)\\}$ and eigenfunctions $\\{u_{j}(s)\\}$ of
$K_{s}$. Quite remarkably $K_{s}$ commutes with a second-order ordinary
differential operator, the prolate spheroidal operator, and so the
eigenfunctions of $K_{s}$ are given by prolate spheroidal functions. The
asymptotics for the $\lambda_{j}(s)$’s were obtained by Slepian [96] in 1965.
The asymptotics of the prolate spheroidal functions are well-known and
classical. However, the behavior of the eigenfunctions for finite values of
$s$, a problem of great practical interest, presents a serious numerical
challenge (a priori, the problem is ill conditioned), which was taken up only
very recently by A. Osipov, V. Rokhlin and H. Xiao in their remarkable text
[84]. A related open problem concerns the Airy operator $A_{s}$ with kernel
$\frac{A_{i}(\lambda)A^{\prime}_{i}(\mu)-A^{\prime}_{i}(\lambda)A_{i}(\mu)}{\lambda-\mu}$
acting on $L^{2}(s,\infty)$ (here $A_{i}(\lambda)$ is the standard Airy
function). Note that $\det(1-A_{s})$ is the probability distribution for the
largest eigenvalue of a GUE matrix in the edge scaling limit. As in the case
of $K_{s}$, the asymptotic behavior as $s\to\infty$ of the eigenvalues
$\\{\tilde{\lambda}_{j}(s)\\}$ and eigenfunctions $\\{\tilde{u}_{j}(s)\\}$ of
$A_{s}$ are of independent interest. For example, it turns out that
$\tilde{\lambda}_{1}(s)$ controls the behavior of the solution of the KdV
equation in the so-called collisioness shock region. The asymptotics of the
$\tilde{\lambda}_{j}(s)$’s were derived recently by T. Bothner [17],
confirming an earlier conjecture of C. Tracy and H. Widom [102]. It is an open
problem to analyze the behavior of the eigenfunctions $\tilde{u}_{j}(s)$ as
$s\to\infty$, and to compute their behavior for finite $s$. We note that the
$\tilde{u}_{j}(s)$’s also solve a second-order ODE (see C. Tracy and H. Widom
[102]).

Finally, in 2005 (see [34]), the idea was raised of a “Painlevé project”,
which would provide a forum to assemble all current information on the
algebraic, analytical and numerical properties of the Painlevé equations. Such
a project was indeed launched in the form of a website maintained at NIST, but
unfortunately the project has not yet met with much success. Given the
importance of the Painlevé equations as the core of modern special function
theory, such a project remains a much needed national resource.

######  Problem 7 (additional problems for integrable systems).

In 2002, P. Deift and X. Zhou [46] analyzed the behavior of solutions of the
Cauchy problem for perturbations of the defocusing NLS equation

$\displaystyle\begin{split}&iu_{t}+u_{xx}-2|u|^{2}u-\epsilon|u|^{\ell}u=0,\\\
&u(x,t=0)=f(x)\to 0\qquad\text{as}\quad|x|\to\infty.\end{split}$ (12)

For $\ell>0$ fixed (for technical reasons they required $\ell>7/2:\ell>2$
should be enough) they showed, quite surprisingly, that the perturbed equation
(12) remained integrable for $\epsilon>0$ and suitably small, and they showed
that as $t\to\infty$ the solution $u(x,t)=u(x,t;\epsilon)$ behaved like
solutions of NLS (i.e., $\epsilon=0$). An open problem of great and basic
interest is to analyze solutions of the Cauchy problem for perturbations of
the focusing NLS equation

$\displaystyle\begin{split}&i\tilde{u}_{t}+\tilde{u}_{xx}+2|\tilde{u}|^{2}u+\epsilon|\tilde{u}|^{\ell}\tilde{u}=0,\\\
&\tilde{u}(x,t=0)=\tilde{f}(x)\to
0\qquad\text{as}\quad|x|\to\infty.\end{split}$ (13)

Of particular interest is the case where $\tilde{f}(x)$ corresponds to a
2-soliton solution for (unperturbed) focusing NLS. The principal new
difficulty in analyzing (13) over (12), is that for solutions $u(x,t)$ of (12)
we have decay in the sup norm

$\displaystyle\sup_{x}|u(x,t)|\leq c/t^{1/2}\to 0\qquad\text{as}\quad
t\to\infty.$

As $\ell>2$, this makes the perturbation term $\epsilon|u|^{\ell}n$ in (12)
small compared to $2|u|^{2}u$ as $t\to\infty$. Such a decay estimate, however,
cannot be true in general for (13) because of the presence of solitons.

In 1991, S. Venakides, P. Deift and R. Oba [111] considered the so-called Toda
shock problem in which a system of particles $x_{k}$, $k\geq 0$, interact with
exponential forces (cf. Problem 5)

$\displaystyle\ddot{x}_{k}=e^{(x_{k-1}-x_{k})}-e^{(x_{k}-x_{k+1})},\qquad
k\geq 1$ (14)

subject to “shock” initial conditions

$\displaystyle x_{k}=k,\qquad\dot{x}_{k}(0)=0,\quad k\geq 1,\qquad
x_{0}(t)=2at,\quad a>0.$

The authors showed that for $a>a_{\text{crit}}=1$, in a frame moving with the
driving particle $x_{0}(t)$, the system executes a binary periodic motion as
$t\to\infty$. Said differently, as $t\to\infty$, the system settles into a
1-gap solution of the doubly infinite periodic Toda lattice. The proof of this
result, however, was not fully rigorous. There are two open problems here. The
first is to use Riemann–Hilbert/steepest-descent methods to prove the results
in Venakides et al. rigorously. The second problem concerns the fact that if
we replace the exponential forces $e^{x}$ by a general force $F(x)$,

$\displaystyle\ddot{x}_{k}=F(x_{k-1}-x_{k})-F(x_{k}-x_{k+1}),\qquad k\geq 1$
(15)

then numerical simulations show that the solution of the shock problem for
this system behaves in the same way as $t\to\infty$ as the Toda system,
provided (15) has a 2-periodic solution

$\displaystyle x_{k}(t+T)=x_{k}(t),\quad T>0,\qquad x_{k+2}(t)=x_{k}(t)$

for $k\in\mathbb{Z}$. Second open problem: Analyze the general shock problem
(15) for suitably small perturbations $F(x)$ of $e^{x}$.

This paper has been about random matrix theory and integrable systems.
Hovering in the background, there is the notion of “integrable probability
theory”. What many people mean by this terminology is that the stochastic
problem at hand has special algebraic property which make it possible to solve
the system explicitly and ultimately express the key statistics for the
problem (e.g., gap probabilities, nearest-neighbor distributions, …) in terms
of solutions of classically integrable Hamiltonian dynamical systems (e.g.,
Painlevé functions). While this interpretation of the terminology certainly
describes accurately the way things work, it ultimately interprets integrable
probability theory in terms of a completely different mathematical field,
viz., dynamical systems. One would like, rather, a definition which is
intrinsic to probability theory. The salient feature of integrable systems
with Hamiltonian $H=H(\xi)$ and $n$ degrees of freedom is that there is a
symplectic change of variables $\xi\to(x,y)$ in which the flow generated by
$H$ is displayed as $n$ independent and elementary motions

$\displaystyle\dot{x}_{k}=y_{k},\qquad\dot{y}_{k}=0,\qquad 1\leq k\leq n,$

yielding

$\displaystyle x_{k}(t)=x_{k}(0)+ty_{k}(0),\qquad y_{k}(t)=y_{k}(0),\qquad
1\leq k\leq n.$

Now coin flips are certainly elementary probabilistic systems and so a
natural, reasonable and intrinsic probabilistic analog of the dynamical notion
of integrability is that there is a change of variables in which key
statistics for the stochastic problem at hand become just a product of
independent coin flips. And the fact of the matter is that this is precisely
the situation for many, if not most, of the known “integrable probabilistic
systems”. For example, for the GUE gap probability $F_{s}=\det(1-K_{s})$
discussed above, the operator $K_{s}$ is trace-class and satisfies the bounds
$0<K_{s}<1$. Hence, $F_{s}$ can be expressed as the product of the eigenvalues
$\\{\lambda_{k}(s)\\}$ of $K_{s}$,

$\displaystyle F_{s}=\prod^{\infty}_{k=1}(1-\lambda_{k}(s)).$ (16)

As $0<\lambda_{k}(s)<1$, (16) displays the gap probability explicitly as a
product of independent flips. Indeed we can imagine that we have a box and a
bag of balls. With probability $\lambda_{1}(s)$ the first ball is placed in
the box. With probability $1-\lambda_{1}(s)$, it is discarded. With
probability $\lambda_{2}(s)$, the second ball is placed in the box, and so on.

Similar considerations apply to all determinantal particle systems, with the
correlation kernel playing a role similar to the role of a generating function
in integrable Hamiltonian theory. This intrinsic definition of an “integrable
probabilistic system”, is just a post facto observation. I don’t know what to
make of it and whether it will be useful, and I am just putting it up there
for the readers’ consideration.

Finally I would like to raise a problem that has been on my mind for many
years. About 120 years ago, in seminal work, A. Sommerfeld [98] showed how to
solve explicitly the diffraction problem for scalar waves at fixed frequency
off a semi-infinite slit $\\{x>0\\}$ in the plane. In the 1930’s P. Morse and
P. Rubenstein [80] showed how to solve explicitly diffraction from a single
finite slit in the plane. However the rigorous solution of the problem of
diffraction from two slits remains open. This problem is of great mathematical
and physical importance, describing as it does Young’s famous double slit
experiment in optics, and the double slit experiment in quantum mechanics. It
is possible that Fokas’ methods (see Problem 5) may be applied to this
problem, giving rise to a Riemann–Hilbert solution from which the relevant
asymptotics of the solution would easily follow. This is a problem of the
first order, whose rigorous and explicit solution would generate great
interest in the mathematical and physics communities.

We end with an extended paraphrase of [34]:

> “There are many other areas, closely related to the problems in the above
> list, where much progress has been made in recent years, and where much
> remains to be done. These include: partition functions for random matrix
> models and their connections to mappings on Riemann surfaces, the
> representation theory of “large” groups such as the infinite symmetric group
> and the infinite unitary group, Macdonald processes, singular value
> statistics of matrix products, statistics for normal matrix ensembles, the
> six vertex model, random matrix models with an external source, non-
> intersecting random paths and the tacnode phenomenon, SVD statistics,
> asymptotics for Toeplitz and Hankel matrices, singular values of $n$
> products of $m\times m$ random matrices as $n$, $m\to\infty$, and many
> others. But as I said in the beginning, my list is not complete or
> definitive.”

## Acknowledgements

This work was supported in part by NSF Grant DMS-1300965. The author would
also like to thank Ivan Corwin for his help with Problem 3, and Jinho Baik and
Ivan Corwin for their unstinting help in preparing this article for
publication. The author is also very grateful for the opportunity to present
this list of problems which continue to intrigue and challenge him.

## References

  * [1] Ajanki O., Erdős L., Krüger T., Stability of the matrix Dyson equation and random matrices with correlations, arXiv:1604.08188.
  * [2] Amir G., Corwin I., Quastel J., Probability distribution of the free energy of the continuum directed random polymer in $1+1$ dimensions, Comm. Pure Appl. Math. 64 (2011), 466–537, arXiv:1003.0443.
  * [3] Baik J., Deift P., Johansson K., On the distribution of the length of the longest increasing subsequence of random permutations, J. Amer. Math. Soc. 12 (1999), 1119–1178, math.CO/9810105.
  * [4] Bakhtin Yu., Correll J., A neural computation model for decision-making times, J. Math. Psych. 56 (2012), 333–340.
  * [5] Bertini L., Cancrini N., The stochastic heat equation: Feynman–Kac formula and intermittence, J. Stat. Phys. 78 (1995), 1377–1401.
  * [6] Bikbaev R.F., Tarasov V.O., Initial-boundary value problem for the nonlinear Schrödinger equation, J. Phys. A: Math. Gen. 24 (1991), 2507–2516.
  * [7] Binder I., Damanik D., Goldstein M., Lukic M., Almost periodicity in time of solutions of the KdV equation, arXiv:1509.07373.
  * [8] Binder I., Damanik D., Lukic M., VandenBoom T., Almost periodicity in time of solutions of the Toda lattice, arXiv:1603.04905.
  * [9] Bleher P., Its A., Double scaling limit in the random matrix model: the Riemann–Hilbert approach, Comm. Pure Appl. Math. 56 (2003), 433–516, math-ph/0201003.
  * [10] Bohigas O., Pato M.P., Randomly incomplete spectra and intermediate statistics, Phys. Rev. E 74 (2006), 036212, 6 pages, cond-mat/0608165.
  * [11] Bornemann F., On the numerical evaluation of Fredholm determinants, Math. Comp. 79 (2010), 871–915, arXiv:0804.2543.
  * [12] Borodin A., Corwin I., Macdonald processes, Probab. Theory Related Fields 158 (2014), 225–400, arXiv:1111.4408.
  * [13] Borodin A., Okounkov A., Olshanski G., Asymptotics of Plancherel measures for symmetric groups, J. Amer. Math. Soc. 13 (2000), 481–515, math.CO/9905032.
  * [14] Borodin A., Olshanski G., The ASEP and determinantal point processes, arXiv:1608.01564.
  * [15] Borodin A., Petrov L., Integrable probability: from representation theory to Macdonald processes, Probab. Surv. 11 (2014), 1–58, arXiv:1310.8007.
  * [16] Borot G., Eynard B., Majumdar S.N., Nadal C., Large deviations of the maximal eigenvalue of random matrices, J. Stat. Mech. Theory Exp. 2011 (2011), P11024, 56 pages, arXiv:1009.1945.
  * [17] Bothner T., From gap probabilities in random matrix theory to eigenvalue expansions, J. Phys. A: Math. Theor. 49 (2016), 075204, 77 pages, arXiv:1509.07159.
  * [18] Bothner T., Deift P., Its A., Krasovsky I., On the asymptotic behavior of a log gas in the bulk scaling limit in the presence of a varying external potential I, Comm. Math. Phys. 337 (2015), 1397–1463, arXiv:1407.2910.
  * [19] Bourgade P., Erdös L., Yau H.-T., Edge universality of beta ensembles, Comm. Math. Phys. 332 (2014), 261–353, arXiv:1306.5728.
  * [20] Bourgade P., Erdös L., Yau H.-T., Yin J., Universality for a class of random band matrices, arXiv:1602.02312.
  * [21] Bourgade P., Yau H.-T., The eigenvector moment flow and local quantum unique ergodicity, Comm. Math. Phys. 350 (2017), 231–278, arXiv:1312.1301.
  * [22] Calabrese P., Le Doussal P., Rosso A., Free-energy distribution of the directed polymer at high temperature, Europhys. Lett. 90 (2010), 20002, 6 pages, arXiv:1002.4560.
  * [23] Charlier C., Claeys T., Thinning and conditioning of the circular unitary ensemble, arXiv:1604.08399.
  * [24] Claeys T., Fahs B., Random matrices with merging singularities and the Painlevé V equation, SIGMA 12 (2016), 031, 44 pages, arXiv:1508.06734.
  * [25] Claeys T., Kuijlaars A.B.J., Universality of the double scaling limit in random matrix models, Comm. Pure Appl. Math. 59 (2006), 1573–1603, math-ph/0501074.
  * [26] Corwin I., The Kardar–Parisi–Zhang equation and universality class, Random Matrices Theory Appl. 1 (2012), 1130001, 76 pages, arXiv:1106.1596.
  * [27] Corwin I., Kardar–Parisi–Zhang universality, Notices Amer. Math. Soc. 63 (2016), 230–239, arXiv:1606.06602.
  * [28] Corwin I., Hammond A., KPZ line ensemble, Probab. Theory Related Fields 166 (2016), 67–185, arXiv:1312.2600.
  * [29] Corwin I., Liu Z., Wang D., Fluctuations of TASEP and LPP with general initial data, Ann. Appl. Probab. 26 (2016), 2030–2082, arXiv:1412.5087.
  * [30] Corwin I., Quastel J., Remenik D., Renormalization fixed point of the KPZ universality class, J. Stat. Phys. 160 (2015), 815–834, arXiv:1103.3422.
  * [31] Damanik D., Goldstein M., On the inverse spectral problem for the quasi-periodic Schrödinger equation, Publ. Math. Inst. Hautes Études Sci. 119 (2014), 217–401, arXiv:1209.4331.
  * [32] Damanik D., Goldstein M., On the existence and uniqueness of global solutions for the KdV equation with quasi-periodic initial data, J. Amer. Math. Soc. 29 (2016), 825–856, arXiv:1212.2674.
  * [33] Davies E.B., Linear operators and their spectra, Cambridge Studies in Advanced Mathematics, Vol. 106, Cambridge University Press, Cambridge, 2007.
  * [34] Deift P., Some open problems in random matrix theory and the theory of integrable systems, in Integrable Systems and Random Matrices, Contemp. Math., Vol. 458, Amer. Math. Soc., Providence, RI, 2008, 419–430, arXiv:0712.0849.
  * [35] Deift P., Its A., Krasovsky I., Zhou X., The Widom–Dyson constant for the gap probability in random matrix theory, J. Comput. Appl. Math. 202 (2007), 26–47, math.FA/0601535.
  * [36] Deift P., Kriecherbauer T., McLaughlin K.T.-R., Venakides S., Zhou X., Uniform asymptotics for polynomials orthogonal with respect to varying exponential weights and applications to universality questions in random matrix theory, Comm. Pure Appl. Math. 52 (1999), 1335–1425.
  * [37] Deift P., Kriecherbauer T., Venakides S., Forced lattice vibrations. Part I, Comm. Pure Appl. Math. 48 (1995), 1187–1249.
  * [38] Deift P., Li L.C., Nanda T., Tomei C., The Toda flow on a generic orbit is integrable, Comm. Pure Appl. Math. 39 (1986), 183–232.
  * [39] Deift P., Li L.C., Tomei C., Matrix factorizations and integrable systems, Comm. Pure Appl. Math. 42 (1989), 443–521.
  * [40] Deift P., Menon G., Olver S., Trogdon T., Universality in numerical computations with random data, Proc. Natl. Acad. Sci. USA 111 (2014), 14973–14978, arXiv:1407.3829.
  * [41] Deift P., Park J., Long-time asymptotics for solutions of the NLS equation with a delta potential and even initial data, Int. Math. Res. Not. 2011 (2011), 5505–5624, arXiv:1006.4909.
  * [42] Deift P., Trogdon T., Universality for the Toda algorithm to compute the eigenvalues of a random matrix, arXiv:1604.07384.
  * [43] Deift P., Trogdon T., Universality for eigenvalue algorithms on sample covariance matrices, arXiv:1701.01896.
  * [44] Deift P., Trogdon T., Menon G., On the condition number of the critically-scaled Laguerre unitary ensemble, Discrete Contin. Dyn. Syst. 36 (2016), 4287–4347, arXiv:1507.00750.
  * [45] Deift P., Venakides S., Zhou X., New results in small dispersion KdV by an extension of the steepest descent method for Riemann–Hilbert problems, Int. Math. Res. Not. 1997 (1997), 286–299.
  * [46] Deift P., Zhou X., Perturbation theory for infinite-dimensional integrable systems on the line. A case study, Acta Math. 188 (2002), 163–262,.
  * [47] Deift P., Zhou X., Long-time asymptotics for solutions of the NLS equation with initial data in a weighted Sobolev space, Comm. Pure Appl. Math. 56 (2003), 1029–1077, math.AP/0206222.
  * [48] Dotsenko V., Bethe ansatz derivation of the Tracy–Widom distribution for one-dimensional directed polymers, Europhys. Lett. 90 (2010), 20003, 5 pages.
  * [49] Dyson F.J., A Brownian-motion model for the eigenvalues of a random matrix, J. Math. Phys. 3 (1962), 1191–1198.
  * [50] Dyson F.J., Statistical theory of the energy levels of complex systems. II, J. Math. Phys. 3 (1962), 157–165.
  * [51] Dyson F.J., Fredholm determinants and inverse scattering problems, Comm. Math. Phys. 47 (1976), 171–183.
  * [52] Erdős L., Universality for random matrices and log-gases, in Current Developments in Mathematics 2012, Int. Press, Somerville, MA, 2013, 59–132, arXiv:1212.0839.
  * [53] Erdős L., Knowles A., Yau H.-T., Yin J., Delocalization and diffusion profile for random band matrices, Comm. Math. Phys. 323 (2013), 367–416, arXiv:1205.5669.
  * [54] Erdős L., Péché S., Ramírez J.A., Schlein B., Yau H.-T., Bulk universality for Wigner matrices, Comm. Pure Appl. Math. 63 (2010), 895–925, arXiv:0905.4176.
  * [55] Erdős L., Yau H.-T., Yin J., Rigidity of eigenvalues of generalized Wigner matrices, Adv. Math. 229 (2012), 1435–1515, arXiv:1007.4652.
  * [56] Fokas A.S., A unified transform method for solving linear and certain nonlinear PDEs, Proc. Roy. Soc. London Ser. A 453 (1997), 1411–1443.
  * [57] Fokas A.S., Lenells J., The unified method: I. Nonlinearizable problems on the half-line, J. Phys. A: Math. Theor. 45 (2012), 195201, 38 pages, arXiv:1109.4935.
  * [58] Fröhlich J., Spencer T., Absence of diffusion in the Anderson tight binding model for large disorder or low energy, Comm. Math. Phys. 88 (1983), 151–184.
  * [59] Fyodorov Y.V., Mirlin A.D., Scaling properties of localization in random band matrices: a $\sigma$-model approach, Phys. Rev. Lett. 67 (1991), 2405–2409.
  * [60] Grava T., Its A., Kapaev A., Mezzadri F., On the Tracy–Widomβ distribution for $\beta=6$, SIGMA 12 (2016), 105, 26 pages, arXiv:1607.01351.
  * [61] Hairer M., Solving the KPZ equation, Ann. of Math. 178 (2013), 559–664, arXiv:1109.6811.
  * [62] Its A., Shepelsky D., Initial boundary value problem for the focusing nonlinear Schrödinger equation with Robin boundary condition: half-line approach, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 469 (2013), 20120199, 14 pages, arXiv:1201.5948.
  * [63] Johansson K., Shape fluctuations and random matrices, Comm. Math. Phys. 209 (2000), 437–476, math.CO/9903134.
  * [64] Johansson K., Transversal fluctuations for increasing subsequences on the plane, Probab. Theory Related Fields 116 (2000), 445–456, math.PR/9910146.
  * [65] Johansson K., Universality of the local spacing distribution in certain ensembles of Hermitian Wigner matrices, Comm. Math. Phys. 215 (2001), 683–705, math-ph/0006020.
  * [66] Johansson K., Discrete polynuclear growth and determinantal processes, Comm. Math. Phys. 242 (2003), 277–329, math.PR/0206208.
  * [67] Johansson K., Two time distribution in Brownian directed percolation, Comm. Math. Phys. 351 (2017), 441–492, arXiv:1502.00941.
  * [68] Kamvissis S., McLaughlin K.D.T.-R., Miller P.D., Semiclassical soliton ensembles for the focusing nonlinear Schrödinger equation, Annals of Mathematics Studies, Vol. 154, Princeton University Press, Princeton, NJ, 2003\.
  * [69] Kardar M., Parisi G., Zhang Y.-C., Dynamic scaling of growing interfaces, Phys. Rev. Lett. 56 (1986), 889–892.
  * [70] Kriecherbauer T., Shcherbina M., Fluctuations of eigenvalues of matrix models and their applications, arXiv:1003.6121.
  * [71] Lax P.D., Levermore C.D., The small dispersion limit of the Korteweg–de Vries equation. I, Comm. Pure Appl. Math. 36 (1983), 253–290.
  * [72] Lax P.D., Levermore C.D., The small dispersion limit of the Korteweg–de Vries equation. II, Comm. Pure Appl. Math. 36 (1983), 571–593.
  * [73] Lax P.D., Levermore C.D., The small dispersion limit of the Korteweg–de Vries equation. III, Comm. Pure Appl. Math. 36 (1983), 809–829.
  * [74] Li X.H., Menon G., Numerical solution of Dyson Brownian motion and a sampling scheme for invariant matrix ensembles, J. Stat. Phys. 153 (2013), 801–812, arXiv:1306.1179.
  * [75] Manakov S.V., Nonlinear Fraunhofer diffraction, Sov. Phys. JETP 38 (1974), 693–696.
  * [76] McKean H.P., Trubowitz E., Hill’s operator and hyperelliptic function theory in the presence of infinitely many branch points, Comm. Pure Appl. Math. 29 (1976), 143–226.
  * [77] Mehta M.L., Random matrices, Pure and Applied Mathematics (Amsterdam), Vol. 142, 3rd ed., Elsevier/Academic Press, Amsterdam, 2004.
  * [78] Miller P.D., Xu Z., On the zero-dispersion limit of the Benjamin–Ono Cauchy problem for positive initial data, Comm. Pure Appl. Math. 64 (2011), 205–270, arXiv:1002.3178.
  * [79] Minami N., Local fluctuation of the spectrum of a multidimensional Anderson tight binding model, Comm. Math. Phys. 177 (1996), 709–725.
  * [80] Morse P.M., Rubenstein P.J., The diffraction of waves by ribbons and by slits, Phys. Rev. 54 (1938), 895–898.
  * [81] Moser J., Three integrable Hamiltonian systems connected with isospectral deformations, Adv. Math. 16 (1975), 197–220.
  * [82] Olver S., Numerical solution of Riemann–Hilbert problems: Painlevé II, Found. Comput. Math. 11 (2011), 153–179.
  * [83] Olver S., Nadakuditi R.R., Trogdon T., Sampling unitary ensembles, Random Matrices Theory Appl. 4 (2015), 1550002, 22 pages, arXiv:1404.0071.
  * [84] Osipov A., Rokhlin V., Xiao H., Prolate spheroidal wave functions of order zero. Mathematical tools for bandlimited approximation, Applied Mathematical Sciences, Vol. 187, Springer, New York, 2013.
  * [85] Pfrang C.W., Deift P., Menon G., How long does it take to compute the eigenvalues of a random symmetric matrix?, in Random Matrix Theory, Interacting Particle Systems, and Integrable Systems, Math. Sci. Res. Inst. Publ., Vol. 65, Cambridge University Press, New York, 2014, 411–442, arXiv:1203.4635.
  * [86] Prähofer M., Spohn H., Scale invariance of the PNG droplet and the Airy process, J. Stat. Phys. 108 (2002), 1071–1106, math.PR/0105240.
  * [87] Ramírez J.A., Rider B., Virág B., Beta ensembles, stochastic Airy spectrum, and a diffusion, J. Amer. Math. Soc. 24 (2011), 919–944, math.PR/0607331.
  * [88] Rumanov I., Beta ensembles, quantum Painlevé equations and isomonodromy systems, in Algebraic and Analytic Aspects of Integrable Systems and Painlevé Equations, Contemp. Math., Vol. 651, Amer. Math. Soc., Providence, RI, 2015, 125–155, arXiv:1408.3847.
  * [89] Rumanov I., Painlevé representation of Tracy–Widomβ distribution for $\beta=6$, Comm. Math. Phys. 342 (2016), 843–868, arXiv:1408.3779.
  * [90] Sasamoto T., Spohn H., One-dimensional KPZ equation: an exact solution and its universality, Phys. Rev. Lett. 104 (2010), 230602, 4 pages, arXiv:1002.1883.
  * [91] Schenker J., Eigenvector localization for random band matrices with power law band width, Comm. Math. Phys. 290 (2009), 1065–1097, arXiv:0809.4405.
  * [92] Shcherbina M., Orthogonal and symplectic matrix models: universality and other properties, Comm. Math. Phys. 307 (2011), 761–790, arXiv:1004.2765.
  * [93] Shcherbina T., On the second mixed moment of the characteristic polynomials of 1D band matrices, Comm. Math. Phys. 328 (2014), 45–82, arXiv:1209.3385.
  * [94] Shcherbina T., Universality of the second mixed moment of the characteristic polynomials of the 1D band matrices: real symmetric case, J. Math. Phys. 56 (2015), 063303, 23 pages, arXiv:1410.3084.
  * [95] Sinai Ya.G., Soshnikov A.B., A refinement of Wigner’s semicircle law in a neighborhood of the spectrum edge for random symmetric matrices, Funct. Anal. Appl. 32 (1998), 114–131.
  * [96] Slepian D., Some asymptotic expansions for prolate spheroidal wave functions, J. Math. and Phys. 44 (1965), 99–140.
  * [97] Sodin S., The spectral edge of some random band matrices, Ann. of Math. 172 (2010), 2223–2251, arXiv:0906.4047.
  * [98] Sommerfeld A., Mathematische Theorie der Diffraction, Math. Ann. 47 (1896), 317–374.
  * [99] Soshnikov A., Universality at the edge of the spectrum in Wigner random matrices, Comm. Math. Phys. 207 (1999), 697–733, math-ph/9907013.
  * [100] Symes W.W., The $QR$ algorithm and scattering for the finite nonperiodic Toda lattice, Phys. D 4 (1982), 275–280.
  * [101] Tao T., Vu V., Random matrices: universality of local eigenvalue statistics, Acta Math. 206 (2011), 127–204, arXiv:0906.0510.
  * [102] Tracy C.A., Widom H., Level-spacing distributions and the Airy kernel, Comm. Math. Phys. 159 (1994), 151–174, hep-th/9211141.
  * [103] Tracy C.A., Widom H., A Fredholm determinant representation in ASEP, J. Stat. Phys. 132 (2008), 291–300, arXiv:0804.1379.
  * [104] Tracy C.A., Widom H., Integral formulas for the asymmetric simple exclusion process, Comm. Math. Phys. 279 (2008), 815–844, Erratum, Comm. Math. Phys. 304 (2011), 875–878, arXiv:0704.2633.
  * [105] Tracy C.A., Widom H., Asymptotics in ASEP with step initial condition, Comm. Math. Phys. 290 (2009), 129–154, arXiv:0807.1713.
  * [106] Trefethen L.N., Embree M., Spectra and pseudospectra. The behavior of nonnormal matrices and operators, Princeton University Press, Princeton, NJ, 2005.
  * [107] Trogdon T., Olver S., Riemann–Hilbert problems, their numerical solution, and the computation of nonlinear special functions, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2016.
  * [108] Trogdon T., Olver S., Deconinck B., Numerical inverse scattering for the Korteweg–de Vries and modified Korteweg–de Vries equations, Phys. D 241 (2012), 1003–1025.
  * [109] Venakides S., The zero dispersion limit of the Korteweg–de Vries equation for initial potentials with nontrivial reflection coefficient, Comm. Pure Appl. Math. 38 (1985), 125–155.
  * [110] Venakides S., The Korteweg–de Vries equation with small dispersion: higher order Lax–Levermore theory, Comm. Pure Appl. Math. 43 (1990), 335–361.
  * [111] Venakides S., Deift P., Oba R., The Toda shock problem, Comm. Pure Appl. Math. 44 (1991), 1171–1242.
  * [112] Witte N.S., Bornemann F., Forrester P.J., Joint distribution of the first and second eigenvalues at the soft edge of unitary ensembles, Nonlinearity 26 (2013), 1799–1822, arXiv:1209.2190.
  * [113] Wu T.T., McCoy B.M., Tracy C.A., Barouch E., Spin-spin correlation functions for the two-dimensional Ising model: exact theory in the scaling region, Phys. Rev. B 13 (1976), 316–374.