Stressed Cooper pairing in QCD at high isospin density: effective Lagrangian and random matrix theory

We generalize QCD at asymptotically large isospin chemical potential to an arbitrary even number of flavors. We also allow for small quark chemical potentials, which stress the coincident Fermi surfaces of the paired quarks and lead to a sign problem in Monte Carlo simulations. We derive the corresponding low-energy effective theory in both p- and ε-expansion and quantify the severity of the sign problem. We construct the random matrix theory describing our physical situation and show that it can be mapped to a known random matrix theory at low baryon density so that new insights can be gained without additional calculations. In particular, we explain the Silver Blaze phenomenon at high isospin density. We also introduce stressed singular values of the Dirac operator and relate them to the pionic condensate. Finally we comment on extensions of our work to two-color QCD.

1 Introduction Understanding the nonperturbative physics of Quantum Chromodynamics (QCD) is one of the central challenges in theoretical physics.QCD in the vacuum is strongly coupled, giving rise to a variety of emergent phenomena such as chiral symmetry breaking, quark confinement, formation of nuclei, and mass gap generation of gluons.Since the seminal work by Banks and Casher [1] it is known that chiral symmetry breaking is associated with the condensation of near-zero eigenvalues of the Dirac operator.The correlations of Dirac eigenvalues on the scale ∼ 1/V 4 Σ, also known as the microscopic domain, strictly obey the predictions of chiral random matrix theory (ChRMT), which corresponds to the leading order of the ε-expansion of chiral perturbation theory (ChPT) [2,3] (see [4,5] for reviews).
Here, V 4 and Σ stand for the volume of Euclidean space-time and the chiral condensate in the chiral limit, respectively.The equivalence between a rather simple Gaussian matrix model with no space-time structure and QCD in a certain limit is truly surprising, but it has been confirmed explicitly again and again through lattice QCD simulations.Not only theoretically intriguing, the equivalence also provides us with a means of extracting low-energy constants in ChPT from lattice QCD data, where Dirac eigenvalues are easily computable.
The dynamics of QCD at nonzero temperature T and/or chemical potential µ is relevant for the physics of the early Universe, relativistic heavy-ion collisions, and compact stars [6][7][8].At high density, the physics is entirely different from that of the vacuum: the celebrated BCS mechanism leads to the condensation of quark pairs, which breaks gauge and chiral symmetries in three-color QCD, a phenomenon referred to as color superconductivity [9,10].However, conventional Monte Carlo simulations based on importance sampling are hindered by the infamous sign problem, which originates from the complex phase of the fermion determinant at nonzero µ [11].While several promising approaches to overcome this obstacle have been proposed [11,12], a feasible way to simulate dense QCD is yet to be found.To gain insights into the physics of dense quark matter, a number of QCD-like theories that have a nonnegative path-integral measure even at nonzero chemical potential have been investigated intensively by many authors, with numerical methods as well as in effective models.Such special theories include QCD with gauge group SU(2) (called two-color QCD) [13], QCD with adjoint fermions [14], G 2 gauge theory [15], and QCD with isospin chemical potential [16,17]. 1 Those theories share many features, such as the existence of light bosons that condense at nonzero chemical potential, and the interested reader is referred to [18,19] for reviews.
In the absence of reliable numerical simulations, analytical first-principle studies are highly valuable.The study of the Dirac spectrum in QCD with nonzero quark chemical potential in the regime µ 2 q 1/ √ V 4 was undertaken in [20][21][22][23][24] on the basis of low-energy effective theories and ChRMT (see [21,25] for reviews).It was found that the sign problem is manifested in an extreme oscillation of the spectral density of the Dirac operator, and that the latter is actually responsible for the fact that observables in QCD (e.g., the chiral condensate) at T = 0 are independent of µ q below roughly one third of the nucleon mass, even though the fermion determinant itself depends on µ q .This is informally called the Silver Blaze phenomenon of QCD [26,27].The baryon-number Dirac spectrum was also studied in [28].
The microscopic Dirac spectrum in QCD and QCD-like theories at high density was investigated in [19,[29][30][31][32]. Through the extension of ChRMT to dense QCD it was shown that the fluctuations of the complex Dirac eigenvalues of order 1/ √ V 4 ∆ 2 (with ∆ the BCS gap of quarks) are universal, i.e., independent of the microscopic details of the QCD interaction and solely determined by global symmetries.The whole analysis was extended to the singular values of the Dirac operator [33].A Banks-Casher-type relation in dense QCD-like theories was also established, which connects the Dirac spectral density at the origin and ∆ 2 [34].
In this paper we consider QCD with an even number N f of flavors at large isospin chemical potential µ I Λ QCD [17,35].For two flavors and zero quark chemical potential ε-regime

Inhomogeneous phase
Single-flavor pairing Figure 1.A schematic phase diagram of N f = 2 QCD with large isospin chemical potential µ I as a function of quark chemical potential µ q at T = 0.The ε-regime will be defined in section 3 assuming that the system is placed in a four-dimensional Euclidean box of linear extent L.
the partition function is given by2 where the Dirac operator D with the property D(µ) † = −D(−µ) is defined in section 2 and the subscript YM implies an average over the gauge fields.At low T , the ground state is dominated by the Fermi sea of u and d quarks plus the condensate uγ 5 d that originates from the attractive interaction between quarks near the Fermi surface. 3This leads to a BCS gap ∆ for quarks.In [19] a low-energy effective theory at energy scale ∆ was constructed for the generalization of (1.2) to N f flavors.Furthermore, the ChRMT describing the spectrum of D(µ) was identified and solved analytically [19].
Then a natural question to ask is what happens if the condition µ q = 0 is loosened.This is a long-standing subject, and a rough physical picture is known at least for asymptotically large µ I where the weak-coupling BCS mechanism is at work (see figure 1).Namely, for small µ q = 0, the pairing between u and d quarks is stressed by the mismatch of Fermi levels, but the ground state at T = 0 is unchanged as long as µ q is too small to compensate for the energy cost of breaking the Cooper pairs.In this region, no quark number is generated and the BCS gap is independent of µ q [36].(This property will be referred to as the highdensity Silver Blaze phenomenon in the rest of this paper, to distinguish it from the original one at low density.)When µ q reaches a threshold µ c q ≈ ∆/ √ 2 (called the Chandrasekhar-Clogston limit [37,38]) the standard BCS pairing is no longer energetically preferable and a first-order phase transition occurs to an inhomogeneous phase (e.g., a Fulde-Ferrell-Larkin-Ovchinnikov phase, where the pair carries a net nonzero momentum) [39].As µ q grows further, the system is expected to undergo yet another phase transition to a state with a single-flavor pairing (u u and dd) [40,41].On the other hand, for low and intermediate µ I the physics is less transparent because the system is strongly coupled; see, e.g.,  for studies on QCD-like theories and [69][70][71] for reviews on possible inhomogeneous phases in A typical situation considered in this paper is shown for N f = 4.The chemical potentials for u quarks (d quarks) are assumed to be slightly perturbed from −µ I (+µ I ).
the phase diagram.We also note that in recent years similar physics has been discussed in the context of imbalanced ultracold atomic Fermi gases [72,73].
The high-density Silver Blaze phenomenon for 0 < µ q < µ c q is puzzling at first sight, since observables are independent of µ q while the fermion determinant det[D(−µ I + µ q ) + m] det[D(µ I + µ q ) + m] in the path-integral measure depends on µ q .In this paper we elucidate the mechanism behind this phenomenon by constructing the low-energy effective theory and the corresponding ChRMT for QCD at large isospin and small quark chemical potential, and by looking into the spectral properties of the Dirac operator.As an idealization we will neglect beta decay and the charge neutrality condition, which would strongly suppress the formation of a pionic condensate [51,56,74].
This paper is organized as follows.In section 2 we summarize basic properties of QCD with isospin and quark chemical potential and fix the notation.In section 3 we construct the corresponding low-energy effective theory in both pand ε-expansion and study the severity of the sign problem.In section 4 we construct the ChRMT corresponding to the leading order of the ε-expansion.By mapping it to a known ChRMT that is applicable at low density we can gain a number of insights at high density.We define stressed singular values of the Dirac operator and relate them to the pionic condensate, and also study the baryon-number Dirac spectrum.In section 5 we briefly comment on two-color QCD and mention which parts of the arguments for QCD with N c ≥ 3 have to be modified for N c = 2.We conclude in section 6.In Appendix A we clarify a potential ambiguity in the effective theory.

QCD with large isospin chemical potential
Assuming even N f , we consider QCD with N f /2 pairs of u and d quarks.We will refer to u quarks as u f and to d quarks as d f with f = 1, . . ., N f /2.We introduce chemical potentials of the form µ u,f = −µ I + μu,f for the u quarks and µ d,f = µ I + μd,f for the d quarks, respectively, where we assume |μ i,f | µ I for i = u, d and all f .In other words, we consider QCD at large isospin chemical potential µ I but allow for small quark chemical potentials on top of µ I .For convenience of notation we define ) An example of a chemical potential distribution for N f = 4 is shown in figure 2.
The partition function of the microscopic theory is given by where4 D(µ) ≡ γ ν D ν − µγ 4 is the Euclidean Dirac operator in the fundamental representation of SU(N c ) for N c ≥ 3,5 which is an analytic continuation of the Minkowski Dirac operator D M (µ) ≡ iγ ν D ν +µγ 0 with x 0 = −ix 4 .In this paper we always work in Euclidean space-time unless stated otherwise.We also assume sufficiently low temperature T ∆ throughout.The special case (1.1) is recovered by setting N f = 2, μu,1 = μd,1 = µ q , and Note that shifting μu → μu − δ μ1 N f /2 and μd → μd + δ μ1 N f /2 simply corresponds to a shift µ I → µ I + δ μ, as is evident from (2.2).We are not interested in such a trivial shift and therefore impose the condition which implies µ I = Tr[µ d − µ u ]/N f .Not imposing this condition leads to an ambiguity in the effective theory that is discussed in appendix A, which should best be read after section 3.1.2.
3 Low-energy effective theory

p-expansion
The purpose of this subsection is to derive a low-energy effective theory of Nambu-Goldstone (NG) bosons for the theory defined in section 2. We are interested in a regime where the coincident Fermi surfaces of u and d quarks are slightly disrupted by nonzero |μ i,f | ∆.Before looking into this case we first consider the limit μu = μd = 0 [17,35] as a starting point.

Effective theory for zero stress
As noted in [17,19,35], the symmetry breaking for (2.2) at µ I Λ QCD (and μu = μd = 0) is driven by the condensate6 u f γ 5 d f (f = 1, . . ., N f /2), resulting in the breaking pattern [19,Sec. 4 where the suppression of the axial anomaly by medium effects is taken into account.This pattern of spontaneous symmetry breaking is consistent with QCD inequalities [19].The breaking pattern (3.1) gives rise to N 2 f /2 NG bosons which we parameterize by U and V and which reside in the coset spaces We can employ the spurion method to determine the form of the effective Lagrangian L eff .Under flavor transformations of quarks, the NG bosons and the quark masses transform as where M u and M d are the mass matrices for the u and d quarks 7 and For later use, let us also insert a source term .The role of this source term is similar to that of the diquark source in two-color QCD and adjoint QCD.This term enables us to derive a Banks-Casher-type relation for the pionic condensate [33].To leave (3.4) invariant under flavor transformations, Ω 1 and Ω 2 should transform as Next we consider the parity transformation P .Recalling U ∼ d L u R and V ∼ u L d R , we have Assuming the "p-regime" counting of this theory to be8 the leading O(p 2 ) effective Lagrangian invariant under (3.3), (3.5), and parity turns out to be To better understand this result we add a few comments: 1.The non-existence of O(M u,d ) terms in L eff (U, V ) points to the fact that the chiral condensate vanishes in this theory owing to the huge energy gap of anti-quarks due to the Fermi sea.It then follows from the first term on the third line of (3.8) that the masses of the NG modes are m 2 NG = O(M 2 u,d ).In contrast, the sources Ω 1,2 appear linearly, as they couple to the condensate in this theory.This implies for the masses of the NG modes that m 2 NG = O(Ω 1,2 ).To perform a consistent low-energy expansion based on a propagator 1/(p 2 + m 2 NG ) it is natural to count M u,d as O(p) and Ω 1,2 as O(p 2 ), which explains (3.7).

Cross terms, i.e., Tr[U
, are suppressed at high density [75] and have been dropped here.Terms with a single derivative, i.e., are also allowed by symmetries, but these are total derivatives that do not contribute to the action.We note in passing that the second line in (3.8)only affects the U(1) part of U and 3. The low-energy constants in (3.8) are defined in the limit and depend on µ I .Φ is proportional to the magnitude of the pionic condensate uγ 5 d + c.c. ∆ is the BCS gap of quarks.F and f are the decay constants of the NG modes, and v and ṽ are the corresponding velocities in the medium.At asymptotically high density we have relations such as ,9 but precise knowledge of these quantities is not needed in the rest of this paper.
4. The coefficient 3N c ∆ 2 /4π 2 of the first term on the third line of (3.8) was determined in [19,34] through matching between high-density effective theory (HDET) [76][77][78] and chiral effective theory (see [75,78,79] for the corresponding analysis in the colorflavor-locked phase).The positive overall sign of this term fixes the parity of the ground state: since the minimum of this term is attained at U = −V ∝ 1 for M u,d real and positive, the ground state is odd under parity [17].If Ω 1 = −Ω 2 (a source for the 0 − condensate) the last term of (3.8) is also minimized by U = −V .However, there will be a competition if 5. The so-called Bedaque-Schäfer terms [80] are not included in (3.8) as they are subleading in the present p-expansion.

Constant terms
are not explicitly shown in (3.8) because they do not affect the dynamics of NG modes and because they are irrelevant for the analysis of microscopic Dirac eigenvalues [34].

Effective theory for nonzero stress
We now incorporate the effects of μu,d into (3.8),assuming that these chemical potentials are much smaller than the gap ∆ and can thus be regarded as low-energy expansion parameters in L eff .For this purpose we again employ the p-counting (3.9) Let us begin with the d quarks.To use HDET we momentarily switch to Minkowski space-time.The fermionic part of the microscopic Lagrangian is then given by In the regime µ I Λ QCD this theory can be treated in the framework of HDET, where we expand in powers of μd,f in a way analogous to [80] (where high-density QCD, rather than high isospin density, was considered).To second order in p the result is then given by where vν = (1, v F ) with Fermi velocity v F , D and / D ⊥ are counted as O(p), and the dots denote higher orders in p.The definitions of the projected modes d i+ (i = R, L) and of / D ⊥ are given in [76,77].The first two terms in parentheses are O(p), while the next two terms are O(p 2 /µ I ), i.e., the expansion parameter is p/µ I .
In (3.10) both μd and i∂ 0 come with γ 0 .Furthermore, in (3.11) the mass matrix and μd appear in the combination μd This implies that L HDET at this order would be invariant under a time-dependent [80].Since according to (3.3) the NG fields U and V transform in the d-quark sector as the effective theory can also be made invariant under the spurious symmetry via the replacements10 Note that the second term in parentheses is suppressed by O(p/µ I ) with respect to the first term.Therefore it can be dropped when we construct the effective Lagrangian to O(p 2 ).
The u-quark sector can be treated in a similar manner.In the end, after analytic continuation to Euclidean space-time ∂ 0 → i∂ 4 , we find the leading O(p 2 ) effective Lagrangian including the effects of μu,d to be given by where11 This completes the derivation of the effective theory in the presence of μu,d .The low-energy constants in (3.14) are the same as those in (3.8).In particular, they are defined in the limit μu = μd = 0. Equation (3.8) follows as a limit of (3.14) if we set μu = μd = 0. Let us recall that, when L eff was constructed in (3.8), terms with a single derivative such as Tr[U † ∂ 4 U ] were dropped as they are total derivatives.Retaining this term and replacing ∂ 4 U by ∇ 4 U according to (3.15a) would result in a non-derivative term, (3.16) The second term vanishes thanks to (2.3), so omission of the single-derivative terms in (3.8) does not influence our current discussion.In appendix A we discuss an ambiguity that appears if the condition (2.3) is not respected.
If we set μu = μd = µ q 1 N f /2 , representing a small common quark chemical potential on top of a large isospin chemical potential, we find that µ q disappears from the covariant derivatives in (3.15), leaving no effect on L eff .This is the high-density analogue of the Silver Blaze phenomenon we mentioned in the introduction.It could have been anticipated from the fact that the NG modes ∼ ud in this theory carry no net baryon number.We expect a nonzero baryon number to emerge only if µ q is greater than µ c q ∼ ∆/ √ 2, at which the isotropic BCS phase gives place to a new phase, but this is beyond the domain of validity of our low-energy effective theory.

ε-expansion
We now move on to the ε-regime [29,81,82].We consider the system to be confined in a 4-dimensional Euclidean box with linear extent L and volume where m NG is the mass scale of the NG fields.The first inequality ensures that the contribution of non-NG modes to the partition function is negligible, while the second inequality implies that the Compton wavelength of the NG fields is much larger than the size of the box.In this limit the partition function is dominated by the zero-momentum modes of the NG fields.This regime can be defined through the "ε-expansion" counting12 Here, ξ(x) represents the nonzero-momentum modes of U and V , which are given by U where U 0 and V 0 denote the zeromomentum modes.
Extracting the leading terms up to O(ε 4 ) from (3.14) and discarding higher-order terms we obtain In deriving (3.19) we omitted several terms at O(ε 4 ) either because they are total derivatives or because they are proportional to Tr( μu − μd ), which vanishes according to condition (2.3).
In the ε-regime, the zero-momentum modes are no longer suppressed as V 4 → ∞, and one has to sum up their contributions nonperturbatively [81].This is in contrast to the p-regime (3.9),where they are counted as O(p) like nonzero-momentum modes and can be treated perturbatively.The kinetic terms for ξ U (x) and ξ V (x) in (3.19)only affect the multiplicative normalization of the partition function and are irrelevant for the dependence of the partition function on μu,d , M u,d , and Ω 1,2 .
We thus find that the finite-volume partition function of QCD for µ I Λ QCD at leading order of the new ε-expansion (3.18) is given by This completes the derivation of the effective partition function in the ε-regime.In section 4 we will show that the expression (3.20) can be reproduced by a certain zero-dimensional random matrix theory.Note that for Ω 1 = Ω 2 = 0, (3.20) can be computed analytically in two limits: If at least one of μu or μd is zero we obtain the Berezin-Karpelevich integral [83,84].If at least one of M u or M d is zero we obtain the Harish-Chandra-Itzykson-Zuber integral [85,86].

Sign problem
Consider N f = 4 QCD with μu = μd = µ q 1 2 , Ω 1 = Ω 2 = 0, and equal mass m, i.e., Z This theory suffers from a sign problem at µ q = 0. Let us denote the complex phase of the fermion determinants inside . . .by e iθ .To estimate the severity of the sign problem it is useful to compare the partition function (3.22) with the phase-quenched (phq) theory, Z phqQCD (µ The change due to the phase quenching is shown schematically in figure 3. Then For a rough estimate it suffices to apply the mean-field approximation by dropping derivative terms in (3.14), which leads us to the microscopic limit (3.20).The result, to leading order in the thermodynamic limit (V 4 ∆ 2 m 2 1 and V 4 F 2 µ 2 q 1), is given by13 Z (4) QCD (µ I ; μu = μd = µ q 1 2 , m) Z QCD (µ I ; μu = μd = µ q τ 3 , m) where in the last step we have evaluated the integral for the configuration for arbitrary θ and ϕ to maximize the exponent in the integrand.Consequently, the sign problem is exponentially hard at any nonzero µ q , e iθ phq ∼ e −4V 4 F 2 µ 2 q . (3.28) This is in marked contrast to QCD without isospin chemical potential, where the sign problem becomes severe only for µ q m π /2 [90,91].The difference stems from the fact that Z phqQCD (µ I ; µ q , m) contains strictly massless NG modes that couple to µ q .Let us recall that at µ q = m = 0 there were eight NG modes in total.At m = 0, four of them acquire masses while the other four remain massless. 14Because two of the four massless modes are charged under the U(1) symmetry to which µ q couples in Z (4) QCD (µ I ; μu = μd = µ q τ 3 , m), they Bose-condense as soon as a nonzero µ q is turned on.This leads to the exponentially severe sign problem (3.28).By contrast, in QCD without isospin chemical potential, all pions are gapped for m = 0.This postpones the onset of the sign problem until µ q = m π /2.
The analysis of this subsection can straightforwardly be generalized to any N f divisible by 4 since in this case we can flip half of the chemical potentials and combine the Dirac determinants pairwise to obtain the absolute value.However, this is no longer possible for N f ≡ 2 mod 4.
4 Random matrix theory and spectral properties

Random matrix model for nonzero stress
A random matrix model that exactly reproduces part B of (3.20) has been constructed in [19].Here we present an extension of this model to incorporate the effects of μu,d and Ω 1,2 , where P and Q are N ×N complex matrices while μu , μd , Mu , and Md are (N f /2)×(N f /2) matrices acting on flavor indices (i.e., we write P − μu instead of P ⊗ 1 N f /2 − 1 N ⊗ μu etc. for brevity).All dimensionless parameters carry a hat to distinguish them from physical variables. 15The inclusion of the chemical potentials in this form was motivated by Stephanov's model [92], which was devised for QCD at low density.There is another wellknown way of incorporating the chemical potential into RMT devised by Osborn [20] where the chemical potential is multiplied by another Gaussian random matrix.We expect such a formulation to belong to the same large-N universality class as (4.1).As we will show shortly, our model (4.1) describes QCD at large isospin chemical potential.Models with a similar structure were investigated in [21,93] with the aim of describing QCD at small isospin chemical potential (called phase-quenched QCD by those authors).These models must not be confused with ours.It is worthwhile to note that refs.[21,93] confirmed through explicit calculation that the two formulations of incorporating µ into RMT lead to an identical quenched microscopic spectral density.This is strong evidence that these two formulations are indeed equivalent in the large-N microscopic limit.
Let us return to the model (4.1).Using standard techniques (see, e.g., [2,19,94]) of fermionization and Hubbard-Stratonovich transformation, we find that in the large-N limit with the scaling (4.1) reduces to a nonlinear sigma model, Comparing (4.3) with (3.20) we find the correspondence with the identifications This proves the equivalence of the partition function for low-energy QCD in the ε-regime and chiral RMT, both at large µ I .Let us add a few comments.
1. Just as the quark mass couples to the Dirac eigenvalues, the pionic source (3.4) couples to the singular values of the Dirac operator [33].Therefore the above correspondence, including the Ω1,2 terms, shows not only the equivalence between QCD and RMT for the Dirac eigenvalue distribution, but also for the singular-value distribution of the Dirac operator.While a complete proof would necessitate partially quenched ChPT [95,96], in this paper we shall be satisfied with the equivalence at the level of the fermionic partition function.
2. Within RMT there is no parameter corresponding to µ I .The effect of µ I is included implicitly in ∆, F , and Φ in (4.5).This is true in RMT for two-color QCD at high density as well [30,32].
3. The two partition functions in (4.4) differ by a factor e N Tr( μ2 u + μ2 d ) .This factor does not affect expectation values in both theories and is irrelevant, unless one is interested in the partition function itself, or in its derivative w.r.t.μu,d .Actually such a discrepancy generally arises when matching QCD and chiral RMT with chemical potential [28,97].[19] that Z (N f ) RMT ( Mu,d ) for μu,d = Ω1,2 = 0 may be cast into the form of the determinant of a certain matrix of dimension N f /2.Using this result one can derive the microscopic spectral density of the Dirac matrix 0 P −Q † 0 analytically for arbitrary masses [19].However, such a simple formula is not known for the current extension to nonzero μu,d .

It was shown in
In order to study the spectral properties of the Dirac matrices the first step would be to derive the eigenvalue representation of the partition function (4.1).However, this is a difficult task even in the limit Mu,d = Ω1,2 = 0, and we postpone this to future work.In the next section we turn to the singular values of these matrices to discuss the Silver Blaze phenomenon of QCD at high density, where we will find that a number of insights can be gained without any additional calculation of spectral correlations.

Mapping high to low density
In the following we set μu = μd = μq 1 N f /2 for simplicity, which satisfies condition (2.3).
From the mapping between RMT and QCD in the ε-regime found in the previous subsection we have the exact correspondence As remarked above, it is technically difficult to compute the eigenvalue correlations of these matrices.However, as will be shown below, one can analytically compute the eigenvalue correlations for the product of these matrices, Note that for µ q = 0 the operator on the LHS equals D(µ I ) † D(µ I ), 16 whose eigenvalues {ξ 2 n } are real and nonnegative.Their positive square roots {ξ n } (with ξ n ≥ 0 for all n) are called the singular values of D(µ I ).As a generalization, we will refer to the positive and negative square roots of the eigenvalues of the operator in (4.8) as the "stressed singular values".They are no longer real for µ q = 0.In the limit µ q → 0 they reduce to {+ξ n } ∪ {−ξ n }, i.e., the singular values of D(µ I ) and their negatives.We will show in section 4.4 that the stressed-singular-value spectrum encodes essential information on the pionic condensate uγ 5 d .
In the remainder of this section we concentrate on the influence of nonzero μq by setting Mu,d = 0.17 Furthermore we assume Ω1 = ω1 1 N f /2 and Ω2 = ω2 1 N f /2 from now on.Then the partition function (4.1) reads where the {±ip k } and {±iq } are the eigenvalues (2N each) of respectively.Since μq enters the Dirac matrices as an anti-Hermiticity-breaking parameter, the spectra {±ip k } and {±iq } spread from the imaginary axis to the entire complex plane, marking the emergence of the sign problem for the weight (4.9).Note that by definition the set {±p k } ∪ {±q k } constitutes the stressed singular values of the Dirac operator.We now notice an interesting fact: the measure in (4.9) consists of two components, each of which is mathematically identical to the massive partition function of RMT for QCD with N f /2 flavors at small quark chemical potential and vanishing isospin chemical potential [20,92], where the subscript ν = 0 implies the restriction to the topologically trivial sector.According to this exact correspondence, the universal microscopic correlation functions for the stressed singular values {±p k } and {±q } with weight (4.9) are precisely given by those of the well-known matrix model (4.11), provided that the pionic sources ω1 and ω2 in (4.9) are identified with the quark masses m in (4.11).The microscopic correlation functions in the model (4.11) have been computed exactly in [20] using orthogonal polynomials and in [21,93] from the replica limit of the Toda lattice equation.
In the following three subsections we present insights that can be gained from earlier works through the mapping from (4.9) to (4.11).

Microscopic stressed-singular-value spectrum
It is well known that the microscopic spectral density of the Dirac operator in QCD with µ 2 q 1/ √ V 4 and µ I = 0 changes its behavior qualitatively as a function of µ q [21,22,24].At µ q = 0 the spectral density is supported only on the imaginary axis, and its value at the origin is proportional to Σ 0 in the chiral limit, as known from the Banks-Casher relation [1].(Here, Σ 0 , F 0 , etc. denote low-energy constants of ChPT in the QCD vacuum.)For 0 < µ q < m π /2 the spectral density is roughly constant on a two-dimensional straight band along the imaginary axis, with 18width ∼ As µ q exceeds m π /2, the spectral density develops an elliptical domain of strong oscillations, with an amplitude that scales exponentially with V 4 and a period that shrinks as 1/V 4 (see figure 4 left).The spectral density is no longer real and positive, signaling the onset of a severe sign problem for µ q > m π /2.Through the mapping explained in section 4.2, these mathematical results carry over to dense QCD with µ I Λ QCD and µ q = 0.For a physical interpretation of the mathematical

Re ξ
Im ξ Figure 4. Left: Sketch of the Dirac spectral density ρ D (µ q ; λ) defined in (4.15) for QCD at small quark chemical potential µ q > m π /2.It is roughly constant in the yellow region and strongly oscillating in the blue elliptical regions (whose boundaries have been computed in [24]).Right: Sketch of the stressed-singular-value density ρ sv (µ I , µ q ; ξ) defined in (4.18) for µ I Λ QCD and µ q > ΦΩ 1,2 /F .It behaves just like ρ D (µ q ; λ), except that the real and imaginary parts are interchanged.
formulas we need to (i) trade the quark masses for the pionic sources Ω 1,2 , (ii) set the number of flavors to N f /2, and (iii) replace the low-energy constants in the QCD vacuum by those in the high-density chiral effective theory (3.14).In particular, the chiral condensate is mapped to the pionic condensate.
Instead of quoting complicated mathematical formulas from earlier works, we would like to discuss the overall structure of the stressed-singular-value spectrum.For µ I Λ QCD and µ q = 0, the square roots of the eigenvalues of the operator −D(−µ I + µ q )D(µ I + µ q ) µq=0 = D(µ I ) † D(µ I ) are the singular values of the Dirac operator D(µ I ), as explained after (4.8).The associated Banks-Casher-type and Smilga-Stern-type relations have been derived in [33,34].For 0 < µ q ΦΩ 1,2 /F ∼ Ω 1,2 ∆/g ,19 the (positive and negative) square roots of the eigenvalues of −D(−µ I + µ q )D(µ I + µ q ), i.e., the stressed singular values, extend to the two-dimensional complex plane and have a density that is roughly constant over a straight band along the real axis, with width ∼ For µ q Ω 1,2 ∆/g a severe sign problem sets in: the flat stressed-singular-value density is invaded by an elliptical domain of strong oscillations that amplify with V 4 as described above (see figure 4 right).In particular, for Ω 1 = Ω 2 = 0 the sign problem sets in as soon as a nonzero µ q is turned on, as we have seen in section 3. 3.
In this manner one can attain a quantitative picture of the microscopic domain of the operator −D(−µ I + µ q )D(µ I + µ q ) by simply translating known formulas for D(µ q ) with µ I = 0 and N f /2 flavors to the high-density regime µ I → ∞ with N f flavors.It may seem surprising that the same formulas apply to the description of two seemingly unrelated operators, in two radically distinct situations.We can interpret this finding as a notable manifestation of the universal applicability of RMT.
In the present treatment we have neglected nonzero quark masses.Understanding their effect on the stressed-singular-value spectrum is an intriguing problem that is left for future work.

Pionic condensate and stressed singular values
Let us begin with µ I = 0.For µ q = 0, the Dirac eigenvalues spread over the complex plane and the Banks-Casher relation ceases to be valid, but the Dirac spectral density is still related to the chiral condensate through the relation with the Dirac eigenvalue density20 where Tr δ(λ − A) is shorthand for i δ(λ − a i ) with a i the eigenvalues of A. At zero temperature, a general thermodynamic argument suggests that observables must be independent of µ q for µ q < µ C q M N /N c , where M N is the nucleon mass.This is referred to as the Silver Blaze phenomenon of dense QCD [26,27].Therefore the chiral condensate (4.14) must also be independent of µ q , despite the fact that ρ D (µ q ; λ) strongly varies as a function of µ q , as illustrated in the last subsection and in figure 4 left.This puzzling situation was investigated mathematically in the microscopic limit [22,23].The authors found that the explanation for the µ q -independent chiral condensate may be attained through properties of suitable orthogonal polynomials in the complex plane, which lead to nontrivial cancellations of oscillating contributions in the integral (4.14).They also realized that it is the whole spectral density, including the flat strip as well as the strongly oscillating domain, that is responsible for the correct behavior of ψψ as a function of m.
Let us see how these findings add to our understanding of high-density QCD.For µ I = 0, the partition function in the chiral limit for It follows from (3.4) that the condensate is with the stressed-singular-value density where the ξ 2 n are the eigenvalues of −D(−µ I + µ q )D(µ I + µ q ).A sketch of ρ sv is given in figure 4 right.
As is clear from (4.17) and (4.14), the relation between the stressed-singular-value density ρ sv (µ I , µ q ; ξ) and the pionic condensate is the same as the relation between the spectral density ρ D (µ q ; λ) and the chiral condensate.In the microscopic domain, the two densities are given by the same functions, as noted in the previous subsection, so all findings for the Dirac spectral density at µ q = 0 apply to the high-isospin-density regime.The µ qindependence of the pionic condensate (i.e., the high-density Silver Blaze phenomenon) at zero temperature and its discontinuity as ω crosses zero can be accounted for by the same mathematical mechanism as found for the chiral condensate in [22,23].The puzzle that the µ q -dependent function ρ sv (µ I , µ q ; ξ) leads to a constant pionic condensate is resolved in this way.
It must be emphasized, though, that the essential physics behind the Silver Blaze phenomena at low and high density is not the same.On the one hand, the QCD vacuum does not respond to small µ q > 0 since it cannot excite a nucleon.On the other hand, dense isospin matter is insensitive to small µ q > 0 because it is not energetically preferable to break the Cooper pairs of u and d quarks.It is intriguing that the same mathematical resolution applies to those two radically different situations.

Baryon-number Dirac spectrum
In this section we discuss the quark-number density n q (µ q ), which is obtained from the partition function as We again begin with the low-density regime with µ I = 0. Expressing the partition function in terms of the Dirac operator D(0) at zero chemical potential, we find Physically one expects n q (µ q ) at T = 0 to vanish for 0 ≤ µ q M N /N c .This property of QCD was discussed in connection with the spectral properties of γ 4 [D(0) + m] in [26,27].Recently this issue was revisited in [28], where ρ q (µ q , m; z) was computed explicitly for m = 0 in the microscopic limit, i.e., for λ ∼ µ q ∼ O(1/ √ V 4 F 0 ).Next we proceed to the dense regime (µ I Λ QCD ).Since the condensate uγ 5 d does not carry net baryon charge, the quark-number density must vanish identically for µ q below a threshold ∼ ∆/ √ 2 at which a first-order phase transition occurs (as reviewed in section 1).For simplicity we will only consider degenerate masses, ignore Ω 1,2 , and set μu = μd = µ q 1 N f /2 .From (2.2) we then obtain Therefore the quark-number density is given by The spectral density R q can be computed in the microscopic domain λ ∼ µ q ∼ O(1/ √ V 4 F ) using RMT.From (4.1) the corresponding random matrix can be read off as It is a challenging task to compute the spectral density of these random matrices.However, the problem simplifies considerably if we take the chiral limit m = 0, as P and Q are then decoupled.The spectral density of the simplified matrices was worked out analytically in [28] in an effort to find ρ q (µ q , m; z) at µ I = 0, cf.(4.21).The mathematical equivalence between [28] and this work enables us to extract information for R q (µ I , µ q ; z) with no additional calculation.Adapting the findings of [28] to our context, we can conclude the following.
1.At µ q = 0, the sign problem is absent and the density R q (µ I , 0; z) is positive definite.
In the macroscopic regime it varies smoothly, and in the microscopic regime it is actually constant: R q (µ I , 0; z) ∼ F 2 ∼ µ 2 I .As µ q increases from zero, a circular domain of radius µ q appears around the origin in which R q (µ I , µ q ; z) shows extremely rapid oscillations with amplitude growing exponentially with V 4 , similarly to what is observed in the Dirac spectral density at µ q = 0.
2. The quark-number density n q (µ q ) follows from R q (µ I , µ q ; z) via (4.24).If the integral is computed using only the constant part of R q , the resulting n q increases monotonically with µ q , in apparent contradiction with the expected Silver Blaze phenomenon of dense isospin matter.However, inclusion of the oscillating part of the spectrum cures this problem, and the resulting n q shows the correct µ q -independence. 21 Comment on two-color QCD While the main body of this paper concentrates on QCD with N c ≥ 3, it seems worthwhile to comment on possible extensions of this work to two-color QCD, because the finitedensity dynamics of the latter has been actively explored in lattice simulations (see, e.g., [13,98,99]).To avoid complications we will only consider the case μu = μd = µ q 1 N f /2 .First and foremost, the Dirac operator for SU(2) gauge group possesses an anti-unitary symmetry Cτ 2 γ 5 D(µ)Cτ 2 γ 5 = D(µ) * , with τ 2 the second generator of SU(2) [82].As a consequence, the partition function of two-color QCD is invariant under the exchange of quark chemical potential and isospin chemical potential [42]: The patterns of symmetry breaking with or without chemical potentials are summarized in table 1.It is notable that, unlike in QCD with N c ≥ 3, the quark chemical potential µ q in two-color QCD enters as a symmetry-breaking external field.
The unique symmetries of two-color QCD can readily be incorporated into RMT by simply replacing the complex random matrices in (4.1) with real random matrices.This prescription was introduced in chiral RMT at zero density in [100] and later generalized to chiral RMT for two-color QCD at high density [19,30,33].After applying this prescription, the mapping of section 4.2 from high to low density is still valid, and the ensuing analysis for the stressed-singular-value density and the baryon-number Dirac spectrum parallels the N c ≥ 3 case, although the actual calculations are technically more difficult [32].
We now briefly highlight some physically distinctive features of two-color QCD.As we will see shortly, the symmetry-breaking pattern essentially depends on whether N f /2 is even or odd. 22As an example for even N f /2, let us take N f = 4 with quarks {u 1 , u 2 , d 1 , d 2 }.For nonzero µ q , the Cooper pairing between u and d becomes energetically costly, so the dominant pairing channels are u 1i u 2i and d 1i d 2i with i = R, L. (Note that these condensates are color singlets for N c = 2.) Thus in this case the unbroken global symmetry that same as µ I = 0, µ q = 0 Table 1.Comparison of the patterns of spontaneous symmetry breaking in two-color QCD and in QCD with N c ≥ 3 with quark and isospin chemical potential in the chiral limit (m = 0).N f is assumed to be even.In the lower two rows the axial anomaly is ignored, as it is irrelevant at high density.In the bottom row, µ q is assumed to be much smaller the other scales (e.g., µ I and ∆) so that µ q = 0 does not disrupt the condensate at µ q = 0.
Residual symmetry (H) Sign problem Table 2. Global symmetries that remain intact after spontaneous symmetry breaking in two-color QCD with µ I = 0 and µ q = 0 in the chiral limit (m = 0).Again, µ q is assumed to be much smaller than the other scales.
leaves these condensates unchanged is [Sp(2)] 4 .It generalizes to [Sp(N f /2)] 4 for general even N f /2, as given in table 2. We note in passing that the high-density Silver Blaze phenomenon does not occur in this case, as the NG modes respond to any small µ q = 0 right away -the pionic condensate transmutes into the diquark condensates, in a way analogous to two-color QCD at low density where the chiral condensate transmutes into the diquark condensate [101].
Next we move on to odd N f /2, focusing on N f = 2 and N f = 6 for illustration.For N f = 2 and at large µ I , the condensate uγ 5 d forms and persists until µ q reaches a threshold µ c q ∼ ∆/ √ 2 (see [55,58] for detailed model analyses of this transition), while at the same time the quark-number density remains zero at T = 0, exhibiting the high-density Silver Blaze phenomenon.This theory shows essentially the same behavior as QCD for N c ≥ 3. The unbroken symmetry is [U(1)] 2 , one of which is the quark-number symmetry and the other is a rotation generated by γ 5 I 3 , with I 3 the third isospin generator.
For N f = 6, µ q = 0 tries to split the coincident Fermi levels of {u 1,2,3 , d 1,2,3 } to two levels, one each for {u 1,2,3 } and {d 1,2,3 }.However, {u 1,2,3 } or {d 1,2,3 } alone involve an odd number of flavors and cannot support an isotropic BCS pairing by themselves.Then it would be energetically more preferable to pair as u 1 u 2 , d 1 d 2 , and u 3 d 3 (up to trivial permutations).The last pairing is stressed by µ q .The residual symmetry in this phase is the product of [Sp(2)] 4 , which leaves u 1i u 2i and d 1i d 2i (i = R, L) unchanged, and [U(1)] 2 , which acts on u 3 and d 3 in the same way as in the N f = 2 case.The symmetry for general odd N f /2 is given in table 2.
The emergence of the sign problem at µ q = 0 also depends on whether N f /2 is even or odd.Since the fermion determinant in two-color QCD is real, the path-integral measure in (5.1) for even N f /2 is nonnegative definite, and therefore no sign problem arises. 23The stressed-singular-value density is a smooth function over the complex plane, unlike for N c ≥ 3 where µ q = 0 inevitably causes strong oscillations (recall figure 4).In contrast, for odd N f /2, the sign fluctuation of the determinant in (5.1) is not completely canceled at µ q = 0. Combining the mapping from high to low density in section 4.2 with the exact spectral densities in two-color QCD at low density [32] we learn that the stressedsingular-value spectrum at high density should exhibit a domain of strong oscillations just as depicted in figure 4. A quantitative study of this phenomenon in two-color QCD is an interesting future direction.

Concluding remarks
In this paper we have studied QCD with large isospin chemical potential µ I for an arbitrary even number of flavors, allowing for a small mismatch of chemical potentials for different flavors.In section 3 we have systematically constructed the low-energy effective theory of Nambu-Goldstone modes which emerge from the symmetry breaking due to the BCS pairing of u and d quarks.After formulating the p-expansion for coincident Fermi surfaces, we have extended the scheme to the case where the BCS pairing is stressed by small µ q = 0, by utilizing the invariance of the high-density effective theory under a spurious temporal gauge transformation involving µ q .We also established counting rules for the ε-expansion at high isospin density and constructed the low-energy effective theory in the leading order of this expansion.Using this effective theory we have estimated the severity of the sign problem showing that, with nonzero stress, the average sign factor becomes exponentially small for large space-time volume.In section 4 we provided a new random matrix theory that reproduces the finite-volume partition function in the ε-regime.We introduced "stressed singular values" of the Dirac operator for nonzero stress and showed that the pionic condensate at large µ I is linked to the near-zero spectrum of the stressed singular values.Moreover, we found that the microscopic correlation functions of the stressed singular values in the chiral limit at large µ I are exactly described by those of the Dirac eigenvalues at µ I = 0 and small µ q , which is a consequence of an interesting equivalence between our RMT at large µ I and the conventional one at µ I = 0 and small µ q .This equivalence also enabled us to elucidate the microscopic mechanism of the high-density Silver Blaze phenomenon: the partition function at T = 0 is independent of µ q although the quark determinant depends on µ q .We found that this is due to a rapidly oscillating part of the stressed-singular-value spectrum.Intriguingly, this feature is mathematically the same as for the Silver Blaze phenomenon at µ I = 0 and µ q = 0. Furthermore, we pointed out that the baryon-number Dirac spectrum, i.e., the spectrum of the operator γ 4 (D(µ I ) + m), can be computed analytically from our new RMT at least in the chiral limit.The extension of the present work to two-color QCD was also discussed.
There are many possible future directions.First, the Dirac eigenvalues with nonzero stress have not been considered in this work.It is an important but challenging task to compute their microscopic correlation functions explicitly in the framework of our new RMT.Second, it would be intriguing to analytically compute the group integrals of the ε-regime partition function (3.20) for the general case of nonzero mass and nonzero stress.Third, we pointed out that it is possible to obtain the stressed-singular-value spectrum and the baryon-number Dirac spectrum in the chiral limit by way of the mapping to low density.However, this mapping does not work for nonzero quark masses, and it deserves further study to understand those spectra in the massive case.Fourth, the meson mass spectrum for nonzero stress can be determined from the effective theory constructed in this paper.Fifth, it would be interesting to look into two-color QCD more thoroughly on the basis of our brief account in section 5. Sixth, the extension of this work to the regime with strong stress (µ q ∼ ∆) is quite important, but for us to formulate a low-energy expansion we must first pin down the correct pattern of symmetry breaking as well as the condensates that form.This is not yet fully resolved in dense QCD, and a lot of elaborate work would be necessary before one can discuss anything about the spectral properties of the Dirac operator.Last but not least, the results of this paper should be checked in future lattice simulations. 24Our analytical predictions are not only of physical relevance, but also offer a nontrivial benchmark test for any computational technique that aims to overcome the sign problem.
.4) into the Lagrangian, where Ω 1 and Ω 2 are (N f /2) × (N f /2) matrices and P R/L = (1 ± γ 5 )/2 are the usual chiral projectors.This allows us to extract the pionic condensate uγ 5 d + c.c. by taking the derivative of log Z