On the extreme value statistics of normal random matrices and 2D Coulomb gases: Universality and finite N corrections

TL;DR

This paper extends orthogonal polynomial methods to analyze extreme eigenvalue statistics of normal random matrices and 2D Coulomb gases, demonstrating universality, finite N corrections, and connections to Hermitian matrix techniques.

Contribution

It introduces a generalized approach for extreme value analysis of non-Hermitian matrices, including finite N corrections and universality results, using methods adapted from Hermitian matrix theory.

Findings

Eigenvalue with largest modulus converges to Gumbel distribution

Universality holds for arbitrary radially symmetric potentials

Finite N corrections are explicitly computed

Abstract

In this paper we extend the orthogonal polynomials approach for extreme value calculations of Hermitian random matrices, developed by Nadal and Majumdar [1102.0738], to normal random matrices and 2D Coulomb gases in general. Firstly, we show that this approach provides an alternative derivation of results in the literature. More precisely, we show convergence of the rescaled eigenvalue with largest modulus of a normal Gaussian ensemble to a Gumbel distribution, as well as universality for an arbitrary radially symmetric potential. Secondly, it is shown that this approach can be generalised to obtain convergence of the eigenvalue with smallest modulus and its universality for ring distributions. Most interestingly, the here presented techniques are used to compute all slowly varying finite N correction of the above distributions, which is important for practical applications, given the slow convergence. Another interesting aspect of this work is the fact that we can use standard techniques from Hermitian random matrices to obtain the extreme value statistics of non-Hermitian random matrices resembling the large N expansion used in context of the double scaling limit of Hermitian matrix models in string theory.

Equations

$$\mu (z):=\underset{N\to \infty }{\lim }{\mu }_N(z)=\frac{1}{\pi }{\mathbb{I}}_{\{|z|\le 1\}}.$$

$${\mathbb{P}}_N(M)\mathrm{d}M=\frac{1}{N!Z_N}{e}^{-N\mathrm{tr}V(M)}\mathrm{d}M,$$

$$Z_N=\frac{1}{N!}\int e^{-N\mathrm{tr}V(M)}\mathrm{d}M$$

$${\mathbb{P}}_N(z_1,\dots ,z_N)=\frac{1}{N!Z_N}\prod _{1\le j<k\le N}|z_j-z_k{|}^2{e}^{-N{\sum }_{i=1}^{N}V(z_i)}.$$

$${\mathbb{P}}_N(z_1,\dots ,z_N)=\frac{e^{-N^2{S}_{\mathrm{eff}}(z_1,\dots ,z_N)}}{N!Z_N}$$

$$S_{\mathrm{eff}}(z_1,\dots ,z_N)=\frac{1}{N}\sum _{i=1}^{N}V(z_i)-\frac{2}{N^2}\sum _{1\le j<k\le N}\log |z_j-z_k|$$

$$V(z)\equiv V(r),z=r{e}^{i\varphi }.$$

$$\underset{r\to \infty }{\lim }\left(V(r)-\log r^2\right)=\infty$$

$$\mu(z)\,\mathrm{d}^{2}z=\frac{1}{4\pi}\,\frac{\mathrm{d}}{\mathrm{d}r}\Big{(}r\,V^{\prime}(r)\Big{)}\,\mathbb{I}_{\{a_{-}\leq r\leq a_{+}\}}\,\mathrm{d}r\,\mathrm{d}\phi.$$

$$V^{\prime }(a_{-})=0$$

$$a_{+}{V}^{\prime }(a_{+})=2.$$

$$⟨p_n|p_m{⟩}_w:=\int w(z)p_n(z)\overline{p_m(z)}{\mathrm{d}}^{2}z=h_n{\delta }_{nm}$$

$$Z_N=\prod _{n=0}^{N-1}{h}_n,$$

$$p_n(z)=z^n.$$

$$h_n=2\pi {\int }_0^{\infty }{r}^{2n+1}w(r)\mathrm{d}r=2\pi {\int }_0^{\infty }{r}^{2n+1}{e}^{-NV(r)}\mathrm{d}r.$$

$$Z_N(\mathcal{D})=\frac{1}{N!}\int \prod _{1\le j<k\le N}|z_j-z_k{|}^2{e}^{-N{\sum }_{i=1}^{N}V(z_i)}\prod _j{\mathbb{I}}_{z_j\in \mathcal{D}}{\mathrm{d}}^2{z}_j,$$

$$\mathbb{P}_{N}(\mbox{all eigenvalues in $\mathcal{D}$})=\frac{Z_{N}(\mathcal{D})}{Z_{N}}.$$

$$F_N(y):={\mathbb{P}}_N(|z_{max}|\le y)={\mathbb{P}}_N(|z_1|\le y,\dots ,|z_N|\le y)=\frac{Z_N(y)}{Z_N(\infty )}.$$

$$\int e^{-NV(z)}{\mathbb{I}}_{\{|z|\le y\}}{p}_n(z;y)\overline{p_m(z;y)}{\mathrm{d}}^{2}z=h_n(y){\delta }_{nm}.$$

$$Z_N(y)=\prod _{n=0}^{N-1}{h}_n(y)$$

$$F_N(y)={\mathbb{P}}_N(|z_{max}|\le y)=\prod _{n=0}^{N-1}\frac{h_n(y)}{h_n(\infty )}.$$

$$h_n(y)=2\pi {\int }_0^y{r}^{2n+1}{e}^{-NV(r)}\mathrm{d}r,$$

$$y=a_{+}+C{N}^{-1/2}({\alpha }_N+{\beta }_{N}Y),$$

$$F_N(y)$$

$$h_n(y)=2\pi {\int }_0^y{e}^{Nf(r)}\mathrm{d}r,\text{mbox}withf(r)=\frac{2n+1}{N}\log r-r^2.$$

$$f(r)=f(r_{0})-\frac{1}{2}\big{|}f^{\prime\prime}(r_{0})\big{|}\,(r-r_{0})^{2}+\mathcal{O}\left((r-r_{0})^{3}\right),$$

$$\frac{h_n(y)}{h_n(\infty )}=\frac{\mathrm{erf}(\sqrt{Nb}{r}_0)+\mathrm{erf}(\sqrt{Nb}(y-r_0))}{1+\mathrm{erf}(\sqrt{Nb}{r}_0)}+...,$$

$$\mathrm{log}\big{(}F_{N}(y)\big{)}\!=\!N\!\int_{0}^{1}g(\xi;y)\,\mathrm{d}\xi+\frac{1}{2}\big{[}g(1;y)-g(0;y)\big{]}\!+\!\mathcal{O}\!\left(N^{-1}\right)\!,$$

$$g(\xi ;y)$$

$$\mathrm{erf}(z)\equiv \frac{2}{\sqrt{\pi }}{\int }_0^z{e}^{-x^2}\mathrm{d}x=1-\frac{e^{-z^2}}{\sqrt{\pi }z}+\mathcal{O}\left(z^{-3}{e}^{-z^2}\right).$$

Full text

1 Introduction

Complex non-Hermitian random matrices were first investigated by Ginibre [1], who analysed matrices whose elements are independent and identically distributed (real or complex) Gaussian variables. He showed that the eigenvalue density $\mu_{N}(z)$ of such $N$ dimensional random matrices in the limit $N\to\infty$ is given by the well-known circular law 111The eigenvalue density is formally defined as $\int_{A}\mu_{N}(z)\,\mathrm{d}^{2}z=\frac{1}{N}\,\mathbb{E}_{N}(\#\{\mbox{eigenvalues in A}\})$ for any $A\subset\mathbb{C}$, where $\mathbb{E}_{N}(\cdot)$ indicates the expectation value with respect to the joint probability density function $\mathbb{P}_{N}$ over the $N$ eigenvalues.

[TABLE]

In other words, the eigenvalues are distributed uniformly over the disk of radius 1 in the complex plane (see Figure 1 (Left)).

In recent years there has been an increasing interest in non-Hermitian random matrices (see e.g. [2, 3, 4, 5, 6] for some examples as well as Chapter 18 of [7] for a review). Random matrices where all elements are chosen independently from some distribution are generally referred to as Wigner matrices. The Gaussian case discussed above is a special case, since it allows for a Coulomb gas formulation describing the joint eigenvalues distribution in terms of an attractive Gaussian potential and a logarithmic, Coulomb repulsion [1]. In general, Wigner random matrices do not allow for such a formulation due to the lack of symmetries. One closely related ensemble is that of normal random matrices [8] where one introduces an additional symmetry that the matrices are normal (they commute with their adjoint). This symmetry allows for a Coulomb gas formulation with arbitrary potential. Another important class of non-Hermitian random matrices are isotropic or bi-unitarily invariant random matrices (see for example [9, 10]).

Normal random matrices have been studied intensively during the last years [11, 12, 13, 14]. In this work we are interested in the extreme value statistics of normal random matrices with general radially symmetric potential. As the eigenvalues can spread over the complex plane (within the support), they are less correlated than in the unitary or orthogonal case. This is reflected by the extreme value statistics. In the above example, considering the distribution of the eigenvalue with the largest modulus, it is known [15, 16] that the rescaled distribution is given by a Gumbel distribution. Recall that the Gumbel distribution is one of three families which describe the extreme value statistics of independent and identically distributed random variables (Weilbull, Gumbel, Fréchet) [17, 18, 19]. While the eigenvalues of non-Hermitian ensembles are not independent, their correlation is much weaker than for Hermitian ensembles. Indeed for Hermitian ensembles the eigenvalues are strongly correlated and the distribution of the largest eigenvalue is given by the Tracy-Widom distribution [20, 21, 22].

Outline and what’s new.

After a brief introduction to normal matrices and their relation to 2D Coulomb gases in Section 2, we extend in Section 3 the orthogonal polynomials approach to extreme value statistics of random matrices from the Hermitian case [23] to non-Hermitian random matrices. In Section 4 we use this approach to show convergence of the probability distribution of the rescaled eigenvalue with largest modulus of normal random matrices to a Gumbel distribution providing a simplified, alternative derivation of result in [15, 16] for Coulomb gases. In Section 5 we show how the here presented approach can be used to calculate finite $N$ corrections to the Gumbel distribution. In particular, we provide an analytical expression including all slowly varying correction terms. In Section 6 and 7 we show that our approach can also be used to show convergence and universality of the distribution of the eigenvalue with largest modulus at the outer edge as well as the eigenvalue with smallest modulus at the inner edge of a ring distribution. In Section 8 we derive finite $N$ results, analogous to those derived in Section 5 for the Gaussian, for the extreme value statistics at the outer and inner edge of arbitrary radially symmetric potentials. We conclude in Section 9 where we also comment on future directions of research. At various places in this article analytical results are complemented with numerical results. The methods to obtain those numerical results are explained in A.

Note that while the results in Section 4 and 6 are an alternative derivation of previously proven theorems [15] and [16], the here presented derivation is novel and interesting for several reasons. Firstly, we can use the same orthogonal polynomial approach as used for Hermitian random matrices in the case of non-Hermitian random matrices. Those techniques are widely used in the field of random matrix theory, thus making the derivation easily accessible. Secondly, it is non-trivial that we can use the same derivational steps as used to derive the Tracy-Widom distribution in [23] to show convergence of the extreme value statistics of normal random matrices to the Gumbel distribution. Furthermore, the derivation resembles the double scaling limit of the free energy in Hermitian matrix models [24].

The results of Section 5, 7 and 8 are new. Particularly interesting and important for practical applications is the analytical expression derived in Section 5 and 8 for the finite $N$ contributions to the extreme value statistics including the Gumbel distribution and all slowly varying correction terms. Regarding the results for the inner edge, Section 7, it is interesting that one can use the same general approach as in Sections 4 and 6. However, alternatively one could have also used the mapping of eigenvalues $z_{i}\to 1/z_{i}$ and the corresponding change in the potential, $V(z)\to V(1/z)+2(1+1/N)\log|z|$, to derive the results of Section 7 from Section 6 or the corresponding result in the literature [16].

2 Normal random matrices and Coulomb gases

Normal random matrices were introduced in [8] and are defined through the measure

[TABLE]

where $M$ are complex $N\times N$ dimensional matrices satisfying the constraint $[M,M^{\dagger}]=0$. Here $\mathrm{d}M=\prod_{ij}\mathrm{d}\Re M_{ij}\,\mathrm{d}\Im M_{ij}$ is the so-called Lebesgue measure and the normalisation

[TABLE]

is also referred to as the partition function.

The eigenvalue statistics of normal random matrices is in many ways similar to that of general non-Hermitian matrices. However, normal random matrices are computationally much more feasible due to their symmetry. In particular, the normality condition $[M,M^{\dagger}]=0$ implies that there exists a unitary matrix which simultaneously diagonalises $M$ and $M^{\dagger}$. This allows one to employ Coulomb gas techniques and find the joint probability measure of the complex eigenvalues, given by

[TABLE]

Here the term $\Delta(\{z\}):=\prod_{1\leq j<k\leq N}(z_{j}-z_{k})$ is called the Vandermonde determinant. Writing

[TABLE]

with

[TABLE]

we can interpret the effective model as that of a gas of $N$ particles in a potential $V(z)$ experiencing a logarithmic repulsion. This is equivalent to a Coulomb repulsion in 2D, thus the name Coulomb gas.

As mentioned above, the (complex or real) Ginibre ensemble where all entries are (complex or real) Gaussian random variables also has a joint probability density function of its eigenvalues given by (4) with Gaussian potential. It is interesting to study normal random matrices as their joint eigenvalue distribution is described by a Coulomb gas with general potential $V(z)$.

In this article we are only interested in radially symmetric potentials

[TABLE]

In this case one can use saddle point techniques or more mathematical results from the theory of harmonic analysis [25] to find the eigenvalue density for a given radially symmetric potential. More precisely, given the condition that

[TABLE]

and $r\,V^{\prime}(r)$ increasing in $\mathbb{R}^{+}$ or $V^{\prime}(r)>0$ and $V$ convex in $\mathbb{R}^{+}$ one finds that the eigenvalue density is given by [25]

[TABLE]

We see that the eigenvalue density has support on a ring with inner radius $a_{-}$ and outer radius $a_{+}$ (or a disk in the case where $a_{-}=0$), where the radii are given by

[TABLE]

and $a_{+}$ is the solution to following equation

[TABLE]

We see that the eigenvalue density either has topology of a disk or of an annulus. This is the famous one ring theorem [26, 27, 25, 6].

As an example, for the Gaussian normal ensemble with potential $V(r)=r^{2}$, we have $a_{-}=0$ and $a_{+}=1$, and the eigenvalues spread over the unit disk in the complex plane. We see that the eigenvalue density is identical to that of the Ginibre ensemble, i.e., the circular law (see Figure 1 (Left)). Another easy but more interesting example is the potential $V(r)=r^{2}+sr$. Here for $s>0$ the topology of the support is a disk while for $s<0$ the topology is that of an annulus (see Figure 1 (Right)). Interesting statistics can be computed for the transition when the topology changes from disk to annulus (see for example [28]).

Besides saddle point analysis, another powerful technique employed in random matrix theory in general and in this article in particular is the method of orthogonal polynomials. The aim of the orthogonal polynomials approach is to define a set of polynomials which are orthogonal with respect to a weight $w(z)$, i.e.,

[TABLE]

where in our case $w(z)=e^{-N\,V(z)}$. Having found such a set of polynomials one can use them to evaluate many quantities in random matrix theory, in particular, the partition function can be obtained as

[TABLE]

where $h_{n}$ are the normalisation constants of the orthogonal polynomials as defined in (12).

For the case of radially symmetric potentials, the orthogonal polynomials become very simple. In fact, for any weight $w(z)=w(|z|)$ it can be readily checked that they are simply given by the monomials

[TABLE]

Thus one has

[TABLE]

For the case of the Gaussian potential with $V(r)=r^{2}$ the normalisation constants are given by $h_{n}=\pi N^{-n-1}\Gamma(n+1)$.

3 Orthogonal polynomial approach for extreme value statistics of random matrices

One of the most well-known examples of extreme value statistics of random matrices is the Tracy-Widom distribution [20, 21, 22] originally obtained for the GUE. The orthogonal polynomial approach to extreme value statistics of random matrices was first introduced in [23] to provide an easy derivation of the Tracy-Widom distribution for the GUE. It has then been extended to arbitrary potentials and multi-critical points [29] as well as large deviations [30]. Here we show how one can use the same formulation for normal random matrices with complex eigenvalues. More precisely we are interested in the statistics of the complex eigenvalue $z_{max}$ with the largest modulus. The basic idea underlying this approach is simple. Let us define a partition function $Z_{N}(\mathcal{D})$ where the eigenvalues are restricted to a domain $\mathcal{D}\subseteq\mathbb{C}$,

[TABLE]

where $\mathbb{I}_{\,C}$ is the indicator function which is 1 if $C$ is true and 0 otherwise. Then $Z_{N}(\mathbb{C})\equiv Z_{N}$, i.e., when integrating over the entire domain we recover the conventional partition function. Given the above, it follows that

[TABLE]

We can use this relation to express extreme value statistics in terms of $Z_{N}(\mathcal{D})$. More precisely, let us choose $\mathcal{D}=\{z:\,|z|\leq y\}$ and denote as shorthand $Z_{N}(y):=Z_{N}(\{z:\,|z|\leq y\})$ (and thus $Z_{N}(\infty)=Z_{N}$), then the cumulative probability function of the eigenvalue with the largest modulus is given by

[TABLE]

Following the analogous approach for the GUE [23], the key idea is now to define orthogonal polynomials for the partition function $Z_{N}(y)$. In particular, we can find the orthogonal polynomials such that

[TABLE]

We have

[TABLE]

and thus

[TABLE]

As above, we are interested in the case of radially symmetric potentials $V(z)=V(|z|)$. In this case, we observe that (19) is equivalent to (12) with a weight $w(z)\equiv w(|z|;y)=e^{-N\,V(|z|)}\,\mathbb{I}_{\{|z|\leq y\}}$. Since $w(z)$ is radially symmetric, we know from the previous section that the orthogonal polynomials are simply given by the monomials $p_{n}(z;y):=z^{n}$. We see that, while in general $p_{n}(z;y)$ is a function of $y$, here all $y$-dependence is in the normalization $h_{n}(y)$. Therefore, using equation (15) we obtain

[TABLE]

where the indicator function in the weight reduces the integration range of $r$ to $[0,y]$.

Once $N$ becomes large the probability density $f_{N}(y)=F_{N}^{\prime}(y)$ will be more and more centered around the edge of the support, $a_{+}$. More precisely, one expects $(|z_{max}|-a_{+})\sim N^{-1/2}$.222Note that $N^{-1/2}$ is the natural distance scaling for $N$ eigenvalues in a finite 2D domain. However, for the Unitary ensembles it is known that the scaling changes at the edge. This is not expected here as the correlation is weaker (which is the reason why one obtains the Gumbel distribution). Mathematically, one can see that this scaling is indeed correct from the terms $\sqrt{N}\,(y-r_{0})$ in later equations, cf. (27). It is however interesting to consider so-called multi-critical models where more derivatives in (26) become zero, potentially leading to different scalings and universality classes, cf. [29] for corresponding results for Unitary ensembles. Thus to obtain convergence to a limiting distribution one has to rescale $y$ to a “continuum” variable $Y$,

[TABLE]

where $C$ is a constant and $\alpha_{N},\beta_{N}$ are slowly varying functions.333Recall that a function $f(N)$ is called slowly varying (at infinity) if $\lim_{N\to\infty}f(aN)/f(N)=1$ for all $a>0$. A typical example is $\log N$. Loosely speaking a slowly varying function grows slower than any power of $N$. The constant $C$ is chosen merely for convenience and could equally be absorbed in $\alpha_{N}$ and $\beta_{N}$. Since $y\geq 0$, we have $Y\geq-\alpha_{N}/\beta_{N}$. For the limiting distribution to have support on $Y\in(-\infty,+\infty)$, we thus require $\alpha_{N}/\beta_{N}\to\infty$ as $N\to\infty$. This is enough to obtain the extreme value statistics related to the eigenvalue with the largest modulus for the Gaussian potential, as we will show in the next section.

4 The Gaussian potential

From the results of the previous section, more precisely (21) and (22), it follows that the cumulative probability function of the eigenvalue with largest modulus for the Gaussian random normal matrix is given by

[TABLE]

in which $P(a,z):=\gamma(a,z)/\Gamma(a)$ is the regularised incomplete gamma function. Note that this result can also be obtained from [31], where it was shown that the moduli of eigenvalues for a Gaussian potential follow the distribution of independent and identically distributed $\chi^{2}$ random variables.

One can try to reduce the large $N$ behaviour of (24) using the asymptotic expansion of the regularised incomplete gamma function, where both arguments become large, however, it is more instructive to calculate the asymptotic behaviour from first principles using a saddle point evaluation. In particular, such a derivation will generalise to arbitrary potentials, where we cannot find a closed form expression such as (24).

In the following we determine the large $N$ expansion of $h_{n}(y)$ using a saddle point expansion. Note that we can also obtain it from the so-called “string equation”, a set of equations describing a recursion relation for the coefficients of the orthogonal polynomials, which is widely used in Hermitian matrix models [24].

To perform the asymptotic expansion using a saddle point evaluation we start by writing

[TABLE]

We can now expand $f(r)$ around its maximum at $r_{0}=\sqrt{(2n+1)/(2N)}$ as

[TABLE]

The integration in (22) then reduces to Gaussian integrals over finite intervals, yielding error functions:

[TABLE]

where b\equiv\big{|}f^{\prime\prime}(r_{0})/2\big{|}=2 for the Gaussian potential. Here “$+...$” refers to the corrections coming from the $\mathcal{O}\left((r-r_{0})^{3}\right)$ term in (26). That those corrections are indeed of lower order can be verified later.

Now using the Euler-Maclaurin summation formula, we look at the first two leading order terms in (21) when both $N$ and $n$ become large, while their ratio is given by $\frac{n}{N}=:\xi\in[0,1]$. One obtains

[TABLE]

where we defined

[TABLE]

in which $r_{0}(\xi)=\sqrt{\xi}+\mathcal{O}\left(N^{-1}\right)$. Recall the asymptotic expansion of the error function

[TABLE]

Note that the argument of the first error function in the denominator of (29) is not necessarily large if $\xi$ is close to zero. However, the asymptotics of (29) is dominated by the second term

[TABLE]

Note that the scaling (23) ensures that the argument of the second error function \sqrt{Nb}\,\big{(}y-r_{0}(\xi)\big{)}\to\infty as $N\to\infty$ as long as the slowly varying function $\alpha_{N}\to\infty$. Furthermore, since $r_{0}(\xi)\to a_{+}$ as $\xi\to 1$, the dominant contribution in (31) comes from $\xi$ close to one, in which case we have

[TABLE]

Therefore, the leading term on the right hand sight of (28) is given by

[TABLE]

Let us now insert the scaling relation (23). We know that $\alpha_{N}$ is a slowly varying function which diverges for $N\to\infty$. Furthermore, we require $\alpha_{N}/\beta_{N}\to\infty$ as $N\to\infty$. It will become clear later that we will have to choose $\beta_{N}\sim 1/\alpha_{N}$ and moreover that it will be convenient to choose $\beta_{N}=1/(2\alpha_{N})$. The factor of $2$ is merely a rescaling of $Y$ which could otherwise be determined at the end of the calculation. We thus have the scaling

[TABLE]

where here $a_{+}=1$. Since we know that the contribution of $g(\xi;y)$ comes from $\xi$ close to one, we introduce the following change of variable in the integration in (33),

[TABLE]

For this choice one has, recalling that $b=2$,

[TABLE]

This relation also justifies the above choices of scalings. In particular,

[TABLE]

is a slowly varying function plus a finite term. To get this it was necessary to choose $\beta_{N}=1/(2\alpha_{N})$. Putting everything together, one then obtains

[TABLE]

For the pre-factor right-hand-side to converge to be 1, one is led to the following choice for the slowly varying function $\alpha_{N}$, to leading order,

[TABLE]

Using this, one arrives at

[TABLE]

Furthermore, under this choice, for the subleading term on the right hand side of (28) we have

[TABLE]

In addition, one can see that the corrections terms to (27) are of the same order.

The final result is thus given by

[TABLE]

which is precisely the Gumbel distribution. This provides an alternative derivation of the results in [15].

Figure 2 shows $F_{N}(Y)$ evaluated for increasing values of $N$. It is clear that for any practical applications, the limiting distribution is not a good approximation. However, the result presented in (24) provides a closed expression for $F_{N}(Y)$ for finite $N$ in terms of the regularised incomplete gamma function which can easily be evaluated numerically by inserting the slowly varying function (39) into (24), as had been done in Figure 2. For large values of $N$, one can benefit from the fact that the sum in (24) is dominated by $n\approx N$. This is equivalent to the statement that the major contribution to the integral (33) comes from the terms with $\xi$ close to one, as given by (35). We can also verify the correctness of this expression by comparing the result of $F_{N}(Y)$ with the empirical cumulative probability function obtained from sampling $m$ representations of a Ginibre matrix of size $N$. This is shown in Figure 3 for $m=1000$, with $N=200$ and $N=1000$ respectively.

5 Finite $N$ corrections for the Gaussian potential

One interesting aspect of the here presented approach is that it follows closely the standard finite $N$ expansion of the free energy of Hermitian random matrices using orthogonal polynomials (see for example Section 2.3 of [24]). This provides a well-established framework for calculating finite $N$ corrections as we discuss in this section. Furthermore, besides this “perturbative” expansion, “non-perturbative” expansions, so-called instantons, can in principle be used to obtain large deviations (for a corresponding result for Unitary ensembles see [30]).

Calculating the next-to-leading order terms of $F_{N}(Y)$ can be important in practice since the convergence to the Gumbel distribution is very slow as can be seen from the order of the subleading term in (40), manifested in Figure 2. We now compute the next-to-leading order terms in (40) and (42) up to, but not including terms which are of order $1/\sqrt{N}$ and potential logarithmic factors (which is the first correction that is not slowly varying). Those are the most important corrections, since those logarithmic terms only fall off very slowly. As we will see later, there is an infinite number of such slowly varying correction terms.

Before presenting an explicit expression for all logarithmic correction terms, we outline the general procedure to obtain any higher-order correction. Firstly, we can extend (28) to any order using the Euler-MacLaurin formula

[TABLE]

where $B_{p}$ are the Bernoulli numbers and

[TABLE]

Using that $f(r)=\frac{2\xi N+1}{N}\;\log r-V(r)$ (for general potential $V(r)$) and demanding that $f^{\prime}(r_{0})=0$ which determines $r_{0}(\xi)$, we can expand

[TABLE]

From this we can perform a Feynman type expansion

[TABLE]

Each integral is Gaussian (given that $f^{\prime\prime}(r_{0})<0$) and results in a term involving an error function which has an expansion

[TABLE]

Inserting back each expansion in the previous one and carefully collecting terms of equal order one can obtain the general large $N$ expansion of $\log F_{N}(y)$.

Here we are only interested in logarithmic corrections for the Gaussian potential. We have already seen in the previous section that the second term in (43) is of order $\mathcal{O}(1/\sqrt{N\log N})$, given our choice of $y$, i.e. (34). This is also true for terms coming from the Taylor expansion of $Nf(r)$, i.e. (45), beyond the quadratic term. Thus we have that 444For the remainder of this section $+...$ refers to correction terms of the order of $\mathcal{O}(N^{-1/2})$ times some logarithmic function. In principle, we would be happy in ignoring anything which is going to zero faster than any slowly varying term, i.e. terms of $\mathcal{O}(N^{-\epsilon})$, for any $\epsilon>0$. It just happens that the most dominant of such terms is $\mathcal{O}(N^{-1/2})$ times some logarithmic function.

[TABLE]

with $g(\xi;y)$ given in (29). For $y$ given in (34), $g(\xi;y)$ is dominated by the second error function in the numerator of (29) which can be used to obtain, for the Gaussian potential,

[TABLE]

where we expanded the logarithm to first order using the fact that higher order terms from the expansion of the logarithm in (29) beyond the linear term go to zero at least as fast as $\mathcal{O}(1/\sqrt{N})$. This can furthermore, be written as:

[TABLE]

Using a change of variable z:=\sqrt{2N}\,\big{(}y-r_{0}(\xi)\big{)} and the explicit form of $r_{0}(\xi)=\sqrt{\xi}+\mathcal{O}(1/N)$ for the Gaussian potential we can write

[TABLE]

In fact, one can integrate this expression which yields555Note that the leading order term cancels out.

[TABLE]

where we substituted in the scaling relation (34) for $y$ and defined

[TABLE]

with $\alpha_{N}$ is given in (39). We observe that (52) still contains terms of order $\mathcal{O}(N^{-1/2})$, dropping those we get

[TABLE]

with

[TABLE]

This provides an analytic expression including the leading order term of $\log F_{N}$, the Gumbel distribution, and all slowly varying correction terms! Figure 4 shows a plot of $F_{N}(Y)$ versus $\exp(\phi_{N}(Y))$ for $N=10$ and $N=100$. We see that already for $N=100$ the analytical expression $\exp(\phi_{N}(Y))$ matches the numerically evaluated $F_{N}(Y)$ perfectly.

To see that (55) is indeed an infinite series of slowly varying terms in $N$, we use the asymptotic expansion of the complementary error function

[TABLE]

and the definition of $\alpha_{N}$, given in (39), to write

[TABLE]

This highlights that (55) indeed constitutes the leading order term plus an infinite number of slowly varying correction terms.

6 Universality at the outer edge

In this section we consider the generalisation of the Gaussian potential to an arbitrary radially symmetric potential, $V\equiv V(r)$. The approach presented in Section 4 translates directly to an arbitrary potential by simply replacing $f(r)$ in (25) with

[TABLE]

The local extremum $r_{0}=r_{0}(\xi)$ determined by $f^{\prime}(r_{0})=0$ is given by the solution to the following equation

[TABLE]

Comparing (59) and (11) shows that in the large $N$ limit, the outer edge and $r_{0}$ coincide, i.e., $r_{0}(1)=a_{+}+\mathcal{O}(N^{-1})$.

Next we Taylor expand $f(r)$ around its extremum at $r_{0}$. To ensure that the extremum is indeed a maximum we compute the second derivative,

[TABLE]

Therefore, we have $f^{\prime\prime}(r_{0})<0$ and thus have a maximum, if the condition

[TABLE]

or equivalently

[TABLE]

holds for the potential $V(r)$ given the condition (8). Note that these conditions are the same as the ones given below (8).

Since we have shown that $f^{\prime\prime}(r_{0})<0$ under the condition (61) or (62), it then follows that (27) and (32) from the Gaussian case remain valid, where now $f^{\prime\prime}(r_{0})$ in $b=|f^{\prime\prime}(r_{0})/2|$ is given by (60). Furthermore, using a Taylor expansion around $\xi=1$ we get to leading order,

[TABLE]

Therefore, the leading order term in the probability distribution of the eigenvalue with the largest modulus is given by

[TABLE]

As in the Gaussian case, we use the general scaling ansatz for $y$ given by (34), i.e.

[TABLE]

where we replaced $b$ by (63) and the slowly varying function $\alpha_{N}^{(+)}$ will be determined later.

We thus also change the integration variable $\xi$ in (64) as in (35) only with an additional constant $\gamma_{+}>0$,

[TABLE]

The value of $\gamma_{+}$ will be determined in the following. We expand $r_{0}(\xi)$ around $\xi=1$:

[TABLE]

where $\delta_{+}$ is given by

[TABLE]

Hence, according to the constraint already imposed in (61), we have $\delta_{+}>0$ and therefore obtain

[TABLE]

which by choosing

[TABLE]

reduces to (36). Thus we have,

[TABLE]

The pre-factor is one to leading order if we choose,

[TABLE]

For the Gaussian potential the last term in the above expression is zero.

We have thus shown universality of the distribution of the eigenvalue with the largest modulus when rescaled around the outer edge of a generic potential $V(r)$ satisfying the condition given in (61) or (62), where the limiting probability distribution is the Gumbel distribution

[TABLE]

with $\alpha_{N}^{(+)}$ given in (72) and $F_{+}$ in (63). This reproduces results from [16] using our general approach.

7 Universality at the inner edge

Now we study the extreme value statistics of a normal matrix ensemble with potential $V=V(r)$ at the inner edge of the eigenvalue support for potentials where $a_{-}>0$. To do so we now choose $\mathcal{D}$ in (16) to be $\mathcal{D}=\{z:\,|z|\geq y\}$ and define

[TABLE]

The probability distribution for the eigenvalue with the smallest modulus $z_{min}$ is then given by

[TABLE]

where now $h_{n}(y)$ is given through

[TABLE]

We want to evaluate the above expression and scale it around the finite inner edge of the support at $|z|=a_{-}>0$. We have

[TABLE]

with $f(r)=(2n+1)/N\;\log r-V(r)$. Then performing a Taylor expansion and following a similar strategy as in the previous section results in

[TABLE]

where again b\equiv\big{|}f^{\prime\prime}(r_{0})/2\big{|} and $r_{0}$ is given as a solution of (59). Thus in the large $N$ limit one has again (28), i.e.,

[TABLE]

but where now

[TABLE]

which using a reasoning similar to the one in the previous sections yields

[TABLE]

where now to leading order we expand around $\xi=0$,

[TABLE]

In fact, the argument is now even simpler since $r_{0}(\xi)$ stays finite and does not go to zero as $\xi$ goes to zero – recall that $V^{\prime}(r_{0}(\xi))$ goes to zero when $\xi$ goes to zero and thus (59) does not require that $r_{0}(\xi)$ becomes zero. This is assured by the requirement that $a_{-}>0$.

Now we scale $y$ around the inner edge of the support, $a_{-}$, in an analogous way as for the outer edge

[TABLE]

Note that as $N$ becomes large $y$ approaches $a_{-}$. Moreover, $r_{0}(\xi)$ with $\xi=0$ also approaches $a_{-}$ as $N$ becomes large, as is clear when comparing (10) with (59). This implies that the leading order contribution of (81) now comes from $\xi=0$. This suggests the following change of variables in (81)

[TABLE]

Now expanding $r_{0}(\xi)$ around $\xi=0$ gives

[TABLE]

where recalling (59) and $V^{\prime}(a_{-})=0$ one finds that $\delta_{-}$ is given by

[TABLE]

We notice that since $a_{-}$ is finite and the potential $V$ is convex, we have $\delta_{-}>0$. We therefore obtain

[TABLE]

which by choosing

[TABLE]

becomes

[TABLE]

We thus have

[TABLE]

Similar to the previous section, the pre-factor is one to leading order, if we choose,

[TABLE]

Following the reasoning from the previous section we thus get

[TABLE]

where $\alpha_{N}^{(-)}$ is given in (91) and $F_{-}$ in (82). This shows universality by proving convergence of the rescaled distribution of the eigenvalue with the smallest modulus to a Gumbel distribution. The derivation assumes that the potential fulfills the condition (61) or (62) and that the inner radius $a_{-}>0$, i.e., the support of the eigenvalue density has topology of a ring.

8 Finite $N$ corrections for outer and inner edge of an arbitrary potential

To obtain all logarithmic corrections for the extreme value statistics at the outer edge of an arbitrary potential we start with the generalisation of (50),

[TABLE]

We can now change the integration to $r_{0}$, where the change in the integration measure can be obtained from (59),

[TABLE]

with

[TABLE]

thus eliminating any $\xi$-dependence. Next we change the integration variable to

[TABLE]

Knowing that for our scaling of $y$ the leading contribution from the error function comes from $r_{0}=a_{+}$ we Taylor expand the measure around this value, leading to

[TABLE]

where in addition we substituted $r_{0}(1)=a_{+}+...$ and $r_{0}(0)=a_{-}+...$. Performing the integration, inserting the scaling for $y$, and following the step of Section 5, we arrive at

[TABLE]

with

[TABLE]

Here we defined, as previously,

[TABLE]

but with $\alpha_{N}^{(+)}$ given by (72) and $F_{+}$ is defined in (63). The above expression provides an infinite series of all logarithmic correction terms to the extreme value statistics of the largest eigenvalue scaled around the outer edge of the eigenvalue support of an arbitrary radially symmetric potential. We see that only the first two derivatives of the potential, $V^{\prime}(r)$ and $V^{\prime\prime}(r)$, at the outer edge enter the expression, but not higher order derivatives.

As for the Gaussian case, we can also give a series expression for (99), highlighting that it is an infinite series of slowly varying terms,

[TABLE]

where the only change is in the form of the slowly varying function.

As an example we look at the potential $V(r)=r^{3}/3$. For this potential, $a_{+}=2^{1/3}$ and $F_{+}=3a_{+}$. In Figure 5 we plot for this potential $\exp(\phi_{N}^{(+)}(Y))$, with $\phi_{N}^{(+)}(Y)$ given in (99), against the numerically evaluated $F_{N}(Y)$ for $N=100$.

The derivation for the inner edge proceeds analogously. Firstly, starting with (79) and (80) and proceeding as in the previous case, we obtain for $F_{N}(y):=\mathbb{P}_{N}(|z_{min}|\geq y)$,

[TABLE]

The derivation now proceeds similarly as for the outer edge, where now, due to the different scaling of $y$, the main contribution of the error function comes from $\xi=0$ which corresponds to $r_{0}=a_{-}$. Proceeding as above one obtains,

[TABLE]

with

[TABLE]

as well as

[TABLE]

with $\alpha_{N}^{(-)}$ for given by (91) and $F_{-}$ is defined in (82).

The infinite expansion of $\phi_{N}^{(-)}(Y)$ in terms of logarithmic corrections is exactly as in (101) with all plus sub/super-scripts replaced by minus sub/super-scripts.

As an illustration for how well $\exp(-\phi_{N}^{(-)}(Y))$ approximates the finite $N$ expression $F_{N}(Y)$ we plot in Figure 6 both for the potential $V(r)=r^{2}/2-r$ for $N=500$. For this potential the inner edge is located at $a_{-}=1$ and $F_{-}=1$.

9 Discussion

The celebrated Tracy-Widom distribution provides the extreme value statistics for Hermitian random matrices. For non-Hermitian matrices the situation is somewhat easier as eigenvalues are only weakly correlated and the extreme value statistics is given by the much simpler Gumbel distribution. However, due to the lack of symmetries, non-Hermitian random matrices are in general harder to deal with. A special case are normal random matrices which are non-Hermitian and at the same time allow for a Coulomb gas formulation for general potential.

In this work we investigate the extreme value statistics of normal random matrices and 2D Coulomb gases for general radially symmetric potentials. This is done by extending the orthogonal polynomial approach, introduced in [23] for Hermitian matrices, to normal random matrices and 2D Coulomb gases. We first analyse the simplest case of Gaussian normal random matrices and show convergence of the eigenvalue with the largest modulus rescaled around the outer edge of the eigenvalue support to a Gumbel distribution. One strength of this approach lies in the fact that it immediately generalises to an arbitrary potential $V=V(r)$ with radial symmetry which satisfies the condition (8) together with the condition (61) or (62). We use this to show universality of the distribution of the eigenvalue with largest modulus when rescaled around the outer edge of the eigenvalue support. This provides an alternative, simplified derivation of results presented in [15, 16]. In addition, it is shown that the approach presented here also generalised to compute convergence of the distribution of the eigenvalue with smallest modulus rescaled around the inner edge of the eigenvalue support with topology of an annulus.

As noted in the introduction, it is interesting that one can extend the orthogonal polynomial approach developed for Hermitian random matrices to non-Hermitian random matrices. Furthermore, it is surprising that one can use the same derivational steps to obtain the extreme value statistics of normal random matrices as are used to obtain the Tracy-Widom distribution for Hermitian random matrices [23, 29]. This universal underlying structure suggests that the here presented approach can easily be generalised to new cases, as we indeed show when calculating the extreme values statistics for the eigenvalue with the smallest modulus for the ring distribution.

The here presented approach can also be used to obtain finite $N$ corrections. In (55) we give an analytical expression for the leading order term and all slowly varying finite $N$ corrections for the Gaussian potential. It is quite interesting that one can obtain such a closed-form expression given that there is an infinite number of such slowly varying correction terms as one can see from (57). Such corrections are important in practice given the extremely slow convergence of the extreme value statistics $F_{N}(Y)$, as illustrated in Figure 2 which can be contrasted with the extremely good approximation given by (55) and illustrated in Figure 4 for much smaller $N$. In (99) and (104) we obtain analogous results for the outer and inner edge of an arbitrary radially symmetric potential. The logarithmic corrections only depend on the first and second moment of the potential evaluated at the corresponding edge, but not on higher moments. Figure 5 and 6 show how well the expression including all logarithmic correction terms approximates the finite $N$ expression for two examples of potentials both at the outer as well as at the inner edge.

An interesting direction of future research is to explore the extreme value statistics of the eigenvalue with smallest modulus at the transition when the inner radius of the eigenvalue support goes to zero as $N$ goes to infinity. This transition was analysed in [28] in the context of the mean radial displacement. It would be interesting to see whether the extreme value statistics can be tuned to a different universality class in this transition by using a double scaling limit in which $a_{-}$ goes to zero as $N$ becomes large.

Acknowledgements – The authors would like to thank John Wheater, Max Atkin and Carlos Tomei for fruitful discussions, as well as the anonymous referee for constructive comments. RE is thankful for the financial support of CCPG (PUC-Rio). SZ was previously supported by Nokia Technologies and Lockheed Martin through the Quantum Optimisation and Machine Learning Project at Oxford University.

Appendix A Numerical methods

A.1 Ginibre matrices

For the Gaussian case we can generate samples from the joint eigenvalues density by first creating a complex matrix of random elements

[TABLE]

and then use standard techniques to numerically diagonalise the matrix.

As a concrete example in R one can write the following code, which generates $m$ samples of the matrix $M_{ij}$, diagonalises it, extracts the eigenvalue with the largest modulus and stores it in an array $\vec{y}$.

y <- rep(NULL, m) for (i in 1:m){ M = matrix(complex(real=rnorm(NN),imaginary=rnorm(NN)),ncol=N) y[i] = sort(Mod(eigen(M, symmetric=FALSE, only.values = TRUE)$values)/sqrt(N*2))[N] }

A.2 Monte-Carlo methods for normal random matrices and Coulomb gases

To obtain samples of the eigenvalues $\vec{z}$ from the joint eigenvalue density (4) for a general (radially symmetric) potential we employ Monte-Carlo techniques. To do so we decompose $\vec{z}=\vec{x}+i\vec{y}$ and initialise the $N$ elements of $\vec{x}$ and $\vec{y}$ with uniform random numbers in the interval $(-1,1)$. Then we iterate over the following steps. First one perturbs a given randomly chosen element of $\vec{x}$ and $\vec{y}$ each by a Gaussian random number, multiplied by a scaling $\eta/N$ where $\eta$ is a number of $\mathcal{O}(1)$. Given the perturbed vector $\vec{z}\mkern 2.0mu\vphantom{z}^{\prime}$ we then evaluate the difference in energy in the Boltzmann weight (5),

[TABLE]

If the new configuration has lower energy, i.e., $\Delta E>0$, we change $\vec{z}$ to $\vec{z}\mkern 2.0mu\vphantom{z}^{\prime}$. If the new configuration has larger energy we still allow to change $\vec{z}$ to $\vec{z}\mkern 2.0mu\vphantom{z}^{\prime}$ if $\exp(\Delta E)$ is larger than a random uniform number $u$ in $(0,1)$. This is the standard Metropolis step. The above iteration is repeated $m$ times where $m$ is usually chosen $\mathcal{O}(N^{3})$.

For completeness we provide the pseudo-code of the algorithm below.

1:procedure Monte Carlo simulation for eigenvalue distribution

2:     initialise: $N,m,\eta$ and the function $V(r)$

3:     initialise: $\vec{x}\leftarrow\mathrm{uniform}(-1,1)$

4:     initialise: $\vec{y}\leftarrow\mathrm{uniform}(-1,1)$

5:     for interation in $1,..,m$ do

6:         // perturb one randomly chosen element of $\vec{x}$ and $\vec{y}$:

7:         $k\leftarrow\mathrm{choose}([1,...,N])$

8:         $x^{\prime}[k]\leftarrow x[k]+\eta\delta_{x}/N$ with $\delta_{x}\sim\mathrm{normal}(0,1)$

9:         $y^{\prime}[k]\leftarrow y[k]+\eta\delta_{y}/N$ with $\delta_{y}\sim\mathrm{normal}(0,1)$

10:         // Calculate energy difference for $\vec{z}=\vec{x}+i\vec{y}$ and $\vec{z}\mkern 2.0mu\vphantom{z}^{\prime}=\vec{x}\mkern 2.0mu\vphantom{x}^{\prime}+i\vec{y}\mkern 2.0mu\vphantom{y}^{\prime}$:

11:         $\Delta E=N^{2}(S_{\mathrm{eff}}(\vec{z})-S_{\mathrm{eff}}(\vec{z}\mkern 2.0mu\vphantom{z}^{\prime}))$

12:         // sample a uniform number and do the Metropolis step:

13:         $u\sim\mathrm{uniform}(0,1)$

14:         if $\exp(\Delta E)>u$ then

15:              $\vec{x}\leftarrow\vec{x}\mkern 2.0mu\vphantom{x}^{\prime}$

16:              $\vec{y}\leftarrow\vec{y}\mkern 2.0mu\vphantom{y}^{\prime}$

17:         end if

18:     end for

19:end procedure