3-dimensional QCD phase diagram with pion condensate in the NJL model

With the isovector coupling constants adjusted to reproduce the physical pion mass and lattice QCD results in baryon-free quark matter, we have carried out rigourous calculations for the pion condensate in the 3-flavor Nambu-Jona-Lasinio model, and studied the 3-dimensional QCD phase diagram. With the increasing isospin chemical potential $\mu_I$, we have observed two nonzero solutions of the pion condensate at finite baryon chemical potentials $\mu_B$, representing respectively the pion superfluid phase and the Sarma phase, and their appearance and disappearance correspond to a second-order (first-order) phase transition at higher (lower) temperatures $T$ and lower (higher) $\mu_B$. Calculations by assuming equal constituent mass of $u$ and $d$ quarks would lead to large errors of the QCD phase diagram within $\mu_B \in (500, 900)$ MeV, and affect the position of the critical end point.


I. INTRODUCTION
Understanding the structure of the phase diagram for the Quantum Chromodynamics (QCD) is one of the main goals of high-energy nuclear physics. Efforts are mostly devoted to exploring the 2-dimensional QCD phase diagram, i.e., in the plane of the temperature and the baryon chemical potential. Due to the sign problem [1][2][3][4][5][6], lattice QCD (LQCD) calculations are unable to provide solid information in the high baryon chemical potential region, where our knowledge on the QCD phase diagram mostly relies on low-energy relativistic heavyion collisions and effective QCD models. Experiments at RHIC-BES, FAIR-CBM, and also at NICA and HIAF, etc., have been or are to be carried out, in order to search for the signal of the QCD critical end point, which was observed in effective QCD models, such as the Nambu-Jona-Lasinio (NJL) model [7][8][9][10], the Dyson-Schwinger approach [11,12], and the functional renormalization group method [13,14].
Our knowledge on the QCD phase diagram can be extended to another degree of freedom, i.e., the isospin [15]. As the isospin chemical potential increases and reaches the mass of a pion, pions can be produced out of the vacuum and a Bose-Einstein condensate is expected to form. The formation of the pion condensate has profound effects on the QCD phase diagram [16][17][18][19][20][21][22][23][24][25][26][27][28][29]. Unlike the case at finite baryon chemical potentials, LQCD calculations do not suffer from the sign problem and can give reliable results at finite isospin chemical potentials. Intuitively, the QCD phase diagram at finite isospin chemical potentials is related to the isovector interaction in quark * Corresponding author: xujun@zjlab.org.cn matter, and the latter has ramifications in both relativistic heavy-ion collisions and nuclear astrophysics. For example, in relativistic heavy-ion collisions induced by neutron-rich nuclei at RHIC-BES, the elliptic flow splitting between π − and π + favors a finite isovector interaction in quark matter [30], and the charge susceptibility is largely affected by the isovector interaction [31]. The isovector quark interaction may also affect properties of strange quark stars [32,33]. It is of great interest to reproduce the LQCD results [34,35] at finite isospin chemical potentials but zero baryon chemical potential by varying the strength of the isovector quark interaction, and then extrapolate the calculations to finite baryon chemical potentials, thus exploring the whole 3-dimensional QCD phase diagram.
In this manuscript, we report such a study based on a 3-flavor NJL model. We reproduce the physical pion mass and LQCD results in baryon-free quark matter by adjusting the coupling constants of the scalar-isovector and vector-isovector interaction. Afterwards, we extrapolate the model to finite baryon chemical potentials, and study on the 3-dimensional QCD phase diagram in the presence of the pion condensate. We note that the constituent masses of u and d quarks are generally set to be equal in previous studies at small or vanishing baryon chemical potentials (see, e.g., Refs. [25,35,36]). In the present study for exploring the whole 3-dimensional QCD phase diagram, we carry out a rigourous calculation without this drawback. Effects from productions of other mesons, e.g., the kaon condensate, are neglected in the present study. Section II gives the main formulas for the 3-flavor NJL model, with the derivations of the Lagrangian in the mean-field approximation, the thermodynamic potential, and the quark condensates and densities detailed in the Appendices. Results in baryon-free and baryon-rich quark matter as well as the 3-dimentional QCD phase diagram are given in Sec. III. We conclude and outlook in Sec. IV.

II. THEORETICAL FRAMEWORK OF 3-FLAVOR NJL MODEL
A. The Lagrangian We start from the Lagrangian density of a 3-flavor NJL model expressed as [37] where are the kinetic term, the scalar-isoscalar term, the vectorisoscalar term, the Kobayashi-Maskawa-t' Hooft (KMT) term, the scalar-isovector term, and the vector-isovector term, respectively. In the above, ψ = (u, d, s) T represents the 3-flavor quark fields with each flavor containing quark fields of 3 colors;m = diag(m u , m d , m s ) is the current quark mass matrix for u, d, and s quarks; λ a (a = 1, ..., 8) are the Gell-Mann matrices in SU (3) flavor space with λ 0 = 2/3I; G S and G V are respectively the scalarisoscalar and the vector-isoscalar coupling constant; G IS and G IV are respectively the scalar-isovector and the vector-isovector coupling constant. Since the Gell-Mann matrices with a = 1, 2, 3 are identical to the Pauli matrices in u and d space, the isovector couplings break the SU (3) symmetry while keeping the isospin symmetry. K denotes the strength of the six-point KMT interaction [38] that breaks the axial U (1) A symmetry, where 'det' denotes the determinant in flavor space, i.e., with Γ = 1 ± γ 5 , and ǫ ijk being the Levi-Civita symbol with q i , q j , and q k representing the u, d, and s quark fields. In the present study, we employ the parameters m u = m d = 3.6 MeV, m s = 87 MeV, G S Λ 2 = 3.6, KΛ 5 = 8.9, and the cutoff value in the momentum integral Λ = 750 MeV/c given in Refs. [7,39], and define R IS = G IS /G S and R IV = G IV /G S as the reduced scalar-isovector and vector-isovector coupling constant, respectively.
B. The Lagrangian density and the thermodynamic potential in the mean-field approximation To study the system at finite chemical potentials and temperature, we introduce the chemical potentials in the Lagrangian density or equivalently where µ B , µ I , and µ S are the baryon, isospin, and strangeness chemical potential, respectively.
Based on the mean-field approximation as detailed in Appendix A, the above Lagrangian density can be written as where is the inverse of the quark propagator S(p) as a function of quark momentum p, with being the gap parameter, and being the condensation energy independent of the quark fields. In the above, ρ q = qγ 0 q and σ q = qq are the net-quark density and the chiral condensate, respectively, with q = u, d, s, and π = ψ iγ 5 λ 1 ψ is the pion condensate. The Dirac effective mass or the constituent mass of quarks can be expressed as Note that M u = M d is used in some previous studies for the 2-flavor system [18,34] or at µ B ∼ 0 [25,35,36], while we consider the most general case in the present study on the QCD phase diagram at high baryon and isospin chemical potentials. The effective chemical potentials can be expressed as Similarly, the effective baryon, isospin, and strangeness chemical potentials arẽ Starting from the partition function as detailed in Appendix B, the thermodynamic potential can be expressed where β = 1/T is the inverse of the temperature, and are the kinetic contribution from light quarks and s quarks, respectively. The quantities in the above are defined as λ ′ k = λ k −μ B 3 with λ k being the quasiparticle energy as detailed in Appendix B, and E ± s = E s ±μ s with E s = M 2 s + p 2 being the single s quark energy. In the present study without considering the color superconductivity, the color degree of freedom contributes a factor of N c = 3.

C. Gap equations
Using the quark propagator as detailed in Appendix C, the expressions of the chiral condensates and the netquark densities for u, d, and s quarks in terms of the phase-space distribution function can be written as The net baryon density ρ B and the isospin density ρ I can be calculated from The expression of the pion condensate π is In the above expressions, is the Fermi-Dirac distribution, the g functions have the form of and they satisfy the following relations Eqs. (18)-(23) can be obtained equivalently from with q = u, d, s being the quark flavor, leading to the relations

D. Equation of state
The energy density can be obtained from the thermodynamic potential through the thermodynamical relation where ε 0 is to ensure that the energy density is zero in vacuum. The pressure of the quark matter is where Ω 0 = ε 0 is the thermodynamic potential in vacuum, to ensure that the pressure is zero in vacuum.

III. RESULTS AND DISCUSSIONS
Based on the theoretical framework of the 3-flavor NJL model, we discuss the behavior of the pion condensate in both baryon-free and baryon-rich quark matter as well as the corresponding 3-dimensional phase diagram. We neglect the vector-isoscalar interaction G V = 0, and the s quark chemical potential is set to be µ s = 0, throughout the study. Reduced pion and chiral condensate (π/2σ0 and σu/σ0) (a), reduced isospin density ρI /m 3 π (b), and reduced energy density ε/εSB (c) as a function of the reduced isospin chemical potential µI /mπ in cold (T = 0) and baryon-free (µB = 0) quark matter. Available results from lattice QCD calculations are compared in panels (b) and (c).

A. Pion condensate in baryon-free quark matter
We start from fitting the physical pion mass and the lattice results by adjusting the isovector coupling constants. By choosing R IS = −0.002, the pion mass is fitted to be m π = 140.9 MeV. Once µ I is larger than this value, the reduced pion condensate π/2σ 0 becomes nonzero in cold and baryon-free quark matter, with σ 0 being the light quark condensate in vacuum, as shown in Fig. 1(a). This behavior can be intuitively understood [40] from the expression of the pion condensate π [Eq. (26)] together with the relation Eqs. (30) and (14). It can be seen that π = 0 is always a solution of Eq. (26), while for µ I > m π there appear nonzero solutions of π. For very large µ I ∼ 2 GeV, the nonzero solutions of π disappear again. Also shown in Fig. 1(a) is the decrease of the reduced chiral condensate of u quarks σ u /σ 0 after the appearance of the pion condensate. This corresponds to a second-order phase transition, where the values of the chiral and pion condensate change continuously while their derivatives have a sudden jump as µ I increases. For cold and baryon-free quark matter, the densities are zero below the threshold isospin chemical potential, and become finite above µ I ∼ m π as a result of the nonzero pion condensate. By choosing R IV = 0.25, the reduced isospin density as a function of the reduced isospin chemical potential reproduces very well the result from LQCD calculations [41], as shown in Fig. 1(b). With the fitted R IS and R IV , which are used throughout this study, the reduced energy density ε/ε SB as a function of the reduced isospin chemical potential is shown in Fig. 1 being the energy density in the Stefan-Boltzmann limit.
Comparing with the LQCD results in Fig. 22 of Ref. [42], our result gives a similar peak position of µ I /m π ∼ 1.3, while the peak value of ε/ε SB decreases with increasing spatial extent L from the LQCD calculations [42]. Similar results but in hot and baryon-free quark matter at T = 50, 100, and 150 MeV are displayed in Fig. 2. Due to the diffuseness of the Fermi-Dirac distribution at finite temperatures in Eqs. (18) and (26), the reduced chiral condensate σ u /σ 0 of u quarks is smaller at small µ I , while the reduced pion condensate π/2σ 0 also decreases, compared to that at T = 0. It is also interesting to see that the threshold isospin chemical potential µ I increases with the increasing temperature. Since at finite temperatures the densities become nonzero even for small chemical potentials, the isospin density ρ I becomes finite and increases with the increasing temperature below the threshold isospin chemical potential. As a consequence, the energy density shows a similar behavior below the threshold isospin chemical potential. It is seen that there is still a peak around µ I = 1.3m π at T = 50 MeV, while such peak disappears at higher temperatures. It is worth noting that ε SB approaches µ I ∼ 0 in the power of µ 4 I . Although ε is zero at T = 0 below the threshold isospin chemical potential, it becomes finite at finite temperatures. This leads to the divergence behavior of ε/ε SB when µ I approaches 0 as shown in Fig. 2(c). The equation of state (EOS) of cold and baryon-free quark matter is displayed in Fig. 3. It is seen that the energy density and the pressure become finite at a much larger isospin chemical potential µ I ∼ 2M , with M being the Dirac mass of light quarks in vacuum, in the absence of the pion condensate, compared to the case with the pion condensate incorporated. This leads to a larger energy density and pressure for a given isospin chemical potential in the presence of the pion condensate. As for the EOS as a function of the isospin density ρ I , the energy density is reduced while the pressure is increased with the pion condensate, compared to the case without considering the pion condensate. This is due to the different relations between ρ I and µ I in the two cases. It can be seen from Fig. 3 the pion condensate will stiffen the P ∼ ε relation. It could thus be further speculated that this might affect the properties of compact stars from solving the Tolman-Oppenheimer-Volkoff equation. However, this generally does not happen [43,44], since it requires a large isospin chemical potential but a not too large baryon chemical potential reached at the same time, related to the discussions in the next subsection. and 5 the reduced pion and chiral condensate in cold (T = 0 MeV) and hot (T = 50) and baryon-rich quark matter. Our rigourous calculation reduces to those in other studies using the assumption M u = M d at µ B = 0, but some differences are expected to appear at large baryon chemical potentials, and results from the two cases are compared in this and next subsections. It is seen from Fig. 4(a) and Fig. 5(a) that the results for smaller µ B are not qualitatively different from those in baryonfree quark matter shown in Fig. 1 (a) and Fig. 2(a). As µ B increases, except that the pion condensate becomes zero again at a slightly smaller µ I , there appears a second nonzero solution (π 2 ) of the pion condensate from Eq. (26) between about µ I ∼ 10m π and µ I ∼ 14m π , as shown in Fig. 4(b) and Fig. 5(b). At an even higher µ B , the occurrence of π 2 is at an even lower µ I . On the other hand, there is a first-order phase transition of the pion condensate at large µ I for both π and π 2 , where the pion condensate changes suddenly from a finite value to zero, as shown by the vertical lines in Fig. 4(c) and Fig. 5(c).
In addition, as shown in the same panel, it is seen that the assumption of M u = M d overestimates significantly the threshold µ I for π 2 , and underestimates the threshold µ I for π, especially at T = 50 MeV compared with those at T = 0 MeV. If the baryon chemical potential is further increased, the occurrence of the pion condensate becomes a first-order phase transition as well, and this is displayed in Fig. 4(d) and Fig. 5(d). It is also seen in Fig. 5(d) that the region of π = 0 and π 2 = 0 is much smaller from the rigourous calculations of M u = M d compared to those by assuming M u = M d . In addition, the chiral condensate generally has a sudden increase (decrease) once the pion condensate has a sudden decrease (increase) with the increasing isospin chemical potential. Figure 5(d) shows that the QCD phase structure can be different at large µ B and higher temperatures from the rigourous calculations comparing with those by assuming equal constituent mass for u and d quarks. Although π and π 2 calculated based on Eq. (26) satisfy the condition ∂Ω/∂π = 0, they are not both stable. For the T , µ B , and 4 typical isospin chemical potential µ I chosen according to Fig. 4(c), we show in Fig. 6 the thermodynamic potential Ω as a function of the reduced pion condensate, after subtracting the contribution Ω 0 in vacuum. At µ I = 141 MeV, it is seen that a local minimum point of Ω is about to appear, leading to the occurrence of π. At µ I = 355 MeV, except for a local minimum point corresponding to π, a local maximum point of Ω is about to appear, leading to the occurrence of π 2 . At µ I = 800 MeV, there are both local maximum and local minimum points of Ω, showing the existence of both π and π 2 states. At µ I = 2437 MeV, both the local maximum and local minimum of Ω are about to disappears, and π and π 2 turn to 0 accordingly. Since the solution π 2 corresponds to a maximum thermodynamic potential, it is an unstable solution, and the system favors π rather than π 2 . It is also argued that the instability of π 2 could be cured by considering the free energy of a system with a fixed baryon density [45] or in a Fermi system with a finite-range momentum-dependent interaction [46]. To understand in more details the properties of π and π 2 , we display in Fig. 7 the dispersion relations and netquark momentum distributions for u and d quarks in the condition of Fig. 6(c) and also at T = 50 MeV. For the dispersion relation, we show one of the 4 solutions of λ ′ k = λ k −μ B 3 with k = 1, 2, 3, 4 and λ k being the quasiparticle energy as detailed in Appendix B. The net-quark momentum distribution is defined as so integrating n u(d) gives the net-quark density, i.e., It is seen from Fig. 7(a) that the quasiparticle energy for π 2 becomes negative when the momentum p is between about 200 and 500 MeV/c. For the other 3 quasiparticle energy solutions for π 2 , they do not change sign as a function of the momentum. The negative quasiparticle energy state corresponds to the so-called Sarma phase [47], where the quasiparticle excitation does not need additional energy. The Sarma phase breaks the pairing of u andd quarks, as can be see from Fig. 7(c) that n d and n u is approximately of a constant value 0 and 1 in the region of negative λ ′ , respectively, compared with the pairing states when n u (p) = −n d (p) is always satisfied. At T = 50 MeV, the dispersion relation of λ ′ is similar to that at T = 0 MeV, while there is no sudden jump in the net-quark momentum distribution, though n u (p) and n d (p) become asymmetric around the momentum region where λ ′ is negative. This behavior corresponds to the 'quasi' Sarma phase at finite temperatures [48,49]. Also shown are the dispersion relation of λ ′ for the π solution corresponding to the local minimum of the thermodynamic potential as shown in Fig. 6. In such case, it is seen that λ ′ does not change sign, and n u (p) = −n d (p) is always satisfied at T = 0 MeV and approximately satisfied at T = 50 MeV. We then move to the discussions of the 3-dimensional QCD phase diagram. Since the appearance of the pion condensate always leads to the decrease of the chiral condensate, here we only discuss the behavior of the pion condensate in the T − µ B − µ I space, where the chiral condensate generally has an opposite behavior. Figure 8 displays the phase diagrams in the T −µ B plane at different isospin chemical potentials, where we show 3 phases as discussed above. Phase I is the normal baryon-rich and isospin-asymmetric quark matter with π = 0. Phase II is the pion superfluid phase with π = 0. Phase III is the phase with both nonzero solutions of π and π 2 , with the latter corresponding to the existence of the Sarma phase as discussed above. It is seen that Phase I generally exists at larger T or larger µ B , while Phase II generally exists at smaller T and µ B . Phase III exists in the area between the solid line and the dash-dotted line. It is seen that the phase transition between Phase I and Phase II, in the absence of Phase III, is always a second-order one, with the phase boundary represented by the dashed lines, so are the phase transition between Phase II and Phase III, with the phase boundary represented by the dash-dotted lines. On the other hand, the phase transition between Phase I and Phase III is always a first-order one, with the phase boundary represented by the solid lines. The critical end point (CEP), which connects the boundaries of the first-order phase transition and the second-order phase transitions, moves to a higher temperature with µ I changing from 200 MeV to 400 MeV, and the increasing trend saturates above µ I = 400 MeV. The assumption of M u = M d leads to a different phase structure at larger µ B , resulting in a CEP at lower temperatures and larger baryon chemical potentials.  Figure 9 displays the phase diagrams in the T − µ I plane at different baryon chemical potentials. It is seen that the normal quark phase (Phase I) generally exists at very small or large isospin chemical potentials, or at high temperatures, while the area of the pion superfluid phase (Phase II) shrinks dramatically with the increasing baryon chemical potential. The phase transitions are always of second-order at smaller baryon chemical potentials, while the first-order phase transition becomes more and more dominate with the increasing baryon chemical potential. Phase III with π 2 = 0 doesn't exist at µ B = 0 (not shown here), but it gradually appears inside Phase II at small baryon chemical potentials, and its area becomes larger and dominate at large baryon chemical potentials. The difference between results from rigourous calculations (M u = M d ) and approximations (M u = M d ) on the phase diagram is seen in panels (b) and (c), mainly exists in the relative area of Phase II and Phase III as well as the position of the CEP.   Figs. 8 and 9. The area of Phase II shrinks with the increasing temperature. Also, the first-order phase transition and Phase III become less dominate at higher temperatures. The deviations of results using the approximation (M u = M d ) from rigourous calculations (M u = M d ) can now be quantitatively seen within µ B ∈ (500, 900) MeV. Figure 10 is of course consistent with Figs. 8 and 9, and they together give a whole picture of the 3-dimensional QCD phase diagram.

IV. CONCLUSION
With the scalar-isovector and vector-isovector coupling constants adjusted to fit the physical pion mass and the lattice QCD results in baryon-free quark matter, we have studied the 3-dimensional QCD phase diagram by considering the pion condensate based on the 3-flavor NJL model. We found that the pion condensate becomes less important at higher temperatures or larger baryon chemical potentials. Thus, although incorporating the pion condensate would stiffen the equation of state of strange quark matter, it generally does not affect the properties of compact star systems where large baryon chemical potentials are reached. Besides the normal solution, we observe the appearance of a second nonzero solution of the pion condensate with the increase of the isospin chemical potential in baryon-rich quark matter, while both solutions disappear at very large isospin chemical potentials. The normal solution corresponds to the local minimum of the thermodynamic potential and represents the pion superfluid phase, while the second solution corresponds to the local maximum of the thermodynamic potential and is related to the Sarma phase. The occurrence or the disappearance of the pion condensate is a second-order phase transition at higher temperatures or smaller baryon chemical potentials, while it becomes a first-order one at lower temperatures or larger baryon chemical potentials. The calculations by assuming equal constituent mass of u and d quarks may introduce large errors in the 3-dimensional QCD phase diagram within µ B ∈ (500, 900) MeV, and affect the extraction of the critical end point, compared with the rigourous calculations in the present study.
To further explore the QCD phase structure, one can incorporate the polyakov loop into the NJL model, and study the interplay among the chiral condensate, the pion condensate, and the polyakov loop [22]. The kaon condensate can be further incorporated by considering systems at large strangeness chemical potentials. In addition, although the pion condensate generally will not affect the equation of state of strange quark matter and thus properties of compact stars, it is of great interest to further incorporate the chiral imbalance [50], the color superconductivity [51], etc., and see their effects on the QCD phase diagram and compact star properties. Further detailed properties of the QCD phase structure, e.g., the Larkin-Ovchinnikov-Fudde-Ferrell phase, are also worth investigating, as shown in Refs. [45,52].

Appendix A: The Lagrangian from mean-field approximation
In the mean-field approximation, it is assumed that deviations due to fluctuations of all quantities A from their thermal average values A are small. Thus, the following relations can be introduced to linearize the Lagrangian with Γ = {1, γ 5 , γ µ , γ 5 γ µ }, and the angular bracket denoting the expectation value from the quantum statistical average. In our previous studies [37], we assumed ψ γ k ψ = ψ γ 5 τ ψ = ψ γ 5 λ a ψ = ψ γ 5 γ µ ψ = 0 due to the parity symmetry in a static quark matter, so the condensates ψ i ψ j with i = j vanishes since it is assumed that the flavor is conserved in the case of µ I < m π . In order to study the 3-dimensional QCD phase diagram, we consider systems at larger isospin chemical potentials where pion condensates may arise, i.e., . In such case, the nonzero expectation value of ūiγ 5 d or d iγ 5 u spontaneously break the U I (1) symmetry, corresponding to the Bose-Einstein condensation of charged pions. The phase θ ud represents the direction of the U I (1) symmetry breaking. Since the thermodynamic potential does not depend on θ ud but depends only on |π ± | 2 or |π| 2 , we can set them to be real values corresponding to θ ud = 0 without losing generality. The Kaon condensate, which could be important at large strangeness chemical potentials, is not considered in the present study.
In the mean-field approximation by using the relations of Eq. (A1), the scalar-isoscalar term can be expressed as where V S = G S σ 2 u + σ 2 d + σ 2 s + GS 2 π 2 is the scalarisoscalar condensate energy, with being the chiral condensates for u, d, and s quarks, respectively, and is the self-energy contributed from the scalar-isoscalar interaction. The scalar-isovector term in the mean-field approximation can be expressed as where V IS = G IS (σ u − σ d ) 2 + G IS π 2 is the scalarisovector condensate energy, and is the self-energy contributed from the scalar-isovector interaction. The KMT term in the mean-field approximation can be expressed as is the condensate energy contribution from the KMT interaction, and is the self-energy contributed from the KMT interaction. Considering only the flavor-singlet state, the vectorisoscalar term in the mean-field approximation can be expressed as is the vector-isoscalar condensate energy, with ρ q (q = u, d, s) being the netquark density, and represents the contribution of the effective chemical potential from the vector-isoscalar interaction, with ρ = ρ u +ρ d +ρ s being the total net-quark density. The vectorisovector term in the mean-field approximation can be expressed as 2 is the vector-isovector condensate energy, and represents the contribution of the effective chemical potential from the vector-isovector interaction.
Appendix B: Thermodynamic potential In this appendix, we obtain the thermodynamic potential from the partition function, which in the grand canonical ensemble is written as where β = T −1 is the inverse of the temperature, µ is the chemical potential, andĤ andN are the Hamiltonian operator and the quark number operator, respectively. The sum a dΨ a is carried out over all states. According to the finite-temperature field theory, the partition function in the mean-field approximation can be expressed in the form of the path integral with the real time t = x 0 converted to the imaginary time τ = it, and the functional integral DψDψ covering all quark species. In the above equation, the condensate energy independent of ψ andψ can be factorized in Z, and after applying the relation the partition function can be simplified as where V is the system volume, and the 4-momentum becomes p = (p 0 , p) = (iω n , p) with ω n = (2n + 1)πT being the Matsubara frequency for a fermi system. The thermodynamic potential of the quark system can be obtained from the partition function through With the form of the quark propagator as Eq. (13) and keeping in mind that p 0 = iω n , we can get the following relation after some algebras with ∆ = (G S + 2G IS − Kσ s ) π being the gap parameter, and E q = M 2 q + p 2 with q = u, d, s being the single-quark energy. Replacingμ u andμ d with the effective baryon and isospin chemical potentialμ B andμ I according to Eq. (16), Eq. (B6) becomes where λ k is the solution of the quartic equation for p 0 + µ B /3, with the coefficients with Ξ 1 = c 2 + 12 e, Ξ 2 = 2 c 3 + 27 d 2 − 72 ce, These roots satisfy the following relations In the special case of M u = M d , the four roots become with E = M 2 u + p 2 = M 2 d + p 2 . Since the sign is degenerate, there are actually only two solutions. In the baryon-free system, the general case always reduces to the special case of M u = M d .
Using the summation formula of the Matsubara frequencies and combining Eqs. (B5) and (B7), we can get the expression of the thermodynamic potential as Eq. (17).

Appendix C: Quark condensate and quark density
In this appendix we obtain the expressions of the chiral condensates, the net-quark densities, and the pion condensate in terms of the phase-space distribution function from the quark propagator. In the absence of the pion condensate, u, d, and s quarks are decoupled, and the quark propagator can be written as In the presence of the pion condensate, u and d quarks are mixed, and the off-diagonal terms appear. The quark propagator is then In the above, the diagonal terms of the quark propagator in the absence of the pion condensate are and E ± q = M 2 q + p 2 ±μ q , for q = u, d, s. The detailed expressions of each quark propagator are S uu (p) = 4 k=1 g uu (λ ′ k ) With the following relations [18]