GUE-chGUE Transition preserving Chirality at finite Matrix Size
We study a random matrix model which interpolates between the the singular values of the Gaussian unitary ensemble (GUE) and of the chiral Gaussian unitary ensemble (chGUE). This symmetry crossover is analogous to the one realized by the Hermitian Wilson Dirac operator in lattice QCD, but our model preserves chiral symmetry of chGUE exactly unlike the Hermitian Wilson Dirac operator. This difference has a crucial impact on the statistics of near-zero eigenvalues, though both singular value statistics build a Pfaffian point process. The model in the present work is motivated by the Dirac operator of 3d staggered fermions, 3d QCD at finite isospin chemical potential, and 4d QCD at high temperature. We calculate the spectral statistics at finite matrix dimension. For this purpose we derive the joint probability density of the singular values, the skew-orthogonal polynomials and the kernels for the -point correlation functions. The skew-orthogonal polynomials are constructed by the method of mixing bi-orthogonal and skew-orthogonal polynomials, which is an alternative approach to Mehta’s one. We compare our results with Monte Carlo simulations and study the limits to chGUE and GUE. As a side product we also calculate a new type of a unitary group integral.
MSC numbers: 15B52, 33C45
Random matrix theory (RMT) is a versatile tool for analyzing spectral statistics of operators like Hamiltonians in quantum chaotic and disordered systems [1, 2, 3, 4], the density operator in quantum information theory [5, 6], and the Dirac operator in Quantum chromodynamics (QCD) [7, 8]. It even allows to compare spectra of completely different systems ranging over many orders of scales. Applications of RMT can also be found beyond physics, like telecommunications, time series analysis in finance, ecology, sociology and medicine, and mathematical topics like algebraic geometry, number theory, combinatorics and graph theory. For more examples, see, e.g., . In QCD the applicability is two-fold. First, it allows to derive analytical relations between low energy constants of the chiral effective theory (non-linear -models) and spectral observables of the Dirac operator. This enables to determine the low energy constants by lattice simulations; see [7, 8]. Second, RMT can be applied to situations where the notorious sign problem impedes lattice simulations, like at finite baryon chemical potential [10, 11, 12, 13] or at finite -angle [14, 15, 16, 17].
The random matrix model we consider here is inspired by a certain type of Dirac operators. Hence, it will be of interest in QCD although we may expect applications in other areas, as well. Especially, Hamiltonians in condensed matter theory sometimes share similar or even the same global symmetries as those of Dirac operators in QCD-like theories. Our model is a Gaussian distributed, chiral, two-matrix model exhibiting statistics corresponding to the Dyson index in the bulk of the spectrum. The random matrix is explicitly of the form
with denoting the set of Hermitian matrices. The coupling parameter can be chosen real in general. However, due to the symmetries and with an arbitrary eigenvalue of , we can reduce its parameter range to . For the model is exactly the one of chGUE while for we have the spectral statistics of the singular values of the GUE. Hence the exact chiral pairs of eigenvalues are at any time present.
The model (1.1) is related to the elliptic complex Ginibre ensemble, for which the primary focus has been on the complex eigenvalue spectrum of the matrix [18, 19, 12]. Its physical application includes the scattering at disordered and chaotic systems , as well as 3d QCD at finite baryon chemical potential [19, 12]. In comparison to these works, we are interested in the singular value statistics of .
There are three applications of (1.1) in QCD. The first application is 4d QCD at high temperature. Since the early 80’s it is understood that at high temperature QCD-like gauge theories undergo dimensional reduction [21, 22, 23, 24, 25, 26]. In this regime the chiral condensate evaporates and RMT loses its validity for the infrared Dirac spectrum . However, by judiciously choosing the boundary condition of quarks along the time-like circle it is possible to avoid chiral restoration up to an arbitrarily high temperature [28, 29]. Then the dimensional crossover should manifest itself particularly strong in the smallest eigenvalues of the Dirac operator because they encode the type of spontaneous symmetry breaking. The Dirac operator of 4d QCD with more than two colors () and quarks in the fundamental representation shares the global symmetries of chGUE [30, 31, 8]. In three dimensions the symmetries are those of GUE [32, 19, 33]. Since chiral symmetry has to be always present for the Dirac operator in the 4d continuum theory, we expect the spectral statistics of (1.1).
The second application can be found in 3d QCD at finite isospin chemical potential [19, 34]. When analyzing the pion condensate ( and are the up and down quarks), one needs to introduce a source variable . Arranging the two quark fields as , the fermionic part of the Lagrangian reads111The authors of [19, 34] did not study the pion condensate . As such, they had no source term and thus had a decoupling of the two quarks.
where is the Euclidean anti-Hermitian 3d Dirac operator, are the quark masses and and are the Pauli matrices in spinor and flavor space, respectively. Let us take the chiral limit for simplicity. The resonances (zeros of the characteristic polynomials of ) in are the eigenvalues of the operator . Nonzero density of the latter at the origin is a necessary condition for the pion condensate formation . By replacing the operator by an anti-Hermitian random matrix (GUE) we arrive at the random matrix model
The relation of this model to (1.1) is similar to the relation between the Stephanov model  and the Osborn model  for 4d QCD at finite baryon chemical potential. This means that for large we have a phase transition of the model (1.4) to a phase where develops a spectral gap about the origin. Such a phase does not exist in the model (1.1). However, in the other phase where the spectral gap is closed, we have shown in  that the hard edge statistics at the origin will be the same for both models.
The third application is 3d lattice QCD for staggered fermions. It is known [37, 38, 39, 40, 41, 42] that symmetries of the staggered lattice Dirac operator do not necessarily agree with those of the continuum theory. Recently a complete classification of such symmetry shift was given for all dimensions . Towards the continuum limit the Dirac operator has to undergo a change of symmetries to reach the correct continuum theory. In  the model (1.1) was proposed as a description of the symmetry crossover of 3d staggered fermions. The comparison of lattice simulations and Monte Carlo simulations of (1.1) in  supports their idea.
Let us mention another model which interpolates between GUE and chGUE, namely of the form
This model was considered in  for the Hermitian Wilson Dirac operator (see also [44, 45] for a related model) , in particular we want to compare the results of our model with those of the Hermitian Wilson Dirac operator in its chiral limit (only then the spectral gap of is closed). The main difference of (1.5) to (1.1) is the loss of chirality. Whenever there are no exact chiral pairs of eigenvalues like . We will see in sections 5 and 6 that this difference has an immediate impact on the behavior of the eigenvalues closest to the origin.
In the present work we will analyze the model (1.1) at finite matrix dimension. For this purpose we first derive the joint probability density of the eigenvalues of (equivalent to the singular values of modulo sign) in section 2. To achieve this, we evaluate a unitary group integral of a new kind in A which generalizes the Leutwyler-Smilga integral . The joint probability density turns out to have a Pfaffian form. This is true also for the model (1.5) in  and several other two-matrix models [47, 48, 49, 50, 51, 52, 53, 54, 55, 56], where this list is by no means exhaustive. This Pfaffian form allows us to exploit general results on the method of mixing bi-orthogonal and skew-orthogonal polynomials . Those results are summarized in B and are used to find the Pfaffian point process (section 3), the kernels (section 5) and the skew-orthogonal polynomials (section 4). Explicit expressions for the skew-orthogonal polynomials are computed via the supersymmetry method [57, 58]. In section 6, we study the limits to GUE () and to chGUE () in more detail. The qualitative and quantitative difference between the two models (1.1) and (1.5) becomes clearer in the limit . We summarize our results in section 7.
2 Joint Probability Distribution
To obtain the eigenvalues of , see Eq. (1.1), we perform the diagonalization with the singular values of and . Thus the singular values of completely determine the eigenvalue spectrum of . We are interested in the joint probability distribution of at finite . For this purpose we express the distribution of in terms of as
To shorten the notation we define
Upon the singular value decomposition the measure transforms as 
where is the Haar measure of and the differential is the product of all independent real differentials of the matrix entries of . Hence we have for the joint probability distribution of
The integral over can be absorbed in the integration over . The remaining group integral can be performed with the formula derived in A.
The joint probability density is cast into the following form
for even and
for odd . The normalization constant is
and the weight functions are
For the definition of we need the integrals
Let us underline that the joint probability density (2.6) for odd can also be written with replaced by . Indeed this would be more natural from the perspective of deriving the joint probability density, see A. The difference of the two representations is that we subtracted the last row and column from the first rows and columns which does not change the Pfaffian. In this way the two-point weight is orthogonal to the constant, i.e. . The reason why we do this is because we want to pursue the ideas in  regarding the construction of the finite results via skew-orthogonal polynomials, especially those for odd . This construction differs from Mehta’s  only for odd . It has the advantage that all pairs of skew-orthogonal polynomials can be derived in the same way regardless of the parity of , while in Mehta’s construction all polynomials are of the same order since they are modified by the one of the highest order, see more in section 4 and in B.
3 Pfaffian Point Process
simplifies drastically. The masses of the bosonic valence quarks must have a non-vanishing real part, , to guarantee the integrability. Usually one sets and chooses the first masses equal to the masses of the dynamical quarks and the remaining and being the valence quark masses which might be complex as it is the case for calculating the -point correlation function.
where the indices take the values and . It is worth noting that this representation is valid both for even and odd . The first fermionic quark masses can be identified with those of the dynamical quarks, . The remaining fermionic quark masses, as well as those of the bosonic quarks, are from valence quarks.
To get the corresponding result for odd , one of the bosonic quark masses has to be taken to infinity yielding a row and a column in the Pfaffian which comprise the partition functions and , only. The Pfaffian structure (3.2) carries over to and their expressions in the hard edge limit are given in .
Another consequence of the Pfaffian form of is that the singular values build a Pfaffian point process . This means that each -point correlation function,
can be represented as a Pfaffian, see B,
The minus sign results from the arrangement of the blocks, namely the upper left corner only comprises the matrix , here. We could also arrange the columns and rows such that each entry consists of a block containing all three kernels which would absorb the overall sign. Let us underline that the three kernels have a different form for even and odd , as shown in section 5, while the structure (3.4) itself does not change.
The normalized level density is given by
We will make use of this relation in section 5.
Again this holds due to the general form of the joint probability density and does not need any detail of the considered model, as can be readily shown by the algebraic rearrangement method proposed in [62, 63, 55].
We can also include dynamical quarks with masses in the -point correlation function (3.4). This would yield a shift in the subscripts of the kernels and, additionally, we would get rows and columns comprising , and . For odd we can introduce an additional mass and send it afterwards to infinity. This would give us a further row and column with and .
Concluding this subsection, due to the very particular structure of the joint probability density of the eigenvalues of the Dirac operator most quantities can be reduced to the knowledge of only a few functions. How the quantities depend on them is independent of the parity of , only the explicit form of these few functions strongly depends on it.
We mention that similar Pfaffian structures have been derived for several other two matrix models [47, 48, 49, 50, 51, 52, 53, 54, 43, 55, 56]. Considering the fact that determinantal point processes can also be rewritten as Pfaffian ones , Pfaffian point processes seem to be more natural than determinantal point processes.
4 Skew-Orthogonal Polynomials
Random matrix ensembles having a probability weight of the form (2.5) and (2.6) can generally be solved with the method of skew-orthogonal polynomials  or a mixed version of bi-orthogonal and skew-orthogonal polynomials  (for odd ), see B. As we have seen, the quenched limit () [7, 12] is enough to consider since the theory with dynamical quarks can be easily constructed from it. Therefore we construct only the polynomials corresponding to the quenched weight.
Let us denote by
the average of a function over a complex matrix . In this definition the two parameters and are independent, which is advantageous at a particular step of the calculation below. The random matrix model (1.1) corresponds to .
with arbitrary constants which can be adjusted appropriately at the end. The polynomials are of order in and the polynomials are of order . When we set we choose the short hand notation and . Moreover we define the skew-symmetric products
for any integrable functions . The subscripts refer to even and odd .
When using the algebraic rearrangement method in , see also B, we notice that the polynomials are proportional to Pfaffians, cf. Eq. (2.7). Due to this Pfaffian structure of the polynomials they satisfy the following orthogonality relations by construction (for any )
This is the foundation of our choice for the skew-orthogonal polynomials in section 5.
Before proceeding let us find explicit representations for the two kinds of polynomials and . We first consider and follow the ideas of the supersymmetry method [57, 58]. We refer to  for a mathematical introduction to supersymmetry. In the first step we rewrite the determinant as a Gaussian integral over a -dimensional complex Grassmann-valued vector ,
We omit the overall constants at the moment since we know that the polynomials are given in monic normalization, In the next step we employ the Hubbard-Stratonovich transformation with a Hermitian matrix
After integration over we obtain
The Gaussian integral over can be performed via the identity
which is valid for any positive definite Hermitian matrix and can be proven by spectral decomposing and then integrating over each matrix entry, separately. It remains to simplify this expression when we set . For this simplification we make use of the identity
with arbitrary matrices and , several times. In the end we arrive at
Note that everything depends on , only. Hence, we can employ the superbosonization formula [67, 68, 69] and replace the integration over by an integration over a phase, . This yields after proper normalization
The contour only encircles the origin counter-clockwise. Changing we can rewrite the polynomial to
When expanding the square root in we can identify the Laguerre polynomials
This yields a more explicit expression in terms of a finite sum,
Here, we have used the floor function yielding the largest integer which is smaller or equal to .
An expression analogous to (4.12) can be derived for . In fact it can be completely derived from . Recalling the definition (4.3) we notice that the term can be pulled out of the integral such that these terms are proportional to . The term with can be generated by a derivative in . However we have, then, also to differentiate the normalization constant but this yields only a shift in the arbitrary constants . With a slight abuse of notation, we denote the new constants also by . We have
When applying this relation to the result (4.12), we find (after shifting again)
In terms of the Laguerre polynomials this expression reads
after an additional shift of the constant from to . Here we used the identities and .
When we define the quotient
with for , each pair satisfies the normalization
This can be readily checked by the Pfaffian representation (2.7) of the polynomials.
The skew-orthogonal polynomials are different for even and odd because the two-point weight changes, cf. Eqs. (2.8) and (2.10). Therefore also the explicit form of the kernels (3.6) will be different. We collect the results for even in subsection 5.1 and for odd in subsection 5.2.
The level density at finite is then
It is notable that for decreasing a discontinuity of the level density is building up at the origin. The reason for this is that the limit is not uniform, see also  where it was observed for the level density of the staggered Dirac operator in three dimensions. This can be understood by the level densities of the GUE and the chGUE. While the level density of the GUE is non-zero at the origin it vanishes linearly for chGUE, see .
Another important point is the approach to the limit compared to the GUE and the chGUE interpolation in [43, 44, 45] where chirality is broken. In  the authors considered the random matrix model (1.5). The particular form of this model implies that regardless of how small is the chirality is broken and one has a finite density at the origin. In our model (1.1) we preserve chirality which implies that we have always a linear drop off at the origin. The level repulsion reflected in this behavior results from the exact chiral pairs of eigenvalues of which feel each other and which is missing in the model (1.5). The regime were the interaction of the chiral pairs will take place is of order for small and shows up in the level density about the origin, see Fig. 1.
The third point we want to emphasize is the merging of eigenvalue peaks of for on the positive and negative line, cf. Fig. 1. The reason is that we have on average only eigenvalues on the positive and negative axis, separately, for GUE. Those are represented by peaks in the level density. For chGUE we have peaks, thus, twice as much. This is also the reason why the width of the level density for is obviously bigger than the one for ; we note that the level density is always normalized to unity. One can interpret this behavior also differently. Since we plot in Fig. 1 the singular values of one has to compare it with the level density of the eigenvalues of GUE which is a direct sum of two independent random matrices, see [33, 70, 71, 72].
In comparison to the limit , the limit seems to be less dramatic. The level density approaches this limit uniformly without any surprising features.