New universality classes of the non-Hermitian Dirac operator in QCD-like theories

In non-Hermitian random matrix theory there are three universality classes for local spectral correlations: the Ginibre class and the nonstandard classes $\mathrm{AI}^\dagger$ and $\mathrm{AII}^\dagger$. We show that the continuum Dirac operator in two-color QCD coupled to a chiral $\mathrm{U}(1)$ gauge field or an imaginary chiral chemical potential falls in class $\mathrm{AI}^\dagger$ ($\mathrm{AII}^\dagger$) for fermions in pseudoreal (real) representations of $\mathrm{SU}(2)$. We introduce the corresponding chiral random matrix theories and verify our predictions in lattice simulations with staggered fermions, for which the correspondence between representation and universality class is reversed. Specifically, we compute the complex eigenvalue spacing ratios introduced recently. We also derive novel spectral sum rules.


I. INTRODUCTION
The profound relationship between random matrix theory (RMT) and natural sciences has been the subject of research over decades [1,2]. The energy-level statistics of quantum systems obeying RMT are considered to be a signature of quantum chaos [3], and the three universality classes of RMT (GOE/GUE/GSE) [4], called the Wigner-Dyson (WD) classes, provide an exhaustive description of energy-level repulsion in non-integrable quantum systems [5]. 1 In the classification of Hermitian RMT [9,10] there are seven non-WD symmetry classes. The three of them that preserve chiral symmetry exactly describe the distribution of near-zero eigenvalues of the Euclidean Dirac operator in gauge theories with spontaneous chiral symmetry breaking [11][12][13][14]. The characteristics of the seven non-WD classes are manifest only in the vicinity of the spectral origin. In the bulk of the spectrum, i.e., far away from both the origin and the spectral edge, their level statistics reduce to those of the WD classes since the chiral symmetry and the particle-hole symmetry that ensure the global pairing λ ↔ −λ of the spectrum have no effect on local correlations in the bulk. Indeed, lattice simulations of QED and QCD have revealed agreement between non-chiral RMT and the spectral statistics of the Dirac operator in the bulk [15,16]. It is therefore necessary to distinguish "symmetry classes" from "universality classes" 2 -the correct statement is that for Hermitian RMT there are 10 symmetry classes and 3 universality classes in total.
Recent years have witnessed a surge of interest in non-Hermitian quantum systems [21][22][23]. The symmetry classification of non-Hermitian RMT has been completed [24,25]. There are 38 distinct symmetry classes. 3 Out of 1 RMT ensembles that interpolate between WD classes have also been discussed [6][7][8]. It turns out that a constant nonzero breaking of antiunitary symmetries of GOE and GSE always leads to GUE statistics in the large matrix-size limit [6]. 2 In mathematics, the universality of RMT refers to the invariance of the spectral statistics for general matrix potentials [17][18][19][20]. 3 Earlier classifications [26,27] contained 43 classes. The arXiv the 38 classes, 22 possess chiral symmetry, and 3 of those describe the spectral statistics of the Dirac operator with baryon chemical potential, as reviewed in [28][29][30]. They are summarized below.

Class
Matrix form Matrix elements chiral real Ginibre 0 A B 0 real chiral complex Ginibre complex chiral symplectic Ginibre quaternion real Despite the diversity of non-Hermitian RMT, only one universality class has been known until recently: the Ginibre universality class [31]. For real, complex and quaternion Ginibre ensembles, the short-range eigenvalue correlations are always determined by the Vandermonde determinant |∆({z k })| 2 = i<j |z i − z j | 2 (see [32] for a summary of results). In [33,34], agreement was found between the bulk spectral statistics of the Dirac operator with baryon chemical potential and the Ginibre ensemble.
It has recently been pointed out [35,36] that there are in fact three universality classes for local level statistics of non-Hermitian RMT: Here, the naming scheme of [36] was adopted, T denotes transposition, and σ 2 is the second Pauli matrix. These three ensembles exhibit distinct nearest-neighbor eigenvalue spacing distributions [35,36]. Complex symmetric random matrices have received little attention from physicists for a long time, with the exception of [37]. In view of the aforementioned developments it is natural to ask the following questions. First, is there any version of [26] has been replaced with a corrected version that contains 38 classes. Hence there is no contradiction between [26] and [24]. application of the 19 nonstandard non-Hermitian chiral RMT classes to high-energy physics? Second, can the non-Ginibre universal statistics be observed in the Dirac spectrum? We answer both questions in the affirmative. Our contributions in this paper are as follows. (1) We show that Dirac operators in two-color QCD and adjoint QCD coupled to a chiral U(1) gauge field (including an imaginary chiral chemical potential as a special case) have special transposition symmetries. (2) We introduce two novel chiral RMTs that share the above symmetries. (3) Through extensive lattice simulations with staggered fermions and using the complex spacing ratio analysis introduced in [38] we find strong numerical evidence that spectral correlations of complex Dirac eigenvalues in these theories match those of the non-Ginibre universality classes of non-Hermitian RMT. (4) We highlight the physical importance of the Dirac operator coupled to a constant chiral gauge field.
Our findings are in stark contrast to the Dirac spectrum at finite baryon chemical potential, which exhibits correlations in the Ginibre class. This work provides a new point of view on quantum chaos in gauge theories.
This paper is organized as follows. In Sec. II we derive transposition properties of the continuum Dirac operator with fermions in real and pseudoreal representations of SU (2) in the presence of a chiral U(1) gauge field. In Sec. III we present the non-Hermitian random matrix ensembles corresponding to these symmetries and derive an integral representation for the fermionic partition function. In Sec. IV we define the staggered lattice Dirac operators to be used in the numerical simulations and derive their transposition symmetries. We also derive novel spectral sum rules for these operators. In Sec. V we present our numerical results, which verify the correspondence between university classes and fermion representations in SU (2). We conclude in Sec. VI. Three appendices are provided for a discussion of constant chiral gauge fields and for technical details.
Let us first consider quarks in the fundamental representation of the gauge group SU (2). The Euclidean Dirac operator coupled to a chiral U(1) gauge field is given by where A a ν and B ν are the SU(2) and U(1) gauge fields, respectively, and the τ a are the generators of SU(2), i.e., the Pauli matrices acting in color space. The U(1) charge assignment should be such that gauge anomalies are absent. As is well known, in the absence of B ν , D possesses the antiunitary symmetry [39] [iD, where C = iγ 4 γ 2 is the charge conjugation matrix and K is the complex conjugation operator. The coupling to B ν breaks this symmetry, but D retains the transposition symmetry As Cτ 2 is a symmetric unitary matrix, Theorem 3 in [40] can be applied, and it follows that D = D T in a suitable basis. It is straightforward to verify that the argument so far holds for an arbitrary pseudoreal representation of the gauge group. Next, we turn to quarks in a real (e.g., adjoint) representation of a non-Abelian compact gauge group. In the absence of B ν , D has the antiunitary symmetry [39] [iD, CK] = 0 .
With the coupling to B ν , the antiunitary symmetry is broken, but D fulfills the transposition symmetry (In a real representation, the generators can always be chosen to be antisymmetric, see, e.g., [41,Eq. (5)].) Hua's decomposition [42,Theorem 7] implies that one can find a unitary matrix U such that C = U Σ 2 U T with where σ 2 acts on the Dirac indices. Performing a similarity transformation D → U * DU T and using C = −C * it follows that D can be cast into a complex quaternion form satisfying D T = Σ 2 DΣ 2 . The following table gives a summary.

Representation
Note that, if there is nonzero baryon chemical potential, the Dirac operator D(µ) ≡ D − µγ 4 no longer satisfies (3) since D(µ) T = Cτ 2 D(−µ)Cτ 2 . In contrast, in the case of an isospin chemical potential, where the t a are the generators of SU(2) flavor. Thus, by a rerun of the argument above, D I (µ) T = Σ 2 D I (µ)Σ 2 in a suitable basis.
If we drop the spatial components of B ν and make the temporal component constant, the latter represents an imaginary chiral chemical potential, We comment on the physical implication of this term in App. A. Recently [43] has investigated the effect of a term analogous to (8) in the NJL model to formulate PTsymmetric chiral symmetry breaking. We note in passing that a real chiral chemical potential has been extensively used in lattice simulations to create chirality imbalance [44][45][46], see [47] for a review.

A. Definition of ensembles
On the basis of the transposition symmetries (7) we propose that the spectral properties of these Dirac operators are described by the following non-Hermitian random matrices:

Representation
Matrix form Matrix elements The matrix elements of V obey independent Gaussian distributions with the same variance and mean zero. The ensemble in the first row corresponds to class AI † with sublattice symmetry S + (equivalent to class D with sublattice symmetry S + ), and the ensemble in the second row to class AII † with sublattice symmetry S + (equivalent to class C with sublattice symmetry S + ), respectively [24, Tables VII, XII and XIII]. Note that if V is real or quaternion real, these ensembles reduce to chiral GOE and chiral GSE, respectively. Every eigenvalue of the ensemble in the second row is doubly degenerate due to the non-Hermitian generalization of Kramers' theorem [24].
In section V we perform extensive numerical analyses to verify our proposal. In this paper, we focus on bulk properties of the spectrum and leave the analysis of spectral statistics near the origin to future work.

B. Fermionic partition function
The fermionic partition function, i.e., the average of the product of characteristic polynomials, is of special interest in RMT. The arguments for pseudoreal and real quarks can be run in parallel, so we shall concentrate on pseudoreal quarks in the following. At first sight it may seem natural to consider the partition function where N f is the number of quark flavors of mass m, but this Z is pathological as it vanishes in the chiral limit m → 0. Instead we must consider K pairs of conjugate quarks [48], where w f and z f are masses. We consider even N . The determinants may be expressed as a Grassmann integral, After integrating out V we introduce auxiliary bosonic fields to perform a Hubbard-Stratonovich transformation. The Grassmann variables can then be integrated out trivially to yield where After a singular-value decomposition the saddle point is given by Ω = 1 2K . This implies spontaneous symmetry breaking U(2K)×U(2K) → U(2K) diag [11]. The soft mode around the saddle point may be parametrized as Ω = U ∈ U(2K). Then This is the most generic result in the large-N microscopic limit. Let us consider the special case With the skew-symmetric matrix This integral has been evaluated analytically in [49], which enables us to obtain where Pf denotes the Pfaffian and I is the modified Bessel function of the first kind.

IV. LATTICE DIRAC OPERATORS
In this section we define the lattice Dirac operators that will be used in the numerical simulations. To keep the notation simple we will use the same symbol D as for the continuum Dirac operator since confusion is unlikely to arise. We work in four dimensions. All dimensionful quantities are given in lattice units, i.e., we set the lattice spacing a to 1.
A. Staggered Dirac operator with chiral U(1) field The massless staggered Dirac operator coupled to a lattice gauge field U µ (x) is given by its matrix elements where µ is a unit vector in one of the four dimensions and the η µ (x) are the usual staggered phases. We consider SU(2) gauge fields but do not specify the representation for the time being. The continuum chiral symmetry is replaced by the remnant chiral symmetry {D, ε} = 0, where the operator ε is given by its matrix elements ε xy = ε(x)δ xy with Since D is anti-Hermitian its eigenvalues are purely imaginary. It is straightforward to show that due to {D, ε} = 0 the eigenvalues of D come in pairs ±λ with eigenvectors ψ and εψ. We now introduce a chiral U(1) gauge field where ϕ µ (x) is a smooth real gauge field. The staggered Dirac operator coupled to θ µ (x) reads The U(1) gauge transformations of the fermions and the gauge field are given by where α(x) is an arbitrary gauge potential. For θ µ (x) = 1, D(θ) is no longer anti-Hermitian so that its eigenvalues move from the imaginary axis into the complex plane. The remnant chiral symmetry {D(θ), ε} = 0 still holds, and therefore the complex eigenvalues still come in pairs ±λ with eigenvectors ψ and εψ.
We now consider the SU(2) gauge field in the fundamental and adjoint representation, denoted by U F and U A , respectively. The latter is obtained from the former by It is straightforward to check that U A is real. All statements made in this section so far hold for both representations. The difference between the representations lies in the transposition symmetries. Using the pseudoreality (U F ) * = τ 2 U F τ 2 it is straightforward to show that As a consequence, the eigenvalues of D F (θ) are twofold degenerate. We show this in App. B. Similarly, for the adjoint representation we use (U A ) * = U A to show that We therefore expect D(θ) in the fundamental or adjoint representation of SU(2) to be in universality class AII † or AI † , respectively. 4 This expectation will be confirmed numerically in Sec. V. We note that the universality classes of the fundamental and adjoint representation are interchanged for staggered fermions compared to the continuum. This is a generic phenomenon [50][51][52].
The eigenvalues of D(θ) satisfy a sum rule that is very useful as a check of the eigenvalues obtained numerically. We show in App. C that where the average is over all links of the lattice, V is the lattice volume and N rep is the dimension of the representation, i.e., for gauge group SU(2) we have N rep = 2 (N rep = 3) in the fundamental (adjoint) representation, respectively.

B. Staggered Dirac operator with chiral chemical potential
To introduce the chiral potential µ 5 we follow Braguta et al. [45] but slightly change their notation to conform with our conventions. We write the staggered Dirac operator in the presence of µ 5 in the form where D is the operator in (18) and D 5 is given by its matrix elements with s(x) = (−1) x2 , δ = (1, 1, 1, 0) and In the continuum limit, the second term in (27) reduces to the µ 5 term in (8) [45]. It is straightforward to show that D 5 is anti-Hermitian and anticommutes with ε. Thus, for µ 5 ∈ R the eigenvalues of D(µ 5 ) are purely imaginary, while for µ 5 / ∈ R they are generically complex. Furthermore, we still have {D(µ 5 ), ε} = 0, and thus the eigenvalues of D(µ 5 ) again come in pairs ±λ with eigenvectors ψ and εψ. Equations (24) and (25) also hold with D(θ) replaced by D(µ 5 ), and hence we expect the same correspondence between representation and universality class. For SU(2) fundamental every eigenvalue of D(µ 5 ) is again twofold degenerate, see App. B.
The eigenvalues of D(µ 5 ) also satisfy a useful sum rule. We show in App. C that where the average is over all sites of the lattice.

V. NUMERICAL RESULTS
The aim of this section is to test our proposal by numerical simulations of the corresponding lattice quantum field theories. As mentioned above, we use staggered fermions. We work in the quenched approximation since it is numerically cheap and sufficient for the point we are making. The inclusion of dynamical fermions does not change the universality class and therefore makes no difference for the quantities we compute.

A. Overview of the simulations
Our simulations were done in the gpt framework [53], which provides a convenient python interface to the Grid library [54]. We used the simple (unimproved) Wilson plaquette action for the gauge fields. The SU(2) gauge field in the fundamental representation was generated using the Creutz-Kennedy-Pendleton heatbath algorithm [55,56]. The adjoint SU(2) field was obtained from the fundamental field via (23). For the generation of the compact U(1) field exp(iϕ µ (x)) we used the Hattori-Nakajima heatbath algorithm [57]. The coupling constants were chosen as β SU(2) = 4/g 2 = 2.0 and β U(1) = 1/e 2 = 0.9, which in both cases corresponds to the confined phase. Numerically, the most expensive operation is the calculation of all eigenvalues of the Dirac operator, which scales like O(V 3 ). Fortunately, for the comparison of lattice data and RMT predictions, a relatively small volume is sufficient since the local spectral correlations have very small finite-volume effects. We found that a lattice size of 8 3 × 16 already shows nearly perfect agreement with RMT and therefore restrict ourselves to this lattice size. In this case the eigenvalue calculation can be done on a single CPU in a very short time, i.e., we do not need to employ HPC resources. The number of configurations for each case is given in Table I.
Scatter plots of the Dirac eigenvalues for a single configuration are shown in Fig. 1 for the cases we considered. The spectral correlations on the scale of the mean level spacing are insensitive to the global features of the spectrum. The main purpose of Fig. 1 is to demonstrate that the distribution of the eigenvalues is smooth, i.e., there U(1) µ5 = i µ5 = 2i SU (2) (21), while the last two columns give the value of µ5 in (27).
are no clustering effects. Note that for a fixed lattice volume the number of independent eigenvalues for the adjoint representation is larger by a factor of 3 compared to the fundamental representation since (i) the dimension of the representation is larger by a factor of 3/2 and (ii) the eigenvalues in the fundamental representation are twofold degenerate.
Once the eigenvalues λ k have been obtained, we follow    Ginibre   FIG. 3. Results for P (r) (left) and P (θ) (right) from staggered lattice data for SU(2) fundamental and from RMT obtained by integrating the distributions P (z) in Fig. 2 over θ and r, respectively. The histograms have been plotted as lines for better visibility (horizontal values correspond to bin centers). The four curves for the AII † ensemble and the lattice data are essentially on top of one another. In the bottom row we plot the same data as in the top row, but relative to the Ginibre ensemble for better visibility. [38] and compute the complex spacing ratios where λ NN k and λ NNN k are the nearest and next-to-nearest neighbor of λ k in the complex plane. We used a k-d tree algorithm to perform the nearest-neighbor search. We included all eigenvalues in the computation of the z k . This causes some boundary effects, see the discussion in [38], but these effects are suppressed in large volumes.
Writing z = r e iθ we have r ≤ 1 by construction. 5 We computed the distributions P (z), P (r) and P (θ) as well as moments of these distributions from the lattice data and compared them with the corresponding RMT predictions. To the best of our knowledge, analytical RMT results for these quantities are not available yet. We therefore generated the RMT predictions numerically by diagonalizing 300 random matrices of dimension 4000 for every ensemble.
Ginibre -U(1) FIG. 5. Jensen-Shannon distance between the distributions P (z) for SU(2) fundamental in Fig. 2, computed in polar coordinates as a function of the number of bins for r and θ.
the figure captions, and therefore we only summarize our observations here. From the heatmap plots in Fig. 2 we cannot draw unambiguous conclusions, but we already see that the lattice data for P (z) look roughly equal for all three cases and agree with the AII † class, while the AI † and Ginibre classes look different. From Figs. 3 and 4 it becomes quite clear that the lattice data correspond to class AII † . An interesting observation in Fig. 3 (and Fig. 7 below) is that the curves for all three ensembles appear to intersect at the same points. In the absence of analytical results we do not know whether this is actually true, but perhaps this can be shown in future work. We also computed the Jensen-Shannon distance between pairs of the two-dimensional distributions shown in Fig. 2, comparing lattice results to RMT predictions as a function of the number of bins. The results are displayed in Fig. 5. The Jensen-Shannon distance is an information-theoretic measure for how similar two prob- ability distributions are. The fact that the distance extrapolates to zero for RMT class AII † strongly supports the agreement of the data with that universality class.

C. Results for SU(2) adjoint
Our numerical results for gauge group SU(2) and fermions in the adjoint representation are summarized in Figs. 6 through 9. The conclusions are completely analogous to those in Sec. V B, except that we now observe agreement with RMT class AI † .

VI. CONCLUSIONS AND OUTLOOK
We have shown that the two nonstandard universality classes AI † and AII † of non-Hermitian RMT are realized in the bulk spectral correlations of the Dirac operator in SU(2) gauge theory coupled to a chiral U(1) gauge field or an imaginary chiral chemical potential. In the continuum, we find class AI † for pseudoreal representations of SU(2) and class AII † for real representations. We have established the corresponding chiral RMTs and verified our predictions in lattice simulations, where the above correspondence between representation and universality class is reversed for the staggered Dirac operator. We have also derived spectral sum rules that are very useful to check the numerical eigenvalues.
There are several avenues for future work. Analytically, it would be interesting to study the integral repesentation (15) in the general case of nonequal masses. Also, one should try to derive exact RMT result for the spectral correlations, which we have generated numerically so far. On the lattice side, one could study the bulk spectral correlations in the deconfined phase, for which the same universality classes are expected. One could also focus on the local spectral correlations of the near-zero eigenvalues, for which the chiral symmetry is relevant.

ACKNOWLEDGMENTS
This work was supported in part by DFG grant SFB TRR 55.

Appendix A: Constant chiral gauge fields
Spatially inhomogeneous bilinear condensates arise in QCD and QCD-like theories under various conditions [58][59][60][61][62][63][64]. To measure a modulated chiral condensate one needs to add a source term to the action, which is the Fourier component of the chiral condensate with a wave vector q. Its nonzero value signals spontaneous breaking of chiral, translational and rotational symmetries. The chiral limit is assumed here. By performing a space-dependent chiral rotation ψ → e −iγ5q·x ψ the exponential factor in (A1) may be absorbed, while a new term arises in the Lagrangian, which has the same form as the chiral vector field B ν in (1). The accumulation of near-zero eigenvalues of the modified Dirac operator D − iγ 5 γ k q k is linked to the formation of modulated condensates.

Appendix B: Twofold degeneracy
It is well-known [65] that the eigenvalues of the staggered Dirac operator in the fundamental representation of SU(2) are doubly degenerate since SU(2) is pseudoreal. We briefly show that this is also true in the presence of the chiral U(1) field θ µ (x) and the chiral chemical potential µ 5 ∈ C, for which D is no longer anti-Hermitian. Our arguments for this non-Hermitian Kramers degeneracy closely parallel those in [24]. For greater generality we formally include both θ µ (x) and µ 5 in the staggered operator. We first consider the massless case and discuss the mass term at the end of this section. We split the massless operator in the form where D 1 /2 and D 5 are given in Eqs. (21) and (28), respectively. Both operators satisfy (24), hence so does D. Now let |ψ be a right eigenvector of D with eigenvalue λ and ϕ| be the corresponding left eigenvector with the same eigenvalue, i.e., we have Taking the transpose of the second equation and using (24) we find D(τ 2 |ϕ * ) = −λ(τ 2 |ϕ * ). Using {D, ε} = 0 then yields Hence λ is doubly degenerate with eigenvectors |ψ and ετ 2 |ϕ * , provided that these two vectors are linearly independent. To show this we use (ετ 2 ) T = −ετ 2 and obtain ϕ|ετ 2 |ϕ * = ϕ|ετ 2 |ϕ * T = − ϕ|ετ 2 |ϕ * , i.e., we have ϕ|ετ 2 |ϕ * = 0. Together with the third equation in (B2) this implies that |ψ and ετ 2 |ϕ * are biorthogonal and hence linearly independent. Adding a mass term m1 to D simply shifts all eigenvalues by m, which preserves the twofold degeneracy.
Appendix C: Derivation of spectral sum rules In this section we give a brief derivation of the sum rules for the squared staggered Dirac operator. Again, for greater generality we include both the chiral U(1) field θ µ (x) and the chiral chemical potential µ 5 and now also a mass term in the operator. The special cases considered in the main text, i.e., the massless staggered operator with either a chiral U(1) field or nonzero µ 5 , are obtained from the general sum rule (C6) by setting m = 0 and either µ 5 = 0 or θ µ (x) = 1.
We split the massive operator D m into three parts, where the unit matrix has dimension N rep V and the operators D 1 /2 and D 5 are given in Eqs. (21) and (28). We trivially have Tr D 1 = 0 and Tr D 5 = 0 because the diagonal elements of these operators are zero. We also have Tr D 1 D 5 = xy Tr(D 1 ) xy (D 5 ) yx = xyµ η µ (x) Tr U µ (x)θ µ (x)δ x+µ,y − U † µ (y)θ µ (y)δ x,y+µ × s(y) Ū δ (y)δ y+δ,x +Ū † δ (x)δ y,x+δ = 0 (C2) since after the summation over y we end up with Kronecker deltas of the form δ x±µ,x+δ , which are always zero because δ = ±µ. 6 Hence We now derive the last two terms. First, where we have used η µ (x)η µ (x ± µ) = 1 and the fact that U is unitary. In the last line, · · · xµ denotes an average over all links of the lattice. Furthermore, Tr D 2 5 = xy s(x)s(y) Tr Ū δ (x)δ x+δ,y +Ū † δ (y)δ x,y+δ × Ū δ (y)δ y+δ,x +Ū † δ (x)δ y,x+δ = x s(x) s(x + δ) TrŪ δ (x)Ū † δ (x) where we have used s(x)s(x±δ) = −1 and · · · x denotes an average over all sites of the lattice. Note thatŪ δ (x) is not unitary. From (C3) we obtain In our numerical simulations µ 5 is purely imaginary. 6 Here and below we assume that the lattice has more than two sites in every direction because otherwise some Kronecker deltas could be nonzero due to the toric nature of the lattice.