Reentrant pion superfluidity and cosmic trajectories within PNJL model

Gaoqing Cao¹ Lianyi He² Pengming Zhang¹ ¹School of Physics and Astronomy, Sun Yat-sen University, Zhuhai 519088, China
²Department of Physics, Tsinghua University, Beijing 100084, China

Abstract

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.

pacs:

11.30.Qc, 05.30.Fk, 11.30.Hv, 12.20.Ds

I Introduction

As we know, the field of high energy nuclear physics (HENP) initiated with the search of quark-gluon plasma (QGP) phase Shuryak:1980tp for quantum chromodynamics (QCD) matter, and the QGP phase was expected to be realized through relativistic heavy ion collisions (HICs) Lee:1982qw ; Heinz:2000bk . 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 Adams:2005zg . Of course, other properties of QGP have also been well explored and one remarkable discovery is that the QGP is a nearly perfect liquid Kovtun:2004de ; Kolb2004 . Nevertheless, a much more sophisticated and challenging mission is to depict the first QCD phase diagram in $T-\mu_{\rm 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_{\rm B}$ according to either lattice QCD simulations Aoki:2006we ; Bhattacharya:2014ara or experimental detections Floris:2014pta ; Adamczyk:2017iwn . 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_{\rm B}$ Luo:2017faz .

Besides, the $T-\mu_{\rm I}$ phase diagram has also been extensively studied and charged pion ( $\pi^{\pm}$ ) superfluidity was expected theoretically in the region with not too large $T$ and large isospin chemical potential $\mu_{\rm I}~{}(>m_{\pi}^{\rm v})$ Son:2000xc ; Kogut:2002zg ; Brandt:2017oyy ; He:2005nk . Specifically, the transition between chiral symmetry breaking or restoration phase and $\pi^{\pm}$ superfluidity was consistently found to be of second-order in chiral perturbation theory Son:2000xc , lattice QCD Kogut:2002zg ; Brandt:2017oyy and effective models such as Nambu–Jona-Lasinio (NJL) model He:2005nk . And with the increasing of $\mu_{\rm I}$ , the Bose-Einstein condensation of $\pi^{\pm}$ was found to smoothly crossover to the Bardeen-Cooper-Schieffer phase Sun:2007fc which then possibly becomes the quarksonic matter Cao:2016ats . However, it’s a pity that the systems in nature with large $\mu_{\rm I}$ , neutron stars, are also large $\mu_{\rm B}$ matters, and the constraint of electric neutrality eventually disfavors $\pi^{\pm}$ superfluidity in cold neutron stars Andersen:2007qv ; Abuki:2008tx .

The hope of finding $\pi^{\pm}$ superfluidity in nature reburned in the explorations of proto-neutron stars Abuki:2009hx and early Universe Middeldorf-Wygas:2020glx ; Vovchenko:2020crk , where the temperature and lepton flavor densities can be much larger. Especially, large lepton flavor densities would also help to stabilize $\pi^{\pm}$ mesons against weak decays thanks to the Pauli blocking effect to the final state Abuki:2009hx . 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 Vovchenko:2020crk . According to the Big Bang theory, temperature drops to $\sim 200\,{\rm MeV}$ in $10^{-6}\,{\rm 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_{\rm B}/s=8.6*10^{-11}$ Ade2016 and $|n_{\rm l}/s|<0.012$ Oldengott2017 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 Wygas:2018otj . Later, the possibility of $\pi^{\pm}$ superfluidity was realized at large lepton flavor asymmetries Middeldorf-Wygas:2020glx ; Vovchenko:2020crk and the lower phase boundary was reasonably depicted by utilizing the criteria $|\mu_{\rm Q}|>m_{\pi^{\pm}}(\mu_{\rm Q},T)$ Vovchenko:2020crk , where $\mu_{\rm Q}$ plays the role of $\mu_{\rm I}$ for charged pions.

In Ref. Vovchenko:2020crk , the temperature effect on $m_{\pi^{\pm}}(\mu_{\rm 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. He:2005nk . So we intend to recheck the phase boundary of $\pi^{\pm}$ superfluidity in the Polyakov loop extended NJL model (PNJL model) Fukushima:2017csk and show how the order parameters and cosmic trajectories evolve across the $\pi^{\pm}$ 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 Fukushima:2017csk .

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.1 and Sec.II.2, respectively. The QEWD sector will be approximated as free gases in Sec.II.3. 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 Meisinger:1995ih ; Fukushima:2003fw ; Ratti:2005jh ; Fukushima:2017csk . 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.

II.1 The QCD sector with two flavors

The Lagrangian density of two-flavor PNJL model with electric charge chemical potential $\mu_{\rm Q}$ and baryon chemical potential $\mu_{\rm B}$ can be given as Fukushima:2017csk ; Klevansky:1992qe ; Hatsuda:1994pi

	$\displaystyle{\cal L}$	$\displaystyle=$	$\displaystyle\bar{\psi}\left[i\not{\partial}-i\gamma^{4}\left(ig{\cal A}^{4}+Q_{\rm q}\mu_{\rm Q}+{\mu_{\rm B}\over 3}\right)-m_{0}\right]\psi$		(1)
			$\displaystyle+G\left[\left(\bar{\psi}\psi\right)^{2}+\left(\bar{\psi}i\gamma_{5}\bm{\tau}\psi\right)^{2}\right]-V(L,L^{*})$		(1)

in Euclidean space. In the NJL model part, $\psi=(u,d)^{T}$ represents the two-flavor quark field and ${\cal A}^{4}=A^{\rm 4c}T^{\rm c}/2$ is the non-Abelian gauge field with $T^{\rm c}$ the Gell-Mann matrices in color space; $m_{0}\equiv m_{\rm 0l}\mathbbm{1}_{2}$ is the current mass matrix, the charge number matrix is

\displaystyle Q_{\rm q}\equiv{\rm diag}(q_{\rm u},q_{\rm d})={1\over 3}{\rm diag}(2,-1),

(2)

and ${\bm{\tau}}$ are Pauli matrices in flavor space. The pure gluon potential is given as a function of the Polyakov loop

L={1\over N_{\rm c}}{\rm tr}\,e^{ig\int{\rm d}x_{4}{\cal A}^{4}}

and its complex conjugate $L^{*}$ by fitting to the lattice QCD data, that is,

	$\displaystyle{V(L,L^{*})\over T^{4}}$	$\displaystyle=$	$\displaystyle-{1\over 2}\left(3.51-{2.47\over\tilde{T}}+{15.2\over\tilde{T}^{2}}\right)\|L\|^{2}-{1.75\over\tilde{T}^{3}}$		(3)
			$\displaystyle\times\ln\left[1-6\|L\|^{2}+4(L^{3}+{L^{*}}^{3})-3\|L\|^{4}\right]$		(3)

with $\tilde{T}\equiv T/T_{0}$ and $T_{0}=0.27\,{\rm GeV}$ Fukushima:2017csk .

To obtain the analytic form of the basic thermodynamic potential, we take Hubbard-Stratonovich transformation with the help of the auxiliary fields $\sigma=-2{G}\bar{\psi}\psi$ and ${\bm{\pi}}=-2{G}\bar{\psi}i\gamma^{5}{\bm{\tau}}\psi$ Klevansky:1992qe and the Lagrangian becomes

	$\displaystyle{\cal L}$	$\displaystyle=$	$\displaystyle\bar{\psi}\!\left[i{\not{\partial}}\!-\!i\gamma^{4}\!\left(\!ig{\cal A}^{4}\!+\!Q_{\rm q}\mu_{\rm Q}\!+\!{\mu_{\rm B}\over 3}\!\right)\!-\!i\gamma^{5}\bm{\tau}\cdot\bm{\pi}\!-\!\sigma\!-\!m_{0}\right]\!\psi$		(4)
			$\displaystyle-{\sigma^{2}+\bm{\pi}^{2}\over 4G}-V(L,L^{*}).$		(4)

For later convenience, we alternatively represent it as

	$\displaystyle{\cal L}\!$	$\displaystyle=$	$\displaystyle\!\bar{\psi}\!\Big{[}i{\not{\partial}}\!-\!i\gamma^{4}\!\!\left(ig{\cal A}^{4}+Q_{\rm q}\mu_{\rm Q}\!+\!{\mu_{\rm B}\over 3}\right)\!\!-\!i\gamma^{5}\!\!\left({\tau_{3}}{\pi^{0}}\!+\!{\tau_{\pm}}{\pi^{\pm}}\right)$		(5)
			$\displaystyle\!-\sigma-m_{0}\Big{]}\!\psi-{\sigma^{2}+{\left({\pi^{0}}\right)^{2}}+2{\pi^{+}}{\pi^{-}}\over 4G}-V(L,L^{*})$		(5)

in the forms of physical particles: $\pi^{0}=\pi_{3}$ and $\pi^{\pm}={1\over\sqrt{2}}(\pi_{1}\mp i\pi_{2})$ , where $\tau_{\mp}={1\over\sqrt{2}}(\tau_{1}\mp i\tau_{2})$ 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

\langle\sigma\rangle=m-m_{0},\ \langle\pi^{0}\rangle=0,\ \langle\pi^{\pm}\rangle=\Pi/\sqrt{2}.

Then, in mean field approximation, the thermodynamic potential can be given in energy-momentum space as

	$\displaystyle\Omega_{\rm 2f}$	$\displaystyle=$	$\displaystyle-{\rm Tr}\ln\left[{\not{k}}\!-\!m\!-\!i\gamma^{4}\!\!\left(ig{\cal A}^{4}\!+\!Q_{\rm q}\mu_{\rm Q}\!+\!{\mu_{\rm B}\over 3}\right)\!-\!i\gamma^{5}\tau_{1}\Pi\right]$		(6)
			$\displaystyle+{(m-m_{0})^{2}+\Pi^{2}\over 4G}+V(L,L^{*})$		(6)

with the trace ${\rm 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

\displaystyle{\rm Det}\!\left[{\not{k}}\!-\!m\!+\!\gamma^{0}\!\!\left(ig{\cal A}^{4}\!+\!Q_{\rm q}\mu_{\rm Q}\!+\!{\mu_{\rm B}\over 3}\right)\!-\!i\gamma^{5}\tau_{1}\Pi\right]=0.

(7)

We get $k_{0}=E^{\rm t}(k)\pm\left({\mu_{\rm Q}+2\mu_{\rm B}\over 6}+i\langle g{\cal A}^{4}\rangle\right)$ with He:2005nk

E^{\pm}(k)=\sqrt{\left(\epsilon(k)\pm{\mu_{\rm Q}\over 2}\right)^{2}+\Pi^{2}},\ \epsilon(k)=\sqrt{k^{2}+m^{2}}.

Finally, in the saddle point approximation Fukushima:2017csk : $L=L^{*}$ , the thermodynamic potential can be given directly as Xiong:2009zz

	$\displaystyle\Omega_{\rm 2f}$	$\displaystyle=$	$\displaystyle V(L,L)+{(m-m_{0})^{2}+\Pi^{2}\over 4G}-2N_{c}\int^{\Lambda}{{\rm d}^{3}k\over(2\pi)^{3}}\sum_{t=\pm}E^{\rm t}(k)-2T\!\!\int\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\!\sum_{t,u=\pm}Fl\left(L,u,E^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\right),$
			$\displaystyle Fl(L,u,x,y)=\log\left[1+3L\,e^{-{1\over T}\left(x-u\,y\right)}+3L\,e^{-{2\over T}\left(x-u\,y\right)}+e^{-{3\over T}\left(x-u\,y\right)}\right],$

where three-momentum cutoff $\Lambda$ is adopted to regularize the divergent vacuum term.

Armed with the equation of state, the coupled gap equations follow the minimal conditions, $\partial_{\rm m}\Omega_{\rm 2f}=\partial_{\rm\Pi}\Omega_{\rm 2f}=\partial_{\rm L}\Omega_{\rm 2f}=0$ , as

$\displaystyle 0$	$\displaystyle=$	$\displaystyle{m-m_{0}\over 2G}-2N_{c}\int^{\Lambda}{{\rm d}^{3}k\over(2\pi)^{3}}\sum_{t=\pm}{m\over\epsilon(k)}{\epsilon(k)\!+\!t{\mu_{\rm Q}\over 2}\over E^{\rm t}(k)}+6\int{{\rm d}^{3}k\over(2\pi)^{3}}\sum_{t,u=\pm}{m\over\epsilon(k)}{\epsilon(k)\!+\!t{\mu_{\rm Q}\over 2}\over E^{\rm t}(k)}dV_{1}\left(L,u,E^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\right),$	(9)
$\displaystyle 0$	$\displaystyle=$	$\displaystyle{\Pi\over 2G}-2N_{c}\int^{\Lambda}{{\rm d}^{3}k\over(2\pi)^{3}}\sum_{t=\pm}{\Pi\over E^{\rm t}(k)}+6\int{{\rm d}^{3}k\over(2\pi)^{3}}\sum_{t,u=\pm}{\Pi\over E^{\rm t}(k)}dV_{1}\left(L,u,E^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\right),$	(10)
$\displaystyle 0$	$\displaystyle=$	$\displaystyle T^{4}\left[-\left(3.51\!-{2.47\over\tilde{T}}\!+\!{15.2\over\tilde{T}^{2}}\right)L\!+\!{1.75\over\tilde{T}^{3}}{12L(1-L)^{2}\over 1\!-\!6L^{2}\!+\!8L^{3}\!-\!3L^{4}}\right]-6T\!\int\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\!\sum_{t,u=\pm}dV_{2}\left(L,u,E^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\right),$	(11)

where we’ve defined two dimensionless auxiliary functions for future use:

	$\displaystyle dV_{1}(L,u,x,y)$	$\displaystyle=$	$\displaystyle{L\,e^{-{1\over T}\left(x-u\,y\right)}+2L\,e^{-{2\over T}\left(x-u\,y\right)}+e^{-{3\over T}\left(x-u\,y\right)}\over 1+3L\,e^{-{1\over T}\left(x-u\,y\right)}+3L\,e^{-{2\over T}\left(x-u\,y\right)}+e^{-{3\over T}\left(x-u\,y\right)}},$		(12)
	$\displaystyle dV_{2}(L,u,x,y)$	$\displaystyle=$	$\displaystyle{e^{-{1\over T}\left(x-u\,y\right)}+e^{-{2\over T}\left(x-u\,y\right)}\over 1+3L\,e^{-{1\over T}\left(x-u\,y\right)}+3L\,e^{-{2\over T}\left(x-u\,y\right)}+e^{-{3\over T}\left(x-u\,y\right)}}.$		(13)

Furthermore, the entropy, electric charge number, and baryon number densities can also be derived analytically according to the thermodynamic relations as

$\displaystyle\!\!\!\!\!\!s_{\rm 2f}$	$\displaystyle=$	$\displaystyle-{\partial\Omega_{\rm 2f}\over\partial T}=2\!\!\int\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\!\sum_{t,u=\pm}\left[Fl\left(\!L,u,E^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\!\right)\!+\!{3\left(\!E^{\rm t}(k)\!-\!u\,{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\!\right)\over T}dV_{1}\left(\!L,u,E^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\!\right)\right]$	(14)
		$\displaystyle+T^{3}\left\{{1\over 2}\left(4\times 3.51-3\times{2.47\over\tilde{T}}+2\times{15.2\over\tilde{T}^{2}}\right)L^{2}+{1.75\over\tilde{T}^{3}}\ln\left[1-6L^{2}+8L^{3}-3L^{4}\right]\right\},$	(14)
$\displaystyle\!\!\!\!\!\!n_{\rm Q}^{\rm 2f}$	$\displaystyle=$	$\displaystyle-{\partial\Omega_{\rm 2f}\over\partial\mu_{\rm Q}}=N_{c}\!\!\int^{\Lambda}\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\!\sum_{t=\pm}t{\epsilon(k)\!+\!t{\mu_{\rm Q}\over 2}\over E^{\rm t}(k)}-\!3\!\!\int\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\sum_{t,u=\pm}t{\epsilon(k)\!+\!t{\mu_{\rm Q}\over 2}\over E^{\rm t}(k)}dV_{1}\left(\!L,u,E^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\!\right)$	(15)
		$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ +\int\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\!\sum_{t,u=\pm}u\,dV_{1}\left(\!L,u,E^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\!\right),$	(15)
$\displaystyle\!\!\!\!\!\!n_{\rm B}^{\rm 2f}$	$\displaystyle=$	$\displaystyle-{\partial\Omega_{\rm 2f}\over\partial\mu_{\rm B}}=2\int{{\rm d}^{3}k\over(2\pi)^{3}}\!\sum_{t,u=\pm}u\,dV_{1}\left(\!L,u,E^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\!\right).$	(16)

II.2 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_{\rm 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 Fukushima:2017csk ; Klevansky:1992qe ; Hatsuda:1994pi :

	$\displaystyle{\cal L}_{\rm NJL}\!$	$\displaystyle=$	$\displaystyle\!-\!V(L,L)\!+\!\bar{\psi}\!\left[i\not{\partial}\!-\!i\gamma^{4}\!\!\left(\!ig{\cal A}^{4}\!+\!Q_{\rm q}\mu_{\rm Q}\!+\!{\mu_{\rm B}\over 3}\!\right)\!-\!m_{0}\right]\!\psi$		(17)
			$\displaystyle+G\sum_{a=0}^{8}\left[(\bar{\psi}\lambda^{a}\psi)^{2}+(\bar{\psi}i\gamma_{5}\lambda^{a}\psi)^{2}\right]+{\cal L}_{\rm tH},$		(17)

where $\psi=(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

	$\displaystyle m_{0}$	$\displaystyle\equiv$	$\displaystyle{\rm diag}(m_{\rm 0u},m_{\rm 0d},m_{\rm 0s}),$
	$\displaystyle Q_{\rm q}$	$\displaystyle\equiv$	$\displaystyle{\rm diag}(q_{\rm u},q_{\rm d},q_{\rm s})={1\over 3}{\rm diag}(2,-1,-1);$		(18)

the interaction index $\lambda^{0}=\sqrt{2/3}~{}\mathbbm{1}_{3}$ and $\lambda^{i}~{}(i=1,\dots,8)$ are Gell-Mann matrices in flavor space. For later use, the ’t Hooft term ${\cal L}_{\rm tH}\equiv-K\sum_{t=\pm}{\rm Det}~{}\bar{\psi}\Gamma^{t}\psi$ can be represented as

\displaystyle\!\!\!\!\!\!\!\!{\cal L}_{\rm tH}\!=\!-{K\over 2}\sum_{t=\pm}\epsilon_{ijk}\epsilon_{imn}(\bar{\psi}^{i}\Gamma^{t}{\psi}^{i})(\bar{\psi}^{j}\Gamma^{t}{\psi}^{m})(\bar{\psi}^{k}\Gamma^{t}{\psi}^{n})

(19)

with the interaction vertices $\Gamma^{\pm}=\mathbbm{1}_{4}\pm\gamma_{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:

\sigma_{\rm f}=\langle\bar{\psi}_{\rm f}\psi_{\rm f}\rangle,\ \Delta_{\pi}=\langle\bar{u}i\gamma^{5}d\rangle,\ \Delta_{\pi}^{*}=\langle\bar{d}i\gamma^{5}u\rangle.

For brevity, we set $\Delta_{\pi}=\Delta_{\pi}^{*}$ without loss of generality in the following. To facilitate the study, we’d like first to reduce ${\cal L}_{\rm 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 six-fermion interaction term Klevansky:1992qe , we immediately find

	$\displaystyle{\cal L}_{\rm tH}^{4}\!$	$\displaystyle=$	$\displaystyle\!-{K}\!\left\{\epsilon_{ijk}\epsilon_{imn}\sigma_{i}\!\left(\bar{\psi}^{j}{\psi}^{m}\bar{\psi}^{k}{\psi}^{n}\!-\!\bar{\psi}^{j}i\gamma^{5}{\psi}^{m}\bar{\psi}^{k}i\gamma^{5}{\psi}^{n}\right)\!+\right.$		(20)
			$\displaystyle\left.\!\!\!\!2\Delta_{\pi}\!\!\left[\bar{s}{s}\!\left(\bar{u}i\gamma^{5}d\!+\!\bar{d}i\gamma^{5}u\!-\!\Delta_{\pi}\right)\!+\!\bar{s}i\gamma^{5}{s}\left(\bar{u}d\!+\!\bar{d}u\right)\right]\right\},$		(20)

where the second term in the brace is induced by $\pi^{\pm}$ condensations. Armed with the reduced Lagrangian density:

	$\displaystyle{\cal L}_{\rm NJL}\!$	$\displaystyle=$	$\displaystyle\!-V(L,L)\!+\!\bar{\psi}\!\left[i\not{\partial}\!-\!i\gamma^{4}\!\!\left(\!ig{\cal A}^{4}\!+\!Q_{\rm q}\mu_{\rm Q}\!+\!{\mu_{\rm B}\over 3}\!\right)\!-\!m_{0}\right]\!\psi$		(21)
			$\displaystyle+G\sum_{a=0}^{8}\left[(\bar{\psi}\lambda^{a}\psi)^{2}+(\bar{\psi}i\gamma_{5}\lambda^{a}\psi)^{2}\right]+{\cal L}_{\rm tH}^{4},$		(21)

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

\displaystyle{\cal L}_{\rm NJL}^{2}\!\!=\!\bar{\psi}\left[i\not{\partial}\!-\!i\gamma^{4}\!\!\left(ig{\cal A}^{4}\!+\!Q_{\rm q}\mu_{\rm Q}\!+\!{\mu_{\rm B}\over 3}\right)\!-\!m_{\rm i}\!-\!i\gamma^{5}\lambda^{1}\Pi\right]\psi,

where the scalar and pseudoscalar masses are respectively

	$\displaystyle m_{i}$	$\displaystyle=$	$\displaystyle m_{0i}-4G\sigma_{i}+2K(\sigma_{j}\sigma_{k}+\Delta_{\pi}^{2}\delta_{i3}),$
	$\displaystyle\Pi$	$\displaystyle=$	$\displaystyle(-4G+2K\sigma_{3})\Delta_{\pi}$		(23)

with $i\neq j\neq k$ . The $G$ and $K$ dependent terms in Eq. (23) are from the $U_{A}(1)$ symmetric and anomalous interactions, respectively. According to Eq.(II.2), $s$ quark decouples from $u,d$ quarks, so the gap equation for $\sigma_{\rm s}$ can be simply given by Klevansky:1992qe :

\displaystyle\!\!\!\!\!\sigma_{\rm s}\!=\!\langle\bar{s}{s}\rangle\!=\!{\rm tr}\left[i\not{\partial}-i\gamma^{4}\!\!\left(ig{\cal A}^{4}\!+\!Q_{\rm q}\mu_{\rm Q}\!+\!{\mu_{\rm B}\over 3}\right)\!-\!m_{\rm s}\right]^{\!-\!1}\!\!.

(24)

However, the $u$ and $d$ light quarks couple with each other through the non-diagonal pseudoscalar mass $\Pi$ . Since $\mu_{\rm B}$ is usually small in early Universe, we can simply set

m_{\rm 0u}=m_{\rm 0d}\equiv m_{\rm 0l},\ \sigma_{\rm u}=\sigma_{\rm d}\equiv\sigma_{\rm 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

\displaystyle\!\!\!\!\!\!\Omega_{\rm bl}

\displaystyle=

\displaystyle\!-\!2N_{c}\!\!\int^{\Lambda}\!\!\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\!\!\left[\sum_{t=\pm}\!E_{\rm l}^{\rm t}(k)\!+\!\epsilon_{\rm s}(k)\right]\!\!-\!2T\!\!\int\!\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\!\!\sum_{\rm u=\pm}\!\!\left[\sum_{t=\pm}\!Fl\!\left(\!L,u,E_{\rm l}^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\!\right)\!\!+\!Fl\!\left(\!L,u,\epsilon_{\rm s}(k),{-\mu_{\rm Q}\!+\!\mu_{\rm B}\over 3}\!\right)\!\right]

(25)

with the energy functions defined by

\displaystyle\!\!\!\!\epsilon_{\rm i}(k)\!=\!\sqrt{k^{2}\!+m_{\rm i}^{2}},\ E_{\rm l}^{\rm t}(k)\!=\!\sqrt{\left[\epsilon_{\rm l}(k)\!+t{\mu_{\rm Q}\over 2}\right]^{2}\!+\Pi^{2}}.

(26)

Eventually, the coupled gap equations follow directly from the definitions of condensates:

	$\displaystyle\sigma_{\rm s}\equiv\langle\bar{s}{s}\rangle={\partial\Omega_{\rm bl}\over\partial m_{\rm s}},\ 2\sigma_{\rm l}\equiv\langle\bar{u}{u}\rangle+\langle\bar{d}{d}\rangle={\partial\Omega_{\rm bl}\over\partial m_{\rm l}},$
	$\displaystyle 2\Delta_{\pi}\equiv\langle\bar{u}i\gamma^{5}{d}\rangle+\langle\bar{d}i\gamma^{5}{u}\rangle={\partial\Omega_{\rm bl}\over\partial\Pi}$		(27)

and the minimal condition $\partial_{\rm L}[V(L,L)+\Omega_{\rm bl}]=0$ as Xia:2013caa

$\displaystyle\sigma_{\rm s}\!$	$\displaystyle=$	$\displaystyle\!-2N_{c}\int^{\Lambda}\!\!{{\rm d}^{3}k\over(2\pi)^{3}}{m_{\rm s}\over\epsilon_{\rm s}(k)}+2N_{c}\int\!\!{{\rm d}^{3}k\over(2\pi)^{3}}{m_{\rm s}\over\epsilon_{\rm s}(k)}\sum_{u=\pm}dV_{1}\left(\!L,u,\epsilon_{\rm s}(k),{-\mu_{\rm Q}\!+\!\mu_{\rm B}\over 3}\!\right),$	(28)
$\displaystyle 2\sigma_{\rm l}\!$	$\displaystyle=$	$\displaystyle\!-2N_{c}\int^{\Lambda}\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\sum_{t=\pm}{m_{\rm l}\over\epsilon_{\rm l}(k)}{\epsilon_{\rm l}(k)+t{\mu_{\rm Q}\over 2}\over E_{\rm l}^{\rm t}(k)}+2N_{c}\int\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\sum_{t,u=\pm}{m_{\rm l}\over\epsilon_{\rm l}(k)}{\epsilon_{\rm l}(k)+t{\mu_{\rm Q}\over 2}\over E_{\rm l}^{\rm t}(k)}dV_{1}\left(\!L,u,E_{\rm l}^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\!\right),$	(29)
$\displaystyle 2\Delta_{\pi}\!$	$\displaystyle=$	$\displaystyle\!-2N_{c}\int^{\Lambda}\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\sum_{t=\pm}{\Pi\over E_{\rm l}^{\rm t}(k)}+2N_{c}\int\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\sum_{t,u=\pm}{\Pi\over E_{\rm l}^{\rm t}(k)}dV_{1}\left(\!L,u,E_{\rm l}^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\!\right),$	(30)

	$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!T^{4}\left[-\left(3.51\!-{2.47\over\tilde{T}}\!+\!{15.2\over\tilde{T}^{2}}\right)L\!+\!{1.75\over\tilde{T}^{3}}{12L(1-L)^{2}\over 1\!-\!6L^{2}\!+\!8L^{3}\!-\!3L^{4}}\right]=6T\!\int\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\!\sum_{u=\pm}\left[\sum_{t=\pm}dV_{2}\left(L,u,E^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\right)\right.$
	$\displaystyle\left.+dV_{2}\left(\!L,u,\epsilon_{\rm s}(k),{-\mu_{\rm Q}\!+\!\mu_{\rm B}\over 3}\!\right)\right].$		(31)

Note that $\Delta_{\pi}=0$ is a trivial solution of Eq.(30), so $\Delta_{\pi}$ or $\Pi$ is still a true order parameter for $I_{3}$ isospin symmetry Son:2000xc in three-flavor case. The total self-consistent thermodynamic potential can be found to be

	$\displaystyle\Omega_{\rm 3f}$	$\displaystyle=$	$\displaystyle V(L,L)+\Omega_{\rm bl}+2G(\sigma_{\rm s}^{2}+2\sigma_{\rm l}^{2}+2\Delta_{\pi}^{2})$		(32)
			$\displaystyle-4K(\sigma_{\rm l}^{2}+\Delta_{\pi}^{2})\sigma_{\rm s}$		(32)

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

$\displaystyle\!\!\!\!\!\!s_{\rm 3f}$	$\displaystyle=$	$\displaystyle 2\!\!\int\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\!\sum_{t,u=\pm}\left[Fl\left(\!L,u,E^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\!\right)\!+\!{3\left(\!E^{\rm t}(k)\!-\!u\,{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\!\right)\over T}dV_{1}\left(\!L,u,E^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\!\right)\right]$	(33)
		$\displaystyle+2\!\!\int\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\!\sum_{u=\pm}\left[Fl\left(\!L,u,\epsilon_{\rm s}(k),{-\mu_{\rm Q}\!+\!\mu_{\rm B}\over 3}\!\right)\!+\!{3\left(\!E^{\rm t}(k)\!-\!u\,{-\mu_{\rm Q}\!+\!\mu_{\rm B}\over 3}\!\right)\over T}dV_{1}\left(\!L,u,\epsilon_{\rm s}(k),{-\mu_{\rm Q}\!+\!\mu_{\rm B}\over 3}\!\right)\right]$
		$\displaystyle+T^{3}\left\{{1\over 2}\left(4\times 3.51-3\times{2.47\over\tilde{T}}+2\times{15.2\over\tilde{T}^{2}}\right)L^{2}+{1.75\over\tilde{T}^{3}}\ln\left[1-6L^{2}+8L^{3}-3L^{4}\right]\right\},$
$\displaystyle\!\!\!\!\!\!n_{\rm Q}^{\rm 3f}$	$\displaystyle=$	$\displaystyle N_{c}\!\!\int^{\Lambda}\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\!\sum_{t=\pm}t{\epsilon(k)\!+\!t{\mu_{\rm Q}\over 2}\over E^{\rm t}(k)}-\!3\!\!\int\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\sum_{t,u=\pm}t{\epsilon(k)\!+\!t{\mu_{\rm Q}\over 2}\over E^{\rm t}(k)}dV_{1}\left(\!L,u,E^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\!\right)$	(34)
		$\displaystyle+\int\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\!\sum_{t,u=\pm}u\,dV_{1}\left(\!L,u,E^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\!\right)-2\int{{\rm d}^{3}k\over(2\pi)^{3}}\!\sum_{t,u=\pm}u\,dV_{1}\left(\!L,u,\epsilon_{\rm s}(k),{-\mu_{\rm Q}\!+\!\mu_{\rm B}\over 3}\!\right),$	(34)
$\displaystyle\!\!\!\!\!\!n_{\rm B}^{\rm 3f}$	$\displaystyle=$	$\displaystyle 2\int{{\rm d}^{3}k\over(2\pi)^{3}}\!\sum_{t,u=\pm}u\,dV_{1}\left(\!L,u,E^{\rm t}(k),{\mu_{\rm Q}\!+\!2\mu_{\rm B}\over 6}\!\right)+2\int{{\rm d}^{3}k\over(2\pi)^{3}}\!\sum_{t,u=\pm}u\,dV_{1}\left(\!L,u,\epsilon_{\rm s}(k),{-\mu_{\rm Q}\!+\!\mu_{\rm B}\over 3}\!\right).$	(35)

II.3 The QEWD sector: free gases

In free gas approximation, the thermodynamic potentials for the QEWD sector can be easily given by Kapusta2006

$\displaystyle\Omega_{\gamma}$	$\displaystyle=$	$\displaystyle 2T\int{{\rm d}^{3}k\over(2\pi)^{3}}\log\left(1-e^{-k/T}\right),$	(36)
$\displaystyle\Omega_{\rm l}$	$\displaystyle=$	$\displaystyle-T\!\!\sum^{\rm i=e,\mu,\tau}_{u=\pm}\!\int\!\!{{\rm d}^{3}k\over(2\pi)^{3}}\!\left\{2\log\left[1+e^{-\left(\epsilon_{\rm i}(k)-u\,(-\mu_{\rm Q}+\mu_{\rm i})\right)/T}\right]\right.$	(37)
		$\displaystyle\left.+\log\left[1+e^{-(k-u\,\mu_{\rm i})/T}\right]\right\},$	(37)

where the degeneracy is one for neutrinos and anti-neutrinos due to their definite chiralities. Then, the corresponding entropy, electric charge number and lepton flavor number densities can be derived directly as

$\displaystyle s_{\gamma}$	$\displaystyle=$	$\displaystyle 2\int{{\rm d}^{3}k\over(2\pi)^{3}}\left[-\log\left(1-e^{-k/T}\right)+{k/T\over e^{k/T}-1}\right],$	(38)
$\displaystyle s_{\rm l}$	$\displaystyle=$	$\displaystyle\sum^{\rm i=e,\mu,\tau}_{u=\pm}\int{{\rm d}^{3}k\over(2\pi)^{3}}\left\{2\log\left[1+e^{-\left(\epsilon_{\rm i}(k)-u\,(-\mu_{\rm Q}+\mu_{\rm i})\right)/T}\right]\right.$	(39)
		$\displaystyle+\log\left[1+e^{-(k-u\,\mu_{\rm i})/T}\right]\!+\!{2\left(\epsilon_{\rm i}(k)\!-\!u\,(\!-\mu_{\rm Q}\!+\!\mu_{\rm i})\right)/T\over 1+e^{\left(\epsilon_{\rm i}(k)-u\,(-\mu_{\rm Q}+\mu_{\rm i})\right)/T}}$
		$\displaystyle\left.+{(k-u\,\mu_{\rm i})/T\over 1+e^{(k-u\,\mu_{\rm i})/T}}\right\},$

$\displaystyle n_{\rm Q}^{\rm l}$	$\displaystyle=$	$\displaystyle 2T\sum^{\rm i=e,\mu,\tau}_{u=\pm}\int{{\rm d}^{3}k\over(2\pi)^{3}}{-u\over 1+e^{\left(\epsilon_{\rm i}(k)-u\,(-\mu_{\rm Q}+\mu_{\rm i})\right)/T}},$	(40)
$\displaystyle n_{\rm i}$	$\displaystyle=$	$\displaystyle-{\partial\Omega_{\rm l}\over\partial\mu_{\rm i}}=T\sum_{u=\pm}\int{{\rm d}^{3}k\over(2\pi)^{3}}\left[{2u\over 1+e^{\left(\epsilon_{\rm i}(k)-u\,(-\mu_{\rm Q}+\mu_{\rm i})\right)/T}}\right.$	(41)
		$\displaystyle\qquad\qquad\ \ \ \left.+{u\over 1+e^{(k-u\,\mu_{\rm i})/T}}\right],\ i=e,\mu,\tau.$	(41)

Now, collecting contributions from both QEWD and QCD sectors, the total entropy, electric charge number and lepton number densities are respectively

\displaystyle\!\!\!\!\!\!s=s_{\gamma}\!+\!s_{\rm l}\!+\!s_{\rm 2f/3f},\,n_{\rm Q}=n_{\rm Q}^{\rm l}\!+\!n_{\rm Q}^{\rm 2f/3f},\,n_{\rm l}=\!\sum_{\rm i=e,\mu,\tau}\!\!n_{\rm i}

(42)

in the QCD epoch. To better catch the expansion nature of early Universe, we define several reduced quantities:

\displaystyle b=n_{\rm B}^{\rm 2f/3f}/s,\ l=n_{\rm l}/s,\ l_{\rm i}=n_{\rm i}/s

(43)

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$ Kolb:1983ni . Furthermore, due to neutrino oscillations at late stage of the Universe, the lepton flavor densities $l_{\rm i}$ are not well constrained at the QCD epoch, thus we will take $l_{\rm e}$ and $l_{\rm\mu}$ 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_{\mu}=113\,{\rm MeV}$ , simply set the electron mass $m_{\rm e}=0$ , and suppress the contribution of heavy $\tau$ leptons for the QEWD sector. The model parameters are fixed for the QCD sector as the following Zhuang:1994dw ; Rehberg:1995kh

$\displaystyle{\rm PNJL_{2f}}\!:$	$\displaystyle\ m_{\rm 0l}\!=\!5\,{\rm MeV},\,\Lambda\!=\!653\,{\rm MeV},\,G\Lambda^{2}\!=\!2.10;$	(45)
$\displaystyle{\rm PNJL_{3f}}\!:$	$\displaystyle\ m_{\rm 0l}\!=\!5.5\,{\rm MeV},\,m_{\rm 0s}\!=\!140.7\,{\rm MeV},$
	$\displaystyle\ \Lambda\!=\!602.3\,{\rm MeV},\,G\Lambda^{2}\!=\!1.835,\,K\Lambda^{5}\!=\!12.36.$

The $T-(l_{\rm e}+l_{\mu})$ phase diagrams of two-and three-flavor PNJL models are illuminated in Fig.1. As we can see, the ratio $l_{\rm e}/l_{\mu}$ 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_{\rm e}+l_{\mu}|\sim 0.1$ in the extrapolated lattice QCD calculations Vovchenko:2020crk , the values are consistently $|l_{\rm e}+l_{\mu}|\sim 0.09$ at their top temperature $T_{\rm pc}=0.16\,{\rm 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\sim 0.21\,{\rm GeV}$ which is much larger than $T_{\rm pc}$ , a well-known drawback of PNJL model Xiong:2009zz . 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_{\rm pc}$ in Ref. Vovchenko:2020crk . By the way, here the threshold lepton flavor asymmetries are $0.08$ and $0.075$ in the two- and three-flavor cases, respectively.

Refer to caption — Figure 1: The $T-(l_{\rm e}+l_{\mu})$ phase diagrams in two- (upper panel) and three-flavor (lower panel) PNJL models. The heavy and light shadows correspond to the pion superfluidity phase for $l_{\rm e}=0$ and $l_{\rm e}=l_{\mu}$ respectively with the total lepton asymmetry $l=-0.012$ . The dashed line in the upper panel is the boundary for the total lepton asymmetry $l=-51/28b$ and $l_{\rm e}=l_{\mu}$ .

Now, we take the more realistic three-flavor PNJL model for example to show the features of cosmic trajectories at $l_{\rm e}+l_{\mu}=-0.2$ , where the early Universe could evolve through the pion superfluidity phase. We compare two cases: $l_{\rm e}=0$ and $l_{\rm e}=l_{\mu}$ , 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_{\rm e}/l_{\mu}$ is only important on the cosmic trajectories of the directly related quantities: $\mu_{\rm e}$ and $\mu_{\mu}$ , and the results almost overlap with each other for other quantities. According to the lower panel of Fig.2, the pion condensate $\Pi$ 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 $|\mu_{\rm 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 Xiong:2009zz .

Following the monotonous feature of $|\mu_{\rm Q}|$ in the upper panel of Fig.3, we find the criterion $|\mu_{\rm Q}|=m_{\pi}^{\rm v}$ to be well satisfied at the lower critical temperature $T_{\rm l}$ ; but that is no longer useful for the exploration of upper critical temperature $T_{\rm u}$ , since the effective $\pi^{\pm}$ mass increases with $T$ . Since the $|\mu_{\rm Q}|$ varies from $m_{\pi}^{\rm v}$ to $\sim 3\,m_{\pi}^{\rm v}$ within the pion superfluidity phase, we can well recognize the BCS-BEC crossover therein Sun:2007fc 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_{\rm u}=m_{\rm d}$ and $L=L^{*}$ . Furthermore, the opposite signs between $\mu_{\rm B}$ and $\mu_{\rm Q}$ and the abrupt jump of $\mu_{\rm B}$ at low temperature both qualitatively fit the findings in the extrapolated lattice QCD study Wygas:2018otj . We note that the exist of pion superfluidity at $T_{\rm l}$ only leaves visible sign in the behavior of $\mu_{\rm Q}$ but the entrance at $T_{\rm u}$ gives rise to important signs in all the chemical potentials. Compared to the signs ${\cal S}(\mu_{\rm e})={\cal S}(\mu_{\mu})={\cal S}(\mu_{\tau})$ for the choice $l_{\rm e}=l_{\mu}=l_{\tau}$ in Ref. Wygas:2018otj , we find ${\cal S}(\mu_{\rm e})={\cal S}(\mu_{\mu})=-{\cal S}(\mu_{\tau})$ due to the choice $|l_{\rm e}+l_{\mu}|\gg|l|$ .

IV Summary

In this work, the possibility of pion superfluidity and the corresponding cosmic trajectories are self-consistently explored by varying the lepton flavor asymmetries within the PNJL model. The effects of strange quarks, total lepton asymmetry and the ratio $l_{\rm e}/l_{\mu}$ are all found to be mild on the phase boundary of pion superfluidity. Following the previous study in Ref. Vovchenko:2020crk , the phase boundary is constrained from both lower and upper sides in our study. At $T_{\rm cp}$ , the lepton flavor asymmetry $|l_{\rm e}+l_{\mu}|$ is $0.09$ in our work, quite consistent with the threshold value $0.1$ obtained in Ref. Vovchenko:2020crk . 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 $\mu_{\rm 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 $|\mu_{\rm Q}|$ at $T_{\rm u}$ is found to be around $3\,m_{\pi}^{\rm v}$ , which is so large that it explains why the extrapolated lattice QCD study was unable to fix $T_{\rm u}$ at all Vovchenko:2020crk .

In Ref.Vovchenko:2020crk , the equation of state of QCD+QEWD matter with different $|l_{\rm e}+l_{\mu}|$ was adopted to study the relic density of primordial gravitational wave; then inversely the observations of GW would help to constrain $|l_{\rm e}+l_{\mu}|$ 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 Coleman:1977py ; Callan:1977pt ; Witten:1984rs . Actually, with strong magnetic field presented in early Universe Vachaspati:1991nm ; Son:1998my ; Grasso:2000wj , 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.

Acknowledgments G.C. is supported by the National Natural Science Foundation of China with Grant No. 11805290. L.H. is supported by the National Natural Science Foundation of China with Grant Nos. 11775123 and 11890712. P.Z. is supported by the National Natural Science Foundation of China with Grant No. 11975320.

References

(1) E. V. Shuryak, Phys. Rept. 61, 71-158 (1980).
(2) T. D. Lee, eConf C8206282, 202-206 (1982) CU-TP-226.
(3) U. W. Heinz and M. Jacob, [arXiv:nucl-th/0002042 [nucl-th]].
(4) J. Adams et al. [STAR], Phys. Rev. Lett. 95, 122301 (2005).
(5) P.K. Kovtun, D. T. Son and A. O. Starinets, Phys. Rev. Lett. 94, 111601 (2005).
(6) P.F. Kolb and U. Heinz, Quark-Gluon Plasma 3, eds. R.C. Hwa and X-N Wang, Singapore: World Sci. (2004).
(7) Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, Nature 443, 675-678 (2006).
(8) T. Bhattacharya, M. I. Buchoff, N. H. Christ, H. T. Ding, R. Gupta, C. Jung, F. Karsch, Z. Lin, R. D. Mawhinney and G. McGlynn, et al. Phys. Rev. Lett. 113, no.8, 082001 (2014).
(9) M. Floris, Nucl. Phys. A 931, 103-112 (2014).
(10) L. Adamczyk et al. [STAR], Phys. Rev. C 96, no.4, 044904 (2017).
(11) X. Luo and N. Xu, Nucl. Sci. Tech. 28, no.8, 112 (2017).
(12) D. T. Son and M. A. Stephanov, Phys. Rev. Lett. 86, 592-595 (2001).
(13) J. B. Kogut and D. K. Sinclair, Phys. Rev. D 66, 034505 (2002).
(14) B. B. Brandt, G. Endrodi and S. Schmalzbauer, Phys. Rev. D 97, no.5, 054514 (2018).
(15) L. y. He, M. Jin and P. f. Zhuang, Phys. Rev. D 71, 116001 (2005).
(16) G. f. Sun, L. He and P. Zhuang, Phys. Rev. D 75, 096004 (2007).
(17) G. Cao, L. He and X. G. Huang, Chin. Phys. C 41, no.5, 051001 (2017).
(18) J. O. Andersen and L. Kyllingstad, J. Phys. G 37, 015003 (2009).
(19) H. Abuki, M. Ciminale, R. Gatto, N. D. Ippolito, G. Nardulli and M. Ruggieri, Phys. Rev. D 78, 014002 (2008).
(20) H. Abuki, T. Brauner and H. J. Warringa, Eur. Phys. J. C 64, 123-131 (2009).
(21) M. M. Middeldorf-Wygas, I. M. Oldengott, D. Bödeker and D. J. Schwarz, [arXiv:2009.00036 [hep-ph]].
(22) V. Vovchenko, B. B. Brandt, F. Cuteri, G. Endrődi, F. Hajkarim and J. Schaffner-Bielich, Phys. Rev. Lett. 126, no.1, 012701 (2021).
(23) P. A. R. Ade et al. (Planck Collaboration), Astron. Astrophys. 594, A13 (2016).
(24) I. M. Oldengott and D. J. Schwarz, Europhys. Lett. 119, 29001 (2017).
(25) M. M. Wygas, I. M. Oldengott, D. Bödeker and D. J. Schwarz, Phys. Rev. Lett. 121, no.20, 201302 (2018).
(26) P. N. Meisinger and M. C. Ogilvie, Phys. Lett. B 379, 163-168 (1996).
(27) K. Fukushima, Phys. Lett. B 591, 277-284 (2004).
(28) C. Ratti, M. A. Thaler and W. Weise, Phys. Rev. D 73, 014019 (2006).
(29) K. Fukushima and V. Skokov, Prog. Part. Nucl. Phys. 96, 154-199 (2017).
(30) S. P. Klevansky, “The Nambu-Jona-Lasinio model of quantum chromodynamics,” Rev. Mod. Phys. 64, 649 (1992).
(31) T. Hatsuda and T. Kunihiro, Phys. Rept. 247, 221 (1994).
(32) J. Xiong, M. Jin and J. Li, J. Phys. G 36, 125005 (2009).
(33) T. Xia, L. He and P. Zhuang, Phys. Rev. D 88, no.5, 056013 (2013).
(34) J. I. Kapusta and C. Gale, Finite-Temperature Field Theory: Principles and Applications, (Cambridge University Press, 2006).
(35) E. W. Kolb and M. S. Turner, Ann. Rev. Nucl. Part. Sci. 33, 645-696 (1983).
(36) P. Zhuang, J. Hufner and S. P. Klevansky, Nucl. Phys. A 576, 525 (1994).
(37) P. Rehberg, S. P. Klevansky and J. Hufner, Phys. Rev. C 53, 410 (1996).
(38) S. R. Coleman, Phys. Rev. D 15, 2929 (1977) Erratum: [Phys. Rev. D 16, 1248 (1977)].
(39) C. G. Callan, Jr. and S. R. Coleman, Phys. Rev. D 16, 1762 (1977).
(40) E. Witten, Phys. Rev. D 30, 272-285 (1984).
(41) T. Vachaspati, Phys. Lett. B 265, 258-261 (1991).
(42) D. T. Son, Phys. Rev. D 59, 063008 (1999).
(43) D. Grasso and H. R. Rubinstein, Phys. Rept. 348, 163-266 (2001).