Index Distribution of Cauchy Random Matrices

Ensembles of matrices with random entries have been extensively studied since the seminal works of Wigner [1], Dyson [2] and Mehta [3]. However, many years before the official birth of Random Matrix Theory (RMT) in nuclear physics in the 1950s, statisticians had already introduced some of the RMT machinery in their studies on multivariate analysis [4, 5, 6]. Restricting our scope to matrices with real spectra, two main classes of ensembles are typically considered, $\iota)$ matrices with independent entries, and $\iota\iota)$ matrices with rotational invariance. While limited analytical insight is generally available for the former, rotationally invariant ensembles of $N\times N$ matrices are generally characterized by a joint probability density function (jpdf) of the $N$ real eigenvalues of the form

$$P_{\beta}({\bm{\lambda}})=\frac{1}{Z_{N,\beta}}\prod_{j<k}|\lambda_{j}-\lambda_{k}|^{\beta}\prod_{j=1}^{N}e^{-\beta V(\lambda_{j})}\;,$$

(1)

where $Z_{N,\beta}$ is a normalization constant, and $\beta=1,2,4$ is the Dyson index of the ensemble, identifying real symmetric, complex Hermitian and quaternion self-dual matrices, respectively. $V(x)$, the potential, is a function suitably growing at infinity that defines the model. For instance, $V(x)=x^{2}/2$ for the Gaussian ensemble, or $V(x)=x-\alpha\ln x$ for the Wishart-Laguerre ensemble.

Armed with (1), a wealth of statistical questions about eigenvalue distributions can be efficiently tackled for both finite and large $N$, such as the average density of states, gap distributions and statistics of extreme eigenvalues. While usually the focus is on typical fluctuations of such random variables, several interesting cases have been lately considered, where the interest lies instead on atypical (rare) fluctuations, far away from the average (see Ref. [7] for a review). To mention just a few, the large deviation probability of extreme eigenvalues of Gaussian [8, 9, 10, 11, 12] and Wishart random matrices [9, 13, 14], the number of stationary points of random Gaussian landscapes [15, 16, 17], the distribution of free energies in mean-field spin glass models [18, 19], conductance and shot noise distributions in chaotic mesoscopic cavities [13, 20], entanglement entropies of a pure random state of a bipartite quantum system [21, 22, 23, 24] and the mutual information in multiple input multiple output (MIMO) channels [25]. In addition, RMT has also proven an invaluable tool in understanding large deviation properties of various observables in the so called vicious walker (or nonintersecting Brownian motion) problem [26, 27, 28, 29, 30]. A powerful tool to deal with such instances is the Coulomb gas technique, originally popularized¹¹1E.P. Wigner had however used it already in 1957 [31] to compute the density of states of the Gaussian ensemble. by Dyson [2], which will be reviewed in Section 2.

Perhaps the simplest and most natural of such questions concerns the random variable $N_{\mathcal{I}}$, defined as the number of eigenvalues contained in an interval $\mathcal{I}$ on the real line. The average value $\langle N_{\mathcal{I}}\rangle$ can be clearly computed as $\langle N_{\mathcal{I}}\rangle=N\int_{\mathcal{I}}dx\rho(x)$, where $\rho(x)$ is the average density of eigenvalues of the ensemble. What about its fluctuations? Dyson  [2] studied the variance of the number of eigenvalues in the “bulk limit”, i.e. when $\mathcal{I}=[-\delta_{N}L/2,\delta_{N}L/2]$, where $\delta_{N}=\pi/\sqrt{2N}$ is the mean spacing in the bulk and $L$ is kept fixed as $N\to\infty$. Clearly $\langle N_{\mathcal{I}}\rangle=L$ and the variance grows logarithmically with $L$,

$$\langle(N_{\mathcal{I}}-L)^{2}\rangle\sim\frac{2}{\beta\pi^{2}}\ln(L)+B_{\beta}\;,$$

(2)

where the constant $B_{\beta}$ was computed by Dyson and Mehta [3]. Therefore the typical scale of fluctuations around the mean is $\sqrt{\ln(L)}$, and the computation of higher moments [32, 33] reveals that on this scale, the random variable $N_{\mathcal{I}}$ has a Gaussian distribution.

Another related observable, which on the contrary has surprisingly escaped a thorough investigation until very recently, is the index $N_{+}(\zeta)=\sum_{i=1}^{N}\theta(\lambda_{i}-\zeta)$, defined as the number of eigenvalues exceeding a threshold $\zeta$. In particular, we will focus on the number $N_{+}:=N_{+}(0)$ of positive eigenvalues. Note that in this case, where $\mathcal{I}$ is the full unbounded interval $[0,\infty)$ the previous result (2), valid on a small symmetric interval around the origin, is no longer applicable.

This random variable $N_{+}$ naturally arises in the study of the stability of a multidimensional potential landscape $V(x_{1},x_{2},\ldots,x_{N})$ [34]. For instance, in string theory $V$ may represent the potential associated with a moduli space [35]. As far as glassy systems are concerned, the point $\{x_{i}\}$ may instead represent a configuration of the system and $V(\{x_{i}\})$ the energy of that configuration [36]. Generally, for disordered systems $V(\{x_{i}\})$ may represent the free energy landscape. Typically this $N$-dimensional landscape displays a complex pattern of stationary points that play a relevant role both in statics and dynamics of such systems [34]. The stability of a stationary point of this $N$-dimensional landscape depends on the $N$ real eigenvalues of the $(N\times N)$ Hessian matrix $M_{i,j}=\left[\partial^{2}V/{\partial x_{i}\partial x_{j}}\right]$ which is symmetric. If all the eigenvalues are positive (negative), the stationary point is a local minimum (local maximum). If some, but not all, are positive then the stationary point is a saddle. The number of positive eigenvalues (the index), $0\leq N_{+}\leq N$, is therefore a crucial indicator of how many directions departing from the stationary point are stable. Given a random potential $V$, the entries of the Hessian matrix at a stationary point are usually correlated. However, often important insights can be obtained by discarding these correlations. In the simplest case, one may assume that the entries of the Hessian matrix are independent Gaussian variables. This then leads to the index problem for a real symmetric Gaussian matrix. This toy model, called the random Hessian model (RHM), has been studied extensively in the context of disordered systems [36], landscape based string theory [37], quantum cosmology [38], random supergravity theories [39] and multifield inflation theories [40].

For Gaussian matrices with $\beta=1$, the statistics of $N_{+}$ was considered by Cavagna et al. [36]. Using replica methods with Grassmann variables, they found that around its mean value $N/2$, the random variable $N_{+}$ has typical fluctuations of $\mathcal{O}(\sqrt{\ln N})$ for large $N$. Moreover, the distribution of these typical fluctuations is Gaussian, i.e.

$$\mathcal{P}^{(G)}_{\beta}(N_{+},N)\approx\exp\left[-\frac{\pi^{2}}{2\ln(N)}(N_{+}-N/2)^{2}\right]\;,$$

(3)

where this form of the distribution is valid over a region of $\mathcal{O}(\sqrt{\ln N})$ around the mean $\langle N_{+}\rangle=N/2$ for large $N$. This implies that variance grows logarithmically with $N$, $\mathrm{Var}(N_{+})\sim\sigma_{G}\ln N$, with $\sigma_{G}=1/\pi^{2}$, for large $N$. For atypically large fluctuations ($N_{+}\gg N/2$), the Gaussian distribution (3) is no longer valid, and the large deviation tails were computed in [41, 42], this time for all $\beta$ using a Coulomb gas method. Setting $\kappa=N_{+}/N$, the probability density of the random variable $N_{+}$ was found to scale for large $N$ as²²2Hereafter, $\approx$ stands for a logarithmic equivalence, $\lim_{N\to\infty}\frac{-\ln\mathcal{P}^{(G)}_{\beta}(N_{+}=\kappa N,N)}{\beta N^{2}}=\psi_{G}(\kappa)$.

$$\mathcal{P}^{(G)}_{\beta}(N_{+}=\kappa N,N)\approx\exp\left[-\beta N^{2}\psi_{G}(\kappa)\right]\;,$$

(4)

where the rate function $\psi_{G}(\kappa)$ was computed exactly [41, 42] over the full range $0\leq\kappa\leq 1$. It is independent of $\beta$ and has a minimum at $\kappa=1/2$ (corresponding to the typical situation, where $\langle N_{+}\rangle=N/2$ i.e. half of the eigenvalues are positive on average). The case $\kappa=1$ (and similarly $\kappa=0$) corresponds to the extreme situation where all eigenvalues are positive (negative). Therefore $\mathcal{P}^{(G)}_{\beta}(N_{+}=N,N)=\mathrm{Prob}[\lambda_{1}>0,\ldots,\lambda_{N}>0]=\mathrm{Prob}[\lambda_{\mathrm{min}}>0]$, i.e. the probability that the smallest eigenvalue $\lambda_{\mathrm{min}}$ is positive. Hence there is a natural connection between the index problem and the distribution of extreme eigenvalues, tackled in [8, 9, 13, 43, 44]. Note that the index problem in the complex plane (i.e. the statistics of the number of complex eigenvalues with modulus greater than a threshold) has also been recently considered [45, 46].

Expanding the rate function $\psi_{G}(\kappa)$ around the minimum, it was found that it does not display a simple quadratic behavior as one could have naively expected. Instead, the quadratic behavior is modulated by a logarithmic singularity, implying that in the close vicinity of $N_{+}=N/2$ over a scale of $\sqrt{\ln N}$ one recovers a Gaussian distribution,

$$\mathcal{P}^{(G)}_{\beta}\left(N_{+},N\right)\approx\exp\left[-\frac{\left(N_{+}-N/2\right)^{2}}{2\left(\mathrm{Var}(N_{+})\right)}\right]\mbox{ for }N_{+}\to N/2\;,$$

(5)

with

$$\mathrm{Var}(N_{+})\sim\sigma_{G}\ln N+C_{\beta}+\mathcal{O}(1/N)\;,\;\sigma_{G}=\frac{1}{\beta\pi^{2}}\;.$$

(6)

Note that for $\beta=1$ one recovers the result in [36]. The constant term $C_{2}$ was found [42] via the asymptotic expansion of a finite-$N$ variance formula conjectured by Prellberg and later rigorously established [47]. Interestingly, the same analysis performed on positive definite Wishart matrices [48] for a threshold at $\zeta$ within the support of the spectral density leads to

$$\mathrm{Var}(N_{+}(\zeta))\sim\sigma_{G}\ln N+\mathcal{O}(1)\;,$$

(7)

where the leading term is independent of $\zeta$ and exactly identical to the Gaussian case. Note that an explicit formula for the full probability of the index for finite $N$ is not available to date in either case.

Given the rather robust large $N$ behavior of the variance which holds both for Gaussian matrices (6) and Wishart matrices (7), we investigate here the index probability distribution of yet another ensemble of random matrices for which we expect a different behavior. We consider indeed the Cauchy ensemble of $N\times N$ matrices $\mathbf{H}$ which are real symmetric ($\beta=1$), complex Hermitian ($\beta=2$) or quaternion self-dual ($\beta=4$). The Cauchy ensemble is characterized by the following probability measure

$$P(\mathbf{H})\propto\left[\det\left({\bm{1}}_{N}+{\mathbf{H}}^{2}\right)\right]^{-\beta(N-1)/2-1}\;,$$

(8)

where ${\bm{1}}_{N}$ is the identity matrix $N\times N$. The definition (8) is evidently invariant under the similarity transformation $\mathbf{H}\to\mathbf{UHU}^{-1}$, where $\mathbf{U}$ is an orthogonal $(\beta=1)$, unitary $(\beta=2)$ or symplectic $(\beta=4)$ matrix. The jpdf of the $N$ real eigenvalues can be then immediately written as

$$P_{\beta}({\bm{\lambda}})\propto\prod_{j=1}^{N}\frac{1}{(1+\lambda_{j}^{2})^{\beta(N-1)/2+1}}\prod_{i<k}|\lambda_{i}-\lambda_{k}|^{\beta}\;.$$

(9)

As in the Gaussian and Wishart cases, we can give an electrostatic interpretation of the jpdf (9), where the $\lambda_{i}$s are positions of charged particles (with say positive unit charge) on the real line and repelling each other via the 2d-Coulomb (logarithmic) interaction. Here they feel an additional interaction with a single particle, with charge $-(N-1+2/\beta)$ placed at the point of coordinate $(0,1)$ in the 2d plane. A closely related ensemble occurs in the context of mesoscopic transport where it represents the scattering matrix of a quantum dot coupled to the outside world by non ideal leads containing $N$ scattering channels [49, 50]. It is also one of the typical examples where free probability theory efficiently applies in the context of random matrices models [51, 52]. Interestingly, the Cauchy ensemble (9) is related to the circular ensemble of RMT via the stereographic projection [53]. Indeed, if one defines the angles $\theta_{k}$ via the relation

$$e^{\mathrm{i}\theta_{k}}=\frac{1-\mathrm{i}\lambda_{k}}{1+\mathrm{i}\lambda_{k}}\;,$$

(10)

then the jpdf of the $\theta_{k}$’s is precisely the one of the eigenvalues in the $\beta$-circular ensemble. This implies in particular that local fluctuations of the eigenvalues in the Cauchy ensemble are described, for large $N$, by the sine-kernel [54]. This connection with the circular ensembles implies also that, in contrast to most other random matrix models, the finite-$N$ spectral density $\rho_{N}(\lambda)$ is independent of $N$, i.e. it coincides for all $N$ with its asymptotic expression $\rho^{\star}_{C}(\lambda)$. This density has fat tails extending over the full real axis, and its expression is given for all $\beta$ by [53, 55]

$$\rho^{\star}_{C}(\lambda)=\frac{1}{\pi}\frac{1}{1+\lambda^{2}}\;.$$

(11)

We will also see below that the Cauchy ensemble, for $\beta=2$, possesses the remarkable property of being exactly solvable for finite $N$ and $\beta=2$, as the suitable orthogonal polynomials can be determined in terms of Jacobi polynomials (see Section 5 for details). Hence, a major difference with Gaussian or Wishart matrices for which the mean spectral density has a finite support is that in the case of Cauchy matrices, the density $\rho^{\star}_{C}(\lambda)$ has no edge (11). It was recently shown in [56] that, for large $N$, the absence of an edge has indeed important consequences on the right large deviations of the top eigenvalue, $\lambda_{\max}=\max_{1\leq i\leq N}\lambda_{i}$, in this ensemble (see also Refs. [57, 58] for the study of $\lambda_{\max}$ for $\beta=2$, though the large deviations were not studied there).

The purpose of this paper is to study the full probability distribution of the index for the Cauchy ensemble, using a Coulomb gas technique. As a bonus, we also obtain the constrained spectral density of a Cauchy ensemble with a prescribed fraction of positive eigenvalues. In the limit $\kappa\to 0$ (where all eigenvalues are negative), we recover the large deviation law for the largest Cauchy eigenvalue derived in [56]. We show that the probability distribution that a Cauchy matrix has a fraction $\kappa=N_{+}/N$ of positive eigenvalues decays for large $N$ as

$$\mathcal{P}^{(C)}_{\beta}(N_{+}=\kappa N,N)\approx\exp\left[-\beta N^{2}\psi_{C}(\kappa)\right]\;,$$

(12)

where the rate function $\psi_{C}(\kappa)$, defined for $0\leq\kappa\leq 1$, is independent of $\beta$ and is calculated exactly (in terms of single integrals that cannot be further simplified) in (53) and (59). The rate function has the following behavior close to the minimum at $\kappa=1/2$,

$$\psi_{C}(\kappa=1/2+\delta)\sim-\frac{\pi^{2}}{2}\frac{\delta^{2}}{\ln|\delta|}\mbox{ for }\delta\to 0,$$

(13)

resulting in a Gaussian distribution of the index around the typical value $\langle N_{+}\rangle=N/2$, albeit with a variance growing with $\ln N$, i.e.

$$\mathcal{P}^{(C)}_{\beta}\left(N_{+},N\right)\approx\exp\left[-\frac{\left(N_{+}-N/2\right)^{2}}{2\left(\mathrm{Var}(N_{+})\right)}\right]\mbox{ for }N_{+}\to N/2\;,$$

(14)

where

$$\mathrm{Var}(N_{+})\sim\sigma_{C}\ln N+\mathcal{O}(1)\;,\;\sigma_{C}=\frac{2}{\beta\pi^{2}}=2\sigma_{G}\;.$$

(15)

This result (15) obtained via a Coulomb gas approach, valid for any $\beta$, is confirmed by an exact finite-$N$ formula, using orthogonal polynomials, for $\beta=2$. An interesting feature of this result (15) is that the prefactor of the leading behavior $\propto\ln N$ of the variance is twice as large here as in the Gaussian (6) and Wishart (7) cases: $\sigma_{C}=2\sigma_{G}$. Given that the local bulk statistics in all these cases is described by the sine-kernel, one may argue that this factor of $2$ is due to the presence of an edge in the density of eigenvalues for Gaussian and Wishart matrices, which is absent for Cauchy matrices. A thorough investigation of this issue is deferred to a separate publication [59].

The plan of the paper is as follows. In section 2, we introduce the Coulomb gas formulation of the index problem, in terms of the minimization of an action which leads to a singular integral equation for the constrained density of eigenvalues. This integral equation is solved, in section 2, using the resolvent method. In section 3 we present the evaluation of the action at the constrained density, from which we compute the rate function $\psi_{C}(\kappa)$ (12) in terms of single integrals. Next, in section 4, we evaluate the asymptotic behavior of $\mathcal{P}^{(C)}_{\beta}\left(N_{+},N\right)$ when $N_{+}$ is close to its average value ${N}/{2}$, extracting the leading behavior of the variance of $N_{+}$ as a function of $N$. Finally, in section 5, we derive an exact finite $N$ formula for the variance of $N_{+}$ in the case $\beta=2$, showing a perfect agreement with the leading trend for large $N$, before concluding in section 6. Some technical details have been confined to A and B.

Let $N_{+}(\zeta)$ the number of eigenvalues of a Cauchy random matrix larger than $\zeta$. The probability density of $N_{+}(\zeta)$ is by definition

$$\mathcal{P}_{\beta}^{(C)}(N_{+}(\zeta),N)=\int d^{N}\lambda P_{\beta}({\bm{\lambda}})\delta\left(N_{+}(\zeta)-\sum_{j=1}^{N}\theta(\lambda_{j}-\zeta)\right),$$

(16)

where $\theta$ is the Heaviside step function and $P_{\beta}({\bm{\lambda}})$ is defined in (9).

We start by writing the jpdf $P_{\beta}({\bm{\lambda}})$ in exponential form,

$$P_{\beta}({\bm{\lambda}})\propto e^{-\beta H[{\bm{\lambda}}]},$$

(17)

where:

$$H[{\bm{\lambda}}]=\left(\frac{N-1}{2}+\frac{1}{\beta}\right)\sum_{j=1}^{N}\ln(1+\lambda^{2}_{j})-\sum_{i>j}\ln|\lambda_{i}-\lambda_{j}|.$$

(18)

In this form, the jpdf (17) mimics the Gibbs-Boltzmann weight of a system of charged particles in equilibrium under competing interactions. Following Dyson [2] we can treat this system for large $N$ as a continuous fluid, described by a normalized density $\rho(\lambda)=(1/N)\sum_{i=1}^{N}\delta(\lambda-\lambda_{i})$. Consequently, the multiple integral in (16) can be converted into a functional integral in the space of normalizable densities. This procedure was first successfully employed to compute the large deviation of maximal eigenvalue of Gaussian matrices [43] and afterwards applied in several different contexts [7, 8, 9, 13, 42].

In the continuum limit, the multiple integral (16) becomes

$$\mathcal{P}_{\beta}^{(C)}(N_{+}(\zeta),N)\propto\int\mathcal{D}[\rho]\int dA_{1}\int dA_{2}\ e^{-\frac{\beta}{2}N^{2}S[\rho]},$$

(19)

where the action $S$ is given by

	$\displaystyle S[\rho]=$	$\displaystyle\operatorname*{\int_{-\infty}^{+\infty}\!\!\!}dx\rho(x)\ln(1+x^{2})-\int\operatorname*{\int_{-\infty}^{+\infty}\!\!\!}dxdx^{\prime}\rho(x)\rho(x^{\prime})\ln|x-x^{\prime}|$		(20)
		$\displaystyle+A_{1}\left(\operatorname*{\int_{-\infty}^{+\infty}\!\!\!}dx\rho(x)-1\right)+A_{2}\left(\int_{\zeta}^{+\infty}dx\rho(x)-\kappa\right)\;,$		(21)

and $A_{1},A_{2}$ are Lagrange multipliers, introduced to enforce the overall normalization of the density, and a fraction $\kappa$ of eigenvalues exceeding $\zeta$. The action $S$ has an evident physical meaning: it represents the free energy of the Coulomb fluid, whose particles are constrained to split over two regions of the real line, a fraction $\kappa$ to the right of $\zeta$ and a fraction $1-\kappa$ to the left of $\zeta$. This free energy scales as $N^{2}$ (and not just $N$) because of the strong all-to-all interactions among the particles. The energetic component $\sim\mathcal{O}(N^{2})$ of this free energy dominates over the entropic part $\sim\mathcal{O}(N)$, making it possible to use a saddle point method (see below). Note that the entropic term has been thoroughly studied and employed to define interpolating ensembles in [60, 61].

As mentioned earlier, the integral (19), where we neglected terms of $\mathcal{O}(N)$, can be evaluated using a saddle point method for large $N$. The constrained density of eigenvalues $\rho^{\star}(x)$ (depending parametrically on $\kappa$ and $\zeta$) is determined by the following variational condition

$$\frac{\delta S}{\delta\rho}=\ln(1+x^{2})-2\int_{-\infty}^{+\infty}dx^{\prime}\rho^{\star}(x^{\prime})\ln|x-x^{\prime}|+A_{1}+A_{2}\theta(x-\zeta)=0\;.$$

(22)

This Eq. (22), valid for $x$ inside the support of $\rho^{\star}(x)$, can be differentiated once with respect to $x$ to give the following singular integral equation

$$\frac{x}{1+x^{2}}+A_{2}\delta(x-\zeta)=\mathcal{P}\int_{-\infty}^{\infty}\frac{\rho^{\star}(y)}{x-y}dy\;,$$

(23)

where $\mathcal{P}$ stands for Cauchy principal part. Solving (23) with the constraint $\int_{\zeta}^{\infty}dx\rho^{\star}(x)=\kappa$ is the main technical challenge. The physical intuition, supported by numerical simulations, points towards a density supported on two disconnected intervals (see Fig. 1): for $\kappa>1/2$, a compact blob with two edges to the left of $\zeta$, and a non-compact blob to the right of $\zeta$, extending all the way to infinity and with a singularity for $x\to\zeta^{+}$ (for $\kappa<1/2$ the situation is reversed, while only for $\kappa=1/2$ the two blobs merge into the unconstrained single-support density $\rho^{\star}_{C}(\lambda)$ (11)). In such a situation, the standard inversion formula [62] for singular integral equation of the type in (23) cannot be applied, as it holds only for single support (“one-cut”) solutions. However, a more general method, which we now present, allows to compute the resolvent (or Green’s function) for this two-cuts problem³³3See [63] for a more sophisticated approach based on loop equation techniques., and from it to deduce $\rho^{\star}$.

We introduce the resolvent

$$G(z)=\int\frac{\rho^{\star}(x)}{z-x}dx\;,$$

(24)

for the Cauchy case. It is an analytic function in the complex plane outside the support of the density. From the resolvent, the density can be computed in the standard way as

$$\frac{1}{\pi}\lim_{\varepsilon\to 0}\text{Im }G(x+\mathrm{i}\varepsilon)=\rho^{\star}(x)\;,$$

(25)

where Im stands for the imaginary part.

Unconstrained case. As a warm-up exercise, we first derive the resolvent equation (using purely algebraic manipulations) for the unconstrained case (corresponding to (23) when $A_{2}=0$), where we expect to recover the density $\rho^{\star}_{C}(x)$ in Eq. (11). First, we multiply both sides in (23) (dropping the principal value) by ${\rho^{\star}(x)}/{(z-x)}$ and we integrate it over $x$, obtaining

$$\int\frac{x}{1+x^{2}}\frac{\rho^{\star}(x)}{z-x}dx=\iint\frac{\rho^{\star}(x)\rho^{\star}(y)}{(x-y)(z-x)}dxdy\;.$$

(26)

Our goal is to express both sides in terms of $G(z)$ and, by doing so, obtain an algebraic equation for $G(z)$. We start by the right hand side (RHS) where we use the simple identity

$$\frac{1}{(z-x)(x-y)}=\left(\frac{1}{z-x}+\frac{1}{x-y}\right)\frac{1}{z-y}\;.$$

(27)

Hence the RHS of (26) can be written as

	$\displaystyle\iint\frac{\rho^{\star}(x)\rho^{\star}(y)}{(x-y)(z-x)}dxdy=$	$\displaystyle\iint\left(\frac{1}{z-x}+\frac{1}{x-y}\right)\frac{\rho^{\star}(x)\rho^{\star}(y)}{z-y}dxdy$
	$\displaystyle=$	$\displaystyle\iint\frac{\rho^{\star}(x)\rho^{\star}(y)}{(z-x)(z-y)}dxdy+\iint\frac{\rho^{\star}(x)\rho^{\star}(y)}{(x-y)(z-y)}dxdy.$		(28)

The second term of the sum in (28) (under the replacement $x\leftrightarrow y$) is the original RHS of (26) with the sign changed. The first term of the sum is just $G(z)^{2}$. Hence, we have that the RHS of (26) is just equal to $(1/2)G(z)^{2}$:

$$\iint\frac{\rho^{\star}(x)\rho^{\star}(y)}{(x-y)(z-x)}dxdy=\frac{1}{2}G(z)^{2}\;.$$

(29)

The left hand side (LHS) of (26) requires a little more algebraic manipulation to be expressed in terms of $G(z)$. We manipulate this expression in two different ways and exploit the equality between the results to get rid of one integral. First, using the identity $x/(1+x^{2})=1/x-1/(x(1+x^{2}))$ one has

$\displaystyle\int\frac{x}{1+x^{2}}\frac{\rho^{\star}(x)}{z-x}dx=\int\frac{1}{x}\frac{\rho^{\star}(x)}{z-x}dx-\int\frac{1}{x(1+x^{2})}\frac{\rho^{\star}(x)}{z-x}dx\;.$

(30)

Note that in this splitting the two integrals in the RHS are individually divergent due to the pole at $x=0$, but the divergence cancels out between the two. We may regularize each individually divergent integral by adding a small imaginary part in the denominator that is removed at the end of the calculation. Using the relation (27), we may express the first term of the sum in (30) as

$\displaystyle\int\frac{1}{x}\frac{\rho^{\star}(x)}{z-x}dx=$

$\displaystyle\int\left(\frac{1}{x}+\frac{1}{z-x}\right)\frac{\rho^{\star}(x)}{z}dx=\frac{1}{z}\underbrace{\int\frac{\rho^{\star}(x)}{x}dx}_{a_{0}}+\frac{1}{z}\underbrace{\int\frac{\rho^{\star}(x)}{z-x}dx}_{G(z)}=\frac{a_{0}}{z}+\frac{G(z)}{z}.$

(31)

The second term of the sum in (30) will not be calculated for now, and will be called $-\alpha(z)$. Using this manipulation (31), we have:

$$\int\frac{x}{1+x^{2}}\frac{\rho^{\star}(x)}{z-x}dx=\frac{a_{0}}{z}+\frac{G(z)}{z}-\alpha(z).$$

(32)

Now, we use a different strategy, using the identity $x/(1+x^{2})=(x-z)/(1+x^{2})+z/(1+x^{2})$, to obtain

	$\displaystyle\int\frac{x}{1+x^{2}}\frac{\rho^{\star}(x)}{z-x}dx=$	$\displaystyle\int\frac{x+z-z}{1+x^{2}}\frac{\rho^{\star}(x)}{z-x}dx$
	$\displaystyle=$	$\displaystyle-\underbrace{\int\frac{\rho^{\star}(x)}{1+x^{2}}dx}_{a_{1}}+z\int\frac{1}{1+x^{2}}\frac{\rho^{\star}(x)}{z-x}dx.$		(33)

The first term in the sum (33) is a constant, which we call $a_{1}$. Now we proceed to manipulate the second term in (33) to obtain

	$\displaystyle z\int\frac{1}{1+x^{2}}\frac{\rho^{\star}(x)}{z-x}dx=$	$\displaystyle z\int\frac{x-z+z}{x(1+x^{2})}\frac{\rho^{\star}(x)}{z-x}dx$		(34)
	$\displaystyle=$	$\displaystyle-z\underbrace{\int\frac{\rho^{\star}(x)}{x(1+x^{2})}dx}_{a_{2}}+z^{2}\underbrace{\int\frac{1}{x(1+x^{2})}\frac{\rho^{\star}(x)}{z-x}dx}_{\alpha(z)}.$		(35)

Therefore, the LHS of (26) can also be written as $-a_{1}-za_{2}+z^{2}\alpha(z)$. We have then two distinct ways of writing the LHS, and we can use them both to cancel $\alpha(z)$:

$$\frac{a_{0}}{z}+\frac{G(z)}{z}-\alpha(z)=\text{(RHS)}=-a_{1}-za_{2}+z^{2}\alpha(z)\;.$$

(36)

Eliminating $\alpha(z)$ and recalling that the RHS is equal to $(1/2)G(z)^{2}$ (29), we find the algebraic equation for $G(z)$

$$-a_{1}+z(a_{0}-a_{2})+zG(z)=\frac{(1+z^{2})}{2}G(z)^{2}\;.$$

(37)

We proceed to determine the constants $a_{0}$, $a_{1}$ and $a_{2}$ using the normalization condition of the density $\rho^{\star}(x)$. From (24) (setting $|z|\to\infty$), it implies that $G(z)$ should asymptotically go as $G(z)\sim 1/z$. Taking the limit $|z|\to\infty$ in equation (37), we find equations for the coefficients

$$-a_{1}+z(a_{0}-a_{2})+z\left(\frac{1}{z}+\mathcal{O}(z^{-2})\right)=\frac{z^{2}}{2}\left(\frac{1}{z}+\mathcal{O}(z^{-2})\right)^{2}\;,$$

(38)

which implies $a_{0}=a_{2}$ and $a_{1}=\frac{1}{2}$. Our algebraic equation, finally, becomes:

$$(1+z^{2})G(z)^{2}-2zG(z)+1=0.$$

(39)

The two solutions read

$$G(z)=\frac{2z\pm\sqrt{4z^{2}-4(1+z^{2})}}{2(1+z^{2})}=\frac{z\pm\mathrm{i}}{1+z^{2}}.$$

(40)

Using (25), the density comes out as expected

$$\frac{1}{\pi}\lim_{\varepsilon\to 0}\text{Im }G(x+\mathrm{i}\varepsilon)=\frac{1}{\pi}\frac{1}{1+x^{2}}=\rho^{\star}_{C}(x)\;.$$

(41)

Constrained case. Now, we consider the full index problem, i.e. with an extra term in the potential as in (23),

$$\frac{x}{1+x^{2}}+A_{2}\delta(x-\zeta)=\mathcal{P}\int_{-\infty}^{\infty}\frac{\rho^{\star}(y)}{x-y}dy,$$

(42)

where the constant $A_{2}$ will be determined by the new normalization condition of the rightmost blob $\int_{\zeta}^{\infty}\rho^{\star}(x)dx=\kappa$. We repeat the same steps as for the unconstrained integral equation, multiplying (23) (without the principal value) by $\frac{\rho^{\star}(x)}{z-x}$ and integrating over $x$. We get an extra term compared to Eq. (26), arising from the Lagrange multiplier

$$\int\frac{x}{1+x^{2}}\frac{\rho^{\star}(x)}{z-x}dx=\iint\frac{\rho^{\star}(x)\rho^{\star}(y)}{(x-y)(z-y)}dxdy-\frac{A_{2}}{z-\zeta}.$$

(43)

We absorb this new term into the RHS and proceed to express, as before, all integrals in terms of $G(z)$. Our new algebraic equation will then be:

$$-a_{1}+z(a_{0}-a_{2})+zG(z)=(1+z^{2})\left(\frac{G(z)^{2}}{2}-\frac{A_{2}}{z-\zeta}\right).$$

(44)

Imposing the condition that $G(z)\sim 1/z$ for $|z|\to\infty$, we get the two conditions $a_{1}=1/2$ and $a_{0}-a_{2}+A_{2}=0$. Calling $A_{2}=B/2$, for $\zeta=0$ we get the equation

$$\frac{B}{2}z+G(z)z-\frac{1}{2}=\left(z^{2}+1\right)\left(\frac{B/2}{z}+\frac{G(z)^{2}}{2}\right),$$

(45)

whose solutions are

$$G(z)=\frac{z^{2}\pm\sqrt{-Bz^{3}-Bz-z^{2}}}{z^{3}+z}\;.$$

(46)

Using (25), it is then easy to derive that the constrained density is:

$$\rho^{\star}(x)=\frac{1}{\pi}\frac{\sqrt{B(x^{3}+x)+x^{2}}}{|x^{3}+x|},$$

(47)

or equivalently

$$\rho^{\star}(x)=\frac{1}{\pi(1+x^{2})}\sqrt{\frac{B(x+\ell_{1})(x+\ell_{2})}{x}}\;,$$

(48)

where the edge points of the leftmost blob (for $\kappa>1/2$) $-\ell_{1},-\ell_{2}$ are determined as a function of $B$

	$\displaystyle\ell_{1}$	$\displaystyle=\frac{1}{2B}(1+\sqrt{1-4B^{2}})\;,$		(49)
	$\displaystyle\ell_{2}$	$\displaystyle=\frac{1}{2B}(1-\sqrt{1-4B^{2}}).$		(50)

The functional form in Eq. (48) holds for $x$ belonging to any of the two supports. The constant $B$ is then determined as a function of $\kappa$ by

$$\int_{0}^{\infty}dx\frac{1}{\pi(1+x^{2})}\sqrt{\frac{B(x+\ell_{1})(x+\ell_{2})}{x}}=\kappa\;.$$

(51)

For $\kappa\to 1/2$ (unconstrained case), the solution is $B\to 0$, and $\rho^{\star}(x)\to\rho^{\star}_{C}(x)$ as expected. This means that we recover the unconstrained Cauchy case if we impose a number of positive eigenvalues exactly equal to $N/2$, and this unconstrained case materializes when $B\to 0$. As $\kappa$ approaches $1/2$, the $\ell_{2}$ edge moves towards the origin, while the other edge $\ell_{1}$ approaches infinity, until exactly at the critical point $\kappa=1/2$ the two blobs merge back together. In Fig. 1, we show a plot of the density (48) for $\kappa=0.9$. We also show standard Monte Carlo simulations of $N=300$ particles distributed according to the Boltzmann weight $\propto e^{-\beta H[{\bm{\lambda}}]}$ under the Hamiltonian $H[{\bm{\lambda}}]$ in (18), with the constraint that $90\%$ of them are forced to stay on the positive semi axis . We observe a nice agreement between our exact formula and the numerics.

Finally, from (19), we obtain the decay of the probability of the index for large $N$ as

$$\mathcal{P}_{\beta}^{(C)}(N_{+},N)\approx\exp\left[-\beta N^{2}\psi_{C}(N_{+}/N)\right]\;,$$

(52)

with

$$\psi_{C}(\kappa)=\frac{1}{2}\left(S[\rho^{\star}]-S[\rho^{\star}]\Big{|}_{\kappa=1/2}\right)\;,$$

(53)

where the additional term $S[\rho^{\star}]\Big{|}_{\kappa=1/2}$ comes from normalization.

In the next section, we simplify the action (21) at the saddle point and express it in terms of two single integrals which are hard to evaluate in closed form. The resulting rate function (53) can anyway be efficiently evaluated numerically with arbitrary precision.

The action (21) evaluated at the saddle point density (48) reads

	$\displaystyle S[\rho^{\star}]=$	$\displaystyle\operatorname*{\int_{-\infty}^{+\infty}\!\!\!}dx\rho^{\star}(x)\ln(1+x^{2})-\int\!\!\!\operatorname*{\int_{-\infty}^{+\infty}\!\!\!}dxdx^{\prime}\rho^{\star}(x)\rho^{\star}(x^{\prime})\ln|x-x^{\prime}|$
		$\displaystyle+A_{1}\left(\operatorname*{\int_{-\infty}^{+\infty}\!\!\!}dx\rho^{\star}(x)-1\right)+A_{2}\left(\int_{0}^{+\infty}dx\rho^{\star}(x)-\kappa\right).$		(54)

Since by construction $\rho^{\star}(x)$ satisfies the normalization conditions, the terms in the second line are automatically zero. We can now replace the double integral with single integrals. We use equation (22) for $\zeta=0$,

$$2\int_{-\infty}^{+\infty}dx^{\prime}\rho^{\star}(x^{\prime})\ln|x-x^{\prime}|=\ln(1+x^{2})+A_{1}+A_{2}\theta(x).$$

(55)

Multiplying this expression by $\rho^{\star}(x)$ and integrating over $x$ we obtain

$$\int\!\!\!\operatorname*{\int_{-\infty}^{+\infty}\!\!\!}dxdx^{\prime}\rho^{\star}(x)\rho^{\star}(x^{\prime})\ln|x-x^{\prime}|=\frac{1}{2}\left(\operatorname*{\int_{-\infty}^{+\infty}\!\!\!}dx\rho^{\star}(x)\ln(1+x^{2})+A_{1}+A_{2}\kappa\right).$$

(56)

Inserting (56) in (54) we obtain

$$S[\rho^{\star}]=\frac{1}{2}\operatorname*{\int_{-\infty}^{+\infty}\!\!\!}dx\rho^{\star}(x)\ln(1+x^{2})-\frac{A_{1}}{2}-\frac{A_{2}}{2}\kappa.$$

(57)

We must now determine the constants $A_{1}$ and $A_{2}$. To do so, we first consider, without loss of generality, the case $\kappa>\frac{1}{2}$, where the density has the shape as in Fig. 1. The left blob has a compact support $[-\ell_{1},-\ell_{2}]$. To determine the relation between $A_{1}$ and $A_{2}$, we study the asymptotic behavior of equation (55) when $x\to+\infty$. Since both $2\int_{-\infty}^{+\infty}dx^{\prime}\rho^{\star}(x^{\prime})\ln|x-x^{\prime}|$ and $\ln(1+x^{2})$ behave like $2\ln(x)$ in the large $x$ limit, we deduce that $A_{1}=-A_{2}$. Evaluating (55) at $x=-\ell_{1}$, we find the value of $A_{1}$, and hence $A_{2}$, in terms of $\ell_{1}$:

$$A_{1}=2\operatorname*{\int_{-\infty}^{+\infty}\!\!\!}dx\rho^{\star}(x)\ln|x+\ell_{1}|-\ln(1+\ell_{1}^{2}).$$

(58)

We may finally write the action at the saddle point as

$$S[\rho^{\star}]=\underbrace{\frac{1}{2}\operatorname*{\int_{-\infty}^{+\infty}\!\!\!}dx\rho^{\star}(x)\ln(1+x^{2})}_{I_{1}}-(1-\kappa)\underbrace{\operatorname*{\int_{-\infty}^{+\infty}\!\!\!}dx\rho^{\star}(x)\ln|x+\ell_{1}|}_{I_{2}}+\frac{(1-\kappa)}{2}\underbrace{\ln(1+\ell_{1}^{2})}_{I_{3}}.$$

(59)

where $\rho^{\star}(x)$ (depending parametrically on $\kappa$) is given in Eq. (48). The single integrals $I_{1}$ and $I_{2}$ do not seem expressible in closed form. However the rate function (53) can be plotted without difficulty (see Fig. 2). One can see that the rate function is symmetric, has a minimum at $\kappa=1/2$, corresponding to the “typical” situation, where half of the eigenvalues are positive and half are negative. In the extreme limit $\kappa=0$ (which is the same as $\kappa=1$), we get

$$\psi_{C}(\kappa=0)=\frac{1}{2}\left(S[\rho^{\star}]\Big{|}_{\kappa=0}-S[\rho^{\star}]\Big{|}_{\kappa=1/2}\right)=\frac{\ln 2}{4}\approx 0.173287...$$

(60)

in agreement with the large deviation law for the largest Cauchy eigenvalue [Eq. (13) in [56] for $w=0$].