Extreme Eigenvalues and the Emerging Outlier in Rank-One Non-Hermitian Deformations of the Gaussian Unitary Ensemble

Abstract

Complex eigenvalues of random matrices $J=\mathrm{GUE}+i\gamma \mathrm{diag}(1,0,\dots ,0)$ provide the simplest model for studying resonances in wave scattering from a quantum chaotic system via a single open channel. It is known that in the limit of large matrix dimensions $N\gg 1$ the eigenvalue density of J undergoes an abrupt restructuring at $\gamma =1$, the critical threshold beyond which a single eigenvalue outlier (“broad resonance”) appears. We provide a detailed description of this restructuring transition, including the scaling with N of the width of the critical region about the outlier threshold $\gamma =1$ and the associated scaling for the real parts (“resonance positions”) and imaginary parts (“resonance widths”) of the eigenvalues which are farthest away from the real axis. In the critical regime we determine the density of such extreme eigenvalues, and show how the outlier gradually separates itself from the rest of the extreme eigenvalues. Finally, we describe the fluctuations in the height of the eigenvalue outlier for large but finite N in terms of the associated large deviation function.

1. Introduction

Rank-one non-normal deformations of the Gaussian and Circular Unitary Ensembles are a useful analytic tool for studying statistics of resonances in quantum scattering from a chaotic domain via a single channel [1,2]. As surveyed in [2,3], these random matrix ensembles are integrable in the sense that the joint probability density of their complex eigenvalues and, in some spectral scaling limits of interest, the eigenvalue correlation functions can be determined in a closed form. Such integrability, which also proves to be useful in other physics contexts, see, e.g., [4], extends to a certain degree to the deformed $\beta –$Gaussian and $\beta –$circular ensembles [5,6], especially to the classical values $\beta =1,4$ [7,8], but is lost if the underlying normal random matrix ensemble (Hermitian or unitary) is not integrable, as is the case with, e.g., finite rank non-Hermitian deformations of Wigner matrices [9,10,11] or band matrices [12]. Still, the latter matrices are found to share, in appropriate parameter ranges, some statistical characteristics of their complex eigenvalues and eigenvectors with their integrable counterparts.

In this paper, we aim to investigate complex eigenvalues with extreme imaginary parts for the rank-one non-Hermitian deformations of the Gaussian Unitary Ensemble (GUE) by exploiting the above-mentioned integrability. The latter feature gives access to the asymptotics of the eigenvalue density in the complex plane on mesoscopic scales and allows us to carry out a quantitative analysis of the separation of the eigenvalue outlier (which is known to exist in this model [9,10]) from the rest of the eigenvalues. Eigenvalue outliers in the complex plane have recently attracted renewed interest [11,13,14,15]. Our analysis refines and complements the existing knowledge about the outliers of nearly Hermitian matrices [9,10,11] albeit for arguably the simplest model of its type. As we will demonstrate, despite the simplicity of the model, its extreme eigenvalues exhibit an interesting transition at a certain value of the deformation parameter, with rich critical behaviour which deserves to be studied in more detail.

The non-Hermitian matrices that we consider are of the form

$$\begin{array}{c}J=H+i\Gamma ,\hfill \end{array}$$

(1)

where H is a GUE matrix and $\Gamma$ is a diagonal matrix with all diagonal entries being zero except the first one,

$$\begin{array}{c}\Gamma =\gamma \mathrm{diag}(1,0,\dots ,0).\hfill \end{array}$$

(2)

Denoting the matrix dimension by N, we fix the global spectral scale by the condition that the expected value of $\mathrm{Tr}{H}^2$ is N. Then the joint probability density function (JDPF) of matrix elements of the GUE matrix H is

$$\begin{array}{c}{f}_N\left(H\right)=const\times \exp \left\{-\frac{N}{2}\mathrm{Tr}{H}^2\right\}.\hfill \end{array}$$

(3)

With this normalisation, the limiting eigenvalue distribution of H, as the matrix dimension is approaching infinity, is supported on the interval $[-2,2]$, and, inside this interval, the eigenvalue density is $\nu \left(X\right)=\frac{1}{2\pi }\sqrt{4-X^2}$.

Note that due to the invariance of the JPDF (3) with respect to unitary rotations $H\to UH{U}^{-1}$ one may equivalently replace $\Gamma$ in (2) with any other rank-one Hermitian matrix. Without loss of generality we may also assume $\gamma$ to be positive. Then the eigenvalues $X_j+i{Y}_j$ of matrices J (1)–(3) are all in the upper half of the complex plane and for N large they all, except possibly one outlier, lie just above the interval $[-2,2]$ of the real line. Whether such an outlier is present or not is determined by the value of $\gamma$. For fixed values of $\gamma <1$, almost surely, for N sufficiently large, all N eigenvalues lie within distance $c_N{N}^{-1}$ from the real line, with $c_N=o\left(N^{ϵ}\right)$ for every $ϵ>0$ [9]. Furthermore, if $\gamma >1$ then the same is true of all but one eigenvalue. This outlier lies much higher in the complex plane: to leading order in N, its imaginary part (the “height”) is $\gamma -{\gamma }^{-1}$ [9,10,14]. For precise statements and proofs we refer the reader to [9,10] where these and similar facts were established for finite rank non-Hermitian deformations of real symmetric matrices with independent matrix entries.

For finite but large matrix dimensions, one would expect to find a transition region of infinitesimal width $\Omega$ about the outlier threshold value $\gamma =1$ which captures the emergence of the outlier from the sea of low lying eigenvalues. Questions about the scaling of $\Omega$ with N and the corresponding characteristic height and distribution of the eigenvalues that lie farthest away from the real line are natural and interesting in this context. These are open questions in the mathematics and mathematical physics literature on the subject.

Apart from the mathematical curiosity, there is also motivation coming from physics. In the physics literature, the eigenvalues of J are associated with the zeroes of a scattering matrix in the complex energy plane, and their complex conjugates with the poles of the same scattering matrix, known as “resonances”. The latter are obviously the eigenvalues of matrices (1)–(2) with $\gamma$ replaced by $-\gamma$. In that context the absolute value of the eigenvalue’s imaginary part is associated with the “resonance width”. The eigenvalues close to real axis are called “narrow resonances” and the outlier is called the “broad resonance”. The use of the Gaussian Unitary Ensemble for H is justified by invoking the so-called Bohigas-Giannoni-Schmidt conjecture [16] describing spectral statistics of highly excited energy levels of some classes of systems whose classical counterparts are chaotic. The resulting ensemble J is then an important ingredient in characterising statistical properties of scattering matrices in systems with quantum chaos and no time-reversal invariance, see [1] for description of the associated framework going back to the pioneering paper [17]. In that framework, the phenomenon of the outlier separation and the simultaneous movement of the rest of the eigenvalues towards the real axis was first discussed, albeit at a heuristic level, already in early theoretical works [18,19], the latter work even establishing the correct asymptotic position of the outlier. Later on, this phenomenon got considerable attention under the name “resonance trapping” and eventually was observed in experiments [20].

Very recently, Dubach and Erdős [11] performed a detailed analysis of the eigenvalue trajectories, with respect to changing the parameter $\gamma$, in the random matrix ensemble $H+i\gamma vs.v^{\ast }$ in the settings when H is assumed to be a Wigner matrix and v a column vector of unit length. It turned out that the evolution of the eigenvalues is governed by a system of deterministic first-order differential equations subject to random initial conditions, with the initial positions and velocities expressed in terms of the eigenvalues and eigenvectors of H. In addition, under suitable conditions on the distribution of matrix entries of H ensuring the validity of the uniform isotropic local law (Theorem 5 in [11]), Dubach and Erdős proved that with high probability the eigenvalue outlier is distinctly separated from the rest of the eigenvalues for all

$$\begin{array}{c}\gamma >1+\frac{N^{\epsilon }}{\sqrt[3]{N}},\epsilon >0.\hfill \end{array}$$

(4)

Moreover, if $\epsilon <1/3$, i.e., if $N^{-1/3+\epsilon }$ is asymptotically small, the outlier’s height is $2N^{-1/3+\epsilon }$ and its real part is in the window of width $N^{-1/3-ϵ/4}$ around the origin, whereas all other eigenvalues are no higher than $N^{-1/3-\epsilon }$. In addition, with high probability, for all

$$\begin{array}{c}\gamma <1-\frac{N^{\epsilon }}{\sqrt[3]{N}},\epsilon >0,\hfill \end{array}$$

(5)

no eigenvalue reaches the heights

$$\begin{array}{c}Y=\frac{m}{\sqrt[3]{N}},m>0.\hfill \end{array}$$

(6)

These findings suggest that the width $\Omega$ of the transition region around $\gamma =1$ scales with as $N^{-1/3}$ for N large. Naturally, for $\gamma$ inside this region one would expect to find several eigenvalues, including the emerging “atypical” outlier, with imaginary parts on the critical scale (6) much exceeding the height $O\left(N^{-1}\right)$ of low lying eigenvalues, as illustrated in Figure 1. One might call such eigenvalues “typical extremes” to emphasise atypicality of the emerging outlier.

To a large extent our paper is motivated by [11] and aims to provide quantitative insights into this picture of the outlier emerging from the cloud of extreme eigenvalues. Whilst the approach of Dubach and Erdős is dynamical (fix matrix H and study eigenvalue trajectories as the magnitude $\gamma$ of the deformation increases), our approach is statistical (fix a scale for $\gamma$ and count the number of eigenvalues on characterisitc spectral scales in the complex plane averaged over the distribution of H which, for technical reasons, we assume to be GUE). Our present approach is limited to the expected values; analysing higher order moments is left as an interesting problem for future investigations. However, even with such a basic tool we are able to develop rather detailed quantitative understanding of the outlier separation and the associated restructuring transition in the spectra of matrices J.

As such, the two approaches complement each other very well. For example, we prove that for

$$\begin{array}{c}\gamma =1+\frac{\alpha }{\sqrt[3]{N}},\alpha \in \mathbb{R},\hfill \end{array}$$

(7)

the expected number of eigenvalues whose height exceeds the level (6) is asymptotically given by the integral ${\int }_m^{\infty }{p}_{\alpha }^{\left(\mathrm{Im}\right)}\left(m^{\prime }\right)d{m}^{\prime }$ with density

$$\begin{array}{c}{\tilde{p}}_{\alpha }^{\left(\mathrm{Im}\right)}\left(m\right)=\frac{1}{2\sqrt{\pi }}\frac{\frac{3}{2m}+{\left(\frac{3m}{2}-\alpha \right)}^2}{m^{3/2}}{e}^{-m{\left(\alpha -\frac{m}{2}\right)}^2},m>0.\hfill \end{array}$$

This density is the average density of the extreme eigenvalues at height (6). Together with findings in [11] this result establishes that the width $\Omega$ of the transition region around $\gamma =1$ indeed scales with $N^{-1/3}$. Similarly, we are able to determine the average density of extreme eigenvalues $Z_j=X_j+i{Y}_j$ of J near the origin in the complex plane in the critical scaling regime when when $q+im=\sqrt[3]{N}Z=O\left(1\right)$. As a function of coordinates q and m, this density, when appropriately rescaled, is given by

$$\begin{array}{c}{\tilde{p}}_{\alpha }(q,m)=\frac{1}{4\pi m}\left[\frac{1}{m}+\frac{q^2}{4}+{\left(\frac{3m}{2}-\alpha \right)}^2\right]e^{-m\left[\frac{q^2}{4}+{\left(\alpha -\frac{m}{2}\right)}^2\right]},q\in \mathbb{R},m>0.\hfill \end{array}$$

It can be verified that ${\int }_{-\infty }^{+\infty }{\tilde{p}}_{\alpha }(q,m)dq={\tilde{p}}_{\alpha }^{\left(\mathrm{Im}\right)}\left(m\right)$, implying that the population of extreme eigenvalues at the critical height (6) which generates the eventual outlier (as $\alpha$ is approaching infinity) is constrained to a narrow vertical strip of width $O\left(N^{-1/3}\right)$ about the origin (the centre of the eigenvalue band of H). Thus, our results both confirm and complement the analysis in [11], and show that it indeed touched the optimal scales in $\gamma$ (7), both along the real and imaginary axes.

We would like to conclude this section with a short description of the structure of our paper. In Section 2 we develop quantitative heuristic analysis of the outlier separation. This analysis elucidates the emerging critical scaling in $\gamma$ and the critical spectral scalings in the complex plane and provides a useful background for rigorous calculations later on. This section also offers our outlook on the universality of the exponent $-1/3$ in (7). Section 3 contains the statement of our main results and discussion. In Section 4 we express the expected density of eigenvalues of J and the density of their imaginary parts at finite matrix dimensions in terms of, respectively, Hermite and Laguerre polynomials. These expressions are then used in Section 5 and Section 6 for asymptotic analysis of eigenvalue densities in various scaling limits. The two appendices contain derivations of technical auxiliary results.

2. Low Lying Eigenvalues and Their Extremes: A Heuristic Outlook

Before presenting our main results in the next Section, we would like to offer our quantitative heuristic insights into the outlier separation elucidating the emerging scalings and mechanisms behind them and providing a useful background for rigorous calculations later on.

With $z_j=X_j+i{Y}_j$ standing for the eigenvalues of matrices $J=H+i\Gamma$, the angular brackets $〈\dots 〉$ standing for averaging over the GUE matrix H (3), and $\delta \left(X\right)$ for the Dirac delta-function, the expected number of eigenvalues of J in domain D can be computed by integrating the mean eigenvalue density

$$\begin{array}{c}{\rho }_N(X,Y)=〈\frac{1}{N}{\displaystyle \sum _{j=1}^N}\delta (X-X_j)\delta (Y-Y_j)〉\hfill \end{array}$$

(8)

over D and multiplying the result by N. For example, the expected number ${\mathcal{N}}_{\gamma }\left(Y\right)$ of the eigenvalues of J which lie above the line $\mathrm{Im}z=Y$ in the complex plane is given by the integral

$$\begin{array}{c}{\mathcal{N}}_{\gamma }\left(Y\right)=N{\int }_{-\infty }^{\infty }{\int }_Y^{\infty }{\rho }_N(X,Y^{\prime })dXd{Y}^{\prime }=N{\int }_Y^{\infty }{\rho }_N^{\left(\mathrm{Im}\right)}\left(Y^{\prime }\right)d{Y}^{\prime },\hfill \end{array}$$

(9)

where ${\rho }_N^{\left(\mathrm{Im}\right)}\left(Y\right)$ is the mean density of the imaginary parts irrespective of the value of the real part,

$$\begin{array}{c}{\rho }_N^{\left(\mathrm{Im}\right)}\left(Y\right)=〈\frac{1}{N}{\displaystyle \sum _{j=1}^N}\delta (Y-Y_j)〉.\hfill \end{array}$$

(10)

Guided by the eigenvalue perturbation theory one can expect that the typical height Y of the eigenvalues whose real part is close to a point $X\in (-2,2)$ in the spectral bulk scales with the mean separation $\Delta ={\left(N\nu \left(X\right)\right)}^{-1}$ between neighbouring real eigenvalues of the GUE matrix H in the limit $N\to \infty$. On a more formal level, introducing the scaled version of ${\rho }_N(X,Y)$ [21]

$$\begin{array}{c}{\tilde{\rho }}_N(X,y):=\frac{1}{\nu \left(X\right)}〈\frac{1}{N}{\displaystyle \sum _{j=1}^N}\delta \left(X-X_j\right)\delta \left(y-2\pi \nu \left(X\right)N{Y}_j\right)〉,-2<X<2,\hfill \end{array}$$

(11)

one finds that such scaled density is well-defined in the limit of large matrix dimensions [1,2,21,22]: for every $y>0$

$$\begin{array}{c}{\displaystyle \tilde{\rho }(X,y):=\underset{N\to \infty }{\lim }{\tilde{\rho }}_N(X,y)=-\frac{d}{dy}\left[e^{-yg\left(X\right)}\frac{\sinh y}{y}\right],g\left(X\right)=\frac{\gamma +\frac{1}{\gamma }}{2\pi \nu \left(X\right)},}\hfill \end{array}$$

(12)

confirming that locally the typical height of low lying eigenvalues scales with $\Delta ={\left(N\nu \left(X\right)\right)}^{-1}$.

Globally, the typical height of low lying eigenvalues scales with $N^{-1}$. Intuitively, this can be understood from the exact sum rule

$$\begin{array}{c}{\displaystyle \sum _{j=1}^N}{Y}_j=\gamma \hfill \end{array}$$

(13)

which follows from the obvious relation $\mathrm{Tr}J=i\gamma +\mathrm{Tr}H$. On a more formal level, consider the expected fraction of the eigenvalues of J which lie above the level $\mathrm{Im}z=Y$, and set $y=NY$. In the limit $N\to \infty$,

$$\begin{array}{c}\frac{1}{N}{\mathcal{N}}_{\gamma }\left(\frac{y}{N}\right)\sim {\int }_{-2}^{2}dX\nu \left(X\right){\int }_{2\pi \nu \left(X\right)y}^{\infty }d{y}^{\prime }\tilde{\rho }(X,y^{\prime })\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & =\frac{e^{-y\left(\gamma +\frac{1}{\gamma }\right)}}{y}{\int }_{-2}^2\frac{dX}{4\pi }\left(e^{y\sqrt{4-X^2}}-e^{-y\sqrt{4-X^2}}\right).\hfill \end{array}$$

(15)

The integral in (15) is the modified Bessel function $I_1\left(2y\right)$. Therefore,

$$\begin{array}{c}\underset{N\to \infty }{\lim }\frac{1}{N}{\mathcal{N}}_{\gamma }\left(\frac{y}{N}\right)=\frac{e^{-y\left(\gamma +\frac{1}{\gamma }\right)}}{y}{I}_1\left(2y\right).\hfill \end{array}$$

(16)

From this,

$$\begin{array}{c}{\tilde{\rho }}^{\left(\mathrm{Im}\right)}\left(y\right):=\underset{N\to \infty }{\lim }\frac{1}{N}{\rho }_N^{\left(\mathrm{Im}z\right)}\left(\frac{y}{N}\right)=-\frac{d}{dy}\left[\frac{e^{-y\left(\gamma +\frac{1}{\gamma }\right)}}{y}{I}_1\left(2y\right)\right]\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill & =\frac{e^{-y\left(\gamma +\frac{1}{\gamma }\right)}}{y}\left[\left(\gamma +\frac{1}{\gamma }-2\right)I_1\left(2y\right)-I_0\left(2y\right)-I_2\left(2y\right)\right].\hfill \end{array}$$

(18)

The density ${\tilde{\rho }}^{\left(\mathrm{Im}\right)}\left(y\right)$ is the mean density of the scaled imaginary parts $y_j=N{Y}_j$ in the limit of large matrix dimensions. Even though it describes low lying eigenvalues it contains some useful information about eigenvalues higher up in the complex plane.

As an example, consider the expected value of the sum of the imaginary parts of low lying eigenvalues. Using definition (10), the sum rule (13) implies that

$$\begin{array}{c}N{\int }_0^{\gamma }Y{\rho }_N^{\left(\mathrm{Im}\right)}\left(Y\right)dY=\gamma .\hfill \end{array}$$

(19)

Upon rescaling $y=NY$, one could naively jump to the conclusion that ${\int }_0^{\infty }y{\tilde{\rho }}^{\left(\mathrm{Im}\right)}\left(y\right)dy=\gamma$. However, by making use of (17) and integral 6.623(3) in [23], one actually finds that

$$\begin{array}{c}{\int }_0^{\infty }y{\tilde{\rho }}^{\left(\mathrm{Im}\right)}\left(y\right)dy={\int }_0^{\infty }\frac{e^{-y\left(\gamma +\frac{1}{\gamma }\right)}}{y}{I}_1\left(2y\right)dy=\frac{\gamma +\frac{1}{\gamma }-\sqrt{{\left(\gamma +\frac{1}{\gamma }\right)}^2-4}}{2}=\left\{\begin{array}{cc}\gamma ,\hfill & \mathrm{if}\gamma <1,\hfill \\ \frac{1}{\gamma },\hfill & \mathrm{if}\gamma >1.\hfill \end{array}\right.\hfill \end{array}$$

Thus, if $\gamma <1$ then the imaginary parts of low lying eigenvalues indeed add up to $\gamma$, in full agreement with the sum rule (19), whereas if $\gamma >1$ they add up only to $\frac{1}{\gamma }<\gamma$. The sum rule deficit $\gamma -\frac{1}{\gamma }$ is exactly the imaginary part of the outlier, and suggests that the rescaled limiting density of low lying eigenvalues, ${\tilde{\rho }}^{\left(\mathrm{Im}\right)}\left(y\right)$, precisely misses the delta-functional mass $\frac{1}{N}\delta \left(y-\left(\gamma -\frac{1}{\gamma }\right)\right)$.

As another example, consider the asymptotic form of ${\tilde{\rho }}^{\left(\mathrm{Im}\right)}\left(y\right)$ when $y\gg 1$. It is markedly different depending on whether $\gamma =1$ or not. In the later case, using in (17) the asymptotic expansion for the modified Bessel function of large argument, $I_p\left(x\right)\sim \frac{e^x}{\sqrt{2\pi x}}(1-\frac{4p^2-1}{8x}+\dots )$ one finds an exponential decay, whilst in the former case the decay is algebraic:

$$\begin{array}{c}{\tilde{\rho }}^{\left(\mathrm{Im}\right)}\left(y\right)=\left\{\begin{array}{cc}{\displaystyle \frac{e^{-y\frac{{(1-\gamma )}^2}{\gamma }}}{2\sqrt{\pi }{y}^{3/2}}\left[\frac{{(1-\gamma )}^2}{\gamma }+\frac{30-3(\gamma +{\gamma }^{-1})}{16y}+O\left(\frac{1}{y^2}\right)\right]}\hfill & \mathrm{if}\gamma \ne 1,\hfill \\ {\displaystyle \frac{3}{4\sqrt{\pi }}\frac{1}{y^{5/2}}+O\left(\frac{1}{y^{7/2}}\right)}\hfill & \mathrm{if}\gamma =1.\hfill \end{array}\right.\hfill \end{array}$$

(20)

It is instructive to return to the unscaled imaginary part Y and take a closer look at the expected number of the eigenvalues of J exceeding the level $\mathrm{Im}z=Y$ in the limit $N\to \infty$. It is evident from (16) that

$$\begin{array}{c}{\mathcal{N}}_{\gamma }\left(Y\right)\sim \frac{e^{-NY\left(\gamma +\frac{1}{\gamma }\right)}}{Y}{I}_1\left(2NY\right),\hfill \end{array}$$

(21)

provided $NY=O\left(1\right)$. Extending this asymptotic relation to large values of $NY$ allows one to get insights, even if only heuristically, about the characteristic scale of the highest placed among the low lying eigenvalues. Along these lines, we define the characteristic scale of the height of typical extreme eigenvalues as such level $Y_e$ that the expected number of eigenvalues with imaginary part exceeding $Y_e$ is of order of unity:

$$\begin{array}{c}{\mathcal{N}}_{\gamma }\left(Y_e\right)=O\left(1\right).\hfill \end{array}$$

(22)

We add the word typical to exclude the atypical eigenvalue (the outlier) which is known to exist when $\gamma >1$. Now, assuming $N{Y}_e$ to be large (but still anticipating $Y_e\ll 1$) one can replace the Bessel function in (21) by its corresponding asymptotic expression and approximate:

$$\begin{array}{c}{\mathcal{N}}_{\gamma }\left(Y_e\right)\approx \frac{e^{-N{Y}_e\frac{{(1-\gamma )}^2}{\gamma }}}{2\sqrt{\pi N}{Y}_e^{3/2}},1\ll N{Y}_e\ll N.\hfill \end{array}$$

(23)

The condition in (22) then leads to two essentially different scenarios depending on the value of $\gamma$. Namely, for every fixed positive $\gamma \ne 1$ the characteristic scale of the typical extreme values is, to leading order in N, $O\left(N^{-1}\ln N\right)$. On the other hand, if $\gamma =1$ then the typical extreme values raise from the sea of low lying eigenvalues to a much higher height of $O\left(N^{-1/3}\right)$. This change of scale for extreme values is easy to trace back to the emerging power-law decay in the vicinity of $\gamma =1$ which is evident in (20).

In fact, as evident from (23), the typical extreme values scale as $Y_e=O\left(N^{-1/3}\right)$ not only at $\gamma =1$, but also as long as $|1-\gamma |\propto N^{-1/3}$. Actually, by setting simultaneously $\gamma =1+\alpha N^{-1/3}$ and $Y_e=m{N}^{-1/3}$ the asymptotic relation (23) is converted into

$$\begin{array}{c}{\mathcal{N}}_{1+\frac{\alpha }{\sqrt[3]{N}}}\left(\frac{m}{\sqrt[3]{N}}\right)\approx \frac{e^{-m{\alpha }^2}}{2\sqrt{\pi }{m}^{3/2}},\hfill \end{array}$$

(24)

an expression that is indeed of order of unity for all fixed values of $\alpha$ and $m>0$. Thus, the width of the transition region about $\gamma =1$ must scale with $N^{-1/3}$. Combined with the existence of a distinct outlier at height $\gamma -{\gamma }^{-1}\gg Y_e$ one may indeed see that our heuristic argument perfectly agrees with the conjecture of Dubach and Erdős about the critical scaling $\gamma =1+O\left(N^{-1/3}\right)$ where the separation of typical and atypical extreme values happens.

Before continuing our exposition of the heuristics behind the restructuring of the density of complex eigenvalues we would like to make two remarks.

Remark 1.

To make further contact with the standard subject of extreme value statistics, it is useful to recourse to the classical theory of extreme values for i.i.d. sequences of random variables$y_1,\dots ,y_N$, a succinct albeit informal summary of which can be found in, e.g., [24]. In that case the probability law of extreme values is characterised by the tail behaviour of the “parent” probability density function (pdf)$p\left(y\right)$ of $y_j$ and is essentially universal in the limit $N\to \infty$. In our context, the pertinent case for comparison is that of non-negative continuous i.i.d. random variables with the parent distribution supported on the entire semi-axis $[0,\infty )$. Then only two possibilities may arise. Those sequences which are characterised by the power-law decaying pdf $p\left(y\right)\sim A{y}^{-(1+\alpha )}$, $\alpha >0$, as $y\to \infty$ have their extreme values scaling with ${(AN/\alpha )}^{1/\alpha }$ and the distribution of their maximum, $y_{max}=\max (y_1,\dots ,y_N)$, after rescaling converges to the so-called Fréchet law in the limit $N\to \infty$. In contrast, if the parent pdf decays faster than any power, e.g., if $\ln p\left(y\right)\sim -y^{\delta }$, $\delta >0$, then, to leading order, extreme values scale with ${\left(\ln N\right)}^{1/\delta }$, and the distribution of the largest value $y_{max}$, converges, after a shift and further rescaling, to the so-called Gumbel law. Although, the imaginary parts of complex eigenvalues in the random matrix ensemble (1)–(3) are not at all independent (as is evident from their JPDF (51) resulting in a non-trivial determinantal two-point and higher order correlation functions at the scale $N^{-1}$, see [22]), our scaling predictions for the typical extreme eigenvalues are in formal correspondence with the i.i.d. picture: a Gumbel-like scaling (with $\delta =1$) if $\gamma \ne 1$ and a Fréchet-like scaling (with $\alpha =3/2$) if $\gamma =1$. This is exactly as would have been implied in the i.i.d. picture by the tail behaviour of the mean eigenvalue densities in the two cases in (20). This fact naturally suggests to conjecture Gumbel statistics for the typical largest imaginary part (excluding possible outlier) for any $\gamma \ne 1$, changing to a Fréchet-like law for $\gamma =1$, with a possible family of $\alpha -$ dependent nontrivial extreme value statistics in the crossover critical regime $\gamma =1+\alpha N^{-1/3}$. Although we are not able to shed light on the distribution of typical extreme eigenvalues in the random matrix ensemble (1)–(3), we will discuss some results in that direction for a somewhat related model at the end of the next section.

Remark 2.

The phenomenon of resonance width restructuring with increasing the coupling to continuum (controlled in the present model by the parameter γ) and the emergence of the broad resonance has many features in common with the so-called super-radiant phenomena in optics. This is well known in the physics literature, see [25] and references therein. Here, we would like to point to a similarity of the spectral restructure in the random matrix ensemble (1)–(3) to a process in a different physics context, the so-called “condensation transition” which occurs in models of mass transport when the globally conserved mass M exceeds a critical value, see, e.g., [26] for a review. In such a regime, the excess mass forms a localised in space condensate coexisting with a background fluid in which the remaining mass is evenly distributed over the rest of the system. A particularly simple case for analysing the condensation phenomenon is when the system has a stationary state such that probability of observing a configuration of masses $m_i$ factorises into the form ${\prod }_{i}f\left(m_i\right)\delta ({\sum }_i{m}_i-M)$. In that context again the tail behaviour of the “parent” mass density $f\left(m\right)$ plays important role. Although we would like to stress again that in our model the imaginary parts of the complex eigenvalues are not independent, the analogy with the condensation phenomenon is quite evident.

Essentially the same heuristic analysis as in the above helps to clarify the numerically observed fact of the outlier emerging mostly close to the origin of the spectrum Re z = 0. From this angle it is instructive to ask what should be the scale of extreme values for eigenvalues satisfying $\left|\mathrm{Re}z\right|<W$, that are sampled in a window of a small widths $W\ll 1$ around the origin (still assuming typically many eigenvalues in the window, so that $W\gg \Delta \sim 1/N$). The total mean number of eigenvalues in the window W whose imaginary parts exceed the level Y (but still formally remain of the order of $1/N$) is now given by

$$\begin{array}{c}{\mathcal{N}}_{\gamma ,W}\left(Y\right)=\frac{e^{-NY\left(\gamma +\frac{1}{\gamma }\right)}}{4\pi Y}\left[T_W\left(NY\right)-T_W(-NY)\right],T_W\left(NY\right)=2{\int }_0^W{e}^{NY\sqrt{4-X^2}}dX.\hfill \end{array}$$

(25)

For $NY\gg 1$ the term $T_W(-NY)$ is exponentially suppressed, while the integral in $T_W\left(NY\right)$ is dominated by $X\ll 1$ and with required accuracy yields the leading-order expression in the form:

$$\begin{array}{c}{\mathcal{N}}_{\gamma ,W}(Y\gg 1/N)\approx \frac{e^{-NY\left(\frac{{(1-\gamma )}^2}{\gamma }\right)}}{2\pi Y^{3/2}}\sqrt{\frac{2}{N}}{\int }_0^{W\sqrt{NY/2}}{e}^{-\frac{t^2}{2}}dt.\hfill \end{array}$$

(26)

Now, let us assume that both the width W of the window and the parameter $\gamma$ scale with N in this non-trivial way as

$$\begin{array}{c}W\sim N^{-1+\kappa },0<\kappa \le 1,\mathrm{and}\gamma =1-\alpha N^{-\delta },0<\delta \le \infty ,\alpha \in \mathbb{R},\hfill \end{array}$$

(27)

and again apply the same heuristic procedure to determine the scale of extreme values $Y_e(\kappa ,\delta )$ in the window as $N\to \infty$ for given values of exponents $\kappa$ and $\delta$. A straightforward computation shows that the arising scale of extreme values very essentially depends on whether the parameter $\delta$ satisfies $0<\delta <1/3$ or $1/3\le \delta <1$. In the former case we find

$$\begin{array}{c}{Y}_e(\kappa ,0<\delta <1/3)\approx \left\{\begin{array}{cc}{\displaystyle N^{-1+\kappa },}\hfill & \mathrm{if}0<\kappa <2\delta ,\hfill \\ {\displaystyle \frac{\kappa -2\delta }{\alpha }{N}^{-1+2\delta }\ln N,}\hfill & \mathrm{if}2\delta <\kappa <1-\delta ,\hfill \\ {\displaystyle \frac{1-3\delta }{\alpha }{N}^{-1+2\delta }\ln N,}\hfill & \mathrm{if}1-\delta <\kappa <1.\hfill \end{array}\right.\hfill \end{array}$$

(28)

whereas in the latter case

$$\begin{array}{c}{Y}_e(\kappa ,1/3\le \delta <1)\approx \left\{\begin{array}{cc}{\displaystyle N^{-1+\kappa },}\hfill & if0<\kappa <2/3,\hfill \\ {\displaystyle N^{-1/3},}\hfill & if2/3<\kappa <1.\hfill \end{array}\right.\hfill \end{array}$$

(29)

One may say that as long as $\delta <1/3$ the system is not fully in the well-developed “critical regime”, and the extreme value scale is growing with the window width, saturating at the Gumbel-like scale $N^{-1+2\delta }\ln N$. At the same time, as long as $\delta$ exceeds the threshold value $\delta =1/3$, the typical extreme values reach the scale $Y_e=O\left(N^{-1/3}\right)$ as long as they are sampled in a window of width exceeding the scale $W_c=O\left(N^{-1/3}\right)$, thus containing $O\left(N^{2/3}\right)$ eigenvalues. This heuristics suggests that only eigenvalues satisfying $\left|X\right|<W_c$ typically have a nonvanishing probability to reach to the maximum height in the complex plane, and eventually to generate an outlier as $\alpha$ increases. It would be also natural to expect the corresponding extreme eigenvalues to follow the Fréchet-type statistics for their imaginary parts, as opposed to the Gumbel statistics in the former case.

We would like to end our heuristic considerations with a brief heuristic outlook on the universality of the scaling factor $N^{-1/3}$ which is key to the correct description of the transition in question. As is evident from (23) the exponent $-1/3$ is implied by the scaling law

$$\begin{array}{c}{\mathcal{N}}_{\gamma =1}\left(Y\right)\propto \frac{1}{N^{1/2}{Y}^{3/2}}\hfill \end{array}$$

(30)

in the limit $NY\gg 1$ for the expected number of eigenvalues exceeding the level line $\mathrm{Im}z=Y$. Thus, to investigate the extent of universality of this exponent one needs to trace the origin of the scaling law (30). This can be readily done by returning to the asymptotic relation (14) and (15) which was used to obtain (30). On evaluating the integral in (15) for large values of $y=Y/N$ by the Laplace method it becomes immediately apparent that the power $Y^{-3/2}$ on the right-hand side in (30) and, hence, the exponent in question stems from the quadratic shape of the limiting GUE eigenvalue density function $\nu \left(X\right)={\left(2\pi \right)}^{-1}\sqrt{4-X^2}$ in the vicinity of its maximum. It is natural to conjecture that had one started from a random Hermitian matrix H taken from the broad class of invariant ensembles characterised by joint probability density function $\propto \exp \left\{-N\mathrm{Tr}V\left(H\right)\right\}$ with a suitable potential $V\left(H\right)$ (or from the class of Wigner matrices with suitable conditions on the iid entries), the asymptotic expression (12) for the scaled eigenvalue density ${\tilde{\rho }}_N(X,y)$ would retain its validity after replacing $\nu \left(X\right)$ in (11) and (12) by the corresponding limiting eigenvalue density of H. For example, as was shown albeit not fully rigorously in [27], such universality of the scaled eigenvalue density near the real line is exhibited by almost Hermitian random matrices which are morally similar to finite rank non-Hermitian deviations as in (1) and (2). Since asymptotic relation (14) and (15) is the immediate corollary of (12), one then concludes that as long as the limiting eigenvalue density of H has a single global parabolic-shaped maximum, an additive rank-one non-Hermitian deformation will demonstrate the same type of critical scaling for its extreme complex eigenvalues, and, most probably, after appropriate rescaling, the same type of critical behaviour of the density of imaginary parts as described in the next section. One can however imagine invariant ensembles where the mean eigenvalue density would have a non-parabolic behaviour close to the maximum point.

From this point of view, the noticed in [11] resemblance of the $N^{-1/3}$ critical scaling in the present model and the edge scaling of extreme real eigenvalues of GUE, which, e.g., manifests itself in the so-called BBP [28] transition under additive rank-one Hermitian perturbation of the GUE, looks to us purely coincidental. Indeed, the latter is known to have its origin in the square root behaviour of the mean density $\nu \left(X\right)$ at the spectral edges where $\nu \left(X\right)$ vanishes, and as such seems to have nothing to do with the behaviour of the same density close to its maximal point.

3. Main Results and Discussion

Our first result concerns the mean density of imaginary parts ${\rho }_N^{\left(\mathrm{Im}\right)}\left(Y\right)$ (10) in the large deviation regime $Y\gg N^{-1}$. We note that no eigenvalue of J has imaginary part equal or greater than $\gamma$. This is a consequence of the sum rule (13). Therefore we only consider the range of values $Y\in [0,\gamma )$.

Theorem 1.

Consider the random matrix ensemble (1)–(3) in the scaling regime

$$\begin{array}{c}{N}^{1-ϵ}Y=y>0,0<ϵ\le 1,N\to \infty .\hfill \end{array}$$

(31)

Then for every fixed $\gamma >0$ and $ϵ\in (0,1]$

$$\begin{array}{c}{\rho }_N^{\left(\mathrm{Im}\right)}\left(Y\right)\sim \frac{1}{\sqrt{N}}{\Psi }_{\gamma }\left(Y\right)\exp \left\{-N{\Phi }_{\gamma }\left(Y\right)\right\},\hfill \end{array}$$

(32)

with

$$\begin{array}{cc}{\Phi }_{\gamma }\left(Y\right)\hfill & =Y(\gamma -Y)-\ln \frac{\gamma -Y}{\gamma }-Y{r}_{\ast }\left(Y\right)+2\ln r_{\ast }\left(Y\right),\hfill \end{array}$$

(33)

$$\begin{array}{cc}{\Psi }_{\gamma }\left(Y\right)\hfill & =\frac{1}{\sqrt{2\pi }}\frac{\gamma }{{(\gamma -Y)}^2}\frac{{\left[1-r_{\ast }\left(Y\right)(\gamma -Y)\right]}^2}{Y^{3/2}{(Y^2+4)}^{1/4}},\hfill \end{array}$$

(34)

and

$$\begin{array}{c}{r}_{\ast }\left(Y\right)=\frac{\sqrt{Y^2+4}-Y}{2}.\hfill \end{array}$$

(35)

The rate function ${\Phi }_{\gamma }\left(Y\right)$ is a smooth non-negative function of Y on the interval $[0,\gamma )$ vanishing at $Y=0$. The rate function is monotone increasing on this interval if $\gamma \le 1$, whereas if $\gamma >1$ then it has two local extrema: a local minimum at $Y_{\ast }=\gamma -{\gamma }^{-1}$ where it vanishes, and a local maximum at $Y_{\ast \ast }=\frac{2(\gamma -{\gamma }^{-1})}{3+\sqrt{1+8{\gamma }^{-2}}}<Y_{\ast }$.

By the way of discussion of the above Theorem a few remarks are in order.

Remark 3.

The two distinct profiles of the rate function are illustrated in Figure 2. If$\gamma >1$, the point$Y_{\ast }=\gamma -{\gamma }^{-1}$where the Large Deviation Rate function${\Phi }_{\gamma }\left(Y\right)$vanishes can be identified as the most probable value of the imaginary part in the region$Y\gg N^{-1}$, converging in the limit$N\to \infty$to (the height of) the outlier, see next comment. At the same time, the other extremal point,$Y_{\ast \ast }$, can be interpreted as the true boundary, along the imaginary axis in the complex plane, between the bulk of eigenvalues and the spectral outlier. This is because the pre-exponential factor${\Psi }_{\gamma }\left(Y\right)$in (32) vanishes at $Y=Y_{\ast \ast }$ too. Hence, ${\rho }_N^{\left(\mathrm{Im}\right)}\left(Y_{\ast \ast }\right)\to 0$ in the scaling limit (31).

Remark 4.

The Large Deviation approximation (32) for $\gamma >1$ describes fluctuations of the imaginary part of the outlier around its most probable value $Y_{\ast }=\gamma -{\gamma }^{-1}$. The law of these fluctuations in the limit $N\to \infty$ can be easily determined from (32). To this end, we first note that for N large the magnitude of fluctuations about $Y_{\ast }$ scales with $1/\left(\sqrt{N}\right|{\Phi }_{\gamma }^{″}\left(Y_{\ast }\right)\left|\right)$. Calculating the second derivative and rescaling the density ${\rho }_N^{\left(\mathrm{Im}\right)}\left(Y\right)$ correspondingly, one finds (in the limit $N\to \infty$) that

$$\begin{array}{c}N{\rho }_N^{\left(\mathrm{Im}\right)}\left(Y_{\ast }+\frac{\sigma u}{\sqrt{N}}\right)\sim \frac{1}{\sqrt{2\pi }}{e}^{-\frac{u^2}{2}},{\sigma }^2=\frac{1}{{\gamma }^2}\frac{{\gamma }^2+1}{{\gamma }^2-1}.\hfill \end{array}$$

(36)

The integral of the rescaled density on the left-hand side over the entire range of values of u counts the expected number of eigenvalues in the $\frac{\sigma }{\sqrt{N}}$-neighbourhood of $Y_{\ast }$. Evidently, this integral is approaching unity as $N\to \infty$, confirming that the rescaled density on the left-hand side in (36) describes the law of fluctuations of a single eigenvalue - the outlier. Thus, we recover one of the results of [10] where laws of outlier fluctuations were established in greater generality than our assumptions (2) and (3). We note that for finite but large values of N the function

$$\begin{array}{c}{p}_N\left(Y\right):=\sqrt{N}{\Psi }_{\gamma }\left(Y\right)\exp \left\{-N{\Phi }_{\gamma }\left(Y\right)\right\}\hfill \end{array}$$

(37)

provides an approximation of the probability density function of the outlier $Y_{max}=\max Y_j$ in the interval $0<\epsilon <Y<\gamma$, $\gamma >1$.

In Figure 3, we plot histograms of the imaginary parts$Y_j$of the eigenvalues and of their maximal value$Y_{max}=\max Y_j$in the random matrix ensemble (1)–(3) and make comparison with the corresponding Large Deviation approximations. Although the value of $N=50$ is only moderately large, one can observe a good agreement. Furthermore, one can observe that the large-N approximation (37) of the probability density of $Y_{max}$ captures well the skewness of the distribution of $Y_{max}$ for finite matrix dimensions. This skewness disappears in the limit $N\to \infty$, see Equation (36).

Remark 5.

Consider now the scales $Y=O\left(N^{-1+\epsilon }\right)$ with $\epsilon \in (0,1)$. The expected number of eigenvalues with $N^{1-\epsilon }Y\in [y_1,y_2]$ is given by the integral

$$\begin{array}{c}N{\int }_{y_1}^{y_2}\frac{1}{N^{1-\epsilon }}{\rho }_N^{\left(\mathrm{Im}\right)}\left(\frac{y}{N^{1-\epsilon }}\right)dy.\hfill \end{array}$$

(38)

The rescaled density in this integral can be found from (32)–(34):

$$\begin{array}{c}\frac{1}{N^{1-\epsilon }}{\rho }_N^{\left(\mathrm{Im}\right)}\left(\frac{y}{N^{1-\epsilon }}\right)\sim \frac{1}{N^{\epsilon /2}}\frac{1}{2\sqrt{\pi }}\frac{{(1-\gamma )}^2}{\gamma }\frac{1}{y^{3/2}}{e}^{-N^{\epsilon }y\frac{{(1-\gamma )}^2}{\gamma }},\epsilon \in (0,1).\hfill \end{array}$$

(39)

Evidently, if $\gamma \ne 1$ then, away from the boundary point $y=0$, the integral in (38) vanishes in the limit $N\to \infty$. Therefore for every fixed $\gamma \ne 1$ and $0<\epsilon <1$ there are no eigenvalues of J whose imaginary part is scaling with $N^{-1+\epsilon }$. On the other hand, according to the heuristics of Section 2, one should expect finite numbers of eigenvalues whose imaginary part is scaling with $N^{-1}\ln N$. These would be the extremes of the eigenvalues with the typical imaginary part $Y=O\left(N^{-1}\right)$.

By formally letting$\epsilon \to 0$in (39) one obtains

$$\begin{array}{c}\frac{1}{N}{\rho }_N^{\left(\mathrm{Im}\right)}\left(\frac{y}{N}\right)\sim \frac{1}{2\sqrt{\pi }}\frac{{(1-\gamma )}^2}{\gamma }\frac{1}{y^{3/2}}{e}^{-y\frac{{(1-\gamma )}^2}{\gamma }}.\hfill \end{array}$$

This relation reproduces the leading order of the asymptotic form of the density of the rescaled imaginary parts $y=NY$ in the region $y\gg 1$, see the top line in (20). Thus, for a fixed value of $\gamma \ne 1$ Theorem 1 describes a crossover of the density of imaginary parts from the characteristic scale of low lying eigenvalues to larger scales, including $Y=O\left(1\right)$ which is the scale of the outlier.

Whereas the picture described by Theorem 1 is quite complete for a fixed $\gamma$, it is not detailed enough to accurately describe the typical extreme eigenvalues in the situation when the parameter $\gamma$ approaches its critical value $\gamma =1$ as N is approaching infinity. For example, from the heuristics of Section 2 we know that both the width of the transition region about $\gamma =1$ and the height of the typical extreme eigenvalues scale with $N^{-1/3}$. The Large Deviation approximation (32), if applied formally in the transition region parametrised by $\gamma =1+\alpha N^{-1/3}$, yields the following approximate expression for the rescaled density of imaginary parts:

$$\begin{array}{c}\frac{1}{N^{1/3}}{\rho }_N^{\left(\mathrm{Im}\right)}\left(\frac{m}{N^{1/3}}\right)\approx \frac{1}{N}\frac{1}{2\sqrt{\pi }}\frac{{\left(\frac{3m}{2}-\alpha \right)}^2}{m^{3/2}}{e}^{-m{(\alpha -\frac{m}{2})}^2}.\hfill \end{array}$$

(40)

Evidently, in the limit of small values of m which corresponds to approaching the scale $Y=O\left(N^{-1}\right)$ from above, this expression does not reproduce the correct power $5/2$ of algebraic decay (20) characteristic of this scale when $\gamma =1$. In contrast, the heuristics based on (21), see the approximations in (23) and (24), do reproduce the correct power. Indeed, by taking the derivative in m of the expression on the right-hand side in (24), one gets

$$\begin{array}{c}\frac{1}{N^{1/3}}{\rho }_N^{\left(\mathrm{Im}\right)}\left(\frac{m}{N^{1/3}}\right)\approx \frac{1}{N}\frac{1}{2\sqrt{\pi }}\frac{\frac{3}{2m}+{\alpha }^2}{m^{3/2}}{e}^{-m{\alpha }^2}.\hfill \end{array}$$

(41)

In the limit of small values of m the expression on the right-hand side agrees with the bottom line in (20). One can also arrive at (41) by making the formal substitution $\gamma =1+\frac{\alpha }{N^{1/3}}$ and $y=NY=m{N}^{2/3}$ in (20).

Our next Theorem is a refinement of Theorem 1 in that it provides an accurate description of the density of the typical extreme eigenvalues in the transition region between the sea of low lying eigenvalues and the eigenvalue outlier.

Theorem 2.

Consider the random matrix ensemble (1)–(3) in the scaling regime

$$\begin{array}{c}\gamma =1+\frac{\alpha }{N^{1/3}},Y=\frac{m}{N^{1/3}},N\to \infty .\hfill \end{array}$$

(42)

Then, for every fixed $\alpha \in \mathbb{R}$ and $m>0$,

$$\begin{array}{c}\frac{1}{N^{1/3}}{\rho }_N^{\left(\mathrm{Im}\right)}\left(\frac{m}{N^{1/3}}\right)\sim \frac{1}{N}\frac{1}{2\sqrt{\pi }}\frac{\frac{3}{2m}+{\left(\frac{3m}{2}-\alpha \right)}^2}{m^{3/2}}{e}^{-m{\left(\alpha -\frac{m}{2}\right)}^2}.\hfill \end{array}$$

(43)

This theorem confirms that the characteristic scale of the height of the typical extreme eigenvalues of matrix J is $O\left(N^{-1/3}\right)$. Indeed, the expected number of eigenvalues with imaginary part exceeding the level $Y=\frac{m}{N^{1/3}}$ is given by

$$\begin{array}{c}N{\int }_m^{\infty }\frac{1}{N^{1/3}}{\rho }_N^{\left(\mathrm{Im}\right)}\left(\frac{m}{N^{1/3}}\right)dm,\hfill \end{array}$$

which is a finite number in the limit $N\to \infty$.

Theorem 2 also describes the density ${\rho }_N^{\left(\mathrm{Im}\right)}\left(Y\right)$ in the cross-over from the characteristic scale of low lying eigenvalues to the Large Deviation regime of Theorem 1. Indeed, for small values of m the asymptotic expression (43) matches the one in (41), whilst in the limit of large values of m it matches (40).

The emerging outlier is captured by (43) when both m and $\alpha >0$ are large. Intuitively this is clear from the comparison of (43) and (40). On a more formal level, one can come to the same conclusion by analysing the limiting density of extreme values

$$\begin{array}{c}{\tilde{p}}_{\alpha }^{\left(\mathrm{Im}\right)}\left(m\right)=\frac{1}{2\sqrt{\pi }}\frac{\frac{3}{2m}+{\left(\frac{3m}{2}-\alpha \right)}^2}{m^{3/2}}{e}^{-m{\left(\alpha -\frac{m}{2}\right)}^2},m>0.\hfill \end{array}$$

(44)

Using Wolfram Mathematica one finds

$$\begin{array}{c}\frac{d}{dm}{\tilde{p}}_{\alpha }^{\left(\mathrm{Im}\right)}\left(m\right)=\frac{e^{-m{\left(\alpha -\frac{m}{2}\right)}^2}}{32\sqrt{\pi }{m}^{7/2}}{Q}_6(\alpha ,m),\hfill \end{array}$$

where

$$\begin{array}{c}{Q}_6(\alpha ,m)=-60-48{\alpha }^{2}m+72\alpha m^2-16{\alpha }^4{m}^2+80{\alpha }^3{m}^3-144{\alpha }^2{m}^4+108\alpha m^5-27m^6.\hfill \end{array}$$

Evidently, $Q_6(\alpha ,m)<0$ for all $m>0$ if $\alpha$ is negative. Therefore, if $\alpha <0$ (subcritical values of $\gamma$) then the limiting density ${\tilde{p}}_{\alpha }^{\left(\mathrm{Im}\right)}\left(m\right)$ is a monotonically decreasing function of m on the entire interval $m>0$. One can interpret this profile as a population of extreme eigenvalues without an obvious “leader”. By continuity, this profile persevere for small positive $\alpha$. Indeed, at $\alpha =0$ the polynomial $Q_6(0,m)$ has three pairs of complex conjugated roots, none are real. Since the roots of polynomials depend continuously on its coefficients, there exists an ${\alpha }_0>0$ such that for all $\alpha \in [0,{\alpha }_0]$ the polynomial $Q_6(\alpha ,m)$ in m will still have no real roots and, hence, will take only negative values, implying that ${\tilde{p}}_{\alpha }^{\left(\mathrm{Im}\right)}\left(m\right)$ is a monotonically decreasing function of m. By computing the roots of $Q_6(\alpha ,m)$ in variable m, we can show that $0.6485<{\alpha }_0<0.649$.

Once $\alpha >{\alpha }_0$, the polynomial $Q_6(\alpha ,m)$ in m acquires real roots. In the limit of large positive $\alpha$ there are two real roots: to leading order these are

$$\begin{array}{c}{m}_1=2\alpha \left(1+\frac{3}{8{\alpha }^3}+o\left(\frac{1}{{\alpha }^3}\right)\right)\mathrm{and}{m}_2=\frac{2}{3}\alpha \left(1+\frac{15}{8{\alpha }^3}+o\left(\frac{1}{{\alpha }^3}\right)\right).\hfill \end{array}$$

The larger root, $m_1$, is the point of local maximum of ${\tilde{p}}_{\alpha }^{\left(\mathrm{Im}\right)}\left(m\right)$, where ${\tilde{p}}_{\alpha }^{\left(\mathrm{Im}\right)}\left(m_1\right)\propto {\alpha }^{1/2}\gg 1$, and the smaller root, $m_2$, is the point of local minimum ${\tilde{p}}_{\alpha }^{\left(\mathrm{Im}\right)}\left(m\right)$, where ${\tilde{p}}_{\alpha }^{\left(\mathrm{Im}\right)}\left(m_2\right)\propto {\alpha }^{-5/2}\ll 1$. In fact, in the limit $\alpha \to \infty$ the larger root is transitioning into $Y_{\ast }$, the most probable value of imaginary parts, and, hence, it can be interpreted as the emerging spectral outlier. At the same time, the smaller root is transitioning into the true boundary $Y_{\ast \ast }$ between the sea of low lying eigenvalues and the outlier. This cross-over can be validated by noticing that in the scaling limit (42) $Y_{\ast }=\gamma -{\gamma }^{-1}\sim 2\alpha$ and $Y_{\ast \ast }=\frac{2(\gamma -{\gamma }^{-1})}{3+\sqrt{1+8{\gamma }^{-2}}}\sim \frac{2}{3}\alpha$.

Further insights into the restructuring of the spectrum of J can be obtained by looking at the $\gamma$-dependence of the expected number of the eigenvalues of J with imaginary parts exceeding the level $Y=m{N}^{-1/3}$. In the scaling limit (42) this number converges to

$$\begin{array}{c}{\tilde{\mathcal{N}}}_{\alpha }\left(m\right)={\int }_m^{\infty }{\tilde{p}}_{\alpha }^{\left(\mathrm{Im}\right)}\left(m^{\prime }\right)d{m}^{\prime }.\hfill \end{array}$$

In Figure 4, we plot ${\tilde{\mathcal{N}}}_{\alpha }\left(m\right)$ as function of $\alpha$ for several values of m. One can observe that for any fixed $m>0$ the population of the extreme eigenvalues of J that exceed the level $Y=m{N}^{-1/3}$ is, on average, growing as $\gamma$ is approaching the critical value $\gamma =1$ from below. For $\gamma$ on the other side of $\gamma =1$, this population peaks a some point and then it starts to decline as $\gamma$ increases further, to a single eigenvalue which is the outlier. All the other extreme eigenvalues are getting closer and closer to the real line with the increase of $\gamma$. One can think of them as being trapped in the sea of low lying eigenvalues. This picture is consistent with the eigenvalue trajectories of Figure 1 and provides a more quantitative description of the “resonance trapping” phenomenon [20] in the framework of random matrix theory.

Our final result aims to clarify the length of the central part of the spectrum of J supporting nontrivial scaling behaviour of the extreme eigenvalues in the vicinity of the separation transition. To this end, let us consider eigenvalues $z_j=X_j+i{Y}_j$ of J in the scaling regime when

$$\begin{array}{c}\gamma =1+\frac{\alpha }{N^{1/3}},X=\frac{q}{N^{1/3}},Y=\frac{m}{N^{1/3}},N\to \infty .\hfill \end{array}$$

(45)

On average, eigenvalue numbers in this regime can be counted using the rescaled density

$$\begin{array}{c}{\tilde{p}}_N(q,m):=〈{\displaystyle \sum _{i=1}^N}\delta \left(q-N^{1/3}{X}_j\right)\delta \left(m-N^{1/3}{Y}_j\right)〉=\frac{N}{N^{2/3}}{\rho }_N\left(\frac{q}{N^{1/3}},\frac{m}{N^{1/3}}\right),\hfill \end{array}$$

where, as before, the angle brackets stand for the averaging over the GUE matrix H in (1) and ${\rho }_N(X,Y)$ is the mean eigenvalue density (8).

Theorem 3.

Consider the random matrix ensemble (1)–(3) in the scaling regime (45). Then, for every fixed $\alpha \in \mathbb{R}$, $q\in \mathbb{R}$ and $m>0$,

$$\begin{array}{c}{\tilde{p}}_{\alpha }(q,m):=\underset{N\to \infty }{\lim }{\tilde{p}}_N(q,m)=\frac{1}{4\pi m}\left[\frac{1}{m}+\frac{q^2}{4}+{\left(\frac{3m}{2}-\alpha \right)}^2\right]e^{-m\left[\frac{q^2}{4}+{\left(\alpha -\frac{m}{2}\right)}^2\right]}.\hfill \end{array}$$

(46)

It is easy to see from (46) that ${\int }_{-\infty }^{\infty }\tilde{p}(q,m)dq={\tilde{p}}_{\alpha }^{\left(\mathrm{Im}\right)}\left(m\right)$. Thus, Theorem 3 confirms the heuristics of Section 2 in that the population of extreme eigenvalues which generates the eventual outlier (as $\alpha$ is approaching infinity) is constrained to a narrow vertical strip of width $O\left(N^{-1/3}\right)$ about the origin.

Our results demonstrate that despite being one of the simplest tools available, the mean eigenvalue density captures the eigenvalue and parameter scales associated with the spectral restructuring in the random matrix ensemble (1)–(3). However, it gives no information about finer details, such as the probability distribution of the extreme eigenvalues during the restructure. Calculating all the higher order eigenvalue correlation functions in the scaling regime (45) would be a significant step towards describing such finer details. Unfortunately, the eigenvalue point process in the random matrix ensemble (1)–(3) is not determinantal at finite matrix dimensions and such a calculation is a considerably more difficult analytic task compared to the mean eigenvalue density.

At this point we want to mention that the probability distribution of extreme eigenvalues can be determined in a related but different random matrix ensemble exhibiting a spectral restructuring not unlike one in (1)–(3). This ensemble consists of subunitary matrices of the form

$$\begin{array}{c}{J}_{CUE}=U\mathrm{diag}(\sqrt{1-T},1,\cdots ,1),\hfill \end{array}$$

(47)

where the matrix U is taken from the Circular Unitary Ensemble (CUE) of complex unitary matrices uniformly distributed over $U\left(N\right)$ with the Haar’s measure and $T\in [0,1]$ is a parameter. The ensemble was originally introduced in [29] and various statistical aspects of their spectra and eigenvectors were addressed in [2,6,30,31] and most recently in [15].

Obviously, if $T=0$ then the matrix $J_{CUE}$ is unitary and all of its eigenvalues lie on the unit circle $\left|z\right|=1$. If $T>0$ and is fixed in the limit $N\gg 1$ then, typically, the eigenvalues of $J_{CUE}$ lie at a distance $O\left(N^{-1}\right)$ from the unit circle with the farthest away being at a distance $O\left(\frac{\log N}{(1-T)N}\right)$ with probability close to one. On the other hand, for $T=1$ one of the eigenvalues becomes identically zero, and the rest are distributed inside the unit circle in the same way as eigenvalues of the so-called “truncated” CUE [32].

The similarity between the random matrix ensembles (47) and (1)–(3) can be exemplified by analysing the mean density of the eigenvalue moduli $r_j=\left|z_j\right|$

$$\begin{array}{c}{\rho }_N\left(r\right)={〈\frac{1}{N}{\displaystyle \sum _{j=1}^N}\delta (r-r_j)〉}_{CUE}\hfill \end{array}$$

in the limit of large matrix dimensions $N\to \infty$. One finds [29] that for every fixed $T\in [0,1]$

$$\begin{array}{c}\underset{N\to \infty }{\lim }N{\rho }_N\left(r\right)=\left\{\begin{array}{cc}\delta \left(r\right)+\frac{2r}{{(1-r^2)}^2},\hfill & \mathrm{if}T=1,\hfill \\ 0,\hfill & \mathrm{if}0<T<1,\hfill \end{array}\right.\hfill \end{array}$$

whereas, on rescaling the radial density near the unit circle [2,32],

$$\begin{array}{c}{\tilde{\rho }}_{CUE}\left(y\right):=\underset{N\to \infty }{\lim }\frac{1}{N}{\rho }_N\left(1-\frac{y}{N}\right)=-\frac{d}{dy}\left[e^{-gy}\frac{\sinh y}{y}\right],\mathrm{with}g=\frac{2}{T}-1.\hfill \end{array}$$

(48)

Equation (48) is identical, with the obvious correspondence

$$\begin{array}{c}\frac{1}{2}\left(\gamma +\frac{1}{\gamma }\right)=\frac{2}{T}-1,\hfill \end{array}$$

(49)

to Equation (12) considered at the centre of the GUE spectrum. In the limit of large values of y,

$$\begin{array}{c}{\tilde{\rho }}_{CUE}\left(y\right)\sim \left\{\begin{array}{cc}{\displaystyle \frac{1-T}{T}\frac{1}{y}{e}^{-2y\frac{1-T}{T}},}\hfill & \mathrm{if}0<T<1,\hfill \\ {\displaystyle \frac{1}{y^2}}\hfill & \mathrm{if}T=1.\hfill \end{array}\right.\hfill \end{array}$$

(50)

The rescaled radial density has an exponentially light tail if $0<T<1$, and it is heavy-tailed if $T=1$ which hints at markedly different behaviour of the extreme eigenvalues in the two cases. Reflecting on (50), one can convince themselves that this change occurs in an infinitesimal region near $T=1$ of width $N^{-1}$. Such a scaling regime was earlier identified and analysed from a somewhat different angle in [15]. The precise relation of our analysis to one in [15] will be given in a separate paper [33].

On setting $T=1-\frac{t}{N}$, $t>0$, one can investigate this transition region in much detail [33]. For example, the smallest eigenvalue modulus of the subunitary matrices $J_{CUE}$,

$$\begin{array}{c}{x}_{min}=\underset{j=1,\dots ,N}{min}\left|z_j\right|,\hfill \end{array}$$

converges in the limit $N\to \infty$ to a random variable X whose cumulative probability distribution function is given by the series

$$\begin{array}{c}\mathrm{Pr}\{X\le x\}={\displaystyle \sum _{n=1}^{\infty }}{(-1)}^{n+1}\frac{x^{n(n-1)}}{{\prod }_{k=1}^n(1-x^{2k})}{e}^{t\left(1-\frac{1}{x^{2n}}\right)},0<x<1.\hfill \end{array}$$

This family of probability distributions interpolates between the Fréchet and Gumbel distributions and is different from the standard family of probability distributions that characterise the extreme values in long sequence of i.i.d. random variables. In the limit of small values of t

$$\begin{array}{c}\underset{t\to 0+}{\lim }\mathrm{Pr}\left\{X<y\sqrt{t}\right\}=e^{-y^{-2}},y>0,\hfill \end{array}$$

whereas

$$\begin{array}{c}\underset{t\to +\infty }{\lim }\mathrm{Pr}\left\{2t(1-X)-\log t+\log \left(\log t\right)<y\right\}=e^{-e^{-y}},y>0.\hfill \end{array}$$

4. Mean Density of Eigenvalues at Finite Matrix Dimensions

Our analysis of various scaling regimes of the random matrix ensemble (1)–(3) is based on finite-N expressions for the mean eigenvalue density and the mean density of imaginary parts in terms of orthogonal polynomials, see Equations (56)–(64). These representations are new and the current Section contains their derivations.

4.1. Joint Eigenvalue Density and Correlation Functions

Our starting point is a closed form expression for the joint density $P_N(z_1,\dots ,z_n)$ of the eigenvalues $z_k=X_k+i{Y}_k$ of J (1)–(3):

$$\begin{array}{l}{P}_N(z_1,\dots ,z_N)=\\ \frac{N^{N^2/2}}{{\left(2\pi \right)}^{N/2}N!G\left(N\right){\gamma }^{N-1}}\exp \left\{-\frac{N}{2}\left({\gamma }^2+{\displaystyle \sum _{k=1}^N}\mathrm{Re}\left(z_k^2\right)\left)\right\}\right.\right.\delta (\gamma -{\displaystyle \sum _{k=1}^N}\mathrm{Im}{z}_k){\displaystyle \prod _{j<k}^N}{\left|z_j-z_k\right|}^2,\end{array}$$

(51)

where G(N) is the Barnes G-function. The expression was derived in [22] (see also [5]) and, for obvious reason, it holds for $(z_1,\dots ,z_N)\in {\mathbb{C}}_{+}^N$, where ${\mathbb{C}}_{+}$ is the upper half of the complex plane ${\mathbb{C}}_{+}=\{z=X+iY:Y\ge 0\}$

The first key fact that makes our analysis possible is that the eigenvalue correlation functions

$$\begin{array}{c}{\displaystyle R_{N,n}\left(z_1,\dots ,z_n\right)={\displaystyle \frac{N!}{\left(N-n\right)!}}\underset{{\mathbb{C}}_{+}^{N-n}}{\int }{P}_N(z_1,\dots ,z_n,z_{n+1},\dots ,z_N){\displaystyle \prod _{k=n+1}^N}d{X}_{k}d{Y}_k,}\hfill \end{array}$$

can be expressed in terms of averages of products of characteristic polynomials of random matrices $J\left(\tilde{\gamma }\right)$ having the same structure as (1)–(3) but of smaller dimension and with a different parameter $\gamma$. The relevance of this to our investigation is in that the mean eigenvalue density ${\rho }_N(X,Y)$ (8) which is the main object of our interest is

$$\begin{array}{c}{\rho }_N(X,Y)=\frac{1}{N}{R}_{N,1}(X+iY).\end{array}$$

(52)

It has been shown in [22] that

$$\begin{array}{c}{R}_{N,n}\left(z_1,\dots ,z_n\right)={\displaystyle \frac{1}{{\left(2\pi \right)}^{n/2}{\gamma }^n}}{\left(1-\frac{{\sum }_{k=1}^n{Y}_k}{\gamma }\right)}^{N-n-1}\frac{N^{\frac{n^2}{2}}{(N-n)}^{Nn-n^2}}{{\displaystyle {\displaystyle \prod _{j=1}^n}\left(N-j-1\right)!}}\prod _{1\le j<k\le n}{\left|z_j-z_k\right|}^2\times \\ \exp \left\{-\frac{N}{2}{\displaystyle \sum _{k=1}^n}{X}_k^2-N{\displaystyle \sum _{k=1}^n}{Y}_k\left(\gamma -Y_k\right)\right\}{〈{\displaystyle \prod _{k=1}^n}{\left|\det \left[{\widehat{z}}_k{\mathbf{1}}_{N-n}-J_{\widehat{\gamma }-{\sum }_{k=1}^n{\widehat{Y}}_k}\right]\right|}^2〉}_{H_{N-n}},\end{array}$$

where

$$\begin{array}{c}\widehat{\gamma }={\left(\frac{N}{N-n}\right)}^{1/2}\gamma ,{\widehat{z}}_k={\left(\frac{N}{N-n}\right)}^{1/2}(X_k+i{Y}_k),{\widehat{Y}}_k={\left(\frac{N}{N-n}\right)}^{1/2}{Y}_k,\hfill \end{array}$$

and $J_{\widehat{\gamma }-{\sum }_{k=1}^n{\widehat{Y}}_k}$ are the random matrices (1)–(3) of dimension $N-n$ with N in (3) replaced by $N-n$ and $\gamma$ in (2) replaced by $\widehat{\gamma }-{\sum }_{k=1}^n{\widehat{Y}}_k$,

$$\begin{array}{c}{J}_{\widehat{\gamma }-{\sum }_{k=1}^n{\widehat{Y}}_k}=H_{N-n}+i\left(\widehat{\gamma }-{\displaystyle \sum _{k=1}^n}{\widehat{Y}}_k\right)\mathrm{diag}(1,0,\dots ,0).\hfill \end{array}$$

The GUE average ${〈\dots 〉}_{H_{N-n}}$ of the product of the characteristic polynomials of $J_{\widehat{\gamma }-{\sum }_{k=1}^n{\widehat{Y}}_k}$ can be performed with the help of the following proposition which we prove in Appendix A.

Proposition 1.

Let

$$\begin{array}{c}{F}_{\gamma }\left(z_1,z_2,\dots ,z_n\right)=〈{\displaystyle {\displaystyle \prod _{j=1}^n}{\left|\det \left(z_j{\mathbf{1}}_N-J_{\gamma }\right)\right|}^2}〉,\hfill \end{array}$$

where$J_{\gamma }$are the rank-one deviations from the GUE of dimension N defined by (1)–(3) and the average is taken over the GUE distribution (3). Then

$$\begin{array}{cc}{F}_{\gamma }(z_1,z_2,\dots ,z_n)\hfill & =\hfill \\ {\displaystyle \frac{1}{2^n}{\left(\frac{N}{\pi }\right)}^{2n^2}\int D\left[S_{2n}\right]\exp \left\{-\frac{N}{2}\mathrm{Tr}{S}_{2n}^2\right\}{\det }^{N-1}\left(Z_{2n}+i{S}_{2n}\right)\det \left(Z_{2n}+i{S}_{2n}-i\gamma L_{2n}\right),}\hfill \end{array}$$

where the integration is over the space of $2n\times 2n$ Hermitian matrices $S_{2n}$, $D\left[S_{2n}\right]$ is the standard volume element in this space and

$$\begin{array}{c}{Z}_{2n}=\mathrm{diag}\left(z_1,z_2,\dots ,z_n,{\overline{z}}_1,{\overline{z}}_2,\dots ,{\overline{z}}_n\right),L_{2n}=\mathrm{diag}\left(1,-1\right)\otimes {\mathbf{1}}_n.\hfill \end{array}$$

Using this Proposition one arrives, after rescaling $S_{2n}={\left(\frac{N}{N-n}\right)}^{1/2}{\widehat{S}}_{2n}$ in the resulting matrix integral, at a useful integral representation for the eigenvalue correlation functions in the random matrix ensemble (1)–(3):

$$\begin{array}{cc}{R}_{N,n}\left(z_1,\dots ,z_n\right)=\hfill & \\ {\displaystyle \frac{c_N}{{\gamma }^n}{\left(1-\frac{{\sum }_{k=1}^n{Y}_k}{\gamma }\right)}^{N-n-1}\exp \left\{-\frac{N}{2}{\displaystyle \sum _{k=1}^n}{X}_k^2-N{\displaystyle \sum _{k=1}^n}{Y}_k\left(\gamma -Y_k\right)\right\}\prod _{1\le j<k\le n}{\left|z_j-z_k\right|}^2\times }\hfill \\ {\displaystyle \int d\left[{\widehat{S}}_{2n}\right]\exp \left\{-\frac{N}{2}\mathrm{Tr}{\widehat{S}}_{2n}^2\right\}{\det }^{N-n-1}\left[Z_{2n}+i{\widehat{S}}_{2n}\right]\det \left[Z_{2n}+i{\widehat{S}}_{2n}-i\left(\gamma -{\displaystyle \sum _{k=1}^n}{Y}_k\right)L_{2n}\right]}\hfill \end{array}$$

(53)

with

$$\begin{array}{c}{c}_{N,n}\left(\gamma \right)={\displaystyle \frac{N^{3n^2/2+Nn}}{{\left(2\gamma \right)}^n{\left(2\pi \right)}^{n/2}{\pi }^{2n^2}{\displaystyle {\displaystyle \prod _{j=1}^n}}\left(N-j-1\right)!}}.\hfill \end{array}$$

4.2. Mean Density of Complex Eigenvalues

Setting $n=1$ and $z_1=X+iY$ in (53) and then shifting the variable of integration by making the substitution ${\widehat{S}}_2=S_2-Y{L}_2$ in the matrix integral, one obtains the following integral representation for the mean density of eigenvalues (52) in the random matrix ensemble (1)–(3):

$$\begin{array}{cc}\hfill & {\rho }_N(X,Y)={\displaystyle \frac{1}{2\gamma {\pi }^2\sqrt{2\pi }}}{\displaystyle \frac{N^{N+1/2}}{\left(N-2\right)!}}{\left(1-\frac{Y}{\gamma }\right)}^{N-2}\exp \left\{-\frac{N}{2}{X}^2-NY\gamma \right\}\times \hfill \\ \hfill & {\displaystyle \int D\left[S_2\right]\exp \left\{-\frac{N}{2}\mathrm{Tr}{S}_2^2+NY\mathrm{Tr}{S}_2{L}_2\right\}{\det }^{N-2}(X{\mathbf{1}}_2+i{S}_2)\det \left(X{\mathbf{1}}_2+i{S}_2-i(\gamma -Y)L_2\right).}\hfill \end{array}$$

(54)

It is convenient to parametrise the hermitian matrix $S_2$ by diagonalising it:

$$\begin{array}{c}{S}_2=U_2{\Sigma }_2{U}_2^{\ast },{\Sigma }_2=\mathrm{diag}\left({\sigma }_1,{\sigma }_2\right),{\sigma }_1\ge {\sigma }_2\in \mathbb{R},\hfill \end{array}$$

where $U_2$ is a $2\times 2$ unitary matrix, which can be parametrised as

$$\begin{array}{c}{U}_2=\left(\begin{array}{cc}\cos \theta & \sin \theta e^{i\varphi }\\ -\sin \theta e^{-i\varphi }& \cos \theta \end{array}\right),\theta \in \left[0,\frac{\pi }{2}\right],\varphi \in \left[0,2\pi \right].\hfill \end{array}$$

Noting that

$$\begin{array}{c}D\left[S_2\right]={\left({\sigma }_1-{\sigma }_2\right)}^2\frac{\sin \left(2\theta \right)}{2}d{\sigma }_{1}d{\sigma }_{2}d\theta d\varphi .\hfill \end{array}$$

one arrives, on making the substitution $S_2=U_2{\mathsf{\Sigma }}_2{U}_2^{\ast }$ in (54), at

$$\begin{array}{l}{\rho }_N(X,Y)=\frac{c_N}{\gamma }{\left(1-\frac{Y}{\gamma }\right)}^{N-2}\exp \left\{-\frac{N}{2}\left(X^2+2Y\gamma \right)\right\}\times \\ {\int }_0^{\frac{\pi }{2}}d\theta {\int }_{-\infty }^{+\infty }d{\sigma }_1{\int }_{-\infty }^{{\sigma }_1}d{\sigma }_2{({\sigma }_1-{\sigma }_2)}^2\sin \left(2\theta \right)\exp \left\{-\frac{N}{2}\left({\sigma }_1^2+{\sigma }_2^2\right)+NY\left({\sigma }_1-{\sigma }_2\right)\cos \left(2\theta \right)\right\}\times \\ {\left(X+i{\sigma }_1\right)}^{N-2}{\left(X+i{\sigma }_2\right)}^{N-2}\left[\left(X+i{\sigma }_1\right)\left(X+i{\sigma }_2\right)+{\left(\gamma -Y\right)}^2-\left(\gamma -Y\right)\left({\sigma }_1-{\sigma }_2\right)\cos \left(2\theta \right)\right],\end{array}$$

where we have introduced

$$\begin{array}{c}{c}_N={\displaystyle \frac{1}{{\left(2\pi \right)}^{3/2}}}{\displaystyle \frac{N^{N+1/2}}{\left(N-2\right)!}}\sim {\displaystyle \frac{N^2{e}^N}{2{\pi }^2}}(N\to \infty ).\hfill \end{array}$$

(55)

The integral over $\theta$ can be performed by the substitution $t=({\sigma }_1-{\sigma }_2)\cos \left(2\theta \right)$. This yields

$$\begin{array}{cc}{\rho }_N(X,Y)=\hfill & {\displaystyle \frac{c_N}{2NY\gamma }}{\left(1-\frac{Y}{\gamma }\right)}^{N-2}\exp \left\{-\frac{N}{2}\left(X^2+2\gamma Y-2Y^2\right)\right\}{J}_N(X,Y),\hfill \end{array}$$

(56)

where

$$\begin{array}{cc}{J}_N(X,Y)=\hfill & {\int }_{-\infty }^{+\infty }d{\sigma }_1{\int }_{-\infty }^{+\infty }d{\sigma }_2{e}^{-\frac{N}{2}({\sigma }_1^2+{\sigma }_2^2)}{(z+i{\sigma }_1)}^{N-2}{(\overline{z}+i{\sigma }_2)}^{N-2}{\displaystyle \frac{(z+i{\sigma }_1)-(\overline{z}+i{\sigma }_2)}{i}}\times \hfill \\ \hfill & \left[(z+i{\sigma }_1)(\overline{z}+i{\sigma }_2)+{(\gamma -Y)}^2+\frac{\gamma -Y}{Ny}-(\gamma -Y){\displaystyle \frac{(z+i{\sigma }_1)-(\overline{z}+i{\sigma }_2)}{i}}\right],\hfill \end{array}$$

with $z=X+iY$.

Further, introducing functions

$$\begin{array}{c}{\pi }_m\left(z\right)={\int }_{-\infty }^{+\infty }d\sigma e^{-\frac{N}{2}{\sigma }^2}{(z+i\sigma )}^{N-m},m=0,1,\dots ,N,\hfill \end{array}$$

(57)

one can rewrite the integral $J_N(X,Y)$ in the following form

$$\begin{array}{cc}{J}_N(X,Y)=\hfill & -i\left[{\pi }_0\left(Z\right){\pi }_1\left(\overline{z}\right)-{\pi }_0\left(\overline{z}\right){\pi }_1\left(z\right)\right]\hfill \\ \hfill & -i\left({(\gamma -Y)}^2+\frac{\gamma -Y}{Ny}\right)\left[{\pi }_1\left(z\right){\pi }_2\left(\overline{z}\right)-{\pi }_1\left(\overline{z}\right){\pi }_2\left(z\right)\right]\hfill & \\ \hfill & +(\gamma -Y)\left[{\pi }_0\left(z\right){\pi }_2\left(\overline{z}\right)+{\pi }_0\left(\overline{z}\right){\pi }_2\left(z\right)-2{\pi }_1\left(z\right){\pi }_1\left(\overline{z}\right)\right],\hfill \end{array}$$

(58)

Now one observes that ${\pi }_m\left(z\right)$ are actually a rescaled version of Hermite polynomials. We have that

$$\begin{array}{c}{\pi }_m\left(z\right)=\sqrt{\pi }{\left({\displaystyle \frac{2}{N}}\right)}^{\frac{N-m+1}{2}}{\tilde{H}}_{N-m}\left(z\sqrt{\frac{N}{2}}\right)=\sqrt{2\pi }{\pi }^{1/4}\sqrt{{\displaystyle \frac{\left(N-m\right)!}{N^{N-m+1}}}}{p}_{N-m}\left(z\sqrt{\frac{N}{2}}\right),\hfill \end{array}$$

(59)

where ${\tilde{H}}_k\left(z\right)$ are the monic Hermite polynomials

$$\begin{array}{c}{\tilde{H}}_k\left(z\right)={\left(-\frac{1}{2}\right)}^k{e}^{z^2}\frac{d}{dz}{e}^{-z^2}\hfill \end{array}$$

and $p_k\left(z\right)$ are the orthonormal Hermite polynomials

$$\begin{array}{c}{p}_k\left(z\right)=\sqrt{{\displaystyle \frac{2^k}{k!\sqrt{\pi }}}}{\tilde{H}}_k\left(z\right)\hfill \end{array}$$

satisfying the orthogonality relations

$$\begin{array}{c}{\int }_{-\infty }^{+\infty }dz{p}_k\left(z\right)p_m\left(z\right)e^{-z^2}dz={\delta }_{k,m}.\hfill \end{array}$$

The polynomials $p_k\left(z\right)$ also satisfy the recurrence relation

$$\begin{array}{c}{p}_{k+1}\left(z\right)=z\sqrt{\frac{2}{k+1}}{p}_k\left(z\right)-\sqrt{\frac{k}{k+1}}{p}_{k-1}\left(z\right).\hfill \end{array}$$

Using the above definitions and the expression for the eigenvalue density ${\rho }_N(X,Y)$ in (56) and with the notation $z=X+iY$ we obtain

$$\begin{array}{cc}{\rho }_N(X,Y)=\hfill & \frac{N-1}{\sqrt{2N}Y\gamma }{\left(1-\frac{Y}{\gamma }\right)}^{N-2}\exp \left\{-\frac{N}{2}{X}^2-NY\left(\gamma -Y\right)\right\}\times \hfill \\ \hfill & \left\{\mathrm{Im}{p}_N\left(z\sqrt{\frac{N}{2}}\right)p_{N-1}\left(\overline{z}\sqrt{\frac{N}{2}}\right)-\left(\gamma -Y\right){\left|p_{N-1}\left(z\sqrt{\frac{N}{2}}\right)\right|}^2\right.\hfill \\ \hfill & +\sqrt{\frac{N}{N-1}}\left[\left({\left(\gamma -Y\right)}^2+{\displaystyle \frac{\gamma -Y}{NY}}\right)\mathrm{Im}{p}_{N-1}\left(z\sqrt{\frac{N}{2}}\right)p_{N-2}\left(\overline{z}\sqrt{\frac{N}{2}}\right)\right.\hfill \\ & \left.\left.+\left(\gamma -Y\right)\mathrm{Re}{p}_N\left(z\sqrt{\frac{N}{2}}\right)p_{N-2}\left(\overline{z}\sqrt{\frac{N}{2}}\right)\right]\right\}.\hfill \end{array}$$

which, by using the recurrence relation, can be further rewritten as

$$\begin{array}{cc}{\rho }_N(X,Y)=\hfill & \frac{1}{Y\gamma }\sqrt{\frac{N}{2}}{\left(1-\frac{Y}{\gamma }\right)}^{N-2}\exp \left\{-\frac{N}{2}{X}^2-NY\left(\gamma -Y\right)\right\}\times \hfill \\ \hfill & \left\{\mathrm{Im}{p}_N\left(z\sqrt{\frac{N}{2}}\right)p_{N-1}\left(\overline{z}\sqrt{\frac{N}{2}}\right)\left[1-\frac{1}{N}+\left(\gamma -Y\right)\left(\gamma +\frac{1}{NY}\right)\right]\right.\hfill \\ \hfill & -{\left|p_{N-1}\left(z\sqrt{\frac{N}{2}}\right)\right|}^2\left[Y{\left(\gamma -Y\right)}^2+\left(\gamma -Y\right)-Y\left(\gamma -Y\right)\right]\hfill \\ \hfill & \left.-{\left|p_N\left(z\sqrt{\frac{N}{2}}\right)\right|}^2\left(\gamma -Y\right)+{\mathrm{Re}}_{pN}\left(z\sqrt{\frac{N}{2}}\right)p_{N-1}\left(\overline{z}\sqrt{\frac{N}{2}}\right)X\left(\gamma -Y\right)\right\}.\hfill \end{array}$$

(60)

4.3. Density of the Imaginary Parts

In this section, we present the derivation of the density for the imaginary parts of the eigenvalues, irrespective of their real parts, as defined in (10). We start with an observation, see integral 7.377 in [23]:

Lemma 1.

Let $\beta \ge \alpha$ be two non-negative integers and $z=X+iY$. Then

$$\begin{array}{c}{\displaystyle {\displaystyle \underset{-\infty }{\overset{\infty }{\int }}}{e}^{-\frac{N}{2}{X}^2}{p}_{N-\alpha }\left(z\sqrt{\frac{N}{2}}\right)p_{N-\beta }\left(\overline{z}\sqrt{\frac{N}{2}}\right)dX=}\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{i}^{\beta -\alpha }\sqrt{\frac{2}{N}\frac{N^{\beta -\alpha }}{\left(N-\alpha \right)\dots \left(N-\beta +1\right)}}{Y}^{\beta -\alpha }{L}_{N-\beta }^{\left(\beta -\alpha \right)}\left(-N{Y}^2\right),\end{array}$$

where $L_M^{\left(\alpha \right)}$ is a standard Laguerre polynomial.

Integrating with respect to X expression for the density ${\rho }_N(X,Y)$ in (60) one gets the probability density of imaginary parts in the form

$$\begin{array}{c}{\rho }_N^{\left(\mathrm{Im}\right)}\left(Y\right)={\displaystyle \frac{1}{Y\gamma }}{\left(1-\frac{Y}{\gamma }\right)}^{N-2}{e}^{-NY\left(\gamma -Y\right)}{F}_N\left(Y\right)\hfill \end{array}$$

(61)

with

$$\begin{array}{ccc}{F}_N\left(Y\right)=\hfill & \frac{N-1}{N}Y{L}_{N-1}^{\left(1\right)}\left(-N{Y}^2\right)-\frac{N-1}{N}\left(\gamma -Y\right)L_{N-1}^{\left(0\right)}\left(-N{Y}^2\right)\hfill & \\ \hfill & +Y\left[{\left(\gamma -Y\right)}^2+\frac{\gamma -Y}{NY}\right]L_{N-2}^{\left(1\right)}\left(-N{Y}^2\right)-\left(\gamma -Y\right)Y^2{L}_{N-2}^{\left(2\right)}\left(-N{Y}^2\right)\hfill & \end{array}$$

(62)

$$\begin{array}{ll}=-2\gamma L_{N-1}^{\left(0\right)}\left(-N{Y}^2\right)+\left(\frac{N-1}{N}3Y+\frac{2\gamma }{N}\right)L_{N-1}^{\left(1\right)}\left(-N{Y}^2\right)& \\ \hspace{1em}+\left[-2Y+Y{(\gamma -Y)}^2\right]L_{N-2}^{\left(1\right)}\left(-N{Y}^2\right)& \end{array}$$

(63)

$$\begin{array}{cc}=\hfill & \frac{N-1}{N}\left(3Y-2\gamma \right)L_{N-1}^{\left(1\right)}\left(-N{Y}^2\right)+\left[2\gamma -2Y+Y{\left(\gamma -Y\right)}^2\right]L_{N-2}^{\left(1\right)}\left(-N{Y}^2\right)\hfill \end{array}$$

(64)

where we systematically used the recursion relations:

$$\begin{array}{c}{L}_{N-1}^{\left(0\right)}\left(-N{Y}^2\right)=L_{N-1}^{\left(1\right)}\left(-N{Y}^2\right)-L_{N-2}^{\left(1\right)}\left(-N{Y}^2\right)\hfill \end{array}$$

and

$$\begin{array}{c}-Y^2{L}_{N-2}^{\left(2\right)}\left(-N{Y}^2\right)=L_{N-2}^{\left(1\right)}\left(-N{Y}^2\right)-\frac{N-1}{N}{L}_{N-1}^{\left(1\right)}\left(-N{Y}^2\right).\hfill \end{array}$$

5. Proof of Theorems 1 and 2

In both proofs we use the following integral representation for the Laguerre polynomials in terms of the modified Bessel functions $I_{\alpha }\left(x\right)$ (see, e.g., Equation 4.19.13 in [34]):

$$\begin{array}{c}{\displaystyle L_{N-k}^{\left(\alpha \right)}\left(-N{Y}^2\right)=\frac{2N^{N-k+1}}{\left(N-k\right)!}\frac{e^{-N{Y}^2}}{{\left|Y\right|}^{\alpha }}{\displaystyle \underset{0}{\overset{\infty }{\int }}}{\tau }^{2N-2k+\alpha +1}{e}^{-N{\tau }^2}{I}_{\alpha }\left(2\tau \right|Y\left|N\right)d\tau (\alpha >-1).}\hfill \end{array}$$

(65)

The integral in (65) can be evaluated in the limit $N\to \infty$ in various scaling regimes for Y using the Laplace method, see Appendix B. The resulting asymptotic expression depends on the scaling of the variable $Y>0$ with N.

Proof of Theorem 1.

Consider the scaling regime (31) with $\gamma >0$ being fixed. In this regime the asymptotic form of the mean density of the imaginary parts can be found using the leading order form of $L_{N-k}^{\left(1\right)}\left(-N{Y}^2\right)$ which can be read from (A8) as

$$\begin{array}{c}{L}_{N-k}^{\left(1\right)}\left(-N{Y}^2\right)\sim \frac{e^{NY{r}_{\ast }}}{\sqrt{2\pi N}}\frac{r_{\ast }{\left(Y\right)}^{-2(N-k+1)}}{Y^{3/2}{\left(Y^2+4\right)}^{1/4}},r_{\ast }\left(Y\right)=\frac{\sqrt{Y^2+4}-Y}{2}.\hfill \end{array}$$

(66)

On substituting (66) into (64) one gets an asymptotic expression for the density (61) precisely in the Large Deviation form (32) with the rate function (33) and the pre-exponential factor in the form

$$\begin{array}{c}\frac{1}{\sqrt{N}}{\Psi }_{\gamma }\left(Y\right)=\frac{1}{\sqrt{2\pi N}}\frac{\gamma }{{(\gamma -Y)}^2}\frac{3Y-2\gamma +r_{\ast }{\left(Y\right)}^2(\gamma -Y)\left(2+Y(\gamma -Y)\right)}{Y^{5/2}{(Y^2+4)}^{1/4}}.\hfill \end{array}$$

(67)

Finally, by exploiting the relation $1-r_{\ast }{\left(Y\right)}^2=Y{r}_{\ast }\left(Y\right)$,

$$\begin{array}{cc}\hfill 3Y-2\gamma +r_{\ast }{\left(Y\right)}^2(\gamma -Y)\left(2+Y(\gamma -Y)\right)=& \hfill Y-(\gamma -Y)\left[2(1-r_{\ast }{\left(Y\right)}^2)-r_{\ast }{\left(Y\right)}^{2}Y(\gamma -Y)\right]\\ \hfill =& Y-(\gamma -Y)\left[2r_{\ast }\left(Y\right)Y-r_{\ast }{\left(Y\right)}^{2}Y(\gamma -Y)\right]\hfill \\ \hfill =& Y\left[1-r_{\ast }\left(Y\right)(\gamma -Y)\left(2-r_{\ast }\left(Y\right)(\gamma -Y)\right)\right].\hfill \end{array}$$

This brings the function ${\Psi }_{\gamma }\left(Y\right)$ in (67) to the form as given in (34).

To analyse the shape of the rate function ${\Phi }_{\gamma }\left(Y\right)$ in (33) it is convenient to parametrise

$$\begin{array}{c}Y=e^{\theta }-e^{-\theta },\theta >0.\hfill \end{array}$$

(68)

In this parametrisation, the rate function transforms to

$$\begin{array}{c}{\tilde{\Phi }}_{\gamma }\left(\theta \right):={\Phi }_{\gamma }\left(e^{\theta }-e^{-\theta }\right)=\gamma \left(e^{\theta }-e^{-\theta }\right)+1-e^{2\theta }-2\theta -\ln \left(1-\frac{e^{\theta }-e^{-\theta }}{\gamma }\right),\hfill \end{array}$$

and its derivative in $\theta$ factorises as follows:

$$\begin{array}{cc}{\tilde{\Phi }}_{\gamma }^{\prime }\left(\theta \right)=\hfill & \gamma \left(e^{\theta }+e^{-\theta }\right)-2\left(e^{2\theta }+1\right)+\frac{e^{\theta }+e^{-\theta }}{\gamma -\left(e^{\theta }-e^{-\theta }\right)}\hfill \\ \hfill =& \left(e^{\theta }+e^{-\theta }\right)\left(\gamma -e^{\theta }\right)\left[1-\frac{e^{\theta }}{\gamma -\left(e^{\theta }-e^{-\theta }\right)}\right].\hfill \end{array}$$

Therefore, the stationary points of ${\tilde{\Phi }}_{\gamma }\left(\theta \right)$ solve the equations

$$\begin{array}{c}{e}^{\theta }=\gamma \hfill \end{array}$$

(69)

and

$$\begin{array}{c}{e}^{\theta }=\gamma -\left(e^{\theta }-e^{-\theta }\right).\hfill \end{array}$$

(70)

These equations yields two stationary points $e^{{\theta }_{\ast }}=\gamma$ and $e^{{\theta }_{\ast \ast }}=\frac{\gamma +\sqrt{8+{\gamma }^2}}{4}$. Correspondingly, the rate function ${\Phi }_{\gamma }\left(Y\right)$ has two stationary points

$$\begin{array}{c}{Y}_{\ast }=\gamma -{\gamma }^{-1}\mathrm{and}{Y}_{\ast \ast }=\frac{3\gamma -\sqrt{8+{\gamma }^2}}{4}=\frac{2\left(\gamma -\frac{1}{\gamma }\right)}{3+\sqrt{1+8/{\gamma }^2}}.\hfill \end{array}$$

It is evident that if $0<\gamma <1$ both stationary points $Y_{\ast }$ and $Y_{\ast \ast }$ are negative. One can easily check that in this case ${\Phi }_{\gamma }\left(Y\right)$ is monotonically increasing on the interval $Y>0$ and is positive on this interval.

If $\gamma >1$ then taking the second derivative in $\theta$ one can easily show that

$$\begin{array}{c}{\tilde{\Phi }}_{\gamma }^{″}\left({\theta }_{\ast }\right)=({\gamma }^2-1)({\gamma }^2+1)>0,{\tilde{\Phi }}_{\gamma }^{″}\left({\theta }_{\ast \ast }\right)=-\frac{\left(e^{2{\theta }_{\ast \ast }}-e^{-2{\theta }_{\ast \ast }}\right)\left(1+\gamma e^{{\theta }_{\ast \ast }}\right)}{{\left[\gamma -\left(e^{{\theta }_{\ast \ast }}-e^{-{\theta }_{\ast \ast }}\right)\right]}^2}<0,\hfill \end{array}$$

so that $Y_{\ast }$ is the point of local minimum of the rate function ${\Phi }_{\gamma }\left(Y\right)$, and $Y_{\ast \ast }$ is the point of local maximum. It is also easy to verify that the rate function ${\Phi }_{\gamma }\left(Y\right)$ vanishes in the limit $Y\to 0$ and also at $Y=Y_{\ast }$, staying positive at all other $Y>0$, so that that the point $Y=Y_{\ast }$ is the point of absolute minimum. Finally, to verify that the pre-exponential factor (34) vanishes at $Y=Y_{\ast \ast }$ it suffices to show that $r_{\ast }\left(Y_{\ast \ast }\right)(\gamma -Y_{\ast \ast })=1$. On noticing that

$$\begin{array}{c}{r}_{\ast }\left(Y\right)=\frac{\sqrt{Y^2+4}-Y}{2}=e^{-\theta }.\hfill \end{array}$$

this relation evidently follows from (68) and (70). □

Proof of Theorem 2.

In the scaling regime (42) the variable Y scales with $N^{-1/3}$. As $NY\gg 1$ in this case, the required asymptotic expressions for Laguerre polynomials can be read from (A8). It turns out that in order to calculate the density of imaginary parts to leading order in this regime, one has to retain the subleading term in the pre-exponential factor as specified in (A8). On substituting $Y=m{N}^{-1/3}$ in (A8) we obtain that with the required precision

$$\begin{array}{cc}{L}_{N-1}^{\left(1\right)}\left(-N^{1/3}{m}^2\right)\hfill & =\frac{e^{N{\mathcal{L}}_0\left(Y\right)}}{\sqrt{2\pi m^3}{\left(Y^2+4\right)}^{1/4}}\left(1-\frac{3}{16}\frac{1}{m{N}^{2/3}}\right)\hfill \end{array}$$

(71)

$$\begin{array}{cc}{L}_{N-2}^{\left(1\right)}\left(-N^{1/3}{m}^2\right)\hfill & =\frac{e^{N{\mathcal{L}}_0\left(Y\right)}}{\sqrt{2\pi m^3}{\left(Y^2+4\right)}^{1/4}}{r}_{\ast }{\left(Y\right)}^2\left(1-\frac{3}{16}\frac{1}{m{N}^{2/3}}\right)\hfill \end{array}$$

(72)

$$\begin{array}{cc}{L}_{N-1}^{\left(0\right)}\left(-N^{1/3}{m}^2\right)\hfill & =\frac{e^{N{\mathcal{L}}_0\left(Y\right)}}{N^{1/3}\sqrt{2\pi m^3}{\left(Y^2+4\right)}^{1/4}}\left(1+\frac{1}{16}\frac{1}{m{N}^{2/3}}\right),\hfill \end{array}$$

(73)

where ${\mathcal{L}}_0\left(Y\right)=Y{r}_{\ast }\left(Y\right)-2\ln r_{\ast }\left(Y\right)$ with $r_{\ast }$ (35) and $\ln r_{\ast }$ expanded in powers of $Y\ll 1$:

$$\begin{array}{c}{r}_{\ast }\left(Y\right)=1-\frac{Y}{2}+\frac{Y^2}{8}+O\left(Y^4\right),\ln r_{\ast }\left(Y\right)=-\frac{Y}{2}+\frac{Y^3}{48}+O\left(Y^4\right).\hfill \end{array}$$

(74)

It is easy to see that the overall exponential behaviour of the mean density (61) will still be given by (33) duly expanded:

$$\begin{array}{c}{\Phi }_{\gamma }\left(Y\right)=Y\left(\gamma +\frac{1}{\gamma }-2\right)-\frac{Y^2}{2}\left(1-\frac{1}{{\gamma }^2}\right)+\frac{Y^3}{3}\left(\frac{1}{{\gamma }^3}-\frac{1}{4}\right)+O\left(Y^4\right).\hfill \end{array}$$

(75)

Putting in here the scaling form $\gamma =1+\frac{\alpha }{N^{1/3}}$ and recalling $Y=\frac{m}{N^{1/3}}$ we find from (75), assuming that the parameters $\alpha \in \mathbb{R}$ and $m>0$ are fixed, that

$$\begin{array}{c}N{\Phi }_{1+\frac{\alpha }{N^{1/3}}}\left(\frac{m}{N^{1/3}}\right)=m{\alpha }^2-m^2\alpha +\frac{m^3}{4}=m{\left(\alpha -\frac{m}{2}\right)}^2:={\Phi }_{\alpha }\left(m\right).\hfill \end{array}$$

This verifies the exponent in (43). To find the pre-exponential terms we find it most convenient to use Equation (63). Substituting there (71)–(73) we first get

$$\begin{array}{cc}{F}_N\left(m\right)=\hfill & \frac{e^{N{\mathcal{L}}_0\left(Y\right)}}{\sqrt{2\pi m}{\left(Y^2+4\right)}^{1/4}}\left\{-\frac{2\gamma }{N^{1/3}}{r}_{\ast }\left(Y\right)\left(1+\frac{1}{16}\frac{1}{m{N}^{2/3}}\right)\right.\hfill \\ \hfill & +\frac{1}{m}\left(1-\frac{3}{16}\frac{1}{m{N}^{2/3}}\right)\left(\frac{3m}{N^{1/3}}+\frac{2\gamma }{N}-\frac{3m}{N^{4/3}}\right)\hfill \\ \hfill & \left.+\left[-\frac{2}{N^{1/3}}+\frac{1}{N^{1/3}}\left({\gamma }^2-2\gamma \frac{m}{N^{1/3}}+\frac{m^2}{N^{2/3}}\right)\right]r_{\ast }{\left(Y\right)}^2\left(1-\frac{3}{16}\frac{1}{m{N}^{2/3}}\right)\right\}.\hfill \end{array}$$

(76)

After rearranging and collecting the relevant terms in the above expression we arrive at

$$\begin{array}{c}{F}_N\left(m\right)=\frac{e^{N{\mathcal{L}}_0\left(Y\right)}}{\sqrt{2\pi m}{\left(Y^2+4\right)}^{1/4}}\left\{\frac{{(r_{\ast }\left(Y\right)\gamma -1)}^2+2(1-r_{\ast }{\left(Y\right)}^2)}{N^{1/3}}-\frac{2\gamma m{r}_{\ast }{\left(Y\right)}^2}{N^{2/3}}\right.\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\left.+\frac{1}{N}\left[-\frac{9}{16m}+m^2{r}_{\ast }{\left(Y\right)}^2-\frac{3}{16}\frac{({\gamma }^2-2)r_{\ast }{\left(Y\right)}^2}{m}+\frac{2\gamma }{m}-\frac{\gamma r_{\ast }\left(Y\right)}{8m}\right]\right\}.\end{array}$$

(77)

The expansion (74) together with $\gamma =1+\frac{\alpha }{N^{1/3}}$ give the relations

$$\begin{array}{cc}\frac{{(r_{\ast }\left(Y\right)\gamma -1)}^2+2(1-r_{\ast }{\left(Y\right)}^2)}{N^{1/3}}=\hfill & \frac{2m}{N^{2/3}}+\frac{1}{N}\left[{\left(\alpha -\frac{m}{2}\right)}^2-m^2\right]\hfill \end{array}$$

(78)

and

$$\begin{array}{cc}-\frac{2\gamma m{r}_{\ast }{\left(Y\right)}^2}{N^{2/3}}=\hfill & -\frac{2m}{N^{2/3}}-\frac{2m(\alpha -m)}{N},\hfill \end{array}$$

(79)

which are exact to the subleading order. We can now see that the leading order terms inside the curly brackets in (77) cancel. This also implies that at the leading order it is enough to replace the factor ${\left(Y^2+4\right)}^{1/4}$ in (77) with $\sqrt{2}$. Finally, adding the leading order contribution from

$$\begin{array}{c}\frac{1}{N}\left[-\frac{9}{16m}+m^2{r}_{\ast }{\left(Y\right)}^2-\frac{3}{16}\frac{({\gamma }^2-2)r_{\ast }{\left(Y\right)}^2}{m}+\frac{2\gamma }{m}-\frac{\gamma r_{\ast }\left(Y\right)}{8m}\right]=\frac{1}{N}\left(m^2+\frac{3}{2m}\right)\hfill \end{array}$$

to the $1/N$ terms in (78) and (79) results in

$$\begin{array}{c}{F}_N\left(m\right)=\frac{e^{N{\mathcal{L}}_0\left(Y\right)}}{2\sqrt{\pi m}}\frac{1}{N}\left[\frac{3}{2m}+{\left(\alpha -3\frac{m}{2}\right)}^2\right],\hfill \end{array}$$

(80)

thus verifying the pre-exponential factors in (43). □

Let us finally present the derivation of the marginal density of imaginary parts () pertinent to keeping the product $y=YN$ fixed as $N\to \infty$. This task is straightforwardly achieved by performing the limit $N\to \infty$ in (61) via substituting the corresponding asymptotics of Laguerre polynomials (A3) into the Formula (63) and using the identity $\frac{d}{dy}{I}_1\left(2y\right)=I_0\left(2y\right)-I_2\left(2y\right)$.

6. Proof of Theorem 3

Proof. 

We will use Equations (56)–(58) which express the mean density of eigenvalues ${\rho }_N(X,Y)$ in terms of the rescaled Hermite polynomials ${\pi }_k(X+iY)$ (59).

Using the integral representation in (57) it can be shown that in the scaling limit

$$\begin{array}{c}z=X+iY,X=\frac{q}{N^{1/3}},Y=\frac{m}{N^{1/3}}>0\hfill \end{array}$$

(81)

the rescaled Hermite polynomials ${\pi }_k\left(z\right)$ are given by the asymptotic equations

$$\begin{array}{c}{\pi }_k\left(z\right)\sim \sqrt{\frac{2\pi }{N(1+{\sigma }_{+}^2)}}{\left(-i{\sigma }_{+}\right)}^k{e}^{-\frac{N}{2}\left(1+iz{\sigma }_{+}+2\ln \left(-i{\sigma }_{+}\right)\right)},\hfill \end{array}$$

(82)

$$\begin{array}{c}{\pi }_k\left(\overline{z}\right)\sim \sqrt{\frac{2\pi }{N(1+{\sigma }_{-}^2)}}{\left(-i{\sigma }_{-}\right)}^k{e}^{-\frac{N}{2}\left(1+i\overline{z}{\sigma }_{-}+2\ln \left(-i{\sigma }_{-}\right)\right)},\hfill \end{array}$$

(83)

where we have introduced the notations

$$\begin{array}{c}{\sigma }_{+}=\frac{iz+\sqrt{4-z^2}}{2},{\sigma }_{-}=\frac{i\overline{z}-\sqrt{4-{\overline{z}}^2}}{2}.\hfill \end{array}$$

(84)

This implies for $J_N(X,Y)$ (58) that

$$\begin{array}{c}{J}_N(X,Y)\sim \frac{2\pi }{N}{e}^{-N}\frac{({\sigma }_{+}-{\sigma }_{-})}{\sqrt{(1+{\sigma }_{+}^2)(1+{\sigma }_{-}^2)}}{e}^{-N\left(\frac{i}{2}\left(z{\sigma }_{+}+\overline{z}{\sigma }_{-}\right)+\ln \left(-{\sigma }_{+}{\sigma }_{-}\right)\right)}\hfill \end{array}$$

(85)

$$\begin{array}{c}\times \left[\left(1-{\sigma }_{+}(\gamma -y)\right)\left(1+{\sigma }_{-}(\gamma -Y)\right)-{\sigma }_{+}{\sigma }_{-}\frac{\gamma -Y}{NY}\right].\hfill \end{array}$$

(86)

We are here interested in the limit of small $\left|z\right|$ (81), and, hence, can use the expansions

$$\begin{array}{c}{\sigma }_{+}=1+\frac{iz}{2}-\frac{z^2}{8}+\dots ,{\sigma }_{-}=-1+\frac{i\overline{z}}{2}+\frac{{\overline{z}}^2}{8}+\dots \hfill \end{array}$$

and, consequently,

$$\begin{array}{cc}1-{\sigma }_{+}(\gamma -Y)=\hfill & 1-\gamma +Y+\frac{Y}{2}\gamma -\frac{iX}{2}\gamma +{O\left(\right|z|}^2),\hfill \\ 1-{\sigma }_{-}(\gamma -Y)=\hfill & 1-\gamma +Y+\frac{Y}{2}\gamma +\frac{iX}{2}\gamma +{O\left(\right|z|}^2).\hfill \end{array}$$

Hence,

$$\begin{array}{c}\left(1-{\sigma }_{+}(\gamma -Y)\right)\left(1+{\sigma }_{-}(\gamma -Y)\right)={\left[1-\gamma +Y\left(1+\frac{\gamma }{2}\right)\right]}^2+\frac{X^2}{4}{\gamma }^2+O\left(X^2+Y^2\right).\hfill \end{array}$$

Setting here $X=\frac{q}{N^{1/3}}$, $Y=\frac{m}{N^{1/3}}$ and $\gamma =1+\frac{\alpha }{N^{1/3}}$ one obtains that to leading order in N

$$\begin{array}{c}\left(1-{\sigma }_{+}(\gamma -Y)\right)\left(1+{\sigma }_{-}(\gamma -Y)\right)=\frac{1}{N^{2/3}}\left[{\left(\frac{3}{2}m-\alpha \right)}^2+\frac{q^2}{4}\right].\hfill \end{array}$$

With the same precision we have

$$\begin{array}{c}-{\sigma }_{+}{\sigma }_{-}\frac{\gamma -Y}{NY}=\frac{1}{N^{2/3}}\frac{1}{m},\hfill \end{array}$$

and, consequently,

$$\begin{array}{c}\left[\left(1-{\sigma }_{+}(\gamma -Y)\right)\left(1+{\sigma }_{-}(\gamma -Y)\right)-{\sigma }_{+}{\sigma }_{-}\frac{\gamma -Y}{NY}\right]=\frac{1}{N^{2/3}}\left[\frac{1}{m}+{\left(\frac{3m}{2}-\alpha \right)}^2+\frac{q^2}{4}\right].\hfill \end{array}$$

(87)

On inspecting (56) and (85), one concludes that the overall exponential factor in (56) is given by $e^{-N{\tilde{\Phi }}_{\gamma }}$, where

$$\begin{array}{c}{\tilde{\Phi }}_{\gamma }=\frac{i}{2}(z{\sigma }_{+}+\overline{z}{\sigma }_{-})+\ln \left(-{\sigma }_{+}{\sigma }_{-}\right)+\frac{X^2}{2}+\gamma Y-Y^2-\ln \left(1-\frac{Y}{\gamma }\right).\hfill \end{array}$$

The leading order form of ${\tilde{\Phi }}_{\gamma }$ can be found by expanding in powers of X and Y, in a similar way as before:

$$\begin{array}{l}\frac{i}{2}(z{\sigma }_{+}+\overline{z}{\sigma }_{-})=-Y-\frac{X^2-Y^2}{2}+\frac{Y}{8}(3X^2-Y^2)+{O\left(\right|z|}^4),\\ \hspace{1em}\ln \left(-{\sigma }_{+}{\sigma }_{-}\right)=-Y+\frac{Y^3}{24}-\frac{Y{X}^2}{8}+{O\left(\right|z|}^4),\\ -\ln \left(1-\frac{Y}{\gamma }\right)=\frac{Y}{\gamma }+\frac{Y^2}{2{\gamma }^2}+\frac{Y^3}{3{\gamma }^3}+O\left(Y^4\right).\end{array}$$

Adding all contributions,

$$\begin{array}{c}{\tilde{\Phi }}_{\gamma }=Y\frac{{(\gamma -1)}^2}{\gamma }+\frac{Y^2}{2{\gamma }^2}(1-{\gamma }^2)+\frac{Y^3}{3}\left(\frac{1}{{\gamma }^3}-\frac{1}{4}\right)+\frac{Y{X}^2}{4}+{O\left(\right|z|}^4).\hfill \end{array}$$

Setting here $X=\frac{q}{N^{1/3}}$, $Y=\frac{m}{N^{1/3}}$ and $\gamma =1+\frac{\alpha }{N^{1/3}}$, one obtains that to leading order

$$\begin{array}{c}{\tilde{\Phi }}_{\gamma }=\frac{m}{N}\left({\left(\frac{m}{2}-\alpha \right)}^2+\frac{q^2}{4}\right).\hfill \end{array}$$

(88)

Combining (88) with (87), and trivially taking into account asymptotic expressions for the remaining multiplicative factors in (56) and (86), one arrives at (46). □