A Asymptotic analysis

Index Distribution of Cauchy Random Matrices


Using a Coulomb gas technique, we compute analytically the probability that a large Cauchy random matrix has positive eigenvalues, where is called the index of the ensemble. We show that this probability scales for large as , where is the Dyson index of the ensemble. The rate function is computed in terms of single integrals that are easily evaluated numerically and amenable to an asymptotic analysis. We find that the rate function, around its minimum at , has a quadratic behavior modulated by a logarithmic singularity. As a consequence, the variance of the index scales for large as , where is twice as large as the corresponding prefactor in the Gaussian and Wishart cases. The analytical results are checked by numerical simulations and against an exact finite formula which, for , can be derived using orthogonal polynomials.

1 Introduction

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, matrices with independent entries, and matrices with rotational invariance. While limited analytical insight is generally available for the former, rotationally invariant ensembles of matrices are generally characterized by a joint probability density function (jpdf) of the real eigenvalues of the form


where is a normalization constant, and is the Dyson index of the ensemble, identifying real symmetric, complex Hermitian and quaternion self-dual matrices, respectively. , the potential, is a function suitably growing at infinity that defines the model. For instance, for the Gaussian ensemble, or 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 , 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 popularized1 by Dyson [2], which will be reviewed in Section 2.

Perhaps the simplest and most natural of such questions concerns the random variable , defined as the number of eigenvalues contained in an interval on the real line. The average value can be clearly computed as , where 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 , where is the mean spacing in the bulk and is kept fixed as . Clearly and the variance grows logarithmically with ,


where the constant was computed by Dyson and Mehta [3]. Therefore the typical scale of fluctuations around the mean is , and the computation of higher moments [32, 33] reveals that on this scale, the random variable has a Gaussian distribution.

Another related observable, which on the contrary has surprisingly escaped a thorough investigation until very recently, is the index , defined as the number of eigenvalues exceeding a threshold . In particular, we will focus on the number of positive eigenvalues. Note that in this case, where is the full unbounded interval the previous result (2), valid on a small symmetric interval around the origin, is no longer applicable.

This random variable naturally arises in the study of the stability of a multidimensional potential landscape  [34]. For instance, in string theory may represent the potential associated with a moduli space [35]. As far as glassy systems are concerned, the point may instead represent a configuration of the system and the energy of that configuration [36]. Generally, for disordered systems may represent the free energy landscape. Typically this -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 -dimensional landscape depends on the real eigenvalues of the Hessian matrix 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), , is therefore a crucial indicator of how many directions departing from the stationary point are stable. Given a random potential , 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 , the statistics of was considered by Cavagna et al. [36]. Using replica methods with Grassmann variables, they found that around its mean value , the random variable has typical fluctuations of for large . Moreover, the distribution of these typical fluctuations is Gaussian, i.e.


where this form of the distribution is valid over a region of around the mean for large . This implies that variance grows logarithmically with , , with , for large . For atypically large fluctuations (), the Gaussian distribution (3) is no longer valid, and the large deviation tails were computed in [41, 42], this time for all using a Coulomb gas method. Setting , the probability density of the random variable was found to scale for large as2


where the rate function was computed exactly [41, 42] over the full range . It is independent of and has a minimum at (corresponding to the typical situation, where i.e. half of the eigenvalues are positive on average). The case (and similarly ) corresponds to the extreme situation where all eigenvalues are positive (negative). Therefore , i.e. the probability that the smallest eigenvalue 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 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 over a scale of one recovers a Gaussian distribution,




Note that for one recovers the result in [36]. The constant term was found [42] via the asymptotic expansion of a finite- 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 within the support of the spectral density leads to


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

Given the rather robust large 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 matrices which are real symmetric (), complex Hermitian () or quaternion self-dual (). The Cauchy ensemble is characterized by the following probability measure


where is the identity matrix . The definition (8) is evidently invariant under the similarity transformation , where is an orthogonal , unitary or symplectic matrix. The jpdf of the real eigenvalues can be then immediately written as


As in the Gaussian and Wishart cases, we can give an electrostatic interpretation of the jpdf (9), where the 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 placed at the point of coordinate 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 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 via the relation


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


We will also see below that the Cauchy ensemble, for , possesses the remarkable property of being exactly solvable for finite and , 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 has no edge (11). It was recently shown in [56] that, for large , the absence of an edge has indeed important consequences on the right large deviations of the top eigenvalue, , in this ensemble (see also Refs. [57, 58] for the study of for , 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 (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 of positive eigenvalues decays for large as


where the rate function , defined for , is independent of 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 ,


resulting in a Gaussian distribution of the index around the typical value , albeit with a variance growing with , i.e.




This result (15) obtained via a Coulomb gas approach, valid for any , is confirmed by an exact finite- formula, using orthogonal polynomials, for . An interesting feature of this result (15) is that the prefactor of the leading behavior of the variance is twice as large here as in the Gaussian (6) and Wishart (7) cases: . Given that the local bulk statistics in all these cases is described by the sine-kernel, one may argue that this factor of 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 (12) in terms of single integrals. Next, in section 4, we evaluate the asymptotic behavior of when is close to its average value , extracting the leading behavior of the variance of as a function of . Finally, in section 5, we derive an exact finite formula for the variance of in the case , showing a perfect agreement with the leading trend for large , before concluding in section 6. Some technical details have been confined to A and B.

2 Coulomb gas formulation and integral equation for the constrained density

Let the number of eigenvalues of a Cauchy random matrix larger than . The probability density of is by definition


where is the Heaviside step function and is defined in (9).

We start by writing the jpdf in exponential form,




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 as a continuous fluid, described by a normalized density . 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


where the action is given by


and are Lagrange multipliers, introduced to enforce the overall normalization of the density, and a fraction of eigenvalues exceeding . The action 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 to the right of and a fraction to the left of . This free energy scales as (and not just ) because of the strong all-to-all interactions among the particles. The energetic component of this free energy dominates over the entropic part , 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 , can be evaluated using a saddle point method for large . The constrained density of eigenvalues (depending parametrically on and ) is determined by the following variational condition


This Eq. (22), valid for inside the support of , can be differentiated once with respect to to give the following singular integral equation


where stands for Cauchy principal part. Solving (23) with the constraint 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 , a compact blob with two edges to the left of , and a non-compact blob to the right of , extending all the way to infinity and with a singularity for (for the situation is reversed, while only for the two blobs merge into the unconstrained single-support density (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 problem3, and from it to deduce .

We introduce the resolvent


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


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 ), where we expect to recover the density in Eq. (11). First, we multiply both sides in (23) (dropping the principal value) by and we integrate it over , obtaining


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


Hence the RHS of (26) can be written as


The second term of the sum in (28) (under the replacement ) is the original RHS of (26) with the sign changed. The first term of the sum is just . Hence, we have that the RHS of (26) is just equal to :


The left hand side (LHS) of (26) requires a little more algebraic manipulation to be expressed in terms of . 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 one has


Note that in this splitting the two integrals in the RHS are individually divergent due to the pole at , 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


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


Now, we use a different strategy, using the identity , to obtain


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


Therefore, the LHS of (26) can also be written as . We have then two distinct ways of writing the LHS, and we can use them both to cancel :


Eliminating and recalling that the RHS is equal to (29), we find the algebraic equation for


We proceed to determine the constants , and using the normalization condition of the density . From (24) (setting ), it implies that should asymptotically go as . Taking the limit in equation (37), we find equations for the coefficients


which implies and . Our algebraic equation, finally, becomes:


The two solutions read


Using (25), the density comes out as expected


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


where the constant will be determined by the new normalization condition of the rightmost blob . We repeat the same steps as for the unconstrained integral equation, multiplying (23) (without the principal value) by and integrating over . We get an extra term compared to Eq. (26), arising from the Lagrange multiplier


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


Imposing the condition that for , we get the two conditions and . Calling , for we get the equation


whose solutions are


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


or equivalently


where the edge points of the leftmost blob (for ) are determined as a function of 


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


For (unconstrained case), the solution is , and as expected. This means that we recover the unconstrained Cauchy case if we impose a number of positive eigenvalues exactly equal to , and this unconstrained case materializes when . As approaches , the edge moves towards the origin, while the other edge approaches infinity, until exactly at the critical point the two blobs merge back together. In Fig. 1, we show a plot of the density (48) for . We also show standard Monte Carlo simulations of particles distributed according to the Boltzmann weight under the Hamiltonian in (18), with the constraint that of them are forced to stay on the positive semi axis . We observe a nice agreement between our exact formula and the numerics.

Figure 1: Constrained density of eigenvalues for a matrix of size where 90% of the eigenvalues are forced to lie on the positive semi axis and the corresponding expected theoretical curve for [ from (51)].

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




where the additional term 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.

3 Computation of the rate function

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


Since by construction 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 ,


Multiplying this expression by and integrating over we obtain


Inserting (56) in (54) we obtain


We must now determine the constants and . To do so, we first consider, without loss of generality, the case , where the density has the shape as in Fig. 1. The left blob has a compact support . To determine the relation between and , we study the asymptotic behavior of equation (55) when . Since both and behave like in the large limit, we deduce that . Evaluating (55) at , we find the value of , and hence , in terms of :


We may finally write the action at the saddle point as


where (depending parametrically on ) is given in Eq. (48). The single integrals and 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 , corresponding to the “typical” situation, where half of the eigenvalues are positive and half are negative. In the extreme limit (which is the same as ), we get


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

Figure 2: Plot of the rate function .

In the next section, we perform a careful asymptotic analysis of the rate function around the minimum , which brings a quadratic behavior modulated by a logarithmic singularity. This is in turn responsible for the logarithmic growth of the variance of the index with (to leading order), and provides the correct prefactor.

4 Asymptotic analysis of in the vicinity of

We have already remarked that the typical situation is realized when the constant . Therefore it is necessary to expand the terms and in the action (59), as well as the integral constraint (51), for . It turns out that this calculation is highly nontrivial, as several cancellations occur in the leading and next-to-leading terms of each contribution (see A for details). Denoting (with ), the final result reads:


In A, we show that the relation between and is (to leading order for )


Inverting this relation, we find to leading order


implying from (61) that




From (52), we then have (for close to )


Restoring in the RHS of (66), we obtain the Gaussian behavior announced in the introduction [Eq. (14)]:




In the next section, we compare the asymptotic behavior of the variance with a closed expression valid for finite and , finding perfect agreement between the trends.

5 Exact derivation for the variance of for finite and

In B, we derive a general formula for the variance of any linear statistics at finite and . Specializing it to the index case, we deduce that


where is the kernel of the ensemble, built upon suitable orthogonal polynomials. It turns out that in the Cauchy case, the orthogonal polynomials satisfying