Enhanced contribution of the pairing gap to the QCD equation of state at large isospin chemical potential

I study QCD at large isospin density, which is known to be in the superfluid state with Cooper pairs carrying the same quantum number as pions. I solve the gap equation derived from the perturbation theory up to the next-to-leading order corrections. The pairing gap at large isospin chemical potential is found to be enhanced compared to the color-superconducting gap at large baryon chemical potential due to the $\sqrt{2}$ difference in the exponent arising from the stronger attraction in one-gluon exchange in the singlet channel. Then, using the gap function, I evaluate the contribution of the condensation energy of the superfluid state to the QCD equation of state. At isospin chemical potential of a few GeV, where the lattice QCD and the perturbative QCD can be both applied, the effect of the condensation energy becomes dominant even compared to the next-to-leading order corrections to the pressure in the perturbation theory. It resolves the discrepancy between the recent lattice QCD results and the perturbative QCD result.


I. INTRODUCTION
The better understanding of QCD at nonzero chemical potential is essential for unraveling the dynamical phenomena involving the strong interaction such as binary neutron star mergers, core-collapse supernovae, and heavy-ion collisions, from the first principles.The QCD equation of state (EoS) is the most important quantity characterizing the thermodynamics of the system.
The QCD EoS can be evaluated reliably in perturbation theory at small temperature T and large quark chemical potential µ [1][2][3][4][5][6][7].The validity range of perturbative QCD (pQCD) is limited to the range µ ≳ 1 GeV.This value corresponds to the baryon density n B ≳ 50 n sat , where n sat ≃ 0.16 fm −3 is the nuclear saturation density.It is far beyond the reach of the core density of heavy neutron stars, so the pQCD cannot be directly applied to the realistic environment although it has an indirect impact [8,9] (see, however, Refs.[10,11]).In the non-perturbative regime, lattice QCD simulations have pinned down the EoS at high T and small µ by Taylor expansion in terms of µ/T [12,13].However, the large-µ region of QCD at µ/T > 1 is generally inaccessible by the current Monte-Carlo algorithm due to the sign problem (see Ref. [14] for review).
The sign problem can be circumvented in two-flavor QCD at nonzero chemical potential by taking the same values of chemical potentials with opposite signs for u and d quarks [15].This specific setup corresponds to setting nonzero chemical potential for the (third component of) isospin, which I denote as µ I , while keeping the baryon chemical potential, µ B , to zero.There have been several lattice QCD studies of multi-pion system at nonzero µ I regarding the phase structure and thermodynamic properties [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34] (see Ref. [35] for review of some of the earlier works).* yfuji@uw.eduA recent lattice QCD calculation [34] provided the QCD EoS at T ≈ 0 and up to µ I ≃ 3 GeV.The authors were able to construct states with the quantum numbers of up to 6144 pions, which correspond to large µ I , on ensembles of gauge configurations based on an ingenious algorithm introduced in Ref. [36,37], and measured the thermodynamic properties from the extracted ground state energies of these systems.They compared their results with the pQCD calculation from Ref. [38], and they found a discrepancy between their lattice results and the pQCD calculation at µ I ≳ 2 GeV.
QCD at nonzero µ I can be regarded as the phasequenched theory, in which the complex phase of the fermion determinant is neglected, of QCD at nonzero µ B .Recently, Moore and Gorda pointed out that phase quenching works for any linear combination of chemical potentials not only for µ I .They claim that the relative difference between the phase-quenched theory and the original theory is O(α 3 s ), where α s is the QCD coupling constant, from the perturbative consideration [39].Nevertheless, the perturbative O(α 3 s ) difference between the phase-quenched theory and the original theory will not account for the discrepancy between the lattice QCD and pQCD calculations mentioned above as the latter discrepancy is noticeably larger than O(α 3 s ).Further, Moore and Gorda proposed that the phasequenched lattice calculation can be used to extract the pressure of original theory up to O(α 4 s ) corrections in the perturbation theory through the rigorous QCD inequality [40].By the use of this inequality, the phase-quenched theory places the bound on the pressure of the original theory (see Ref. [41] for an application).
The purpose of this work is to demonstrate that the EoS of the phase-quenched theory has an exponentially enhanced nonperturbative correction compared to the original theory.To this end, I take QCD at nonzero µ I as an example.In this specific theory, the nonperturbative correction to the EoS related to the phase quenching is embodied as the BCS condensation energy.
The ground state of the isospin QCD is expected to  [42] as observed in the recent lattice QCD result [34] (see Ref. [43] for a review).In the isospin QCD, u and d quarks fill up the Fermi sea, and form the color singlet quark-antiquark condensate, while in two-flavor QCD at nonzero µ B and zero µ I , u and d quarks fill up the Fermi sea, and form the color non-singlet diquark condensate leading to the color superconductivity.The former (q q channel) has stronger one-gluon attraction compared to the latter (qq channel), thus the pairing gap is exponentially enhanced at nonzero µ I .
The BCS condensation energy correction to the pQCD calculation also accounts for the discrepancy between the lattice QCD and pQCD results.The gap parameter in the isospin QCD is evaluated by using a method similar to those used to derive the diquark gap parameter in the color-superconducting phase of QCD at nonzero µ B (see Refs. [44,45] and reference therein).I note that the exponential enhancement of the BCS gap [42] and the significance of the condensation energy contribution to the EoS in the pQCD calculation [46] (see also Ref. [47] for the demonstration that the attractive interaction near the Fermi surface stiffens the EoS) were previously pointed out in the literature although there was no reliable evaluation of these effects prior to this work.
The paper is organized as follows.In the next section, I review the QCD EoS at nonzero µ.I review the path integral representation and QCD partition function, including the notion of the phase quenching and the perturbative expansion of the QCD EoS.In Sec.III, I discuss the Cooper pairing in QCD at nonzero µ I .After an overview of the Cooper pairing is given, the derivation and solution of the gap equation are presented.Section IV constitutes the main result of this paper.Based on the pairing gap obtained in the previous section, I calculate the condensation energy and show the numerical results.In Sec.V, I mention the possible implications of this work to the quark-hadron crossover.Finally, I end the paper with a summary and discussion.
Throughout the paper, I mainly follow the notations in Ref. [45], and take the number of flavors N f as N f = 2 unless otherwise stated.

II. REVIEW OF QCD EQUATION OF STATE AT NONZERO CHEMICAL POTENTIAL
In this section, I first review the Euclidean path integral representations of the QCD partition function with nonzero chemical potential along with the notion of phase quenching.Then, I discuss the perturbative expansion of the partition function and the phase-quenching effect in the perturbation theory.The problem setting in this work and the possible resolution are briefly mentioned.

A. Path integral representation of QCD partition function
Here, I introduce the notion of phase quenching, as discussed in detail in Refs.[39,40].The Dirac operator D f (µ f ) for a quark of flavor f at nonzero chemical potential µ f is defined as where the covariant derivative, / D ≡ γ µ ∂ µ + igγ µ T a A a µ , is a skew-Hermitian operator, i.e. / D † = − / D. Integrating out the fermionic degrees of freedom, the grand canonical partition function in the path integral representation is expressed as where S G is the Euclidean action of QCD in the gauge sector.In general, the fermion determinant det D f (µ f ) is complex-valued at nonzero µ f .This complex phase is the source of the sign problem preventing the Monte-Carlo simulation on the lattice.I define the phase-quenched theory by discarding the phase of the fermion determinant as Because the complex phase is absent, the phase-quenched theory is free from the sign problem.From the relation one can rewrite the phase-quenched fermion determinant as Therefore, the phase-quenched theory is equivalent to a theory with an equal number of fermions with opposite chemical potentials.The fermion determinant appears with a fractional power 1/2, so this theory is unphysical.The phase-quenched theory becomes physically sensible in two flavors with mass-degenerate u and d quarks.In this theory, u and d quarks have opposite chemical potentials µ u = µ and µ d = −µ.Usually, the quark chemical potential µ is relabeled as µ I ≡ 2µ, where µ I is a conjugate variable to the third component of isospin I 3 and it is called isospin chemical potential.The corresponding partition function is defined as From the last line, one can indeed verify this is the phasequenched version of a theory at nonzero baryon chemical potential µ B ≡ N c µ and zero µ I :

B. Perturbative expansion
Hereafter, instead of the partition function, I consider the pressure P = (T /V ) ln Z.I quote the perturbation expansion of pressure up to next-to-next-to-leading order (NNLO) in the massless limit [1,2] regularized in the MS scheme [4,48] and the coefficients at each order in α s are where β 0 ≡ (11/3)N c − (2/3)N f is the first coefficient of the QCD beta function with N c being the number of colors.I also define the pressure of the ideal quark gas as For the coupling constant α s ( Λ), I use the expression at the NNLO, and take into account the running of α s ( Λ), which is evaluated at the renormalization scale Λ.The MS scale is fixed as Λ MS ≃ 330 MeV, which is the value suggested from the N f = 2 lattice-QCD data [49,50].I set Λ = 2µ in the following calculation as 2µ is a typical hard interaction scale in the system, but there is an ambiguity in the choice of Λ.As in the conventional prescription, uncertainties associated with this ambiguity are evaluated by varying Λ by a factor 2, namely by taking 1/2 ≤ Λ/(2µ) ≤ 2. Below, I review how the effect of the phase quenching appears in the perturbative expansion as described by Moore and Gorda in [39].One can construct a Feynman rule for each determinant in the square root in the phase-quenched fermion determinant (5) in analogy to the ordinary perturbation theory.In the Feynman rules of the phase-quenched theory, an additional factor 1/2 arises from the following relation Therefore, the only difference between the perturbative expansion of the phase-quenched theory and the original theory is that the sign of the µ is reversed for half of the fermions.
The phase-quenched theory and the original theory have the same perturbative expansion up to O(α 2 s ).One can explicitly verify that by setting µ → −µ in Eq. ( 8); it does not change the expression of P .
The effect of µ → −µ appears at O(α 3 s ).This can be described schematically by going into the Nambu-Gorkov basis where ψ C = C ψ⊤ is the charge-conjugate spinor and C ≡ iγ 2 γ 0 is the charge conjugation operator.In this basis, free fermion propagators are where where The upper and lower component of Ψ gives the equivalent description of the theory.Reversing the sign of µ in the upper component of Ψ gives the equivalent description in the lower component of Ψ with the original value of µ but with the quark-gluon coupling in the conjugate representation.Therefore, flipping the sign µ → −µ is equivalent to changing T a → −T ⊤ a (= −T * a ) in the quarkgluon vertex while keeping the original value of µ in the propagator, and the effect of µ → −µ only appears in the color factor of diagrams.The difference in the diagrammatic expansion appears only at O(α 3 s ).This is due to the difference in color factors with the fundamental and antifundamental representations [39] tr This color factor appears in the diagram with three gluons attached to two fermion loops involving different flavors.

III. COOPER PAIRING IN QCD AT LARGE ISOSPIN DENSITY
In this section, I first review the Cooper pairing in the isospin QCD.Then, I derive and solve the gap equation perturbatively.

A. Overview
At large µ I > 01 , u and d quarks fill up the Fermi sphere with the radius of µ I /2 in the ground state.The Cooper instability leads to that these u and d quarks form Cooper pairs in the color-singlet, pseudoscalar, and 1 S 0 channel: where the greek letters α, β are the color indices in the fundamental representation and ∆ is the gap parameter.
Note that it has the same quantum number as π + , and this pattern of pairing is favored from the QCD inequality [42].
The gap ∆ can be calculated in a similar setup as in the diquark condensation at nonzero µ B .In the weak coupling expansion, ∆ on the Fermi surface has the form where I have truncated the perturbation series at O(1) and O(ln g).
as first pointed out in Refs.[42].Note that this value is √ 2 first pointed out in Ref. [51] (see also Refs.[52][53][54][55][56][57][58][59].Therefore the magnitude of ∆ is exponentially enhanced.For later convenience, I rewrite the remaining terms as The factor b arises from the gluon sector, in which magnetic and electric gluon exchange occurs at a large angle, so it is independent of the color structure of the condensate and whether the chemical potential being µ B or µ The factor b ′ 0 arises from the wave-function renormalization.In this work, I compute this factor for the first time at nonzero µ I , and I find and the numerical value is e −b ′ 0 ≃ 0.420.At nonzero µ B , it was calculated in Refs.[55,60], and they found with the numerical value e −b ′ 0 ≃ 0.177.Summarizing these results, the superfluid gap ∆ is concisely summarized as The gap function in QCD at µ I ̸ = 0 is exponentially enhanced compared to the color-superconducting gap in QCD at µ B ̸ = 0.This is because the attraction arising from the one-gluon exchange is stronger in the color singlet q q channel compared to that in the color antitriplet qq channel [42].
In the following two subsections, I will derive and solve the gap equation for ∆, and clarify the difference from the diquark condensation at nonzero µ B .

B. Gap equation
In this subsection, I derive the gap equation from the perturbation theory.The gap equation is very similar to that of the two-flavor color superconductor, so I will follow the formalism developed in the context of color superconductivity, in which the Nambu-Gorkov basis is employed (14).The major difference is that the Cooper pairing occurs in q q-channel, not in the diquark channel.Consequently, the Nambu-Gorkov basis ( 14) is replaced with the isospin basis [61], namely, I write down the gap equation in the isospin basis.In this basis, the free quark propagators are The inverse free propagator is expressed in the momentum space as where µ = µ I /2 is the quark chemical potential, and Λ is the energy projector The quark-gluon coupling ΨΓ µ a ΨA a µ in this basis is characterized by the following matrix Notice the difference between the quark-gluon vertex in the isospin basis (27) and the Nambu-Gorkov basis (17); in the latter case, the lower component of Γ µ a is in the antifundamental representation −T ⊤ a .The quark part of the two-particle-irreducible (2PI) action Γ can be written as the following functional of the full quark propagator S [62-65] where Γ 2 is the sum of the 2PI skeleton diagrams, and it also depends on the full gluon propagator D. Notice that this expression does not include a factor 1/2 as in Eq. ( 14) of Ref. [45], which is required to cancel the double counting in the Nambu-Gorkov basis.The ground state is given by the stationary point of Γ. From the stationarity condition δΓ[S]/δS = 0, I obtain Schwinger-Dyson equation where Σ is the quark self-energy and defined by the functional derivative of Γ 2 at the stationary point Here, I follow the conventional approximation for Γ 2 in which one truncates the infinite sum of the 2PI skeleton diagrams up to the two-loop order; this two-loop approximation for Γ 2 corresponds to a one-loop approximation for Σ.Then, the gap equation (34) becomes where Q ≡ T ωn d 3 q (2π) 3 denote the sum over the Matsubara modes and integration in the momentum space.I write the quark self-energy in the isospin basis as and the quark full propagator as The off-diagonal elements of the self-energy are related via Φ − = γ 0 (Φ + ) † γ 0 .Through the Schwinger-Dyson equation (33), one can express the diagonal and anomalous propagators in terms of the free propagator and the self-energy as The diagonal elements of the quark self-energy are calculated as [58,66] where ḡ ≡ g/(3 √ 2π) and M 2 ≡ (3π/4)m 2 g with the gluon mass being m 2 g ≡ N f g 2 µ 2 /(6π 2 ).For the off-diagonal part of the quark self-energy, I use the following ansatz for the gap matrix given the pairing pattern in Eq. ( 19) and M is the matrix in the color space By substituting all these equations in Eqs. ( 38) and (39), I obtain where the wave-function renormalization is defined as and The quasiparticle energy is defined as Now, I consider the (2, 1)-component of the gap equation (35) where D ab µν is the gluon propagator.The gluon propagator is assumed to have the color structure D ab ∝ δ ab , where the roman indices a, b, . . .are the color indices in the adjoint representation.By substituting Φ + (41) and F + (44) in the equation above, multiplying M † γ 5 Λ (+) k to both hand sides of the equation, and taking the trace, I obtain I neglected the antiparticle contribution.Hereafter, I suppress the superscript (+) as there is only a quasiparticle contribution.The color factor C F appearing in front is This color factor is associated with the 1 channel, which is the most attractive among the available color channels in the one-gluon exchange interaction between q q, and the available color channels are given by the decomposition 3⊗ 3 = 1⊕8.Note that this is twice as large as that in the 2SC phase of the two-flavor color superconductor at nonzero µ B , in which the color structure of the gap is antisymmetric M αβ = ϵ αβ3 and the color factor is This color factor is associated with the 3 channel, which is the most attractive among the available color channels in the one-gluon exchange interaction between qq, and the available color channels are given by the decomposition 3 ⊗ 3 = 3 ⊕ 6.
After splitting the gluon propagator into the longitudinal and transverse component, taking the Matsubara sum, which is the same procedure as in the nonzero-µ B case [57], the gap equation becomes where I denote εq ≡ Z(ϵ q )ϵ q , and ∆ k ≡ ∆(ε k , k).The factor b is defined in Eq. ( 23).Note that an additional prefactor 2 = N c − 1 arises in the isospin QCD by replacing the color factor (50) in the baryonic QCD by (49).I elaborate on the solution to this equation in the next subsection.

C. Solution of the gap equation
In this subsection, I clarify how the additional factor two in Eq. ( 51) modifies the solution of the gap equation by explicitly solving it.I follow the calculation at nonzero µ B presented in Refs.[57,60] (the same results can be derived using different formalisms as in Refs.[55,58,59,67]).In the equations below, all the modifications arising at nonzero µ I are underlined.Namely, if all the underlined coefficients are discarded, one recovers the calculation at nonzero µ B in Ref. [60].
I solve the gap equation ( 51) at zero temperature to obtain the gap function at the Fermi surface ∆ * ≡ ∆ q=µ .The thermal factor becomes tanh[ε q /(2T )] = 1.As for the dressed energy εq in logarithms, I approximate as ln( bµ/|ε 2 q − ε2 k |) ≃ ln( bµ/|ϵ 2 q − ϵ 2 k |) and Z(ε q ) ≃ Z(ϵ q ), which are valid at the leading order.Then the gap equation becomes In Ref. [51], Son observed that the logarithm can be replaced by max{ln( bµ/ϵ k ), ln( bµ/ϵ q )} at this order, so I make a further approximation and introduce the variables [57,60] x ≡ ḡ ln Since ∆ * ∼ µ exp(−1/ḡ), these variables scale in the ḡ expansion as The gap equation up to O(ḡ) is written in terms of these variables as ∆(x) ≃ 2x To determine the functional form of ∆(x), I take the second derivative of the gap equation ( 56) to convert the integral equation into the differential equation where the ′ symbol denotes the derivative with respect to the argument of the function.By changing the independent variable from x to z ≡ −2 1/3 (2ḡ) −2/3 (1 − 2ḡx), the equation ( 57) becomes the Airy equation The solution to this equation is given by a linear combination of the Airy functions Ai and Bi.With arbitrary constants C 1 and C 2 , ∆(z) is expressed as For later convenience, the Airy functions are decomposed into the phase and modulus as and also for the derivative of the Airy functions, I decompose as The coefficients C 1 and C 2 are fixed by the boundary conditions ∆(z * ) = ∆ * and ∆ ′ (z * ) = 0, where The solution ∆(z) and its derivative are Now, I have determined the functional form of ∆(z).For the remaining part, I need to evaluate the undetermined constant ∆ * .To this end, I set x = x * in Eq. ( 56) and change the integration variable as y → w ≡ −2 1/3 (2ḡ) −2/3 (1 − 2ḡy), then I arrive at The function w∆(w) in the integral can be replaced with ∆ ′′ (w) using the Airy equation (58).Then, the integration by part gives the following relation where z 0 = −2 1/3 (2ḡ) −2/3 (1 − 2ḡx 0 ).By substituting Eqs. ( 62) and (63) into the above equation, I obtain I expand this equation up to the next-to-leading order in terms of ḡ.In the above equation, I use the following asymptotic formulae, which are valid at weak coupling, and I get I note that I used the fact that x 0 is one order higher than x * based on the scaling behavior (55).Also, x 0 is an arbitrary scale and far from the Fermi surface, so it is expected that x 0 dependence cancels in the final expression.Indeed, the x 0 -dependence terms cancel in the equation above, and results become independent of x 0 up to O(ḡ) in the perturbative expansion.From this relation, and its expansion owing to the relation arctan By solving this quadratic equation up to O(ḡ), and using the relation x * ≡ ḡ ln(2 bµ/∆ * ), I finally arrive at

IV. PAIRING GAP CONTRIBUTION TO THE EQUATION OF STATE
In this section, I calculate the condensation energy of the BCS state up to O(g).And, I numerically evaluate the magnitude of such a correction and compare it with the lattice QCD data.

A. Calculation of the condensation energy
The physical pressure is obtained as the value of Γ[S] at its extremum.The stationarity condition implies the relation (34), and formally this relation can be rewritten as Γ 2 [S] = 1  2 tr(ΣS).By substituting this relation into the expression of Γ[S] (32) with the use of the Schwinger-Dyson equation, I obtain the pressure By substituting the bare quark propagator ( 28) and ( 29), and the full quark propagator ( 37), (43), and (44) into the above expression, and then by taking the Matsubara sum, the pressure is obtained (the detailed derivation is presented in Sec.2.4 in Ref. [68]) At zero temperature, the pressure becomes where I suppressed the superscript (+) of the dressed quasiparticle energy εk since only the (+) component has the contribution from the pairing gap ∆; the (−) component is the same as in the unpaired vacuum.Henceforth, the logarithm can be approximated as Z(ε k ) ≃ Z(ϵ k ) so that The condensation energy is defined as Now, I limit the integral around the Fermi surface where −δ ≤ k − µ ≤ δ, and neglect the momentum dependence of the gap ∆.Around the Fermi surface, the density of states is µ 2 /(2π 2 ), and I use the integration parity to limit the range of integration to [0, δ].Then the condensation energy is reduced to where ϵ The leading-order contribution to δP in the expansion in terms of the coupling constant g is where in the last line, I only kept the term that is leading in the expansion in terms of ∆/δ.The next-to-leading-order contribution to δP of O(g) is Again, in the last line, I kept the term that is leading in the expansion in terms of ∆/δ, which is found to be proportional to ∆ 2 .Among the terms in the square bracket in the last line, I only keep the last term − ln[∆/(2 bµ)].
From the scaling in Eq. ( 55), this is the only term at O(1) and the other terms are O(g).This is because of ∆ ∼ e −1/g , so ln ∆ ∼ 1/g reduces the power of g by one.One can also think of this as absorbing the ∼ ln δ term into the definition of ∆ 2 in the prefactor since ∆ has an implicit δ-dependence neglected in the derivation of ∆ above; although it is not confirmed whether the δ-dependence in ∆ and the ∼ ln δ-term match or not.Therefore, the NLO contribution at O(g) to δP is Summarizing the result, the condensation energy δP I will plot the numerical value of this term in the next subsection.

B. Numerical results
In Fig. 1, I plot the isospin matter pressure normalized with the ideal gas value (12).At µ ∼ 1 GeV (equivalently, µ I ∼ 2 GeV), it is expected that the pQCD works since the typical interaction scale in the system is large enough compared to the QCD scale Λ MS .Meanwhile, the lattice QCD calculation provides the EoS up to µ I ∼ 3 GeV [34].As one can see in Fig. 1, there is a discrepancy between these two calculations.
On the one hand, the pQCD calculation predicts P/P id ≃ 1 − 2α s /π + O(α 2 s ) < 1 and an increase in P/P id with increasing µ I (the green band in Fig. 1).For the pQCD calculation, I use the expression (8), and evaluate the scale variation uncertainty by taking 1/2 ≤ Λ/µ I ≤ 2. The lines in the green band from the bottom to the top correspond to the value Λ/µ I = 1/2, 1, and 2.
On the other hand, the lattice QCD data [34] dictates P/P id > 1 and P/P id decreases with increasing µ I (the blue and orange bands in Fig. 1).The blue and orange shaded bands marked with Lattice QCD A and Lattice QCD B are the results sampled from different ensembles; the lattice geometry is L 3 × T = 48 3 × 96 and 64 3 × 128 for the ensemble A and B, respectively.
As explained in Sec.II B, there can be a difference of O(α 3 s ) between these two within the perturbation theory, however, the discrepancy is clearly larger than O(α 3 s ).Now, I add the condensation energy of the BCS state δP (84) to the pQCD pressure.The red lines in the figure show the perturbative estimate with the Cooper pairing taken into account.For the gap function, I use the expression in Eq. (26).The red lines from the bottom to the top correspond to the value Λ/µ I = 2, 1, and 1/2.Note that this is in reverse order of the pQCD pressure without the pairing correction.The smaller Λ/µ I corresponds to the larger value of α s , and the gap is large for the larger value of α s .
One can see that the discrepancy is resolved by adding the condensation energy.I stress that this result is without any fine-tuning, and the only ambiguity in the calculation is the choice of the renormalization scale Λ.
In Fig. 2, I compare the relative magnitudes of different contributions to the pressure.Note that the x-axis is the quark chemical potential µ, not the isospin chemical potential µ I .I show the perturbative corrections to P at each order, NLO (10) and NNLO (11) relative to the ideal quark gas pressure P id (12).I also plot O(α 3 s ) contribution with α 3 s P id mimicking the higher order corrections as the full contribution at O(α 3 s ) is yet incomplete [7].This is the order of magnitude expected for the difference between the phase-quenched theory and the original theory because this difference does not contain any logarithmically enhanced contribution ∼ α s ln(α s ) as reported in Ref. [39].
I plot the condensation energy (84) by the red lines.The condensation energy contribution surpasses the dominant corrections in the pQCD at the NLO around µ ∼ 1-5 GeV, depending on the renormalization scale.Note that the scale variation uncertainty becomes larger for pQCD with the pairing contribution.
I also overlay the pairing gap contribution diquark gap in the 2SC phase in Fig. 2. The expression of the 2SC gap is As expected, this contribution is exponentially suppressed and does not contribute to bulk thermodynamics.

V. QUARK-HADRON CROSSOVER
The results presented above may also suggest the quark-hadron crossover from the BEC phase to the BCS phase.At the equation level, one can see that the pion condensate in the BEC phase changes into a BCS condensate.
At low density, from the chiral perturbation, the leading contribution to the pressure is ∝ µ 2 I , and it reads [42,69] I call it the BEC term.From the one-loop correction in the chiral perturbation theory, the term ∝ µ 4 I arises with a small prefactor ∼ 1/(4π) 2 [69].As the density becomes larger, the µ 4 I -term becomes more relevant compared to the µ 2 I -term.From the calculations above, even at µ I ≃ 2 GeV, I find there is a substantial contribution from the pairing in the BCS regime, which has the form I call it the BCS term.Although there is a µ-dependence in ∆, it varies slowly with increasing µ, so this µdependence is mild, and ∆ can be regarded roughly as a constant.Therefore, it has the same structure as the BEC term (86).Namely, the BCS term can be rewritten as ∝ f 2 ∆ µ 2 I /2 by defining the prefactor f ∆ as In the pQCD, there also arises a term with ∝ µ 4 I (12).
In Fig. 3, I show the relative magnitude of the prefactor of the BCS term, f ∆ , to the pion decay constant f π , which is the prefactor of the BEC term.I find that the value of f ∆ stays close to f π and varies slowly with increasing µ I .Note that the band in Fig. 3 corresponds to the scale-variation uncertainty.
It may also be interesting to see the behavior of these terms in the large-N c limit as discussed in Refs.[42] (see also Refs.[70,71]).At small µ I , the BEC term (86) scales as O(N c ) since f 2 π ∼ O(N c ).At large µ I , the bulk thermodynamics is dominated by P id = N c N f µ 4 /(12π 2 ) (12), which also scales as O(N c ).Therefore, the BEC and BCS regime is continuous in the N c scaling, which is similar to that at nonzero µ B [72] although the physics origin at small chemical potential is different; in the nonzero-µ I case, the N c dependence arises from the chiral symmetry breaking through the Gell-Mann-Oakes-Renner relation, while in the nonzero-µ B case, it arises from the interaction among baryons.Now I turn to the large-N c limit of the BCS term.As I will explain shortly, it scales differently from the other terms in the pressure.Unlike the BCS term in the colorsuperconducting phase, which is suppressed in the limit N c → ∞ with fixed 't Hooft coupling λ = g 2 N c [73], the BCS term at nonzero µ I is nonvanishing in the large-N c limit as this is a color singlet.The naive estimate of the gap in the large-N c by taking this limit in the expression (26) gives With this, the BCS term (87) scales as O(N 6 c ) although it is parametrically small compared to the bulk µ 4 -term for a small value of λ.The prefactor (N c /λ) 5/2 in the above expression arises from the Debye screening and the Landau damping effects in the gluon exchange.Therefore, the BCS term (87) in the large-N c limit scales differently from the pion BEC term (86) as well as the bulk µ 4 -term.
However, these effects responsible for the large prefactor (N c /λ) 5/2 are suppressed by O(N ) in the large-N c limit, therefore the naive estimate above may be modified and gives the consistent N c -counting also for the BCS term.I will justify this by the hand-waving estimate.The gap equation takes the form [45] ∆ where θ is the angle between the momenta k and q (see Eq. ( 47) for definition) and δ is a cutoff scale for the collinear divergence.At finite N c , δ arises from the Landau damping δ ∼ (∆m 2 g ) 1/3 .Here, if we assume that the confinement persists at large N c , which is the fundamental assumption in Quarkyonic matter [72], δ can be taken as the QCD scale Λ QCD , which is the characteristic scale for the confinement.Then, the solution of the gap equation becomes ∆/µ ∼ exp(−Λ QCD /λ) without the N 5/2 c factor in front.Even though the parametric dependence of the gap on λ is different, it still survives in the large-N c limit.This will give the same large-N c scaling for the BCS term as the BEC term, so the quark-hadron crossover is implied from the large-N c limit.
I also note that the gap parameter here wins over the gap of the chiral density wave, such as of the Deryagin, Grigoriev, and Rubakov (DGR) type [74][75][76] (see also [77]).It is natural to expect so as the chiral density wave uses only the part of the phase space near the Fermi surface.The detailed analysis will be reported elsewhere.

VI. SUMMARY AND DISCUSSION
In this work, I studied QCD at nonzero isospin chemical potential µ I ̸ = 0 as a specific example of the phasequenched theory.I calculated the gap parameter and the condensation energy associated with it up to the next-toleading order in the expansion in terms of the coupling constant g.I found an exponential enhancement in this contribution, and the inclusion of this nonperturbative correction to the equation of state explains the discrepancy between the lattice QCD and naive pQCD results without any fine-tuning, as shown in Fig. 1.
This implies that when extracting the perturbative coefficients of O(α 4 s ) from the phase-quenched lattice simulation, one needs to take into account the nonperturbative correction arising from the Cooper pairing.In the physical sense, the effect of the phase quenching of the fermion determinant in the partition function is interpreted as an enhancement in the pairing gap (it is evident in Fig. 2).This correction cannot be treated within the perturbation theory explained in Sec.II B. Each fermion determinant in the square root of Eq. ( 5) corresponds to quarks with positive and negative chemical potentials, and these quarks are treated separately in the perturbation theory.However, in the actual lattice QCD calculation, there is a large contribution to the thermodynamics from the mixing between these quarks with positive and negative chemical potentials.
This can be exemplified clearly by relabeling the quarks with positive and negative chemical potentials in Eq. ( 5) as u and d quarks, respectively.Then, the Cooper pair condensation ⟨ dγ 5 u⟩ mixes u and d quarks, and one can have a diagram with u and d quarks running inside already at a one-loop level.One can generalize the results in this paper from N f = 2 to even N f straightforwardly.For odd N f , the generalization is more non-trivial, namely, the gap equation and its solution presented in Sec.III should be modified to include an additional factor of 1/2 in Eq. ( 13) in the perturbation theory.It would also be interesting to see the effect of flavor mixing of different origins, for example from the instanton-induced interaction (see, e.g.[78]).
To the best of my knowledge, this is the first crossvalidation of the perturbative QCD at large density confronted with the lattice simulation.This has several implications and impacts on QCD at nonzero baryon chemical potential, and possibly on neutron star physics.
I also verified that the calculation of the pairing gap as a solution of the gap equation derived in the perturbation theory is reliable in the regime where the perturbative expansion of the partition function is valid.It is highly plausible that the evaluation of the pairing gap is reliable as well in QCD at nonzero baryon chemical potential.Posit that the weak-coupling calculation of the pairing gap is still valid around µ ∼ 1 GeV, the pairing gap at nonzero baryon chemical potential is exponentially small compared to the isospin-QCD counterpart.As a consequence, this fact implies that the color-superconducting gap does not affect the bulk properties such as the equation of state at least in the perturbative regime.It further implies that the behavior of the trace anomaly introduced in Ref. [79] is very different between QCD at nonzero µ I and µ B .Namely, the former case shows the large negative value for the trace anomaly as shown in Ref. [34] (see also Ref. [46]), which is mainly caused by the pairing gap term in thermodynamics, while in the latter case, the trace anomaly can still be positive owing to the absence of the large pairing gap term as conjectured in Ref. [79].
As a future extension of this work, changing the number of colors N c is an interesting direction, particularly taking N c = 2 and N c → ∞.In two-color QCD, one can use technology very similar to the present work to calculate the equation of state.Two-color QCD at nonzero chemical potential is a theory free from sign problem, so one can confront the lattice data (to date, there are several lattice equation of state data available, e.g.[80][81][82][83][84]).One may also expect the large diquark gap as well because the representations 2 and 2 are equivalent in SU(2) due to the pseudoreality.As I mentioned partially in the text, the expression for the pairing gap may be different in the large-N c limit and may hint at the existence of a phase transition as a function of N c .This would also deserve further investigation.

FIG. 1 .
FIG.1.Comparison of the pQCD pressure with the lattice QCD data.The bands for the pQCD results correspond to uncertainties arising from the ambiguity in the renormalization scale Λ, and dotted lines in the middle denote the common choice Λ = µI .The pressure is normalized by the ideal-quarkgas value(12).

FIG. 2 .
FIG.2.The magnitude of the perturbative and pairing corrections to P relative to the ideal quark gas value P id .The central line corresponds to the renormalization scale Λ = 2µ, and the band around it represents Λ = µ and 4µ.Note that the horizontal axis is quark chemical potential µ.

FIG. 3 .
FIG.3.Relative magnitude of the prefactor of the BCS term compared to that of the BEC term.The vertical dashed line is the value of the pion decay constant in the vacuum, fπ ≃ 93 MeV.The band shows the scale-variation uncertainty.
a Bose-Einstein condensate (BEC) phase with pion condensation at µ I > m π .As µ I increases, it undergoes crossover to the Bardeen-Cooper-Schrieffer (BCS) state with Cooper pairs carrying the same quantum number as pions arXiv:2312.11443v3 [hep-ph] 26 Mar 2024 be