Individual eigenvalue distributions of crossover chiral random matrices and low-energy constants of SU(2)$\times$U(1) lattice gauge theory

We compute individual distributions of low-lying eigenvalues of a chiral random matrix ensemble interpolating symplectic and unitary symmetry classes by the Nystr\"om-type method of evaluating the Fredholm Pfaffian and resolvents of the quaternion kernel. The one-parameter family of these distributions are shown to fit excellently the Dirac spectra of SU(2) lattice gauge theory with a constant U(1) background or dynamically fluctuating U(1) gauge field, which weakly breaks the pseudo-reality of the unperturbed SU(2) Dirac operator. Observed linear dependence of the crossover parameter with the strength of U(1) perturbations leads to precise determination of the pseudo-scalar decay constant, as well as the chiral condensate in the effective chiral Lagrangian of AI class.

Individual eigenvalue distributions of crossover chiral random matrices and low-energy constants of SU (2) We compute individual distributions of low-lying eigenvalues of a chiral random matrix ensemble interpolating symplectic and unitary symmetry classes by the Nyström-type method of evaluating the Fredholm Pfaffian and resolvents of the quaternion kernel. The one-parameter family of these distributions is shown to fit excellently the Dirac spectra of SU(2) lattice gauge theory with a constant U(1) background or dynamically fluctuating U(1) gauge field, which weakly breaks the pseudoreality of the unperturbed SU(2) Dirac operator. The observed linear dependence of the crossover parameter with the strength of the U(1) perturbations leads to precise determination of the pseudoscalar decay constant, as well as the chiral condensate in the effective chiral Lagrangian of the AI class.

Introduction
The grounds for universality, i.e., insensitivity to details of the system of concern, of local correlation of energy levels of stochastic and quantum-chaotic Hamiltonians have been well uncovered by now, in terms of ten-fold classification of symmetric superspaces on which spectral nonlinear σ models describing spontaneous symmetry breaking reside [1], and of semiclassical equivalence between periodic orbits and the aforementioned σ models [2]. On the other hand, the presence of explicit symmetry breaking perturbations is known to induce a crossover between different universality classes [3], in such a way that is also insensitive to the systems' details. An example of this universality crossover is the Gaussian orthogonal ensemble (GOE)-Gaussian unitary ensemble (GUE) transition in a disordered or chaotic system under a magnetic field [4,5]. In the realm of lattice gauge theory, where Dirac operators play the rôle of stochastic Hamiltonians [6,7], the crossover between the chiral Gaussian unitary ensemble (chGUE) itself and the chGUE-GUE crossover, associated with an imaginary isospin chemical potential [8] and a finite lattice-spacing effect in the Wilson Dirac operator [9], respectively, have been utilized to determine the pion decay constant and the Wilsonian chiral perturbation constants from relatively small lattices. The aim of this work is to apply this spectral approach to the determination of low-energy constants in another setting, namely SU(2) gauge theory under U(1) perturbations, either in the form of a constant imaginary chemical potential [10,11] or a dynamically fluctuating one. In contrast to these preceding works which used n-level correlation functions or the smallest eigenvalue distribution, our strategy in this paper is to employ multiple spectral observables which allows for precise fitting of lattice data, namely, individual distributions of the kth smallest Dirac eigenvalues [12] inclusively. The practical advantages of our method will be proved in the precision of the low-energy constants determined, as preliminarily reported in Ref. [13]. This paper is composed mainly of two parts: Sect. 2 is devoted to analytic treatments in random matrix theory, and Sect. 3 to its application to Dirac eigenvalue distributions measured in lattice simulations. In Sect. 2 we start by reviewing the established results of the chiral Gaussian symplectic ensemble (chGSE)-chGUE crossover, namely, the derivation of the quaternion kernel K [14][15][16]. Then we apply the Nyström-type method [17,18] to that kernel and compute individual eigenvalue distributions in the form of the Fredholm Pfaffian and resolvents ofK [19]. The relationship between the chiral Lagrangian in the ε regime and the nonlinear σ model from random matrices interpolating chGSE-chGUE leads to identification between parameters in each theory. In Sect. 3 we introduce our models of lattice gauge theory and explain our strategy of fitting the Dirac spectra using individual eigenvalue distributions of chGSE and of the chGSE-chGUE crossover. Optimally fitting parameters (mean level spacings ∆ and crossover parameters ρ) will be exhibited in Tables 1-6, leading to very precise determination of the low-energy constants (chiral condensate Σ and pseudo-scalar decay constant F ) as presented in Table 7. Our conclusions, including a possible direction of study, will be summarized in Sect. 4.

chGSE-chGUE crossover 2.1. Crossover random matrix ensemble
Let N and N be even positive integers. Let A and B be (N/2) × (N /2) quaternion matrices, which can be represented as ordinary N × N matrices A and B as * Here the four units of the quaternion field H are represented by the 2 × 2 unit matrix and Pauli matrices, {σ µ } = (I, −iσ 1 , −iσ 2 , −iσ 3 ). Let these matrix elements belong to i.e., A is quaternion-real and B is not (i.e., B is a generic N × N complex matrix). We consider A jk to be independent random variables, distributed according to the Gaussian distributions e − 1 2 tr AA † and e −tr BB † , respectively. We define an ensemble of (N + N ) × (N + N ) Hermitian matrices H of the form where τ is a real parameter, initially introduced by Dyson as a fictitious time for the Brownian motion of eigenvalues [3]. This ensemble is called a "chiral" random matrix ensemble because it enjoys the chiral symmetry {H, γ 5 } = 0 with γ 5 = diag(I N , −I N ). This anticommutation relation implies that the eigenvalues of H consist of min(N , N ) pairs of generically nonzero eigenvalues of equal magnitude and opposite signs (i.e., (±1)× singular values of C), and ν = |N − N | zero eigenvalues. The presence of B violates the quaternion-reality of C and the self-duality of H (i.e., H ji ⊗ σ † µ , i, j = 1, . . . , N + N ), and lifts the Kramers degeneracy of all nonzero eigenvalues of H (i.e., singular values of A alone). Accordingly, this random matrix ensemble interpolates between two limiting cases, the chiral Gaussian Symplectic Ensemble (chGSE) at τ = 0 and the chiral Gaussian Unitary Ensemble (chGUE) at τ → ∞, depending on a single parameter τ (at a fixed, finite N and N ).

Joint eigenvalue distribution
In order to make the paper self-contained, below we sketch the derivation of the probability distribution of singular values of C, and refer the reader to Refs. [14][15][16] for rigorous proofs of the relevant formulas. We start from the unnormalized probability measure of the matrix elements of C, where dA = j,k,µ dA jk . Without loss of generality we assume N ≤ N . We employ the singular value decomposition A = SM S † and C = U ΛU † with S ∈ USp(N ), S ∈ USp(N ), U ∈ U(N ), U ∈ U(N ) and parametrize the singular values as M = diag (µ 1 , . . . , µ N ) 0 N ×ν , Λ = diag (λ 1 , . . . , λ N ) 0 N ×ν . Kramersdegenerate pairs of singular values of A are ordered such that µ i+N/2 = µ i (i = 1, . . . , N/2). The measures dA and dB on quaternion-real matrices and complex matrices take the following respective forms (in what follows we suppress all constant factors in the measures): Here d(S, S ) and d(U, U ) denote the invariant measures on the respective angular degrees of freedom, and denote the Vandermonde determinants N/2 (µ 2 ) := . The probability measure of the singular values {λ i } of C follows from Eqs. (4) and (5) by integrating out the unitary matrices (U, U ). The integrations over symplectic matrices (S, S ) decouple after redefining the unitary matrices by S † U → U and S † U → U , leading to the expression 3/28 We employ the Berezin-Karpelevich formula [20][21][22] for the integration over (U, U ) and take the pairwise confluent limit µ i+N/2 → µ i for all i = 1, . . . , N/2: Here I ν denotes the Bessel function of the pure imaginary argument, I ν (z) = J ν (iz). By substituting Eq. (7) into Eq. (6) and performing a change of variables from the singular values of rectangular matrices C and A to the eigenvalues of Wishart matrices CC † and AA † , x i = λ 2 i and y i = µ 2 i , the probability measure becomes proportional to Here we have introduced the Laguerre weight w(x) := x ν e −x , and the symmetric function With the multiplicative constant chosen as the above, the function g(x, y) admits an alternative interpretation as a one-particle Green's function for the Brownian motion at time τ , with the norm given by h k = (k + ν)!/k! and the "one-particle energy" by γ k = 2k + ν + 1. Using a lemma by Mehta (A.17 of Ref. [23]), the (N/2)-fold integral of an N × N determinant in Eq. (8) can be decomposed into an N × N Pfaffian of single integrals:

Quaternion determinant
In this subsection we summarize the procedure presented in Ref. [23], Sect. 14. We introduce a set of arbitrary monic polynomials {R k (x)} k=0,1,... and arbitrary positive numbers {r k } k=0,1,... , and define functions {ψ k (x)} by and their F -convolutions {φ k (x)} by 4/28 Next we introduce functions D(x, x ), S(x, x ), and I(x, x ), which are bilinear combinations of {ψ k (x)} and {φ k (x)}: and the corresponding N × N matrices D N , S N , and I N by Since a 2N × 2N antisymmetric matrix D N S N t −S N −I N is a product of two rectangular matrices of size 2N × N and N × 2N : On the other hand, Using Eqs. (19) and (20), the probability measure (11) now reads Since its upper diagonal block is the transpose of the lower diagonal block and both off-diagonal blocks are antisymmetric, K N can be regarded a 2N × 2N complex matrix representative of an N × N quaternion self-dual matrix, which we call K N = [K (x i , x j )] i,j=1,...,N . Using Dyson's lemma qdet Φ = Pf (ZΦ) for the quaternion determinant (qdet) [24] of a quaternion self-dual matrix Φ, we finally obtain

Skew-orthogonal polynomials
We would like to choose {R k (x)} and {r k } (which are so far arbitrary) in such a way that the 2 × 2 matrix K(x, y) = S(x, y) J(x, y) D(x, y) S(y, x) representing the quaternion kernel K (x, y) and is correctly normalized: These two relationships would yield a crucial property that the restricted quaternion matrix K n = [K (x i , x j )] i,j=1,...,n satisfies the recursion relation ∞ 0 dx n qdet K n (x 1 , . . . , x n ) = (N − n + 1) qdet K n−1 (x 1 , . . . , x n−1 ), from which a k-level correlation function is expressed as qdet K k = Pf(ZK k ). It is well established that the properties (23), (24) are fulfilled by requiring {R k (x)} to be skeworthogonal with respect to the skew inner product , and {r k } to be their skew-norms, Under this choice, F (x, x ) itself is expressed as Thus the matrix element of due to the definition (17).

Microscopic crossover scaling limit
Now we concentrate on the case in which the Kramers degeneracy is weakly broken by the small parameter τ 1. Then the spectral density ρ(λ) of H in the large-N limit is identical to that of chGSE (τ = 0), i.e., Wigner's semicircle ρ(λ) = π −1 √ 4N − λ 2 . We magnify the vicinity of the origin by introducing the rescaled variables † x i := λ i /∆ which measure the eigenvalues in units of the mean level spacing at the origin, ∆ = 1/ρ(0) = π/ √ 4N . Moreover, in order to realize a nontrivial crossover behavior we take the triple-scaling limit N → ∞, λ → 0, τ → 0 while keeping the combination ρ = √ τ /∆ and x i fixed finite. In this limit, sums over k turn to integrals over v := k/N and Laguerre polynomials reduce to Bessel functions, w(z)L ν k (z) ∼ k ν/2 J ν (2 √ kz) as k → ∞. Accordingly, the quaternion kernel elements (15), (16), (29) reduce to (see footnote † ) [14] S Thus the correlation function of n positive rescaled eigenvalues {x i } of H in the vicinity of the origin is finally expressed as (36) † Here we have abused the notation slightly: Up to Sect. 2.4, x i = λ 2 i denotes squared eigenvalues of H, whereas after Sect. 2.5 x i = λ i /∆ denote microscopically rescaled eigenvalues of H. Accordingly the function symbols K(x, y), S(x, y), etc. are also used with two different meanings. 7/28

Individual eigenvalue distributions
It is well known that, for a determinant process in which n-point correlation functions are expressed in terms of a scalar kernel K(x, y) as R n (x 1 , . . . , x n ) = det [K(x i , x j )] n i,j=1 , the probability E (I) for an interval I to contain exactly points is expressed as a Fredholm determinant (Det) of K over I [23]: . (37) . This argument directly carries over to our case of the quaternion determinant process in which n-point correlation functions are expressed in terms of a quaternion kernel K (x, y) as Tr log(I − ξK I ) in ξ, the first few E (s) are expressed in terms of the Fredholm determinant and the resolvents of the operatorK I , After specializing to I = [0, s) and abbreviating E (s) := E ([0, s)), the probability distribution p k (s) of the kth smallest positive eigenvalue is given in terms of E 0 (s), . . . , E k−1 (s) as This relationship follows from a simple observation that, for a joint of two intervals [0, s + ds) = [0, s) ∪ [s, s + ds) := I ∪ dI, the probability that the narrower interval dI contains 8/28 more than one eigenvalue is of order O(ds 2 ), so to the order O(ds 1 ) one has Subtracting Eq. (41) from the definition of E (s) gives which is equivalent, in the limit ds 0, to with p 0 (s) = 0 understood. Summing over = 0, . . . , k − 1 gives Eq. (40). An efficient way of numerically evaluating the Fredholm determinant of a trace-class operatorK I is the Nyström-type discretization [17,18] is an M × M matrix evaluated with a quadrature rule consisting of a set of M points {x i } ∈ I and associated weights As the order M of the quadrature increases, the RHS of Eq. (44) is proven to converge to its LHS uniformly and exponentially fast in M [17,18]. We also need to evaluate the resolvents in Eq. (39), which are likewise approximated as Obviously these formulae hold for a 2 × 2-matrix-valued kernel as well.  x R1 x Fig. 1 First four eigenvalue distributions p 1 (s), . . . , p 4 (s) (left) for 0.04 ≤ ρ ≤ 0.70 (step 0.01, purple to red) and the microscopic spectral density R 1 (x) (right) for 0.01 ≤ ρ ≤ 1.00 (step 0.01, purple to red) for the chGSE (black) to chGUE (gray) crossover, at ν = 0. 9/28 For comparison, the spectral densities that comprise the former, R 1 (x) = ∞ k=1 p k (x) = S(x, x) in Eq. (33), are plotted in Fig. 1 (right). The practical advantage of adopting individual eigenvalue distributions over n-level correlation functions (including R 1 (x)) for fitting is clear from the figures: As the oscillation of the latter consists of overlapping multiple peaks, the characteristic shape of each peak is inevitably smeared, resulting in a rather structureless curve for which an accurate fit is difficult. On the other hand, the shape of the former is clearly distinguishable and is extremely sensitive to the ρ parameter, because the ratio of the two p k (s) of the chGUE-chGSE crossover at different ρ grows as exp(const.s 2 ) for large s. Therefore, the p k (s) are expected to admit very precise one-parameter fitting of data by the least-squares method (as will be shown in Tables 3-6 of Sect. 3).

Effective theory and low-energy constants
Now we shall relate the crossover random matrix ensemble to the nonlinear σ models originating from gauge theory. In continuum, QCD-like theories with N F flavors of quarks in a real representation have pseudoreal Dirac operators, as does our case of SU(2) lattice gauge theory with fundamental staggered fermions (Eq. (52) in Sect. 3). The low-energy effective Lagrangian of these theories is universally determined by the spontaneous breaking of the Pauli-Gürsey extended flavor group SU(2N F ) down to its vector subgroup SO(N F ) and takes the form in the leading order of the p-expansion. Here Q(x) is a symmetric SU(2N F ) matrix-valued Nambu-Goldstone field (called a nonlinear σ model of class AI), m the degenerate quark mass, andM = σ 1 ⊗ I NF . It contains two phenomenological constants: F the pseudo-scalar decay constant and Σ = ψ ψ /N F the chiral condensate (both measured in the chiral and zero-chemical potential limit m, µ → 0). The effect of introducing the quark number chemical potential µ to the fundamental theory is unambiguously incorporated in L eff through flavor covariantization in the symmetric rank-2 tensor representation [27], If the theory is in a finite volume V = L 4 and the Thouless energy defined as E c ∼ F 2 /(ΣL 2 ) is much larger than m (called the ε-regime), the path integral is dominated by the zero mode only and takes the form The above form (not the concrete symmetric space on which Q takes values) of action is in fact common for all classes of nonlinear σ models (in even or odd dimensions, with Dyson indices β = 1, 2, 4) in which the global symmetry is broken by the chemical potential [28].
On the other hand, the characteristic polynomial det(λ − H) NF of a random matrix H (3) in the microscopic crossover scaling limit explained in Sect. 2.5 can be evaluated by exponentiating the determinant with N F flavors of Grassmannian vectors and using a standard technique of Hubbard-Stratonovich transformations (see, e.g., Ref. [29]) for matrices A and B. One can easily show that the outcome is the same nonlinear σ model as Eq. (48) with 10/28 parameters replaced by respectively. Note that the above identification can also be read off from the exponents in the quaternion kernel elements (33)-(35).

Dirac spectrum
In this section we shall fit the probability distributions analytically derived in the previous section to the Dirac operator spectra measured from two types of lattice gauge simulations: (a) SU(2) gauge theory with the imaginary chemical potential, and (b) SU(2) × U(1) gauge theory. In either case the pseudoreality of the staggered SU(2) Dirac operator is weakly violated by the U(1) field, which is applied as a fixed background [30] or dynamically fluctuating.

Simulation details
The phase 2πϕ can be regarded as the Aharonov-Bohm (AB) flux [4] penetrating the temporal circle and is gauge-equivalent to the imaginary chemical potential µ = iµ I = 2πiϕ/L, i.e., a fixed U(1) background B ν = (2πϕ/L)δ ν,4 . We chose antiperiodic/periodic boundary conditions in the temporal/spatial directions, respectively, and consider a small twisting along the temporal direction (ϕ 1).
(b) SU(2) × U(1) gauge theory: Following Ref. [31], non-compact U(1) link variables B ν (x) are generated under the Coulomb gauge-fixing condition (with an additional constraint for the B 4 (x)) and at the unit coupling constant, and are multiplied to SU(2) link variables with e denoting the bare U(1) coupling constant. Fermions are quenched both for the SU(2) and U(1) gauge fields. For the pure SU(2) case e = 0 (or ϕ = 0), the staggered SU(2) Dirac operator in the fundamental representation is pseudoreal, i.e., satisfies (C denotes complex conjugation) for either choice of periodic or antiperiodic boundary condition in each direction, and we impose periodic boundary conditions in all four directions and consider the U(1) part as a small perturbation (i.e., e 1).

11/28
As the presence of the AB flux ϕ or the U(1) coupling e breaks the pseudoreality (53), they parametrize antiunitary symmetry breaking in Dirac operators and are anticipated to be the lattice gauge theory counterpart of the crossover parameter ρ in random matrix ensemble interpolating chGSE and chGUE. We note that while the effect of the ICP µ I on the low-energy effective Lagrangian is completely dictated on the symmetry ground [32,33] and is related to the pseudo-scalar decay constant F as for the real chemical potential [27], the effect of e cannot be directly related to F , because integrating over the dynamical U(1) gauge field in the effective Lagrangian would lead to nonlocal self-couplings of pseudo-scalar mesons.  Then we define χ 2 from the measured frequency F b = #{λ 2i ∈ I b } and its analytic pre- The statistical error δ∆ i is estimated as a deviation from the optimal ∆ i at which χ 2 /d.o.f. increases by unity. The combined value of the mean level spacing at the origin∆ is obtained as the weighted average of (∆ i , δ∆ i ), i = 1, . . . , 4. We have confirmed that these four data are always mutually consistent, so that their combination helps to improve the statistical error in ∆ as compared to the previous method of using the smallest eigenvalue only [10,11]. 12/28 (ii) Determination of the crossover parameter ρ. Next we switch on the U(1) perturbation (a) or (b) and measure the Dirac spectra {λ k } for N conf = O(10 4 ) independent configurations. The effect of such perturbations on ∆ (i.e., the chiral condensate) is negligible in the lowest order in the ε-expansion that we are working on. U(1) perturbations split once-Kramers-degenerate pairs of eigenvalues, λ 2i−1 < λ 2i . We take the first two pairs of Dirac eigenvalues (λ 2i−1 , λ 2i ), i = 1, 2, and define the unfolded eigenvalues as (s 2i−1 , s 2i ) = (λ 2i−1 /∆ i , λ 2i /∆ i ). Then by using the same strategy as in Step (i), their histograms P k (s k ) are fitted to the analytic results p k (s) (39), (40) of the crossover random matrices, with the crossover parameter ρ k being varied. The statistical error δρ k of the crossover parameter is again estimated as a deviation from the optimal ρ k at which χ 2 /d.o.f. increases by unity. The crossover parameterρ corresponding to a particular choice of ϕ or e is eventually determined as the weighted average of (ρ k , δρ k ), k = 1, . . . , 4. The advantage of using both once-degenerate eigenvalue pairs (over using, e.g., three low-lying eigenvalues) is now evident: Histograms of s 2i−1 and s 2i always shift in the opposite directions under the Kramers-breaking perturbation, so that the effect of small errors in the determination of the overall scale ∆ i in Step (i) (which shifts both histograms in the same direction) is expected to be canceled in the final value ofρ obtained by combining ρ 2i−1 and ρ 2i .

Low-energy constants.
Due to the correspondence (49), two low-energy constants contained in the chiral Lagrangian, the chiral condensate Σ and the pseudo-scalar decay constant F , are directly related to ∆ and ρ measured in Steps (i) and (ii), the former by the Banks-Casher relation Σ = π ∆V ; (54) the latter for (a) the SU(2)+ICP model by and for (b) the SU(2)×U(1) model by Note that the combination "F 2 µ 2 I " in the LHS of Eq. (56) is to be regarded as a single coefficient of the −trBQ †B Q term in the chiral Lagrangian (48). As this quantity should be proportional to e 2 , our aim here is to determine the unknown proportionality constant set by the dynamics. Thus we shall check, for each β in either case of (a) or (b), the stability of the ratio R(ϕ) = ρ/ϕ or R(e) = ρ/e as ϕ or e is varied, and determine its mean value as a weighted average of {R(ϕ), δR(ϕ)} or {R(e), δR(e)}. Tables 1 and 2 we exhibit optimal values of the mean level spacings ∆ i of SU(2) Dirac spectra (a) under the antiperiodic boundary condition on the temporal direction (to be used for SU(2)+ICP), and (b) under the periodic boundary conditions on all directions (to be used for SU(2) × U(1)), respectively. At each β, we adopt 13/28 only the ∆ i that pass the χ 2 test of the fitting: χ 2 /d.o.f. < 2.0. In Fig. 2 (top) we exhibit sample plots of histograms of the four smallest nondegenerate Dirac eigenvalues P i (λ i ) versus the corresponding individual eigenvalue distributions p i (s) of chGSE, each being optimally rescaled by ∆ i . The mutual consistency of such ∆ i observed in Fig. 2 (bottom) justifies the use of the weighted average, listed in the seventh columns in Tables 1 and 2. This sample figure illustrates that the error in ∆ 1 (i.e., in Σ), which is relatively larger than that in the other three ∆ i , is considerably improved by a factor of 5-10 and is down to O(10 −4 ) by the use of the weighted average of the four. Note that histograms at β = 1.75 on V = 4 4 and at β = 2.1 on V = 6 4 (marked by * in Tables 1 and 2) failed to be fitted into the chGSE predictions in the above criterion as χ 2 /d.o.f. exceeds 3. Thus we chose to relax it by cutting off the tails of the distributions for which p 1 (s) < 0.2 from fitting, which leads safely to 14/28 Table 1 Mean level spacings of the pure SU(2) Dirac spectrum on V = 4 4 and 6 4 with the antiperiodic boundary condition on the temporal direction, in units of 10 −2 a −1 .  Table 2 Mean level spacings of the pure SU(2) Dirac spectrum on V = 4 4 and 6 4 with periodic boundary conditions on all four directions, in units of 10 −2 a −1 . Next we exhibit the optimal values of the crossover parameters ρ k (k = 1, . . . , 4), determined by fitting individual Dirac eigenvalue histograms to the chGSE-chGUE predictions, (a) for the SU(2)+ICP model in Tables 3 and 4 Tables  5 and 6. Again we adopted only the cases that passed the χ 2 test with χ 2 /d.o.f. < 2.0, and the exceptional treatment of cutting off the tails of the smallest eigenvalue distributions for which p 1 (s) < 0.2 from fitting is applied to β = 1.75 on V = 4 4 and β = 2.1 on V = 6 4 (marked by * in Tables 3-6). Whenever the histogram of the unfolded eigenvalues s 2i−1 = s 2i of the pure SU(2) Dirac operator is well fitted into the chGSE prediction in our criterion, so is the corresponding pair {P 2i−1 (s 2i−1 ), P 2i (s 2i )} of the perturbed Dirac operator to the chGSE-chGUE crossover prediction, as expected. Figure 3 shows samples of the histograms {P 1 (s 1 ), P 2 (s 2 )} of the perturbed Dirac operator eigenvalues and best-fit distributions of the chGSE-chGUE crossover. As the perturbation ϕ or e is increased, the eigenvalue distributions are clearly seen to follow the random matrix curves. They respond more rapidly to the perturbation for smaller β (left panels) than for larger β (right panels), indicating that F is a decreasing function of β. For the sake of graphical visibility, we also exhibit samples of histograms of {P 1 (s 1 ), . . . , P 4 (s 4 )} and best-fit distributions of the chGSE-chGUE crossover, at fixed β = 1.75 and with increasing perturbations ϕ (Fig. 4) and e (Fig. 5), respectively. Figure 6 illustrates the crossover parameters ρ k determined from P k (s k ) for the SU(2)+ICP model and for the SU(2)×U(1) model at β = 1.75 and V = 6 4 , at each value of perturbation ϕ or e. We notice that ρ 1 and ρ 2 (ρ 3 and ρ 4 ) have a tendency to counter-move across the weighted average of the four, as anticipated in Sect. 3.1. Combined use of P 2i−1 (s 2i−1 ) and P 2i (s 2i ) is indeed seen to reduce the errors of the best-fit parametersρ in both panels. 16 (2) (2) (2) Tables 1 and 2. In Table 7 we list these values in the third column (SU(2)+ICP) and in the fifth column (SU(2)×U(1)), at each β and on V = 4 4 and 6 4 lattices. In addition, the values of the chiral condensate in the thermodynamic limit V = ∞, extrapolated from V = 4 4 and 6 4 , are listed for the coupling range 0 ≤ β ≤ 1.75. ‡ The pseudo-scalar decay constant F in Eqs. (55) and (56) is essentially the proportionality constant between the perturbation strength µ I or e and the crossover parameter ρ. In Fig. 7 we exhibit the ratiosρ/µ I for the SU(2)+ICP model (left) and theρ/e for SU(2)×U(1) model (right), both on V = 6 4 and at various β. For each model we observe excellent stability of the ratios as the perturbation µ I or e is varied. This justifies the use of combiningρ/µ I or ρ/e for all available µ I or e at each β, so that errors in the combined ratios are significantly reduced. Pseudo-scalar decay constants determined in this way are exhibited in the fourth column (SU(2)+ICP) and in the sixth column (SU(2)×U(1)) of Table 7. For the SU(2)+ICP model, Σ and F 2 are determined with 10 −4 and 10 −3 precision, respectively, on our small lattices of 4 4 and 6 4 . This observation, first advocated in Ref. [8] in the context of SU(3) gauge theory with isospin ICP, enables us to extrapolate their value to the thermodynamic limit V → ∞, listed in the bottom rows of Table 7. On the other hand, for the SU(2)×U(1) ‡ Although the fitting range of P k (s) at β = 1.75 on V = 4 4 is compromised as compared to β ≤ 1.5, we dare to estimate the thermodynamic limit of the chiral condensate and pseudo-scalar decay constant using these data. 23/28 model the combination F 2 µ 2 I /e 2 is found to scale linearly with the lattice volume V . Thus we listed the extrapolated values of F 2 µ 2 I /(e 2 V ) in the thermodynamic limit. Finally, in Fig. 8 we exhibit SU(2)-coupling dependences of the low-energy constants on V = 4 4 , 6 4 , and ∞ in the SU(2)+ICP model (top) and in the SU(2)×U(1) model. The error bars are so small that they are almost obscured by the symbols. Table 7 Low-energy constants Σ and F derived from the SU(2)+ICP model and SU (2)

Conclusions
We have analytically evaluated the kth smallest eigenvalue distributions p k (s) for a random matrix ensemble interpolating chGSE and chGUE using a Nyström-type method applied to the Fredholm Pfaffian and resolvents of the quaternion kernel. These random matrix results are applied to fit the spectra of fundamental, staggered Dirac operators of SU(2) gauge theory with imaginary chemical potential and of SU(2)×U(1) gauge theory on small lattices, from the strong-coupling to the near-scaling regions. Combined use of the first four nondegenerate Dirac eigenvalue distributions in place of the spectral density or the smallest eigenvalue distribution of unperturbed SU(2) gauge theory enables us to determine the chiral condensate with O(10 −4 ) precision. Excellent one-parameter fitting of χ 2 /d.o.f. < 2 between non-cumulative individual distributions and eigenvalue histograms is achieved for almost all cases of U(1) perturbations. Combined use of the first four eigenvalue distributions also contributed to a reduction of the errors in the crossover parameter ρ. The acute sensitivity of p k (s) on ρ, and the observed linear dependence of ρ on the perturbation strength (AB flux ϕ or U(1) coupling e) resulted in determination of the pseudo-scalar decay constant F (i.e., the coefficient of the pseudoreality-breaking term) with O(10 −3 ) precision.
Our method of determining F in QCD-like theories, which has proved to be feasible on relatively small-sized lattices, is clearly advantageous over the conventional method of using axial current correlators, which inevitably requires a large temporal dimension. A possible application of our method would be towards technicolor candidate gauge theories with fermions in (pseudo)real representations, such as SU(N ) gauge theory with two adjoint flavors [35], which corresponds to the chGSE class if simulated with an overlap Dirac operator [7]. Hyper-precise determination of its "low-energy" constants (i.e., Higgs couplings) from lattice simulations using our ICP method would, upon comparison with knowledge from collider experiments, contribute to single out a credible scenario from such BSM candidates.