Unquenched QCD Dirac Operator Spectra at Nonzero Baryon Chemical Potential

The microscopic spectral density of the QCD Dirac operator at nonzero baryon chemical potential for an arbitrary number of quark flavors was derived recently from a random matrix model with the global symmetries of QCD. In this paper we show that these results and extensions thereof can be obtained from the replica limit of a Toda lattice equation. This naturally leads to a factorized form into bosonic and fermionic QCD-like partition functions. In the microscopic limit these partition functions are given by the static limit of a chiral Lagrangian that follows from the symmetry breaking pattern. In particular, we elucidate the role of the singularity of the bosonic partition function in the orthogonal polynomials approach. A detailed discussion of the spectral density for one and two flavors is given.

The microscopic spectral density of the QCD Dirac operator at nonzero baryon chemical potential for an arbitrary number of quark flavors was derived recently from a random matrix model with the global symmetries of QCD. In this paper we show that these results and extensions thereof can be obtained from the replica limit of a Toda lattice equation. This naturally leads to a factorized form into bosonic and fermionic QCD-like partition functions. In the microscopic limit these partition functions are given by the static limit of a chiral Lagrangian that follows from the symmetry breaking pattern. In particular, we elucidate the role of the singularity of the bosonic partition function in the orthogonal polynomials approach. A detailed discussion of the spectral density for one and two flavors is given.

I. INTRODUCTION
One of the features that makes QCD at nonzero baryon chemical potential both elusive and interesting is that the Euclidean Dirac operator does not have any Hermiticity properties. The non-Hermiticity occurs in an essential way because a finite baryon number density is obtained by promoting the propagation of quarks in the forward time direction and inhibiting the propagation in the backward time direction. Because of the non-Hermiticity the quark determinant attains a complex phase, which prevents the analysis of the QCD partition function by means of probabilistic methods. An approximation that is often used at zero chemical potential is to ignore the fermion determinant altogether. However, this approximation fails dramatically at nonzero chemical potential [1]. In particular, the critical chemical potential at zero temperature is determined by the pion mass instead of the nucleon mass. Therefore, the main source of information for QCD at nonzero baryon density is based on simplified models or on perturbative expansions at very high densities which may never be accessed experimentally. Lattice QCD simulations are feasable in a region where the ratio of the chemical potential and the temperature is sufficiently small so that extrapolation from an imaginary chemical potential [2,3] or from zero chemial potential [4] becomes reliable (for a critical review of the field and additional references, we refer to [5]).
In this paper, we will investigate an observable for which nonperturbative results can be obtained for unquenched QCD at nonzero baryon density. This observable is the spectral density of the QCD Dirac operator. At nonzero chemical potential its support is a two-dimensional domain. The naive argument that observables do not depend on the chemical potential at zero temperature does not apply to the spectral density of the Euclidean Dirac operator. This puzzling fact can be understood as follows. To define an eigenvalue density we need both the eigenvalues and the complex conjugate eigenvalues or the Dirac operator and its complex conjugate. However, complex conjugation is equivalent to changing the sign of the chemical potential. Therefore, the generating function for the QCD Dirac spectrum has to contain both quarks and conjugate quarks in the valence sector that have an opposite baryon charge. This opens the possibility that the low energy limit of this theory contains Goldstone bosons made out of quarks and conjugate anti-quarks which have a nonzero baryon number [6,7]. More formally, the chemical potential in the generating function breaks the flavor symmetry resulting in a chiral Lagrangian that becomes dependent on the chemical potential. The region of the Dirac spectrum we wish to analyze determines the mass of the valence quarks [8,9]. In order to determine the eigenvalue density of the Dirac operator at small eigenvalues the masses of the valence quarks are equally small. Consequently, the mass of the charged Goldstone modes is small and the generating function depends on the chemical potential even at small values of µ.
The low energy effective theory for the generating function of the Dirac spectrum is a theory of Goldstone modes determined by spontaneous breaking of chiral symmetry. This effective theory is a version of a chiral Lagrangian [10] which takes into account the presence of the conjugate quarks. The allowed terms in the effective theory are determined by the flavor symmetries and the way they are broken by the quark masses and the chemical potential: the effective theory must break the flavor symmetries in precisely the same way as in the generating function we started from. To be specific, the Goldstone field, U , takes values in the coset U ∈ SU (N f + 2n), where n is the number of additional quark -conjugate quark pairs. The leading order chiral Lagrangian in the standard counting scheme is the usual non-linear sigma model [11]. The dependence of the chemical potential is completely fixed by the flavor symmetries. It enters as the zeroth component of an external vector field through a shift of the Euclidean time derivative (see e.g. [10,12,13,14,15]) where µ is the chemical potential and B is the charge matrix. Since the baryon charge of conjugate quarks is opposite to the standard quarks [6] the charge matrix B is not proportional to the identity and the commutator is non-vanishing. Although the eigenvalue spectrum of the Dirac operator is not a direct physical observable, we hope that detailed knowledge of this observable in a nonperturbative domain where no other information is available, will ultimately lead to a better understanding of the problems that hinder numerical simulations at nonzero baryon chemical potential. At zero temperature these problems are manifest unless µ 2 F 2 π V ≪ 1, where F π is the pion decay constant. In this paper, we will focus on a scaling regime where the product µ 2 F 2 π V is fixed as the volume, V , is taken to infinity. The fixed combination µ 2 F 2 π V can take any value and our results will show the effect of the fermion determinant on the Dirac spectrum. The quark masses will be taken such that m π scales like µ with m 2 π F 2 π V kept fixed. To be specific, we consider the range [8,16] as the volume L 4 = V is taken to infinity, and Λ is the scale of the lightest non-Goldstone particle. This regime is sometimes known as the ergodic domain or as the epsilon regime 1 of chiral perturbation theory [17]. In this domain the zero modes of the Goldstone fields dominate the partition function which reduces to a static integral over the Goldstone manifold [14,17]. The static integral is completely determined by the flavor symmetries of the QCD partition function. This implies that any theory with the same flavor symmetries and flavor symmetry breaking is described by the same static integral. The only memory of the underlying theory is two coupling constants, namely the chiral condensate Σ and the pion decay constant F π . Hence the partition function in this limit is universal in the sense that all theories with a given symmetry breaking pattern and a mass gap will have a partition function given by the same static integral over the Goldstone manifold provided that m π , µ ≪ 1/L ≪ Λ [9,18,19,20]. The simplest theories in this class are invariant random matrix theories in the limit of large matrices. Because of the large invariance group it is sometimes simpler to analyze the random matrix theory rather than to evaluate the integral over the Goldstone manifold directly. In order to derive the microscopic eigenvalue density of the Dirac operator, which is defined by rescaling its eigenvalues as z k V Σ, we start from the microscopic partition function with N f + 2n flavors. Then we take derivatives with respect to the mass of the additional n quarks and n conjugate quarks, and finally we remove the additional flavors by taking the limit n → 0. This procedure is known as the replica trick [21]. Since the microscopic partition functions are only known for integer n there is no guarantee that a correct nonperturbative answer is obtained this way. After a two decade long discussion [22] two closely related methods have emerged that result in the correct nonperturbative results: the replica limit of the Painlevé equation [23,24] and the replica limit of the Toda lattice equation [25,26]. Both the Painlevé equation and the Toda lattice equation are well-known in the theory of exactly solvable systems (see for example [27,28]). The Painlevé equation is a complicated nonlinear differential equation, whereas the Toda lattice equation is a simple two step recursion relation. For that reason, it is much simpler to work with the Toda lattice equation. The advantage of the Painlevé equation is that nonperturbative results can be obtained from fermionic partition functions only [23]. Recently, the consistency of the replica limit of the Toda lattice equation [29] and the supersymmetric method [30] was established. A necessary ingredient for the applicability of the Toda lattice method is that the QCD partition function satisfies the Toda lattice equation. This is the case if the QCD partition function is a τ -function [28]. For the low-energy limit of QCD at µ = 0 this was shown for fermionic partition functions in [31,32]. These results were extended to bosonic and supersymmetric partition functions as well as to QCD at nonzero chemical potential in [25,26,29].
In the replica limit of the Toda lattice equation, generating functions with bosonic quarks appear in addition to generating functions with fermionic quarks. While the fermionic generating functions for the spectral density at nonzero baryon chemical potential are known from the integral over the Goldstone manifold [26,33], only the simplest case with bosonic quarks is known [26]. The relevant supersymmetric generating functions, for which the Goldstone manifold is a supermanifold, have proved to be quite challenging to compute, and no explicit expressions have been obtained so far. In the present paper we derive these generating functions from the random matrix model introduced in [34]. The flavor symmetries and their spontaneous and explicit breaking are exactly the same in the large N limit of this random matrix model and in full QCD. Since both theories have a mass gap, in the microscopic limit, they are both given by the same integral over the Goldstone manifold.
In the framework of random matrix models, spectral correlation functions can also be computed by means of (bi-)orthogonal polynomials in the complex plane. The complex orthogonal polynomial approach was developed in [35,36,37,38,39,40] and was applied to unquenched [34] QCD Dirac spectra at nonzero chemical potential. The outcome of our analysis here is that both the Toda lattice and the orthogonal polynomials approach produce equivalent results. The divergence of the bosonic generating functions also sheds some light on universality of non-Hermitian random matrix models in general.
Spectra of the non-Hermitian Dirac operator at nonzero baryon density have been obtained from lattice QCD in the quenched case [1,41,42,43] and for QCD with two colors [44,45] and have been compared successfully to random matrix theory [41,42,43,45]. The microscopic spectral density of quenched lattice QCD at nonzero baryon density was first analyzed in [43] where quantitative agreement with analytical predictions [26,37] was found in an asymptotic domain where the results derived from the chiral Lagrangian [26] agree with the expression obtained in [37]. The two flavor phase quenched partition function which is calculated in this paper does not suffer from a sign problem either and could be simulated on a lattice (see [46,47] for recent results). A first analytical prediction for the unquenched QCD Dirac spectrum at nonzero baryon chemical potential in the non-perturbative regime was obtained in [34]. These results and the work presented in this paper could in principle be compared to numerical simulations in the region where the sign problem sets in. In particular, any proposal for a solution of the sign problem on the lattice can be tested against these predictions.
For the spectral density with one dynamical flavor we show both analytically and numerically that large fluctuations occur for increasing value of the chemical potential. These findings further illustrate the difficulties encountered when one wishes to reproduce these results in lattice QCD simulations.
In section II we explain the derivation of the spectral density of the QCD Dirac operator within the replica framework. The required generating functions are all given in terms of group integrals that represent low-energy effective QCD-like partition functions. The bulk of this paper is devoted to the calculation of the group integrals corresponding to bosonic quarks using the large N limit of a random matrix model. After introducing the relevant random matrix model for QCD with chemical potential in section III, we compute the necessary generating functions by means of complex orthogonal polynomials in section IV. In section V we show that these partition functions satisfy the Toda lattice equation. This allows us to calculate the spectral density from the replica limit of the Toda lattice equation. We show in general terms that the results agree with a direct computation using orthogonal polynomials [34]. In section VI we present our results for one, two and any number of dynamical flavors and compare them to the quenched and phase quenched results. It also contains a discussion of the thermodynamic limit of our results. This section is self contained and can be read independently of the previous sections. Our main findings are summarized in the conclusions, and several technical details are worked out in two appendices.

II. THE SPECTRAL DENSITY FROM GENERATING FUNCTIONS
At nonzero baryon chemical potential the Dirac operator is non-Hermitian so that its spectrum has a twodimensional support in the complex plane. In this section we start with the definition of the spectral density for complex eigenvalues. Then we show how it can be obtained from a generating function with additional flavors.
We consider the eigenvalues of the Dirac operator given by the eigenvalue equation Since (D η γ η ) † = −D η γ η and (µγ 0 ) † = µγ 0 the Dirac operator is neither Hermitian nor anti-Hermitian and the eigenvalues are complex. Because of the axial symmetry {D η γ η + µγ 0 , γ 5 } = 0 the nonzero eigenvalues come in pairs with opposite sign, ±z j . The density of eigenvalues in the presence of N f flavors is defined as the vacuum expectation value over a sum of delta functions in the complex plane at the positions of the eigenvalues, vanishing and non-vanishing, where µ is the chemical potential, N f is the number of flavors and {m f } = m 1 , . . . , m N f are the quark masses. For later convenience we have not included the zero eigenvalues in our definition of the spectral density. The average over gauge fields in a fixed topological sector, ν, is defined by Since the phase of the product of the fermion determinants is non-vanishing, such expectation values of operators will in general be complex. In particular we will find that the unquenched spectral density is complex rather than real non-negative. The two-dimensional δ-functions appearing in the definition of the eigenvalue density can be obtained from the derivative with the resolvent defined by the average where δ z = zd/dz and δ z * = z * d/dz * . Since the replica limit of the right hand side of this equation is the spectral density 2 we find the remarkably simple expression [25,26,29] In (16) and elsewhere in this paper it is our notation to separate the fermionic and bosonic quark masses in the argument of the partition function by a vertical line. In addition to the fermionic partition functions, we will need to evaluate Z N f ,n ν for n = −1, i.e. for one bosonic quark and one bosonic conjugate quark. In the quenched case (N f = 0), this partition function is given by the integral [26,49] where dU θ(U )/det 2 U is the integration measure on positive definite 2 × 2 Hermitian matrices and It was found that in order to obtain a finite limit ǫ → 0, the normalization constant has to be chosen C ǫ ∼ 1/ log ǫ. The appearance of this logarithmically diverging term is not surprising since the inverse fermion determinants are regularized as Zero eigenvalues of D η γ η do not lead to zero eigenvalues of the operator in the r.h.s. of this equation. Therefore, generically its eigenvalues λ k are nonzero real and occur in pairs ǫ ± λ k , and we expect that the ǫ dependence enters in a similar way to the mass dependence of the QCD partition function with one bosonic quark at zero topological charge. In the microscopic limit this is given by K 0 (m) which diverges logarithmically in m.
To derive the partition functions (14) and (17) we have only used the global symmetries of the underlying QCD partition functions. In particular, this implies that any microscopic theory with the same global symmetries and spontaneous breaking thereof will have the same zero momentum effective theory. The simplest theory in this class is chiral Random Matrix Theory at nonzero chemical potential which is obtained from the QCD partition function by replacing the matrix elements Dirac operator by (Gaussian) random numbers. The partition function (17) was also derived explicitly starting from a chiral random matrix model instead of only using symmetry arguments [26]. Another advantage of using a random matrix model is that one can easily perform numerical simulations. For example, the quenched spectral density was calculated numerically [26] and was found to be in agreement with (16).
Below we will show that the replica limit of the Toda lattice equation (16) can also be used to obtain the unquenched spectral density using a similar regularization of the partition function with bosonic quarks. In this case we need Z N f ,n=−1 ν with N f = 0. This involves an integral over the supergroupĜl(N f |2). This integral is not known explicitly for N f = 0. For that reason we will derive Z N f ,n=−1 ν directly from the corresponding random matrix model instead of the low-energy effective partition function based onĜl(N f |2).

III. RANDOM MATRIX MODEL FOR QCD AT NONZERO CHEMICAL POTENTIAL
In this section we define the random matrix model which we will use to calculate the generating functions. Random matrix models for QCD originally [8] focused on the quark mass dependence at zero chemical potential. Later [7], a random matrix model including the chemical potential successfully explained why quenched lattice QCD at zero temperature has a phase transition at a chemical potential of half the pion mass. However, a disadvantage of this model is that no eigenvalue representation is known which is required for the use of orthogonal polynomial methods. The random matrix model that will be used in the present paper was introduced in [34]. It differs from the model in [7] by a different form of the chemical potential term. Its partition function can be reduced to an eigenvalue representation and allows for an explicit solution in terms of orthogonal polynomials [34]. Because this model captures the correct global symmetries of the QCD Dirac operator, its low energy limit is also given by (14). A related model defined in terms of a joint eigenvalue distribution was introduced in [37]. Although this model is not in the universality class of QCD partition functions, the complex orthogonal polynomial methods that were developed for the derivation of the spectral density are applicable to the random matrix model introduced in [34]. In the section V we will comment on the relations between the different models.
The random matrix partition function with N f quark flavors of mass m f and n replica pairs of regular and conjugate quarks with masses y and z * each is defined by where the non-Hermitian Dirac operator is given by Negative numbers of flavors N f , n < 0 will denote insertions of bosons given by inverse powers of determinants. To allow for the presence of both bosonic and fermionic flavors we will distinguish between N b and N f in later formulas.
Here Φ and Ψ are complex (N + ν) × N matrices with the same Gaussian weight function The matrix D(µ) replaces the QCD Dirac operator plus chemical potential in (3). It has exactly ν zero modes which identifies ν as the absolute value of the topological charge. The model (20) has the same flavor symmetries as QCD which are broken by the quark masses and the chemical potential in exactly the same way as in QCD. The only difference with the model [7] is that the chemical potential is multiplied by a complex matrix Ψ with Gaussian distributed matrix elements. In the microscopic limit, where µ 2 N , m f N and zN are fixed as N → ∞ this partition function will match the partition function (14) after fixing the scale of the parameters appropriately. The latter is uniquely determined in terms of the flavor symmetries and their breaking, and thus universal. As was already done in the previous section, the microscopic limit of the partition function will be denoted by Z instead of the notation Z which we use both for the QCD partition function and the random matrix theory partition function. Generally, the overall normalization of Z and the microscopic limit of Z is different. The form of the Gaussian weight w G (Φ) respects the flavor symmetries, but these symmetries do not exclude traces of higher powers of Φ † Φ in the exponent as well as other non-invariant terms. At zero chemical potential it was shown [50,51,52,53,54,55,56,57] that partition functions and eigenvalue correlations in the microscopic limit are independent of the form of this weight, which, in random matrix theory, is known as universality. The only condition for the probability density is that the spectral density near the origin is nonvanishing in the thermodynamic limit. In the derivation of the joint eigenvalue density of (20) it is essential that the two weight functions are Gaussian. However, we believe that this is only a technical requirement and higher order invariant terms will not alter the microscopic limit of the joint probability distribution (for more discussion see section V).
In [34] it was shown that the partition function (20) has an eigenvalue representation after choosing an appropriate representation for the matrices Φ and Ψ yielding where the integration extends over the full complex plane and the joint probability distribution of the eigenvalues is given by The Vandermonde determinant is defined as and the weight function reads where K ν is a modified Bessel function. The parameters of the weight function, {m f }, y, z * and µ are collectively denoted by a. In the quenched case, N f = n = 0, the weight function will be denoted by w(z k , z * k ; µ). We have included a normalization factor 1/µ 2N in eq. (24). This way it reduces to the joint eigenvalue density of the chGUE times a product of delta-functions in the limit µ → 0. We expect that the terms that are nonvanishing in the microscopic limit are universal and are required for the derivation of the microscopic spectral density.
In the asymptotic limit, N |z k | 2 /µ 2 ≫ 1, the modified Bessel function in the weight is well approximated by its leading order asymptotic expansion. The corresponding joint eigenvalue distribution was introduced in [37] as a model for QCD at nonzero chemical potential. Let us look closer at the condition N |z k | 2 /µ 2 ≫ 1. In terms of the half-width, x max ∼ µ 2 , of the cloud of eigenvalues and the spacing, ∆ ∼ 1/N , of the imaginary part of the eigenvalues it can be rewritten as For x max > ∆ this implies that |z k | ≫ ∆ where we expect to find the weight function of the Gaussian Unitary Ensemble perturbed by an anti-Hermitian matrix which indeed has a Gaussian weight function (see [36]). In other words, the appearance of the modified Bessel function K ν is a direct consequence of the chiral symmetry of the problem.

IV. GENERATING FUNCTIONS FROM RANDOM MATRIX MODELS
The aim of this section is to derive explicit expressions for the generating functions introduced in section II using the random matrix model introduced in the previous section. In particular, we will present new expressions for the corresponding partition functions with any given number of fermions and one pair of conjugate bosonic quarks, as they are needed for the replica limit of the Toda lattice equation. We will find that the partition functions with bosonic replicas have to be regularized. This divergence will lead to a factorization into a partition function containing only fermions and a purely bosonic one. After taking the microscopic limit, they can be expressed entirely in terms of known partition functions (14) and (17). In the following section we consider partition functions with N f flavors and n fermionic replicas to show that the random matrix model partition functions satisfy the Toda lattice equation (15).

A. Orthogonal polynomials and Cauchy transforms
As is the case for Hermitian random matrix models, the partition functions (23) can be expressed in terms of orthogonal polynomials, and the integrals to obtain the spectral correlation function can be performed by means of the orthogonality relations. In this section we present the polynomials corresponding to the weight function (26) for The orthogonal polynomials are defined as solutions to the following orthogonality relation where the integral is over the full complex plane with measure d 2 z = dRez dImz, and r k denotes the squared norms of the polynomials. Since the quenched weight w(z, z * ; µ) appearing in this relation is real positive the polynomials form a complete set with positive squared norms r k . The polynomials p k (z) depend only on the variable z and not its complex conjugate, as is indicated by the notation. They are also functions of the chemical potential µ, the size N of the random matrix, as well as of the topological charge ν. We suppress the dependence on µ, N and ν.
In monic normalization, p k (z) = z 2k + . . ., the solution of (28) is given by [34] Since the Vandermonde determinant only contains squared variables only polynomials in those squared variables will appear. The same happens for the polynomials on the real line when mapping the Laguerre ensemble from the positive real numbers to the full real line. For the weight (26) the norms are given by Because the Laguerre polynomials L We note that h 0 (y) has a nontrivial y-dependence, unlike the polynomial p 0 (y). In contrast to random matrix models with eigenvalues on a curve in the complex plane, the integral over the pole in the complex plane is always well defined. However, if we take the limit µ → 0 of an anti-Hermitian Dirac operator we have to make sure that the argument y lies outside the support of eigenvalues. The leading order asymptotic large-y behavior of h k (y) is obtained by expanding 1/(y 2 − z 2 ) in a geometric series. Because of orthogonality all coefficients up to order 1/y 2k vanish resulting in the asymptotic behavior This result is useful for checking the identities derived below.

B. Partition functions with one flavor
The polynomials and their Cauchy transform enjoy the following relations to one flavor partition functions. The former is related to a single fermionic flavor partition function through the so-called Heine formula, Here we have written the fermion determinant as an expectation value with respect to the quenched partition function at N f = n = 0. (Vanishing replica (or flavor) indices will be commonly suppressed throughout the following.) For µ < 1 the partition function is manifestly real and positive for real masses z, as the orthogonal polynomials are Laguerre polynomials of argument −N z 2 /(1 − µ 2 ) with zeros on the imaginary axis. For µ > 1 the zeros are located on the real axis. Remarkably, the Yang-Lee zeros [58] of the model in [7] behave quite different from the model (20) discussed in this paper. The polynomials can also be interpreted as characteristic polynomials. The relation (33) can also be written down for p l =N (z), expressing it as an average over l variables. However, in order to be precise we would have to chose the weight function in (26) to be N -independent as in [38,39,40]. Since we are mainly interested in taking the large-N limit we prefer to explicitly keep the N or volume dependence inside the weight. For l ≈ N we can safely ignore this technical subtlety.
Turning to the Cauchy transform of p N (z), it is given by the partition function with one bosonic quark flavor [39,40], The same remarks about positivity and the N -dependence of the weight function made before also apply here.
Our program is to express spectral correlation functions in terms of partition functions. This has the advantage that the microscopic limit is given by a group integral based on the symmetries of the partition functions. To achieve our goal we also need all possible two flavor partition functions. At zero chemical potential the program of expressing spectral correlation functions in terms of multiflavor partition functions was carried out in [59,60].

C. Partition function with two flavors
It is known from Hermitian random matrix models [61] that the expectation value of two characteristic polynomials is proportional to the kernel of the orthogonal polynomials (28). As was shown in [38], for an arbitrary weight function, similar results can be derived for non-Hermitian random matrix models. Since we are now dealing with non-Hermitian matrices we have to distinguish between a characteristic polynomial (or mass term) and its Hermitian conjugate. For a single pair of one fermionic quark with mass z and one conjugate fermionic quark (with complex conjugated eigenvalues) with mass u * we have the relation [38] Here we have defined the bare kernel K N (z, u * ) which differs from the full kernel by the square root of the weight function at each argument. The full kernel appears in the computation of correlation functions from orthogonal polynomials (see section V). In particular, for z = u we obtain the result for one pair of fermionic replicas n = 1 This quantity is manifestly real and positive. On the other hand, if we compute the partition function with two flavors with the same eigenvalues, we obtain [38] For Hermitian random matrix models the two results (35) and (38) are related through the Christoffel-Darboux formula, which in general does not hold in the complex plane (see also Appendix A). We now turn to the most important two-flavor partition function, containing a pair of a bosonic quark and its conjugate. In the replica limit of the Toda lattice equation this partition function appears when computing the quenched density (see (16)). In the microscopic limit this partition function was calculated in [26], both by starting from a random matrix model, and by integrating explicitly over the Goldstone manifold in the chiral Lagrangian (17). It turned out that this partition function is singular and has to be regularized according to the prescription given in (19). The partition function is obtained in the limit of vanishing regulator, ǫ → 0, where it diverges as ∼ log(ǫ). This diverging constant can be absorbed in the normalization of the partition function. Below we will establish the nature of this singularity in the orthogonal polynomial approach. In fact, it does not enter in the calculation of the spectral density for the random matrix model (20) by means of complex orthogonal polynomials [34]. However, we will show that it occurs even at finite N in the random matrix model partition function with one bosonic quark and one conjugate bosonic quark of the same mass.
Generalizing the results in [39] to the chiral case we find that the following relation holds where Equation (40) (41) is logarithmically divergent, just like the divergence encountered in (17). In order to regularize this divergence we cut out a circle C(x, ǫ) of radius ǫ around the pole at z = x and one around the pole at z = −x.
To evaluate the singular part of the regularized Q ǫ (x * , x; µ) we expand the weight function around the singular point and keep only the leading part. For |x| ≫ ǫ we find for the pole at z = x We find the same contribution from the pole at z = −x. Thus we have obtained the following interesting relation for the singular part of the bosonic partition function We are not aware that such a relation between the weight function (26) of complex orthogonal polynomials and the bosonic partition function at a degenerate mass pair was previously known. It translates easily to non-chiral non-Hermitian random matrix models. Relation (43) shows that the singularity that was found in [26] for the microscopic limit of Z n=−1 N (x, x * ; µ) even persists for finite size random matrices and can be expressed simply in terms of the weight function. Only weak assumptions have been made in the derivation of (43), and its validity extends to a wide class of weight functions. For x ∼ O(ǫ) the above derivation does not apply. Instead we find that lim x→0 Q ǫ (x, x * ; µ) ∼ 1/ǫ 2 . As we will explain in the next section, the logarithmic singularity we have found for the two-point function (43) will also occur in more general partition functions with at least one pair of conjugate bosonic quarks.
The last two-point function we need below is the partition function with one fermion and one boson. This partition function reads [39] Here, we have defined the bare kernel N N −1 (x, y) consisting of orthogonal polynomials and Cauchy transforms. This partition function is nonsingular for any values of the arguments and correctly normalizes to unity at equal arguments y = x. The remaining two flavor partition functions combining (conjugate) fermionic quarks and (conjugate) bosonic quarks are given in Appendix A.

D. Partition Functions with an arbitrary number of flavors and one conjugate flavor
General expressions for the expectation value of ratios of characteristic polynomials valid for an arbitrary complex weight function were given in [38,39,40]. We will apply these results to partition functions with N f flavors plus one pair of regular and conjugate fermionic (n = +1) or bosonic (n = −1) flavors which both enter in the expression for the spectral density (16) after taking the replica limit n → 0.
The partition function with N f fermionic quarks with masses m f is given by [38] Next we add a single pair of a fermion and its conjugate. This partition function is given by [38] Z where K N +1 (x, y) is the bare kernel introduced in (36). In the replica limit of the Toda lattice equation, the partition function with one pair of conjugate bosons instead of conjugate fermions enters as well. It is given by [39] (see also Appendix A) = .
If we take the limit of degenerate bosonic masses, y → x, the function Q(x * , y; µ) inside the matrix element A N +N f −2 (x * , y) becomes singular, while all other matrix elements remain finite. Expanding the determinant with respect to the last row this singular part ∼ Q ǫ (x * , x; µ) gives the dominant contribution, Using the relations (40) and (46) we thus arrive at the following factorized result, Both sides are regularized by replacing Q(x, x * ; µ) by Q ǫ (x, x * ; µ) as discussed above (42). This is the main result of this section. Let us stress that the factorization we obtained is not a consequence of the suppression of nonplanar diagrams in the large N limit. In (50) it occurs at finite-N and is strictly due to the singularity for coinciding arguments. This completes the derivation of all generating functions required for the derivation of the spectral density from the replica limit of the Toda lattice equation (16). If we were to compute not only the spectral density but also higher order correlation functions from the Toda lattice approach (see e.g. in [26]), we would have to add more pairs of replicated bosons, each with different masses y j and x * j . Taking the degenerate limit y j → x j would lead to the same type of factorization into ∼ j Z nj =−1 N (x j , x * j ; µ) and a remaining fermionic partition function.

V. TODA LATTICE EQUATION FROM ORTHOGONAL POLYNOMIALS
The purpose of this section is threefold. First we show that the generating functions with n pairs of fermionic replicas satisfy the Toda lattice equation. Since the microscopic limit of the random matrix partition function, Z N f ,n N , is given by the effective partition function (14) [26,33], this result also shows that the Toda lattice equation (15) for N f = 0 [26] can be extended to arbitrary N f . Second, we show that the spectral density obtained from the replica limit of the Toda Lattice equation agrees with the result computed from orthogonal polynomials [34]. This extends the agreement between the two approaches to cases where the weight in the partition function is no longer positive definite. We close this section with some remarks about universality.
A. The spectral density from the Toda lattice equation We start by showing that the Toda lattice equation (15) can be modified to hold for random matrix model partition functions at finite N . The starting point of the evaluation is the expression for the partition function with N f + n ordinary flavors and n conjugate fermionic flavors generalizing (47). While this result is given in [38] for non-degenerate masses we have to take the limit where the replicated masses become degenerate. This leads to successive derivatives in these variables of the polynomials and of the kernels (36). For simplicity we only display the result for N f = 1, omitting constant factors. Here the derivatives are defined as We observe that this expression depends on the matrix size explicitly through the combination N + n and implicitly through factors that originate from the N -dependence in the exponent of the weight function. Using the Sylvester identity as in [26] while keeping only the N dependence that enters as N +n we easily derive the Toda lattice equation, again ignoring constant factors. This equation is valid at finite N provided that the factor N that appears explicitly in the exponent of the weight function is not changed. The subscript of the partition function thus denotes the total number of complex integration variables. For large N this distinction can be ignored. It is not difficult to generalize the Toda lattice equation to arbitrary N f . The result is We still have to fix the overall normalization constant of the partition functions which in principle depends on N, N f , as well as on the replica index n. Since the l.h.s. of (54) is linear in n we have to choose the normalization constants in the r.h.s. to retain this n-dependence. With the appropriately adjusted constants the random matrix model generating functions satisfy the equality (15).
Since the microscopic limit of the random matrix partition functions Z N f ,n N ({m f }, z, z * ; µ) is given by the group integrals (14), we conclude that (15) is satisfied for any number of flavors. Explicit expressions for these partition functions in terms of Bessel functions will be given in section VI where we also give results for the spectral densities.
Before doing so let us step back to discuss what we have achieved for the spectral density using the replica limit of the Toda lattice equation (54). Combining equations (15) and (12) with the factorization (50) and canceling one of the partition functions in the denominator of the Toda lattice equation, we arrive at the following result for the spectral density, where we have explicitly displayed the divergent normalization factor −(π/2) log ǫ which cancels the divergent factor that was obtained in (43). The ν dependence of the spectral density is contained implicitly in the partition functions. This is the central result of the Toda lattice approach applied to the random matrix model (20). The spectral density factorizes into QCD like partition functions which, at low energy, can be expressed as group integrals of the form (14) and (17).
Let us discuss some properties of the spectral density. First, it is symmetric under z → −z (together with z * → −z * ), because the partition function Z N f ,n=1 N contains only squared variables and Z n=−1 N is even as well. This is a consequence of the axial symmetry of the partition function which requires that the nonzero eigenvalues occur in pairs of opposite sign. Inserting the result (47) for Z N f ,n=1 N we see that the prefactor is canceled from the Vandermonde (for an explicit expression in terms of determinants see (57) below). Replacing the polynomials and kernels in Z N f ,n=1 N by their expressions in terms of partition functions we notice one further symmetry (for real quark masses), that is taking the complex conjugate and at the same time exchanging z and z * . These symmetries follow directly if we write the spectral density as an integral over the joint eigenvalue distribution. Thus the real part is symmetric about the real axis in the complex eigenvalue plane while the imaginary part is anti-symmetric. Therefore, the eigenvalue density on the real axis is real. However, away from the real axis it is, in general, neither real nor positive and therefore does not define a probability density. The lack of reality directly relates to the asymmetric way the variables z and z * enter in the polynomials and kernel respectively. The conjugation symmetry is related to the invariance of the partition function under changing the sign of µ, and is valid for QCD beyond the ergodic domain (2) as well.

B. Equivalence of the Toda lattice and orthogonal polynomials approach
In order to compare the result (55) from the replica limit of the Toda lattice equation to the result from orthogonal polynomials [34] it is instructive to rewrite (55) in terms of the weight function, polynomials and kernels. Reinserting (43), (47) and (46) we obtain again after normalizing appropriately. As we will explain next, this is precisely the result that was obtained in [34]. It is well-known from random matrix theory that the spectral density can be written in terms of a kernel of orthogonal polynomials [62]. This also holds for complex weight functions. However, if the weight function is not a symmetric function of z and z * , we have to use bi-orthogonal polynomials in the complex plane defined by [39,63], assuming the pseudo norms r (N f ) k to be nonvanishing. (Bi-orthogonal polynomials in general no longer form a scalar product with positive definite norms.) In contrast to orthogonal polynomials on the real line, the p k (z) and q k (z) are in general different polynomials. Defining their kernel by the following expression for the correlator of k complex eigenvalues holds [62] In particular, for the spectral density we obtain In [34] it was shown that the kernel for N f flavors can be expressed as a determinant of kernels for zero flavors which immediately leads to the result (57).
This establishes the equivalence between the Toda lattice approach and the bi-orthogonal polynomial approach at the level of the spectral density. Starting from the observation that also in the case of generating functions for multipoint correlation functions the partition function Z n=−1 N (z l , z * l ; µ) factors out, we are confident that the equivalence between the Toda lattice approach and the bi-orthogonal polynomial approach can be established in that case as well.

C. Universality
After having computed the spectral density from the replica limit of the Toda lattice equation and compared to that from the orthogonal polynomial approach, we would like to use the insight gained to discuss the issue of universality. First of all, we can distinguish between two different ways to approach universality. The first way is based on the observation that different theories with the same pattern of spontaneous symmetry breaking and a mass gap have the same low-energy limit. For example, gauge field theories and random matrix theories with the same global symmetries are described by the same effective partition function in the microscopic scaling limit. For QCD with three colors and fundamental quarks at nonzero chemical potential this universal partition function is given by the group integral (14). In QCD, a mass gap exists because of confinement, and in Random Matrix Theory, a mass gap appears in the large N limit. Therefore, if we change the action without affecting the global symmetries and the existence of a mass gap, the microscopic scaling limit will remain the same. A second way to understand universality is based on a direct calculation of correlation functions for a deformed probability distribution which does not necessarily have the unitary invariance of the Gaussian model. In the case of deformations that respect the unitary invariance, universality then follows from the asymptotic behavior of the modified orthogonal polynomials. In particular, for chiral random matrix theory at µ = 0 this approach has been very successful [51,55].
Let us also comment on the role of the weight function in non-Hermitian random matrix theories. From the replica limit of the Toda lattice equation it is clear that the eigenvalue density will be proportional to the partition function with one pair of conjugate bosonic flavors which, according to this work, is proportional to the weight function. In general, this weight function contains both universal and non-universal parts, and in the microscopic scaling limit only the universal parts survive. If we allow for deformations of the Gaussian weight by higher order powers the resulting weight function, and therefore the spectral density, will be modified. It is only in the microscopic scaling limit that the non-universal parts of the weight function disappear and a universal spectral density will be recovered. In Hermitian matrix models the weight function at finite N is also non-universal but reduces to a universal result in the microscopic scaling limit. Non-Hermitian random matrix models for QCD at µ = 0 represent a continuous one parameter deformation of the Hermitian model at µ = 0, governed by µ 2 N . In this case, the weight function contains an additional universal piece in the microscopic scaling limit that develops into a delta function ∼ δ(Im(z)) for µ 2 N → 0. The additional universal piece occurs because an eigenvalue representation can only be obtained after a nontrivial integration over the similarity transformations that diagonalize the non-Hermitian matrices. In the Hermitian case it is trivial to obtain the eigenvalue representation, and no additional universal piece is generated. For technical reasons we have not been able to integrate out the similarity transformations when higher powers are added to the Gaussian probability distribution of the model (20). We mention in passing that for non-chiral non-Hermitian random matrix models a universality proof by deforming the random matrix potential could be given [64].

VI. RESULTS AND EXPLICIT EXAMPLES
A. The partition functions in the microscopic limit Having obtained analytical expressions for the partition functions of the random matrix model (20) at finite N we now consider the large N limit. This limit will be taken according to the counting scheme where m f N and µ 2 N are kept fixed as N → ∞ and is referred to as the microscopic limit. The rescaling of the non-Hermiticity parameter µ 2 with N was first introduced in [65] as the concept of weak non-Hermiticity. As discussed in the previous section, in the microscopic limit the random matrix partition function is given by an integral over the Goldstone manifold, i.e. it is uniquely determined by the spontaneous breaking of the flavor symmetries. The random matrix model can be viewed as an alternative way to calculate the integral over the Goldstone manifold, which in the case of nonzero chemical potential and both bosonic and fermionic quarks has not been accomplished in any other way.
The dependence of the quark mass and the chemical potential in the integral over the Goldstone manifold is only through the combinations mΣV and µ 2 F 2 π V , cf. (14) and (17) [14]. The dimensionful scales can be recovered from the random matrix model by making the replacements To make contact with the physical content of the theory we will express our results in terms of the dimensionful scales in the remainder of this section.
In the case we have only either fermionic or bosonic quarks and no conjugate quarks, the Goldstone bosons are not charged with respect to the chemical potential, and we should find the microscopic partition functions at zero chemical potential. One easily observes that the effective partition function (14) does not depend on the chemical potential in this case. In the microscopic limit of (33) we thus find the partition function (14) with one fermionic flavor [10,66,67,68] An overall scaling factor e −µ 2 F 2 π V has been removed. This unphysical factor can be removed in the random matrix model by scaling the Gaussian potential in the weight (22) by a factor of 1 − µ 2 [34]. We did not include an explicit scale factor in the weight function, but will assume this procedure has been done in order to arrive at the correct expressions in the microscopic limit.
When there is one pair of conjugate fermionic quarks, the microscopic partition function is given by the result obtained by performing the integral (14) over the Goldstone manifold [26] Z n=1 For the universal partition function with a pair of conjugate bosonic quarks we find using (43) where we have divided out the divergent factor −π log(ǫ)/2. This is again in agreement with the result [26] obtained by an explicit integration over the Goldstone manifold (17). All other microscopic partition functions required to calculate the spectral density can be obtained from the relation (50) B. The microscopic spectral density In this section we write out several different examples for the eigenvalue density of the QCD Dirac operator in the microscopic limit. We emphasize that these results are nonperturbative analytical predictions for the QCD Dirac spectrum at nonzero chemical potential which should be reproduced by lattice QCD simulations. In order to get more compact expressions we will express our results in terms of the partition function (64). The normalization of the quenched density is chosen such that there are ΣV /π eigenvalues per unit length along the imaginary axis. The density of the projection of the eigenvalues on the imaginary axis is therefore equal to the eigenvalue density at µ = 0. The other densities reduce to the quenched density when the quark masses are taken to infinity (m f ΣV ≫ µ 2 F 2 π V ) and we use this limit to fix the normalization.
The general expression for the spectral density is given by where the microscopic limit of Z n=−1 N (z, z * ; µ) is given in the previous subsection. The microscopic limit of Z where x k = m k V Σ. Note that with only regular fermionic quarks, µ has no effect on the microscopic partition function. The partition function Z N f ,n=1 N ({m f }, z, z * ; µ) is given in terms of the polynomials p N (m) and the bare kernel K N +1 (x, y) in (47). The microscopic limit of this partition function is obtained by expressing the polynomials and the kernel in terms of partition functions that are known in the microscopic limit. The polynomials p N (m) can be expressed as the partition functions with one fermionic flavor (see (33)), which in the microscopic limit does not depend on µ and is given by I ν (mV Σ) in (63). In (37) the kernel K N +1 (x, y) was shown to be equal to Z n=1 ν which is given in (64). In the remaining sections we will set ν = 0 for simplicity.

The quenched spectral density
The quenched microscopic eigenvalue density of the Dirac operator at nonzero baryon chemical potential is given by [26] This result has been checked by direct numerical simulations [26] within a random matrix model, and successfully compared to results from lattice QCD [43]. Plots of the quenched spectral density are given in Fig. 1. By using the leading order asymptotic forms of the Bessel functions, one easily shows that in the limit µ 2 F 2 π V ≫ 1 the support of the quenched spectrum is inside the strip [26,49] |Re(z)|Σ 2µ 2 F 2 π < 1 (70) with a density that is equal to In Fig. 1 we observe that this plateau is already present for µF π √ V = 2.5. This step in the spectral density is associated with a phase transition of the generating function. The theory with a pair of conjugate fermionic flavors of mass x, described by Z n=1 (x, x; µ), is in the normal phase for xΣ/F 2 π = m 2 π /2 > 2µ 2 and in a Bose condensed phase for xΣ/F 2 π = m 2 π /2 < 2µ 2 [14]. For the x-integrated quenched spectral density we find (z = x + iy) consistent with our choice of normalization.

The spectral density for one flavor
The microscopic eigenvalue density with one dynamical quark of mass m and with baryon chemical potential µ is given by [34] This result is shown in Fig. 2 for µF π √ V = 0.1 and mV Σ = 10, and in Fig. 3 for µF π √ V = 2.5 and mV Σ = 5. It is our first example of a spectral density that is complex due to the sign problem.
Since the unquenched spectral density is proportional to the quenched spectral density its support cannot go beyond the support of the quenched spectral density. In particular for µ 2 F 2 π V ≫ 1 the unquenched spectral density is also confined inside the strip given by (70). As we can observe from Fig. 4 the asymptotic behavior is much more complicated than in the quenched case. An asymptotic form of the integral (64) was derived in [37]. For µ 2 F 2 π V ≫ 1 and (x + m)Σ/(4µ 2 F 2 π ) < 1, it is well approximated by This approximation breaks down if (x + m)Σ/(4µ 2 F 2 π ) > 1 because then the saddle point of the t-integration will be outside the interval [0, 1]. The exponential functions in the asymptotic behavior of the Bessel functions no longer compensate each other in the additional contribution to the unquenched spectra density. Instead of a plateau for N f = 0 we find an oscillatory contribution with an amplitude that increases exponentially with the volume. For the total spectral density a plateau is still visible in the region where the quenched contribution dominates. As an illustration we show in Fig. 4 the real part of the eigenvalue density for µ 2 F 2 π V = 100 and mV Σ = 100. The period of the oscillations is of the order of the level spacing at µ = 0 whereas the amplitude is of the order exp(10) for our choice of parameters. With these values the real part of the density deviates substantially from the quenched density. Note that the real part of the density changes sign at z = m. In this case the scale of the imaginary part is comparable to that of the real part.
(75) The severe sign problem manifest itself in the strongly oscillating region of the eigenvalue density which has amplitudes of the order exp (10). The full size of the peaks has been clipped for better illustration.
In Fig. 5 we show the real and imaginary parts of the spectral density for µF π √ V = 2.5 and m 1 ΣV = m 2 ΣV = 5. The spectral density is neither real nor positive.

The spectral density with one pair of conjugate flavors
The microscopic spectral density of the Dirac operator in QCD with one pair of conjugate quarks (with mass m and m * ) is found to be (see (87) and [38]) From this expression it follows that the density is real, positive, and has a zero at z = m. This is illustrated by Fig.  6, where we plot the spectral density for mV Σ = 5 and µF π √ V = 2.5. Since the quark mass, m, is real, the determinant of the conjugate quark is identical to the determinant of a quark with a chemical potential of the same magnitude but with opposite sign, cf. (13). That is, the partition function in this case is equivalent to a theory with two degenerate flavors and nonzero isospin chemical potential, µ I . This version of QCD does not have a sign problem [69] as has been explored by lattice QCD simulations [46]. 5. The density with one ordinary flavor and one pair of conjugate flavors As a final example we write out the microscopic eigenvalue density in QCD with 3 light flavors with real masses m u = m d = m and m s and chemical potentials µ I = µ s = µ. This theory has a sign problem due to the extra quark that isn't matched by a conjugate quark. The phase diagram as a function of µ I and µ s was determined in [70]. The eigenvalue density is (88) For z = m the two top rows are identical and the first and third column are identical in the 3 × 3 determinant. Hence the density is zero at z = m but does not change sign. While for z = m s the first and second column in the 3 × 3 determinant are identical, so the density changes sign at z = m s . A plot of this density is shown in Fig. 7.  . Note that the density bounces off the real axis at z = mu and changes sign at z = ms as expected.

VII. CONCLUSIONS
We have analyzed the spectrum of the QCD Dirac operator at nonzero baryon chemical potential. In the microscopic limit this spectrum is uniquely determined by the global symmetries of the QCD partition function. This is true both for the quenched and the unquenched theory. In both cases the spectral density of the Dirac operator can only be obtained after the introduction of a complex conjugate pair of valence quarks resulting in a nontrivial baryon charge matrix and Goldstone modes with nonzero baryon number. The microscopic limit of the generating function for the QCD Dirac spectrum can therefore be represented in two different ways. First, as an integral over the Goldstone manifold which is determined by the pattern of symmetry breaking. Second, as the large-N limit of a partition function of an ensemble of non-Hermitian matrices with the symmetries of the QCD partition function. In recent work, the Dirac spectrum was derived directly from the random matrix model by means of complex orthogonal polynomials. In this paper we have obtained the Dirac spectrum from the replica limit of the Toda lattice equation. This explains that the spectral density factorizes into the product of a fermionic partition function and the partition function for one conjugate pair of bosonic valence quarks. The bosonic factor does not depend on the dynamical quarks and therefore goes beyond what is implied by the Toda lattice structure. This arises as a consequence of a singularity in the partition functions that contain a pair of conjugate bosonic quarks with the same mass. This singularity persists even for finite size matrices. In fact, it is present for an ensemble of 2 × 2 random matrices. Also the Toda lattice structure is valid for finite size matrices, and one can easily verify that the correct spectral density is obtained for an ensemble of 2 × 2 matrices.
As a new result, we have also obtained explicit expressions for the spectral density when the absolute value of the fermion determinant is present in the QCD partition function. In this case a direct check from lattice simulations is possible for all values of the scaled chemical potential, as has already been done in the fully quenched case.
The calculation of the unquenched microscopic spectral density has been a highly nontrivial test of the Toda lattice approach. We expect that it will be equally successful for all remaining β = 2 cases. For example, the method could be applied to the calculation of GUE S-matrix fluctuations and parametric correlations [71] as well as to other non-Hermitian random matrix ensembles [72]. The extension of this approach to Dyson index β = 1 and β = 4 is still an open problem.
The spectral density can be obtained from the replica limit of a Painlevé equation that has been derived from fermionic partition functions only. On the other hand, if a family of partition functions satisfies a Painlevé equation, a Toda lattice equation can be derived from the Bäcklund transformations of the Painlevé system. Therefore, all systems in the same universality class have the same microscopic limit for the fermionic partition functions, the bosonic partition functions and the microscopic spectral density.
Another important difference between Hermitian and non-Hermitian Gaussian random matrix theories is that the probability density of non-Hermitian theories do not factorize into the joint eigenvalue density and the distribution of the eigenvectors. The joint eigenvalue density is only obtained after a nontrivial integration over the similarity transformations that diagonalize the non-Hermitian matrices. This results in the appearance of a nontrivial universal factor in the joint eigenvalue distribution, in our case a modified K-Bessel function. The appearance of this function is a direct consequence of the chiral symmetry of the problem.
The microscopic spectral density derived in this paper is valid for small temperatures and all values of the dimensionless parameters µ 2 F 2 π V and mΣV . Within this range of parameters, where the sign problem develops and becomes severe, the effect on the spectral density can be followed analytically. A first manifestation of the sign problem is that the spectral density in the complex eigenvalue plane is no longer real. This occurs already in a parameter range where lattice simulations might be feasible. For larger values of µ 2 F 2 π V the analytical expressions for the eigenvalue density can serve as tests for numerical methods that apply to the nonperturbative regime of QCD at nonzero baryon chemical potential. In the limit where µ 2 F 2 π V ≫ 1 the spectral "density" shows oscillations with a period of the order of the level spacing and an amplitude that diverges exponentially with µ 2 F 2 π V . This behavior, which is not present in the quenched case, should be responsible for the breaking of chiral symmetry. This issue will be addressed in a future publication.
denominator can be written solely as determinants of polynomials and Cauchy transforms. We find and All remaining combinations of bosons and fermions and their conjugates can be obtained from the given ones by complex conjugation.
Next we wish to comment on the relation between these formulas and in particular on the additional pole terms which occur in the two kernels containing Cauchy transforms, equations (40) and (45). For random matrix models with real eigenvalues the kernel of the polynomials (36) satisfies the Christoffel-Darboux identity, For orthogonal polynomials in the complex plane (28) this relation is generally not satisfied so that (35) and (38) are different. In the Hermitian limit µ → 0 they are equal, due to the relation (80). The Eqs. (39) and (79) or (44) and (78) can be related if the following Christoffel-Darboux identities are valid: This is the case for Hermitian random matrix models, but in general there are no such Christoffel-Darboux identities for non-Hermitian random matrix theories. This identity follows from the observation that in general the orthogonal polynomials and their Cauchy transforms obey the same three-step recursion relation. Only the recursion involving the lowest Cauchy transform h k=0 (y) is different from that involving the polynomial p k=0 (y), which results in the extra integral Q(s, t) over poles in (40), or the single pole 1/(s 2 − t 2 ) in (45), respectively. Thus in the Hermitian limit (39) and (79), and (44) and (78) become equal due to their Christoffel-Darboux identities (81). This explains the presence of pole terms in (39) and (44). In (45) this term can also be understood differently, as it ensures the correct normalization of the partition functions (44) to unity at equal arguments. Let us now turn to more general partition functions. We will briefly explain here how our generating function (48) follows from the results of [39] (generalized to the chiral ensemble). There the following expectation value of complex eigenvalues is computed as In order to simplify the denominator the determinant can be expanded with respect to the next to last row. The remaining determinant only contains polynomials p k (m i ) and the kernel N N f −2 (m i , y). The invariance of determinants under adding and subtracting rows and columns can be used to first eliminate all sums in the kernels N N f −2 (m i , y). Then all polynomials can be reduced to monomials, leading to the following Vandermonde like determinant