Reentrant pion superfluidity and cosmic trajectories within PNJL model

In this work, we self-consistently explore the possibility of charged pion superfluidity and cosmic trajectories in early Universe under the framework of Polyakov-Nambu--Jona-Lasinio model. By taking the badly constrained lepton flavor asymmetries $l_{\rm e}$ and $l_\mu$ as free parameters, the upper boundaries of pion superfluidity phase are consistently found to be around the pseudocritical temperature at zero chemical potentials. So the results greatly support the choice of $T=0.16~{\rm GeV}$ as the upper boundary of pion superfluidity in the previous lattice QCD study. Take $l_{\rm e}+l_\mu=-0.2$ as an example, we demonstrate the features of pion condensation and the associated cosmic trajectories with the evolution of early Universe. While the trajectory of electric chemical potential reacts strongly at both the lower and upper boundaries of reentrant pion superfluidity, the trajectories of other chemical potentials only respond strongly at the upper boundary.


I. INTRODUCTION
As we know, the field of high energy nuclear physics (HENP) initiated with the search of quark-gluon plasma (QGP) phase [1] for quantum chromodynamics (QCD) matter, and the QGP phase was expected to be realized through relativistic heavy ion collisions (HICs) [2,3]. Up to now, the existence of QGP at high temperature becomes a consensus in HENP -a direct evidence is the observation of number-of-constituent quark scaling of the elliptic flows for mesons and baryons in large center of mass HICs [4]. Of course, other properties of QGP have also been well explored and one remarkable discovery is that the QGP is a nearly perfect liquid [5,6]. Nevertheless, a much more sophisticated and challenging mission is to depict the first QCD phase diagram in T − µ B plane with the help of HICs. At early time, people didn't find it a real phase transition from QGP to hadron phase at small baryon number density n B according to either lattice QCD simulations [7,8] or experimental detections [9,10]. Recently, the STAR group is carrying out lower energy collisions in their BES II experiments with the hope of catching the critical end point by increasing n B [11].
Besides, the T − µ I phase diagram has also been extensively studied and charged pion (π ± ) superfluidity was expected theoretically in the region with not too large T and large isospin chemical potential µ I (> m v π ) [12][13][14][15]. Specifically, the transition between chiral symmetry breaking or restoration phase and π ± superfluidity was consistently found to be of second-order in chiral perturbation theory [12], lattice QCD [13,14] and effective models such as Nambu-Jona-Lasinio (NJL) model [15]. And with the increasing of µ I , the Bose-Einstein condensation of π ± was found to smoothly crossover to the Bardeen-Cooper-Schieffer phase [16] which then possibly becomes the quarksonic matter [17]. However, it's a pity that the systems in nature with large µ I , neutron stars, are also large µ B matters, and the constraint of electric neutrality eventually disfavors π ± superfluidity in cold neutron stars [18,19].
The hope of finding π ± superfluidity in nature reburned in the explorations of proto-neutron stars [20] and early Universe [21,22], where the temperature and lepton flavor densities can be much larger. Especially, large lepton flavor densities would also help to stabilize π ± mesons against weak decays thanks to the Pauli blocking effect to the final state [20]. And it is interesting that the QCD phase would impact primordial gravitational wave (GW) and generation rate of black holes quite well in early Universe [22]. According to the Big Bang theory, temperature drops to ∼ 200 MeV in 10 −6 s of the Big Bang and the early Universe enters the QCD epoch where strong interaction dominates particle scattering. In the QCD epoch, the baryon and lepton number U (1) symmetries were already violated, and the tracing back of the observations in present Universe constrains the corresponding number densities as n B /s = 8.6 * 10 −11 [23] and |n l /s| < 0.012 [24] with s the entropy density. Under these and electric neutrality constraints, the cosmic trajectories of several chemical potentials were explored by combining hadron resonance gas model, lattice QCD and free quark gas model [25]. Later, the possibility of π ± superfluidity was realized at large lepton flavor asymmetries [21,22] and the lower phase boundary was reasonably depicted by utilizing the criteria |µ Q | > m π ± (µ Q , T ) [22], where µ Q plays the role of µ I for charged pions.
In Ref. [22], the temperature effect on m π ± (µ Q , T ), the order parameter of their effective mass model, was only taken into account through the ideal gas part. This might be the reason why the lattice QCD inspired model couldn't predict the upper phase boundary. And it is also a pity that they evaluated the phase boundary without demonstrating the evolutions of the true order parameter: charged pion condensates. The advantage of NJL model is that chiral symmetry breaking and restoration, pion superfluidity, and the corresponding pion masses can be self-consistently studied by solving the gap equations and the zero points of pion propagators, all of which can be derived analytically, see Ref. [15]. So we intend to recheck the phase boundary of π ± superfluidity in the Polyakov loop extended NJL model (PNJL model) [29] and show how the order parameters and cosmic trajectories evolve across the π ± superfluidity phases. In principle, the PNJL model is able to mimic QCD more realistically by counting the deconfinement effect to quarks and gluon contributions to the total entropy [29].
The paper is organized as follows. In Sec.II, we develop the overall formalism for the study of the QCD and quantum electroweak dynamics (QEWD) matter in the QCD epoch. For the QCD sector, the two-and three-flavor PNJL models will be laid out explicitly in Sec.II A and Sec.II B, respectively. The QEWD sector will be approximated as free gases in Sec.II C. Then, we present numerical results in Sec.III and finally summarize in Sec.IV.

II. THE OVERALL FORMALISM
In the QCD epoch of early Universe, both QCD and QEWD sectors are relevant: the elementary degrees of freedom are quarks and gluons in the QCD sector and leptons and photons in the QEWD sector. To study the QCD sector self-consistently, we adopt the chiral effective Polyakov-Nambu-Jona-Lasinio (PNJL) model, where quarks contribute through the NJL model part and the contributions of gluons are given in terms of Polyakov loop (PL) according to the lattice QCD simulations [26][27][28][29]. For the QEWD sector, the coupling constants are usually small, so we are satisfied to utilize free gas approximation for the involved particles. The following sections are devoted to developing detailed formalisms for both the QCD matter with two or three flavors and the free QEWD matter.

A. The QCD sector with two flavors
The Lagrangian density of two-flavor PNJL model with electric charge chemical potential µ Q and baryon chemical potential µ B can be given as [29][30][31] in Euclidean space. In the NJL model part, ψ = (u, d) T represents the two-flavor quark field and A 4 = A 4c T c /2 is the non-Abelian gauge field with T c the Gell-Mann matrices in color space; m 0 ≡ m 0l ½ 2 is the current mass matrix, the charge number matrix is and τ are Pauli matrices in flavor space. The pure gluon potential is given as a function of the Polyakov loop and its complex conjugate L * by fitting to the lattice QCD data, that is, withT ≡ T /T 0 and T 0 = 0.27 GeV [29].
To obtain the analytic form of the basic thermodynamic potential, we take Hubbard-Stratonovich transformation with the help of the auxiliary fields σ = −2Gψψ and π = −2Gψiγ 5 τ ψ [30] and the Lagrangian becomes For later convenience, we alternatively represent it as in the forms of physical particles: π 0 = π 3 and π ± = 1 √ is the lowering/raising operator in flavor space. To our present interests, it is enough to assume the expectation values of the auxiliary fields to be Then, in mean field approximation, the thermodynamic potential can be given in energy-momentum space as with the trace Tr over the energy-momentum, spinor, flavor, and color spaces. To derive the explicit form of the trace term, we need to solve the quark dispersions from the zero points of their inverse propagator in Minkowski space, that is, from Finally, in the saddle point approximation [29]: L = L * , the thermodynamic potential can be given directly as [32] where three-momentum cutoff Λ is adopted to regularize the divergent vacuum term. Armed with the equation of state, the coupled gap equations follow the minimal conditions, where we've defined two dimensionless auxiliary functions for future use: Furthermore, the entropy, electric charge number, and baryon number densities can also be derived analytically according to the thermodynamic relations as B. The QCD sector with three flavors In order to explore the properties of QCD matter more realistically, we adopt the three-flavor PNJL model where more low-lying mesons are involved and the QCD U A (1) anomaly has been properly taken into account through the 't Hooft term. In saddle point approximation, the corresponding Lagrangian density can be given by [29][30][31]: where ψ = (u, d, s) T is now the three-flavor quark field. Similar to the two-flavor case, the current mass and electric charge number matrices of quarks are respectively m 0 ≡ diag(m 0u , m 0d , m 0s ), the interaction index λ 0 = 2/3 ½ 3 and λ i (i = 1, . . . , 8) are Gell-Mann matrices in flavor space. For later use, the 't Hooft term L tH ≡ −K t=± DetψΓ t ψ can be represented as with the interaction vertices Γ ± = ½ 4 ± γ 5 for right-and left-handed channels, respectively. Here, one should note the Einstein summation convention for the flavor indices i, j, k, m, n and the correspondences between 1, 2, 3 and u, d, s.
To our main concerns, we choose the following scalar and charged pseudoscalar condensates to be nonzero: For brevity, we set ∆ π = ∆ * π without loss of generality in the following. To facilitate the study, we'd like first to reduce L tH to an effective form with four-fermion interactions at most. By applying the Hartree approximation to contract a pair of quark and antiquark in each sixfermion interaction term [30], we immediately find where the second term in the brace is induced by π ± condensations. Armed with the reduced Lagrangian density: the left calculations can just follow the two-flavor case in principle. By contracting quark and antiquark pairs once more in the interaction terms of Eq.(21), we find the quark bilinear form as where the scalar and pseudoscalar masses are respectively with i = j = k. The G and K dependent terms in Eq. (23) are from the U A (1) symmetric and anomalous interactions, respectively. According to Eq.(22), s quark decouples from u, d quarks, so the gap equation for σ s can be simply given by [30]: However, the u and d light quarks couple with each other through the non-diagonal pseudoscalar mass Π. Since µ B is usually small in early Universe, we can simply set m 0u = m 0d ≡ m 0l , σ u = σ d ≡ σ l in order to further carry out analytic derivations. Then, by following a similar procedure as the previous section, the explicit thermodynamic potential can be worked out for the bilinear terms as with the energy functions defined by Eventually, the coupled gap equations follow directly from the definitions of condensates: and the minimal condition ∂ L [V (L, L) + Ω bl ] = 0 as [33] Note that ∆ π = 0 is a trivial solution of Eq.(30), so ∆ π or Π is still a true order parameter for I 3 isospin symmetry [12] in three-flavor case. The total self-consistent thermodynamic potential can be found to be by utilizing the definitions of condensates and their relations to scalar and pseudoscalar masses, refer to Eqs. (27) and (23). And the entropy, electric charge number and baryon number densities can be given according to the thermodynamic relations as C. The QEWD sector: free gases In free gas approximation, the thermodynamic potentials for the QEWD sector can be easily given by [34] where the degeneracy is one for neutrinos and antineutrinos due to their definite chiralities. Then, the corresponding entropy, electric charge number and lepton flavor number densities can be derived directly as Now, collecting contributions from both QEWD and QCD sectors, the total entropy, electric charge number and lepton number densities are respectively s = s γ +s l +s 2f/3f , n Q = n l Q +n by following the conventions. According to the introduction, l is not so well constrained as b from the observations, but the standard picture well predicts that l = −51/28 b [35]. Furthermore, due to neutrino oscillations at late stage of the Universe, the lepton flavor densities l i are not well constrained at the QCD epoch, thus we will take l e and l µ as free variables in the following.

III. NUMERICAL RESULTS
To carry out numerical calculations, we get the muon mass from the Particle Data Group as m µ = 113 MeV, simply set the electron mass m e = 0, and suppress the contribution of heavy τ leptons for the QEWD sector. The model parameters are fixed for the QCD sector as the following [36,37]  The T −(l e +l µ ) phase diagrams of two-and three-flavor PNJL models are illuminated in Fig.1. As we can see, the ratio l e /l µ and the effect of strange quarks are all not so important for determining the phase boundaries of pion superfluidity, especially the upper ones. And the calculations with the standard lepton asymmetry in the upper panel indicate that the uncertainty of l plays a negligible role in the exploration of phase boundary. Compared to the threshold lepton flavor asymmetry |l e + l µ | ∼ 0.1 in the extrapolated lattice QCD calculations [22], the values are consistently |l e + l µ | ∼ 0.09 at their top temperature T pc = 0.16 GeV in our evaluations. So recalling the effectivenesses of PNJL model and their criterion for the phase boundary, the agreement is remarkable! In advance, we obtain the upper boundaries of pion superfluidity to be consistently T ∼ 0.21 GeV which is much larger than T pc , a well-known drawback of PNJL model [32]. Nevertheless, the upper boundaries are almost the corresponding pseudo-critical temperatures at zero chemical potentials, which then supports the setting of the upper boundary around T pc in Ref. [22]. By the way, here the threshold lepton flavor asymmetries are 0.08 and 0.075 in the two-and three-flavor cases, respectively. Now, we take the more realistic three-flavor PNJL model for example to show the features of cosmic trajec-tories at l e + l µ = −0.2, where the early Universe could evolve through the pion superfluidity phase. We compare two cases: l e = 0 and l e = l µ , and demonstrate the order parameters and chemical potentials in Fig.2 and Fig.3, respectively. As we can see, the effect of the ratio l e /l µ is only important on the cosmic trajectories of the directly related quantities: µ e and µ µ , and the results almost overlap with each other for other quantities. According to the lower panel of Fig.2, the pion condensate Π shows a reentrant feature with T : although the decreasing at higher temperature can be easily understood as isospin symmetry restoration, the increasing at lower temperature is not trivial. Actually, the latter is due to the enhancement of |µ Q | with T under the constraint of electric neutrality, see the upper panel in Fig.3. In Fig.2, it is also interesting to notice the consistency between the critical temperatures of I 3 isospin symmetry restoration and deconfinement transition at finite charge chemical potential [32]. Following the monotonous feature of |µ Q | in the upper panel of Fig.3, we find the criterion |µ Q | = m v π to be well satisfied at the lower critical temperature T l ; but that is no longer useful for the exploration of upper critical temperature T u , since the effective π ± mass increases with T . Since the |µ Q | varies from m v π to ∼ 3 m v π within the pion superfluidity phase, we can well recognize the BCS-BEC crossover therein [16] as the early Universe cooled down. As expected, the baryon chemical potential is small except for very low temperature, which justifies our assumptions in the QCD sector: m u = m d and L = L * . Furthermore, the opposite signs between µ B and µ Q and the abrupt jump of µ B at low temperature both qualitatively fit the findings in the extrapolated lattice QCD study [25]. We note that the exist of pion superfluidity at T l only leaves visible sign in the behavior of µ Q but the entrance at T u gives rise to important signs in all the chemical potentials. Compared to the signs S(µ e ) = S(µ µ ) = S(µ τ ) for the choice l e = l µ = l τ in Ref. [25], we find S(µ e ) = S(µ µ ) = −S(µ τ ) due to the choice |l e + l µ | ≫ |l|.

IV. SUMMARY
In this work, the possibility of pion superfluidity and the corresponding cosmic trajectories are selfconsistently explored by varying the lepton flavor asymmetries within the PNJL model. The effects of strange quarks, total lepton asymmetry and the ratio l e /l µ are all found to be mild on the phase boundary of pion superfluidity. Following the previous study in Ref. [22], the phase boundary is constrained from both lower and upper sides in our study. At T cp , the lepton flavor asymmetry |l e +l µ | is 0.09 in our work, quite consistent with the threshold value 0.1 obtained in Ref. [22]. However, with the pseudocritical temperature in PNJL model much larger than that from lattice QCD, we find that the threshold values shift to 0.08 and 0.075 for two-and three-flavor cases, respectively.
According to the three-flavor example, the pion condensation shows a non-monotonous or reentrant feature as should be for the existences of both upper and lower second-order phase boundaries. While the sign of lower critical temperature is only visible in µ Q , the signs of upper one show up in all the chemical potentials. So in principle, the phase transitions to and from pion superfluidity during the evolution of early Universe can be identified through the non-analytic features in the cosmic trajectories. Moreover, the critical |µ Q | at T u is found to be around 3 m v π , which is so large that it explains why the extrapolated lattice QCD study was unable to fix T u at all [22]. In Ref. [22], the equation of state of QCD+QEWD matter with different |l e + l µ | was adopted to study the relic density of primordial gravitational wave; then inversely the observations of GW would help to constrain |l e + l µ | and thus indirectly indicate whether pion superfluidity had happened or not in the QCD epoch. For a second order phase transition, we don't expect it to leave direct relics on GW, such as the pion superfluidity in this work; but for a first order one, the transition itself will give rise to a specific GW spectrum [38][39][40]. Actually, with strong magnetic field presented in early Universe [41][42][43], the transitions relevant to pion superfluid, which is also superconductor, could be shifted from second order to first order. We will explore such interesting situation in more detail in our coming work. Hopefully, the advanced GW detectors, such as LIGO, SKA, LISA, and Tianqin, could help to capture the signals in future.