Symmetry Crossover Protecting Chirality in Dirac Spectra

We consider a random matrix model in the hard edge limit (local spectral statistics at the origin in the limit of large matrix size) which interpolates between the Gaussian unitary ensemble (GUE) and the chiral Gaussian unitary ensemble (chGUE). We show that this model is equivalent to the low-energy limit of certain QCD-like theories in the epsilon-regime. Moreover, we present a detailed derivation of the microscopic level density as well as the partially quenched and unquenched partition functions. Some of these results have been announced in a former letter by us. Our derivation relies on the supersymmetry method and is performed here step by step. Additionally, we compute the chiral condensate and the pion condensate for the quenched as well as unquenched settings. We also investigate the limits to GUE and chGUE and confirm our conjecture that the non-uniformity of the GUE limit would carry over to the hard edge limit.


Introduction
Random matrix theory (RMT) is an extremely versatile tool in the statistical description of spectra in physical systems [1]. This is especially true in QCD-like systems where RMT has been applied since the early 90's. Verbaarschot et al. [2][3][4] have shown that non-linear sigma models emerge from RMT in the large-N limit as a low-energy effective theory. Such sigma models exactly match those that arise in the ε-regime of QCD under appropriate conditions. For instance, the Compton wavelength of the Nambu-Goldstone modes have to exceed the system size. Then the path integral is dominated by their zero-momentum modes [5,6] and the contribution from the kinetic term factorizes.
Usually the terms in the QCD chiral effective theory can be derived by spurion analysis and invoking local flavour symmetries as in [7,8]. In the same way one can create the corresponding random matrix models. To each quantity like quark mass, chemical potential or lattice spacing one can introduce a dimensionless counterpart in RMT. The ε-regime of the partition function, then, uniquely fixes and identifies both sets, the physical variables in QCD and the dimensionless variables in RMT. This way one can derive an infinite number of spectral sum rules for the QCD Dirac operator along the lines of [6] in the QCD vacuum and [9][10][11] at finite chemical potential, as long as the matrix model is in the same universality class as the considered physical QCD system.
We consider a chiral random matrix model that interpolates between the Gaussian unitary ensemble (GUE) and the chiral Gaussian unitary ensemble (chGUE) statistics. This model has been first proposed in [12] to describe the continuum limit of three dimensional staggered fermions. This is one of three possible applications which we discuss in detail in Sec. 2. Beside this application, we also point out the possible use of our model to 3d continuum QCD with isospin chemical potential and to 4d gauge theories at high temperature. Furthermore, the considered model is also related to quantum chaos [13,14] where the complex eigenvalues of the off-diagonal matrix block of the chiral matrix have been studied. In this topic the model is known as the elliptic complex Ginibre ensemble [13,[15][16][17][18]. In this paper we study its singular value statistics and, thus, a different aspect of this model.
Apart from analysing why our model might fit to QCD-like theories, we derive the low-energy effective partition function in Sec. 3, which is known as the hard edge scaling limit in RMT, that has been given by us in the letter [19]. By exploiting the fact that the considered random matrix model exhibits a Pfaffian point process [20], we concentrate on the partition functions of one and two flavours, either bosonic or fermionic, in Sec. 3.1, because they are the building blocks of any spectral correlation.
Another quantity, which we have already reported without proof in the letter [19], is the miscroscopic level density at the origin. Its derivation is outlined in Sec. 4.3, and the lengthy details of the calculation are given in Appendix A.4. These computations are based on the supersymmetry method and we refer the reader to [21,22] for a pedagogical introduction. Note that this approach is different to the supersymmetric spurion analysis in [23,24]. In this section we also analyse the limits to GUE and chGUE and identify quantities which seem to be ideal to measure some low energy constants. Moreover, we study a quantity in Sec. 5, that resembles the chiral condensate in 4d continuum QCD. Indeed it can be identified with the chiral condensate when considering the applications of 3d staggered fermions or 4d QCD with twisted boundary conditions. However for 3d continuum QCD with isospin chemical potential it is essentially the pion condensate, see Sec. 2.3.
Additionally to the sections pointed out above, we summarize our results in Sec. 6 and give details to several computations in Appendix A.

Motivations of the Model
We are interested in the spectral statistics of the chiral random matrix [12,19,20] D = 0 iW iW † 0 , W ≡ H 1 + iµH 2 , H 1 , H 2 ∈ Herm(N ) and µ ∈ R (2.1) drawn from the distribution Tr (H 2 1 + H 2 2 ) . (2. 2) The set of Hermitian N × N matrices is denoted by Herm(N ). Our analysis starts with the partition function of N f quarks, with masses m f and additional source variables j a (f = 1, . . . , N f and a = 1, 2, 3). The source variables are helpful for calculating the observables in section 5. The Pauli matrices τ j are embedded in the 2N -dimensional space as 1 N ⊗ τ j though we omit the tensor notation.
When considering the spectral statistics of the complex matrix W only, our model is also known as the complex elliptic ensemble [13,[15][16][17]. The spectrum of W is generically complex and its support is given by an ellipse for large N , thus the name. The complex eigenvalues play an important role in the scattering at disordered or chaotic systems [13][14][15].
Let us point out another model which interpolates between GUE and chGUE. It is of the form D 5 = 0 W W † 0 + µH, W ∈ C N ×N and H ∈ Herm(2N ). (2.4) This model describes the Hermitian Wilson-Dirac operator [26], see also [27] for an equivalent realization. There is a crucial difference between the models (2.4) and (2.1). While the eigenvalues of D come always in chiral pairs ±iΛ n with Λ n ≥ 0 (n = 1, . . . , N ), it is only the case for D 5 when µ = 0. As already explained in [19,20] this difference is crucial for the behaviour of the eigenvalues close to the origin. When we have chiral pairs of eigenvalues we find a level repulsion from the origin regardless of how small µ is. This behaviour carries over to the microscopic limit about the origin when taking N → ∞ as we will see below.
In particular we are interested in the scaling µ 2 ∝ 1/N because it is where the symmetry crossover sets in. For the scaling µ 2 ∝ 1 we will see that we always obtain the statistics of chGUE.
Before we come to the spectral statistics in the limit N → ∞, let us point out some applications of this random matrix model in QCD. In Sec. 2.1 we discuss the relation to staggered fermions and their continuum limit. Also the high temperature limit of some quantum field theories can be modelled by D, see Sec. 2.2. The third application is presented in Sec. 2.3 and deals with 3d QCD at finite isospin chemical potential.

Symmetries and the Continuum Limit
Lattice Dirac operators of naive and staggered fermions [28] do not necessarily satisfy the same global symmetries as the continuum Dirac operator. Several lattice simulations [12,29,30] have shown this in three and four dimensions for different numbers of colors and fermions in the fundamental or adjoint representation. The good thing is that QCD in the Standard Model (4 dimensions, fermions in the fundamental representation and three colors) does not suffer from this behaviour. Only the rooting and the lack of a well-defined topological charge are problematic in this particular theory, which we will not address here. An attempt to understand this shift of symmetries in detail has been made in twodimensional QCD-like theories [31]. A detailed discussion of this phenomenon for gauge theories in general space-time dimension d ≥ 2 was done in [32] where a Bott-periodicity was revealed, not only in the space-time dimension but also in the number of directions with an even number (partition) of lattice sites. Let us briefly recount the situation.
For each direction, which has an even parity of lattice sites, a "chirality" operator Γ j (Γ j = Γ † j = Γ −1 j and Tr Γ j = 0) can be defined. This operator assigns to an odd-lattice site (counted only in this direction) a "+1" and to an even lattice site a "−1". As a consequence, this operator anti-commutes with the naive Dirac operator D and, hence, generates an additional chiral symmetry. The collection of all these chiral operators Γ j build a Clifford algebra, i.e. [Γ k , Γ l ] + = 2δ kl with [., .] + the anti-commutator. This has two consequences. First and foremost, the naive Dirac operator may become (highly) degenerate. Moreover, the Dirac operator changes its global symmetries along the Bott periodicity [33].
What does this mean for staggered fermions? Staggered fermions are essentially naive fermions with an even parity in each direction. Due to the Bott periodic shift of global symmetries, it was shown in [31] that staggered Dirac operators share always the global symmetries of the corresponding eight-dimensional continuum theory. This explains why QCD in the Standard Model does not suffer from this problem because the global symmetries of the four-dimensional and eight-dimensional continuum theory are the same. This is not true for 3d QCD as well as in any other odd dimension. Staggered fermions of QCD with the gauge group SU(N c ≥ 3) in the fundamental representation yield always a chiral, complex, anti-Hermitian Dirac operator and, thus, shares the symmetries with the eightdimensional and, thence, four-dimensional continuum theory. The question is whether the correct global symmetries are recovered when the continuum limit is taken.

Matrix Model for the Symmetry Shift
In [12] a mechanism of such a change from symmetries of even to odd dimensional SU(N c ≥ 3) gauge theory in the fundamental representation was proposed. They considered the model (2.1) and fitted Monte-Carlo simulations of this random matrix model (2.1) to lattice simulations for three-dimensional staggered fermions in the quenched theory for several coupling constants and lattice sizes. The comparison seems to look surprisingly good despite the fact that the degeneracy of the eigenvalues does not fit with the number of doubler fermions of unrooted staggered fermions. Without rooting the number of flavours of staggered fermions should be enhanced from N f to 4N f in three dimensions. However the Dirac operator (2.1) is for µ = 0 only doubly degenerate. We underline that the continuum limit is given by µ = 0 in this model. A more appropriate model for the described situation would be , H ∈ Herm(N ), W ∈ C 2N ×2N and µ ∈ R, (2.5) which obviously leads to more terms in the chiral Lagrangian compared to the one in the partition function (3.6) derived for the model (2.1). We briefly show this, here. For this purpose, we want to consider the partition function into a sum of Hermitian N × N matrices W R,j and W I,j with χ 0 = 1 2 and χ 1,2,3 the three Pauli matrices. The distribution is chosen to be Gaussian (2.8) The choice of different standard deviations c R,j , c I,j takes into account that the lattice is not necessarily invariant when interchanging its axes. Thus, our model is more in the spirit of [34] where a random matrix model for 4d staggered fermions has been proposed. Next, we introduce an N × 4N f rectangular matrix V whose matrix entries are independent complex Grassmann variables with a complex conjugation of the second kind, i.e. [35] for an introduction to superanalysis and superalgebra. Then, we can rewrite the product of determinants as a Gaussian integral over V , i.e.
The Pauli matrices τ j originate from the chiral structure in D stag and are embedded in the 4N f -dimensional space as a tensor τ j ≡ 1 N f ⊗ τ j ⊗ 1 2 ; similarly we embed M ≡ M ⊗ 1 2 ⊗ 1 2 and χ j ≡ 1 N f ⊗ 1 2 ⊗ χ j . The average over the random matrices yields where in the last line we employed the bosonization formula [36][37][38] and replaced V † V → √ N U . The measure dµ(U ) is the Haar measure of the unitary group U(4N f ).
In the end we take N → ∞ while keeping M = √ N M and µ 2 = N µ 2 fixed, which tells us that the integral concentrates on the manifold U = U τ 3 U † with U ∈ U(4N f ). See the detailed discussion in Sec. 3.1 since the saddle point equation is exactly the same. Thence, we end up with (2.11) Our model only describes the first two terms of the chiral Lagrangian, cf. (3.6). The inclusion of all terms makes the model certainly more sophisticated but, as a drawback, it becomes also more analytically involved, in particular the finite N discussion [20] would have not been possible anymore. Nonetheless, the good qualitative agreement [12] of the spectral statistics of the Dirac operators of the model (2.1) and of three-dimensional staggered fermions shows that the random matrix model (2.1) is worthwhile to study. In particular one might conjecture, that the missing terms proportional to c 2 R,j>0 and c 2 I,j>0 in the chiral Lagrangian play a sub-leading role, at least in the quenched theory. What happens for the unquenched theory has to be still perused.

Gauge theories at high temperature
Euclidean quantum field theories generally undergo dimensional reduction at high temperature T = 1/β [39][40][41][42][43][44][45] where β is the circumference of S 1 . A rough perturbative picture of this phenomenon is as follows. Let us consider QCD-like theories with Dirac fermions ψ on spacetime R 3 ×S 1 . The temporal part of the Lagrangian of ψ reads ψγ 4 [∂ 4 + iA 4 (x)]ψ. If we substitute A 4 (x) by a constant A 4 via a mean field approximation, the eigenvalues of ∂ 4 + iA 4 are given by i[(2nπ + θ)/β + A 4 ] with n ∈ Z, where θ specifies the boundary condition, ψ(x 4 + β) = e iθ ψ(x 4 ). In a thermal phase with a trivial Polyakov loop P = 1, one can set θ = π and A 4 = 0 so that the smallest (in magnitude) eigenvalue is π/β. Therefore in the limit π/β Λ QCD fermions decouple from the low-energy dynamics. In particular, the chiral condensate evaporates and chiral symmetry is restored, inhibiting applications of chiral random matrix theory to this hot phase.
What happens in other phases is more interesting: when the mass gap in units of T , which is min n∈Z |2nπ + θ + β A 4 |, is small, the chiral symmetry breaking tends to persist up to higher T . Let us consider a pure gauge theory, where the situation is simplest. Recall that in hot SU(N c ) pure gauge theory there are N c distinct vacua having Tr P ∝ e iφ with φN c /2π ∈ Z. It has been known that ψψ of valence quarks obeying θ = π exhibits a strong dependence on φ [46][47][48]. In particular, when N c is even, ψψ seems to remain nonzero up to arbitrarily high T in the vacuum with φ = π. Indeed, substituting θ = π and A 4 = π/β yields min n∈Z |2nπ +θ +β A 4 | = 0, implying that quarks acquire no perturbative mass gap at high T . 1 Such "high-temperature chiral symmetry breaking" can also be explained on the basis of topological excitations of the gauge field, called instanton-monopoles or dyons, that carry fermion zero modes due to the index theorem [49][50][51]. We conjecture that the near-zero region of the quenched Dirac spectrum in this phase would undergo a dimensional crossover from chGUE in d = 4 to non-chiral GUE in d = 3. This transition can be studied via the model (2.3) with N f = 0. In this case the parameter µ signifies the effective size of the fourth dimension. This conjecture can be tested in Monte Carlo simulations. We emphasize that chiral symmetry of the Dirac operator is not explicitly broken throughout the symmetry crossover.
Whether such a smooth dimensional crossover is possible or not in the presence of dynamical quarks is a highly nontrivial issue. We refer to [52,53] for recent works on this subject, which were inspired by the idea of adiabatic continuity in [54][55][56].

3d QCD at Finite Isospin Chemical Potential
Let us consider again QCD in three dimensions with the gauge group SU(N c ≥ 3) in the fundamental representation. For two flavours the ground state will accommodate a pion condensate ud = 0 at T = 0 when the isospin chemical potential µ I is large enough, that entails a spontaneous breaking of the U(1) isospin symmetry, see [11,[57][58][59][60] for a similar discussion in four dimensions.
The fermionic part of the Euclidean Lagrangian with a source term j π for the pion condensate (similar to the diquark source term in two-color QCD [11,[61][62][63]) is given by with the covariant derivative D ν = ∂ ν + iA ν (ν = 1, 2, 3), the gauge vector field A µ ∈ su(N c ≥ 3) and the two quark fields V = (u, d). In this section, we use Einstein's summation convention. Let us recall that the x 3 -direction is the imaginary time-direction in the 3 dimtheory. Furthermore the Pauli matrices σ ν act on the spinor space while the Pauli matrices τ l act on the flavour index. 1 The same phenomenon occurs for any Nc if we instead impose θ = 0 (i.e., periodic boundary condition) and look at the sector with φ = 0. For more details on the interplay of chiral symmetry breaking and the fermionic boundary condition, we refer the reader to [64][65][66][67].
Taking the derivative of Eq. (2.13) with respect to j π at j π = 0 and averaging over the gauge field configurations, one obtains a Banks-Casher-type relation that links the pion condensate to the density of the smallest singular values of σ ν D ν + m q + µσ 3 [11,60,68]. Note that one has to set m u = −m d = m in order to derive the Banks-Casher-type relation. Statistical fluctuations of these singular values can be analyzed with the model (2.1) for one flavour N f = 1. To see this we briefly show that the chiral Lagrangian of the physical system in the ε-regime and with degenerate quark masses agrees with the one following from the random matrix model (2.1).
After the quark fields are integrated out, the partition function comprises the determinant Comparing with the random matrix model (2.1) one can identify σ ν D ν ↔ iH 1 , σ 3 ↔ −H 2 and j π ↔ m f for N f = 1. Note that j π plays the role of the mass here. The random matrix model, which can be naturally associated with this determinant, can be chosen as with H a 2N × 2N dimensional Hermitian matrix. Let us underline that the quantities j and µ are not equal to their physical counterparts in (2.12) but need a rescaling with the low energy constants in the ε-regime. The same also holds for the quark masses M q , which will be represented by the dimensionless diagonal 2 × 2 matrix M q in the random matrix setting.
We now proceed as in Sec. 2.1 and derive the corresponding chiral Lagrangian in the ε-regime, which corresponds to the partition function (2.15) Anew we introduce a Grassmann valued rectangular matrix V , which is of dimension 2N ×2, to rewrite the determinant as a Gaussian integral (2.16) The integration over the random matrix H can now be readily done and we obtain (2.17) A Hubbard-Stratonovich transformation [69,70] has been applied in the second step, where we introduced the integral over the 2 × 2 Hermitian matrix Q, and the Grassmann variables were integrated out in the last one. When keeping j = √ N j, M q = √ N M q , and µ = µ fixed, we find the saddle point equation Q 2 = 1 2 in the limit N → ∞. The solutions Q = ±1 2 are algebraically suppressed compared to Q = U τ 3 U † with U ∈ U(2). Thence, the asymptotics of the partition function (2.17) is This is the same result when applying the spurion analysis for the physical QCD-model (2.12) where the chiral Lagrangian at the leading order is equal to Σ is the condensate uu − dd , F is the pion decay constant, and W is another low energy constant related to the pion condensate ud . Comparing (2.17) and (2.19), we notice that the dimensionless quantities from the random matrix model are given by the physical quantities as 2 M q = V ΣM q , 2 j = V W j π , and µ 2 = V F 2 µ 2 iso /2, that entails also the physical scaling M q , j π , µ 2 iso ∼ 1/V . For degenerate quark masses m u = m d the chiral Lagrangian reduces to the one, which results from the random matrix model (2.5) considered by us, see Sec. 3.1, though, the case m u = −m d for the Banks-Casher-type relation [11,60,68] can be readily obtained from a slight modification of our calculations in Appendix A.2.
Indeed, at some point one sees the difference between the models (2.1) and (2.14). When choosing the parameter µ ∼ √ N on the level of random matrix theory or V µ 2 iso ∝ V the spectrum of D iso develops a spectral gap about the origin, whereas D of (2.1) does not. Thus, the relation between D and D iso is similar to the relation of the Osborn model [71] and the Stephanov model [72] for the baryon chemical potential.

Effective Lagrangians of the Pseudo-Scalar Mesons
In Sec. 3.1 we derive the effective Lagrangian of the model (2.1) in the hard edge limit, which is the RMT counterpart of the ε-regime in QCD. In this way we establish the intimate connection between our model and the applications delineated in Sec. 2. As a cross check we study the limit to the GUE (3d QCD) and the chGUE (4d QCD) result in Sec. 3.2. Furthermore we derive the macroscopic level density of the model (2.1) in Sec. 3.3, to underline that one has to be careful in comparing the hard edge results of RMT with Monte Carlo simulations at finite matrix dimension N .

Effective Lagrangian
The connection of RMT with QCD is given via the non-linear σ-model which is the chiral Lagrangian in the 4-dimensional continuum QCD. The reason why these two very different theories are related is due to same spontaneous breaking of global symmetries. Thus we consider the partition function (2.3) for N f ≥ 1 in the large-N limit. We first map this partition function to a dual matrix space whose dimension only depends on N f . This method is called the supersymmetry method and introductions can be found in [21,22]. In Secs. 2.1 and 2.3, we have shown two kinds of its procedure, namely the superbosonization and the Hubbard-Stratonovich approach, respectively. Here, we pursue the superbosonization approach [36][37][38].
In the first step, we introduce a rectangular matrix V of dimension N × 2N f whose matrix entries consist of independent complex Grassmann variables. We rewrite the product of determinants as The denominator on the right hand side correctly normalizes the Gaussian integral over V . Let us emphasize that M and the Pauli matrices τ j act on two different components of a tensor space; we embed them as M ≡ M ⊗ 1 2 and τ j ≡ 1 N f ⊗ τ j though this is an abuse of the notation.
The average over H 1 and H 2 in Eq. (3.1) can be readily performed and yields .

(3.3)
We recall that we get additional signs in the exponential function due to the anti-commuting nature of the matrix entries of V . Now we are ready to apply the superbosonization formula [36][37][38] and replace the nilpotent The scaling factor √ N is needed to perform the saddle point analysis. As a price of this exchange, we not only replace the flat Berezin measure dV by the Haar measure dµ(U ) on U(2N f ), but also get a factor det −N U , as it has been the case in Eq. (2.10). The additional term reflects the nature of the integral over the Grassmann valued matrix V , that picks out only the highest term of a Taylor expansion in V † V . The group integral with the term det −N U is a multidimensional contour integral selecting the correct terms of this Taylor expansion. This way we arrive at where the denominator is only a constant depending on N and N f but nothing else. The double scaling we are looking into is given by M = √ N M and µ = √ N µ fixed when N → ∞. This scaling is obviously different from the Stephanov-type model (2.14) and originates from the fact that the parameter µ is the prefactor of a random matrix H 2 with fully occupied matrix entries, cf. (2.1), while for the model (2.14) µ stands in front of a fixed diagonal matrix.
To get a finite result in this limit, we need to renormalize the partition function by a factor C N,N f that only depends on N and The explicit form of C N,N f is irrelevant for physics, because it depends on the random matrix model and is, hence, not universal. We assume that it is chosen such that the resulting partition function is given by an integral with the normalized Haar measure, cf. (3.6). We also absorb the denominator in Eq. (3.4) in this constant.
Taking N → ∞ we perform a saddle point approximation. We first have to solve the saddle point equation Hence the eigenvalues of U τ 1 are ±1. Diagonalizing U τ 1 yields a Jacobian proportional to the modulus square of the Vandermonde determinant of the eigenvalues of U τ 1 . This Vandermonde determinant implies that all solutions with Tr U τ 1 = 0 are algebraically suppressed by factors of 1/ √ N to those solutions which have an equal number of eigenvalues +1 and −1. All of these contributing solutions are unitarily equivalent such that the general solution is given by When plugging this result into Eq. (3.4), we finally arrive at the effective partition function cf. Eqs. (2.11) and (2.19). This equation is the main result of this section. Actually, the integration is effectively only over the coset It is the manifold for the Nambu-Goldstone bosons of the spontaneous symmetry breaking

Limits of the Effective Lagrangians
There are two interesting limits of the effective partition function (3.6). First and foremost, we can take µ → 0 followed by a π/4-rotation U → e −iπτ 2 /4 U . Then, the effective partition function is given by This is exactly the finite volume partition function of QCD in three dimensions [4,73,74]. The second limit is given by µ 1, which is slightly more involved. The ground state corresponds to the minimum of Tr (U τ 3 U † τ 3 ) 2 . To determine it, we consider the Hermitian The equality holds if and only if U τ 3 U † τ 3 = −τ 3 U τ 3 U † . Hence U τ 3 U † has to have a chiral form. In addition, it is unitary and Hermitian, i.e.
We plug this into Eq. (3.6) and find This is the conventional ε-regime partition function of QCD in four dimensions [5,6]. As schematically shown in Fig. 1, the parameter µ = N µ 2 effectively controls the symmetry in the low-energy theory. The coset manifold that represents the Nambu-Goldstone field also evolves with µ = N µ 2 accordingly; with increasing µ 2 = N µ 2 , some of the Nambu-Goldstone modes acquire a mass gap and gradually decouple from the low-energy physics.
which are the two Nambu-Goldstone manifolds corresponding to the two limits. We emphasize that the symmetry crossover in Fig. 1 is different from the one realized by Wilson fermions [26,27], which break chiral symmetry explicitly.

Macroscopic Level Density
Let us turn to the macroscopic level density despite the lack of direct connection to continuum QCD. The reason why we want to calculate it comes from an observation while we have performed Monte Carlo simulations of the model (2.1) (see Ref. [19] as well as Sec. 4.3) that the magnitude of the level density at the origin varies with µ when µ becomes large. Hence, one has to be careful when simulating the random matrix model (2.1) and comparing it with our results in the hard edge limit (microscopic spectral statistics about the origin).
To understand this behaviour let us take up the fermionic partition function (3.1) for N f = 1 and take the limit N → ∞ at fixed M/ √ N = ix − = ix + and µ. It is known that the macroscopic level density R 1 (x) is obtained as Note that the order of the limits and the derivative is important to select the correct solution. Moreover, we want to emphasize that this relation only works for the global spectral statistics. The calculation of the microscopic spectral density will be postponed until Sec. 4.3. The quotient Z is the skew-orthorgonal polynomial corresponding to the random matrix ensemble (2.1), see Ref. [20] where the following formula has been derived Thus the level density is proportional to For large N one has to solve the saddle point equation which yields the two solutions with sign(x) = x/|x| the sign of x, survives the limit N → ∞ due to the finite increment > 0 so that the macroscopic level density is or after proper normalization we find the Wigner semi-circle, of radius 2 1 + µ 2 . It is this µ-dependence of the radius which changes the height of the level density about the origin which is 1/[π 1 + µ 2 ]. The relation between the macroscopic level density R 1 (x) of the eigenvalues and the microscopic level density ρ(λ) follows from the relation x = λ/N between the global and local coordinates of the eigenvalues: We want to point out that quite often the relation between the two regimes is also chosen to be x = λ/(2N ), meaning divided by the whole matrix dimension of D, which leads to the 1/π asymptotics for chGUE , see [3,75].
As a conclusion from the above discussion, in the hard edge limit, which is the one related to the ε-regime of QCD, the coupling parameter µ scales as 1/ √ N and, hence, the asymptotic height is lim λ→∞ ρ(λ) = 2/π. On the other hand, when µ is large enough we have to take into account a correction when comparing the microscopic limit results in Sec. 4.3 with finite N Monte Carlo simulations. For example for µ = √ N µ = 2 with N = 50 it is still of about 5% and, thence, visible. Therefore one has to go to much larger matrix sizes when increasing the rescaled parameter µ; noting that N has to scale quadratically with µ.

Partition Function with Valence Quarks
Now we turn to partition functions with valence quarks, i.e.
The bosonic valence quarks have a non-vanishing real part, Re κ j,b = 0, to guarantee the integrability. Usually one sets k b = k f − N f = k and chooses the first N f masses κ f,j equal to the masses of the dynamical quarks and the remaining κ f,j and κ b,j being the valence quark masses which might be complex as it is the case for calculating the k-point correlation function.
In [20], we have shown for the model (2.1) at finite N that the whole spectral statistics are completely described by the partition function of one and two bosonic or/and fermionic flavours. The reason is that the spectral statistics of the singular values of (2.1) exhibit a Pfaffian point process, see [76] for the definition. This carries over to the hard edge limit and, hence, to the physical QCD systems summarized in Sec. 2. Therefore we especially concentrate on those quantities when giving explicit results.
In the following discussion, we pursue the same strategy as in Sec. 3.1, but extend the ideas to superspace. First we arrange the masses in a diagonal ( . . , κ f,k f ). Here we use the notation that the boson-boson block is in the upper left block of a supermatrix and the fermion-fermion block in the lower right block. For the preparation of our calculation, we additionally need to keep the integrability throughout the calculation.
In the first step we introduce a complex rectangular matrix V which is this time supersymmetric with dimensions N × (2k b |2k f ). The supersymmetric counterpart of the is Str σ = Tr σ bb − Tr σ ff and the related superdeterminant is given as Sdet σ = det(σ bb − σ bf σ −1 ff σ fb )/ det σ ff . Now we are able to integrate over the two Hermitian matrices H 1 and H 2 and find Here we have chosen the normalization of an integral over a Grassmann variable η as dη = 0 and ηdη = 1 and the order in its measure is dηdη * .
In the final step we replace the supermatrix V † V Lτ 1 by the supermatrix √ N U ∈ Gl C (2k b |2k f )/U(L bb τ 1 |2k f ) with the help of the superbosonization formula [36][37][38]. What does this explicitly mean, especially the first part of the notation U(L bb τ 1 |2k f )? The matrix U has explicitly the form with U ff ∈ U(2k f ) being an ordinary unitary matrix and U bb L bb τ 1 being a positive definite Hermitian matrix. Since the embedding of the non-compact, ordinary group U(k b , k b ) in the supergroup U(L bb τ 1 |2k f ) is non-trivially given via the matrix L bb τ 1 we have highlighted it in our notation. The off-diagonal block η is a 2k f ×2k b rectangular matrix with independent complex Grassmann variables as matrix entries and η † is its Hermitian adjoint. After exploiting the superbosonization formula, we arrive at The (non-normalized) Haar measure dμ(U ) can be expressed in terms of the flat measure dU (products of the differentials of independent matrix entries) as dμ(U ) = Sdet (U ) 2(k f −k b ) dU . The normalization constant is which can be checked for κ → ∞, where the partition function becomes Z Starting from Eq. (4.6) we perform the hard edge limit for particular cases. In Secs. 4.1 and 4.2 we consider the partition functions with either only fermionic or bosonic quarks, respectively. The partially quenched case is analyzed in Sec. 4.3 where we also derive the microscopic level density.
The case with only fermionic quarks has been already discussed in Sec. 3.1. Here, we keep track of all constants and give very explicit results for the one and two flavour partition function. Especially the latter two functions are sufficient to construct all fermionic partition functions. For instance for an even number N f ∈ 2N of flavours, it is [20] In the case of an odd N f one introduces an additional mass m N f +1 and applies Eq. (4.8) for N f + 1 masses. Eventually, one lets m N f +1 → ∞ yielding in the last row and column the one-flavour partition function. In appendix A.2 we will derive Eq. (4.8) directly from the large-N partition function (3.6).
Before we arrive at the hard edge result (3.6), we want to give a detailed list of contributions which result from the saddle point approximation. Let us take Eq. (4.6) with no bosonic degrees of freedom, k b = 0. In the hard edge limit, we expand the unitary matrix U as U = U zU † = U (τ 3 + diag(iδz 1 , −iδz 2 )/ √ N )U † with δz 1 and δz 2 two real diagonal k f -dimensional matrices parametrizing the massive modes and U ∈ U(2k f ) a unitary matrix distributed by the normalized Haar measure dµ(U ). The different signs in front of δz 1 and δz 2 result from the opposite crossings of the contours of the eigenvalues of U through the saddle points ±1. This change of variables yields the following approximation of the Vandermonde determinant: . The differentials transform as dz = N −k f dδz 1 dδz 2 . Additionally, we get the constant 1/(2k f )! 2k f −1 j=0 (π j /j!) that is essentially the volume of the coset where S 2k f is the symmetric group permuting the diagonal elements of z. The latter has to be compensated since we ordered the signs in front of δz, which produces an additional factor of (2k f )!/(k f !) 2 . Moreover, we obtain the sign (−1) N k f at the saddle points from the determinant Sdet N U = det −N U .
Collecting everything, we have for the fermionic partition function with an arbitrary k f with Z (k f ) given as in Eq. (3.6). In the particular cases of one and two flavours this reads and respectively, where we used Stirling's formula for the Gamma function in the limit of a large argument.
Comparing Eqs. (4.8) and (4.9), it is apparent that most of the group integral over U ∈ U(2k f ) can be performed, explicitly. In particular it yields a Pfaffian determinant. The remaining integrals in the one-and two-flavour partition function in Eqs. (4.10) and (4.11) are simplified further in Appendix A.2 so that we end up with and (4.13) These two results are reminiscent of the Bessel kernel [75], which is deformed now. Indeed the limit µ → ∞ can be readily checked yielding the Bessel-kernel for vanishing topological charge, because the integrals (4.12) and (4.13) are evaluated at the saddle points x = 0 and x 1 = x 2 = 0, respectively. Confirming the correct limit for µ → 0 in Eqs. (4.12) and (4.13) is not that simple though one can set µ = 0 in the integrand. To see that the partition functions become indeed those of the sine-kernel [25], one has to reverse the Berezin-Karpelevich integral (A.11) and integrate over U instead of U τ 3 U † in Eqs. (4.10) and (4.11) which is the Harish-Chandra-Itzykson-Zuber integral [77,78]. The result is the sine-kernel result which is for two or four flavours (masses always appear as chiral pairs ±κ j ), and , one can alternatively calculate this result more directly by employing the Taylor expansion of the Bessel function in Eq. (4.12) and then integrating each term separately.

The case
As for the fermionic partition functions, we want to determine the exact normalization constants as well as the partitions function of one and two bosonic flavours as those partition functions also satisfy [20] for an even number k b of bosonic flavours. The counterpart of this equation for odd k b can be obtained in the same way as for fermionic flavours via introducing an additional auxiliary bosonic quark with mass κ b,k b +1 and sending this mass to infinity in the end.
To begin with, we again calculate the components of the saddle point approximation. The hard edge limit works the same way as in the fermionic case, in particular the matrix U can have only the eigenvalues ±1 at the saddle point. However, we now have a restriction for the contours that is reflected by the fact that U L bb τ 1 is a positive definite Hermitian matrix. Hence, we can parametrize U and δz being a 2k b × 2k b matrix which satisfies the following conditions. The matrix δzL bb τ 1 is Hermitian and fulfills the commutation relation [δz, L bb τ 1 ] = 0. We want to emphasize that in the saddle point manifold we cannot arbitrarily permute the entries of L bb τ 1 because of the nature of the non-compact group U(L bb τ 1 ) 2 which is composed of disjoint parts. Therefore a combinatorial factor does not appear this time, the saddle point is unique. In spite of this, for all choices of L bb the noncompact group U(L bb τ 1 ) is unitarily equivalent to the non-compact unitary group U(k b , k b ) because Tr L bb τ 1 = 0. More precisely, we can bring all matrices in U(L bb τ 1 ) to a standard form of matrices in U(k b , k b ) by a single unitary transformation.
As U(L bb τ 1 ) is a non-compact group, there is no normalizable Haar measure. Consequently, the normalization constant cannot be computed in the traditional way with the volumes of the groups, but we have to stick with the measure which is given by the pseudo-Riemannian length element, g(dU , dU ) = Re Tr (dU ) 2 = 1 N Re Tr δz 2 + Re Tr U −1 dU, L bb τ 1 + δz √ N 2 ≈ 1 N Re Tr δz 2 + Re Tr U −1 dU, L bb τ 1 2 . (4.17) We denote the volume element resulting from the invariant length element by dμ(U ) which is given by the standard formula dμ(U ) = √ det g dU . Thence, we have for the flat measure The measure dμ(U ) is the pseudo-Riemannian volume element of the invariant length element Re Tr U −1 dU, L bb τ 1 2 .
We put again everything together and find for a general bosonic partition function This intermediate result is more suitable after performing a π/4-rotation, in particular (4.20) Then the 2k b × 2k b matrix U changes its symmetries to U −1 = τ 2 U † τ 2 so that it becomes independent of L bb . The cases of one and two flavours are given by As before the integral (4.20) over the coset U(1 k b ⊗τ 2 )/[U(k b )×U(k b )] can be evaluated explicitly, even when it is of a non-compact type, because Eqs. (4.16) and (4.20) must agree. Therefore, also this integral yields a Pfaffian determinant whose matrix entries are essentially given by Eqs. (4.21) and (4.22). In Appendix A.3 we simplify these two partition functions, which are explicitly

The case
The quenched partition function of one bosonic and one fermionic flavour is equal to where we have defined the Lagrangian This function should not be confused with the diagonal matrix L consisting of signs defined right before Eq. (4.2). The calculation of the large N -limit is quite lengthy and is deferred to Appendix A.4. Here we state the final result, Despite this cumbersome result, the chGUE result can be readily regained in the limit µ → ∞. Then, only the second and third term contribute because the integration variables scale like x, y ∝ 1/ µ. Therefore the quenched partition function becomes which indeed agrees with the result of chGUE. For µ = 0, the integrals can be performed explicitly as well, Its simple form in this limit is again reminiscent of the sine-kernel structure of the GUE. After we gained the quenched partition function (4.29) we can calculate the microscopic level density by differentiating in κ f , thereafter setting κ b = κ f = L − iλ, and eventually taking the real part in the limit → 0, i.e. (4.32) For this purpose we employ the identities [81] for λ > 0. Therefore the microscopic level density of the quenched system is This is the main result of this section and has been published in the letter [19], but without the detailed calculation in Appendix (A.4). 3 The microscopic level density (4.34) has also been compared with Monte Carlo simulations of the random matrix model (2.1) in [19]. It has been very surprising for us how fast it converges to the chGUE limit [75] ( µ → ∞), since it is almost perfectly achieved for µ = 2, see [19] as well as the right plot in Fig. 2. In contrast, the level density in the GUE limit ( µ → 0) has always a small but persistent peak close to the origin, see [19], while it should be only a constant lim µ→0 ρ(λ) = 2/π. This has been already seen for the quenched staggered Dirac operator in three dimensions [12]. To understand this non-uniform behaviour we have to scale the spectrum on the scale µ where this peak lives. Thence we rescale λ → µλ in Eq. (4.34) and take the limit µ → 0. Before we can commute the limit with the integrals we have to rescale the integration variable y → y/ µ, as well, i.e. × 4y 4 + y 2 − λ 2 y 2 − 1 4 |λ|J 0 (2λy) + 4y 2 + 3 2 λ 2 yJ 1 (2|λ|y) We show the behaviour of this functions as well as how it is approached by the random matrix model (2.1) in the left plot of Fig. 2. The limit seems to be almost perfectly achieved for µ < 0.1. The persistent peak can be interpreted as the smallest singular value which is enforced to lie on the positive real line. The eigenvalue of the GUE closest to the origin has not this restriction, yet the smallest singular value of the model (2.1) shall become exactly this eigenvalue. This obvious conflict is the origin of the non-uniform behaviour of the level density in the limit µ → 0.

Chiral Condensate/Pion Condensate
In this section we analyze the quantity where the prefactor √ N correctly rescales the eigenvalues for the hard edge scaling limit. This quantity is either proportional to the chiral condensate Σ( κ = m) = uu−dd with m = V Σm the dimensionless quark mass (either of a sea or a valence quark) in the application to the staggered Dirac operator in three dimensions, see Sec. 2.1, as well as in particular high temperature four dimensional quantum field theories, see Sec. 2.2, or it is proportional to the pion condensate ud with κ = j = V W j π with j the corresponding dimensionless source variable in 3d QCD at finite isospin chemical potential as discussed in Sec. 2.3. Since the computation of O( κ) differs between the quenched theory and a theory with dynamical quarks, we discuss it in separate Secs. 5.1 and 5.2.
One last general comment on O( κ): the asymptotics for large argument is determined by the asymptotics of ρ(λ) which is lim λ→∞ ρ(λ) = 2/π. Hence the limit for O( κ) is This behaviour is also true with dynamical quarks because our random matrix D only models the topologically trivial phase (topological charge is ν = 0).

Quenched theory
In the quenched case the observable (5.1) can be computed via a derivative of the quenched partition function Z After plugging the result (4.29) into this equation, we find Interestingly, both limits are reached extremely quickly. For µ → 0 differences can be hardly seen at µ = 0.1, whilst for µ → ∞ it is already in good agreement with µ = 0.9. Thus the convergence is even faster than that for the microscopic level density (cf. Sec. 4.3), which can be actually expected because integrals usually show an improved convergence. The impressive convergence for "large" µ can be understood by the Gaussian factor in the integral (5.4) that suppresses very rapidly any deviation from the saddle point x = y = 0. For "small" µ it helps that everything only depends on its square. Therefore the corrections of µ = 0.1 is only of about 1%.

One and Two flavours
When we have dynamical quarks we can set the auxiliary variable κ equal to one of the quark masses, say m 1 , so that we can express it in terms of a derivative of the partition function Z (N f ) , in particular we have Inserting the results (4.12) and (4.13) into this relation we obtain for one and two flavours and , (5.9) respectively, where the dependence of O (N f =2) on m 2 is made implicit for brevity. The function for the two-flavour case is (5.10) Applying the derivative ∂ m 1 in the two-flavour case explicitly yields a lengthy expression. This statement is even true for its limits µ → 0 and µ → ∞; therefore we omit these expressions.
The behaviour for the one-flavour case is shown in Fig. 4(a). It is apparent that the GUE limit lim as well as the chGUE limit are uniformly approached and, surprisingly, well-achieved for moderate values of µ, namely µ = 0.5 and µ = 5, respectively. This has to be seen in contrast to the quenched theory in Sec. 5.1 where no uniform limit to GUE exists. The dynamical quark mass pushes the spectrum away from the origin so that the region on the scale µ is almost void of eigenvalues. This behaviour can be also seen for the two-flavour case, which can be observed for instance in Fig. 4(c) where we have chosen m 2 . The slope at the origin of O( m 1 ) flattens out when the second mass decreases, cf. Fig. 4(b) for vanishing mass m 2 . Considering the duality of flavour and topological charge, the case m 2 = 0 can be identified with the ν = 1 case of chGUE. Then the observable (5.13) -25 - This expression has the asymptotic behaviour for the chGUE limit and for the GUE limit. Also illustrated in Fig. 4 Again the results for the values µ = 5 and µ = 0.5 almost perfectly match the asymptotics. As for the quenched theory, the Gaussian term in the integrals (5.8) and (5.9) as well as the dependence on µ 2 help the convergence. Furthermore, the repulsion from the dynamical quark masses shifts the region where we have a good a agreement to larger values of µ. This can be understood from the analytical results (5.8) and (5.9) because only modified Bessel functions of the first kind are involved which grow exponentially. In the quenched case (5.4), those exponential contributions are compensated by the modified Bessel functions of the second kind.

Conclusions
We continued our study [19,20] of the spectral statistics of a Gaussian random matrix model that interpolates between the GUE and the chGUE while preserving chirality. In the present work we investigated the statistics on the scale of the local mean level spacing about the origin (hard edge scaling) which agrees with the Dirac spectra in several QCD-like theories in the ε-regime. We gave a detailed account of approximations involved in showing such correspondences.
Additionally, we gave full details of the derivation of the microscopic level density (4.34), that we already presented in the letter [19]. We even carried further the analysis of the level density and identified a scale-free function (4.36) in the limit of vanishing coupling constant µ → 0. This limit is a persistent deviation from the level density of the GUE which is a constant on the local scale. This deviation results from the protected chirality and essentially reflects the fact that the smallest singular value λ > 0 always feels a residual interaction with its "mirror charge" −λ < 0. The scale-free function (4.36) shows a single peak that can be interpreted as this singular value and has also been observed for the staggered Dirac operator in three-dimensional QCD. Therefore we are confident that our results may help in fixing the low energy constant corresponding to the term Tr (U τ 3 U † τ 3 ) 2 in the effective Lagrangian of the Nambu-Goldstone bosons, cf. Eq. (2.11).
Moreover, we anlytically evaluated an observable (5.1) which corresponds to quark bilinear condensates in QCD-like theories. We derived explicit expressions of this quantity for the quenched theory and the unquenched theory with one and two flavours. We observed that in the presence of dynamical quarks the spectrum is pushed away from the origin so that both limits, GUE ( µ → 0) and chGUE ( µ → ∞), are approached uniformly in the unquenched theory. In contrast, the quenched theory exhibits a non-uniform asymptotics to the GUE result. As it has been the case for the quenched microscopic level density, one can make out a scale free function (5.6) on the scale µ also for this observable.

Acknowledgments
MK acknowledges support by the German research council (DFG) via the CRC 1283: "Taming uncertainty and profiting from randomness and low regularity in analysis, stochastics and their applications".

A Details of Some Calculations
In this appendix we show detailed computations of several results presented in the main text. The normalization (4.7) of the supersymmetric integral (4.6) is derived in Appendix A.1. The one-and two-flavour partition functions for fermions and bosons are evaluated in Appendices A.2 and A.3, respectively. In these two sections we also consider the partition functions with an arbitrary number of flavours and show how the Pfaffian forms (4.8) and (4.16) can be directly obtained from the large N results (3.6) and (4.20). In Appendix A.4 we explain the main steps of the calculation from Eq. (4.27) to Eq. (4.29).

A.1 Derivation of the Normalization (4.7)
To obtain the second line of Eq. (4.7), we first integrate over the complex Grassmann variables η which only appear in the measure and the superdeterminant: In the first step we rescaled η → U ff ηU bb while keeping η † fixed, which is possible for Grassmann variables. The integrals over U in the third line of Eq. (A.1) are obtained by expressing the determinant det(1 2k b − η † η) as The measure dµ(U ) can be chosen as the normalized Haar measure. Afterwards, we have integrated over η because it is simply a Gaussian. We have got the last line of Eq. (A.1) via the Selberg-type integral [25] U(n) a a,b=1,...,n = n−1 j=0 j! (l + j)! (A.3) with n, l ∈ N. We recall that the Vandermonde determinant is given as The integrals over U bb and U ff follow from similar formulas as (A.3), i.e.
Putting everything together we arrive at the second line of Eq. (4.7).
A.2 Derivation of Eqs. (4.12) and (4.13) The next goal of this appendix is to evaluate the group integral involved in the fermionic partition function (3.6). Here we consider the general case of an arbitrary number k f of flavours. We achieve this by using the parametrization with U 1 , U 2 ∈ U(k f ) and ϑ = diag(ϑ 1 , . . . , ϑ k f ) ∈ [0, π] k f . The Haar measure can be calculated with the help of the invariant length element up to normalization. The normalization constant K can be calculated as follows, (A.10) The unitary matrices U 1 and U 2 drop out in the quadratic term which is proportional to µ 2 , see Eq. (3.6). Thus we are left with a Berezin-Karpelevic integral [82,83] U(k f ) with I ν the modified Bessel function of the first kind. The fermionic partition function is then The Vandermonde determinants can be combined via the Schur Pfaffian identity [84] ∆ 2 where the index a labels the rows and b the columns. Here we normalize the Pfaffian as With the help of de Bruijn's integration theorem [85] we already see at this step that the fermionic partition functions with an arbitrary number k f of flavours can be reduced to those with k f = 1, 2. In particular we can read off from Eq. (A.12) the one-flavour result (4.12) and the two-flavour one (4.13).

A.3 Derivation of Eqs. (4.23) and (4.24)
As for the fermionic partition function, we can calculate the group integral in Eq. (4.20) for a general number of bosonic flavours k b and then read off the result for k b = 1, 2. The integration over U can be done after choosing a parametrization very similar to Eq. (A.7), namely The division with respect to [U(1)] k b is due to the invariance under G → ΦG with Φ an arbitrary unitary, diagonal matrix. One can readily check that the matrix H is indeed Hermitian, positive definite and Hτ 2 is self-inverse. The invariant length element becomes Neglecting the normalization for the moment, the partition function is The normalization of the remaining non-compact coset integral is fixed in the limit L bb κ = κ 0 1 k b → ∞ at µ = 0. Before coming to this problem we want to calculate the non-compact coset integral. The integral over G is the non-compact counterpart of the Berezin-Karpelevic integral (A.11). Indeed one can show that it satisfies the same differential equation as the compact one, i.e.
regardless of whether U 1 and U 2 are unitary matrices or G = U 1 = U † 2 is a complex matrix. The two matrices V and W and independent and k b × k b dimensional and ∆ V = a,b ∂ V ab ∂ V * ba is the Laplace operator with respect to V . After performing a singular value decomposition of V with λ = diag(λ 1 , . . . , λ k b ) its singular values, the singular value part of the Laplace operator becomes such that the differential equation factorizes. There are two fundamental solutions of the eigenvalue equations of the operator ∂ 2 λ j +λ −1 j ∂ λ j , either the modified Bessel function of the first kind I 0 or the modified Bessel function of the second kind K 0 . Since Eq. (A.17) decays exponentially for large κ, we can exclude the first case and we have for the non-compact coset integral We plug this result into Eq. (A.17) and arrive at √ κ 0 + δH 2 /(2κ 0 ) with δH = δH † = −τ 2 δHτ 2 . Therefore we can write δH = δH 1 τ 1 + δH 3 τ 3 with two independent Hermitian k b × k b matrices δH 1 and δH 3 . Here we used again the abbreviated tensor notation H j τ j ≡ H j ⊗ τ j . The metric is then and, hence, the measure becomes with dδH the flat Lebesgue measure. Therefore we have the asymptotic behaviour of the partition function This result has to be compared with the asymptotics of the intermediate result (A.21). Suppose C is the constant we are looking for. Then Now we can identify C by setting this result equal to Eq. (A.24). The final result of this subsection is We start our calculation with the second line of Eq. (4.27). Before we take the limit N → ∞ we integrate over the Grassmann variables. For this purpose we introduce a unitary matrix S ∈ U(2) as ) dS exp N Tr S − N Tr SU (A.28) The first determinant can be expanded as follows.