This post will provide an accessible introduction to work contained in my joint paper with Thomas Bothner. In this work we find integrable structure in the elliptic Ginibre ensemble and find interesting connections to the theory of finite temperature Airy processes. Our paper also performs a nonlinear steepest descent analysis of the cumulative distribution function of the rightmost particle, but this is left out of this post to keep the discussion nontechnical.

Background and motivation

Suppose you have an ecosystem with nn species, with populations x1,,xnx_1, \dots, x_n. Is this ecosystem stable? This question was raised by Robert May in his article “Will a Large Complex System be Stable?” (1972). Suppose as a first approximation that each species is isolated and the only thing limiting population is competition over some non-organic resources. You could model this by the system of uncoupled ODEs

x˙i=μxi\dot{x}_i = - \mu x_i

for μ>0\mu > 0. We have shifted the populations so that xi=0x_i = 0 is the equilibrium. Now if we wanted to model interactions between species we could write

x˙i=μxi+j=1nAijxj\dot{x}_i = - \mu x_i + \sum_{j=1}^n A_{ij} x_j

where AA is some interaction matrix. AijA_{ij} represents the effect of the population of species jj on the growth rate of species ii. (Note that there is no reason to think this matrix is symmetric.) This could be written in vector form

x˙=(μI+A)x.\dot{\mathbf{x}} = (- \mu \mathbb{I} + A)\mathbf{x} .

It is clear that this system is stable when μ>maxiλi\mu > \max_i \Re \lambda_i and unstable when μ<maxiλi\mu < \max_i \Re \lambda_i, where λ1,,λn\lambda_1, \dots, \lambda_n are the eigenvalues of AA. Given that we expect nn to be very large and the interactions between the species very complex, we cannot expect to model AA exactly. Instead we take AA to be random. Note that the probability of stability, P(maxiλi<μ)\mathbb{P}(\max_i \Re \lambda_i < \mu), is exactly the cumulative distribution function for maxiλi\max_i \Re \lambda_i.

In our work we are interested in AMn(C)A \in \mathcal{M}_n(\mathbb{C}), so it models stability of a linear system over the field of complex numbers. I will list three candidate distribution functions for AA. This list is by no means exhaustive – you can give many others – however these three cases are “exactly solvable” in that one can compute their correlation functions explicitly.

The (complex) Ginibre Ensemble

This ensemble was introduced by Ginibre in 1965. This is an ensemble of n×nn \times n matrices with independent, identically distributed matrix elements, each given by complex Gaussians. One could write the density as

PGin(X)dX=πn2etrXXdX.P_{\mathrm{Gin}}(X) \, \mathrm{d}X = \pi^{-n^2} \mathrm{e}^{-\mathrm{tr}\, X X^\ast} \, \mathrm{d}X .

The eigenvalues have the joint distribution

ϱn(z1,,zn)=Cnek=1nzk21i<jnzizj2.\varrho_n(z_1, \dots, z_n) = C_n \mathrm{e}^{-\sum_{k=1}^n |z_k|^2} \prod_{1 \leq i < j \leq n } |z_i - z_j|^2 .

The Ginibre ensemble obeys a “circular law.” If we define the 1-point density as

ϱ(x)=1n(C)n1ϱn(x,z2,,zn)d2z2d2zn\varrho(x) = \frac{1}{n} \int_{(\mathbb{C})^{n-1}} \varrho_n(x, z_2, \dots, z_n) \, \mathrm{d}^2 z_2 \dots \mathrm{d}^2 z_n

then ϱ(nx)1πχx<1\varrho(\sqrt{n}x ) \to \frac{1}{\pi} \chi_{\lvert x \rvert < 1} as nn \to\infty in the weak-\ast sense. That is, the spectral density tends towards the unit disk. Furthermore, the eigenvalue with largest real part is asymptotically Gumbel distributed. More precisely

P(maxiλin+γn4+t4γn)eet\mathbb{P}\left( \max_i \Re \lambda_i \leq \sqrt{n} + \sqrt{\frac{\gamma_n}{4}} + \frac{t}{\sqrt{4 \gamma_n}} \right) \to \mathrm{e}^{-\mathrm{e}^{-t}}

as nn \to \infty, where

γn=12(lnn5lnlnnln(2π4)).\gamma_n = \frac{1}{2}(\ln n - 5 \ln \ln n - \ln (2\pi^4)) .

This result was originally due to Bender (2010) and has recently been given a simple proof by G. Cipolloni, L. Erdős, D. Schröder, and Y. Xu (2022). Furthermore, we show in work forthcoming on the arXiv, that the real parts of the eigenvalues converge locally, in the bulk of the spectrum, to a Poisson point process as nn\to \infty.

The Gaussian Unitary Ensemble

The Gaussian Unitary Ensemble (GUE) is an ensemble of n×nn \times n complex Hermitian random matrices. The “unitary” in the name comes from the fact that the ensemble has a unitary symmetry. The density is formally very similar to the Ginibre case

PGUE(X)dX=Cnetr(X2)dX.P_{\mathrm{GUE}}(X) \, \mathrm{d}X = C^\prime_n \mathrm{e}^{-\mathrm{tr}(X^2)} \, \mathrm{d}X .

The difference is that this density is restricted to the subspace of matrices such that X=XX=X^\ast. In this case all the eigenvalues lie on R\mathbb{R} and their joint density is given by

ϱn(x1,,xn)=Cnek=1nxi21i<jnxixj2.\varrho_n(x_1, \dots , x_n) = C^{\prime \prime}_n \mathrm{e}^{-\sum_{k=1}^n x_i^2} \prod_{1\leq i < j \leq n }|x_i - x_j|^2 .

This looks practically identical to the case of Ginibre, but in fact the restriction to R\mathbb{R} rather than C\mathbb{C} changes things considerably. As is well known, the 1-point density converges to a semicircular distribution. Furthermore, its rightmost eigenvalue asymptotically obeys a Tracy-Widom (1994) law. That is

P(maxi=1,,nλi2n+t2n16)exp(t(st)q(s)2ds)\mathbb{P}\left( \max_{i=1, \dots, n} \lambda_i \leq \sqrt{2n} + \frac{t}{\sqrt{2} n^\frac{1}{6}} \right) \to \exp\left(-\int_t^\infty (s-t) q(s)^2 \, \mathrm{d}s\right)

as nn \to \infty, and where qq satisfies Painlevé II,

q(t)=tq(t)+2q(t)3q^{\prime \prime}(t) = t q(t) + 2 q(t)^3

with boundary condition q(t)Ai(t)q(t) \sim \mathrm{Ai}(t) as t+t \to +\infty. This Tracy-Widom distribution is a highly universal object appearing in many different contexts such as random matrix theory, Ulam’s problem, the KPZ equation, asymmetric exclusion processes, Aztec diamonds etc. This result is also celebrated because it connects random matrix theory to integrable systems (though it was not the first to do so).

The Elliptic Ginibre Ensemble

We start with the following observation. A Ginibre matrix can be generated by

X=12(H1+iH2)X = \frac{1}{\sqrt{2}}(H_1 + i H_2)

where H1,H2H_1, H_2 are independently sampled GUE matrices. Motivated by this, consider the random matrix

X=1+τ2H1+i1τ2H2X = \sqrt{\frac{1+\tau}{2}}H_1 + i \sqrt{\frac{1-\tau}{2}} H_2

for τ[0,1]\tau \in [0,1]. τ=0\tau = 0 yields the Ginibre ensemble, τ=1\tau = 1 yields the GUE. Thus we can interpolate between these two ensembles. This is called the “elliptic Ginibre ensemble.” Its density is given by

PeGin(X)dX=Cne11τ2tr(XXτ(X2))dX.P_{\mathrm{eGin}}(X)\,\mathrm{d}X = C_n^{\prime \prime \prime} \mathrm{e}^{- \frac{1}{1-\tau^2}\mathrm{tr}\, (X X^\ast - \tau \Re (X^2))} \, \mathrm{d}X .

The ensemble derives its name because, up to a n\sqrt{n} scaling, the spectral density tends (as nn \to \infty) towards a constant on the ellipse

{(x,y)R2:x2(1+τ)2+y2(1τ)2<1}R2C\left\{ (x,y) \in \mathbb{R}^2 \, : \frac{x^2}{(1+\tau)^2} + \frac{y^2}{(1-\tau)^2} < 1 \right\} \subset \mathbb{R}^2 \simeq \mathbb{C}

and zero outside.

In order to see something interesting for the extremal eigenvalue we need to take a double scaling limit, where τn=1σ2n13\tau_n = 1 - \frac{\sigma^2}{n^\frac{1}{3}}, for σ>0\sigma > 0 a fixed parameter. This limit where the ellipse is very “flat” is called “weak non-Hermicity.” We are interested in the distribution of the largest real part in the limit. This scaling (due to Bender, 2010) is carefully chosen. Let τn\tau_n tend to 11 too fast and local correlations look like the GUE, too slow and local correlations look like Ginibre.

Results

Bender (2010) defines a rescaled point process at the edge of the ellipse. One zooms in on the edge of ellipse in a specified way.

xj:=λjcnan,x_j := \frac{\Re \lambda_j - c_n}{a_n}, yj:=λjbn.y_j := \frac{\Im \lambda_j}{b_n}.

The scalings an,bn,cna_n, b_n, c_n can be found in Bender’s paper, and I will just state the result: that the rescaled point process {zj:=(xj,yj)}j=1n\{ z_j := (x_j , y_j) \}_{j=1}^n converges as nn\to \infty. This yields a determinantal point process parametrised by σ>0\sigma > 0. We now arrive at our main result.

Theorem (Bothner-L.):

Fσ(t):=limnP(maxj=1,,nxjt)=exp(t(st)[Rpσ(s,y)2dνσ(y)]ds)\boxed{F_\sigma(t) := \lim_{n \to \infty} \mathbb{P}\left( \max_{j=1,\dots, n} x_j \leq t \right) = \exp\left( - \int_t^\infty (s-t) \left[ \int_\mathbb{R} p_\sigma(s,y)^2 \, \mathrm{d}\nu_\sigma(y) \right] \mathrm{d}s \right)}

where dνσ(λ)=1σπe(λσ)2dλ\mathrm{d}\nu_\sigma(\lambda) = \frac{1}{\sigma \sqrt{\pi}} \mathrm{e}^{-\left( \frac{\lambda}{\sigma}\right)^2}\, \mathrm{d}\lambda and where pσp_\sigma satifies integro-differential Painlevé II,

2t2pσ(t,y)=[t+y+2Rpσ(t,λ)2dνσ(λ)]pσ(t,y)\boxed{\frac{\partial^2}{\partial t^2} p_\sigma(t,y) = \left[ t+y+2\int_\mathbb{R} p_\sigma(t,\lambda)^2 \, \mathrm{d}\nu_\sigma(\lambda) \right] p_\sigma(t,y)}

with boundary condition pσ(t,y)Ai(t+y)p_\sigma(t,y) \sim \mathrm{Ai}(t+y) as t+t \to +\infty, pointwise in yRy \in \mathbb{R}. \triangle

How does one see that this generalises the Tracy-Widom result? If one takes σ0\sigma \downarrow 0 one expects to reduce to the GUE. One sees that dνσ(λ)δ(λ)dλ\mathrm{d}\nu_\sigma(\lambda) \to \delta(\lambda) \, \mathrm{d}\lambda, and so

limσ0Fσ(t)=exp(t(st)pσ(s,0)2ds)\lim_{\sigma \downarrow 0} F_\sigma(t) = \exp\left( - \int_t^\infty (s-t) p_\sigma(s,0)^2 \, \mathrm{d}s \right)

and our integro-differential equation, when evaluated at y=0y =0, reduces to Painlevé II, with the right boundary condition.

On the other hand, recovering the Gumbel distribution as σ\sigma \to \infty is actually harder and is done in our paper. This is best done by studying the asymptotics of the corresponding Fredholm determinant rather than working with the integro-differential equation. We also look at asymptotics of Fσ(t)F_\sigma(t) under various scalings of tt and σ\sigma by means of nonlinear steepest descent techniques.

Remark: There is work left to be done here, especially for the left tail tt \to -\infty, since we only obtain asymptotics in a scaling régime where σ\sigma is very small.

Sketch of proof

As the point process defined by Bender is determinantal, gap probabilities may be expressed in terms of Fredholm determinants.

Fσ(t)=det(1KAiσ)L2((t,)×R)F_\sigma(t) = \det(1- K^{\sigma}_{\mathrm{Ai}})_{L^2((t,\infty)\times \mathbb{R})}

The kernel of KAiσK^{\sigma}_{\mathrm{Ai}} (found by Bender, 2010) is complicated and is given in Equation 1.10 of our paper. We make the following observation. Let F:L2(R+×R)L2(R+×R)\mathcal{F} : L^2(\mathbb{R}_+ \times \mathbb{R}) \to L^2(\mathbb{R}_+ \times \mathbb{R})

Ff(x,η)=12πeiyηf(x,y)dy\mathcal{F}f(x,\eta) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty \mathrm{e}^{-i y \eta} f(x,y) \, \mathrm{d}y

be the Fourier transform in the vertical (imaginary) direction. Then we have

Kσ=FKAiσF1K_\sigma = \mathcal{F} K^{\sigma}_{\mathrm{Ai}} \mathcal{F}^{-1}

where KσK_{\sigma} is an integral operator with kernel

Kσ(z1,z2)=0ϕσ(x1+s,y1)ϕσ(x2+s,y2)dsK_{\sigma}(z_1, z_2) = \int_0^\infty \phi_\sigma(x_1 + s, y_1 ) \phi_\sigma(x_2 + s, y_2 ) \, \mathrm{d}s

where ϕ(x,y)=π14ey22Ai(x+σy)\phi(x,y) = \pi^{-\frac{1}{4}} \mathrm{e}^{- \frac{y^2}{2}}\mathrm{Ai}(x+\sigma y) and zi(xi,yi)z_i \equiv (x_i, y_i). By the invariance of the determinant, we then have

Fσ(t)=det(1Kσ)L2((t,)×R).F_\sigma(t) = \det(1- K_{\sigma})_{L^2((t,\infty)\times \mathbb{R})} .

Notice that Kσ(z1,z2)K_{\sigma}(z_1, z_2) takes the form of a Hankel composition, and so one can use the method given in my previous post to derive a Tracy-Widom-like result.

Finite temperature Airy processes

Our result is curious because the integro-differential equation we derive also appears in the study of finite temperature Airy processes on R\mathbb{R}, whereas we are looking at a point process in R2\mathbb{R}^2. What is the relationship?

Begin by transferring the tt dependence to the operator.

Fσ(t)=det(1Kσt)L2(R+×R)F_\sigma(t) = \det(1- K_{\sigma}^t)_{L^2(\mathbb{R}_+\times \mathbb{R})}

where

Kσt((x1,y1),(x2,y2))=Kσ((x1+t,y1),(x2+t,y2))K_{\sigma}^t((x_1, y_1),(x_2, y_2)) = K_{\sigma}((x_1+t, y_1),(x_2+t, y_2))

Now define the operator Pt,σ:L2(R+×R)L2(R+)P_{t,\sigma} : L^2(\mathbb{R}_+ \times \mathbb{R}) \to L^2(\mathbb{R}_+) with kernel

Pt,σ(s,(x,y))=π14e12y2Ai(s+x+σy+t).P_{t,\sigma}(s, (x,y)) = \pi^{-\frac{1}{4}} \mathrm{e}^{-\frac{1}{2}y^2} \mathrm{Ai}(s+x+\sigma y + t) .

Then a short calculation shows that

Pt,σPt,σ=KσtP_{t,\sigma}^\ast P_{t,\sigma} = K_\sigma^t

Now consider Pt,σPt,σ:L2(R+)L2(R+)P_{t,\sigma} P_{t,\sigma}^\ast : L^2(\mathbb{R}_+) \to L^2(\mathbb{R}_+). A short calculation shows that this operator has kernel

Pt,σPt,σ(s1,s2)=RΦ(yσ)Ai(s1+y+t)Ai(s2+y+t)dy=:Nσt(s1,s2)P_{t,\sigma} P_{t,\sigma}^\ast(s_1, s_2)=\int_\mathbb{R} \Phi\left( \frac{y}{\sigma} \right) \mathrm{Ai}(s_1+y+t) \mathrm{Ai}(s_2+y+t) \, \mathrm{d}y =: N_\sigma^t(s_1, s_2)

where Φ(x)=1πxet2dt\Phi(x) = \frac{1}{\sqrt{\pi}}\int_{-\infty}^x \mathrm{e}^{-t^2}\, \mathrm{d}t. Undoing the shift by letting Nσ(s1+t,s2+t)=Nσt(s1,s2)N_\sigma (s_1 +t , s_2 + t) = N_\sigma^t (s_1 , s_2) and using Sylvester’s identity we find

Fσ(t)=det(1Nσ)L2(t,).F_\sigma(t) = \det(1-N_\sigma)_{L^2(t,\infty)} .

This means that the largest particle of the finite temperature Airy process with function Φ\Phi is identically distributed to the largest real part in the elliptic edge process.

Remark: Indeed, this observation may be generalised. Let us consider any finite temperature Airy process, with kernel

N(x,y)=Rψ(s)Ai(x+s)Ai(y+s)dsN(x,y) = \int_\mathbb{R} \psi(s) \mathrm{Ai}(x+s) \mathrm{Ai}(y+s) \, \mathrm{d}s

where ψ(x)0\psi(x) \geq 0 is an increasing, continuously differentiable function such that ψ()=0\psi(-\infty) = 0 and ψ(+)=1\psi(+\infty) = 1. By an identical argument to that presented above we have the identity

det(1N)L2(t,)=det(1K)L2((t,)×R)\det(1-N)_{L^2(t,\infty)} = \det(1-K)_{L^2((t,\infty) \times \mathbb{R})}

where

K((x1,y1),(x2,y2))=0ϕ(x1+s,y1)ϕ(x2+s,y2)dsK((x_1,y_1), (x_2,y_2)) = \int_0^\infty \phi(x_1 +s, y_1) \phi(x_2 +s, y_2) \, \mathrm{d}s

where ϕ(x,y)=ψ(y)Ai(x+y)\phi(x,y) = \sqrt{\psi^\prime(y)}\mathrm{Ai}(x+y). Thus KK has the form of a Hankel composition and so the Fredholm determinant has an associated integro-differential representation (by the method of my previous post). That such Fredholm determinants have an integro-differential representation is already known in the literature but, to the author’s knowledge, this method of relating finite temperature kernels to Hankel compositions is new. \triangle

Remark: Here we make a remark that will be well known to experts, namely the relationship between finite temperature kernels and Riemann-Hilbert problems. It is clear that if we let

A:L2(R)L2(t,)A : L^2(\mathbb{R}) \to L^2(t,\infty)

with kernel

A(x,x)=Ai(x+x)ψ(x),A(x,x^\prime) = \mathrm{Ai}(x+x^\prime) \sqrt{\psi(x^\prime)} ,

then clearly AA=N:L2(t,)L2(t,)A A^\ast = N : L^2(t,\infty) \to L^2(t,\infty). A simple calculation shows that Mt:=AA:L2(R)L2(R)M_t := A^\ast A : L^2(\mathbb{R}) \to L^2(\mathbb{R}) has kernel

Mt(x,y)=ψ(x)ψ(y)Ai(x+t)Ai(y+t)Ai(x+t)Ai(y+t)xy.M_t (x,y) = \sqrt{\psi(x)}\sqrt{\psi(y)} \frac{\mathrm{Ai}(x+t)\mathrm{Ai}^\prime (y+t) - \mathrm{Ai}^\prime (x+t)\mathrm{Ai}(y+t)}{x-y} .

This is an integrable Its-Izergin-Korepin-Slavnov type operator (see this nice introduction to integrable operators by Deift). The IIKS theory relates Fredholm determinants of such operators to Riemann-Hilbert problems. Thus

det(1N)L2(t,)=det(1Mt)L2(R)\det(1-N)_{L^2(t,\infty)}= \det(1-M_t)_{L^2(\mathbb{R})}

may be related to a Riemann-Hilbert problem. This then allows for a Deift-Zhou steepest descent analysis.

Remark: As a final remark, these factorisations are useful at proving that the corresponding operator is trace class, since one need only show that the factors are respectively Hilbert-Schmidt.

The bulk of the spectrum (addendum)

The techniques developed above can be used to derive a similar integro-differential representation for gaps between real parts in the bulk of the spectrum (the reader is referred to our paper on the arXiv for details). In particular, the correlation kernel in the bulk also satisfies a curious factorisation property. Recall that we are considering the régime where τn1\tau_n \to 1 and so “the bulk” is approximately the set (2,2)(-2,2).

“Weak non-Hermiticity” in the bulk requires a somewhat different scaling than at the edge. In particular, if λ0(2,2)\lambda_0 \in (-2,2) is the point we’re going to “zoom in” on and ρ1(x)=1π(1x24)+\rho_1(x) = \frac{1}{\pi}\sqrt{\left(1-\frac{x^2}{4}\right)_+}, then we must take

τn=11n(σρ1(λ0))2\tau_n = 1- \frac{1}{n} \left( \frac{\sigma}{\rho_1(\lambda_0)}\right)^2

for σ0\sigma \geq 0. Then, after a suitable rescaling, the point process around the point λ0(2,2)\lambda_0 \in (-2,2) converges to a determinantal point process with kernel given by

Ksinσ(z1,z2)=1σπey12+y222σ212πππe(σu)2cos(u(z1z2))duK_{\mathrm{sin}}^{\sigma}(z_1,z_2) = \frac{1}{\sigma\sqrt{\pi}} \mathrm{e}^{-\frac{y_1^2 + y_2^2}{2\sigma^2}} \frac{1}{2\pi}\int_{-\pi}^\pi \mathrm{e}^{-(\sigma u)^2} \cos(u(z_1 - \overline{z_2})) \, \mathrm{d}u

where yk=zky_k = \Im z_k for k=1,2k=1,2. This kernel describes a determinantal point process in the plane R2C\mathbb{R}^2 \simeq \mathbb{C}.

Now let Jt=(t,t)×RR2J_t = (-t,t)\times \mathbb{R} \subset \mathbb{R}^2 and suppose we are interested in the gap probability given by

det(1Ksinσ)L2(Jt).\det(1-K_{\mathrm{sin}}^{\sigma})_{L^2(J_t)} .

This corresponds to looking at gaps between real parts (note also that the point process is horizontally translation invariant). We now observe the following factorisation. Let

Aσ:L2(Jt)L2(π,π)A_\sigma : L^2(J_t) \to L^2(-\pi,\pi)

with kernel

Aσ(a,z)=12π32σexp(y22σ212(σa)2iaz).A_\sigma(a,z) = \frac{1}{\sqrt{2\pi^\frac{3}{2}\sigma}} \exp\left(-\frac{y^2}{2\sigma^2} - \frac{1}{2}(\sigma a)^2 - ia \overline{z}\right) .

A simple calculation then shows that

Ksinσ=AσAσ:L2(Jt)L2(Jt)K_{\mathrm{sin}}^{\sigma} = A_\sigma^\ast A_\sigma : L^2(J_t) \to L^2(J_t)

But then, by Sylvester’s identity, we have

det(1Ksinσ)L2(Jt)=det(1AσAσ)L2(π,π)=det(1Sσt)L2(t,t)\det(1-K_{\mathrm{sin}}^{\sigma})_{L^2(J_t)} = \det(1-A_\sigma A_\sigma^\ast)_{L^2(-\pi,\pi)} = \det(1-S_\sigma^t )_{L^2(-t,t)}

where SσtS_\sigma^t is a rescaled version of AσAσA_\sigma A_\sigma^\ast (so that its domain is now L2(t,t)L^2(-t,t)). A simple calculation (see our paper) shows that one can write the kernel of SσtS_\sigma^t as

Sσt(a,b)=0[Φ(tσ(z+1))Φ(tσ(z1))]cos(π(ab)z)dz.S_\sigma^t(a,b) = \int_0^\infty \left[ \Phi\left( \frac{t}{\sigma}(z+1) \right) - \Phi\left( \frac{t}{\sigma}(z-1)\right) \right]\cos(\pi(a-b)z) \, \mathrm{d}z .

This is exactly a finite temperature sine kernel. In our paper we show that for any determinantal point process on R\mathbb{R} with kernel of the form

K(a,b)=0w(z)cos(π(ab)z)dzK(a,b) = \int_0^\infty w(z)\cos(\pi(a-b)z) \, \mathrm{d}z

for w:R+[0,1)w: \mathbb{R}_+ \to [0,1) and tending to 00 exponentially fast at ++\infty, the gap probability det(1K)L2(t,t)\det(1-K)_{L^2(-t,t)} can be represented in terms of a solution to an integro-differential Painlevé V equation. The generalises the famous result of Jimbo-Miwa-Môri-Sato (1980).