Heavy-tailed chiral random matrix theory

We study an unconventional chiral random matrix model with a heavy-tailed probabilistic weight. The model is shown to exhibit chiral symmetry breaking with no bilinear condensate, in analogy to the Stern phase of QCD. We solve the model analytically and obtain the microscopic spectral density and the smallest eigenvalue distribution for an arbitrary number of flavors and arbitrary quark masses. Exotic behaviors such as non-decoupling of heavy flavors and a power-law tail of the smallest eigenvalue distribution are illustrated.


Introduction
Random matrix theory has flourished as a versatile tool in theoretical and mathematical sciences over decades [1][2][3]. In Quantum Chromodynamics (QCD), the development of chiral random matrix theory (chRMT) [4][5][6] (also called the Wishart-Laguerre ensemble) that extends the traditional Wigner-Dyson classes has helped us gain a profound understanding of the link between dynamical mass generation of fermions and spectral statistics of the Dirac operator [7]. The chRMT has also been used as a simplistic Ginzburg-Landau-type model of QCD at finite temperature and density [8,9]. Exact spectral correlations of a non-Hermitian Dirac operator at nonzero chemical potential were also worked out [10][11][12][13][14].
In RMT there are various choices for the probabilistic weight of random matrix elements. While the independent Gaussian distribution is the simplest from a mathematical point of view, it often turns out that distributions that deviate from Gaussian lead to the same spectral correlations in the limit of large matrices. This robustness of RMT is known as universality [22]. However, when the deformation of the weight is strong enough, results begin to differ from those of the Gaussian ensemble. Such random matrix ensembles with heavy-tailed weights have found applications to disordered conductors and financial statistics, as reviewed in [23]. In general, when the rotational invariance of matrices is broken as in the Lévy matrix ensemble [24], the models tend to be analytically intractable. Alternatively one can also consider heavy-tailed matrix ensembles with rotational invariance, at the expense of losing statistical independence of matrix elements. Both directions have been actively pursued [25][26][27][28][29][30][31][32][33][34][35], revealing a plethora of exotic behaviors not seen in Gaussian RMTs.
So far, applications of chRMT to QCD have been mostly limited to the hadronic phase with ψψ = 0. Extensions to the high-density regime where the diquark condensate ψψ = 0 preponderates were explored in [13,14,36,37]. While spontaneous symmetry breaking in all these cases can be characterized by a nonvanishing fermion bilinear condensate, symmetry breaking in general can also be triggered by higher-order condensates. In QCD it was stressed by Stern that chiral symmetry breaking with F π = 0 does not necessitate ψψ = 0 [38,39]. Indeed one can imagine a situation where chiral condensate is forbidden by an anomaly-free discrete subgroup of U(1) A and the spontaneous breaking SU(N f ) R × SU(N f ) L → SU(N f ) V is driven by a quartic condensate. (This pattern of symmetry breaking was studied by Dashen long time ago [40].) While this exotic phase that we call the Stern phase is ruled out by rigorous QCD inequalities at vanishing baryon density [41], there are arguments in favor of the Stern phase at finite density. First, color superconducting phases in dense QCD are examples of the Stern phase due to the fact that the leading gauge-invariant order parameter that breaks chiral symmetry is provided by four-quark condensates [42,43]. Secondly, in phases with spatially modulated chiral condensates, the phonon fluctuations associated with translational symmetry breaking wipe out the spatial order and lead to a phase with quartic condensates [44]. Other related arguments can be found in [45][46][47][48][49][50].
In this paper, we introduce a new heavy-tailed chRMT that corresponds to the Stern phase. To be precise, we show that our chRMT with N × N random matrices reproduces, in the large-N limit, the finite-volume partition function of the Stern phase with K > 4 in the ε-regime. (Here we label the Stern phase with an index K that specifies the unbroken subgroup of U(1) A [50].) This implies that all infinitely many sum rules for the Dirac eigenvalues in the Stern phase [45,50] are obeyed by microscopic eigenvalues of random matrices in the new chRMT. In the chiral limit, our chRMT coincides with the model previously considered by Akemann and Vivo [33]. Here we solve the model at large N with arbitrary quark masses and analytically obtain the microscopic spectral density and the smallest eigenvalue distribution in dependence of quark masses. We will not compute the macroscopic large-N limit or spectral densities at finite N because they are not universal quantities.
This paper is structured as follows. In Section 2 we define the model and discuss its relevance to QCD. Then we solve the model analytically in the large-N limit and obtain the microscopic spectral density and the smallest eigenvalue distribution, for an arbitrary number of flavors and arbitrary quark masses. Section 3 is devoted to conclusions and outlook.
2 Random matrix theory for the Stern phase of QCD

Definition of the model and the large-N limit
The matrix model considered in this paper is defined by the partition function where X is a complex N × N random matrix and dX denotes the flat Cartesian measure. This integral converges for arbitrary N ≥ 1 and N f ≥ 0. The weight (2.1) evidently has three interesting properties: (i) it is invariant under unitary rotations X → V 1 XV 2 for V 1,2 ∈ U(N ), (ii) the matrix elements are statistically correlated, and (iii) the distribution is heavy tailed, i.e., it does not decay exponentially for large matrix elements. This random matrix ensemble can be seen as an unquenched generalization of previous RMTs [28-30, 32, 33, 35] that had a heavy-tailed weight similar to (2.1) but with no determinants. If the weight for X in (2.1) is replaced with a Gaussian distribution, the model reverts to the standard β = 2 chRMT called the chiral Gaussian unitary ensemble (chGUE) [4,5]. In the large-N limit, our chRMT enjoys a sigma-model representation. To see this, we rewrite (2.1) up to a trivial multiplicative constant as where in the last step we introduced W ≡ zX. The final expression (2.2) is akin to the normal chGUE except that the mass term is multiplied by another Gaussian random variable z of order 1/ √ N . This implies that we need a large-N limit withm f = O(1/ √ N ), which is different from the conventional large-N limit withm f = O(1/N ). Following the standard route of bosonization [4], one can easily obtain where dU denotes the Haar measure andM ≡ diag(m 1 , . . . ,m N f ). Equation (2.4) exactly coincides with the ε-regime finite-volume partition function of QCD in the Stern phase with K > 4 [41,45,50]. Notably, the sigma model (2.4) has no term at O(M ) in the exponent, in contradistinction to the standard chiral Lagrangian. 1 This reflects that the chiral condensate ψψ → 0 in the chiral limit. We find it intriguing that the sigma model structure has changed in the absence of any changes in symmetries of the underlying random matrix. The new chRMT may serve as a toy model for spontaneous symmetry breaking driven by higher-order condensates. It follows from the coincidence of mass dependence between chRMT and QCD that infinitely many spectral sum rules of the Dirac operator in the Stern phase [45,50] can be reproduced exactly from chRMT. This suggests that the universal behavior of small Dirac eigenvalues originating from chiral symmetry breaking could be probed by using this chRMT, which shares the same pattern of symmetry breaking as the Stern phase of QCD but is much simpler and analytically tractable. We end this subsection with two supplementary remarks.
• The partition function (2.4) has no dependence on the gauge-field topology. In the Stern phase with K > 4, topological sectors with nonzero winding numbers are suppressed in the leading order of the ε expansion [50]. This deprives us of a physical motivation to study the model (2.1) with rectangular X. Nevertheless, it could be mathematically interesting to investigate such extensions in future work. 2 • The sign of the leading term in the exponent of (2.4) was fixed unambiguously by chRMT, despite that both signs are allowed by symmetries. Actually, the sign of the leading term can be flipped if we modify the fermion determinant in (2.2) as However the Dirac operator is now Hermitian! This means that the anti-Hermiticity of the Dirac operator imposes a nontrivial constraint on the sign of the low energy constant. A similar observation was made in chRMT for Wilson fermions [53].

Microscopic spectral density
When ∀m f ∈ R, the partition function (2.2) becomes with the Dirac operator Our primary interest is in the spectral statistics of D S on the scale 1/ √ N . In this "microscopic domain", the eigenvalue density and eigenvalue correlations are expected to be universal in the sense that it is solely determined by the pattern of global symmetry breaking, with no dependence on specific details of UV interactions in QCD that cause the symmetry breaking. In the following, we derive the microscopic spectral density in the large-N limit, first in the chiral limit (Sec. 2.2.1) and then for arbitrary quark masses (Sec. 2.2.2), by making use of a formal similarity of (2.6) to the chGUE. Our notation is fixed as follows.
Obviously, λ S n = λ n |z| for every n.

Chiral limit
Let us begin with the chiral limit. The spectral density of D S at finite N is defined as where we have introduced the spectral density in massless chGUE (2.12) The microscopic limit of (2.12) was derived in [5]. Now we use this result to obtain the microscopic spectral density for the new chRMT, where we have introduced the microscopic spectral density for massless chGUE [5] 16) ρ S N f (ζ) The integral in (2.15) can be performed analytically, resulting in a compact expression Here I n (z) is the modified Bessel function of the first kind. In Fig. 1 we show ρ S N f (ζ) for several values of N f . It converges to 1/ √ π = 0.564 . . . as ζ → ∞. As N f increases, the density of eigenvalues near zero is depleted because of the determinant in the measure.
is flat and monotonic (except for N f = 0), showing no oscillatory behavior typical of ρ N f (ζ). Let us recall that, in the chGUE, the oscillation is produced by peaks in the density of individual small eigenvalues. By contrast, as we will see in Sec. 2.3, the density of individual eigenvalues near zero in our chRMT is so broad that their superposition smears out each peak completely. Aside from this difference, ρ N f (ζ) and ρ S N f (ζ) look similar, but once again we emphasize that ρ S N f (ζ) is the density at the scale λ S n ∼ 1/ √ N , whereas ρ N f (ζ) is the density at the scale λ n ∼ 1/Nthese two regimes are totally different. It was originally pointed out by Stern [38,39] that chiral symmetry could be spontaneously broken when near-zero Dirac eigenvalues scale as 1/ √ V 4 instead of 1/V 4 in the thermodynamic limit V 4 → ∞. Our finding within chRMT is fully consistent with Stern's perspective. While the Stern phase is ruled out by rigorous QCD inequalities in QCD at zero density [41], its realization in zero-dimensional chRMT is not prohibited. As a side remark, we mention that a universal behavior of Dirac eigenvalues on the scale 1/ √ V 4 is known in strongly non-Hermitian chRMTs corresponding to the BCS regime of dense QCD-like theories [13,14,21,37,54].
Usually one associates the approach of ρ N f (ζ) to a constant value at ζ → ∞ with a nonzero chiral condensate through the Banks-Casher relation [7]. It must be noted, however, that the same behavior of ρ S N f (ζ) does not imply a nonvanishing chiral condensate. The reason is that, in this model, the height of the macroscopic spectral density R S N,N f (λ) at the origin scales as √ N for N 1, implying that the chiral condensate lim An important remark on a preceding work is in order. In [33], with a mathematical motivation, Akemann and Vivo studied a deformed Wishart-Laguerre ensemble, which is essentially equivalent to (2.1) with N f = 0. They derived the microscopic spectral density in the large-N limit analytically, not only for square X but also for rectangular X of size (N + ν) × N . In addition, their analysis exhausted all the three symmetry classes with Dyson index β = 1, 2 and 4. Intriguingly, one can show that their model with ν > 0 and β = 2 exactly coincides with our unquenched model (2.1) with N f = ν and ∀m f = 0. This is a manifestation of what is known in chRMT as the duality between flavor and topology [17]. As a result, (2.15) above can be obtained from [33, eq.(4.7)] by letting ν = N f and sending α → 0 there. We confirmed that our results agree with [33]. Nonetheless we have presented the full derivation above, firstly because the compact expression (2.17) is new, and secondly because the computation in the chiral limit is a useful preliminary step for the generalization to the case of arbitrary nonzero quark masses, which is a genuinely new result of this paper and will be worked out in the next subsection.

Nonzero masses
Reinstating quark masses in (2.10) and replacing z by z/ √ N , we obtain The microscopic limit is achieved by sending N to infinity with

Extending (2.13) to nonzero masses, we obtain
where we have introduced the partition function for chGUE (cf. [55][56][57]) whereas the spectral density for massive chGUE is given by . (2.27) The microscopic limit of (2.27) was computed in [58,59] as . In the chiral limit (2.29) reduces to (2.16). Now, substituting (2.25) and (2.28) into (2.22), we finally arrive at the microscopic spectral density of D S with N f massive flavors, Let us examine the simplest case closely. For N f = 1, (2.30) can be simplified to In the limit µ → 0 (2.32) reproduces (2.15) for N f = 1, as it should. We plot ρ S N f =1 (ζ, µ) in Fig. 2 for several values of µ. Clearly ρ S N f =1 (ζ, µ) increases with µ. The asymptotic value at ζ 1 depends on µ, and appears to diverge as µ → ∞. This means that a heavy flavor does not decoupleρ S N f =1 (ζ, µ) does not reduce to the quenched density at large µ. 3 This is quite unusual compared to what is known for standard chRMT [58,59], where ρ S 2 (ζ, {µ 1 , µ 2 }) the microscopic spectral density approaches 1/π asymptotically for any number of flavors and any masses, and where the decoupling of heavy flavors holds in the sense that, when some of the masses are sent to infinity, the massive spectral density reduces to that for a reduced number of flavors. By contrast, Fig. 2 reveals that neither property persists in the Stern phase. For comparison, we also display the massive spectral density for N f = 2 in Fig. 3, which exhibits a similar mass dependence to N f = 1. We look into this curious behavior in more detail in the next subsection.

Large-mass limit
Why does a heavy flavor fail to decouple? Let us examine what happens to the spectral density when one of the masses is made large compared to the others. Our starting point is the ε-regime partition function (2.4) of the Stern phase for N f light flavors, If µ N f is by far the largest among µ f 's, the fluctuation of U over U(N f ) would be effectively restricted to U(N f − 1), hence U Ũ 0 0 1 withŨ ∈ U(N f − 1). By plugging this into (2.33) and introducing a reduced mass matrix µ r ≡ diag(µ 1 , . . . , µ N f −1 ), we get which is nothing but the partition function of chGUE with N f − 1 flavors [4]. Hence one cannot recover Z S N f −1 from Z S N f by sending one of the masses to infinity; this is how our naive expectation of decoupling fails. Instead, one ends up with the conventional chiral Lagrangian with an O(M ) term whose coefficient is set by µ N f . This implies that a large explicit mass µ N f induces large chiral condensates ψ f ψ f for the other N f − 1 flavors. Generation of such induced condensates has been discussed in [46] for large-N c QCD, and our analysis based on chRMT is totally consistent with [46].
As there is a generic correspondence between sigma models and spectral statistics, one can expect that the spectral density for the Stern phase at large µ N f would reduce to that of chGUE whose Gaussian distribution parameter is set by µ N f . In fact, when This can be shown from (2.30) by using the Laplace expansion of a determinant and approximating the modified Bessel function by its asymptotic form. The relation (2.36) explicitly provides a novel link between the spectral density in the Stern phase and that in chGUE.
To assess the accuracy of (2.36), we display ρ S N f =1 (ζ, µ) for µ = 5 and 10 in Fig. 4, together with the RHS of (2.36). While the agreement is good for small ζ, deviations emerge for ζ 1/µ. An oscillatory behavior not present in the chiral limit gradually sets in as µ increases.

Smallest eigenvalue distribution
Next we turn to the smallest eigenvalue distribution in the large-N microscopic limit of the chRMT for the Stern phase. As in the previous sections, we work with the rescaled masses µ f = √ Nm f . We first define the so-called gap probability  38) with P N,N f the smallest eigenvalue distribution. Now, by applying the method of [59], it is somewhat tedious but straightforward to show where S N,N f and R N,N f are the partition function and the spectral density of chGUE, respectively, as defined in (2.23) and (2.27). Then it is easy to take the large-N microscopic limit by exploiting (2.25) and (2.28), with the result where . . . . . .
This is a new result. For small N f , integrals in (2.43) can be carried out analytically and yield simple expressions: They are correctly normalized to 1 when integrated over 0 ≤ ζ ≤ ∞. The result for P 0 agrees with [33]. In P 2 the masses were set equal for simplicity. A salient feature of (2.46) is that they decay only polynomially (∝ ζ −3 ) at large ζ, in contrast to a Gaussian decay in chGUE [59][60][61]. This long tail of P N f (ζ) could be a signal of weak eigenvalue repulsion in this model. Actually, the decay ∼ ζ −3 can be shown for any N f and any masses, on the basis of (2.40). If we rescale the variable as x → x/ζ in the numerator of (2.40), we get an additional overall factor ζ −4 while the rest of the integral tends to a well-defined large-ζ limit. Combined with ζ at the head of (2.40), the prefactor becomes ζ −3 .
In Fig. 5, P N f (ζ) for N f = 1 and 2 are plotted and compared to the microscopic spectral density (2.30). P N f nicely fit the near-zero part of the spectral density. They tend to be more localized near the origin and represent a peak in the density when the masses are increased, as anticipated from the reduction to chGUE discussed in Sec. 2.2.3.
We expect that the extension of our analysis to the k-th smallest eigenvalue distribution for general k ∈ N would be straightforward along the lines of [62].

Conclusions and outlook
In this paper we proposed an unorthodox chiral random matrix model with a heavy tail. The model is a one-parameter reweighting of the standard chGUE and can be solved exactly. We discussed potential relevance of this model to the Stern phase of QCD, where chiral condensate is zero but chiral symmetry is broken by higher-order condensates. We analytically obtained the microscopic spectral density and the smallest eigenvalue distribution in the large-N limit and discussed their dependence on the number of flavors and quark masses. Our model is not only useful as a conceptual toy model for the Stern phase but may also help a numerical evaluation of low-energy constants in future lattice simulations through fitting to the lattice Dirac spectrum.
There remain several issues that call for further investigation. While we only solved the model with unitary symmetry (β = 2), it would be technically straightforward to generalize it to β = 1 and 4; in fact, this has already been done for the case of N f = 0 in [33]. One can also build a non-Hermitian extension of this model in the spirit of [10] by replacing the matrix X in (2.1) by a sum of two independent random matrices. While these extensions are intriguing at least from a mathematical point of view, their physical meaning remains to be clarified. Another unanswered question is whether the present model could be applied to QCD in three dimensions. To the best of our knowledge, there is so far no study of a Stern-like phase in 2 + 1 dimensions. We wish to address some of these issues in future work.