Limits of spiked random matrices II

Introduction

Johnstone (2001) proposed the spiked population model for simple trends in high dimensional data. One takes a data matrix $X$ whose columns are i.i.d. vectors with (population) covariance a fixed rank perturbation of the identity, and studies the behaviour of the largest eigenvalues of the sample covariance matrix $XX^{*}$ when both the dimension and the size of the sample are large. Baik, Ben Arous and Péché (2005) (hereafter BBP) discovered a very interesting phase transition phenomenon in the complex Gaussian setting. Small spikes do not affect the asymptotic behaviour of the top eigenvalues, which display the usual Tracy–Widom fluctuations around the upper edge of the Marchenko–Pastur law; large spikes, however, lead to outliers with Gaussian fluctuations. New structure emerges near the transition point with near-critical spikes deforming the soft edge limit. Understanding this transition regime in the real case remained open for some time. There is a parallel development for fixed rank additive perturbations of Wigner matrices.

In Bloemendal and Virág (2013) (hereafter Part I), we considered rank one spiked real/complex/quaternion Wishart matrices and additive rank one perturbations of the Gaussian orthogonal, unitary and symplectic ensembles. Our approach is based on the continuum operator limit at the general beta soft edge developed in Ramírez, Rider and Virág (2011) (hereafter RRV). We introduced general $\beta$ analogues of the rank one spiked models, modifying the tridiagonal ensembles of Dumitriu and Edelman (2002) and extended the RRV technology to describe the soft-edge scaling limit in terms of the stochastic Airy operator

We went on to characterize the limit laws in terms of the diffusion from RRV and in terms of an associated second-order linear parabolic PDE. We further showed that at $\beta=2,4$ the PDE is related to known Painlevé II representations originating in Baik and Rains (2000) and gave new proofs of these, finally recovering those of the undeformed Tracy–Widom laws.

Even the existence of limiting distributions in the critical regime was in general new for $\beta\neq 2$, though see the prior work of Wang (2008) on the rank one $\beta=4$ case at $w=0$, as well as the subsequent work of Mo (2012) offering a more standard treatment of the rank one $\beta=1$ case. Forrester (2013) comments on all three works and gives an alternative interpretation and construction of our general $\beta$ rank one spiked model.

Here, we deal with $r$ “spikes”, or general bounded-rank perturbations of Gaussian and Wishart matrices. To do so, we introduce a new “canonical form for perturbations in a fixed subspace”, a $(2r+1)$-diagonal band form that has a purely algebraic interpretation. It generalizes the Dumitriu–Edelman forms and is able to handle rank $r$ perturbations. We then develop a generalization of the methods of RRV and Part I to a matrix-valued setting: block tridiagonal matrices converge to a half-line Schrödinger operator with matrix-valued potential, the spikes once again appearing in the boundary condition. We treat the real, complex and quaternion ($\beta=1,2,4$) cases simultaneously. Once again, even the existence of a near-critical soft-edge limit is new off $\beta=2$. Unlike in Part I, however, we do not define a general $\beta$ version of either matrix model, nor of the limiting operator; in Section 2, we will see that the higher rank versions of these objects do not readily admit a $\beta$-generalization.

Dyson’s Brownian motion makes a surprise appearance, providing nice SDE and PDE characterizations of the limit laws—new $r$ parameter deformations of Tracy–Widom($\beta$)—in which $\beta$ reappears as a simple parameter. The derivation makes use of the matrix-valued version of classical Sturm oscillation theory and the Riccati transformation. In a short final section, we report on preliminary evidence that at $\beta=2$ the PDE can be connected with a Painlevé II representation of Baik (2006) for these distributions (which appeared originally in BBP in the form of Fredholm determinants).

We highlight two more features of our approach beyond the novelty of bypassing formulas for joint eigenvalue densities and handling $\beta=1,2,4$ together. First, we treat the perturbation as a parameter. By this, we mean that all perturbations in a fixed subspace are considered jointly (on the same probability space); this picture is carried through to the limit, which is therefore a family of point processes parameterized by an $r\times r$ matrix. Second, we allow more general scalings than those considered in BBP. Most importantly, in the Wishart case we do not require the two dimensional parameters $n,p$ to have a positive limiting ratio but rather allow them to tend to infinity together arbitrarily.

To state our results, we introduce some objects and notation that will be used throughout the paper.

The result is a $(2r+1)$-diagonal matrix with positive outer diagonals. For Gaussian and null Wishart ensembles, the change of basis interacts well with the Gaussian structure; this observation goes back to Trotter (1984) in the $r=1$ case. In the GO/U/SE case, we take $v_{1},\ldots,v_{r}$ to be the initial coordinate basis vectors, while in the Wishart case we use the initial rows of the data matrix $X$. As in Part I, the key observation is then that the perturbations commute with the change of basis.

For the (unperturbed) Gaussian ensembles, the band form looks like

In Section 4, we apply this result to the band forms just described, proving a process central limit theorem for the potential and verifying the required tightness assumptions. The limiting operator turns out to be a multidimensional version of the stochastic Airy operator, which we now describe.

where $B^{\prime}_{x}$ is “standard matrix white noise”, the derivative of a standard matrix Brownian motion, and $rx$ is scalar. (Here, again $\beta$ is restricted to the classical values, as the noise term lacks a straightforward $\beta$-generalization.) The potential is thus the derivative of a continuous matrix-valued function; rigorous definitions will appear in Section 3 in a more general setting.

For now it is enough to know that, together with a general self-adjoint boundary condition

For concreteness, we record that the eigenvalues $\Lambda_{0}\leq\Lambda_{1}\leq\dots$ and corresponding eigenfunctions $f_{0},f_{1},\dots$ of $\mathcal{H}_{\beta,W}$ are given, respectively, by the minimum and any minimizer in the recursive variational problem

Here, candidates $f$ are only considered if the first integral and boundary term are finite; the stochastic integral can then be defined pathwise via integration by parts. The eigenvalues and eigenfunctions are thus jointly defined random processes indexed over $W$.

We note one important property of the eigenvalue processes, namely the pathwise monotonicity of $\Lambda_{k}$ in $W$ with respect to the usual matrix partial order. This is immediate from the variational characterization and the fact that the objective functional is monotone in $W$. (For the higher eigenvalues, it is most apparent from the standard min–max formulation of the variational problem.)

then, jointly for $k=1,2,\dots$ in the sense of finite-dimensional distributions,

where $\Lambda_{0}\leq\Lambda_{1}\leq\dots$ are the eigenvalues of $\mathcal{H}_{\beta,W}$. Convergence holds jointly over $\{P_{n}\},W$ satisfying the condition.

then, jointly for $k=1,2,\dots$ in the sense of finite-dimensional distributions,

where $\Lambda_{0}\leq\Lambda_{1}\leq\dots$ are the eigenvalues of $\mathcal{H}_{\beta,W}$. Convergence holds jointly over $\{\Sigma_{n,p}\},W$ satisfying the condition.

for $k=0,1,\ldots.$ Write simply $F_{\beta}=F_{\beta}^{0}$ for the ground state distribution (limiting largest eigenvalue law). Once again, the generalization from Part I is not straightforward. The proofs are contained in Section 5.

Let $\operatorname{\mathbf{P}}_{x_{0},(w_{1},\ldots,w_{r})}$ be the measure on paths $(p_{1},\ldots,p_{r}):[x_{0},\infty)\to(-\infty,\infty]^{r}$ determined by the coupled diffusions

with initial conditions $p_{i}(x_{0})=w_{i}$ and entering into $\{p_{1}<\cdots<p_{r}\}$, where $b_{1},\ldots,b_{r}$ are independent standard Brownian motions; particles $p_{i}$ may explode to $-\infty$ in finite time whereupon they are restarted at $+\infty$. Then

We describe the diffusion more carefully in Section 5, asserting that it determines a law on paths valued in an appropriate space. Probabilistic arguments lead to the following reformulation in terms of its generator.

Furthermore, $F_{\beta}$ is “continuous to the boundary” as one or several $w_{i}\to+\infty$. For subsequent eigenvalue laws $F_{\beta}^{k}(x;w_{1},\ldots,w_{r})$, (8) is replaced with the recursive boundary condition

At $\beta=2$, these distributions were obtained in BBP in the form of Fredholm determinants of finite-rank perturbations of the Airy kernel. Baik (2006) derived Painlevé II formulas, and by a symbolic computation with a computer algebra system we were able to verify that the latter satisfy the PDE (6) for $r=2,3,4,5$; details are described in Section 6. A pencil-and-paper proof for all $r$ was found since the initial posting [Bloemendal and Baik (2013)].

We make two final remarks. From the finite $n$ matrix models it is clear that the “rank $r$ deformed” limiting distributions $F_{\beta,r}(x;w_{1},\ldots,w_{r})$ reduce to those for a lower rank $r_{0}<r$ in the following way:

Unfortunately, this reduction relation is not readily apparent from any of our characterizations (operator, SDE or PDE).

Lastly, the SDE and PDE characterizations seem to make sense for all $\beta>0$ (although one has to be careful for $\beta<1$). It would be interesting to find natural “general $\beta$ multi-spiked models” at finite $n$, interpolating between those studied here at $\beta=1,2,4$ and generalizing those introduced in Part I for $r=1$. At $\beta=2$, perhaps one could discover a relationship with formulas of Baik and Wang (2013).

A canonical form for perturbation in a fixed subspace

In Part I, we observed that the tridiagonal models of Gaussian and Wishart matrices were amenable to rank one perturbation. In this section, we introduce a banded (also block tridiagonal) generalization amenable to higher-rank perturbation. We first describe it as a natural object of pure linear algebra; we then show how it interacts with the structure of Gaussian and Wishart random matrices to produce the band forms displayed in the \hyperref[sec1]Introduction.

We begin with a geometric, coordinate-free formulation.

Proof of Theorem 2.1 We give an explicit inductive construction. Along the way, we will see that the uniqueness condition holds precisely when the choice is forced at each step.

It is convenient to restate the properties of the orthonormal basis in the theorem in the following equivalent way: for $r+1\leq i\leq n$, we have $\langle v_{i},Tv_{i-r}\rangle\geq 0$ and $Tv_{i-r}\in\operatorname{span}\{v_{1},\ldots,v_{i}\}$. Suppose inductively that$v_{1},\ldots,v_{k-1}$ have been obtained for some $r+1\leq k\leq n$, satisfying the preceding conditions for $r+1\leq i\leq k-1$. Let $w=Tv_{k-r}$; we must choose $v_{k}$ so that $\langle v_{k},w\rangle\geq 0$ and $w\in\operatorname{span}\{v_{1},\ldots,v_{k}\}$. There are two cases to consider. If $w\notin\operatorname{span}\{v_{1},\ldots,v_{k-1}\}$ then $v_{k}$ must be a multiple of $w^{\prime}=w-\sum_{i=1}^{k-1}\langle v_{i},w\rangle v_{i}$; the positivity condition further forces $v_{k}=w^{\prime}/|w^{\prime}|$, which gives $\langle v_{k},w\rangle=|w^{\prime}|>0$. If $w\in\operatorname{span}\{v_{1},\ldots,v_{k-1}\}$, then any $v_{k}\in\{v_{1},\ldots,v_{k-1}\}^{\perp}$ will do, and in this case $\langle v_{k},w\rangle=0$.

When uniqueness holds, as is generically the case, the basis may also be obtained by applying the Gram–Schmidt process to the first $n$ vectors of the sequence

We now state and prove a concrete matrix formulation in which the first $r$ coordinate vectors play the role of $v_{1},\ldots,v_{r}$. The point of the second proof is that it emphasizes the resulting band matrix rather than the change of basis; the algorithm will be used in the next subsection.

Furthermore, if strict positivity holds in (11) then $U$ and $B$ as such are unique.

Proof of Theorem 2.3 We prove existence by giving an explicit algorithm; it generalizes the Lanczos algorithm, which applies in the case $r=1$.

Stop when $k=n-r$. Let $U=U_{n-r}\cdots U_{1}$ and $B=B_{n-r}=UAU^{\dagger}$.

2 Perturbed Gaussian and spiked Wishart models

The change of basis described above interacts very nicely with the Gaussian structure in Gaussian and Wishart random matrices. The $r=1$ case of this observation is due to Trotter (1984), who described the tridiagonal forms explicitly. His forms fall into the framework of Theorem 2.1 by taking the initial vector to be fixed in the Gaussian case, and taking it to be the top row of the data matrix in the Wishart case. As we observed in Part I, the change of basis commutes with rank one additive perturbations for the Gaussian case and with rank one spiking for the Wishart case. We now extend the story to the $r>1$ setting.

In the Gaussian case, we will be perturbing in a fixed (nonrandom) subspace; without loss of generality this may be taken as the initial $r$-dimensional coordinate subspace, and so we take the basis of Theorem 2.1 that begins with the first $r$ standard basis vectors. We can therefore obtain the band form by a direct application of the algorithm from the proof of Theorem 2.3. The Wishart case is a little more complicated; here we want to perturb in the random subspace spanned by the first $r$ rows of the data matrix. Our new basis will begin with the Gram–Schmidt orthogonalization of these initial rows. As in the $r=1$ case, it is most transparent to construct a lower band form of the data matrix first, afterward realizing the band Jacobi form as its multiplicative symmetrization. In both the Gaussian and the Wishart cases, we will see that the uniqueness condition of Theorem 2.1 holds almost surely.

Let $A$ be an $n\times n$ GOE matrix. Applying the algorithm from the proof of Theorem 2.3 while keeping track of the distribution of the matrix $B_{k}$ at each step—the key of course being the unitary invariance of standard Gaussian vectors—yields the following band Jacobi random matrix $G=UAU^{\dagger}$:

As expected the perturbation shows up undisturbed in the upper-left $r\times r$ corner of $G$.

Continue in this way until the rows and columns both run out (stop alternating if one runs out before the other).

The resulting $L=V_{n\wedge(p-r)}\cdots V_{1}XU_{1}\cdots U_{n\wedge p}$ has $n\wedge p$ nonzero columns and $(n+r)\wedge p$ nonzero rows, which can be described as follows:

where we have ignored the issue of truncation in the final $r$ rows and columns ($g^{\prime}s$ and $\chi^{\prime}s$ with indices beyond the allowed range should simply be zero). The change of basis is thus $U_{1}\cdots U_{n\wedge p}$; a little thought shows that, as claimed earlier, the new basis begins with the orthogonalization of the first $r$ rows of $X$. Since the form (16) satisfies the uniqueness condition of Theorem 2.1 a.s., the basis is indeed the one given by the theorem.

where $L=VXU$ and $L_{0}=VX_{0}U$. The point is that same change of basis works in the rank $r$ spiked case, and by the lower band structure of $L_{0}$, the perturbation shows up in the upper-left $r\times r$ corner:

Viewed in terms of the algorithm used to produce $L$, the point is that the first $r$ rows of $X$ are never “mixed” together or with the lower rows, but only “rotated” within themselves.

Limits of block tridiagonal matrices

The banded forms of Section 2 may also be considered as block tridiagonal matrices with $r\times r$ blocks. In this section, we give general conditions under which such random matrices, appropriately scaled, converge at the soft spectral edge to a random Schrödinger operator on the half-line with $r\times r$ matrix-valued potential and general self-adjoint boundary condition at the origin. In Section 4, we verify these assumptions for the two specific matrix models we consider.

Proposition 3.7 establishes that the limiting operator is a.s. bounded below with purely discrete spectrum via a variational principle. The main result is Theorem 3.9, which asserts that the low-lying states of the discrete models converge to those of the operator limit.

The scalar $r=1$ case of Part I, based in turn on RRV, serves as a prototype. Care is required throughout to adapt the arguments to the matrix-valued setting, and we give a self-contained treatment.

respectively; the sub-diagonal process is of course the conjugate transpose of the super-diagonal process. (We could have absorbed $W_{n}$ into $Y_{n,1}$ as an additive constant, but keep it separate for reasons that will soon be clear. Note also that the upper-left block has $m_{n}^{2}$ rather than $2m_{n}^{2}$.) We refer to $H_{n}$ as a rank $r$ block tri-diagonal ensemble.

As in RRV and Part I, convergence rests on a few key assumptions on the potential and boundary terms just introduced. By choice, no additional scaling will be required. The role of the convergence in the first and third assumption below will be clear as soon as we define the continuum limit. The growth and oscillation bounds of the second assumption (and the lower bound implied by the third) ensure tightness of the low-lying states; in particular, they guarantee that the spectrum remains discrete and bounded below in the limit.

with respect to the compact-uniform topology (defined using any matrix norm).

(so $\triangle_{n}Y_{n,i}=\eta_{n,i}+\triangle_{n}\omega_{n,i}$) with $\eta_{n,i;j}\geq 0$ (as matrices), such that for some deterministic scalar continuous nondecreasing unbounded functions $\overline{\eta}(x)>0$, $\zeta(x)\geq 1$ not depending on $n$, and random constants $\kappa_{n}\geq 1$ defined on the same probability spaces, the following hold: the $\kappa_{n}$ are tight in distribution, and for each $n$ we have almost surely

2 Reduction to deterministic setting

In the next subsection, we will define a limiting object in terms of $Y(x)$ and $W$; we want to prove that the discrete models converge to this continuum limit in law. We reduce the problem to a deterministic convergence statement as follows. First, select any subsequence. It will be convenient to extract a further subsequence so that certain additional tight sequences converge jointly in law; Skorokhod’s representation theorem [see Ethier and Kurtz (1986)] says this convergence can be realized almost surely on a single probability space. We may then proceed pathwise.

In detail, consider (1)–(25). Note in particular that nonnegativity of the $\eta_{n,i}$ and the upper bound of (23) give that for $i=1,2$ the piecewise linear process $\{\int_{0}^{x}\eta_{n,i}\}_{x\geq 0}$ is tight in distribution, pointwise with respect to the spectral norm and in fact compact-uniformly. Given a subsequence, we pass to a further subsequence so that the following distributional limits exist jointly:

for $i=1,2$, where convergence in the first two lines is in the compact-uniform topology. We realize (3.2) pathwise a.s. on some probability space and continue in this deterministic setting.

We will assume this subsequential pathwise coupling for the remainder of the section.

3 Limiting object and variational characterization

Formally, the limiting object is the eigenvalue problem

We thus have a completely general homogeneous linear self-adjoint boundary condition. We refer to $\operatorname{span}\{u_{i}:i>r_{0}\}$ as the Dirichlet subspace and the corresponding $f_{i}$ as Dirichlet components; they will require special treatment in what follows.

where the Dirichlet part of the last term is interpreted as zero. Formally, the form $\mathcal{H}(\cdot,\cdot)$ is just the usual one $\langle\cdot,\mathcal{H}\cdot\rangle$ associated with the operator $\mathcal{H}$; the potential term has been integrated by parts and the boundary condition “built in”. See also Remark 3.5 below.

The regularity and decay conditions naturally associated with this form are given by the following weighted Sobolev norm:

where the positive part of $W$ is defined as $W^{+}=\sum_{i=1}^{r}w_{i}^{+}u_{i}u_{i}^{\dagger}$ with $w^{+}=w\vee 0$. [Define the negative part similarly with $w_{i}^{-}=-(w\wedge 0)$, so that $W=W^{+}-W^{-}$.] We refer to $\|\cdot\|_{*}$ as the $L^{*}$ norm and define an associated Hilbert space $L^{*}$ as the closure of $C_{0}^{\infty}$ under this norm. (The formal Dirichlet terms are again interpreted to be zero, but they can also be thought of as imposing the Dirichlet condition.) We record some basic facts about $L^{*}$.

Any $f\in L^{*}$ is uniformly $\mbox{H\"{o}lder}(1/2)$-continuous and satisfies $|f(x)|^{2}\leq 2\|f^{\prime}\|\|f\|$ $\leq\|f\|_{*}^{2}$ for all $x$; furthermore, $f_{i}(0)=0$ for $i>r_{0}$.

We have $|f(y)-f(x)|=|\int_{x}^{y}f^{\prime}|\leq\|f^{\prime}\||y-x|^{1/2}$. For $f\in C_{0}^{\infty}$ we have $|f(x)|^{2}=-\int_{x}^{\infty}2\operatorname{Re}f^{\dagger}f^{\prime}\leq 2\|f\|\|f^{\prime}\|\leq\|f\|_{*}^{2}$; an $L^{*}$-bounded sequence in $C_{0}^{\infty}$, therefore, has a compact-uniformly convergent subsequence, so we can extend this bound to $f\in L^{*}$ and also conclude the behaviour in the Dirichlet components.

Every $L^{*}$-bounded sequence has a subsequence converging in the following modes: (i) weakly in $L^{*}$, (ii) derivatives weakly in $L^{2}$, (iii) uniformly on compacts and (iv) in $L^{2}$.

(i) and (ii) are just Banach–Alaoglu; (iii) is the previous fact and Arzelà–Ascoli again; (iii) implies $L^{2}$ convergence locally, while the uniform bound on $\int\overline{\eta}|f_{n}|^{2}$ produces the uniform integrability required for (iv). Note that the weak limit in (ii) really is the derivative of the limit function, as one can see by integrating against functions $\mathbf{1}_{[0,x]}$ and using pointwise convergence.

By the bound in Fact 3.1 with $x=0$, the boundary term in (32) could be done away with. It is natural to include the term, however, when considering all $W$ simultaneously and viewing the Dirichlet case as a limiting case. More importantly, it clarifies the role of the boundary terms in the following key bound.

For every $0<c<1/\kappa$ there is a $C>0$ such that, for each $b>0$, the following holds for all $W\geq-b$ and all $f\in C_{0}^{\infty}$:

In particular, $\mathcal{H}(\cdot,\cdot)$ extends uniquely to a continuous symmetric bilinear form on $L^{*}\times L^{*}$.

For the first three terms of (31), we use the decomposition $Y=\int\eta+\omega$ from the previous subsection. Integrating the $\int\eta$ term by parts, (27) easily yields

Break up the $\omega$ term as follows. The moving average $\overline{\omega}_{x}=\int_{x}^{x+1}\omega$ is differentiable with $\overline{\omega}_{x}^{\prime}=\omega_{x+1}-\omega_{x}$; writing $\omega=\overline{\omega}+(\omega-\overline{\omega})$, we have

By (28), $\max(|\omega_{\xi}-\omega_{x}|,|\omega_{\xi}-\omega_{x}|^{2})\leq C_{\varepsilon}+\varepsilon\overline{\eta}(x)$ for $|\xi-x|\leq 1$, where $\varepsilon$ can be made small. In particular, the first term above is bounded absolutely by $\varepsilon\|f\|_{*}^{2}+C_{\varepsilon}\|f\|^{2}$. Averaging, we also get $|\overline{\omega}_{x}-\omega_{x}|\leq(C_{\varepsilon}+\varepsilon\overline{\eta}(x))^{1/2}$; Cauchy–Schwarz then bounds the second term absolutely by $\sqrt{\varepsilon}\int_{0}^{\infty}|f^{\prime}|^{2}+\frac{1}{\sqrt{\varepsilon}}\int_{0}^{\infty}{(C_{\varepsilon}+\varepsilon\overline{\eta})|f|^{2}}$ and thus by $\sqrt{\varepsilon}\|f\|_{*}^{2}+C_{\varepsilon}^{\prime}\|f\|^{2}$. Now combine all the terms and set $\varepsilon$ small to obtain a version of (33) with the boundary terms omitted (from both the form and the norm).

We break the boundary term in (31) into its positive and negative parts. For the negative part, Fact 3.2 gives $|f(0)|^{2}\leq(\varepsilon/b)\|f^{\prime}\|^{2}+(b/\varepsilon)\|f\|^{2}$; $W^{-}\leq b$ then implies that

which may be subtracted from the inequality already obtained. For the positive part $f(0)^{\dagger}W^{+}f(0)$, use the fact that $c\leq 1\leq C$ to simply add it in. We thus arrive at (33).

For the $L^{*}$ bilinear form bound, begin with the quadratic form bound $|\mathcal{H}(f,f)|\leq C_{c,b}\|f\|^{2}_{*}$; it is a standard Hilbert space fact that it may be polarized to a bilinear form bound [see, e.g., Section 18 of Halmos (1951)].

We say $f\in L^{*}$ is an eigenfunction with eigenvalue $\Lambda$ if $f\neq 0$ and for all $\varphi\in C_{0}^{\infty}$ we have

Note that (34) then automatically holds for all $\varphi\in L^{*}$, by $L^{*}$-continuity of both sides.

This definition represents a weak or distributional version of the problem (3.3). As further justification, integrate by parts to write the definition

(For a Dirichlet component $f_{i}$ the restriction on test functions implies that $\langle\varphi^{\prime}_{i},1\rangle=0$, so the first boundary term on the right-hand side is replaced with an arbitrary constant.) Now (35) shows that $f^{\prime}$ has a continuous version, and the equation may be taken to hold everywhere. In particular, $f$ satisfies the boundary condition of (3.3) classically. [For a Dirichlet component, we just find that the arbitrary constant is $f_{i}^{\prime}(0)$.] One can also view (35) as a straightforward integrated version of the eigenvalue equation in which the potential term has been interpreted via integration by parts. This equation will be useful in Lemma 3.6 below and is the starting point for the development in Section 5.

We now characterize the eigenvalues and eigenfunctions variationally. As usual, it follows from the symmetry of the form that eigenvalues are real (and eigenfunctions with distinct eigenvalues are $L^{2}$-orthogonal). The $L^{2}$ part of the lower bound in (33) says the spectrum is bounded below. The rest of (33) implies that there are only finitely many eigenvalues below any given level: a sequence of normalized eigenfunctions with bounded eigenvalues must have an $L^{2}$-convergent subsequence by Fact 3.2. At a given level, more is true.

By linearity, it suffices to show a solution of (35) with $f^{\prime}(0)=f(0)=0$ must vanish identically. Integrate by parts to write

which implies that $|f^{\prime}(x)|\leq C(x)\int_{0}^{x}|f^{\prime}|$ with some $C(x)<\infty$ increasing in $x$. Gronwall’s lemma then gives $|f^{\prime}(x)|=0$ for all $x\geq 0$.

There is a well-defined $(k+1)$st lowest eigenvalue $\Lambda_{k}$, counting with multiplicity. The eigenvalues $\Lambda_{0}\leq\Lambda_{1}\leq\ldots$ together with an orthonormal sequence of corresponding eigenvectors $f_{0},f_{1},\ldots$ are given recursively by the variational problem

in which the minimum is attained and we set $f_{k}$ to be any minimizer.

Since we must have $\Lambda_{k}\to\infty$, $\{\Lambda_{0},\Lambda_{1},\ldots\}$ exhausts the spectrum and the resolvent operator is compact. We do not make this statement precise.

Proceed inductively, minimizing now over the orthocomplement $\{f\in L^{*}:\|f\|=1,f\perp f_{0},\ldots,f_{k-1}\}$. Again, $L^{2}$-convergence of a minimizing sequence guarantees that the limit remains admissible; as before, the limit is in fact a minimizer; conclude by applying the arguments of the previous paragraph with $\varphi,g$ also restricted to the orthocomplement.

4 Statement

Let $H_{n}$ be a rank $r$ block tr-diagonal ensemble as in (19) satisfying Assumptions 1–3, and let $\lambda_{n,k}$ be its $(k+1)$st lowest eigenvalue. Define the associated form $\mathcal{H}$ as in (31) and let $\Lambda_{k}$ be its a.s. defined $(k+1)$st lowest eigenvalue. In the deterministic setting of subsequential pathwise coupling, $\lambda_{n,k}\to\Lambda_{k}$ for each $k=0,1,\ldots.$ Furthermore, a sequence of normalized eigenvectors corresponding to $\lambda_{n,k}$ is precompact in $L^{2}$ norm, and every subsequential limit is an eigenfunction corresponding to $\Lambda_{k}$. Finally, convergence holds uniformly over possible $W_{n},W\geq-b>-\infty$. One recovers the corresponding distributional tightness and convergence statements for the full sequence, jointly for $k=0,1,\ldots$ in the sense of finite-dimensional distributions and jointly over $W_{n},W$.

The proof will be given over the course of the next two subsections.

5 Tightness

with the nonnegative part $W_{n}^{+}$ defined as before.

When considering just a single $W_{n},W$, the boundary term in (3.5) is really only required when the limit includes Dirichlet terms; it is simpler, however, not to distinguish the two cases here. More importantly, including this term clarifies the role of the boundary term in the following key bound. Note that the original case considered in RRV has $W_{n}=m_{n}$ in our notation. (The $H_{n}$ form and $L^{*}_{n}$ norm there contained a term $m_{n}|v_{0}|^{2}$, though it is hidden in the fact that, in our notation, they use $\triangle_{n}$ in place of $D_{n}$.)

The potential term $\langle v,Vv\rangle=\int_{0}^{\infty}v^{\dagger}Vv$, defined in (18), is analyzed according to (22):

Together with $|D_{n}v|^{2}$, the $\eta$-terms provide the structure of the bound as we now show. Afterward we will control the $\omega$-terms and lastly deal with the boundary term.

Recall (23) and that $\eta_{i}\geq 0$. For an upper bound, rearrange $(v-Tv)^{\dagger}\eta_{2}(v-Tv)\geq 0$ to

Now $\int\overline{\eta}|Tv|^{2}=\int(T^{\dagger}\overline{\eta})|v|^{2}\leq\int\overline{\eta}|v|^{2}$ since $\overline{\eta}$ is nondecreasing, and we obtain

Toward a lower bound, we use the slightly tricky rearrangement $0\leq(\frac{1}{2}v+Tv)^{\dagger}\eta_{2}(\frac{1}{2}v+Tv)=3\operatorname{Re}v^{\dagger}\eta_{2}Tv+(Tv-v)^{\dagger}\eta_{2}(Tv-v)-\frac{3}{4}v^{\dagger}\eta_{2}v$. With (24), we get

We handle the $\omega$-terms with a discrete analogue of the decomposition used in the continuum proof. Consider the moving average

which has $\triangle\overline{\omega}_{i}=(m/\lfloor m\rfloor)(T^{\lfloor m\rfloor}-1)\omega_{i}$; it is convenient to extend $\omega_{i}(x)=\omega_{i}(\lceil n/r\rceil/m_{n})$ for $x>\lceil n/r\rceil/m_{n}$. Decompose $\omega_{i}=\overline{\omega}_{i}+(\omega_{i}-\overline{\omega}_{i})$. For the $\omega_{1}$-term,

By (25) and Cauchy–Schwarz, the first term is bounded absolutely by $(C_{\varepsilon}+\varepsilon\overline{\eta})|v|^{2}$ and its integral by $\varepsilon\|v\|_{*}^{2}+C_{\varepsilon}\|v\|^{2}$. The second term calls for a summation by parts:

The averaged bound $|\overline{\omega}_{1}-\omega_{1}|\leq(C_{\varepsilon}+\varepsilon\overline{\eta})^{1/2}$ and Cauchy–Schwarz bound the integrand

and its integral by $\sqrt{\varepsilon}\|v\|_{*}^{2}+C_{\varepsilon}^{\prime}\|v\|^{2}$. One thus obtains a similar bound on $|\langle v,(\triangle\omega_{1})v\rangle|$.

There are corresponding bounds for the $\omega_{2}$-terms. For the $\overline{\omega}_{2}$-term, use $2|v||Tv|\leq|v|^{2}+|Tv|^{2}$. For the $(\omega_{2}-\overline{\omega}_{2})$-term, modify the summation by parts:

Incorporating all the $\omega$-terms into (39), (40) and setting $\varepsilon$ small, we obtain (37) but with the boundary terms omitted (from both the form and the norm).

We break the boundary term in (38) into its positive and negative parts. A discrete analogue of a bound from Fact 3.1 will be useful:

It gives $|v(0)|^{2}\leq(\varepsilon/b)\|Dv\|^{2}+(b/\varepsilon)\|v\|^{2}$, and then $W^{-}\leq b$ implies that

which may be subtracted from the inequality already obtained. The positive part may simply be added in using that $c\leq 1\leq C$. We thus arrive at (37).

If the $W_{n}$ are not bounded below then the lower bound in (37) breaks down: in fact, the bottom eigenvalue of $H_{n}$ really goes to $-\infty$ like minus the square of the bottom eigenvalue of $W_{n}$. This is the supercritical regime.

6 Convergence

We begin with a simple lemma, a discrete-to-continuous version of Fact 3.2.

Let $f_{n}\to f$ be as in the hypothesis and conclusion of Lemma 3.15. Then for all $\varphi\in C_{0}^{\infty}$ we have $\langle\varphi,H_{n}f_{n}\rangle\to\mathcal{H}(\varphi,f)$. In particular, $\mathcal{P}_{n}\varphi\to\varphi$ in this way and so

Since $\varphi$ is compactly supported, we have $R_{n}\varphi=\varphi$ for $n$ large and the $R_{n}$s may be dropped. By assumption $D_{n}f_{n}$ is $L^{2}$ bounded and $D_{n}f_{n}\to f^{\prime}$ weakly in $L^{2}$, so by the preceding observations $D_{n}\varphi\to_{L^{2}}\varphi^{\prime}$ and

For the potential term, we must verify that

converges to $-\langle\varphi^{\prime},Yf\rangle-\langle\varphi,Yf^{\prime}\rangle$. Recall by Assumption 1 (1) and (3.2) that $Y_{n,i}\to Y_{i}$ compact-uniformly ($i=1,2$) and $Y=Y_{1}+\frac{1}{2}(Y_{2}+Y_{2}^{\dagger})$. Writing $Y_{n}=Y_{n,1}+\frac{1}{2}(Y_{n,2}+Y^{\dagger}_{n,2})\to Y$ (and disregarding the notational collision with $Y_{i}$), we first approximate $V_{n}$ by $\triangle Y_{n}$:

which converges to the desired limit by the observations preceding the lemma together with the assumptions on $f_{n}$ and the fact that $T_{n}\varphi\to_{L^{2}}\varphi$ in $L^{2}$ since $m_{n}\|T_{n}\varphi-\varphi\|=\|D_{n}\varphi\|$ is bounded. The error in the above approximation comes as a sum of $T_{n}$ and $T_{n}^{\dagger}$ terms. Consider twice the $T_{n}$ term:

Finally, for the boundary terms Assumption 3 gives

where in the Dirichlet case $i>r_{0}$ the left-hand side vanishes for $n$ large because $\varphi_{i}$ is supported away from 0.

Turning to the second statement, we must verify that $\mathcal{P}_{n}\varphi\to\varphi$ as in Lemma 3.15. The uniform $L^{*}_{n}$ bound on $\mathcal{P}_{n}\varphi$ follows from the following observations: $\|(\mathcal{P}_{n}\varphi){\sqrt{1+\overline{\eta}}}\|=\|\mathcal{P}_{n}\varphi{\sqrt{1+\overline{\eta}}}\|\leq\|\varphi{\sqrt{1+\overline{\eta}}}\|$; for $n$ large enough that $R_{n}\varphi=\varphi$ we have $\|D_{n}\mathcal{P}_{n}\varphi\|=\|\mathcal{P}_{n}D_{n}\varphi\|\leq\|D_{n}\varphi\|\leq\|\varphi^{\prime}\|$ (Young’s inequality); for the boundary term note that $(\mathcal{P}_{n}\varphi)_{i}(0)$ is bounded if $i\leq r_{0}$ and in fact vanishes for $n$ large if $i>r_{0}$. The convergence is easy: $\mathcal{P}_{n}\varphi\to\varphi$ compact-uniformly and in $L^{2}$, and for $g\in L^{2}$ we have $\langle g,D_{n}\mathcal{P}_{n}\varphi\rangle=\langle\mathcal{P}_{n}g,D_{n}\varphi\rangle\to\langle g,\varphi^{\prime}\rangle$.

We finish by recalling the argument to put all the pieces together. A technical point: unlike in previous treatments we do not assume that the eigenvalues are simple.

Proof of Theorem 3.9 We first show that for all $k$ we have $\underline{\lambda}_{k}=\liminf\lambda_{n,k}\geq\Lambda_{k}$. Assume that $\underline{\lambda}_{k}<\infty$. The eigenvalues of $H_{n}$ are uniformly bounded below by Lemma 3.13, so there is a subsequence along which $(\lambda_{n,1},\ldots,\lambda_{n,k})\to(\xi_{1},\ldots,\xi_{k}=\underline{\lambda}_{k})$. By the same lemma, corresponding orthonormal eigenvector sequences have $L^{*}_{n}$-norm uniformly bounded. Pass to a further subsequence so that they all converge as in Lemma 3.15. The limit functions are orthonormal; by Lemma 3.16 they are eigenfunctions with eigenvalues $\xi_{j}\leq\underline{\lambda}_{k}$ and we are done.

We proceed by induction, assuming the conclusion of the theorem up to $k-1$. For $j=0,\ldots,k-1$ let $v_{n,j}$ be orthonormal eigenvectors corresponding to $\lambda_{n,j}$; for any subsequence we can pass to a further subsequence such that $v_{n,j}\to_{L^{2}}f_{j}$, eigenfunctions corresponding to $\Lambda_{j}$. Take an orthogonal eigenfunction $f_{k}$ corresponding to $\Lambda_{k}$ and find $f_{k}^{\varepsilon}\in C_{0}^{\infty}$ with $\|f_{k}^{\varepsilon}-f_{k}\|_{*}<\varepsilon$. Consider the vector

The $L^{*}_{n}$-norm of the sum term is uniformly bounded by $C\varepsilon$: indeed, the $\|v_{n,j}\|_{*n}$ are uniformly bounded by Lemma 3.13, while the coefficients satisfy $|\langle v_{n,j},f_{k}^{\varepsilon}\rangle|\leq\|f_{k}^{\varepsilon}-f_{k}\|+\|v_{n,j}-f_{j}\|<2\varepsilon$ for large $n$. By the variational characterization in finite dimensions and the uniform $L^{*}_{n}$ form bound on $\langle\cdot,H_{n}\cdot\rangle$ (by Lemma 3.13) together with the uniform bound on $\|\mathcal{P}_{n}f_{k}^{\varepsilon}\|_{*n}$ (by Lemma 3.16), we then have

where $o_{\varepsilon}(1)\to 0$ as $\varepsilon\to 0$. But (41) of Lemma 3.16 provides $\lim\langle\mathcal{P}_{n}f_{k}^{\varepsilon},\break H_{n}\mathcal{P}_{n}f_{k}^{\varepsilon}\rangle=\mathcal{H}(f_{k}^{\varepsilon},f_{k}^{\varepsilon})$, so the right-hand side of (3.6) is

Now letting $\varepsilon\to 0$, we conclude $\limsup\lambda_{n,k}\leq\Lambda_{k}$.

Thus, $\lambda_{n,k}\to\Lambda_{k}$; Lemmas 3.13 and 3.15 imply that any subsequence of the $v_{n,k}$ has a further subsequence converging in $L^{2}$ to some $f\in L^{*}$; Lemma 3.16 then implies that $f$ is an eigenfunction corresponding to $\Lambda_{k}$. Finally, convergence is uniform over $W_{n},W\geq-b$ since the bound 3.13 is.

CLT and tightness for Gaussian and Wishart models

We now verify Assumptions 1–3 of Section 3 for the band Jacobi forms of Section 2, and thus prove Theorems 1.2 and 1.3 via Theorem 3.9.

The technical tool we use to establish (1) is a functional central limit theorem for convergence of discrete time processes with independent increments of given mean and variance (and controlled fourth moments) to Brownian motion plus a nice drift. Appearing as Corollary 6.1 in RRV, it is just a tailored version of a much more general result given as Theorem 7.4.1 in Ethier and Kurtz (1986). We record it here.

uniformly for $j/m_{n}$ on compact sets as $n\to\infty$. Then $y_{n}(x)=y_{n,\lfloor m_{n}x\rfloor}$ converges in law, with respect to the compact-uniform topology, to the process $h(x)+ab_{x}$ where $b_{x}$ is a standard Brownian motion.

Since the limit is a.s. continuous, Skorokhod convergence (the topology used in the references) implies uniform convergence on compact intervals [see Theorem 3.10.2 in Ethier and Kurtz (1986)] and we may as well speak in terms of the latter.

As usual, this soft-edge scaling can be predicted as follows. Centering $G_{n}$ by $2\sqrt{n}$ gives, to first order, $\sqrt{n}$ times the discrete Laplacian on blocks of size $r$. With space scaled down by $m_{n}$, the Laplacian must be scaled up by $m_{n}^{2}$ to converge to the second derivative. Finally, the scaling $m_{n}=n^{1/3}$ is determined by convergence of the next order terms to the noise and drift parts of the limiting potential.

Decompose $H_{n}$ as in (19), (3.1). The upper-left block is

we want the boundary term $W_{n}$ to absorb the “extra” $m_{n}^{2}$ (the 2 in the right-hand side “should be” a 1) and the perturbation in order to make $Y_{n,1;0}$ small just like the subsequent increments of $Y_{n,i}$. We therefore set

With this choice Assumption 3 is an immediate consequence of the hypotheses of Theorem 1.2. The processes $Y_{n,1},Y_{n,2}$ are determined and it remains to verify Assumptions 1 and 2.

Proof of (1), Gaussian case Define scalar processes $y_{k,l}$ for $1\leq l\leq r$ and $l\leq k\leq l+r$ by

Note that the $y_{k,l}$ are independent increment processes that are mutually independent of one another. With the usual embedding $j=\lfloor n^{1/3}x\rfloor$, Proposition 4.1 together with standard moment computations for Gaussian and Gamma random variables—in particular

for $\alpha$ large [valid since we consider $j=O(n^{1/3})$ here]—leads to the convergence of processes

Noting that the two Brownian motions in each entry are independent and that the entries on and below the diagonal are independent of each other, we conclude that this limiting matrix process is distributed as $Y(x)$ in $(\ref{YSA})$.

We turn to Assumption 2. Here, we need bounds over the full range $0\leq j\leq\lceil n/r\rceil-1$. Recall that we can extend the $Y_{n,i}$ processes beyond the end of the matrix arbitrarily ($R_{n}$ takes care of the truncation), and it is convenient to “continue the pattern” for an extra block or two by setting $\chi_{\alpha}=0$ for $\alpha<0$. For the decomposition (22), we simply take $\eta_{n,i}$ to be the expectation of $\triangle Y_{n,i}$ and $\triangle\omega_{n,i}$ to be its centered version; the components of $\eta_{n,i}$ are then easily estimated and those of $\omega_{n,i}$ become independent increment martingales. We further set $\overline{\eta}(x)=rx$.

Proof of (23)–(25), Gaussian case From (46), we have $\eta_{n,1;j}=0$ and

for some fixed $c$, which yields the matrix inequalities

and verifies (23) with $\overline{\eta}(x)=rx$. Separately, we have the upper bound (24):

are tight over $n$. Squaring, bounding the outer supremum by the corresponding sum, and then taking expectations gives

where we have used the $L^{p}$ maximum inequality for martingales [see, e.g., Proposition 2.2.16 of Ethier and Kurtz (1986)]. To bound the latter expectation, expand the fourth power to obtain $O(m_{n}^{2})$ nonzero terms that are $O(m_{n}^{-2})$ with constants independent of $x$ and $n$. It follows that the entire sum is uniformly bounded over $n$, as required.

2 The Wishart case

See Part I for detailed heuristics behind the scaling; written in this way, it allows that $p,n\to\infty$ together arbitrarily, that is, only $n\wedge p\to\infty$. It is useful to note that

Decompose $H_{n,p}$ as in (19), (3.1). The upper-left block is

Once again, Assumption 3 follows immediately from the hypotheses of Theorem 1.3.

We must still deal with the perturbed term in $Y_{1;0}$ and show that

in probability. We defer this to the end of the proof of Assumption 1, to which we now turn. As in the Gaussian case, $Y$ is given by (43).

Proof of (1), Wishart case By the preceding paragraph it suffices to treat the null case $\Sigma=I$ and afterward check (50). Define processes $y_{k,l}$ for $1\leq l\leq r$ and $l\leq k\leq l+r$ by (44) as in the Gaussian case. From (16) with the centering and scaling of (48) and (3.1), we obtain

where the $O(1)$ terms stand in for the interior Gaussian sums of (16), all of whose moments are bounded uniformly in $n,p$. Since $m^{1+k}/(np)^{k/2}\leq m^{1-2k}=o(1)$ for $k\geq 1$, these terms are negligible in the scaling of Proposition 4.1 in the sense that the associated processes converge to the zero process. Next, use that expressions of type $\chi_{n}-\sqrt{n}$ are $O(1)$ in the same sense, and that $\sqrt{n}-\sqrt{n-j}=O(j/\sqrt{n})=O(m/\sqrt{n})=o(1)$ since we consider $j/m$ bounded here (and similarly for $p$), to write

jointly for $1\leq k,l\leq r$. After the dust clears, we thus arrive at exactly the same limiting process as in the Gaussian case, namely (43).

Turning to Assumption 2, we may continue the processes $Y_{n,i}$ past the end of the matrix for convenience just as in the Gaussian case. The Wishart case presents an additional issue at the “end” of the matrix: recall that the final $r$ rows and columns of $S$ in (16) may have some apparently nonzero terms set to zero. However, these changes are easily absorbed into the bounds that follow. For (22), we once again take $\eta_{n,i}$ to be the expectation of $\triangle Y_{n,i}$ and $\triangle\omega_{n,i}$ to be its centered version. We also set $\overline{\eta}(x)=rx$ as before.

Proof of (23)–(25), Wishart case This time we have

Using (47) one finds, for some constant $c$, that

Alternative characterizations of the laws

In this section, we derive the SDE and PDE characterizations, proving Theorems 1.5 and 1.6.

For each noise path $B_{x}$, the eigenvalue equation $\mathcal{H}_{\beta,W}f=\lambda f$ can be rewritten as a first-order linear ODE with continuous coefficients. We begin with the formal second-order linear differential equation

Now let $g=f^{\prime}-{\sqrt{2}}Bf$. The equation becomes

In other words, the pair $(f(x),g(x))$ formally satisfies the first-order linear system

Since $B_{0}=0$, $g$ simply replaces $f^{\prime}$ in the initial condition (53). If one prefers, this condition can be written in the standard form

2 Matrix oscillation theory

The matrix generalization of Sturm oscillation theory goes back to the classic work of Morse Morse (1932) [see also Morse (1973)]. Textbook treatments of self-adjoint differential systems include that of Reid (1971). Our reference will be the paper of Baur and Kratz (1989), which allows sufficiently general boundary conditions.

We first consider the eigenvalue problem on a finite interval $[0,L]$ with Dirichlet boundary condition $f(L)=0$ at the right endpoint. In the scalar-valued setting, the number of eigenvalues below $\lambda$ is found to coincide with the number of zeros of $f$ (the solution of the initial value problem) that lie in $(0,L)$. The correct generalization to the matrix-valued setting involves tracking a matrix whose columns form a basis of solutions, and counting the so-called “focal points”.

The idea is that focal points are isolated and move continuously to the left as $\lambda$ increases. For sufficiently negative $\lambda$, there are no focal points on $(0,L]$; each time $\lambda$ passes an eigenvalue, a new focal point is introduced at $L$.

We indicate how the proposition follows from the results of Baur and Kratz (1989). Note that Conditions (A1), (A2) on page 337 are satisfied by our coefficients, and that (A3) on page 340 is satisfied by our boundary conditions. Theorem 1 on page 345 thus applies. See (3.5) on page 341 for the definition of $\Lambda(\lambda)$; the Dirichlet condition at $L$ gives the particularly simple result that the right-hand side of (4.1) vanishes, so the quantity $n_{2}(\lambda)$ is constant. Theorem 2 applies as well, and we obtain $n_{1}(\lambda)-n_{1}=n_{3}(\lambda)$. Here, $n_{1}(\lambda)$ is the number of focal points in $[0,L)$, $n_{1}=\lim_{\lambda\to-\infty}n_{1}(\lambda)$ and $n_{3}(\lambda)$ is the number of eigenvalues below $\lambda$. To finish, we consult Theorem 3 on page 353; noting that (A4′) is satisfied by Section 7.2, page 365, to find that $n_{1}$ is simply the multiplicity of the focal point at 0. The oscillation result follows. For the assertion about the spectrum, we apply Theorem 4, noting that (A5), page 358 holds by (i) there, and (A6), page 359 also holds.

We conclude the following for our matrix system.

A soft argument now recovers an oscillation theorem for the original half-line problem.

The variational problem for $\mathcal{H}_{L}$ simply minimizes over the subset of $L^{*}$ functions that vanish on $[L,\infty)$; the Dirichlet condition is important here. It follows immediately that $\Lambda_{L,k}\geq\Lambda_{k}$, using the min–max formulation of the variational characterization. Proceed by induction, assuming that $\Lambda_{L,j}\to\Lambda_{L}$ for $j=0,\ldots,k-1$.

Let $f_{L,j}$ be orthonormal eigenvectors corresponding to $\Lambda_{L,j}$. By the induction hypothesis, the variational characterization for $\mathcal{H}$ and the finite-dimensionality of its eigenspaces, every subsequence has a further subsequence such that $f_{L,j}\to_{L^{2}}f_{j}$, eigenvectors corresponding to $\Lambda_{j}$. Let $f_{k}$ be an orthogonal eigenvector corresponding to $\Lambda_{k}$ and take $f_{k}^{\varepsilon}$ compactly supported with $\|f_{k}^{\varepsilon}-f_{k}\|_{*}<\varepsilon$. Let

For large $L$, the inner products are at most $2\varepsilon$, so $\|g_{L}-f_{k}\|_{*}\leq c\varepsilon$. Noting that $g_{L}$ is eventually supported on $[0,L]$, the variational characterization gives

and the right-hand side tends to ${\mathcal{H}(f_{k},f_{k})}/{\langle f_{k},f_{k}\rangle}=\Lambda_{k}$ as $\varepsilon\to 0$.

3 Riccati SDE: Stochastic airy meets dyson

Let $(F,G)$ be a conjoined basis for (54) as defined in the previous subsection. Then, on any interval with no focal points, the matrix $Q=GF^{-1}$ is self-adjoint and satisfies the matrix Riccati equation

Now let $P=F^{\prime}F^{-1}$. While $P=Q+{\sqrt{2}}B$ is not differentiable, by (56) it certainly satisfies the integral equation

if $[x_{1},x_{2}]$ is free of focal points. In other words, $P$ is a strong solution of the Itô equation

Consider the eigenvalues $p_{1},\ldots,p_{r}$ of $P$. The main point is that the drift term in (58) is unitarily equivariant and passes through the usual derivation of Dyson’s Brownian motion [Dyson (1962)]. The eigenvalues therefore evolve as an autonomous Markov process.

To describe the law on paths we need a space, and there are two issues: it will be necessary to keep the eigenvalues ordered but also allow for explosions/restarts. We therefore define a sequence of Weyl chambers $C_{k}\subset(-\infty,\infty]^{r}$ by

Writing $X=\dot{A}(0)$ and $\nabla_{X}$ for the directional derivative, and taking$v_{1}(0),\ldots,v_{r}(0)$ to be the standard basis, we find

Returning to (58), at each fixed time $x$ we can change to the diagonal basis for $P_{x}$ because the noise term is invariant in distribution and the drift term is equivariant. Itô’s lemma amounts to formally writing $dp_{i}=\nabla_{dP}p_{i}+\frac{1}{2}\nabla_{dP}^{2}p_{i}$ and using that $dB_{ii}$ are jointly distributed as $\sqrt{2/\beta}\,db_{i}$ for $i=1,\ldots,r$ while $|dB_{ij}|^{2}=dt$ for $j\neq i$. We thus arrive at (59).

Recall that the evolution of $P$ through a focal point is still described by an SDE, after changing coordinates. The same is therefore true of $\mathbf{p}$ through an explosion; the form (60) is obtained from (59) by an application of Itô’s lemma.

Just as with the usual Dyson’s Brownian motion, the $p_{i}$ are almost surely distinct at all positive times: $\mathbf{p}(x)\in\mathcal{C}$ for all $x>0$. One can show this “no collision property” holds for any solution of (59), (60), even with an initial condition $\mathbf{p}(0)\in\partial C_{0}$. (Technically, one defines an entrance law from $\partial\mathcal{C}$ by a limiting procedure.) Since the coefficients are regular inside $\mathcal{C}$, this suffices to prove uniqueness of the law. See Anderson, Guionnet and Zeitouni (2010), Section 4.3.1 for a detailed proof in the driftless case.

Proof of Theorem 1.5 Explosions of $\mathbf{p}$ as in Theorem 5.4 correspond to focal points of $F$ for each $\lambda$. By Theorem 5.3, the total number of explosions $K$ is equal to the number of eigenvalues strictly below $\lambda$. (Notice that $\mathbf{p}$ ends up in $C_{K}$.) For a fixed $\lambda$, translation invariance of the driving Brownian motions $b_{i}$ allows one to shift time $x\mapsto x-\lambda/r$ and use (3) started at $x_{0}=-\lambda/r$. Putting $a=-\lambda$ we have $\operatorname{\mathbf{P}}(-\Lambda_{k}\leq a)=\operatorname{\mathbf{P}}(\Lambda_{k}\geq\lambda)=\operatorname{\mathbf{P}}_{a/r,\mathbf{w}}(K\leq k)$ as required.

4 PDE and boundary value problem

We now prove the PDE characterization, Theorem 1.6. We will need two properties of the eigenvalue diffusion.

Let $\mathbf{p}:[x_{0},\infty)\to\overline{\mathcal{C}}$ have law $\operatorname{\mathbf{P}}_{x_{0},\mathbf{w}}$ as in (3) and let $K$ be the number of explosions. Then the following hold:

Given $x_{0},k$, $\operatorname{\mathbf{P}}_{x_{0},\mathbf{w}}(K\leq k)$ is increasing in $\mathbf{w}$ with respect to the partial order $\mathbf{w}\leq\mathbf{w}^{\prime}$ given by $w_{i}\leq w_{i}^{\prime}$, $i=1,\ldots,r$.

$\operatorname{\mathbf{P}}_{x_{0},\mathbf{w}}$-almost surely, $p_{1},\ldots,p_{r}$ remain bounded below in $C_{K}$ (after the last explosion), or equivalently in $C_{0}$ on the event $\{K=0\}$.

Part (i) is a consequence Theorem 1.5 and Remark 1.1, the pathwise monotonicity of the eigenvalues $\Lambda_{k}$ as a function of the boundary parameter $W$ with respect to the usual matrix partial order. It can also be seen from the related fact that the matrix partial order is preserved pathwise by the matrix Riccati equation (58), which implies that a solution started from $W$ explodes no later than one started from $W^{\prime}\geq W$. This fact holds for the $P$ evolution if it holds for the $Q$ evolution (56), and for the latter it is Theorem IV.4.1 of Reid (1972).

Part (ii) follows from the stronger assertion that $p_{i}\sim\sqrt{rx}$ as $x\to\infty$. In the $r=1$ case, this is Proposition 3.7 of RRV. Heuristically, the single particle drift linearizes at the stable equilibrium $\sqrt{rx}$ to $2\sqrt{rx}(\sqrt{rx}-p_{i})$; even with the repulsion terms one expects fluctuations of variance only $C/\sqrt{x}$. We omit the proof.

For $F^{k}$, there is the following more general picture. Consider the PDE in $\overline{C}_{0}\cup\cdots\cup\overline{C}_{k}$, defined across the seams by changing coordinates as in (60). Put the boundary condition (7) on all the chambers and (8) on the bottom of $\overline{C}_{k}$. Then the solution is $F^{k}$ in $\overline{C}_{0}$; the reason is the same as for $F=F^{0}$, but now using (5) and the hitting event “at most $k$ explosions”. Similarly, the solution is $F^{k-1}$ in $\overline{C}_{1}$ and so on down to $F^{0}$ in $\overline{C}_{k}$. Continuity holds across the seams and (1.6) follows after permuting coordinates.

Connection with Painlevé II

In Part I, we used the PDE characterization to give new proofs of certain Painlevé II formulas for the single-parameter (rank one deformed) distribution functions $F_{\beta}(x;w)$ in the cases $\beta=2,4$, in particular recovering the Painlevé II representations for the corresponding undeformed Tracy–Widom distributions by taking $w\to\infty$. The Painlevé formulas appeared originally in Baik and Rains (2000; 2001) in a different context; in the random matrix theory setting, Baik (2006) derived them from the BBP result in the case $\beta=2$ but they are new for $\beta=4$ when $w\neq 0$ [see Wang (2008)].

Baik (2006) also derives a Painlevé II formula for the multi-parameter distribution function $F_{2}(x;w_{1},\ldots,w_{r})$. While we do not have a full independent proof at present, we used the computer algebra system Maple to verify symbolically that it does indeed satisfy our PDE (6) at $\beta=2$ for $r=2,3,4,5$. Since this article was first posted, a pencil-and-paper proof for all $r$ was found [Bloemendal and Baik (2013)]. We first state Baik’s formula and then briefly describe the symbolic computation.

Let $u(x)$ be the Hastings–McLeod solution of the homogeneous Painlevé II equation

where $\operatorname{Ai}(x)$ is the Airy function. Put

Equation (64) is one member of the Lax pair for the Painlevé II equation. The other member of the pair is

Our symbolic verification for small values of $r$ consisted of the following steps. The differential relations given by (61)–(65) were encoded as formal substitution rules. The determinant in (66) was expanded (this step becomes problematic for larger $r$…!) and the result plugged into our PDE (6). The substitution rules were then applied repeatedly. Finally, the result was factored using Maple’s built-in command. Each time, the output contained the factor