The Complex Laguerre Symplectic Ensemble of Non-Hermitian Matrices

We solve the complex extension of the chiral Gaussian Symplectic Ensemble, defined as a Gaussian two-matrix model of chiral non-Hermitian quaternion real matrices. This leads to the appearance of Laguerre polynomials in the complex plane and we prove their orthogonality. Alternatively, a complex eigenvalue representation of this ensemble is given for general weight functions. All k-point correlation functions of complex eigenvalues are given in terms of the corresponding skew orthogonal polynomials in the complex plane for finite-N, where N is the matrix size or number of eigenvalues, respectively. We also allow for an arbitrary number of complex conjugate pairs of characteristic polynomials in the weight function, corresponding to massive quark flavours in applications to field theory. Explicit expressions are given in the large-N limit at both weak and strong non-Hermiticity for the weight of the Gaussian two-matrix model. This model can be mapped to the complex Dirac operator spectrum with non-vanishing chemical potential. It belongs to the symmetry class of either the adjoint representation or two colours in the fundamental representation using staggered lattice fermions.


Introduction
The first complex random matrix models date back to Ginibre [1] where the three classical Wigner-Dyson ensembles were generalised. In recent years we have seen a revival of complex matrix models, due to new applications as well as new insights and solutions. Todays interest in such models ranges from Quantum Chromodynamics (QCD) with chemical potential [2], fractional Quantum-Hall effect [3], Coulomb plasma [4] and growth processes [5] to scattering in Quantum chaos reviewed in [6]. New insights include the discovery of the regime of weak non-Hermiticity [7], interpolating between correlations of real eigenvalues of the Wigner-Dyson ensembles and those correlations computed by Ginibre. Furthermore, we know that there are many more complex matrix model symmetry classes than those with real eigenvalues [8].
The aim of this work is to solve such a new complex matrix model, by generalising a chiral version of one of the Wigner-Dyson ensembles. Although this is motivated by a specific application to QCD we hope that our result will be useful in other areas as well, because of the wide range of applicability of matrix models. A major problem in QCD is that the Dirac operator becomes non-Hermitian when introducing a chemical potential for the quarks. In general this spoils the positivity of the action and thus its interpretation as a Boltzmann weight when performing numerical simulations of QCD on an Euclidean space-time lattice.
A way out of this dilemma is to study properties of field theories similar to QCD, where Dirac eigenvalues come in complex conjugate pairs. There are two symmetry classes different from QCD, with the Dirac operator having an orthogonal or symplectic symmetry compared to unitary. This feature holds at both zero [9] and non-zero chemical potential [10]. In this work we will construct a complex matrix model with symplectic symmetry and solve it analytically. It has the virtue to be testable in lattice simulations including dynamical fermions. These are absent in the so-called quenched approximation, suppressing the Dirac operator determinants, and such matrix model predictions [11,12] have already been confirmed in quenched lattice simulations of QCD [13] and the symplectic symmetry class [14], where preliminary results of this work were used. Let us emphasise that unquenched matrix model predictions for QCD with chemical potential do exist [15,16].
The most important ingredient of matrix models is the way they incorporate the repulsion between eigenvalues or levels. The form of the repulsion term is dictated by the Jacobian when changing from matrices to eigenvalues. For real eigenvalues it is given by the absolute value of a Vandermonde determinant of the eigenvalues to the power of the Dyson index β = 1, 2, or 4, labelling the number of independent matrix elements in the Orthogonal, Unitary or Symplectic Ensembles respectively (OE, UE or SE). The same is true for the squared singular values of the so-called chiral ensembles (chOE, chUE or chSE) which have an additional block structure. When employing the method of orthogonal polynomials [17] all eigenvalue correlations of the corresponding Gaussian ensembles (GUE etc.) can be obtained using Hermite polynomials, or Laguerre polynomials in the chiral classes.
What is known about the level repulsion of complex eigenvalues and the corresponding correlations functions? For the complex GUE one obtains again the modulus squared of the Vandermonde determinant of complex eigenvalues [1]. For the corresponding chiral ensemble it was conjectured in [11] to have the same Jacobian of squared complex eigenvalues, which was later obtained from a two-matrix model realisation [15] 1 . While for β = 4 a Vandermonde to the power 4 can be constructed using normal matrices [18,5], the Jacobian of the complex GSE of generic non-Hermitian quaternion real matrices is more complicated [1]. It shows an additional level repulsion from the real line as pointed out in [19] (see also [20]), where this ensemble was solved for an arbitrary weight function. A similar feature was observed numerically in [10], where the level repulsion for all three complex chiral ensembles were compared.
Our task here is to derive the corresponding Jacobian for the complex chGSE using a two-matrix model of non-Hermitian quaternion real matrices, in analogy to the approach [15]. Because of the QCD application we allow for exact zero-eigenvalues, using rectangular matrices. The result for the Jacobian turns out to be again that of the non-chiral complex GSE of squared variables (generalised to rectangular matrices), hence having an additional level repulsion from the imaginary axis. We solve our model introducing orthogonal polynomials in the complex plane, where we encounter complex Laguerre instead of Hermite polynomials [3,21,19] in non-chiral ensembles. Other techniques include supersymmetry, e.g. in [22,21,20], the replica method [23,12] or a fermionic mapping [24] which we do not discuss here. Our computation completes the task of solving the complex extensions of the chiral and non-chiral Wigner-Dyson β = 2 and 4 ensembles. The difficulty with the remaining β = 1 ensembles is that they have a finite fraction of real eigenvalues, with the corresponding probabilities to be determined. Only the spectral density at weak non-Hermiticity of the complex GOE is know so far [22].
The article is organised as follows. In the following section 2 the complex chSE is defined for a general weight function, and all correlation functions are given at finite-N in terms of the corresponding skew orthogonal polynomials. The subsections 2.1 and 2.2 distinguish between weights without and with explicit insertions of characteristic polynomials or mass terms respectively, where we have to impose a two-fold mass degeneracy. In subsection 2.3 two examples for weight functions and their corresponding skew-orthogonal polynomials are given, a Ginibre type weight and a Gaussian combined with a K-Bessel function. The latter originates from the Gaussian two-matrix representation of our model discussed in detail in section 5. The orthogonal polynomials of the K-Bessel function weight are Laguerre polynomials in the complex plane, and we prove both their orthogonality and skew orthogonality, with the help of appendix A 2 . In section 3 we take the large-N limit at weak non-Hermiticity for the ensemble of complex skew orthogonal Laguerre polynomials. Explicit examples for the quenched and massive spectral density as well as for massive partition functions are given in subsections 3.1 and 3.2 . Section 4 is devoted to the strong non-Hermiticity limit, giving the same set of examples. Technical details for the derivation of the limiting kernel at strong non-Hermiticity are collected in the appendix B. In section 5 we relate the Gaussian two-matrix model representation of our complex model to the Dirac spectrum with chemical potential. We also provide a matrix representation of characteristic polynomials. The derivation of the Jacobian is deferred to appendix C. Our findings are summarised in section 6.

Correlation functions of the complex chSE at finite-N
In this section we define the complex extension of the chSE for a general class of weight functions. The diagonalisation of the matrix model representation in terms of two non-Hermitian quaternion real matrices is discussed in more detail in section 5 and appendix C. This Gaussian matrix model leads to a specific weight function containing K-Bessel functions. The reason is the nontrivial decoupling of eigenvalues and eigenvectors which is not present for ensembles with real eigenvalues. The two-matrix model representation determines the Jacobian we write down below as the interaction term for the eigenvalues.
For the calculation of the eigenvalue correlation functions in terms of skew orthogonal polynomials we proceed in two steps. We first treat the so-called quenched case where no explicit characteristic polynomials (determinants) are present in the weight function. Although the expressions we give 2 The Gaussian weight first introduced in [11] for complex Laguerre polynomials L ν N agrees only asymptotically or at special values of ν with the correct K-Bessel weight. The latter was introduced in [15] and conjectured to satisfy orthogonality.
for the complex eigenvalue correlations in the first step, eq. (2.7) below, are true for general weight functions, it is easier to include determinants as mass terms in a second step. Namely we can reexpress the corresponding massive correlation functions in terms of the quenched ones, following the same logic as in [25,26].

Complex eigenvalue correlation functions
The matrix model partition function of the complex Gaussian chSE is defined as the following matrix integral where Φ and Ψ are two matrices with quaternion real elements of rectangular size (N + ν)× N without further symmetry properties. In order to transform to eigenvalues we parametrise the following linear combinations as being equivalent to independent Schur decompositions of CD and DC. U and V are symplectic matrices and the diagonal matrices X and Y contain the complex conjugate eigenvalue pairs (x k , x * k ) and (y k , y * k ), respectively. After integrating out the second set we obtain in terms of the variables the following expression for the partition function For more details we refer to section 5 and appendix C for the Jacobian. The integrals in z k ∈ C extend over the whole complex plane with d 2 z ≡ dℜe(z)dℑm(z). Here w(z, z * ) denotes a general positive definite weight function defined in complex plane. It depends on both z and z * . The only restriction we want to impose here is the existence of all positive moments, d 2 z w(z, z * ) z 2k < ∞, k ∈ N , and similarly for the complex conjugated ones. For the Gaussian matrix model eq. (2.1) the weight is given by (see section 5) a second example is given below in eq. (2.27) (see also [27]). The general weight w(z, z * ) may depend on extra parameters such as the non-Hermiticity parameter µ > 0, or masses m f as we will discuss later. The former will allow us to take limits to either real eigenvalues by letting µ → 0, or to maximal non-Hermiticity at µ → 1, given for example by the chiral Ginibre type weight |z| exp[−zz * ]. Eq. (2.4) can be derived in a different way, inserting squared arguments into the Jacobian of the non-chiral, complex SE [1]. This is in analogy to the construction of the chiral from the non-chiral complex UE [11]. Below we compute all eigenvalue correlations in terms of skew orthogonal polynomials in the complex plane. This generalises the construction [19] using complex skew orthogonal polynomials for the non-chiral complex SE. There, complex Hermite polynomials were used for the special case of a Gaussian potential. In our chiral model with the special weight of the Gaussian two-matrix model we encounter complex Laguerre polynomials, as expected from the complex chUE (see [11,15]).
For a general weight function the k-point correlation functions of complex eigenvalues can be defined in the standard way [17] Being real positive functions they depend also on the complex conjugated variables z * 1 , . . . , z * k . We suppress this dependence in the following for simplicity. In analogy to the non-chiral model [19] the eigenvalue correlations show a depletion of the eigenvalues from the real axis, being identically zero there. This property was found also numerically for a chiral complex one-matrix model with chemical potential [10]. Since our Jacobian or interaction term is given in terms of squared variables compared to [19], there is the same depletion also from the imaginary axis, due to . Such a depletion has been already detected in lattice data for the corresponding symmetric class [14].
Following [17,19] all k-point correlation functions are determined through a quaternion determinant [17] of a matrix valued kernel, where the kernel is given by . (2.8) Here we have introduced the pre-kernel κ N (z 1 , z * 2 ), 3 Examples for such correlation functions are the spectral density given by 10) and the two-point function The pre-kernel eq. (2.9) is expressed in terms of skew-orthogonal polynomials q k (z). These are polynomials of degree k in the complex variable z ∈ C squared, q k (z) = z 2k + O(z 2k−2 ), where we have chosen the monic normalisation (because of the form of eq. (2.4) we only need polynomials of squared arguments, see also [11]). The skew-orthogonal polynomials are defined with respect to the following antisymmetric scalar product by obeying For a general weight they enjoy the following representation which holds in complete analogy to [19]. Here c N is an arbitrary constant that can be set to zero. Due to this representation we can observe that the polynomials have real coefficients, as q 2N (z) * = q 2N (z * ) and q 2N +1 (z) * = q 2N +1 (z * ). A similar representation of the skew orthogonal polynomials also holds for the SE with real eigenvalues [28]. In section 5 we will also give a matrix representation of the q 2N (z) in terms of a characteristic polynomial. The expectation values are to be taken with respect to the partition function eq. (2.4) (not to be confused with the antisymmetric scalar product). The skew orthogonal polynomials q 2N (z) and q 2N +1 (z) both depend only on the variable z. While eq. (2.7) for the k-point correlation functions holds for general weight functions it is convenient to explicitly split off determinants or mass terms from the weight. We shall do this in the next subsection, allowing us to compute characteristic polynomials more general than in the example eq. (2.14) above.

Massive complex correlation functions
The massive matrix model partition function of the complex chSE we wish to compute is given by where we allow for complex masses m 2 f = m * 2 f is this subsection. Compared to the partition function eq. (2.4) we have explicitly introduced mass terms in complex conjugated pairs as extra sources. The new weight function including the masses is still positive definite 4 . The mass terms can also be written as a determinant of quaternion real non-Hermitian matrices as we will show in section 5. Therefore these partition functions can be interpreted as characteristic polynomials, and we will give explicit formulas for them below.
In the presence of a non-Hermiticity parameter µ we can make contact with the chSE by taking the Hermitian limit µ → 0. This leads to 4-fold degenerate (real) masses, as indicated by the superscript 5 . Our aim is to express the partition function as well as all complex eigenvalue correlation functions in the presence of mass terms through the eigenvalue correlations without mass terms at N f = 0. For that reason we need to impose such a fourfold mass degeneracy as in [25]. There, massive correlation functions were expressed through quenched ones for Hermitian self-dual matrices [25], and we will follow the same strategy here. Similar results were also obtained for complex non-Hermitian matrices in the non-chiral case β = 2 [26], where only a twofold degeneracy is needed due to the different Jacobian k>l |z 2 k − z 2 l | 2 . To start we compute the massive partition function eq. (2.16) in terms of quenched correlation functions. When writing the mass terms as z 2 j + m 2 f = z 2 j − (im f ) 2 they can be absorbed into a larger Jacobian of N + N f eigenvalues z 1 , . . . , z N , im 1 , . . . , im N f , after multiplying with constant factors as The massive partition function is thus proportional to a quenched N f -point correlation function at i times the complex masses. Filling in all explicit factors from the definition eq. (2.6) we thus obtain (2.17) Dividing by Z N +N f and switching back to parameters u f = im f this gives an explicit expression for products of L pairs of characteristic polynomials times their complex conjugates It extends the corresponding results [27] for the corresponding β = 2 class. Next we use the same trick to compute the massive eigenvalue correlation functions, obtaining Inserting the above relation (2.17) for the massive partition function all pre-factors cancel and we obtain the very elegant expression . (2.20) The same relation was obtained previously for real eigenvalues in the chSE [25] and for the non-chiral complex UE [26]. It is easy to see that it also holds for the non-chiral complex SE, replacing squared by un-squared variables everywhere. Together with eqs. (2.7) and (2.8) this gives all eigenvalue correlation functions for our massive model eq. (2.16). Due to the different numbers of eigenvalues occurring, N compared to N +N f , these relations are only exact for an N -independent weight function w(z, z * ). When we consider N -dependent weights as in the examples eqs. (2.27) and (2.21) below our results will be still valid for asymptotically large-N as then N ≈ N + N f in the exponent. We have to add a word of caution to applying eqs. (2.20) and (2.17). While these expressions are true for general complex masses m f the limit of real masses has to be taken with care. As we have explained after eq. (2.6) our model shows a depletion of eigenvalues both from the real and imaginary axis. Therefore both the numerator and the denominator in eq. (2.20) vanish when setting the masses m f to be real or purely imaginary. We thus have to take the limit of the imaginary part of m f going to zero using the rule of de l'Hôspital. This leads to a finite expression for the massive correlation functions. We will illustrate this explicitly in examples in sect. 3 and 4 when taking the large-N limit.

Examples for skew orthogonal polynomials in the complex plane
In order to illustrate the finite-N expressions for the eigenvalue correlation functions we give two explicit examples of weight functions together with their corresponding skew orthogonal polynomials eqs. (2.14) and (2.15) as they appear inside the pre-kernel. We will use some results for non-skew orthogonal polynomials with respect to the same weight functions which are derived in appendix A.
The first example is from the Gaussian two-matrix model representation of the complex chSE (see sect. 5), Here µ is a non-Hermiticity parameter. Compared to the QCD symmetry class [15] the parameter ν measuring the rectangularity of the matrices is shifted, ν → 2ν, as for the real eigenvalue model. The real limit can be taken by letting µ → 0, leading to a delta-function in the imaginary part. This will be discussed in more detail in sect. 3 below. We only mention here that in this Hermitian limit the last term of the Jacobian in eq. (2.4) will provide two more powers of the real part, (ℜe(z)) 2 , arriving at the standard power of zero-modes |ℜe(z)| 4ν+3 for the real chSE (see eq. (5.10)). Let us consider a special case in which the K-Bessel function simplifies. For indices ν = ± 1 4 it exactly coincides [29] with its asymptotic value, K 2ν (x) ∼ e −x π 2x , and we obtain up to a constant This is of the same form as the weight in the corresponding non-chiral ensemble [19]. More generally speaking the K-Bessel function at half-integer index simplifies to an exponential times a finite sum over half-integer powers of its argument. In the non-chiral complex SE [19] complex Hermite polynomials were used to construct the skew orthogonal polynomials eq. (2.13) for a Gaussian weight function as eq. (2.22). This is consistent with the relation between odd or even Hermite polynomials and Laguerre polynomials L ± 1 2 N , respectively (see e.g. [29]). Following the same lines as in [19] we construct our skew orthogonal polynomials here using complex Laguerre polynomials, which are show to be orthogonal with respect to the weight eq. (2.21) in appendix A. We obtain in monic normalisation q 2k+1 (z) = (z 2 ) 2k+1 + . . . and q 2k (z) = (z 2 ) 2k + . . . . The fact that this choice satisfies eq. (2.13) can be seen as follows. In the case where both indices in eq. (2.13) are equal, the second line for k = l, the scalar product vanishes trivially due to antisymmetry. If both indices are odd we now show that the scalar product vanishes because of the orthogonality of the Laguerre polynomials from appendix A. To see that we use the three-step recursion relation 25) in the squared variable Z 2 = N z 2 1−µ 2 . Consequently the factor (z * 2 − z 2 ) in the scalar product eq. (2.12) can only lower or raise the index of the Laguerre polynomials by one unit. The terms where indices are raised or lowered vanish from orthogonality, having different parity. The middle term in the recursion eq. (2.25) only contributes at k = l and the corresponding terms cancel (as we have seen already from antisymmetry above). The same logic applies if both indices in eq. (2.13) are even. Only those terms in the recursion where the index is unchanged contribute, and they also cancel after using orthogonality, for all values of k and l.
It remains to analyse the scalar product in eq. (2.12) with one index even and one odd, q 2k+1 , q 2l S . For k > l the expression again vanishes trivially due to orthogonality, using eq. (2.25). So far the coefficients in the linear combination of even Laguerre polynomials in the definition (2.24) could have been arbitrary. It is the requirement q 2k+1 , q 2l S = r k δ kl for k ≤ l that fixes them uniquely in a given normalisation. One can easily check that the coefficients given in the definition (2.24) imply vanishing for k < l, upon using the recursion and the norms of the Laguerre polynomials eq. (A. 19). The overall constants in eqs. (2.23) and (2.24) are fixed by choosing the monic normalisation. This choice allows us to take the limit µ → 1 of maximal non-Hermiticity on the polynomials. The remaining coefficients r K k in eq. (2.13) at k = l follow to be This determines all quantities in the pre-kernel eq. (2.9) and thus all k-point correlation functions through eq. (2.7). We now turn to the second example given by a Gaussian, chiral Ginibre type weight In can be obtained from the above weight eq. (2.21) at maximal non-Hermiticity µ = 1, when taking the asymptotic limit of large arguments (or at ν = ± 1 4 ). However, at small arguments they differ, and so far no matrix representation is known for the weight eq. (2.27). As for all rotational invariant measures the associated orthogonal polynomials are the monomials z k and their complex conjugates, from which we can construct the skew orthogonal polynomials. We obtain the following result in monic normalisation. The skew orthogonality eq. (2.13) can be checked along the same lines as described above for the Laguerre polynomials, with a much simpler recursion of course. We only give the final result for the coefficients r G k , using the normalisation of the orthogonal polynomials computed in the appendix A eq. (A.22), (2.30) Let us also give an example for the simplest correlation function at finite-N , the spectral density. Taking the simplest case without determinants and the weight eq. (2.21) of the two-matrix model we can explicitly write out eq. (2.10): It is manifestly real and vanishes along the real and imaginary axis.

The large-N limit at weak non-Hermiticity
The limit of weak non-Hermiticity was first introduced in [7] for unitary invariant complex non-Hermitian matrices. It is defined by simultaneously taking the Hermitian limit µ → 0 and N → ∞ such that the following combination is kept constant. In this limit the macroscopic spectral density collapses onto the real axis, while the microscopic correlation functions still extend into the complex plane, on a stripe of the order of α 2 . This limit allows us to extrapolate between the correlation functions of real eigenvalues, by taking α → 0, and those of complex eigenvalues at strong non-Hermiticity to be introduced in the next section 4, by taking α → ∞. We define the following rescaling of the complex eigenvalues at the origin where we take over the convention from µ = 0. There eigenvalues are rescaled as V Σx = ξ where 2N = V is the volume, and the inverse variance or chiral condensate proportional to the macroscopic density at the origin is given by Σ = 1/ √ 2 in the matrix model eq. (5.1) at µ = 0. The pre-kernel and the resulting microscopic correlation functions have to be rescaled accordingly: In the following we derive these quantities explicitly for the weight eq. (2.21). We start with the large-N weak non-Hermiticity limit of the skew orthogonal Laguerre polynomials eqs. (2.23) and (2.24). Here, we can use some of the results from [11], which we repeat for completeness. From the asymptotic of the Laguerre polynomials, we obtain for the odd Laguerre polynomials of eq. (2.23) Here we have introduced the scaling variable s ∈ [0, 1]. For the even polynomials eq. (2.24) we have to replace the sum by an integral, k j=0 → k 1 0 dt, introducing a second scaling variable t = j/k ∈ [0, 1]. In terms of the integration variable t we obtain for the even Laguerre polynomial inside the integrand The resulting expression for the sum over Laguerre polynomials eq. (2.24) is then (3.8) Here we have used the fact that the µ-dependent pre-factors turn into an exponential, and similarly for powers of k. In order to obtain the pre-kernel eq. (2.9) it remains to evaluate the k-dependent pre-factors of both skew orthogonal polynomials eq. (2.24) and (2.23), divided by the norm r K k , in the large-N limit, Having collected the asymptotic of all building blocks we can write down the limiting expression for the pre-kernel eq. (2.9) (3.11) The asymptotic of the even and odd skew orthogonal polynomials eqs. (2.24) and (2.23) follow from setting s = 1 in eqs. (3.8) and (3.5) times e ∓α 2 , respectively. The former will be used later in subsection 3.2 as it enjoys a direct interpretation as a characteristic polynomial or massive partition function.
In order to give the spectral density and all higher correlation functions we also need the limiting form of the weight eq. (2.21): Inserting this together with the limiting form of the pre-kernel eq. (3.11) into eqs. (2.8) and (2.7) rescaled according to eq. (3.3) we obtain all quenched microscopic correlation functions, and through eq. (2.20) also all massive microscopic correlations functions. Since the massive correlation functions at real masses involve taking limits we give explicit examples in subsection 3.2 below, after first discussing the simplest example, the spectral density without masses.

Example: the quenched spectral density
Here we give the quenched (N f = 0) microscopic spectral density for an arbitrary topological charge ν as well as a comparison to its real limit α → 0. It is given by inserting eqs. (3.11) and (3.12) into eq. (2.10) We note that in addition to the ν-dependence of the J-Bessel functions the weight is also explicitly topology dependent. This is in contrast to the spectral density of real eigenvalues, see eq. (3.14) below. It is only in the limit of small α 2 or large arguments, |ξ| 2 /(4α 2 ) ≫ 1, that the K-Bessel function becomes ν-independent, K 2ν (x) ∼ e −x π 2x , leading to a Gaussian weight asymptotically (see also eq. (2.22)).
In Fig. 1 we have plotted the microscopic density for ν = 0 and 2 at the same value of weak non-Hermiticity parameter α = 0.4. It is instructive to look at the density of real eigenvalues x for the chSE [30] ρ real (x) = 2x 2 (3.14) Compared to the real density the density in the complex plane keeps the oscillatory structure, with the maxima and minima at the same places. However, due to the repulsion of eigenvalues from the real axis the peak structure is doubled, with a valley along the real axis. On the right of Fig. 2 the additional zero eigenvalues at ν = 2 repel the density further away from the origin, just as they do in the real case. The double peak structure is very different from the complex chUE [11,12,15] and has been detected in lattice simulations [14]. We can also perform an analytical check by taking the limit α → 0 on the microscopic density eq. (3.13). In this limit the pre-factor of the integral becomes  where we denote ξ = x + iy. Taylor expanding the integrand in the density with respect to y = ℑm(ξ), we obtain to leading order after using a Bessel identity. For the limiting spectral density we thus obtain where in the last step we substituted s = u 2 and t = v 2 . Compared to the Hermitian limit of the chiral complex ensemble [11] the density is not proportional to δ(y) but to −yδ(y) ′ . This reflects the repulsion from the real axis. Still it gives unity when integrated over the imaginary part y, after an integration by parts. To illustrate our results further we also give the complex density for larger values of the non-Hermiticity parameter α. As it is seen in Fig. 3 the oscillations are rapidly damped, with the density starting a to develop a plateau as at strong non-Hermiticity (see Fig. 7). For the latter as well as for an analytic check of this relation by letting α → ∞ we refer to section 4. We also remind the reader that the overall height of the density varies with α. This is because of the normalisation to a delta-function in the limit α → 0 eq. (3.15).

Massive correlation functions at weak non-Hermiticity
After discussing in detail the quenched spectral density we now turn to the presence of massive flavours. We present the spectral density for one pair of degenerate massive flavours 4N f = 4 explicitly, as well as some characteristic polynomials or massive partition functions as examples. Higher order correlation functions and more flavours can then easily be obtained from the given building blocks eq. (3.11), following the same procedure.
At finite-N the simplest spectral density with non-zero flavour content has 4N f = 4 masses. It is given by While the first term gives exactly the quenched spectral density the second term is a correction term, just as for zero chemical potential [25]. When taking the large-N limit we have to rescale the masses in the same way as the complex eigenvalues in eq. (3.2), As we have mentioned at the end of subsection. 2.2, eq. (3.19) includes a limit when setting the masses to be real, m f = m * f . This is because of the level repulsions both from the real and imaginary axis: the denominator as well as the numerator of the second term in eq. (3.19) vanishes. The same feature also holds for more general correlation functions eq. (2.20). We illustrate this by Taylor expanding the limiting pre-kernel κ weak (iη, −iη * ) in the denominator in the imaginary part of the mass ǫ → 0, The imaginary argument rotates J-to I-Bessel functions, where we could have obtained the same result from eq. (3.16), setting x = iη and y = iǫ there. Making the same expansion for the other terms, κ weak (z, ±iη ( * ) ) and its conjugates in the numerator, we arrive at the following expression for the microscopic limit of the density eq. (3.19): with the first term from eq. (3.13) and the correction term to the quenched density being given by Here we have again called η x = η which is now real. Before discussing this main result of the subsection we also briefly give 2 examples for massive partition functions. The simplest case with 2N f = 2 flavours follows from the even skew-orthogonal polynomial at imaginary mass z = im, comparing eq. (2.14) and (5.8): We can set the mass to be real here without taking limits, as q 2N (im) itself does not display any level-repulsion from the axis. We obtain the remarkably simple expression in the weak limit, where we have rotated the result eq. (3.8) from the previous subsection to an imaginary argument, skipping all constants (including N ). As a second example we give the partition function for 4N f = 4 using our formula eq. (2.17), Upon using the same expansion as for the density eq. (3.21) we immediately obtain (3.27) where we have again suppressed all constants. Due to the presence of conjugate pairs of complex eigenvalues the dependence on α cannot be factored out, similar to the same situation in the unitary ensembles [27,12,16]. Such a simplification only occurs in the absence of conjugate flavours, in the QCD matrix model with a sign problem [31,15,16], where the partition function can be normalised to be µ-independent. As a check both expressions eqs. (3.25) and (3.27) reduce to the know results [32,25] in the limit α → 0.   Let us now discuss the massive density eq. (3.22). In Fig. 4 we plot the microscopic, massive spectral density for two different values of masses η = 8.74 and η = 4.26, both at weak non-Hermiticity parameter α = 0.4 and ν = 0, and compare it the corresponding expression for real eigenvalues [25].
We have deliberately chosen the same values for the mass parameter as there, where numerical data for lattice simulations with dynamical fermions at µ = 0 [33] were compared. Similar to having ν = 2 exact zero eigenvalues in Fig. 1 the eigenvalues are pushed further away from the origin through the masses. Below we compare the corresponding expression for real eigenvalues at the same values of the mass parameters [25]. As for the quenched case the correlation functions keep their oscillatory structure, with the maxima located at the same positions as for real eigenvalues, at the given value of α. Furthermore we observe that for larger masses the density moves toward the quenched density in  On the other hand smaller masses moving towards the origin push the eigenvalues further out. In the zero-mass limit the density should approach the quenched density at topological charge ν = 2 approximately, as can be seen in Fig. 5 left. Although for non-zero chemical potential there is no more flavour-topology duality [34] because of the explicit ν-dependence in the K-Bessel functions in the weight, this duality holds approximately, when the K-Bessel function reaches its asymptotic, ν-independent form. For those values of |ξ 2 |/(2α 2 ) in Fig. 6 left, where the density is non-vanishing, this is indeed the case. As a final example we display the massive density at a larger value of α = 2.5 in Fig. 6 right. As seen already for the quenched density the oscillations are washed out (compared to Fig. 4 right). A further feature can be seen here. For the chosen value of α and η the ξ-value where the mass factor |ξ 2 + η 2 | 2 vanishes, ±ℑm(ξ) = η, now lies on the imaginary axis inside the shown picture. It can be seen very clearly that the eigenvalues are pushed away from the very location of the mass, in addition to the level repulsion from the y-axis.
In the previous subsection we performed an analytical check by taking the Hermitian limit α → 0 of the spectral density and comparing it to the known expression for the density of real eigenvalues. We have performed the same check here where we refer to [25] for the corresponding expression for the massive spectral density of real eigenvalues. Without going through the details we just comment of the structure to be obtained.
Since our example eq. (3.18) involves the two-point function it contains all four elements of the matrix kernel eq. (2.8), the pre-kernel κ N at different (complex conjugated arguments). In contrast to that for real eigenvalues the two-point or higher correlation functions are given by a quaternion determinant of a two by two matrix kernel with different matrix elements: a pre-kernel which is basically given by eq. (3.14), its derivative and its integral (see [30]). How can we obtain these three different matrix elements from one pre-kernel? Taking the limit α → 0 invokes a further derivative with respect to the imaginary part y = ℑm(ξ), in addition to the expansion in eq. (3.23) due to real masses. As a result we will get the real limit of the pre-kernel eq. (3.11), its first and second derivative, correctly reproducing all terms in [25] for the real massive density.

The strong non-Hermiticity limit
The regime of strong non-Hermiticity in defined by taking the limit N → ∞ while keeping µ ∈ (0, 1] fixed. The rescaling of the eigenvalues is modified compared to eq. (3.2), keeping finite 6 . Hence the definition of the large-N correlation functions is also altered, reading The difference in scaling can be easily understood as follows. At weak non-Hermiticity the eigenvalues remain close to the real axis, with their imaginary part being small compared to the real part. In order to obtain a macroscopic density at the origin of the order ρ macro (0) ∼ N their average distance has to be O(1/N ). In contrast to that at strong non-Hermiticity the eigenvalues spread out in the complex plane, having a comparable real and imaginary part. To form again a constant macroscopic density ρ macro (0) ∼ N the average distance between nearest neighbours has to be O(1/N 2 ). Here we measure distance with respect to the two-dimensional metric. A clear separation of scales between the two different limits can be defined as follows. At weak non-Hermiticity the microscopic density describes the universal fall off to zero of the eigenvalues away from the imaginary axis. In direction of the real axis along the maxima we obtain as usual lim ℜez→∞ ρ weak (z) = ρ macro (0) 7 . In contrast to that at strong non-Hermiticity the microscopic density attains the constant value of the macroscopic one in all directions, lim |z|→∞ ρ strong (z) = ρ macro (0), as can be seen below in Figs. 7 and 8. This is of course except along the real and imaginary axis, due to the repulsion from the Jacobian in eq. (2.4). The fall-off of the macroscopic density to zero in any direction (which is not seen in the microscopic regime here) can then be expected to be given by another microscopic function at the edge, in analogy to the Airy-kernel for real ensembles at the edge of the spectrum. When numerically solving the matrix model at finite but large-N care has to be taken to ensure a clear separation between micro-and macroscopic density when comparing to the analytic formulas.
The computation of correlation functions at strong non-Hermiticity is more involved, as is already the case in the non-chiral ensemble [19]. In fact we will follow the same strategy as there, by first deriving a differential equation for the large-N pre-kernel at maximal non-Hermiticity µ = 1. A solution of this equation can be obtained either directly or by taking the limit α → ∞ of the prekernel at weak non-Hermiticity, keeping ξ weak /α fixed. Since the two methods match it is suggestive to recover the full range of strong non-Hermiticity µ ∈ (0, 1] by inserting a µ-dependent parameter through ξ → ξΣ(µ). However, we can also argue that such a rescaling can always be absorbed by redefining the scaling limit eq. (4.1). Details of the derivation of the differential equation and its solution are given in the appendix B.
We begin by taking the limit of maximal non-Hermiticity µ → 1 on the skew orthogonal polynomials eqs. (2.23) and (2.24). In this limit only the highest powers survive and they become monic or a sum over monomials, respectively. We thus obtain the following expression for the pre-kernel eq. (2.9) in the microscopic large-N limit eq. (4.2): It can be shown that the function κ strong (ξ, ζ * ) obeys the following set of inhomogeneous second order differential equations where we refer to the appendix B for details. Due to the antisymmetry of the pre-kernel the right hand side changes sign. The solution of this equation can be constructed as explained in the appendix B. The result is given by In order to obtain all correlation functions from eq. (4.5) we next give the weight function eq. (2.21) in the limit of strong non-Hermiticity eq. (4.1): Here we display the weight for any value of the strong non-Hermiticity µ. At µ = 1 the expression simplifies, the exponential factor cancels with the kernel and the argument inside the K-Bessel function reduces to unity times |ξ| 2 . Inserting this into eq. (2.10) we obtain the following result for the microscopic spectral density at maximal non-Hermiticity As an analytic check we have performed the limit α → ∞ of the microscopic spectral density at weak non-Hermiticity eq. (3.13), where we refer to the end of appendix B for details. We obtain exactly the same functional form of the spectral density, eq. (B.20) (up to a different constant normalisation), in terms of the variable ξ S = √ N z/µ √ 2. Together with eq. (4.6) this leads us to conjecture that the µ-dependence at strong non-Hermiticity with µ < 1 can be restored as follows, replacing everywhere in eq. (4.7). We have shown it to be true for µ = 1 and for small µ ≪ 1. However, if this is true we could then redefine our strong non-Hermiticity limit eq. (4.1) to completely reabsorb this µ-dependence, defining As an illustration of our results we show the microscopic spectral density eq. (4.7) and its dependence on ν at maximal non-Hermiticity µ = 1 in Fig. 7. One can clearly see that for ν = 2 the additional zero eigenvalues push the density away from the origin. For the numerically evaluation of the integral in eq. (4.7) necessary to plot it an equivalent representation of the integral is very convenient. Using the identity derived in the appendix eq. (B.23) we have where the exponential pre-factor has cancelled out. This form of the density allows to dicuss its rotation and reflection symmetries most clearly. The fact that it is an even function of (ξ 2 − ξ * 2 ) and otherwise only depends on the modulus |ξ| 2 leads to the following properties. First, one can obviously rotate the density by π/2, π, or 3π/2, multiplying ξ by i , (−1), or −i without changing it. Second, the density is invariant under reflections along the x-, y-, (x = y)-and (x = −y)-axis where ξ = x + iy.

Massive correlation functions at strong non-Hermiticity
In this subsection we give the simplest nontrivial result including one pair of degenerate massive flavours, 4N f = 4, where we will follow closely the corresponding subsection 3.2 at weak non-Hermiticity. Higher order correlations as well as more flavours follow along the same lines. We start by defining how to rescale the masses √ N m f ≡ η f , f = 1, . . . , N f , (4.12) in analogy to the eigenvalues in eq. (4.1). In order to derive the massive spectral density we could simply take the limit α → ∞ of the expressions in eqs. (3.13) and (3.23) at weak non-Hermiticity, together with the results from appendix B. In order to give more compact expressions we derive results using the alternative form of the strong density eq. (4.10) and the corresponding kernel. The massive spectral density follows again from eq. (3.19) at finite-N . The kernel at strong non-Hermiticity we use is obtained from eq. (4.5) by inserting eq. (B.23): Expanding again in the imaginary part of the mass ǫ → 0, setting η = η x + iǫ we easily obtain for the denominator of eq. (3.19) (4.14) The same expansion can be done for κ strong (ξ, ±η ( * ) ) and its conjugates, leading to the following expression for the massive spectral density at strong non-Hermiticity: with the quenched density given by eq. (4.10). The correction term reads ∆ρ (4) (4.16) We have again called η = η x which is now real. Before discussing the massive spectral density let us also give some examples for partition functions. Taking the strong non-Hermiticity limit of eq. (3.26) at 4N f = 4 we simply have to insert the expansion eq. (4.14) to obtain up to suitable normalisation constants. Since we did not compute the asymptotic even skew-orthogonal polynomials at generic strong non-Hermiticity 0 < µ ≤ 1 we cannot directly use eq. (3.24) for 2N f = 2.
We can circumvent this problem by taking the infinite α limit on eq. (3.25), leading to the following Let us come back to the spectral density eq. (4.16). Below we show a few examples for the influence of the mass parameter at maximally strong non-Hermiticity. Compared to the quenched density in  is an additional level repulsion from the mass, when ±ℑm(ξ) = η. Non-zero topology leads to an additional level repulsion from the origin as can be seen in Fig. 9 right, and both the effects of mass and ν can be clearly separated (see also Fig. 7 right).
For smaller values of the mass parameter the density approaches approximatively the quenched density at ν = 2, as the mass acts as additional zero eigenvalues (comparing Fig. 10 left and Fig. 7  right). On the other hand going to larger masses we get back the quenched density as can be seen in Fig. 10 right. Of course there is still level repulsion at the large mass value, in this case η = 8, as displayed in Fig. 11. However, this does no longer affect the density at the origin. Furthermore, there are no strong oscillations observed here at the location of the mass, as seen in unquenched QCD [16]. Here, the situation is very similar to phase quenched QCD also reported in [16] on the matrix model in the β = 2 symmetry class. These findings underline once more the importance of the presence or absence of the sign problem on the level of the microscopic complex eigenvalue correlations of the Dirac operator. We note in passing that because of the analytic relationship between the correlations at  Fig. 10 (right) showing the first quadrant for a larger argument ℑm(ξ). The level repulsion at the mass value ξ = 8i is clearly visible. strong non-Hermiticity and weak non-Hermiticity at large α similar figures could have obtained in this section by plotting the corresponding pictures at large enough α. We have checked this numerically in several cases for both the quenched and massive density, with α = 5 (see also Fig. 6). However, such a density at weak non-Hermiticity will always vanish for large enough ±ℑm(ξ), at any given α. In contrast to that the density at strong non-Hermiticity remains constant.

Relation to the Dirac operator spectrum with chemical potential
The purpose of this section is twofold. First, we give a more detailed discussion of the matrix representation of our complex eigenvalue model eq. (2.4) in terms of two non-Hermitian quaternion real matrices, including mass terms. This justifies to speak of a matrix model ensemble. It also gives an explicit realization of the joint probability distribution in eq. (2.4) by computing the Jacobian that leads to complex eigenvalues. In particular it yields the special weight function eq. (2.21) containing a K-Bessel function. The second purpose is to relate the model to the Dirac operator spectrum of QCD-like gauge theories in the presence of a chemical potential. To do this matching we need to compare the global symmetries of such Dirac operators to a random matrix representation.
The Dirac operator as it appears in four dimensional gauge theories can be represented by an anti- , with non vanishing entries only on the off diagonal due to chiral symmetry. For different gauge groups and representations the off-diagonal blocks W have either a orthogonal, unitary, or symplectic symmetry, as was pointed out in [9]. Adding a chemical potential µ for the quarks amounts to shifting the Dirac operator by / D gauge → / D gauge + µγ 0 where γ 0 is a Dirac-γ matrix. In a given basis it is represented as a block matrix as well, with unit elements on the off diagonals only [2]. Being Hermitian this addition destroys the anti-Hermiticity property, where the sum / D gauge + µγ 0 now has complex eigenvalues. The generalisation of the classification [9] to µ = 0 we use here was made in [10].
If we replace the Dirac operator by a random matrix we have to satisfy two criteria. First, / D gauge may have exact zero eigenvalues corresponding to different sectors of topological charge through an index theorem. For a fixed number ν of zero modes, which we shall always consider here, we thus have to replace the blocks W by rectangular matrices, of size (N + ν) × N . Second and more importantly, these rectangular matrices have to be in the same symmetry class as the corresponding gauge theory, consisting of real, complex or quaternion real elements [9,10]. The latter case which we consider here, also labelled by the Dyson index β = 4, describes SU (N c ) gauge theories in the adjoint representation. When representing the gauge theory on a lattice using staggered fermions the corresponding symmetry is an SU (N c = 2) gauge theory in the fundamental representation.
The addition of a chemical potential in the matrix model can be done in different ways. It was first suggested [2,10] to add it as a constant times γ 0 , the off diagonal unity matrix, assuming this term to be diagonal in the random matrix space. While this model successfully explained the failure of the quenched approximation [2] an analytical treatment of all complex eigenvalue correlations from orthogonal polynomials similar to the µ = 0 case was not achieved to date, apart from the quenched spectral density obtained using the replica method [12]. Therefore a different matrix model representation of the chemical potential term was suggested very recently in [15] for the QCD symmetry class β = 2. It assumes that the chemical potential term is non-diagonal in matrix space and is represented by a second, uncorrelated matrix which has the same symmetries as the one modelling the µ = 0 part. Making the model more complicated at first sight it has the remarkable feature of leading to a complex eigenvalue representation, which then allows for a treatment using (bi-)orthogonal polynomials [15,16]. A different model given only in terms of complex eigenvalues was first suggested in [11] and tested in quenched QCD lattice simulations [13]. Although it is lacking the matrix symmetries of / D gauge it asymptotically agrees with the results of [12,15] (see the weight eq. (2.22) vs. eq. (2.21)). In the following we will adopt the strategy of [15] and apply it to the complex β = 4 symmetry class. Our corresponding matrix model is thus given in terms of two rectangular (N + ν) × N matrices, Φ and Ψ, with quaternion real elements without further symmetry properties: where 1 is the quaternion unity element. Compared to eq. (2.1) we have now included the quark masses. As a first step we define the linear combinations with a trivial constant Jacobian. In terms of these new matrices the partition function reads In this form we have only convergent integrals for µ ∈ [0, 1] at finite-N , as can be seen from eq. (5.8) below. We will not discuss the possibility of phase transitions here, where for lattice simulations of this situation we refer to [35]. The Dirac matrix inside the determinant has ν zero eigenvalues and N complex eigenvalues which come in pairs of opposite sign as well as complex conjugate pairs. The former is due to chiral symmetry and the latter is due to the symplectic structure as the eigenvalues of a quaternion real matrix come in complex conjugate pairs. It is this symplectic structure that is responsible for the fact that the determinant will remain positive definite, without having a sign problem as unquenched QCD. It makes this symmetry class an ideal testing ground for lattice simulations with dynamical fermions in the presence of a chemical potential.
In the next step we parametrise the matrices C and D as follows, which is equivalent to making an independent Schur decomposition of CD and DC. Here U and V are symplectic matrices. The matrices R and S are upper triangular with real quaternion elements, and we have split off their diagonal parts X and Y which are also quaternion real, without containing the second and third quaternion unit (see e.g. [17] for explicit representations). They contain the complex eigenvalues of C and D which come in N pairs (x k , x * k ) for the matrix X and (y k , y * k ) for the matrix Y respectively. The parametrisation is not unique as under an additional diagonal unitary transformation the matrices remain triangular. It only changes the upper triangular matrices R and S which we will integrate out afterwards. This symmetry can be used to restrict the matrices to be V ∈ Sp(N )/U (1) N and U ∈ Sp(N + ν)/(Sp(ν) × U (1) N ). The parametrisation eq. (5.4) as well as the counting of degrees of freedom is discussed most explicitly in appendix C. The relation to the original eigenvalues of the Dirac matrix is after a suitable ordering of the independent variables. We can now integrate out the symplectic matrices U and V , as well as the Gaussian integrals over R and S to obtain Here we have included the mass dependent factor ∼ |m f | 2ν in the normalisation, to the power of mass degeneracy 2 in 2N f . Although we have now achieved an eigenvalues representation in terms of x k and y k we still have twice as many degrees of freedom that what we aimed at, a single set of N complex eigenvalues z k . If we substitute y k by z k : y k = −z 2 k /x k we can integrate out the variables x k in analogy to [15], after decomposing them into polar coordinates, This gives us a complex eigenvalue model as in eq. (2.4) with the specific weight function eq. (2.21) plus mass terms, as well as a matrix representation of the q 2N (z) eq. (2.14) for N f = 1. We note that compared to the QCD symmetry class [15] Let us stress that it is not this imposed degeneracy that is responsible for the absence of a phase in the weight. The non-degenerate partition function eq. (5.8) has already a positive definite weight. Eq. (5.9) provides a determinant representation for the mass terms in eq. (2.16) for which we have been able to compute all complex correlation functions. The computation of partition functions without 2-fold degenerate masses eq. (5.8) remains an open problem for 2N f > 2, where the special case 2N f = 2 is given by eq. (3.24). In the limit µ → 0 the partition function eq. (5.9) reduces to using eq. (3.15). This is the partition function of the chGSE with real eigenvalues computed in [25,32] and successfully compared to lattice simulations with dynamical fermions at µ = 0 in [33,25]. It matches the corresponding chiral Lagrangian in the ǫ-regime calculated for equal masses in [36].

Conclusions
In this paper we have solved the complex extension of the chiral or Laguerre Symplectic Ensemble to non-Hermitian matrices. The solution was shown to be expressible in terms skew orthogonal polynomials for general weight functions in the complex plane. We gave explicit expression for finite-N for all correlations functions of eigenvalues in the absence and presence of twofold degenerate mass terms, as well as for characteristic polynomials. Two examples of a Gaussian Ginibre type weight and a weight with a K-Bessel function were given. We proved that the (skew-)orthogonal polynomials corresponding to the latter are Laguerre polynomials in the complex plane. We investigated the large-N limit of complex skew orthogonal Laguerre polynomials, first treating weak non-Hermiticity, and derived all correlation functions in terms of the limiting kernel and weight function. These finding were illustrated by examples of the microscopic spectral density without and with mass terms as well as by partition functions, where we discussed the influence of both exact zero eigenvalues and the masses. The same analysis was then repeated for strong non-Hermiticity where we had to solve two inhomogeneous differential equations to determine the limiting kernel. The virtue of the weak non-Hermiticity limit is that it allows to extrapolate in terms of the weak non-Hermiticity parameter α between real eigenvalue correlations at α → 0, and strong non-Hermiticity in the limit α → ∞. This provided several analytical checks.
These mathematical results may have applications in various physical systems with complex operators, and we have focused on the application to the Dirac operator spectrum in QCD-like gauge theories. We have identified the correct global symmetries of our ensemble by finding a two-matrix representation that matches with the symmetries of the Dirac operator with chemical potential in SU (N c ) gauge theories in the adjoint representation, or SU (N c = 2) gauge theories in the fundamental representation using staggered fermions in a lattice discretisation. Our two-matrix model representation led to the particular form of a K-Bessel weight function of complex eigenvalues, which is identical to that of the corresponding model in the QCD symmetry class. However, the solution we presented in terms of complex eigenvalues holds for more general classes of weight functions as well.
It would be very interesting to verify our predictions in lattice simulations with dynamical fermions in the corresponding theories with chemical potential. In contrast to QCD these do not suffer from a sign problem and thus can be solved using standard techniques. This project is work in progress.
holds for a > b when integrating over the full complex plane , where we have introduced Note that throughout this appendix we use the shifted topological index 2ν → ν as it appears in the complex chUE [15], compared to the weight eq. (2.21). In the following we can restrict ourselves to the case k ≥ l in eq. (A.1), as then k ≤ l follows from complex conjugation. To show eq. (A.1) it is thus sufficient to prove The remaining case k = l in eq. (A.1) is trivially non-vanishing, as we integrate over a positive quantity. The corresponding norm of the k-th polynomial is computed below as well. To set up our orthogonality proofà la Gram-Schmidt we first compute the following general integral, After introducing polar coordinates we can integrate out the angle to obtain the modified I-Bessel function. Changing variables r 2 = s the remaining integral can be found in terms of the hypergeometric function [29] where convergence is guaranteed due to a > b. We note that all the moments c kl are real. The case k ≤ l can thus simply be obtained from interchanging k and l: We first consider the orthogonality of the zero-th polynomial to L ν k . For l = 0 the hypergeometric function in eq. (A.4) reduces to a simple rational function, F (n, β; β; z) = (1 − z) −n , and we obtain Using the fact that the Laguerre polynomials have an explicit representation, we can simply write down the orthogonality relation between the k-th and 0-th polynomial as Here we have used the definitions (A.2) as well as N (1−µ 2 ) = a 2 −b 2 2b and a 2 − b 2 = N 2 /µ 2 . The insertion of the coefficient c m0 eq. (A.6) leads to the binomial sum (1 − 1) k that vanishes for all k > 0. For k = 0 we obtain the norm of the 0-th polynomial reading In the general case eq. (A.3) the strategy of our the proof will be similar. We shall try to reduce the expressions to vanishing binomial sums. However, there will be remaining terms, which we have to show to cancel as well. As a preliminary step we can use a transformation formula [29] to show that the hypergeometric functions in the coefficients c kl eq. (A.4) terminate, Because of ν ≥ 0, and k ≥ l ≥ 0 being integers the hypergeometric function on the right hand side terminates, containing l + 1 terms. Writing down the definitions we obtain for eq. (A.3) where we have split the sum into two contributions for m < l and m ≥ l in c ml . In the sequel we manipulate the Σ 1 and Σ 2 separately by rearranging the hypergeometric functions, showing that they indeed cancel for k > l. We start with the simpler second sum, using eq. (A.10): In the next step we reduce the individuals sums sorted by powers of b 2 /a 2 = (1 − µ 2 ) 2 /(1 + µ 2 ) 2 to binomial sums by shifting the summation indices: In the last step we have used k−l+1+n , replacing incomplete binomial sums. In particular we need k > l for the first term. Also we have multiplied in the factor (1+µ 2 ) 2l (1−µ 2 ) l . The generic n-th term in the sum is explicitly displayed. In this form we can compare to the first sum Σ 1 , again after doing some manipulations: In contrast to eq. (A.12) the length of each hypergeometric sum now increases from term to term. Let us therefore write them out most explicitly: In order to see that Σ 1 and Σ 2 cancel we now have to compare individual powers a 2 b 2 = (1+µ 2 ) 2 (1−µ 2 ) 2 . For the for highest power, a 2(l−1) b 2(l−1) = (1+µ 2 ) 2(l−1) (1−µ 2 ) 2(l−1) , it is easy to see that this is indeed the case, with only the first term in the last line of Σ 1 in eq. (A.15) contributing. In contrast to that for the lowest power a 0 b 0 the last term in each hypergeometric sum in eq. (A.15) does contribute, leading to 16) This clearly cancels the last term in eq. (A. 13). To see that all other terms cancel too we now write down all contributions to the generic power (1+µ 2 ) 2l−2n−2 (1−µ 2 ) 2l−2n−2 to be compared with that of Σ 2 in eq. (A.13). In Σ 1 we thus have to pick all terms a 2q−2p b 2q−2p with 2q − 2p = 2l − 2n − 2. Because of n ≤ l − 1 and p ≤ q only those hypergeometric sums with 2q ≥ 2l − 2n − 2 start contributing. Eliminating p = q − l + n + 1 we obtain the following contribution where we have shifted summation variables in the last step. This term clearly cancels the corresponding term in Σ 2 in the last but one line of eq. (A.13). This ends our proof of eq. (A.3) and therefore of the orthogonality of the Laguerre polynomials in the complex plane. What remains to be computed are the norms of the k-th Laguerre polynomial, which can be obtained as follows. The only step where we used k > l was when producing a zero in the first term of Σ 2 in eq. (A.13), deducing (1 − 1) k−l = 0. For k = l this is no longer true and this term does contribute, leading to the only non-vanishing contribution from Σ 1 + Σ 2 . Because of the orthogonality only L ν k (z 2 )z * 2k will contribute to ||L ν k || 2 . To obtain the norm for the polynomials in monic normalisation, we still have to normalise eq. (A.11) appropriately, leading to which is our final result on orthogonal Laguerre polynomials in the complex plane. In order to clarify the issue of orthogonality of complex Laguerre polynomials with respect to different weight functions we have explicitly checked that L ν 0 and L ν 1 are not orthogonal with respect to the exponential Gaussian weight in eq. (2.22) times |z| 2ν+1 for a general value of ν, as was claimed in [11]. The integral is given by a linear combination of Legendre functions of index ν + 1 2 and ν + 5 2 which does not vanish in general. Only in the special case of ν = ± 1 2 where the Laguerre polynomials are related to Hermite the orthogonality holds. This does not come as a surprise now in view of the simplification of the K-Bessel weight eq. (2.21) to eq. (2.22) in that case.
For comparison we shall also give the norms for the Gaussian weight eq. (2.27). It trivially holds that for all rotational invariant weight functions w(z, z * ) = w(|z|) the monic powers z k form a set of orthogonal polynomials (A.20) For different weight functions only the corresponding norms differ, which then enter the kernel, or the construction of the skew orthogonal polynomials as in section 2.3. For the Gaussian weight eq. (2.27), with shifted 2ν → ν one easily obtains the following result, after changing to polar coordinates B Differential equation for the kernel at strong non-Hermiticity B.1 Construction of the differential equation In order to derive the differential equation we first rewrite the coefficients inside the pre-kernel eq. (4.3) in a more symmetric way. Using the following identity for the Γ-function called doubling formula [29], It is convenient to change variables, and to define the following functional which is antisymmetric in the variables u, v, This leads to where we have shifted the sum over j in the first term and denote (−1)!! = 1. This compares to multiplication by v reading as follows where we have shifted the sum over k in the second term. Subtracting this from eq. (B.6) we obtain In the first step we have again used the relation (B.1) and the last step follows upon the series representation of the I-Bessel function. Going back to the kernel with eq. (B.5) in the changed variables eq. (B.3) the differential equation (4.4) follows, The corresponding equation in the variable u = ξ 2 2 is obtained in the same way, or by simply using the antisymmetry of the function σ(u, v) B.2 Solution of the differential equation Next we seek for a solution of these two inhomogeneous differential equations. The homogeneous equation can be easily solved, as it holds in both variables u and v. The second independent solution in terms of K-Bessel functions is excluded because of the regularity of the pre-kernel at the origin. We therefore conclude that the following product solves both homogeneous equations. However, this solution is symmetric in both arguments, in contrast to the antisymmetric pre-kernel. Each factor of the homogeneous solution has the following integral representation that will be useful below. In terms of the original variables it is given by eq. (B.3) and similarly for ζ * . Instead of making an Ansatz that solves the full set of inhomogeneous differential equations (B.9) and (B.10) we construct a solution by rewriting the right hand side as an integral and manipulating it until we obtain the correct differential operator acting on it. Omitting the constant pre-factor we start from the following integral identity, Here we first integrated by parts e −q 2 J 2ν (2qζ * ), with a vanishing surface term. Next we replace q∂ q J 2ν (2qξ) = ξ∂ ξ J 2ν (2qξ), and then use that J 2ν (2qξ) satisfies a Bessel differential equation: We can now take the differential operator out of the integral and interchange it with the ξ-dependent pre-factor, giving us precisely the desired differential operator. This provides a solution of the second differential equation (B.10), after multiplying with the appropriate constant. A solution of the corresponding equation (B.9) for ζ * follows along the same lines (B.16) However, none of two solutions is antisymmetric in the variables ξ and ζ * to be proportional to the pre-kernel, nor do they obviously solve both differential equations with the different signs simultaneously. We will therefore add and subtract a multiple of the homogeneous equation eq. (B.12) in its integral representation eq. (B.13) to eq. (B.14) and from (B.16), respectively, in order to obtain an antisymmetric solution of both differential equations. The requirement of antisymmetry of the pre-kernel fixes the solution uniquely. It can be seen that the solution is given as follows The first line is antisymmetric in its arguments ξ and ζ * . In the subsequent lines we have made an integration by parts of either e −q 2 J 2ν (2qζ * ) or −e −q 2 J 2ν (2qξ) in order to show, that this expression solves both differential equations, (B.14) and (B.16) plus or minus a solution of the homogeneous equation. The full solution for the pre-kernel including all factors thus reads As a check and alternative way to derive this result we take the limit α → ∞ of the pre-kernel at weak non-Hermiticity eq. (3.11). Substituting first t by p 2 = 2stα 2 and then s by q 2 = 2sα 2 we obtain from eq. (3.11) lim α→∞ κ weak (ξ, ζ * ) = lim If we then replace the first integral α √ 2 0 dq → ∞ 0 dq, keeping ζ * α √ 2 and ξ α √ 2 fixed, we reproduce the strong kernel eq. (B.18) up to constants and the exponential pre-factor. Therefore we cannot directly identify the kernels, but have to take into account the weight function as well when mapping correlation functions from the weak to the strong limit. From eq. (3.12) we thus obtain exactly reproducing the strong density eq. (4.7) obtained from eq. (B.18). Here we have introduced the scaling variable with ξ being the scaling variable in the weak limit eq. (3.2). ξ S scales as the scaling variable in the strong limit eq. (4.1). The explicit factor 2N in front of the density is added since the weak and strong density are rescaled with different powers of N , see eqs. (3.3) and (4.2). The explicit µ-dependent factor in front can be absorbed by defining a µ-dependent strong scaling limit by ξ S , and rescaling the density accordingly. It can be easily seen that all higher eigenvalue correlation functions at weak and strong non-Hermiticity can be matched in the same way when taking α → ∞. We end this appendix by noting that another integral representation can be obtained from the α → ∞ limit of the weak pre-kernel eq. (3.11). First we substitute s = u 2 and t = v 2 . Next we interchange the integrals with respect to u and v and substitute u by y 2 = u 2 α 2 (1 + v 2 ) to arrive at In the last step we have substituted r = 2v/(1 + v 2 ). It is this last form which is very convenient for a numerical evaluation in order to plot the spectral density at strong non-Hermiticity.

C Jacobian of the two-matrix model
In this appendix we compute the Jacobian for the parametrisation eq. (5.4) resulting from Schur decompositions of CD and DC. To this aim we proceed in two steps. First, we provide a labelling of the independent degrees of freedom to obtain a block triangular form of the Jacobi matrix in terms of quaternion matrix elements. Here we follow some analogy to [17] appendix A.35 where a similar decomposition was given for a single quadratic complex non-Hermitian matrix. Second, we calculate the Jacobi determinant, evaluating the diagonal blocks. In the same way as above we define for the diagonal elements Q l,l ≡ q l,l = y l 0 0 y * l , l = 1, . . . , N , (C. 6) containing the eigenvalues of D. The symplectic matrices U and V of size (N + ν) 2 and N 2 a priori carry 2(N + ν) 2 + (N + ν) and 2N 2 + N d.o.f., respectively. However, the parametrisation eq. (C.1) is not unique. Replacing the quaternion elements in all matrices by 2 × 2 blocks one can see that a diagonal unitary matrix U ′ of size 2N leaves the matrices T and Q triangular, keeping the eigenvalues unchanged: whereŨ ′ is U ′ extended to size 2(N + ν) by the unity matrix. Therefore 2N real parameters are redundant in eq. (C.1). Furthermore, the whole lower right ν × ν sub-block in U can be projected out, subtracting a symplectic matrix with 2ν 2 + ν d In oder to label independent d.o.f. for the computation of the Jacobian it is most convenient to define the following differentials: where dH ≡ +idU † U = −iU † dU dJ ≡ −iV dV † = +idV V † (C.10) Both dH and dJ are anti-self dual [17]. Their (upper N) diagonal elements are given by dH l,l ≡ dh l,l = 0 ds l ds * l 0 , l = 1, . . . , N dJ l,l ≡ dj l,l = 0 dp l dp * l 0 , l = 1, . . . , N . (C.11) where we have used the 2N parameters of the symmetry eq. (C.7) to set the real parameters on the diagonal to zero. The remaining parameters on the upper and lower triangle are related due to the anti-self duality, and we only keep the latter in the following set of independent variables: dH : {dh l,l , dH n,l |n > l; n = 1, . . . N + ν, l = 1, . . . N } dJ : {dj l,l , dJ n,l |n > l; n, l = 1, . . . N } , (C. 12) where in dH we have also projected out the ν 2 -block. They label 4N (N − 1)/2 + 4N ν + +2N = 2N 2 + 4N ν and 4N (N − 1)/2 + 2N = 2N 2 real variables respectively, matching those of U, V in eq. (C.8).
We can now arrange the independent variables {dC, dD} and {dT, dQ, dH, dJ} in such a way that the Jacobian is block triangular. It is most efficient to use quaternion real matrix elements for this ordering. In a second step we then have to compute the determinant of the Jacobi matrix by taking the product of the determinants of the diagonal blocks of quaternion elements. The variables {dC, dD} are ordered row by row, starting with the lowest row of dC, until we reach its quadratic part. We then alternate dC and dD element by element, completing the non-quadratic part of the row with dD only. In detail the ordering looks as follows: The independent variables {dT, dQ, dH, dJ} can also be put into two full matrices without zeros, dA and dB of size (N + ν) × N and N × (N + ν) respectively, and then ordered in the same fashion as dC and dD. The matrix dA consists in its strictly lower triangular part of dH. The diagonal part is composed of dh l,l and dt l,l giving full quaternion real matrix elements: dT l,l ≡ dh l,l + dt l,l = dx l ds l ds * l dx * l , (C.14) and its strictly upper triangular elements are completed by dT . Similarly we define a matrix dB by putting dJ in its strictly lower triangular part, the sum of dj l,l and dq l,l on its diagonal: dQ l,l ≡ dj l,l + dq l,l = dy l dp l dp * l dy * l , (C. 15) and dQ on its strictly upper triangular part. These matrices dA and dB are then ordered exactly as in eq. (C.13), or most explicitly as: For the matrix elements related through anti-self duality we have used the relation dH n,k = −dH k,n , where¯denotes the conjugate quaternion [17]. Of course it holds that dH k,n /dH k,n is non-vanishing. However, we will not need to compute these matrix elements as they appear on the off-diagonal of the Jacobi matrix only. Let us distinguish 5 cases to show the block-triangularity of the Jacobi matrix. The differentiation of quaternion elements by quaternion elements used here is to be understood in a symbolical sense, in order to see the block structure of the Jacobi matrix. In a second step we will differentiate independent matrix element by independent matrix element.
(C. 19) The variable dC k,l only depends on dH k,l : dC k,l /dH k,l = it l,l , (C. 20) and not on variables to the right of dH k,l in the ordering eq. (C.16), that is dH k,n>l , dH n<k,m , or any of the dJ, dQ, dT .
ii) k > l; k = N . . . , 2, l = 1, . . . , N : We have both differentials dC k,l an dD k,l simplifying to This gives the following 2 × 2 entry on the diagonal: dH k,l dJ k,l dC k,l it l,l −it k,k dD k,l −iq k,k iq l,l (C. 22) None of the differentials depends on variables further to the right of dH k,l and dJ k,l in eq. (C.16).
iii) k = l; k, l = 1, . . . , N : dC k,l = dt k,k + i k−1 n=1 dH k,n T n,k + idh k,k t k,k − it k,k dj k,k − i N n=k+1 T k,n dJ n,l , dD k,l = dq k,k + i k−1 n=1 dJ k,n Q n,l + idj k,k q k,k − iq k,k dh k,k − i N +ν n=k+1 Q k,n dH n,l . (C.23) The diagonal block entry reads: {dt k,k , dh k,k } {dj k,k , dq k,k } dC k,l 1 it k,k −it k,k 0 dD k,l 0 −iq k,k iq k,k 1 , (C. 24) with no dependence on variables further to the right of dT k,k and dQ k,k . iv) k < l; k, l = 1, . . . , N : dC k,l = dT k,l + i with dD k,l /dQ k,l = 1 , (C. 29) and no dependence on variables to the right of dQ k,l .
leading to the matrix det ∂(a, d, b, c, a ′ , d ′ , b ′ , c ′ ) ∂(dx k , dx * k , ds k , ds * k , dy k , dy * k , dp k , dp * k ) This result is obtained similarly to the previous case, with the contribution to the Jacobian now reading It is easy to see that the remaining cases iv) and v) simply give a contribution unity. Collecting eqs. (C.32), (C. 35) and (C.38) we obtain as a final result the conjectured Jacobian in eq. (5.6).