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

TL;DR

This paper reviews open problems and recent progress in random matrix theory and integrable systems, highlighting ongoing research challenges and advancements discussed at a 2015 conference and in earlier work.

Contribution

It compiles a list of open problems in the field and reports on progress made since a previous related publication, advancing understanding of key issues.

Findings

Progress on open problems in random matrix theory

Identification of key challenges in integrable systems

Updated research directions from recent conference

Abstract

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 Mathematiques, Montreal, 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]).

Equations

$$u_t+u{u}_x+u_{xxx}=0$$

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

$$V(x)={\gamma }_k{x}^{2k}+{\gamma }_{k-1}{x}^{2k-1}+\cdots +{\gamma }_0,{\gamma }_k>0$$

$$S_n=x_1+x_2+\cdots +x_n\to {\hat{x}}_1+x_2+\cdots +x_n\to {\hat{x}}_1+{\hat{x}}_2+x_3+\cdots +x_n$$

$$\to \cdots \to S_n={\hat{x}}_1+{\hat{x}}_2+\cdots +{\hat{x}}_n$$

$$\bullet \text{(isospectral)}$$

$$\bullet \text{(convergence)}$$

$$X_0=H\text{given, converge to a diagonal matrix X∞ as k→∞},$$

$$\bullet \text{(isospectral) }$$

$$\bullet \text{(convergence)}$$

$$\text{converges to a diagonal matrix X∞,X(t)→X∞ as t→∞}.$$

$${\tau }_{ϵ,n,\mathcal{A},\mathcal{E}}(H)\equiv \frac{T_{ϵ,n,\mathcal{A},\mathcal{E}}(H)-⟨T_{ϵ,n,\mathcal{A},\mathcal{E}}⟩}{{\sigma }_{ϵ,n,\mathcal{A},\mathcal{E}}},$$

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

$$\frac{\log {ϵ}^{-1}}{\log n}\ge \frac{5}{3}+\frac{\tau }{2}.$$

$$\underset{n\to \infty }{\lim }\mathrm{Prob}\left(\frac{T^{(1)}}{C_v^{2/3}{2}^{-2/3}{N}^{2/3}\left(\log {ϵ}^{-1}-\frac{2}{3}\log n\right)}\right)=F_{\beta }^{\text{gap}}(t).$$

$$F_{\beta }^{\text{gap}}(t)=\underset{n\to \infty }{\lim }\mathrm{Prob}\left(\frac{1}{C_v^{2/3}{2}^{-2/3}{n}^{2/3}({\lambda }_n-{\lambda }_{n-1})}\le t\right),$$

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

$$E(t)\equiv \sum _{j=2}^n|X_{ij}(t){|}^2=\sum _{j=1}^n({\lambda }_j-X_{11}(t){)}^2|u_{1j}(t){|}^2,$$

$${\lambda }_j-X_{11}(t)=\sum _{i=1}^{n-1}({\lambda }_j-{\lambda }_i)|u_{1i}(t){|}^2,$$

$$|u_{1j}(t)|=\frac{{\beta }_j{e}^{{\lambda }_{j}t}}{{\left(\sum _{i=1}^n{\beta }_i^2{e}^{2{\lambda }_{i}t}\right)}^{\frac{1}{2}}},$$

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

$$\mathrm{Prob}(R_{n,s})=1+o(1)$$

$$n^{-1/2}\{{\beta }_n,{\beta }_{n-1},{\beta }_{n-2}\}$$

$$2^{-2/3}{n}^{2/3}(b_v-{\lambda }_n,b_v-{\lambda }_{n-1},b_v-{\lambda }_{n-2})$$

$$i{u}_t=u_{xx}\pm 2|u{|}^{2}u$$

$$u(x,0)=f(x),$$

$$u(0,t)=g(t),$$

$${\ddot{x}}_k=e^{x_{k-1}-x_k}-e^{x_k-x_{k+1}},k\ge 1$$

$$x_k(0)={\dot{x}}_k(0)=0,k\ge 1,$$

$$x_0(t)=2at+h(\gamma t).$$

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

$${\gamma }_1>{\gamma }_2>\cdots >{\gamma }_j>\cdots >0,{\gamma }_j\to 0\text{as}j\to \infty$$

$$x_k(t)=c_{1}k+Y_1({\beta }_{1}k+\gamma t),1\ll k\ll t$$

$$x_k=c_{j}k+Y_j({\beta }_{1}k+\gamma t,{\beta }_{2}k+\gamma t,\dots ,{\beta }_{j}k+\gamma t),1\ll k\ll t$$

Full text

[TABLE]

with initial data

[TABLE]

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.

[TABLE]

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$,

[TABLE]

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.

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.

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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].

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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$:

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.

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.,

[TABLE]

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

[TABLE]

as $n\to\infty$.

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

[TABLE]

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

[TABLE]

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.

[TABLE]

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

[TABLE]

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

[TABLE]

and initial/boundary value conditions

[TABLE]

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,

[TABLE]

Moreover, there is a sequence of thresholds

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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.

[TABLE]

has the form

[TABLE]

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

[TABLE]

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

[TABLE]

$\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.

[TABLE]

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

[TABLE]

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

[TABLE]

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)

[TABLE]

subject to “shock” initial conditions

[TABLE]

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)$,

[TABLE]

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

[TABLE]

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

[TABLE]

yielding

[TABLE]

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}$,

[TABLE]

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.