Janossy densities for chiral random matrix ensembles and their applications to two-color QCD

We compute individual distributions of low-lying eigenvalues of massive chiral random matrix ensembles by the Nystr\"om-type quadrature method for evaluating the Fredholm determinant and Pfaffian that represent the analytic continuation of the Janossy densities (conditional gap probabilities). A compact formula for individual eigenvalue distributions suited for precise numerical evaluation by the Nystr\"om-type method is obtained in an explicit form, and the $k^{\rm{\small th}}$ smallest eigenvalue distributions are numerically evaluated for chiral unitary and symplectic ensembles in the microscopic limit. As an application of our result, the low-lying Dirac spectra of the SU(2) lattice gauge theory with $N_F=8$ staggered flavors are fitted to the numerical prediction from the chiral symplectic ensemble, leading to a precise determination of the chiral condensateof a two-color QCD-like system in the future.


Introduction
Random matrix theory (RMT) has served as fundamental tool for analysing quantum spectra of classically chaotic systems. Universality of the level statistics of invariant RMTs provides a basis upon which the system-specific information, due e.g. to the presence of short periodic orbits or to the weak localization effect, may be encoded. In the application of RMT to QCD or gauge theories in general, the focus is on the distributions of several smallest eigenvalues of chiral RM ensembles, as they describe the spectral statistics of gauge-covariant Dirac operators in the broken phase of chiral symmetry. This relation is particularly useful with lattice simulations. If a gauge theory is in the chirally broken phase and not in the conformal window, its low-energy excitations are unambiguously described by the chiral Lagrangian on one of the Riemannian symmetric spaces (Nambu-Goldstone manifolds) M. In that case, (i) the low-lying Dirac eigenvalues 0 ≤ λ 1 ≤ λ 2 ≤ · · · measured on lattices of different volumes V will, after prescribed unfolding x k = ΣV λ k and scaling of quark masses µ f = ΣV m f , with a constant Σ independent of the volumes, obey a single statistical distribution p k (x; {µ f }) = δ(x − x k ) , and (ii) this distribution will be identical to the one from the RMT that is equivalent to the zero-mode part of the chiral Lagrangian on M [1]. If the theory is in the symmetric phase of the chiral symmetry, no such scaling with the volume, which collapses the distributions of λ k 's from different volumes onto a single function, would appear. Previously this criterion was applied to QCD around the critical temperature, and the inconsistency with RMT (including non-scaling of unfolded Dirac eigenvalues with volumes) was considered as a sign of chiral symmetry restoration [2]. In addition, if the theory is conformal, no scale should appear so that the chiral condensate Σ should disappear in the chiral limit and description with RMT is not applicable.
In the proposal of the walking technicolor model [3], the choice of the gauge group of techni-gluons and the representation of techni-quarks are rather open (as long as the one-loop beta function coefficient is negative and small), since these particles would be confined under the energy scale of several hundred TeV and would escape direct detection. This spurred extensive numerical searches of the conformal window (where β(g * ) = 0) and the walking regime (where β(g) < 0 but small) on various lattice settings with choices of colors/flavors/representations. Summaries of recent activities with lattice simulations are found in [4][5][6]. In an attempt to identify the chirally broken phase below the conformal window for the SU(3) N F = 4 and 8 systems, Fodor et al. [7] fitted the Dirac spectra of these gauge theories to the analytic results from the chiral GUE (Dyson index β = 2). Subsequently, one of the present author (I.K.) and others tried a similar comparison of the Dirac spectrum of the SU(2) N F = 8 system (see e.g. [8][9][10] for the current situation of this system) to the chiral GSE (β = 4) [11].
For the above approach of fitting Dirac spectra to the corresponding RMT predictions to be practically useful, it is highly desirable to single out individual distributions of each of the ordered RM eigenvalues p k (x) from the spectral density ρ(x; {µ}) = k δ(x − x k ) = k≥1 p k (x; {µ}), as the latter becomes rather structureless after a couple of oscillations ( Fig. 1). The standard technique to access such individual eigenvalues is to use the spectral kernel. Once the spectral kernel is obtained, one can give an analytic expression of the distribution. Moreover, by combining Nyström-type (quadrature) evaluation of Fredholm determinants and Pfaffians, one can numerically evaluate the distribution of individual eigenvalues. Damgaard and one of the authors (S.M.N.) have previously derived analytic expressions of such individual eigenvalue distributions for chiral RM ensembles at three Dyson indices β and with scaled quark mass parameters {µ f }, initially by the shift-ofvariable method [12] and later by the Nyström-type evaluation of Fredholm determinants and Pfaffians of the spectral kernels [13]. There, technical difficulties have prevented us from obtaining analytic formulas for the chiral GSE (β = 4) with even numbers of massless flavors and for the chiral GOE (β = 1) with even values of the topological charge. Especially, the former restriction is frustrating, as it obstructs applications to the SU(2) systems with N F = 8 and 12 staggered flavors that are popular lattice settings of walking technicolor candidates. Because of this reason, the Monte Carlo method with finite-size matrices was used in [11] to generate the spectral distribution of the RM side in their analysis of SU(2) N F = 8 system. The purpose of this paper is to lift this restriction by providing an analytic formula for the conditional gap probability, a.k.a. the Janossy density, that interpolates the ordinary determinantal or Pfaffian formula for the k-point correlation function and the Fredholm determinant/Pfaffian expression for the gap probability. Then our formula is numerically evaluated very efficiently by the Nyström-type method. As an application of our result, the low-lying Dirac spectra of the SU(2) lattice gauge theory with N F = 8 staggered flavors are fitted to the derived RM prediction. 1 This paper is organized as follows. In Sect. 2 we start by reviewing known formulas on the spectral-statistical distributions of chiral RMTs and their Janossy densities. In Sect. 3 we present a formula for the individual eigenvalue distributions suited for precise numerical evaluation by the Nyström-type method. Specifically, we shall provide numerical data of p 1 (x), . . . , p 4 (x) for the chiral GSE with N F = 4 and 8 degenerate massive flavors. In Sect. 4 we determine the values of chiral condensate of the SU(2) system with N F = 8 the first eigenvalue distribution of the corresponding chiral GSE. Conclusions and discussions on feasible applications of our results are presented in Sect. 5. In order to avoid plethora of formulas in RMT and of lattice details in the main text, some of them are relocated to the Appendices.
2 Fredholm determinants and Pfaffians for chiral Gaussian random matrix ensembles In this section, we will summarize some necessary ingredients about the chiral random matrix ensembles, and derive our main formulae for the Fredholm determinants and Pfaffians of Gaussian chiral random matrix ensembles.

Gaussian chiral random matrix ensembles and the microscopic limit
Each ensemble is labelled by the Dyson index β = 1, 2, 4, respectively. The non-negative integer ν denotes the corank of the matrix H = 0 W W † 0 and will correspond to the topological charge when H is interpreted as modelling Dirac operator of a gauge theory [14]. Let Z N,β,ν ({m a }) be the partition function for the Gaussian chiral random matrix ensembles with α parameters m a (a = 1, . . . , α), which will correspond to quark masses, such that (2.2) [7]. Our lattice data is obtained with the unimproved staggered fermion action and suffers from large taste breaking effects. Consequently, we do not observe the 4-fold degeneracy characteristic of the staggered tastes, and the lightest of these corresponds to 2 flavors. Moreover, due to the Kramers degeneracy of the SU(2) Dirac operator, the degeneracy of the lightest fermion modes are 4-fold, to which we must compare the prediction of RMT with NF = 4 instead of NF = 8.
Likewise the p-level correlation function R (p) N,β,ν (λ 1 , . . . , λ p ; {m a }) of the Hermitian matrix H is defined by [15,16] , Here we introduce variables z j 's such that The p-level correlation functions for β = 2 are rewritten as the determinant of the spectral kernel K(z i , z j ) [17][18][19][20]: N (x 1 , . . . , x p−α ; {m a }) is given by the determinant of the scalar kernel [21]. For β = 1, 4, the skew-orthogonal polynomial method involves the quaternionic determinant qdet [22] of the quaternionic kernel [15,16]. In particular, p-level correlation functions are given by (p + α) × (p + α) quaternionic determinants of the quaternionic kernel, which is rewritten by a 2(p + α) × 2(p + α) Pfaffian of its C-number 2 × 2 representative (denoted by the same K(z i , z j ) for notational simplicity), where Z = iσ 2 ⊗ I p+α stands for the skew-unit matrix Z 2 = −I 2(p+α) . Now we will consider the asymptotic limit: This limit corresponds the microscopic limit of the QCD-like theory on a box of volume V such that where Σ stands for the chiral condensate in the chiral limit.
2 In [21], an alternative representation of the p-level correlation function R (p) (ζ1, . . . , ζp; {µa}) is also found for general mass parameters. (See eq. (B.5) in Appendix B.1.) 3 An explicit formula for the p-level correlation function is known as well for the chiral GOE (β = 1) [16], but the convergence of the Nyström-type discretization of the Fredholm Pfaffian is not guaranteed due to the discontinuity of sgn(ζ − ζ ) in its kernel elements. To avoid such analytical difficulty, we will focus on the study of the Fredholm Pfaffian for the chiral GSE, and leave discussions of the chiral GOE for the future work.

Individual eigenvalue distributions
We now focus on the individual distribution of the k th smallest eigenvalue for the chiral random matrix ensembles. There are various techniques to analyze the gap probabilities [23,24] such as linear differential equations [25] or Painlevé transcendental equations á la Tracy-Widom [26]. An alternative method to find individual distribution of the k th smallest eigenvalue in the asymptotic limit (2.9) has also been developed in [12,27]. The procedure of this method consists of three steps [13]: 1. Relate the joint distribution of the first k eigenvalues to the partition function with βk + β(ν + 1)/2 − 1 additional masses and a fixed topological charge 2/β + 1.

Integrate over the scaled variables
On actual implementation of the above method, the numerical integration over k scaled variables in the third step becomes resource-consuming. To circumvent such technical issue, we will consider Fredholm determinants and Pfaffians for the chiral random matrix ensembles with α mass parameters as the generating function of the joint distribution of the first k eigenvalues, and utilize the quadrature method [30] to evaluate them numerically [31][32][33][34]. In this section, we will derive a compact formula 4 of Fredholm determinants and Pfaffians which will be efficient for numerical computations.
Let P N,β,ν (x 1 , . . . , x N ) be the distribution of the probability for all eigenvalues of the rank N matrix, 15) and P N,β,ν (x 1 , . . . , x N ) obeys the normalization condition.
where χ I (x) stands for the characteristic function on I ⊂ R. If I is a line segment [a, b] (a < b) or a semi-infinite line, the characteristic function is given by (2.17) If I consists of one point {y}, Consider the joint probability E(k; I; {m a }) that one finds exactly k eigenvalues on an interval I along the real axis and α eigenvalues in R <0 such that Such a joint probability E(k; I; {m a }) is known as an analytic continuation of the Janossy density [19,37,38]. (See Appendix C for the definition of the Janossy density.) The cumulative distribution F k (s) and the probability distribution p k (s) of the k th smallest positive eigenvalue are expressed by In the next subsections, we shall show that the generating function τ (z; I; {m a }) of the probability E(k; I; {m a }) given by is rewritten as a block-decomposed Fredholm determinant or Pfaffian of the spectral kernels in (2.7) or (2.8).

Fredholm determinant for chiral Gaussian unitary ensemble
We start by sketching the proof for the simplest case β = 2, α = 1, m 2 1 = −y: To rewrite correlation functions σ (k) N,β=2,ν in terms of the spectral kernel (2.25), we will prepare some notations such as (2. 23) In addition, we assume that the quadrature discretization of the Riemann integral on I to be well-defined in the continuum limit M → ∞ (which is always implicit below), (2.24) We further introduce following notations for the block decomposition of the spectral kernel integrated over I.
Adopting eq. (2.7) and these notations, one can rewrite the Fredholm determinant τ (z; I; √ −y) in terms of the block-decomposed scalar kernel as follows: Reorganizing summations, one finds Thus we obtain a compact expression of τ (z; I; √ −y) in terms of the Fredholm determinant.
The generalization to the case with α eigenvalues lying at y a (a = 1, . . . , α) proceeds in the same way as the derivation of eq. (2.26), leading to where the notation for the block decomposition of kernels (2.25) is generalized as (2.28) The numerator K(z) clearly interpolates the (ordinary) determinantal form for the k-level correlation function det κ in the case I → ∅ (for which k, K → 0) and the Fredholm determinantal form det (I − zK) for the generating function of the gap probability in the 'quenched' limit y a → ∞ (for which κ → I and k → 0). Finally, changing the eigenvalue variables back to the chiral Gaussian and taking the asymptotic limit (2.9), eq. (2.27) leads to with the kernel elements given by their scaled forms (2.12), (2.30)

Fredholm Pfaffian for chiral Gaussian symplectic ensemble
Generalization of the result of the previous subsection to the chiral GOE and GSE is straightforward: one finds the quaternionic determinant formula simply by replacing K with the quaternionic kernel and "det" with "qdet" simultaneously, because the quaternionic determinant shares the same linear algebraic properties which are utilized in the derivation of the determinant formula (2.29). In particular for the chiral GSE, one can use the explicit formulae of the correlation functions and spectral kernels for N F = 4α and N F = 2α in [16]. (See Appendix A.) Indeed, applying the correlation functions R (p) in Appendix A to eq. (2.22) and repeating the same steps leading to zeq. (2.29), one finds the following Pfaffian formula For the quadruply degenerated case N F = 4α, z • K with µ a (a = 1, . . . , α) is given by where Z = iσ 2 ⊗ I α+M , and For the doubly degenerated case N F = 2α, z • K for even α with y a = −µ 2 a (a = 1, . . . , α) is given by and z • K for odd α is by Matrix elements of z • K in the asymptotic limit (2. In case that some of the masses µ a 's are degenerated, one should adopt the confluent limit of the spectral kernel. Some details of the confluent limit of the spectral kernel is discussed in Appendix B.

Numerical evaluation of the Janossy density via the Nyström-type discretization
In evaluating the Fredholm determinant (2.27) and Pfaffian (2.31) numerically, the Nyströmtype discretization proves to be a highly efficient method [39,40]. This numerical method is based on the quadrature rule (see a brief summary in Appendix E), and in seminal works by F. Bornemann, it is shown that the Nyström-type discretization of Fredholm determinants of integral operators of trace class (i.e. for unitary and symplectic ensembles) converges exponentially as the order of the discretization grows. In the following, we employ the Gauss-Legendre quadrature rule of order M with the nodes ζ i and the weights w i = dζ i (i = 1, . . . , M ) given in eq. (E.5).

Chiral GUE with doubly degenerated masses N F = 2α
The Nyström-type discretization of the Fredholm determinant for the individual eigenvalue distribution with β = 2 and N F = 2α is given as follows: where the matrix elements K AB are found in eq. (2.12).

Chiral GSE with quadruply degenerated masses N F = 4α
For β = 4 and with N F = 4α, the Nyström-type discretization of the Fredholm Pfaffian is given by In order to verify the numerical computation with the Nyström-type discretization, we closely looked at the difference F 1 (s)−F HMC 1 (s) for matrix ranks N = 250, 1000, 2000. Fig. 3 right shows that the computational results of the hybrid Monte Carlo simulation indeed converge to the Nyström-type discretization as N grows, confirming that the numerical evaluation of F 1 (s) by the Nyström-type discretization at M = 50 is good enough.  • N F = 8

The generalized gap probability
Numerical plots of F k (s) and p k (s) = ∂ s F k (s) for N F = 8 (α = 2) with 8-fold degenerated mass µ = µ 1 = µ 2 , computed at M = 128, are depicted in Figs. 5. The computed numerics of F k (s) are appended as Supplementary Material because this case is practically important within our application to the two-color lattice QCD with staggered quarks; the case with N F = 4 (α = 1) doubtlessly has its chiral symmetry broken as in the ordinary QCD, and those with N f ≥ 12 (α ≥ 3) have negative 1-loop β-function coefficients β 0 = (11N C − 2N F )/3 and are IR free. Accordingly, N F = 8 is the only case which evokes the question of whether its nature is either QCD-like, conformal, or walking (which would nominate the model as a possible candidate for the technicolor scenario), and motivates us to compare its Dirac spectrum to the massive chiral GSE prediction so as to confirm or exclude if it is QCD-like. In the upper panels, mass parameters are chosen at µ = 0 (black, N F = 0 with ν = 4), 0.5, · · · (step 0.5), 10, · · · (step 1), 20, · · · (step 2), 30, · · · (step 5), 60, · · · (step 10), 100, 200 (red to purple), ∞ (gray, N F = 0 with ν = 0). The lower panels are interpolations of the upper ones.

Chiral GSE with doubly degenerated masses N F = 2α
For β = 4 and with N F = 2α, the quaternionic kernel for Janossy density of the β = 4 ensemble is treated independently for even and odd α.
For the case of odd α, the Nyström-type discretization of τ (z; [0, s]; {µ a }) yields  • The confluent limit for N F = 2 + 2 + 2 + 2 The next example is the chiral GSE for N F = 2 + 2 + 2 + 2 (i.e. α = 4) in the complete confluent limit. The Nyström-type discretization of E 0 (s) of order M is given by where the matrix elements S As the N grows, the HMC results converge to the result from Nyström-type discretization. Compared with the previous cases, however, its convergence is slower.
By applying the explicit expressions in Appendix D, we can evaluate E k (s)'s numerically. Plots for F k (s) in 0 ≤ s ≤ 6 are depicted in Fig. 8, and we find a good agreement with the computations of the hybrid Monte Carlo simulation with N = 4000.

Application: chiral condensate from lattice data
As an application of our RMT results, we use the Dirac eigenvalues of the SU(2) gauge theory with N F = 8 quarks in the fundamental representation. A partial analysis of this system has been presented in [11], where the Monte Carlo method is used to generate the RMT data. Full analyses using the current RMT result will appear elsewhere [41]. As stated in the Introduction, we should use the chiral GSE with N F = 4, because due to the taste breaking effect, the 4-fold degeneracy for the staggered fermions is totally broken so the number of lightest flavors is in fact N F = 2. Furthermore, the pseudo-reality of the SU(2) gauge group yields an additional 2-fold degeneracy yields an additional 2-fold degeneracy, by which N F = 2 is promoted to N F = 4.
The microscopic eigenvalue density is related to the Dirac spectrum through where λ i denotes the eigenvalue of the Dirac operator, Σ the chiral condensate, V the 4volume, and m f the quark masses. We relate the smallest Dirac eigenvalue distribution from lattice simulation through 5 The parametersV ,Σ andm f are the dimensionless 4-volume, the chiral condensate and the fermion mass of the SU(2) gauge theory in the lattice unit, respectively. Dimensionful quantities are λ 1 =λ 1 /a, V = a 4V , Σ =Σ/a 3 , and m f =m f /a, where a is the lattice spacing. The distribution of the smallest eigenvalue p latt. 1 (λ 1 ;m f ) is determined from lattice simulation and its normalization is fixed by As the sole undetermined quantity in eq. (4.2) is the chiral condensate, we can use this relation to best-fit the value ofΣ. If the fit does not work, that is, if eq. (4.2) is not numerically satisfied by any choice ofΣ, it implies that the chiral symmetry is restored and the RM description is not applicable. Note that ζ 1 and µ are dimensionless so they are directly related to quantities in the lattice unit: ζ 1 = λ 1 V Σ =λ 1VΣ and µ = m f V Σ = m fVΣ . An integrated version of eq. (4.2) is where s =ŝVΣ. We use I(ŝ) in the fitting process. Our lattice setting is the following. We have three different lattice sizes, (T /a) × (L/a) 3 = 8 × 8 3 , 12 × 12 3 and 16 × 16 3 . In this paper, we use fermion massm f = am f = 0.010. We use several values of the bare gauge coupling β = 4/g 2 , for which we use β = 1.1-1.475. These values are almost the same ones as used in [11]. See Table 2 in Appendix F for the details of the lattice data. The topological charge ν is calculated with the APE smeared [42] configuration with order-a improved (i.e., "clover") field strength. Note that this gluonic definition does not give an integer value on a lattice. The obtained values, however, cluster around integer values so that we can identify configurations with ν = 0. Eigenvalues and topological charges are calculated for every 10 trajectories.
The details of our fitting procedure is as follows: We divide a given lattice eigenvalue distribution into N bin = 25 bins, whose support covers from 0 to 1.3 times the largest value in the distribution. In addition to the average value and the error in each bin, we estimate the correlation matrix C between bins by using the jackknife method. Since a naive estimation of the correlation matrix causes unstable fitting, we use an improved estimation of the inverse, C −1 imp. . See appendix G for the details. The value of the chiral condensateΣ is determined by minimizing the correlated χ squared: where To estimate p RMT 1 (ζ 1 , µ) with arbitrary ζ 1 and µ, which is needed to calculate I RMT (s i ;Σ) for a givenΣ, we use interpolations in both ζ 1 and µ. We first interpolate in µ and then in ζ 1 , with the 4-point interpolation is used for both. Near the boundary of the available points where the 4-point interpolation is not possible, an interpolation with 3 points or an extrapolation with 2 points is used as well. Fig. 9 is a typical example of a good fit (indicating the chirally broken phase) and a bad fit (chirally symmetric phase). In the broken phase, the RMT well describes the smallest eigenvalue distribution from the lattice data, with a reasonably small value of χ squared. On the other hand, in the broken phase, the RMT curve can by no means describe the lattice data. In the figure, the plotted curve is the result with the best value ofΣ = a 3 Σ. The value of χ squared, however, indicates that the quality of the fit is poor in the right panel and the RMT result is rejected as fitting ansatz. µ=3.20 (04) 1-E 0 (s) s RMT lattice Figure 9. Typical example of a good fit (left) and a bad fit (right). The horizontal scale for the RMT curve is determined by the the best value of the chiral condensate, which is denoted in the plot.
It is interesting to note that even though the fit result is unreliable in the symmetric phase, the obtained value of the chiral condensate is small and consistent with zero, as should be in the symmetric phase. This is clearly seen in Fig. 10. We observe that the larger the bare coupling β = 4/g 2 is, the smaller the obtained chiral condensate becomes and eventually the fit becomes unreliable near the vanishing of the chiral condensate at around β = 1.4-1.5. In this Figure, the unreliable data points, for which χ 2 par degrees of freedom exceeds 1, are plotted with pale colored symbols. Such behavior is also reported in [11], where the HMC with N = 400 is used to obtain the RMT result.

Conclusions and discussions
We have numerically evaluated the k th smallest eigenvalue distributions of chiral random matrix ensembles with multiple flavors using the Nyström-type method applied to the Fredholm determinant and Pfaffian describing the Janossy densities. Adopting the compact determinant formulas (2.27) and (2.31) for the Fredholm determinant for the Janossy densities, we performed numerical computations for the chiral GUE and GSE in the asymptotic limit (2.9). One of our goals of these analyses is an application to the two-color QCD with N F fundamental staggered flavors. For the system of N F = 8 flavors in the fundamental representation of SU(2), the distribution of eigenvalues of the Dirac operators is being studied through the lattice simulation [11]. In the simulation we used, the taste symmetry of the staggered fermions is completely broken due to the finite lattice spacing, so that the remaining flavor symmetry is merely N F = 2. In addition to this flavor symmetry, due to the pseudo-reality of the fundamental representation of the SU(2) gauge group, all the eigenvalues of the Dirac operator are doubly degenerated. As a result, the distribution of the Dirac eigenvalues can fit with the chiral GSE with quadruply degenerated masses N F = 4 in the broken phase Σ = 0.
As shown in the left panel of Fig. 9, we observed that the fitting with the chiral GSE works out very nicely in the broken phase. As the bare coupling β = 4/g 2 grows the chiral condensate becomes smaller and eventually the fitting becomes unreliable at around β = 1.4-1.5 (Fig. 10). This implies that the chiral condensate vanishes and the symmetry is restored at β 1.45. A detailed analysis with more lattice data is currently ongoing [41]. We note that even with large values of the scaled quark masses µ f , fitting with the quenched chiral GSE is valid as long as the magnitude of the eigenvalue is much smaller than F 2 π /(Σ √ V ). Although the value of F π , the pion decay constant, is not available from the current lattice data, it is natural to assume that the smallest of the Dirac eigenvalues satisfies this condition in the broken phase. The Banks-Casher relation tells us that the smallest eigenvalue is small enough to give a non-zero eigenvalue density around the origin.
Finally we will list some directions for the future research. Firstly, the numerical computations developed in this article could also be applied to the two-color QCD with N F = 8-12 fundamental flavors. Among such systems, the existence of the conformal window is strongly expected, and the technology of their lattice simulations is developping remarkably in recent years. We anticipate that the RMT analysis of the spectral statistics of the Dirac operators would discriminate the (near-)conformality of the QCD-like systems and unveil some novel aspects of the conformal window.
The Fredholm Pfaffian for the Janossy density of the chiral Gaussian orthogonal ensemble will deserve a future study direction; the chiral GOE describes the distributions of the Dirac eigenvalues for QCD-like systems with staggered fermions in adjoint representation of SU(N c ) [14]. It is known that the exponential convergence of the Nyström-type discretization of the Fredholm Pfaffian for the orthogonal ensemble is not guaranteed due to the infinite oscillations originating from the discontinuity of the quaternion kernel elements. Even though such hard problem resides, we may still be able to apply the Nyström-type discretization for the practical purpose if the error can be estimated appropriately, and use it to estimate the value of the chiral condensate Σ for the adjoint QCD-like system.
In [35,36], an exact analysis of the Janossy density for the unitary ensemble is done on a basis of the Painlevé II transcendent and its associated isomonodromic system. Generalization of such an exact analysis to the symplectic and orthogonal ensembles could be an interesting direction yet to be studied, and it can be compared with our numerical results.
Recent years, the (0 + 1)-dimensional fermionic model with all-to-all random interactions referred to as Sachdev-Ye-Kitaev (SYK) model [43,44] has been studied very actively in the context of the nonequilibrium quantum many-body systems and its application to the AdS/CFT correspondence (see references in a recent review article [45]). The level statistics of the SYK model was numerically examined, and good agreements with the RMT have been observed. It would be interesting to explore how the Fredholm determinant or Pfaffian expression for the Janossy density of the chiral random matrices appears in the level statistics of the supersymmetric SYK Hamiltonian [46,47].

A Quaternionic kernels for the chiral Gaussian symplectic ensemble
In this appendix, the explicit forms of the spectral kernel for the chiral GSE are summarized for quadruply degenerated masses N F = 4α and doubly degenerated masses N F = 2α.

A.1 Spectral kernel for quadruply degenerated masses
The scaled correlation function of the β = 4 chiral RMT with quadruply degenerated masses N F = 4α in the scaling limit (2.9) is given in [16].
where elements of block matrices are

A.2 Spectral kernel for doubly degenerated masses
The p-level correlation function for β = 4 chiral RMT with N F = 2α doubly degenerated masses in the scaling limit (2.9) is given in [16].
For even α, the kernel ZK (p) is given as follows: where S AB 's, D AB 's, and I AB 's are the same as N F = 4α in Appendix A.1, and T stands for the transposition of the block matrix.
For odd α, the kernel ZK (p) is given as follows: where S AB 's, D AB 's, and I AB 's are the same as N F = 4α in Appendix A.1, and Q + (ζ) = 2 ζζ Let Z β=2,ν (x 1 , . . . , x n ) be the partition function which is obtained as the scaling limit (2.9) of the chiral Gaussian unitary ensemble with N f = 2n mass parameters x a = m a /∆ [21]. . (B.1) To consider the confluent limit x i → x j of this partition function [48], we will use the l' Hôpital's rule given as follows.
Let f, g be differentiable functions on an interval I ∈ R. Assume that for c ∈ I, One finds that the confluent limit x i → x 1 = x (i = 1, . . . , n) of the partition function Z β=2,ν (x 1 , . . . , x n ) by adopting the l' Hôpital's rule (B.2) repeatedly.
Lastly, for the odd α case in addition to the above replacements, the matrix elements Q − 's in eq. (A.7) are also replaced by . (B.14)

B.2.1 Quadruply degenerated kernels in the confluent limit
For the chiral GSE with N F = 4α quadruply degenerated masses, we can use the spectral kernel given in Appendix A.1. In the complete confluent limit, some matrix elements in eqs. (A.2)-(A.4) are replaced as follows. (We also choose ν = 0 for simplicity.) C Janossy density C.1 Janossy density for the determinantal random point process Below we shall overview the definition of the Janossy density for the determinantal random point process [49][50][51]. Consider an ensemble of N particles on Z with the joint distribution (see (1) in Fig.11) given by with the kernel K = [K(n, m)] n,m∈Z obeying the projective condition: Then the k-point function R k (n 1 , . . . , n k ) is given by Consider the probability J k,I (n 1 , . . . , n k ) of finding no particle in an interval I ⊂ Z except for k designated point. (See (2) in Fig.11.) J k,I (n 1 , . . . , n k ) is called Janossy density [52], which is given by the restricted kernel K I = [K(n, m)] n,m∈I on I for the determinantal point process. Here we denote the restricted kernel by K I (n, m) = n|K I |m with the orthonormal complete basis {|n | n ∈ I} and its dual { n| | n ∈ I}. The first line of (C.4) is quoted e.g. from [53] (π(X) on page 341), and the second line is by the identity Generalization to the probability J p,k,I (n 1 , . . . , n k ) of finding exactly p particles in I except for k designated points is straightforward (see (3) in Fig.11). Just as in the case of the ordinary gap probability (k = 0), we merely introduce the spectral parameter z so that J p,k,I (n 1 , . . . , n k ) is given by For the continuous determinantal random point process on X ⊂ R with the measure µ, the Janossy density J k,I (x 1 , . . . , x k )µ(dx 1 ) · · · µ(dx k ) for the distribution of the particles in a subset I ⊂ X is defined as the probability density of finding exactly k particles in I and one at each of the k infinitesimal intervals (x i , x i + dx i ) ⊂ I. J k,I (x 1 , . . . , x k ) is given by the Fredholm determinant det(I − K I ) and the determinant of L I := K I (I − K I ) −1 such that The partition function Z N,β,ν ({m a }) of the massive chiral Gaussian ensemble with N F = βn fermions is given by It can further be rewritten as an N + n eigenvalue integral in the following form: This partition function is regarded as that of the determinantal random point process for 2) around z = 1 is found as combinations of the functional traces T n 's as follows. (The same expansions for the quenched (α = 0) ensembles are given in [13].) The functional traces consist of operators K (n) 's are given by Using these operators, one can showT n (k = 1, . . . , 7) as follows.
The quadrature rule is an efficient method to perform the numerical evaluation for the integral of the smooth function. The quadrature formula for the integral over the interval is represented as [30] Q where w i and x i denote the weight and nodes, respectively, determined by the prescription of the quadrature rule. There are several kinds of quadrature rules. The most basic method is the Gauss-Legendre rule and more efficient one is the Clenshaw-Curtis rule. In the following, we will summarize the Gauss-Legendre rule.

F Details of the lattice result
In Table 2, we list the result of the fitting of lattice data.

G Estimation of the correlation matrix
An element of the correlation matrix is given where y i = I(ŝ i ) andŝ i is the upper end of the i-th bin. The bracket · represents the average over lattice configurations which belong to ν = 0 sector. Since the correlation matrix is an average of fluctuation, one needs to use a resampling method like jackknife or bootstrapping to estimate. In this analysis, we use the jackknife method. What we need in the fitting is not the correlation matrix itself but its inverse. As the estimate of C contains some error, we need some care to invert it. If the bin width is too fine, neighboring bins may give (almost) the same value which causes zero mode (or almost zero mode) of the correlation matrix. If eigenvalue of C is too small, the relative error of the eigenvalue becomes large, which makes estimation of C −1 unreliable. Note that the smallest eigenmode gives the largest contribution to the inverse.
We therefore employ the following steps. First of all, some of the bins do not have eigenvalues of the Dirac operator in it (the largest several bins and sometimes the first bin(s)). Let us suppose that i-th bin has no eigenvalue. Then, i-th column/row of the correlation matrix, C ij and C ji for arbitrary j becomes zero as y i is always 1 (or always 0). This obviously reduces the rank of C. We therefore replace the diagonal element C ii = 0 with the upper bound of the estimate, 1/n 3 , where n is number of independent configurations we use 6 . The off-diagonal elements are kept zero. After this modification of the correlation matrix, which is now denoted as C , we still may have very small eigenvalues. Numerically, we even may observe (small) negative eigenvalue of C 7 . We therefore truncate the correlation matrix by cutting small eigenmodes in inverting the matrix to give an improved estimate of the inverse of the correlation matrix C −1 imp. . The cutoff c cut we use is 0.1 times smallest diagonal element, c cut = 0.1/n 3 . That is,

H Hybrid Monte Carlo (HMC) for RMT
A hybrid Monte Carlo simulation technique [54] is applicable to finite N random matrix theory. By introducing ζ i = √ 8N x i and µ a = √ 8N m a as eq. (2.9), the partition function (2.2) becomes where C represents irrelevant normalization factor and the action is The dynamical variables here is the eigenvalue ζ i . The Hamiltonian for the HMC is where p i the conjugate momentum to ζ i . It is straightforward to write down the equation of motions and apply the HMC algorithm. For the molecular dynamical time evolution, we use a leapfrog integrator. The only non-trivial part is ordering of the variables. We assume that 0 < ζ 1 < ζ 2 < · · · < ζ N . Since there is a divergence in the potential at ζ i = 0 and ζ i = ζ j (i = j), if the initial configuration satisfies this ordering, a smooth molecular dynamical evolution keeps the configuration satisfy the same constraint. Discrete time evolutions, however, can break the constraint so that we use the so called retry trick. We check whether the trial configuration satisfies the constraint before the metropolis test. If it does not, rerun the molecular dynamics with the same random momentum but a finer time step, δτ → δτ /2. If the constraint is still broken after several reductions of the time step (our limit is 6 times), the trial configuration is rejected. For β = 4, the frequency of the retry is order 0.01% and we did not encounter rejections for this reason. As β becomes smaller, the effect of the potential barrier becomes weaker. In fact, more frequent retries are needed for β = 2, and some trial configurations are rejected in the end. Note that β = 1 and ν = 0, the potential barrier at ζ i = 0 disappears.
Here is some parameters we used in β = 4 case. The trajectory length between Metropolis test is τ = 1. We keep the acceptance ratio rather high, typically 0.96-0.97, to reduce the frequency of retries. To avoid the auto correlation, we measure ζ i every 500 trajectories. In making the distribution in Figs 3-8, we check the integrated auto correlation, which is 2τ int 1.2 and used every 2 measurements. I Data of k th smallest eigenvalue distributions for chiral GSE with N F = 8 Numerical data of F k (s; µ) (k = 1, 2, 3, 4) for the chiral GSE with N F = 8 degenerate flavors, in the range 0 ≤ s ≤ 20 and 0 ≤ µ ≤ 100 are appended as a Mathematica Notebook "F1234_chGSE_NF8.nb".