Pseudo-Hermitian Random Matrix Models: General Formalism

Pseudo-hermitian matrices are matrices hermitian with respect to an indefinite metric. They can be thought of as the truncation of pseudo-hermitian operators, defined over some Krein space, together with the associated metric, to a finite dimensional subspace. As such, they can be used, in the usual spirit of random matrix theory, to model chaotic or disordered PT-symmetric quantum systems, or their gain-loss-balanced classical analogs, in the phase of broken PT-symmetry. The eigenvalues of pseudo-hermitian matrices are either real, or come in complex-conjugate pairs. In this paper we introduce a family of pseudo-hermitian random matrix models, depending parametrically on their metric. We apply the diagrammatic method to obtain its averaged resolvent and density of eigenvalues as explicit functions of the metric, in the limit of large matrix size $N$. As a concrete example, which is essentially an ensemble of elements of the non-compact unitary Lie algebra, we choose a particularly simple set of metrics, and compute the resulting resolvent and density of eigenvalues in closed form. The spectrum consists of a finite fraction of complex eigenvalues, which occupy uniformly two two-dimensional blobs, symmetric with respect to the real axis, as well as the complimentary fraction of real eigenvalues, condensed in a finite segment, with a known non-uniform density. The numbers of complex and real eigenvalues depend on the signature of the metric, that is, the numbers of its positive and negative eigenvalues. We have also carried thorough numerical analysis of the model for these particular metrics. Our numerical results converge rapidly towards the asymptotic analytical large-$N$ expressions.


I. Introduction
PT-symmetric quantum mechanics (PTQM) [1] and its broader applications (see [2] for recent reviews) have been at the focus of intensive and prolific research activity during the past twenty-five years or so. Generally speaking, the Hilbert space of a PTsymmetric quantum mechanical model is endowed with a non-trivial metric operator, with respect to which the hamiltonian is hermitian. Hamiltonians of PTQM models with proper positive metrics are sometimes referred to as quasi-hermitian [3,4]. They are diagonalizable, and their spectrum is essentially real, because they are similar to a conventionally hermitian hamiltonian. In contrast, hamiltonians of PTQM models with indefinite metrics are referred to as pseudo-hermitian [5,6], and their eigenvalues are either real or come in complex-conjugate pairs. Pseudo-hermitian hamiltonians describe systems with broken P T −symmetry [2]. For a mathematically precise summary of the nomenclature of quasi-hermiticity and pseudo-hermiticity see [7].
Quasi-hermitian (QH) matrices can be thought of as truncated quasi-hermitian linear operators into finite dimensional space. More precisely, we say that the N × N matrix φ (the "hamiltonian") is quasi-hermitian with respect to the N × N hermitian proper (non-negative) metric B, if φ and its adjoint fulfill the intertwining relation (1.1) The relation (1.1) simply means that φ is hermitian in a vector space endowed with a non-trivial metric B. (For B = 1 it reduces of course to ordinary hermiticity.) Following [8,9,10], we should make further distinction between QH matrices and strictly-quasi-hermitian (SQH) matrices. SQH matrices are hermitian with respect to positive definite (and therefore invertible) metrics B, in contrast to merely QH matrices, whose associated non-negative metrics may be non-invertible.
SQH random matrices arise, for example, in studying the spectrum of stable small oscillations of large mechanical systems with large connectivity [8,9]. The low-frequency part of the spectrum of such systems was found in [8,9] to be universal, with the phonon density of states tending to a non-zero constant at zero frequency. In other words, phonons in such systems have universal spectral dimension d s = 1.
An interesting SQH random matrix model was introduced in [11]. These authors fixed a metric B, and took the hamiltonian φ as random, with the aim of studying numerically the dependence of the average density of eigenvalues and level spacing statistics on the metric. Yet another interesting example of an SQH random matrix model, akin to the Dicke model of superradiance, was provided by [12], in which a numerical study of the level spacing distribution was carried out.
The difference of the two sides of (1.1) is an anti-hermitian matrix. Thus, given a generic invertible hermitian metric B, (1.1) amounts to a system of N 2 independent real equations for the 2N 2 real parameters of φ, resulting in N 2 linearly independent solutions. Indeed, as was pointed out in [13] (and in the physics literature in [11]), the general solution of (1.1) for φ is φ = AB, (1.2) where A = A † is a hermitian matrix. Therefore, given B, there are N 2 linearly independent matrices which are SQH with respect to B. An important corollary of (1.2) is that for a positive definite metric, φ is similar to a hermitian matrix, with the similarity matrix being essentially the (invertible) square root of B. (The latter is well defined and hermitian, because B is positive definite.) Consequently, SQH matrices are diagonalizable and have purely real eigenvalue spectrum. This can be seen easily when we write Upon truncation to finite vector spaces, pseudo-hermitian operators beget pseudohermitian (PH) matrices. Such matrices are still defined by the intertwining relation (1.1), but now the metric B has indefinite signature.
The dimensionality of the space of solutions of (1.1) for φ is clearly independent of the signature of the metric B, as long as it is invertible. Thus, since the number of independent solutions of the form (1.2) saturates this dimensionality, it is the general solution of (1.1) regardless of the signature of B. Consequently, any matrix φ which is PH with respect to an indefinite invertible metric B can be written in the form (1.2), with A being a hermitian matrix.
A PH matrix, in contrast to an SQH matrix, need not 2 be similar to a hermitian matrix. In this case, since according to (1.1) φ is similar to its hermitian adjoint, it follows that the characteristic polynomial of φ has real coefficients: (det(z − φ)) * = det(z * − φ).
Thus, the eigenvalues of φ are either real, or come in complex-conjugate pairs. See [14] for a recent discussion of (real asymmetric) PH random matrices. For the sake of completeness, let us briefly discuss the case of a singular metric B. Thus, assume rank(B) = p < N . Then, by transforming to a basis in which B is diagonal, it is easy to see that (1.1) amounts to a system of N 2 − (N − p) 2 independent real equations for the 2N 2 real parameters of φ, resulting in N 2 + (N − p) 2 linearly independent solutions. In contrast, it is easy to check that solutions of the form (1.2) comprise in this case only a subspace of real dimension N 2 −(N −p) 2 . The complementary subspace of solutions, which cannot be expressed in the form (1.2), has therefore real dimension 2(N − p) 2 , namely, the dimension of an (N − p) × (N − p) complex matrix.
This means that not all complex eigenvalues of φ in this case would come in complex conjugate pairs. We shall not pursue here the study of such matrices any further, since in this paper we focus exclusively on invertible metrics B.
We should also mention in passing that some special PH matrices may be nondiagonalizable, as the reader can easily check by working out all the matrices which are PH with respect to the metric B = σ 3 = diag (1, −1). In this case, matrices with degenerate eigenvalues and whose diagonal elements are not equal, admit only a Jordan form and are not diagonalizable. Such non-diagonalizable matrices clearly form a set of zero measure in the space of all PH matrices with respect to a given indefinite metric, and need not concern us in our statistical analysis of the random matrix models we study in this paper.
In this paper we introduce and study a family of pseudo-hermitian random matrix models, which depend parametrically on a fixed (deterministic), invertible, indefinite metric B. The hermitian matrix A = φB −1 in (1.2) is uniquely determined, and is the source of randomness. Thus, the probability ensemble of the PH matrices φ is induced from that of the hermitian matrices A. We are free to choose the latter at will, and a natural choice is to draw A from the Gaussian Unitary Ensemble (GUE) where m > 0 is a parameter and Z N is the normalization constant. Thus, (1.4) induces a PDF on the 2N 2 real parameters of φ, 5) whereZ N is another normalization constant. Therefore, an N × N complex matrix φ is drawn with the Gaussian weight in (1.5), and then the remaining factor δ φ † B − Bφ , which is supported over the N 2 independent real and imaginary parts of the antihermitian matrix φ † B − Bφ, filters in those matrices which are PH with respect to the metric B. This is the matrix model we focus on, and our main goal in this paper is to determine explicitly the ensemble-averaged density of the eigenvalues of φ on the real axis and in the complex eigenvalues in the limit N → ∞.
We have also studied numerically the spectral statistics of φ for ensembles of A other than (1.4). Our results [15] indicate that many attributes of the spectral statistics of φ are universal, and depend, perhaps not surprisingly, only on the variance of fluctuations of A.
We have recently introduced and previewed the model (1.5) in [10], and the purpose of the present paper is to provide detailed exposition of our analysis. We derive explicit equations for determining the averaged resolvent and density of eigenvalues as explicit functions of the metric, in the limit of large matrix size N , by means of the diagrammatic 4 method. As a concrete example, we set the metric to be diagonal with entries +1 or −1, and compute the resulting resolvent and density of eigenvalues in closed form.
The spectrum in this case consists of a finite fraction of complex eigenvalues, which occupy uniformly two two-dimensional blobs of known shape, symmetric with respect to the real axis, as well as the complimentary fraction of real eigenvalues, condensed in a finite segment, with a known non-uniform density. The numbers of complex and real eigenvalues depend on the signature of the metric, that is, the numbers of its positive and negative eigenvalues, both of which are assumed to be finite fractions of N . The average spectrum thus obtained for a particular such metric, is essentially the spectrum of a randomly chosen generator of the non-compact unitary group with the same signature. Rotating this spectrum by ninety degrees in the complex plane gives us the spectrum of a randomly chosen generator of the corresponding non-compact orthogonal group.
We have also carried thorough numerical analysis of the model for these particular metrics. Our numerical results converge rapidly towards the asymptotic analytical large-N expressions. This is perhaps a good place to mention that products of the form (1.2), with A drawn from the GUE and a deterministic hermitian matrix B, appear also in random matrix theory with an external source [16,17], in studying the eigenvalue statistics 3 of the hermitian matrix H = A + B. This is in contrast to the spectral statistics of the non-hermitian (or more precisely, PH) random matrix φ = AB studied in the present paper.
The rest of the paper is organized as follows: In Section II we apply the method of hermitization and the diagrammatic method to deriving the self-consistent gap equations in the large-N , planar limit. We then analyze the gap equations in Section III and determine from them the phase structure of our model in the complex plane, which consists of the so-called non-holomorphic phase in a two-dimensional region D ⊂ C in which the eigenvalues of our PH matrix φ condense in the large-N limit, and a holomorphic phase, in the complementary domain of the complex plane. In particular, we derive explicit expressions for the Green's function of φ in the holomorphic phase (see Eqs. (3.25) and (3.28)), and in the non-holomorphic phase (see Eq. (3.38)). The discontinuity of this function in the holomorphic phase across the real axis determines the density of real eigenvalues of φ, while the density of complex eigenvalues of φ can be determined from that Green's function in the non-holomorphic phase by applying Gauss' law, as explained in (3.3). Section III is very detailed and contains many analytical results and interrelations between various quantities for generic metrics B. (See, e.g., (3.57) and (3.71).) Finally, in Section IV we apply our general formalism to the concrete case of the metric B = diag(1, . . . , 1, −1, . . . , −1) with k ones and N − k minus ones, as described above, and provide explicit analytical and numerical results for the density of eigenvalues on the real axis and in the complex plane (see Eqs. (4.10) and (4.18), (4.21) respectively). We then discuss the relevance of these results to the non-compact Lie algebras su(k, N − k) and so(k, N − k). Many technical and mathematical details are relegated to the four appendices at the end of this paper.

II. Hermitization and Large-N Analysis
In practice, it is more convenient to work with the representation φ = AB of the pseudo-hermitian matrix in (1.2), rather than with the singular PDF (1.5). The averaged density of eigenvalues of φ in the complex-w plane can be calculated from the Green's function (resolvent) as we explain in Section III. Here, and from here on, angular brackets denote averaging of A over the GUE ensemble (1.4).
Averaging over A becomes simpler if one can avoid the product of matrices AB by the method introduced in [18], following which we introduce the 2N × 2N block matrix which puts A and B in two different blocks. The resolvent of (2.2) is readily computed as in which the upper and lower diagonal blocks are essentially the resolvents of φ = AB and φ † = BA, evaluated at w = z 2 . Thus, we arrive at the Green's function where in the last step we have used equality of the traces of the two diagonal blocks of (2.3) due to isospectrality of φ = AB and φ † = BA. Therefore, we can deduce the desired Green's function (2.1) from (2.4), which is more amenable to diagrammatic expansion, in a straightforward manner.
Since H is typically non-hermitian, it might have complex eigenvalues. In Appendix A.1 we discuss the symmetries the spectrum of H in the complex plane. In the large-N limit, these eigenvalues can be dense in a two-dimensional domain D ⊂ C, rendering G 6 not an analytic function 4 of z, depending both on z and z * . As such, G, or for that matter, any of the averaged blocks of (2.3), cannot be determined everywhere in the complex plane C from the moments of H, which amounts to expanding them in inverse powers of z. The latter so-called perturbative Born series would converge only in the analyticity domains of (2.4), which excludes D. This is in contrast with resolvents of hermitian random matrices, which are always analytic functions of z everywhere off the real axis, where they can be determined, in principle, from their perturbative Born series.
As is well known, this difficulty in computing perturbatively averaged resolvents of non-hermitian matrices by resumming their Born series, is overcome by employing the Method of Hermitization [19,20,21,22]. According to this method we can reduce the difficult problem of averaging (2.3) over the random matrix A to the more familiar problem of averaging the resolventĜ of the 4N × 4N hermitian matrix .
(2.6) Therefore, off the real axis in the complex η-plane, the averaged resolvent Ĝ is an analytic function of the spectral parameter η. This will allow us to compute Ĝ by expanding it diagrammatically in powers 1/η, and then resum the series. This series would then converge to Ĝ everywhere off the real axis in the complex-η plane.
After resumming the series, we could let η → 0. The lower left block of (2.6) would then clearly converge to the desired average 1 z−H of (2.3), while the upper right block of (2.6) would converge to its adjoint. In contrast, the diagonal blocks would suffer a jump discontinuity in the limit, depending on whether one approaches the real η axis from above or from below. More specifically, if we set η = is, s ∈ R, then as s → 0, Thus, the diagonal blocks of (2.7) converge in the limit to two isospectral anti-hermitian matrices, while the limiting off-diagonal blocks are hermitian conjugate of each other. More detailed analysis of these diagonal blocks (before averaging over A) is given in Appendix A.2.
In our subsequent analysis we will have to refer to the individual N × N blocks of Let us therefore denote these blocks aŝ G αβ , α, β := 1, 2, 3, 4 (2.9) Matrix elements within a given block will be indexed by Latin indices as Ĝ αβ ij , i, j := 1, . . . , N (2.10) Following 't Hooft, we will use Feynman diagrams with double-lined matrix ('gluon') propagators [23] in the large-N planar limit [24]. To this end we expand the resolvent Ĝ (η; z, z * ) in powers of the bare 'quark' propagator and rearrange the perturbative expansion in terms of the self-energy (the sum over onequark-irreducible graphs)Σ =Σ , (2.12) in the following way: This rearrangement of the perturbation series (valid to all orders in the large-N expansion) is the diagrammatic representation of the relation Ĝ = 1 between Ĝ andΣ, which is equivalent to the definition ofΣ as the sum over one-quarkirreducible graphs. As was explained above, this rearrangement is allowed because Ĝ is an analytic function of η, and therefore the corresponding Born series is convergent. This would not have been possible before the procedure of Hermitization.
We can express [19,25,26]Σ in terms of the connected cumulants of the distribution of A and the full propagator Ĝ . Since A is drawn from the GUE ensemble (1.4), there is only one connected cumulant, namely, the quadratic one (2.14) In the limit N → ∞ only planar diagrams survive, and we thus arrive at The block traces (2.16) are scalars, leading to the 4×4 matrix factor on the right-hand side of Eq. (2.15). Each block ofΣ is proportional to the N -dimensional unit matrix 1 N after averaging over the unitary-invariant ensemble (1.2). The zero-rich block texture of Σ is very helpful practically, and it arises due to the decoupling of the random blocks from the deterministic blocks, which was the motivation for introducing (2.2) in the first place. Thus,Σ is determined exclusively by the four block traces 44, 41, 14 and 11.
Substitution of (2.15) back in (2.13) leads to the so-called gap equation which determines all the blocks Ĝ αβ (η; z, z * ) explicitly. (Note that in the matrix to be inverted on the right hand side of (2.17), we have suppressed explicit factors of the unit matrix 1 N in the appropriate blocks.) 9 This substitution of (2.15) back in (2.13) amounts to resummation of the perturbation series. Thus, at this stage of solving the gap equation, we can safely set η = is → 0 everywhere. Detailed subsequent calculations of the gap equation (2.17) at η = 0 are given in Appendix A.4.
As we can see from (2.7), in the limit s → 0, As can be seen from (A.45) in Appendix A.4, which we copy here a subset of the gap equations determines the four block traces 11, 14, 41 and 44 selfconsistently, which in turn determine all the remaining blocks of Ĝ .
Here and in what follows, in order to avoid cluttering of our equations, we have defined (with some minimal abuse of matricial notations) As was discussed in Appendix A.4, by equating block traces on both sides of (2.24) and with the definition (2.22) in mind, we obtain the set of self-consistent equations for 11, 14, 41 and 44, or equivalently, the quantities a, b and c, as The first two equations in (2.26) suggest that a(z, z * ) and c(z, z * ) be equal, a proposition we could not prove for the corresponding diagonal block traces 11(is; z, z * ) and 44(is; z, z * ) in Appendix A.2, before taking averages. Let us now prove that An important observation is that the two distinct matrices under the trace in these expressions are positive definite, due to (2.23). Thus, each of the traces appearing in (2.28) is positive. Compute now the two sums of traces appearing in (A.18) is a positive matrix, rendering tr P > 0. According to (A.18), the two sums on the LHS of both equations in (2.29) are equal, thus proving (2.27). Combining (2.27) and (A.14), which implies that a(z, z * ) = c(z * , z), we conclude also that a(z, z * ) = a(z * , z). (2.31) As was mentioned following (A.18) and (A.29) in Appendix A.2, the sums of diagonal block traces on the LHS of (2.29) serve as order parameters indicating the location of the two-dimensional support D ⊂ C of the density of eigenvalues of H in the large-N limit [19]. After averaging over A, this role is clearly played by a(z, z * ) = c(z, z * ), which vanish only outside the domain of eigenvalues D.
Finally, after solving the self-consistent equations for a = c and b explicitly, we can then substitute them into the block (cf (A.49)) and read off from it the resolvent of φ = AB. Thus, we obtain the desired Green's function (2.1) as

III. Solution of the Gap Equations and the Density of Eigenvalues in the Complex Plane and on the Real Axis
Before delving into the technical details of computing the density of eigenvalues of the PH matrix φ from (2.33), let us get oriented by recalling the useful analogy between the formulation of two-dimensional electrostatics in the complex plane and the problem of determining the density of eigenvalues of non-hermitian random matrices in the complex plane [27].

III.1. The Two-Dimensional Electrostatic Analog of Spectra of Pseudo-Hermitian Matrices
Consider an N × N non-hermitian matrix X with eigenvalues λ 1 , . . . , λ N , taken from some probability ensemble. In the large-N limit, these eigenvalues typically become dense on average 7 in a two-dimensional domain D in the complex plane, with continuous two-dimensional eigenvalue distribution Thus, ρ (2) (x, y) is supported in the two-dimensional domain D in the complex plane 8 of w = x + iy, and the associated Green's function (analogous to (2.1)) is which one readily recognizes as the planar electric field generated by the charge density ρ (2) (x, y). By applying Gauss' law to the electric field (3.2) we thus recover the charge density, namely, the eigenvalue density indicating that G 2 (w, w * ) is not a holomorphic function of w in the spectral domain D.
Its domain of analyticity in w is rather the complementary domain D c = C/D of the complex plane, which does not contain any eigenvalues. It also follows from In contrast, the eigenvalues of hermitian matrices X = X † , with real eigenvalues λ i , would typically condense along a one-dimensional segment (or several segments) σ of the real axis R, giving rise to a continuous one-dimensional density supported along σ, with Green's function namely, the planar electric field generated by a charged one-dimensional wire σ placed along the real axis. G 1 (w) is manifestly a holomorphic function of w everywhere off the real axis, and it suffers a discontinuous jump across the real axis. The singular onedimensional charge density ρ (1) (x) which causes this jump, can be recovered from this discontinuity by means 9 of the well-known relation In contrast, as was mentioned above, G 2 (w, w * ) is continuous everywhere, and in particular, across the real axis.
Here we come to the crux of our preliminary discussion: PH matrices are sort of a hybridization of hermitian and non-hermitian matrices, in the sense that their eigenvalues are either real, or come in complex-conjugate pairs. A large N × N PH matrix φ would typically have macroscopic amounts, that is, finite fractions of N , of both real and pairs of complex conjugate eigenvalues. In the large-N limit, the complex pairs of eigenvalues would condense in a two-dimensional domain D, which is symmetric with respect to the real axis R, giving rise to a continuous two-dimensional distribution ρ (2) (x, y), while the real eigenvalues would condense along a segment σ of the real axis, giving rise to a continuous one-dimensional distribution ρ (1) (x). These eigenvalue densities are normalized to the corresponding fractions of N of eigenvalues they account for, that is, with 0 ≤ ν ≤ 1. The combination of these two-dimensional and one-dimensional distributions then gives rise to the Green's function (2.1) of our PH matrix φ, namely, Due to the reflection symmetry of the domain D with respect to the real axis, originating from isospectrality of φ and φ † , we can see that Moreover, it follows from (3.8) that G PH (w, w * ) is continuous everywhere, except in an infinitesimal sliver along the real axis, across which it jumps discontinuously as in (3.6), from which, in combination with (3.9), we can determine the one-dimensional density component of (3.8) according to 10 (3.10) is an analytic function of w in the complementary domain D c (and off the real axis), and evidently, has asymptotic behavior as w → ∞. (This asymptotic regime obviously always belongs in D c if the spectral domain D ∪ σ of φ is compact.) Finally, by applying Gauss' law to (3.8), we obtain the density of eigenvalues as with ρ (2) (x, y) supported throughout D and ρ (1) (x) supported along σ. Thus, ρ (2) (x, y) can be extracted from (3.12) simply by subtracting the singular piece proportional to δ(y). In the simpler case where D ∩ σ = ∅, we can avoid this subtraction by restricting w to D, and determine ρ (1) (x) directly from (3.10). We shall discuss such a case as an example in Section IV.

III.2. Solution of the Gap Equations and Phase Structure
Let us now substitute a = c from (2.27) and w = z 2 in (2.26). Thus, we are instructed to solve the self-consistent gap equations (supplemented by the complex conjugated last equation) for the purely imaginary quantity a(w, w * ), that is, and complex b(w, w * ).
In this section we shall analyze these equations assuming a generic invertible metric B. As a concrete example, in Section IV we shall outline the explicit full solution of these equations for a particular metric.
These equations are algebraic polynomial equations for the unknown functions, and therefore have multiple roots, which typically become degenerate at branch point singularities in the complex-w plane 11 . These roots should be then sewn together into unique single-valued continuous solutions for a and b, defined globally throughout the properly cut complex w-plane. Note that a and b may certainly depend on both w and w * , and are therefore not globally holomorphic functions.
At a given point w (away from any branch point) and in some small neighborhood around it, the "physical" solution of (3.13) should be clearly unique and continuous.
From the first equation in (3.13) we see that at a given point w, the solution for a(w, w * ) is either a(w, w * ) = 0 or a(w, w * ) = 0. On the basis of continuity, whichever of these two possibilities that holds at w, should also be valid in a neighborhood of that point. In fact, following the discussion below (2.31) in the previous section, a(w, w * ) = 0 throughout the two-dimensional spectral domain D of our averaged PH matrix φ = AB, and vanishes in the complementary domain 12 D c = C/D. Thus, a(w, w * ) serves as a local order parameter indicating whether w ∈ D or not.
For reasons that should become clear below, we shall refer to the solution with a(w, w * ) = 0 throughout D as the non-holomorphic phase, and to the complementary solution with a(w, w * ) = 0 throughout D c as the holomorphic phase.
We have established the role of a(w, w * ) as the order parameter telling apart the two phases. The other function b(w, w * ) also has a physical interpretation, as sort of generalized self-energy. We argue as follows: Let us recall the resolvent of φ from (2.32). Thus, we can rewrite the second equation in (3.13) as In the simplest case B = 1 N (discussed in Section III.3 below), this equation reads which we readily recognize from the diagrammatic equation (2.15), considered for the plain hermitian GUE matrices A of (1.4), as the relation between the self-energy and the Green's function. Eq. (3.16) merely generalizes this relation to PH matrices, with arbitrary metric B.
Outside the spectral domain D, a(w, w * ) = 0 identically. Thus, the second equation Alluding to the comments made below (3.14), Eq. (3.19) is an algebraic polynomial equation for b, of degree d + 1, where d is the number of distinct eigenvalues of B. A subset of these roots, subjected to appropriate boundary conditions, such as continuity and asymptotic behavior as w → ∞, should then be sewn together as branches of a single-valued holomorphic function b(w) defined throughout the complementary domain D c .
This multibranched structure of b(w) indicates that the complementary domain D c need not necessarily be a simply connected, or even a connected set. However, if the spectral domain D of φ is compact, D c should clearly have a connected subset D c ∞ which contains w → ∞. Thus, by considering (3.19) in the domain D c ∞ , we can read off the asymptotic behavior of b(w), assuming tr B = 0. This asymptotic behavior is consistent, of course with (2.18) and the definition (2.22), which imply that If, however, tr B = 0, we have to expand the RHS of (3.19) to the next order in b/w, and obtain which leads to a contradiction, unless b(w) vanishes identically in a neighborhood of w → ∞, and therefore must vanish identically throughout D c ∞ : Enforcing this result on the expansion of (3.21) in inverse powers of w ∈ D c ∞ , we come to the non-trivial conclusion, that if tr B = 0, then all moments (and in particular, all even moments) tr(Bφ n ) = tr(B(AB) n ) = 0, n ≥ 0. (3.23) Having solved for b(w), we then substitute it (together with a(w) = 0) in (2.32), to obtain the resolvent and thus, from (2.33), the desired Green's function in the holomorphic phase throughout D c . Comparing this equation with the equation defining the self-energy Σ herm (w) for ensembles of hermitian matrices, the interpretation (3.16) of −b as self-energy for PH matrices φ, alluded to in (3.16), becomes manifest in the holomorphic phase: By multiplying the numerator and denominator under the trace in (3.25) by B −1 and then using (3.19), we can rewrite (3.25) more compactly as The asymptotic behavior (3.20) of b(w) implies that in accordance with the definition of the resolvent of φ and the assumed compactness of D, and of course, with (3.11).
In the special case (3.22) we have identically throughout D c ∞ . Finally, the presence of a one-dimensional density component of real eigenvalues given by the discontinuity across the cut as (3.31) in accordance with (3.6).
Throughout the spectral domain D the order parameter a(w, w * ) = 0 and is pure imaginary, and we can rewrite (3.13) throughout this domain as By using the first equation in (3.32) in the second one, we can rewrite the latter as which leads to a contradiction, unless b(w, w * ) ≡ 0, and Thus, b(w, w * ) is pure imaginary throughout D, like a(w, w * ). These two functions are therefore manifestly not holomorphic functions of w in D.
One can easily verify that the two real equations (3.34), together with the first equation in (3.32), are completely equivalent to the original coupled equations (3.32).Therefore, we can simplify (3.32) further, and rewrite them as two real equations for two unknown real functions as 1 N m 2 tr Having solved for a(w, w * ) = iα(w, w * ) and b(w, w * ) = iβ(w, w * ), we then substitute them in (2.32), to obtain the desired Green's function (2.1) in the non-holomorphic phase as We can simplify this expression further into  .12).

III.2.3. The Phase Boundary Line
By definition, for probability ensembles with finite moments of the matrix elements of A, such as our Gaussian probability ensemble (1.4), the various averaged blocks Ĝ αβ and their traces (2.16) should be continuous functions of w and w * . As we cross the phase boundary from the non-holomorphic phase into the holomorphic phase, a(w, w * ) changes from being a non-zero pure imaginary function iα(w, w * ) into an identically vanishing function. This should happen continuously at the boundary line between the two phases, and therefore the phase boundary line Γ should be the solution of the equation That is, the phase boundary line Γ is the line of zeros of α(w, w * ). Due to the structure of the gap equations (3.13), as we cross Γ from the nonholomorphic phase into the holomorphic phase, b would then change continuously from being a non-analytic pure imaginary function iβ(w, w * ) into a holomorphic function b(w). Let's assume for example that β(w, w * ) is known explicitly in some subset of D, which borders a subset D c 1 of the holomorphic domain D c along an arc γ ⊂ Γ. In order to decide which branch of the solution of (3.19) we should assign to D c 1 we then use the fact that b must be continuous across the phase boundary. Thus, our knowledge of b(w, w * ) = iβ(w, w * ) along the boundary arc γ should fix the appropriate holomorphic branch of b in a unique way.
If it so happens that D c 1 forms a holomorphic island inside the non-holomorphic domain, so that γ is actually a closed curve, then we can compute the holomorphic b(w) inside D c 1 by invoking Cauchy's theorem as Let us derive this relation explicitly from the holomorphic gap equation (3.19), which in the present case simplifies into the quadratic equation from which it is also clear that we should pick the root which behaves asymptotically in accordance with (3.20). By comparing (3.45) with (3.29) we conclude (3.42), as required, because b(w) is analytic in w throughout the complex plane, save for a cut along the real axis. Thus, −b(w) = Σ GUE (w), which means that (3.25) is simply the statement that that is, is indeed the Green's function of the GUE ensemble, which leads to the semicircular eigenvalue density supported in the segment |x| ≤ 2 m , as required.

III.4. The Continuum Limit of the Density of Eigenvalues of the Metric B
As can be seen from Eqs. (3.19) and (3.25) in the holomorphic phase, and from Eqs.
of the metric. The main motivation for rewriting Eqs. (3.19), (3.25), (3.35) and (3.38) explicitly in terms of G B (w) is that in this form these equations are amenable to taking the continuum limit for ρ B (µ). That is, the limit in which the N → ∞ eigenvalues of B condense in a finite segment (or segments) on the real axis, rendering ρ B (µ) a continuous function in this domain 13 . Yet another advantage of such reformulation of the gap equations is that it leads to general results and also to simplification of known results such as (3.57) and (3.70) below, which would be difficult to derive otherwise.
By definition (3.49) is normalized according to resulting in the asymptotic behavior are the moments of ρ B (µ). Note also the obvious reflection property which will be useful for us below.

III.4.1. The Holomorphic Phase
In the holomorphic phase, we can clearly rewrite (3.19) where in the last step we have used (3.51). The integral in the last equation in (3.54) can be expressed in terms of (3.50), and therefore can finally rewrite (3.19) in terms of . (3.55) It can be easily checked, based on (3.52), that (3.55) is consistent with the asymptotic condition (3.20).
To summarize, given the metric B, its Green's function G B (w) is known in principle, and therefore (3.55) constitutes a functional equation for determining b(w).
Similarly, the Green's function G(w) (3.25) of φ in the holomorphic phase can be written as which means that the combination is invariant under the conformal mapping w → −w/b(w). By substituting (3.57) in (3.55) we obtain m 2 b 2 (w) + 1 = wG(w), which is nothing but (3.28). 23

III.4.2. The Non-Holomorphic Phase
After multiplying the numerator and denominator of the coupled non-holomorphic gap equations (3.35) by B 2 , we can rewrite them as where w = x + iy, and with α(w, w * ) and β(w, w * ) defined, respectively, in (3.33) and (3.36).
Let us now multiply the first equation in (3.58) on both sides by α 2 + β 2 . Then by adding and subtracting appropriate terms in its numerator, and after using the second equation in (3.58) once, we can recast it into the form where we have used (3.51), and where the real quantity ξ is defined as in terms of which we can rewrite the denominator in (3.59) as (3.63) and the factor multiplying ρ B (µ) in (3.61) as .
(3.65) By recalling (3.53) we can simplify these two coupled real equations further into where we used the definition (3.60) in the first equation. Alternatively, from (3.65), we can represent these two real equations as a single complex gap equation where we used (3.60) and (3.62). Note that the RHS of this equation is purely real. Thus, to summarize, given the metric B, G B (w) is known in principle. The latter is evaluated at ζ, which is a function of the non-holomorphic unknowns α and β. Therefore, the two coupled real gap equations (3.66), or equivalently, the single complex gap equation (3.68), determine α and β.
Once the non-holomorphic gap equations are solved for α and β, we could then use them to determine the resolvent G(w, w * ) of φ in the non-holomorphic phase. The latter is given by (3.38) as a trace over a function of the metric B. We can rewrite G(w, w * ) in terms of G B (ζ) by carrying algebraic manipulations similar to those carried above and obtain

III.4.3. The Unified Form of the Gap Equation
We can rewrite the holomorphic gap equation (3.55) as Comparison of (3.72) and its non-holomorphic counterpart (3.68) suggests that we extend ζ, defined in the non-holomorphic phase by (3.62) and (3.60), into the entire complex w-plane in a natural way asζ where it is understood that a = 0 and b(w) are holomorphic functions in the holomorphic regime w ∈ D c , and a = iα(w, w * ), b = iβ(w, w * ) in the non-holomorphic regime w ∈ D.
In this way, once the gap equations are solved for the appropriate quantities a and b, the two expressions (3.57) and (3.71) for the Green's function of the PH matrix φ in the holomorphic and non-holomorphic regimes, respectively, are unified into the single form wG(w, w * ) =ζG B (ζ) = 1 + m 2 (a 2 + b 2 ). (3.75) Thus, the combination wG(w, w * ) =ζG B (ζ) is invariant under the partly holomorphic, partly real-analytic mapping w →ζ.

III.4.4. The Phase Boundary Line Γ, Yet Again
Let us see what our reformulation of the gap equations in terms of G B (w) entails for the behavior of our solutions at the phase boundary between the holomorphic and non-holomorphic phases discussed in Section III.2.3. Approaching the boundary from within the non-holomorphic phase we know that α → 0 as well as b = iβ. Substituting these relations in (3.70) we immediately see that at the boundary wG(w, w * ) = 1 + m 2 b 2 , (3.76) 26 which is just (3.28) of the holomorphic phase, evaluated at the boundary. This should be expected, because the various holomorphic quantities should match continuously their non-holomorphic counterparts at the phase boundary. In this sense, (3.70) is quite natural, as it is a rather simple expression which crosses over continuously between the two phases at the phase boundary.
We can derive (3.76) in yet another (but related) way, which underlines the continuity ofζ in (3.73) across the phase boundary line, as follows: At the boundary, as can be seen from (3.60) and (3.62). By substituting (3.77) in (3.66), we can combine these two equations into in which both sides are evaluated at some point on the phase boundary line, where also b = iβ. Thus, we can rewrite (3.78) as which by (3.57), evaluated at the phase boundary coming from the holomorphic phase, reproduces (3.76).

III.4.5. Some Checks and the Continuum Limit
In the latter part of Section III.3 we have verified that our formalism reproduced the GUE results for B = 1 1 N . For this metric, of course, G B (w) = 1 w−1 , and one can easily verify that substituting this expression in (3.55) results in our previously derived quadratic equation (3.44) for b(w) in the GUE. Similarly, substituting this G B (w) in (3.57) reproduces the Green's function (3.46) of GUE.
A non-trivial check is offered by applying our results in this section to the metric (4.1) studied in the next section, for which we readily write We have verified that indeed, our results (3.55) and (3.57) in the holomorphic phase, and (3.66) and (3.70) in the non-holomorphic phase, reproduce all the analytic results reported in Section IV. As was mentioned in the Introduction, in this paper we focus on invertible metrics B. That is, we only discuss in this paper such densities ρ B (µ) which vanish in some neighborhood containing µ = 0. As a non-trivial example of such a continuous density, which is a simple generalization of (4.1), consider a metric whose positive and negative eigenvalues follow flat distributions along the positive and negative axis, namely, for some µ 1 > L − > 0 and µ 2 > L + > 0. For this ρ B (µ) we obtain where the logarithms are defined in the cut plane, with a cut emanating from each branch point and running along the real axis in the negative direction. With this assignment of the cuts, the imaginary part of (3.82) will have the correct discontinuity as in (3.6) along the support of (3.81). Substituting with 1 < k < N . Our discussion in this section of the metric (4.1) will be somewhat telegraphic, avoiding detailed derivation of some of the analytical results and presenting only a limited set of numerical results. We shall dedicate a subsequent paper [15] to providing the full details of the relevant analytical derivations, support them by results of ample numerical computations, including numerical demonstration of the statistical universality of this model, which is sensitive in the large-N limit essentially only to the second cumulant (2.14) of the random matrix A.
With the limit N → ∞ in mind, let us define the fraction of the +1's on the diagonal of B as λ = k N .  For the particular choice of metric B given in (4.1), the holomorphic gap equation (3.19) reduces to the cubic equation where according to (3.20), we must pick that root of (4.4) which behaves asymptotically We can rewrite (4.4) as Let us evaluate this cubic equation for w = x + i0+ with x ∈ R, i.e. when w approaches the real axis from above. Then for fixed x, the coefficients of the cubic equation (4.6) are real and its three possible solutions for b are either all real or one of them is real and the other two come as a complex conjugated pair, depending on the sign of the discriminant −∆/(27m 4 ), where We see from (4.7) that there are three real roots when |x| is large. Furthermore, ∆(x) changes its sign at x = ±x 0 , where (4.8) We thus conclude that along the real axis (4.6) has one real and a pair of complex conjugate solutions for b(x + i0+) when x ∈ [−x 0 , x 0 ], and three real solutions otherwise. For the particular metric (4.1), the general expression (3.25) for the Green's function in the holomorphic phase yields By substituting the appropriate root of (4.4) in (4.9), we can obtain the average density ρ (1) (x) of real eigenvalues of φ from the discontinuity of this function across its cut along the real axis according to (3.31) in the large N limit. After a tedious but straightforward calculation, we find that ρ (1) (x) is supported along the interval [−x 0 , x 0 ] where it is given by where ∆ and ξ have been defined in (4.7) and x 0 in (4.8). In particular, at the center of the band of real eigenvalues the density is given by (4.11) In Fig. 1 we show normalized histograms obtained from numerically generated samples for various values of λ. For each value of λ, we have used 2000 samples of matrix size N = 1024. The histograms fit well with the solid black lines from the corresponding theoretical limits N → ∞ given by Eq. (4.10). A small deviation due to finite N effect can be seen only in the histogram in magenta showing the result for λ = 63/128, which is close to the degenerate case λ = 1/2, for which the theoretical density (4.10) predicts that ρ (1) (x) ≡ 0. As λ decreases from 1/2, the density of real eigenvalues increases. The blue histogram corresponds to having λ close to zero. For λ = 0 we obtain of course Wigner's semi-circular distribution, because then φ = A, which has been drawn from the GUE.
According to (3.11), the fractions of all real and all complex eigenvalues sum up to one. Integrating ρ (1) (x) yields the fraction 1 − ν of real eigenvalues as a function of λ, in accordance with (3.7). Carrying out this integral explicitly is a rather difficult task, which we can avoid by the much easier calculation of the complementary fraction of complex eigenvalues in (4.23) in the next section IV.2. Thus, we find that the fraction  (4.12) which is consistent with the symmetry (4.3) of the model. In Fig. 2 this result has been compared with numerical simulations for matrix sizes N = 128, N = 1024 and N = 8192 averaged over many samples. For λ away from 1/2 the convergence is exponentially fast. Near 1/2, we can see a small deviation for finite N . The theoretical large N prediction, actually is a lower bound for the number of real eigenvalues, even for finite N and for each sample. This follows from a special case of a purely algebraic theorem [13].

IV.2. The Non-Holomorphic Phase, a = 0
In this phase a(w, w * ) = iα(w, w * ) and b(w, w * ) = iβ(w, w * ) are pure imaginary. By substituting the metric B given by (4.1) in the gap equations (3.32), or equivalently (3.35), we thus obtain (4.13) or equivalently (4.14) Subtraction of these two equations yields a linear equation for β from which we determine where we substituted w = x + iy. Finally, feeding (4.15) back in (4.14) we determine According to (3.39), the boundary line between the holomorphic and non-holomorphic phases is the curve along which (4.16) vanishes, which is conveniently represented in polar coordinates x = r cos θ, y = r sin θ by where sin θ 0 = |2λ − 1|. This curve is made of two closed loops, one in the upper halfplane, and its mirror image with respect to the real axis in the lower half-plane. Each loop is made of two arcs (the roots of (4.17)) where sin 2 θ ≥ sin 2 θ 0 . r + is the part of the boundary farther from the origin, and r − the closer one. The two arcs are sewn together at sin 2 θ = sin 2 θ 0 to make the loops. These two loops enclose two two-dimensional blobs which comprise the domain D of the non-holomorphic phase that contains the condensed pairs of complex conjugate eigenvalues of φ in the large-N limit.
We can also derive the boundary equation (4.17) coming from the holomorphic phase: As we approach the phase boundary curve, the holomorphic solution b(w) of (4.6) should become purely imaginary. Thus, by substituting the boundary value b = iβ(x, y) in (4.6) and separating the resulting cubic into its real and imaginary parts, we obtain the two coupled equations By comparison with (4.15), we see from the second equation that the holomorphic b(w) coincides with the non-holomorphic solution at the boundary, as it should, by continuity. Then, by substituting this boundary value of β in the first equation in (4.19) we simply obtain the Cartesian form of (4.17) (multiplied by 2λ − 1). The particular metric (4.1) has the special property that B 2 = 1 N , which upon substitution in (3.38) and use of the first gap equation in (4.13), leads to the particularly simple non-holomorphic Green's function valid throughout D. At the phase boundary line (4.18) the non-holomorphic Green's function (4.20) has to cross-over continuously to the holomorphic Green's function (4.9). Note in passing that (4.20) implies wG(w, w * ) = m 2 w * w. By comparing this with (3.70) we thus conclude that α 2 + β 2 = 1 m 2 − x 2 − y 2 , which is consistent with (4.16). By applying Gauss' law to (4.20), as in (3.12), we see that the density of complex eigenvalues is uniform in D, and also independent of λ.

33
The area of the two blobs of D, enclosed by the two loops (4.18), is given by By multiplying this area and the density (4.21), we find that the fraction ν of complex eigenvalue is given by The scatter plots of complex eigenvalues in Figs. 3 and 4 illustrate numerically the results for various values of λ. The blue dots correspond to complex eigenvalues of φ and the red dots correspond to real eigenvalues. In Fig. 3 we used 64 samples of matrix size N = 128. Fig. 4, on the other hand, was obtained from a single sample of matrix size N = 8192. That a single large sample is enough to represent reliably the averaged density is the result of the self-averaging phenomenon, typical of the eigenvalue density of large random matrices. It demonstrates very well how uniformly the eigenvalues are distributed in the domain D. The solid black lines in both figures represent the theoretical boundaries given in Eq. (4.18). Finite-N corrections are pronounced in Fig. 3 as one can clearly see complex eigenvalues outside of the theoretical domain D. In contrast, for N = 8192 in Fig. 4, one can hardly see any outliers. Scatter plots for N = 32786 can be found in [10].
On the lower right plot in Fig. 3 and 4, we see the degenerate case λ = 1/2, for which the theoretical large-N prediction (4.10) (see also (4.12)) gives null density of real eigen- values. Finite-N corrections for λ = 1/2 will result in small amounts of real eigenvalues, as is demonstrated very well by the corresponding plates in Figs. 3 and 4, and in accordance with Fig. 2. In all other plots in Figs. 3 and 4 we see pronounced distributions of real eigenvalues, with the red dots appearing to form a continuous line. This is so because they are very dense as there are |N − 2k| = N |1 − 2λ| real eigenvalues packed in the finite segment [−x 0 , x 0 ]. For fixed λ, the typical distance between eigenvalues on the real lines is of order N −1 while the mean distance between nearest neighbors of complex eigenvalues in the bulk is of order N −1/2 . For each sample, the eigenvalue distribution is symmetric with respect to the real axis. This is due to the fact that the characteristic polynomial of det(z −φ) has real coefficients and the complex eigenvalues appear in conjugate pairs. It is a manifestation of (broken) PT symmetry. The eigenvalue distribution has an additional reflection symmetry with respect to the imaginary axis, but only after averaging over all samples, because due to (1.4) (or (1.5)), φ and −φ are equi-probable.
As was mentioned above (following (4.23)), at N = 8192 we can already see in Fig. 4 that the complex eigenvalues are spread quite uniformly across the blobs, which agrees with our finding in (4.21) that the theoretical density there is constant. In Fig. 4 of [10] one can see two dimensional histograms which show that the density is very flat in the bulk. Only at the edges of those histograms one can notice some deviation from the constant density which is a typical finite N size effect, that appears also in the Ginibre ensemble or in normal matrix models [28,29].
From (4.18) we find that the distance between the complex domain and the real axis (on both sides) is r − (π/2) = m −1 sin(θ 0 /2). Thus, the two blobs kiss and touch the real axis when λ approaches 1/2. The latter is the degenerate case which we have mentioned above, and in which the complex eigenvalues fill uniformly a disk, like in the (real) Ginibre ensemble. This is of course the case tr B = 0 alluded to in (3.22), with b ≡ 0 and G(w) = 1/w throughout the holomorphic region r > 1/m.

IV.3. The Non-Compact
Algebras su(k, N − k) and so(k, N − k), and Pseudo-and Anti-Pseudo-Hermiticity The metric (4.1) is of special importance, because it appears in the definition of the non-compact Lie groups U (k, N −k) and O(k, N −k). By definition, the group U (k, N −k) preserves the inner product u † Bv of vectors u, v ∈ C N , that is, , which mix the two sub-spaces. In addition, due to their "chiral" off-diagonal rectangular structure, they also possess at least |N − 2k| "kinematical" zero modes which are typical of rectangular matrices. In "physical" terms, think of the subspace orthogonal to such a complexified boost, which remains invariant. In other words, the corresponding boost generator annihilates these vectors. The full matrices φ, which dominate the ensemble and whose spectra was discussed in Sections IV.1 and IV.2, are the direct sum of the hermitian compact generators and anti-hermitian non-compact generators. Thus, in essence, the formalism presented in this paper solves this non-trivial addition problem of the direct sum of block-diagonal hermitian random matrices and block-off-diagonal antihermitian random matrices.
The random variable t = 1 N tr φ = 1 N tr AB, which is a weighted sum of the Gaussian variables in (1.4), follows itself the Gaussian distribution (4.28) Therefore (assuming a finite limit for u), and t becomes deterministically zero in the large-N limit. Thus, more precisely, in this limit, our PH random matrix model is effectively an ensemble of the traceless part of u(k, N − k), and our results in this section describe the spectral properties of elements of the Lie algebra su(k, N − k) in the large-N, k limit. The cases k = 0 and k = N , that is B = ±1, correspond to φ being a GUE matrix. Indeed, GUE can be thought of as an ensemble of elements of randomly chosen generators of the algebra su(N ).
In a similar manner, by definition, the group O(k, N − k) preserves the inner product u T Bv of vectors u, v ∈ R N , that is, in terms of a real matrix φ, an element of the Lie algebra so(k, N − k), we can rewrite (4.30) as e φ T Be φ = B, (4.32) which is equivalent to the intertwining relation meaning that φ is anti-pseudosymmetric (APS) with respect to B.
Similarly to the su(k, N − k) case discussed above, one can easily check that such APS matrices φ are the direct sum of block-diagonal real antisymmetric matrices and block-off-diagonal real symmetric matrices. The antisymmetric ones have diagonal blocks of dimensions k × k and (N − k) × (N − k), which generate compact SO(k) ⊗ SO(N − k) rotations within the two invariant subspaces. Their eigenvalues are purely imaginary, of course. The real symmetric matrices have rectangular blocks of dimensions k×(N −k) and (N − k) × k, with real eigenvalues. They generate non-compact transformations (boosts), which mix the two sub-spaces. In addition, due to their "chiral" off-diagonal rectangular structure, they also possess at least |N − 2k| "kinematical" zero modes, because boosts do not affect directions perpendicular to them. The full matrices φ, which dominate the ensemble, are the direct sum the antisymmetric compact generators and real symmetric non-compact generators.
If φ is APS, the pure imaginary matrix iφ is PH and satisfies the intertwining relation (1.1).
In [15] we shall bring strong numerical evidence in support of universality of our results for the spectral density of the ensemble of PH matrices discussed in this section. In particular, we shall demonstrate that replacing the GUE ensemble (1.4) with a properly rescaled GOE ensemble, that is, taking the matrices A to be random real symmetric matrices, does not affect the large-N results for the eigenvalue density of the real matrices φ = AB. In this case, φ T B = Bφ are pseudo-symmetric with respect to B.
Arguing along the lines of universality of the average spectra of real and complex PH hermitian matrices φ in the large-N limit, we should expect similar universality in the spectra of anti-PH real and complex matrices as well. Any complex anti-PH matrix can be written as iφ, with φ a complex PH matrix. Thus the large-N spectrum of complex anti-PH matrices is given by our results for PH matrices, rotated by 90 degrees in the complex plane. Spectral universality thus tells us that the average spectrum of randomly chosen generators of so(k, N − k) in the large-N, k limit should be just the spectrum obtained in this section, rotated by 90 degrees in the complex plane. 37

V. Discussion
In this paper we have introduced a family of pseudo-hermitian random matrices as a new concept in random matrix theory. The conceptual novelty in these models is that PH matrices are hermitian with respect to a given indefinite metric. Major applications of pseudo-hermiticity include (and is certainly not limited to) the studying P T -symmetric quantum mechanical systems in the phase of broken P T -symmetry, as well as their classical gain-loss balanced analogs. Thus, in the usual spirit of Random Matrix Theory [31], randomness of the PH matrices discussed in the present paper can model complicated or chaotic systems, or systems which are disordered to begin with. We shall study such applications of PH random matrix theory in future publications.
In this paper we have applied the diagrammatic method to derive the self-consistent gap equations in the planar limit, from which we determined the phase structure of the model in the complex plane. In particular, we have derived analytical expressions for the average Green's function of the model in the holomorphic and non-holomorphic phases as functions of the indefinite metric, which can be used to determine the average density of eigenvalues of such matrices both on the real axis and in the complex plane. As a concrete example, we have applied our formalism to studying a family of PH random matrix models for the generators of the classical non-compact Lie algebra su(k, N − k) in the large-N, k limit, with finite λ = k/N . We have calculated explicitly their average density of eigenvalues on the real axis and in the complex plane. We have also carried meticulous numerical analysis of these matrices, part of which was presented here. The numerical results agree very well with our analytical predictions. More analytical and numerical details pertaining to this model will be given in a forthcoming paper [15], which will also demonstrate universality of our results for the density of eigenvalues.
A preview of the main idea and results of this paper was published in [10], which also reviewed our results on quasi-hermitian random matrices and their application to the vibrational spectra of mechanical systems with large connectivity [8,9], as well as analytical and numerical results pertaining to yet another PH random matrix model with very rich and interesting phase structure [32].
The results presented here and in [8,9,10] are just the tip of the iceberg. Quasiand pseudo-hermitian random matrices are clearly a promising new direction in the vast ocean of Random Matrix Theory.
We see that each eigenvalue w i of φ leads to a pair of eigenvalues ± √ w i of H. Moreover, as was mentioned in the Introduction, isospectrality of φ and φ † also means that their complex eigenvalues come in complex conjugate pairs. Thus, if w i = 0, then w * i is also an eigenvalue of φ (distinct from w i ), and therefore ± w * i are eigenvalues of H as well. Therefore, while the PH nature of φ renders its spectrum symmetric under reflection with respect to the real axis, the spectrum of H has four-fold reflection symmetry, with respect to both the real and imaginary axes.
Another way to deduce this symmetry, without resorting directly to the spectral properties of φ, is to depart from the chiral block structure of (2.2) which implies is the chirality matrix. Thus, if ψ is a (right) eigenvector of H with eigenvalue z, then Γψ is also an eigenvector of H, with eigenvalue −z. This accounts for the reflection symmetry z ↔ −z of the spectrum of H through the origin. Moreover, as can be seen from (A.2), the LHS of (A.3) is invariant under interchanging A ↔ B, which from (2.2) is equivalent to interchanging H ↔ H † . Thus, the spectrum of H is invariant under complex conjugation as well, and therefore under reflection with respect to the real axis. 39 In the large-N limit, the eigenvalues of H typically become dense in a two-dimensional domain D in the complex plane. More precisely, as it turns out, due to the PH nature of φ, the spectral domain of H also contains a one-dimensional density component of eigenvalues concentrated in finite segments along the real and imaginary axes, the union of which we denote byσ. We have shown here that the domain D ∪σ must have four-fold reflection symmetry with respect to both the real and imaginary axes. In this section we will establish a certain functional relation between the diagonal blocks of the resolvent (2.6), which hold true for every realization of the random matrix A, that is, before taking averages. To this end we evaluate (2.6) at η = is, with s ∈ R, that is,Ĝ . (A.6) Let us now substitute H from (2.2) in (A.6) and start by computing the hermitian and positive definite inverse and find the corresponding N × N blocks as P (s; z, z * ) = P † (s; z, z * ) = s 2 + |z| 2 + A 2 − (zB + z * A) 1 As should be clear from the structure of the second matrix in (A.7), by interchanging A ↔ B we have P (s; z, z * ) ↔ R(s; z, z * ) and Q(s; z, z * ) ↔ Q † (s; z, z * ). (To see that interchanging A and B in the last expression in (A.8) has the effect of transforming Q to its adjoint requires a little bit of reshuffling of the various matrix factors.) Similarly, we find that is, the result (A.7) after simultaneously interchanging z ↔ z * and A ↔ B. By collecting our results from equations (A.7)-(A.9), substituting them in (A.6) and using the definition of the N × N blocks (2.9), we thus conclude that Upon tracing each N × N block in these equations, we also obtain the corresponding expressions for the unaveraged block traces αβ(is; z, z * ) = 1 N trĜ αβ (is; z, z * ), α, β = 1, 2, 3, 4 (A.11) (whose averaged versions appear in (2.16)) in terms of traces over P , Q and R.
Clearly, the sums of diagonal block traces on the LHS of (A.17) do not vanish if and only if z is an eigenvalue of H (and therefore z * an eigenvalue of H † ). Thus, they serve as sort of order parameters indicating the location of the support D ⊂ C of the density of eigenvalues of H in the large-N limit [19]. The discussion in (section 3 of) [19] focused on non-hermitian matrices H with purely two-dimensional spectral domain D. For example, if H is drawn from Ginibre's ensemble [30], normalized in the large-N limit to have the unit disk as its spectral domain, then which should hold true to high precision even before averaging over the ensemble, due to the self-averaging property of large random matrices. Going beyond Ginibre's ensemble to general rotationally invariant ensembles of non-hermitian matrices, the Single Ring Theorem [25] asserts that the spectral domain D of H is either a disk or an annulus, and for such models as well, one can compute the diagonal trace on the LHS of (A.19) (see Section 4 in [25]) and show that it vanishes only outside D. Thus, this diagonal trace indeed fulfills its role as order parameter indicating the location of the two-dimensional spectral domain. However, as was mentioned at the end of Section A.1, the spectral domain of H in (2.2) also contains a one-dimensional eigenvalue density componentσ, located symmetrically along segments of the real and imaginary axes. We now show that this one-dimensional spectral component has a negligible effect on the order parameter, producing a benign removable discontinuity, invisible in the large-N limit. To get oriented, take as a start H to be a hermitian matrix -the opposite extreme of a Ginibre-type matrix -whose eigenvalues x 1 , . . . , x N condense in the large-N limit in some finite segmentσ along the real axis, with density ρ (1) (x) as in (3.4). In this case, it is straightforward to show that As a function of y, this diagonal trace is null, except a removable discontinuity at y = 0, and fails to serve as an order parameter indicating the position ofσ. This conclusion holds true also in the more generic case, when the spectral domain of H is a union of both a one-dimensional partσ and a two-dimensional part D.
In order to proceed, let us carry a singular value decomposition of the matrix z − H In this case, we immediately see that V † (z k )ê k is an eigenvector on the right of H, belonging to the eigenvalue z k , and U(z k )ê k is the corresponding eigenvector on the left (withê k the kth Cartesian unit vector). By substituting (A.21) in (A.17), we see that the two equal diagonal traces on the RHS of (A.17) When we set, for example, y = 0 and go off the real axis, purely real eigenvalues make no contribution to the trace, while for y = 0, in the large-N limit, they only give rise to a removable discontinuity, because then their contribution converges to a smooth function, as in (A.20). Obviously, a similar conclusion holds for the purely imaginary component inσ. Finally, note that for H a normal matrix, and in particular, a hermitian matrix, clearly Λ k (z) = |z − z k |, and therefore c k = 1 ∀k, which is consistent in the large-N limit with (A.20). Thus, to summarize the discussion of the few previous paragraphs, the sums of diagonal block traces on the LHS of (A.17) serve as order parameters indicating the location of the two-dimensional support D ⊂ C of the density of eigenvalues of H in the large-N limit, and are effectively insensitive to the one-dimensional partσ of the spectral domain.
A corollary of this assertion is that both these sums must vanish (modulo the removable discontinuity alongσ) outside D, that is, in the complementary domain D c = C/ D. All the αα(is = 0; z, z * )'s in (A.17) are of the same sign (depending on whether s → 0+ or 0−). Therefore, vanishing of such a sum implies vanishing of each of its terms separately. We conclude that all four αα(is = 0; z, z * ) must vanish in D c .
In contrast, at this point, before taking averages, we cannot determine whether each of the four diagonal traces αα(i0; z, z * ) = 0 separately for z ∈ D. All we can say at this point is that the two equal sums in (A.18) do not vanish ∀z ∈ D. (The identities (A.14) and (A.15) we were able to establish, merely interchange z ↔ z * in (A.18).) We will be able to answer this question in the affirmative (see (2.27) in the main text) only after averaging over A.

A.3. Derivation of the Self Energy in the Planar Limit
Let us start by deriving a useful elementary algebraic identity. To this end we introduce the standard basis vectors in matrix space (ê ij ) kl = δ ik δ jl , (A. 31) namely, the N × N matrix whose all entries are null, except the ij-th element. Then, it is straightforward to show that (repeated indices are summed over) (ê ij Mê ji ) pq = δ pq trM, (A. 32) for any N × N matrix M.
In order to proceed, we write the inverse of (2.8) aŝ where the subscript c indicates the connected second order cumulant (2.14). Thus, (A.37) Multiplying through, and displaying the blocks (2.9) of Ĝ on the way, we thus obtain