Symmetry Transition Preserving Chirality in QCD: A Versatile Random Matrix Model

We consider a random matrix model which interpolates between the chiral Gaussian unitary ensemble and the Gaussian unitary ensemble while preserving chiral symmetry. This ensemble describes flavor symmetry breaking for staggered fermions in 3d QCD as well as in 4d QCD at high temperature or in 3d QCD at finite isospin chemical potential. Our model is an Osborn-type two-matrix model which is equivalent to the elliptic ensemble but we consider the singular value statistics rather than the complex eigenvalue statistics. We report on exact results for the partition function and the microscopic level density of the Dirac operator in the $\varepsilon$-regime of QCD. We compare these analytical results with Monte Carlo simulations of the matrix model.

We consider a random matrix model which interpolates between the chiral Gaussian unitary ensemble and the Gaussian unitary ensemble while preserving chiral symmetry. This ensemble describes flavor symmetry breaking for staggered fermions in 3d QCD as well as in 4d QCD at high temperature or in 3d QCD at finite isospin chemical potential. Our model is an Osborn-type twomatrix model which is equivalent to the elliptic ensemble but we consider the singular value statistics rather than the complex eigenvalue statistics. We report on exact results for the partition function and the microscopic level density of the Dirac operator in the ε-regime of QCD. We compare these analytical results with Monte Carlo simulations of the matrix model. Introduction. Random Matrix Theory (RMT) has been applied in several areas of physics. To name only a few: quantum chaos, condensed matter theory of disordered systems, and quantum information. Introductions to those and other applications can be found in [1]. There are also applications beyond physics as in economics (finance), engineering (telecommunications), mathematics, and statistics in general. The wide applicability of RMT is referred to as universality. Empirically the spectral statistics of physical systems is quite robust against variations in the distribution of matrix elements. The global symmetries of the system have the biggest impact to these statistics. There are just ten symmetry classes of Hermitian operators (including RMT) as proven by Altland and Zirnbauer [2]. When symmetries of a system are either destroyed or enhanced, the symmetry class usually changes continuously which has been a subject of active research for a long time [3]. In this work we give a complete analytical solution to such a transition, with a strong emphasis on various applications of this result to high-energy physics.
RMT was successfully applied to Quantum Chromodynamics (QCD) in the ε-regime of QCD since the 90's [4][5][6]. The reasons for the applicability of RMT in this field are that, first, QCD shares the same global symmetries with certain random matrix models and, second, the infrared eigenmodes and eigenvalues of the Dirac operator in the chirally broken phase are completely governed by these global symmetries. The two continuum theories of 3d and 4d QCD with N c > 2 colours and quarks in the fundamental representation can be understood from this point of view. The corresponding random matrix model for the former is the Gaussian Unitary Ensemble (GUE) [6] while the one for the latter is the chiral Gaussian Unitary Ensemble (chGUE) [4,5].
Two main ideas have driven the RMT approach in QCD. First, the random matrix models are simple enough to provide analytical formulas relating observ-ables like the microscopic level density of the Dirac operator to low energy constants in QCD, e.g., see [7,8]. Those observables can then be measured with the help of lattice QCD simulations and in this way one can fix the low energy constants accurately. Secondly, RMT allows us to study situations where lattice simulations are not available, e.g., because of the sign problem. In this way it was provided a way to understand the sign problem at finite baryon chemical potential [9][10][11][12] and at finite θ-angle [13][14][15].
The applicability domain of RMT with chiral symmetry is not limited to QCD but includes all areas where fermions with chirality emerge. For example, Dirac and Weyl fermions appear in the low-energy effective field theories of graphene, topological insulators, semimetals and d-wave superconductors, see the reviews [16], which can be combined with the RMT aproach to condensed matter and disordered systems [17]. Also non-relativistic fermions with quadratic dispersion may be described by chiral RMT [18].
In the present work we wish to consider the chiral random matrix distributed as Here Herm(N ) denotes the set of Hermitian N × N matrices. The partition function with N f dynamical quarks with masses m 1 , . . . , m N f is then given by 2 The differential D is the product of all real independent differentials of H 1 and H 2 . The normalization ensures that Z (0) N = 1. For varying µ the level statistics of D interpolates between GUE [3,6] at µ = 0 and chGUE [3][4][5] at µ = 1. Since these two limits are characterized by different patterns of flavor symmetry breaking [i.e., U( , the interpolation is far from trivial. The random matrix D replaces the Euclidean Dirac operator in QCD and its matrix dimension N plays the role of the spacetime volume which will be eventually taken to infinity, keeping N m 2 f and N µ 2 fixed. There are at least three major connections between the model (1) and QCD. First, (1) was proposed in [19] as a model of the 3d Dirac operator for staggered fermions. The authors of [19] employed lattice simulations to show that the microscopic spectrum of the 3d staggered Dirac operator obeys chGUE rather than GUE due to discretization effects. The possible convergence to GUE in the continuum limit can then be modeled by the model (1) with µ playing the role of the lattice spacing. The authors of [19] only compared Monte Carlo simulations of (1) with lattice simulations, finding good agreement of the spectra. The microscopic level density of D has not been derived to date. We hope that one can estimate the lattice artefacts of the staggered discretization with the help of our analytical results.
Another application of (1) is to 4d QCD at high temperature with twisted fermionic boundary conditions. It is well known [20,21] that gauge theories at high temperature undergoes dimensional reduction. In this regime the chiral condensate evaporates and RMT loses its validity for the infrared Dirac spectrum [22]. However, by judiciously choosing the boundary condition of quarks along S 1 it is possible to avoid chiral restoration up to an arbitrarily high temperature [23,24], which has a simple explanation based on instanton-monopoles [25]. We conjecture a crossover transition of the spectral statistics of the smallest eigenvalues of the Dirac operator from chGUE (4d) to GUE (3d). The kind of transition will be similar if not even exactly as in our model since chirality is preserved in this situation. A comparison of our analytical results with simulations would resolve this.
The third application of the model (1) can be found for 3d QCD at finite isospin chemical potential. The Euclidean Lagrangian is given by where D 3d is the 3d Dirac operator which is anti-Hermitian, ψ the two-flavor quark fields, µ I the isospin chemical potential, m u/d the quark masses, j the source term for the pion condensate, and τ j and σ 3 the Pauli matrices in flavour space and spinor space, respectively. One can model this system by the two-matrix model (3) with N f = 1, in close analogy to the 4d case [9]. (Note that j in (4) corresponds to m f in (3).) The matrix model can describe the near-zero singular values of the operator D 3d + µ I σ 3 whose nonzero density at the origin is necessary to support nonzero pion condensate [26]. The random matrix model for the Lagrangian (4) without the source j was introduced in [27]. The application of our model to this system is a classical situation of applying RMT to QCD. The comparison of lattice simulations with our results should fix the low energy constants and thus the size of the nonzero pion condensate.
Before closing the introduction, let us briefly review preceding works that are closely related to ours. The random matrix W is also known as the elliptic Ginibre ensemble [11,27,28] and was studied in the scattering at disordered/chaotic systems [29] as well as in 3d QCD at finite baryon chemical potential [11,27]. The authors of these works were interested in the complex eigenvalues of W while we need the singular value statistics of W .
The singular values of W are bijectively related to the eigenvalues of D and their statistics were not calculated before. Let us mention another random matrix model for the Hermitian Wilson Dirac operator [30,31] which also describes a transition between chGUE and GUE. The difference with our model is that the chirality of the random matrix is destroyed in [30,31]. This has dramatic consequences for the statistics of those eigenvalues closest to the origin. When chirality is destroyed the level repulsion from the origin becomes weaker and weaker, see [30,31]. This is not the case in our model where we preserve chirality. The microscopic level density will always drop to zero at the origin as long as µ = 0, see Fig. 1.
Mapping to Superspace. To derive the chiral Lagrangian of (1) in the ε-regime, in particular the partition function with N f flavors and the quenched microscopic level density, we employ the supersymmetry method [32]. To this aim, we consider the partition function involving k bosonic quarks and N f + k fermionic quarks with k = 0 for the ordinary fermionic partition function and k = 1 and N f = 0 for the quenched microscopic level density. The source variables are κ . , m N f ) with m j the masses of the dynamical quarks, L = ±1 the sign of the regularization ε > 0 and λ b/f the source variables for generating the level density. In the full version of this calculation [33] we consider a more general partition function where k and N f are chosen arbitrarily. In the first step we rewrite the ratio of determinants in (5) by a Gaussian integral over an N × (2k|2N f + 2k) rectangular supermatrix V . Then the integral over H 1 and H 2 is Gaussian which can be carried out explicitly, yielding For the definition of the supertrace and other objects in superspace like the superdeterminant we refer the reader to [34]. Note that we suppress the tensor notation with identity matrices meaning κ has to be understood in the (2k|2N f + 2k) dimensional superspace as κ ⊗ 1 2 and the Pauli matrices τ j as 1 k|N f +k ⊗ τ j .
In the next step we apply the superbosonization formula [35,36] and replace where the the boson-boson block U bb satisfies LU bb τ 1 = (LU bb τ 1 ) † > 0 and the fermion-fermion block U ff is unitary. The 2k × (2k|2N f + 2k) off-diagonal block η comprises only independent complex Grassmann variables.
In terms of this supermatrix, the partition function is up to a global constant equal to Here we have already chosen the rescaling from the εregime, namely κ = √ N κ and µ 2 = N µ 2 which are fixed in the large-N limit.
Chiral Lagrangian. First we want to consider the partition function with N f dynamical quarks, i.e. k = 0 [40]. There is no supersymmetry and the whole supermatrix consists of the fermion-fermion block, only, U = U ff . For N → ∞ we have to solve the saddle point equation This implies that all eigenvalues of U are ±1.
When diagonalizing U to its eigenvalues e iφ = diag(e iφ1 , . . . , e iφ N f ) we would obtain the Vandermonde term |∆ N f (e iφ )| 2 = a<b |e iφ b − e iφa | 2 as the Jacobian. This term algebraically suppresses all solutions of (9) in 1/N which do not have an equal number of ±1. Thus we integrate only over the coset which is the same as for µ = 0 (GUE, 3d QCD), see [6]. This coset can be parametrized as U = U τ 3 U † with U ∈ U(2N f ). The chiral Lagrangian contains a perturbation to the µ = 0 result [6,37] and the partition function in the large-N limit is given as The symmetry crossover is controlled by the second term. It gives nonzero masses to splitt off the Nambu-Goldstone modes, reducing the coset down to U(N f ) for large µ. We expect that the result (10) reproduces the leading order of the ε-expansion in QCD for the three areas of applications discussed above. In [33] we have even carried out the integral over U explicitly, which yielded either an N f × N f or an (N f + 1) × (N f + 1) Pfaffian determinant, depending on whether N f is even or odd, respectively. This Pfaffian structure was also obtained in finite N computations by us in [38].
Microscopic Level Density. Next we want to derive the microscopic level density in the quenched limit. Its relation to the partially quenched partition function for k = 1 and N f = 0 is In the large-N limit the saddle point equation (9) still applies though with a different dimension in the unit matrix. For this reason we use the parametrization U bb = Le −ϑ1τ3/2 e −ϑ2τ2/2 τ 1 e z112+z2τ1 e ϑ2τ2/2 e ϑ1τ3/2 , U ff = e −iϕ1τ3/2 e −iϕ2τ2/2 τ 1 e i(z312+z4τ1) e iϕ2τ2/2 e iϕ1τ3/2 , for the two diagonal blocks of U . The angles take the values z 1 , z 2 , ϑ 1 , ϑ 2 ∈ R, ϕ 1 ∈ [−π, π], and ϕ 2 , z 3 , z 4 ∈ [−π/2, π/2]. From the saddle point equation and the properties of U bb and U ff we have to expand the variables z j as follows: (z 1 , z 2 ) = (δz 1 / √ N , δz 2 / √ N ) and (z 3 , z 4 ) = (δz 3 / √ N , δz 4 / √ N ), (z 3 , z 4 ) = (π/2 + δz 3 / √ N , π/2 + δz 4 / √ N ). The latter expansion for z 3 , z 4 yields an algebraically 1/N -suppressed term. The expansion to leading order in 1/N and integration over η and δz j is quite lengthy and we omit it here. It is done in detail in [33]. Afterwards we integrate over ϕ 1 and ϑ 1 which yields the modified Bessel functions of first order, I j , and of second order, K j , respectively. Then we are able to use the relation (11). One important identity for the Bessel functions is the following [39, Eq. (9.6.1)], lim ε→0 Im(ε−iλ) ν K ν (2(ε−iλ) cosh ϑ 2 ) = π 2 λ ν J ν (2λ cosh ϑ 2 ) (13) which simplifies the result a lot. Substituting x = sin ϕ 2 and y = sinh ϑ 2 , we eventually obtain an expression for the microscopic level density of the quenched system, This is a new result. We normalized it to the asymptotics lim λ→∞ ρ(λ) = 2/π. This density is plotted in Fig. 1 and compared with Monte Carlo simulations of the random matrix model (1) for several µ. One immediately notices the different behaviour of this result compared to the one in [30,31], where the transition between chGUE and GUE was done in a way which breaks chirality, reflecting the form of the Wilson Dirac operator. Due to the preservation of chirality in our model, the GUE limit ( µ → 0) is not uniform about the origin while it is in [30,31] when the quark mass vanishes. As we have already pointed out, an exact chirality yields for a complex matrix a linear level repulsion from the origin, though the range where this is happening is of order O( µ) for small µ. However the two models also have something in common. The oscillations are increasingly suppressed when µ decreases, whereas the limit µ → ∞ to the chGUE result is uniformly approached.
Conclusions. We sketched the derivation of the chiral Lagrangian, see Eq. (10), and the microscopic level density, see Eq. (14), in the ε-regime of QCD for the random matrix model (1) which interpolates between chGUE (= 4d continuum QCD) and GUE (= 3d continuum QCD). The details of these computations are done in [33] and fi-nite N results are given in [38]. Again, we underline that the crucial difference with the related models in [30,31] is the preservation of chirality which has drastic consequences on the eigenvalues close to the origin. However we do not expect any differences for the spectral statistics in the bulk and at the soft edges. The infra-red physics of QCD in three and four dimensions, as we mentioned earlier, is directly related to the Dirac eigenvalues closest to the origin and we believe our results can be used to quantify effects like the discretization via staggered fermions in 3d or of a high temperature in 4d with twisted boundary conditions. The quantification should work in the standard way by comparison of our analytical results with lattice simulations.
The results we presented here can be readily generalized to the situation of N f dynamical quarks as well as to higher order correlation functions. In [38] we have shown that the spectral statistics describe a Pfaffian point process meaning that all k-point correlation functions as well as averages over ratios of characteristic polynomials of the random matrix D are given by Pfaffian determinants which only comprise three functions. This simplifies the analysis greatly and it carries over from finite matrix dimension N to infinite matrix size. In particular the statistics of the infra-red eigenvalues of the physical Dirac operators in QCD (say 3d staggered Dirac operators or 3d continuum theory with isospin chemical potential) should also follow this Pfaffian point process and all the consequences related to it.
Given the connection of RMT to various fields of physics we anticipate that the solution presented in the present work would be of value outside of QCD as well, especially in the field of condensed matter systems with emergent Dirac fermions.