CHAPTER: ZERO TEMPERATURE ANALYSIS

In this chapter
we see the results of the numerical calculation of quark-meson model with vanishing temperature and finite density.
We have used the following vacuum parameters as the input:

\begin{gather}
m_\pi^{\rm vac} = 140 \MeV, ~~~ m_\sigma^{\rm vac.} = 600 \MeV, ~~~ f_\pi^{\rm vac} = 90 \MeV, ~~ M_q^{\rm vac} = 300 \MeV.
\end{gather}

This parameter set corresponds to the following coupling constants for mesons

\begin{gather}
g^{\rm vac} \simeq 3.333, ~~~ \lambda^{\rm vac} \simeq 126.1.
\end{gather}

The value of $\lambda^{\rm vac}$ is apparently too large to use the naive perturbation theory,
but this choice is the conventional one
and consistent with other analysis <cit.>.

For comparison to the lattice data in Ref.<cit.>,
we also examine $m_\pi^{\rm vac}= 170 \MeV$ case
keeping other parameters $(m_\sigma, f_\pi, M_q)$ unchanged.

§ EVOLUTION OF MICROSCOPIC QUANTITIES

we see the evolution of condensates as the function of $\mu_I$.
Fig.$\,$<ref> shows the constituent quark mass $M_q$ and the quark mass gap $\Delta$ associated with the pion condensate.
The results of tree-level are illustrated by the dashed lines and the one-loop results by the solid lines.
There are no illuminating differences in $M_q$ between tree-level and one-loop,
but the one-loop results of the quark mass gap $\Delta$ are significantly different from the tree-level ones;
the $\Delta$ increases linearly in $\mu_I$ at tree-level,
whereas it saturates $\Delta\sim 300 \MeV$ at one-loop.
This drastic change of $\mu_I$ dependence is beyond the range of the naive perturbation theory.
The reason should lie in the physical picture of mesons in the existence of quarks rather mere perturbative calculation.

The constituent quark mass $M_q$ and the quark mass gap $\Delta$ as the function of $\mu_I/m_\pi$.
The location of the second order phase transition $\mu_I = 70 \MeV$ is shown by the vertical dashed line.

At tree-level,
mesons are regarded as the fundamental degrees of freedom.
But adding the quark degrees of freedom,
they are no longer elementary and treated as composite particles made of quark and anti-quark.
Although the dynamical effects of quarks to meson fluctuations are neglected in our large $N_c$ approximation,
the quark degrees of freedom still have the critical roles on the condensed mesons.
Thus the evolution of the pion condensate is suppressed by the Pauli blocking of quarks.
Summed up,
the quark degrees of freedom couples to the condensed mesons
and change the quantitative behavior of the meson condensates,
by switching the fundamental degrees of freedom from mesons to quarks.

Now we closely examine the behavior of the quark mass gap $\Delta$
by considering the large chemical potential limit $\mu_I \to \infty$.
In this limit,
we assume chiral symmetry is asymptotically restored $M_q\to 0$.
Then the potential has the form of

\begin{align}
V_0 \to -\frac{2\mu_I^2}{g^2}\Delta^2 + \frac{\lambda}{24g^4}\Delta^4.
\end{align}

After solving the gap equation,
this form immediately leads that $\Delta$ behaves as $\Delta\sim \mu_I$ in this limit.
In the high density limit only relevant scale is $\mu_I$
and the potential should satisfy the conformal scaling as $P\sim \mu_I^4$.
Therefore this scaling of $\Delta \sim \mu_I$ is inevitable by $\mu_I^2\Delta^2$ and $\Delta^4$ terms.

The linearity of $\Delta$ at the tree-level changes to the saturation at one-loop.
This is because the one-loop potential disfavors the large $\Delta$.
The one-loop potential in the same limit is given by

\begin{align}
V_{\rm 1-loop} = V_0^{\msbar}  + \frac{2N_c}{(4\pi)^2}\qty[\qty(\Delta^4 - 4\mu_I^2\Delta^2)\ln\frac{\Lambda^2}{\Delta^2} + \frac{3}{2}\Delta^4] + V_q^{\rm fin}.
\end{align}

The $V_0^{\msbar}$ is the tree-level potential written in the renormalized parameters
and thus the structure of $\mu_I$ and $\Delta$ are unchanged;
the $\mu_I$ dependence is $V_0^{\msbar} \sim \mu_I^2 \Delta^2$.
The $V_q^{\rm fin}$ is the finite part of the quark single particle energy contribution.
We soon find that the coefficient of $\ln \Lambda^2/\Delta^2$ is negative
and that the one-loop potential disfavors the large $\Delta$.
In addition,
in the case of the tree-level the behavior $\Delta \sim \mu_I$ is necessary to satisfy the conformal scaling $P \sim \mu_I^4$,
as the quark single particle energy $V_q^{\rm fin}$ directly provides the $\mu_I^4$ term
when kinetic energy of quark dominates the potential $E_{\rm kin} \sim \mu_I \gg \Delta$,
the behavior $\Delta \to \rm{const}$ is allowed.

§ EQUATION OF STATE AND SOUND VELOCITY

Next we see the thermodynamic quantities.
Thermodynamic potential is given by the minimum of the thermodynamic potential $V_{\rm eff}$ as

\begin{align}
P(\mu_I) = -V_{\rm eff}(\mu_I;{M_q}_*(\mu_I),\Delta_*(\mu_I))
\end{align}

with the solution of gap equation ${M_q}_*$ and $\Delta_*$.
The isospin density and energy density are given by

\begin{align}
n_I = \pdv{P}{\mu_I}, ~~~ \epsilon = -P + \mu_I n_I.
\end{align}

We also see the sound velocity

\begin{align}
c_s^2 = \pdv{P}{\epsilon} = \frac{n_I}{\mu_I \chi_I}, ~~~ \chi_I = \pdv[2]{P}{\mu_I}.
\end{align}

From now on
we compute the numerical results with LQCD data in Refs.<cit.>.
The former is the results of $2+1$-flavor lattice QCD
which masses are at the physical point $m_\pi = 135 \MeV$.
The pion decay constant is $f_\pi \sim 92-96 \MeV$ (the  definition of $f_\pi$ differs from ours by factor $\sqrt{2}$ and the value above is converted to our definition).
The latter is the recent result of multi-pion correlation function
and the pion mass is $m_\pi \sim 170 \MeV$.
This calculation reaches to the high density region up to $\mu_I \sim 1 \GeV$ or $n_I \sim 400 n_0$
(again the definition of isospin chemical potential differs by factor $2$ and we corrected to ours).
Both of two lattice results are consistent to ChPT prediction at low density.

Fig,$\,$<ref> shows the pressure as a function of $\mu_I$.
Fig.$\,$<ref> shows the isospin density as a function of $\mu_I$
and Fig.$\,$<ref> the pressure as a function of energy density $\epsilon$.
We also presents the results of tree-level and LQCD.
For comparison to LQCD results,
we also show the results of $m_\pi = 170 \MeV$ case.

In the former figure,
we can see the both tree-level and one-loop results have the consistency to LQCD results at low density.
But the tree-level results overestimate the density and have large discrepancy with LQCD and one-loop results at high density $n_I \sim 50 n_0$.
This too rapid growth of the density at tree-level is the consequence of the $\Delta$ behavior $\Delta \sim \mu_I$ at high density.
In contrast,
the one-loop result with the pion mass of $m_\pi = 170\MeV$ matches to the LQCD results <cit.> in wide range of density.

Pressure as a function of $\mu_I$.
We also show the LQCD results in Refs.<cit.>.
The one-loop result of quark meson model with the pion mass $m_\pi = 170 \MeV$ matches to the LQCD results <cit.> in wide range of chemical potential.
The tree-level result increases rapidly
and has large discrepancy with one-loop and LQCD results.
The upper is for the global behavior
and the lower is for the low density behavior.

Isospin density as a function of $\mu_I$.
We also show the LQCD results in Refs.<cit.>.
The one-loop result of quark meson model with the pion mass $m_\pi = 170 \MeV$ matches to the LQCD results <cit.> in wide range of density.
The tree-level result increases rapidly at high density region
and has large discrepancy with one-loop and LQCD results.
The upper is for the global behavior
and the lower is for the low density behavior.

In the latter figure of pressure as a function of energy density, or EOS,
the one-loop pressure are slightly larger than that of tree-level.
Both of them are consistent to the LQCD results in wide range of density
in contrast to the isospin density.

Pressure as a function of $\varepsilon$.
We also show the LQCD results in Refs.<cit.>.
The upper is for the global behavior
and the lower is for the low density behavior.

In fig:zeroT-pressure-parameter,
we illustrated the EOS at the low energy with different model parameters
$(m_\sigma,f_\pi)=(450,90)$, $(600,90)$, $(450,100)$ and $(600,100)$ MeV.
We can see the variation of the model parameters does not change the results qualitatively.
Especially the $m_\sigma$ does not change the results much.
But the larger $m_\sigma$ and $f_\pi$, or stronger chiral symmetry breaking in vacuum, reduce the pressure.
We get back this point in Sec.$\,$<ref>.

The pressure as a function of $\varepsilon$ with different model parameters $(m_\sigma,f_\pi)$.
The results of $f_\pi = 90\MeV$ is shown by blue lines and those of $f_\pi = 100 \MeV$ by red lines.
The results of $m_\sigma = 450 \MeV$ is shown by dashed lines and those of $m_\sigma = 600 \MeV$ by solid lines.
With larger $m_\sigma$ and $f_\pi$ the pressure reduces.

For more detailed examination of the EOS,
we see the behavior of the squared sound velocity $c_s^2$ in Fig.$\,$<ref>.
The $c_s^2$ increases rapidly in low density region $n_I < 2 n_0$ and exceeds the conformal value $c_s^2 = 1/3$ around $c_s^2 \simeq 2 n_0$,
forms the peak at $n_I \simeq 5n_0$ and approaches to the conformal value from above at large density.
These qualitative features are consistent to the LQCD results in Refs.<cit.>.

However the detailed results in the BCS region are slightly different.
The peak position of $c_s^2$ is $\mu_I/m_\pi \sim 1.3$ in our results which is consistent to the LQCD results in <cit.>,
while it is $\mu_I/m_\pi \sim 0.8$ in <cit.>.
The height of the peak has the same tendency;
it is $c_s^2 \sim 0.5$ in our results and $c_s^2 \sim 0.3$ in <cit.>
and the results of <cit.> has agreement with our result.

the rapidity with which the $c_s^2$ decreases after making a peak is different between our results and <cit.>.
In our results the $c_s^2$ decreases slowly and large $c_s^2 \simeq 0.5$ remains at high density region,
while in Ref.<cit.> the $c_s^2$ decreases rapidly.
This rapid decrease to the conformal value indicates that the transition from hadronic to quark matter rapidly takes place.
The origin of this discrepancy in the sound velocity is not clear at present.

The squared sound velocity $c_s^2$ as a function of $n_I/n_0$.
The qualitative behavior of $c_s^2$ is consistent to the both LQCD results,
but the peak position and the height of the peak are slightly different with <cit.>.
Comparing to Ref. <cit.> our results are consistent in wide range of density.

We should notice that the $c_s^2$ at tree-level also has the peak structure in the crossover region
and relax to the conformal value from above.
As we have seen in Sec.$\,$<ref>
the term $a_2\mu_I^2$ in the pressure with positive $a_2$ in addition to $\mu_I^4$ term
concludes that the $c_s^2$ converges to $1/3$ from above.
Solving the gap equation at tree-level,
we can find that the tree-level pressure has the form

\begin{align}
P_0(\mu_I) = \frac{m_\pi^4 f_\pi^2}{8\mu_I^2} + \frac{f_\pi^2 \qty(4\mu_I^2 - m_0^2)^2}{2(m_\sigma^2 - m_\pi^2)}.
\end{align}

This pressure has the $\mu_I^2$ term with positive coefficient,
thus the $c_s^2$ can exceed the conformal value $1/3$.
In addition the chiral restoration term $\sim \mu_I^{-2}$ suppresses $c_s^2$ at low density.
Combining these two effects,
the $c_s^2$ has the peak structure even at tree-level.
Even though the quantitative behavior is similar with one-loop level,
the peak structure at tree-level is physically fictitious.
This is because the peak structure of $c_s^2$ in middle density region is just the result of the competition between the dynamics of low and high density region.
As we have mentioned in Sec.$\,$<ref>,
the underlying dynamics of high density region is not correctly described at tree-level.
Thus we conclude that the peak structure at tree-level is not physical.

The convergence to the conformal value from above is different tendency with pQCD
and we discuss this point in Sec.$\,$<ref>.

§ QUARK SATURATION

As mentioned,
the low density physics is governed by meson degrees of freedom
and the internal structure of mesons is not appeared.
But at high density the distance of these mesons becomes short
and the quark degrees of freedom appears in its dynamics.
The physics of intermediate region is determined by the competition between these two degrees of freedom at low and high density.

To estimate the density where the quark degrees of freedom becomes important,
we introduce the overlap density of pion $n_I^{\rm overlap}$.
This is the density where the typical distance of pions becomes comparable to the pion size.
The experimental determination based on the $\pi e$ scattering and the $e^+e^- \to \pi^+\pi^-$ process <cit.>
gives the estimate $\braket{r^2}_V = 0.434(5) \fm^2$ <cit.>, or

\begin{align}
r_\pi^V = \sqrt{\braket{r^2}_V} \sim 0.66 \fm,
\end{align}

which has been well reproduced by the lattice QCD calculation <cit.>.
Then the overlap density is given by

\begin{align}
n_I^{\rm overlap} = \qty(\frac{4}{3}\pi r_\pi^3)^{-1} \sim 0.83 {\rm fm^{-3}} \sim 5.2 n_0.
\end{align}

Fig.$\,$<ref> shows this density corresponds to

\begin{align}
\mu_I \sim 1.2 m_\pi \sim 168 \MeV.
\end{align}

We should notice that the quark wavefunction is percolated outside of $r_\pi^V$
and the internal quarks starts interacting with other quarks in other mesons at lower density.
This suggests that there should be more suitable measure which contains the information of quark wavefunction.
Following Refs.<cit.>,
we interpret the peak structure as the signature of the quark saturation
or the onset of the quark matter.
This quark saturation is characterized by the occupation probability of the quark matter;
the quark states are occupied with the probability $\sim 1$ from low momentum
and the newly added quarks by increasing the density must fill a state on top of the already occupied states.

the occupation probability of the quark matter is calculated through the test quark propagator,
as shown in fig:zeroT-occupation-cartoon.
Because of the Pauli blocking,
the distribution of test quark propagator reflects the occupation probability of the quark matter.
With the pion condensate in the background,
the low momentum mode propagation of quark is suppressed
while that of the high momentum mode is not.
This means that the background condensate consists of the low-momentum quark
and suppresses the low-momentum mode of the quark propagator through the Pauli blocking.

The test quark propagator in the pion condensed phase and the occupation probability $f(p)$.
The $f(p)$ for the propagator reflect the occupation probability of the quark matter.

The quark occupation probability in the pion condensed phase can be computed in the standard Nambu-Gor'kov formalism.
The derivation is reviewed in Appendix <ref>.
The occupation probability of $u, ~ d, ~ \bar u$ and $\bar d-$ quarks are given by

\begin{align}
f(p) &= f_{u,\bar d}(p) = \frac{1}{2}\qty(1 + \frac{\mu_I - E_D(p)}{E(p)}) \\
\bar f(p) &= f_{d,\bar u}(p) = \frac{1}{2}\qty(1 + \frac{\mu_I + E_D(p)}{E(p)}).
\end{align}

In the absence of mass gap $\Delta$,
$u-$ and $\bar d-$ quarks are occupied up to $p = \mu_I$
and the $d-$ and $\bar u-$ quarks are empty.
When the mass gap $\Delta$ is turned on,
the distribution of $u-$ and $\bar d-$ quarks becomes smoothly reducing distribution around $p\sim\mu_I$.

The Fig.$\,$<ref> shows the quark occupation probability as a function of momentum $p$
with isospin density from $n_I/n_0 = 0.2$ to $10$ in $0.2$ increments for gray curves and $1.0$ for red curves.
The blue dots are the location of the surface of the distribution
where $f(p)$ has the half maximum value; $f(p_{\rm surf}) = f(0)/2$.

Quark occupation probability as a function of momentum $p$.
The gray curves are the results of $n_I/n_0 = 0.2 - 10.0$ in $0.2$ increments.
The red curves are the results of $1.0$ increments.
The blue dots are the location of the surface of the distribution
where $f(p)$ has the half maximum value.

To discuss the transition of the effective degrees of freedom in the context of the quark saturation,
it is useful to decompose the evolution of $f(p)$ into two modes.
The first mode is the vertical evolution (Fig.$\,$<ref>, Left)
in which the distribution $f(p)$ increases in monotonic way as $f(p) \sim n_I \varphi_\pi^{\rm vac}(p)$ with the quark distribution $\varphi_\pi^{\rm vac}(p)$ in a single pion.
This regime corresponds to the ideal pion gas regime at low density
where the pions are not so close to each other
and the quarks inside pions are largely unaffected by other quarks in other pions.
In this region $\epsilon/n_I$ is almost constant
and the pressure $P=n_I^2\partial(\epsilon/n_I)/\partial n_I$ is very small.
Here the quarks have the contributions to energy density but can not increase the pressure.
The second is the horizontal evolution (Fig.$\,$<ref>, Right)
in which the distribution $f(p)$ increases and the surface of $f(p)$ extends to the large momentum region.
This region corresponds to the ordinary quark matter region
where the quark interactions and Pauli blocking become important.
In this region the $\epsilon/n_I$ increases with $n_I$.
Thus the pressure can increase and the sound velocity becomes large.

Schematic figures for the evolution of the occupation probability.
Left: The vertical evolution;
Right: The horizontal evolution.

As we can see in Fig.$\,$<ref>,
the calculated distribution of $f(p)$ is the mixture of these two modes in our quark meson model.
In the low density region from $0$ to $\sim 3.0 n_0$,
the magnitude of $f(p)$ at $p=0$ grows rapidly and reaches the magnitude $\sim 0.8$,
which is considered as the result of the vertical evolution.
Even in this region we can see that
the vertical evolution is gradually decelerating
and the horizontal evolution is slightly increasing.
If the momentum of increase near $n_I=0$ continues,
the value of $f(p=0)$ exceed $1$ at $n_I \simeq 2.0 n_0$,
as shown in Fig.$\,$<ref>.
But the Pauli blocking does not allow $f(p)$ to exceed $1$,
the distribution must extend into the large momentum region.
Thus the horizontal evolution becomes dominant as increasing the density
while the vertical evolution is still remaining.
This smooth transition from the vertical to horizontal evolution suggests the smooth transition of the effective degrees of freedom <cit.>.

The evolution of the occupation probability for $p = 0$ (red line).
If the momentum of increase near $n_I=0$ continues (blue dashed line),
the value of $f(p=0)$ exceed $1$ at $n_I \simeq 2.0 n_0$.

In order to observe more closely the transition from the BEC phase to the BCS phase,
we examine how much the horizontal evolution of $f(p)$ is realized in the intermediate region $n_I \sim 2.0 - 5.0 n_0$.
This is illustrated in Fig.$\,$<ref>
by taking the ratio $f(p)_{n_I}/f(p)_{n_I=0.10n_0}$.
The vertical evolution keeps this ratio flat over all range of momentum.
In contrast the horizontal evolution increases this ratio at the large momentum region.
To characterize the evolution of the Fermi surface,
we trace $p_{\rm peak}$,
the position of the maximum of $f(p)/\varphi_\pi^{\rm vac}(p)$.
In practice,
we take $f(p)_{n_I=0.1n_0}$ instead of $\varphi_\pi^{\rm vac}$
expecting $n_I = 0.1n_0$ is low enough to be in the vertical evolution regime;
$f(p)_{n_I=0.1n_0} = 0.1n_0 \varphi_\pi^{\rm vac}(p)$.
The distribution of $f(p)_{n_I=0.1n_0}$ is shown in Fig.$\,$<ref>
and described well by a simple monopole Ansatz or Breit-Wigner form

\begin{align}
\varphi_\pi(p) \sim \frac{1}{1 + p^2/a^2}.
\end{align}

The value of $a$ is $\simeq 87.6$ MeV as determined by fit.
This form suggests the pion wavefunction in the coordinate space has the exponentially-decaying form.

The ratio $f(p)_{n_I}/f(p)_{n_I=0.10n_0}$ as a function of momentum $p$.
The gray and red curves are the same as Fig.$\,$<ref>.
The blue dots are the location of the peak of $f(p)_{n_I}/\varphi_\pi^{\rm vac}(p)$.

The occupation probability $f(p)$ at $n_I = 0.1 n_0$.
The red line is the numerical result
and the blue dashed line is the fit by the monopole Ansatz.
This wavefunction can be regarded as the occupation probability of the vacuum $\varphi_\pi(p)$
and is the reference distribution for the ratio $f(p)_{n_I}/\varphi_\pi(p)$ in Fig.$\,$<ref>.

Fig.$\,$<ref> shows the position $p_{\rm peak}$ is almost constant in the low density region up to $n_I \sim 2.0 n_0$.
Beyond this density the $p_{\rm peak}$ starts to increase and the horizontal evolution becomes dominant.
The density region of horizontal evolution $n \gtrsim 2.0n_0$ is consistent with substantial increase of the sound velocity $c_s^2 \gtrsim 1/3$.
The onset of the horizontal evolution and the quark degrees of freedom in dynamics is relatively early comparing to the overlap density of pion $n_I^{\rm overlap}$.

We also examine the transition of the effective degrees of freedom by the half maximum value of $f(p)$
as illustrated in Fig.$\,$<ref>
As increasing the density or chemical potential,
the width of the distribution becomes wider.
But in the vertical evolution the maximum value of $f(p)$ also increases.
Thus the value of $p_{\rm surf}$ does not increase much.
On the contrary,
in the vertical evolution where the low-momentum mode of the quarks are almost saturated,
the value of $f(p_{\rm surf})$ does not change
while the position of $p_{\rm surf}$ increases.
The transition point of these two evolution is around $n_I \simeq 2.0 n_0$,
which is the consistent result with the previous analysis.

This transition density is much smaller than the pion overlap density $n_I^{\rm overlap} \simeq 5.2n_0$.
This is because the pions does not have the solid boundary and the quark wavefunction is percolated outside of $r_\pi^V$.

§ DISCUSSION

In this section we have short discussion on several issues which are not fully mentioned in the previous sections.
First we discuss the connection of the chiral symmetry restoration and softening of EOS at high density.
Next we compare our results to the pQCD calculation and discuss the possible origin of the discrepancy of the asymptotic behavior of sound velocity at high density region.
Finally we discuss the trace anomaly
$\Delta_{\rm tr} = 1/3 - P/\varepsilon$ as the another measure of the conformality.

§.§ Chiral symmetry restoration and softening

In the Sec.$\,$<ref> we have seen that larger $f_\pi$ or $m_\sigma$ make the EOS softer at high density.
Here we discuss the connection of this softening to the chiral symmetry breaking in the vacuum and its restoration at high density region.
The stronger chiral symmetry breaking corresponds to larger $f_\pi$ and $m_\sigma$ in our model.

The energy density as a function of the chiral effective mass $M_q$.
After chiral restoration the minimum energy is realized at $M_q = 0$,
but with the broken chiral symmetry,
the minimum is realized at $M_q \neq 0$ and the energy is smaller than that of $M_q = 0$.
This gap in the zero-point energy density is the bag constant.

The chiral symmetry breaking in the vacuum shifts the zero-chemical-potential value of thermodynamic potential as

\begin{align}
B \equiv \Omega(M_q=0) - \Omega(M_q = M_q^{\rm vac}),
\end{align}

where the first term is the thermodynamic potential of the trivial vacuum
and the second term is that of the chiral symmetry broken vacuum.
The $B$ is the energy difference between these two vacuum which is called the bag constant.
As the bag constant is defined by the difference of the thermodynamic potential with symmetry-broken and unbroken vacuum,
its magnitude increases with stronger chiral symmetry breaking (fig:zeroT-bagconst2).

In our model the bag constants $B(m_\sigma,f_\pi)$ in the $m_q^4$ unit are given as

\begin{align}
B(600, 90) = 0.389, ~~~ B(600, 100) = 0.483 , \notag\\
B(800, 90) = 0.394, ~~~ B(800, 100) = 0.489,
\end{align}

In our model the bag constants $B(m_\sigma,f_\pi)$ in the $m_\pi^4 = (140 \MeV)^4$ unit are given as

\begin{align}
B(450, 90) = 0.874, ~~~ B(450, 100) = 1.035 , \notag\\
B(600, 90) = 1.439, ~~~ B(600, 100) = 1.698, % in mpi^4 unit (mpi=140 MeV)
\end{align}

from which we can see stronger chiral symmetry breaking increase the bag constant $B$.

The energy density with different strength of the chiral symmetry breaking.
The bag constant is larger for the stronger chiral symmetry breaking.

To see how the bag constant affects the stiffness of EOS,
it is useful to consider the EOS in the bag model with perturbative corrections.
Since EOS is normalized by $P=0$ at $\mu=0$ and $T=0$ for non-perturbative vacuum,
the perturbative evaluation of EOS should be normalized by the bag constant $B$ as

\begin{align}
P_{\rm pert}^{\rm normalized} = P_{\rm pert} - B, ~~~ \varepsilon_{\rm pert}^{\rm normalized} = \varepsilon_{\rm pert} + B.
\end{align}

This shows that the bag constant $B$ associated with the chiral symmetry breaking in the vacuum
contributes to the EOS as the negative pressure and the positive energy density,
which makes the EOS softer.
The similar results on the sound velocity has been obtained the $2+1$ flavor NJL model <cit.>.

§.§ Comparison to pQCD with power corrections

Now we get back the asymptotic behavior of the sound velocity at high density region.
In our calculation of QM model,
the $c_s^2$ converges to the conformal value $1/3$ from above as increasing density.
This behavior is qualitatively different from the pQCD prediction,
which is $c_s^2 \to 1/3$ from below.
The origin of this discrepancy is the existence of non-perturbative pressure $\sim \Lambda^2\mu_I^2$ in our model.
The QCD scale $\Lambda_{\rm QCD}$ has the non-analytic form for QCD coupling constant $\alpha_s$ as

\begin{align}
\Lambda_{\rm QCD} = Q \exp\qty(-\frac{2\pi}{\beta_0 \alpha_s(Q^2)})
\end{align}

The $\mu_I^2$ term in the pressure need the energy scale
but in the high density region the only relevant scale is $\Lambda_{\rm QCD}$.
Thus pQCD can not include the pressure $\sim \Lambda^2\mu_I^2$.
And as illustrated in Sec.$\,$<ref>,
this $\mu_I^2$ term makes the sound velocity $c_s^2$ converge to $1/3$ from above.

To examine in detail,
we calculate the pressure in the pQCD with power corrections
and how large power corrections are needed to change the asymptotic behavior of $c_s^2$.

The pQCD pressure up to $\rmO(\alpha_s^2)$ for single quark flavor $f$ is given as <cit.>

\begin{align}
P_0^f  &= \frac{N_c}{ 12\pi^2 } \qty[|\mu_f| u_f  \qty( \mu_f^2 - \frac{ 5 }{2} M_q^2 ) + \frac{ 3 }{2} M_q^4 \ln \frac{ |\mu_f| + u_f }{ M_q } ] ,\\
P_1^f  &= - \frac{ \alpha_s N_{\rm G} }{ 16 \pi^2 }
\begin{multlined}[t]
\bigg[ 3 \qty( M_q^2 \ln \frac{ |\mu_f| + u_f }{ M_q }  - |\mu_f| u_f )^2 - 2 u_f^4 \\
+ M_q^2 \qty( 6 \ln \frac{ \Lambda_{\rm reno} }{ M_q } + 4 ) \qty( |\mu_f| u_f - M_q^2 \ln  \frac{ |\mu_f| + u_f }{ M_q } ) \bigg]
\end{multlined}
\end{align}

with $u_f = \sqrt{\mu_f^2 -M_f^2}$ and $N_{\rm G} = N_c^2 -1$.
The $M_q$ is the current quark mass for u/d quarks $M_q = 5.0 \MeV$.
The $\Lambda_{\rm reno}$ is the renormalization scale which commonly taken as $0.5\mu_f$ to $2\mu_f$.
In the isospin symmetric case the each flavor contributes equally to the pressure; $P^u = P^d$.

The running coupling constant $\alpha_s$ is given as <cit.>

\begin{align}
\begin{split}
\alpha_s (\Lambda_{\rm reno}) & = \frac{4 \pi}{\beta_0 L} \bigg[ 1 - 2 \frac{\beta_1}{\beta_0^2} \frac{\ln L}{L} %
+ \frac{\beta_1^2}{\beta_0^4 L^2} \qty( \ln^2 L -\ln L - 1 + \frac{\beta_2\beta_0}{\beta_1^2} ) \\
& \quad + \frac{\beta_1^3}{\beta_0^6 L^3} \qty( - \ln^3 L + \frac{5}{2}\ln^2 L + 2\ln L %
- \frac{1}{2} - 3 \frac{\beta_2 \beta_0}{\beta_1^2} \ln L + \frac{\beta_3 \beta_0^2}{2\beta_1^3} ) \bigg]
\end{split}
\end{align}

with $L = 2 \ln( \Lambda_{\rm reno}/\Lambda_{\msbar})$, $\beta_0 = 11-2N_f/3$, $\beta_1 = 102-19N_f/3$,

\begin{align}
\beta_2 {} &= \frac{ 2857 }{ 2 } - \frac{ 5033 }{ 18 }N_f + \frac{ 325 }{ 54 }N_f^2, \\
\begin{split}
\beta_3 {} &=
\bigg( \frac{ 149753 }{ 6 } + 3564 \xi(3) \bigg) %
- \bigg( \frac{ 1078361 }{ 162 } + \frac{ 6508 }{ 27 }\xi(3) \bigg) N_f \\
& \quad + \bigg( \frac{ 50065 }{ 162 } + \frac{ 6472 }{ 81 }\xi(3) \bigg) N_f^2 %
+ \frac{ 1093 }{ 729 } N_f^3
\end{split}
\end{align}

and $\Lambda_{\msbar} \simeq 340  \MeV$.

Notice that this formula of $\alpha_s$ is valid only for $\Lambda_{\rm reno} > \Lambda_{\msbar}$.
We extrapolated $\alpha_s$ for smaller $\Lambda_{\rm reno}$ by using the Gaussian formula <cit.>

\begin{align}
\alpha_s^{\rm low} = \alpha_s^{\rm low} (0) e^{-Q^2/ 4\kappa^2}
\end{align}

with $\alpha_s^{\rm low}(0) \simeq 1.22$ and $\kappa \simeq 0.51$.
In the aim of calculating the sound velocity and thus the second derivative of the pressure,
we interpolate $\alpha_s$ by the smooth function $\alpha_s^{\rm mid}$ as

\begin{align}
\alpha^{\rm low}_s(Q^2) = \theta(t_{\rm low}-Q^2) \alpha^{\rm mid}_s(Q^2) + \theta(t_{\rm high}-Q^2) \theta(Q^2-t_{\rm low}) + \alpha^{\rm high}_s(Q^2) \theta(Q^2 -t_{\rm high})
\end{align}

with $t_{\rm low}^{1/2} = 0.3 \GeV$ and $t_{\rm high}^{1/2} = 1.1 \GeV$.
The boundary condition requires that up to second derivative of $\alpha_s$ are smooth, which is

\begin{align}
\left.\frac{\partial^n \alpha_s^{\rm low/high}}{(\partial Q^2)^n}\right|_{Q^2 = t_{\rm low/high}} = \left.\frac{\partial^n \alpha_s^{\rm mid}}{(\partial Q^2)^n}\right|_{Q^2 = t_{\rm low/high}}
\end{align}

for $n=0,1,2$.
And we choose $\alpha_s^{\rm mid}$ the polynomial function as

\begin{align}
\alpha_s^{\rm mid}(Q^2) = \sum_{m=0} c_m \mu_I^m.
\end{align}

Six coefficients $c_m$ are fully determined by six boundary conditions.

Fig.$\,$<ref> illustrates the running coupling constant $\alpha_s$,
both the original pQCD formula and our interpolation.
We set $Q^2 = \Lambda_{\rm reno}^2$ and $\Lambda_{\rm reno} = X\mu_I$
and the scale factor $X$ are varied as $0.5, 1.0, 2.0$.
The pQCD pressure has the artificial rapid reduction in the low density region,
but freezing coupling keeps the pressure finite in the IR region
as shown in Fig.$\,$<ref>.

The running coupling constant $\alpha_s$ as a function of $\mu_I$.
The solid curves are the results of interpolated expression.
The dashed curve is the original pQCD formula.
We took $X = 0.5, 1.0$ and $2.0$.

The perturbative pressure with the coupling constant $\alpha_s$ in Fig.$\,$<ref>.
The pressure is normalized by the Stefan-Boltzmann limit $P_{\rm SB} = \pi^2 n_I^2/90$.
The notation of the curves is same as Fig.$\,$<ref>.

Now we examine the impact of the power corrections to the sound velocity.
The power correction comes from the gap in the pion condensed phase.
The phase space factor $\sim 4\pi p_f^2 \Delta$ multiplying the gap $\Delta$ dividing $(2\pi)^3$ for each $u, \dbar$-quarks
gives the naive estimate of the power correction to the pressure

\begin{align}
P_{\rm power} = C \frac{\mu_I^2 \Delta^2}{\pi^2}
\end{align}

with overall factor $C = \calO(1)$.
The similar examination has been done in Ref.<cit.>,
which consider the diquark condensate $\Delta$ in 3-color QCD using the functional renormalization group.

We can compute the gap $\Delta$ analytically in the weak-coupling and high density limit.
In the case of the chiral symmetric gap with $J^P = 0^+$
associating with the two-flavor color-superconductor (2SC) pairing,
the form of the gap is given as <cit.>

\begin{align}
\Delta \sim \mu_I g^{-5} \exp\qty(-\frac{3\pi^2}{\sqrt{2}g_s})
\end{align}

with running coupling constant $g_s = g_s(|\mu_I|)$.
In the case of the BCS gap at vanishing temperature
we have the following relation with gap $\Delta_{\rm BCS}$ and transition temperature $T_{\rm BCS}$

\begin{align}
T_{\rm BCS} \simeq 0.57 \Delta_{\rm BCS}.
\label{eq:zeroT-BCStemperature}
\end{align}

In our calculation the gap is about $\Delta\simeq 300 \MeV$ at high density
which corresponds to $T_{\rm BCS} \simeq 170 \MeV$.
This value has a good agreement with the LQCD result <cit.> $T_c^{\rm lattice} \simeq 161 \MeV$.
The pQCD pressure with ad-hoc power correction is shown in Fig.$\,$<ref>.
We choose the size of gap $\Delta = 0, 200$ and $300 \MeV$.
With the original pQCD running coupling constant,
the pressure decreases to zero around $\mu_I \sim 0.3 \GeV$ as lowering the density.
This is the result of the IR divergence of the running coupling constant.
The freezing coupling constant keeps the pressure non-zero in the low density region.
The finite gap $\Delta$ makes the pressure seemingly diverge at low density,
while the pressure with vanishing gap $\Delta = 0$ is finite.

The perturbative pressure with ad-hoc power correction as a function of isospin density $n_0$.
The scale factor $X$ is fixed to $1.0$.
The solid and dashed line represents the freezing coupling and original pQCD formula, as in Fig.$\,$<ref>.

The sound velocity $c_s^2$ is shown in Fig.$\,$<ref>.
We varied the scale factor $X = 0.5 - 2.0$ which is shown in the band.
With the vanishing gap $\Delta = 0$,
we can see the sound velocity increases monotonically and converges to $1/3$ from below.
Adding the small gap $\Delta = 200 \MeV$ makes the sound velocity larger than $c_s^2 = 1/3$,
but the convergence from below is still alive at high density $n_I > 100 n_0$.
In contrast, with the large gap $\Delta = 300 \MeV$,
the sound velocity changes the trend and converges to $1/3$ from above.

The sound velocity $c_s^2$ with ad-hoc power correction as a function of isospin density $n_0$.
The band shows the variation of the scale factor $X = 0.5 - 2.0$.
The solid line is the result of $X = 1.0$.

As increasing the density
the power correction becomes relatively small comparing to the perturbative pressure.
This can be easily seen by taking ratio of these two terms as

\begin{align}
\frac{P_{\rm power}}{P_{\rm pQCD}} \sim \qty(\frac{\Delta}{\mu_I})^2.
\end{align}

Around $\mu_I \simeq 1 \GeV$ or $n_I \simeq 80 n_0$,
the power correction has the small contribution in the total pressure

\begin{align}
\frac{P_{\rm power}}{P_{\rm pQCD}} \sim \qty(\frac{300 \MeV}{1 \GeV})^2 = 0.09,
\end{align}

but change the qualitative behavior of the sound velocity.
This suggests that the power correction has the essential role even at the high density regime
and can be the main source of the discrepancy of the asymptotic behavior of the sound velocity.

§.§ Trace anomaly

We have discussed the sound velocity $c_s^2$ to closely examine the EOS especially its stiffening.
Trace anomaly is another important quantity to characterize the EOS
and gets attention in the context of not only the NSs <cit.>
but also in the mechanical properties of hadrons <cit.>.
The sound velocity is defined as $c_s^2 = \pdv*{P}{\varepsilon}$ and contains only the information of slope of EOS.
But in the context of the NSs,
the trace anomaly is the relation of $P$ and $\varepsilon$
thus contains the magnitude of EOS.

The trace anomaly is the indicator of the deviation from the conformal limit at high density.
It is defined as

\begin{align}
\braket{T^\mu_{~\mu}} = \varepsilon - 3P,
\end{align}

Then the conformal limit concludes $\Delta_{\rm tr} = 0$.
It is useful to introduce the normalized trace anomaly divided by $3\varepsilon$ as

\begin{align}
\Delta_{\rm tr} \equiv \frac{1}{3} - \frac{P}{\varepsilon}.
\end{align}

The trace anomaly $\Delta_{\rm tr}$ is conjectured to be positive in Ref.<cit.>,
but the positivity is broken with the non-perturbative effects as well as the sound velocity.
To examine the effect of the power correction to the trace anomaly,
we again use the parametric pressure as the function of $\mu_I$

\begin{align}
P_\textrm{with powers} = a_0 \mu_I^4 + a_2 \mu_I^2.
\label{eq:zeroT-parametric-pressure}
\end{align}

But this time we consider the running coefficients $a_0$ and $a_2$ as a function of $\mu_I$,
which comes from the $\mu_I$-dependence of  running coupling constant $\alpha_s(\mu_I)$.
Then the energy density $\varepsilon$ is derived as

\begin{align}
\varepsilon = 3a_0\mu_I^4 + 2a_2\mu_I^2 + \pdv{a_0}{\ln \mu_I}\mu_I^4 + \pdv{a_2}{\ln \mu_I}\mu_I^2
\end{align}

and the trace anomaly is given as

\begin{align}
\braket{T^\mu_{~\mu}}_\textrm{with powers} &= -2a_2\mu_I^2 + \pdv{a_0}{\ln \mu_I}\mu_I^4 + \pdv{a_2}{\ln \mu_I}\mu_I^2.
\end{align}

The original pQCD pressure favors the positive trace anomaly,
but as is the case for the sound velocity
the attractive power correction $a_2 > 0$ can make the trace anomaly negative.

Our numerical result of the trace anomaly as a function of $n_I/n_0$ is shown in Fig.$\,$<ref>
with the LQCD result <cit.>.
The trace anomaly starts with the positive value at $n_I \simeq 0$
and decreases to under zero around $n_I \simeq 3 n_0$.
At the high density around $n_I \simeq 10 n_0$,
$\Delta_{\rm tr}$ starts to increase and converges to zero from below.

The trace anomaly $\Delta_{\rm tr}$ as a function of isospin density $n_I / n_0$.
The qualitative behavior of the trace anomaly is consistent with the LQCD results <cit.>.

The trace anomaly calculated by the pQCD adding the power correction is shown in Fig.$\,$<ref> and Fig.$\,$<ref>.
along with our model calculations.
Fig.$\,$<ref> shows the trace anomaly without the power correction with different renormalization scale $X$,
which has the positive value in the whole range of $\mu_I$, independent on $X$.
But adding the power correction with $\Delta \sim 200 \MeV$ or larger makes the trace anomaly negative
as shown in Fig.$\,$<ref>.
These results suggests that the sign of the trace anomaly is sensitive to the power correction
and can be useful indicator of the non-perturbative effects.

The trace anomaly $\Delta_{\rm tr}$ without gap $\Delta$ as a function of isospin density $n_I$.
The notation of the curves is same as Fig.$\,$<ref>.

The trace anomaly with the power correction $\Delta = 0, 200, 300$ and $400  \MeV$.
The notation of the curves is same as Fig.$\,$<ref>.
Nonzero gap favors the negative trace anomaly.