# Sound velocity peak and conformality in isospin QCD

Ryuji Chiba<EMAIL_ADDRESS>Department of Physics, Tohoku
University, Sendai 980-8578, Japan  Toru Kojo<EMAIL_ADDRESS>Department of Physics, Tohoku University, Sendai 980-8578, Japan

###### Abstract

We study zero temperature equations of state (EOS) in isospin QCD within a
quark-meson model which is renormalizable and hence eliminates high density
artifacts in models with the ultraviolet cutoff (e.g., NJL models). The model
exhibits a crossover transition of pion condensations from the Bose-Einstein-
Condensation regime at low density to the Bardeen-Cooper-Schrieffer regime at
high density. The EOS stiffens quickly and approaches the quark matter regime
at density significantly less than the density for pions to spatially overlap.
The sound velocity develops a peak in the crossover region, and then gradually
relaxes to the conformal value $1/3$ from above, in contrast to the
perturbative QCD results which predicts the approach from below. In the
context of QCD computations, this opposite trend is in part due to the lack of
gluon exchanges in our model, and also due to the non-perturbative power
corrections arising from the condensates. We argue that with large power
corrections the trace anomaly can be negative. In quantitative level, our EOS
is consistent with the lattice results in the BEC regime but begins to get
stiffer at higher density. The sound velocity peak also appears at higher
density. The BCS gap in our model is $\Delta\simeq 300$ MeV in the quark
matter domain, and naive application of the BCS relation for the critical
temperature $T_{c}\simeq 0.57\Delta$ yields the estimate $T_{c}\simeq 170$
MeV, in good agreement with the lattice data.

###### pacs:

##  I Introduction

The quantum chromodynamics (QCD) with a large isospin chemical potential
($\mu_{I}$) can be studied in lattice Monte-Carlo simulations and hence has
been a useful laboratory to test theoretical conceptions in dense matter [1,
2, 3]. In this theory, the positive isospin chemical potential favors the
population of up-quarks and of down-antiquarks. A matter with finite isospin
density starts with a Bose-Einstein-Condensation (BEC) phase of charged pions
as composite particles. The dilute regime is well-described by chiral
effective theories (ChEFT) for pions [4, 5, 6, 7, 8]. As density increases the
quark substructure of pions should become important and the system transforms
into a Bardeen-Cooper-Schrieffer (BCS) phase with a substantial quark Fermi
sea. This BEC-BCS transition is crossover (for BEC-BCS crossover, see,
e.g.,[9, 10, 11]), as confirmed by model studies and lattice Monte-Carlo
simulations. We study this crossover in the context of quark-hadron continuity
[12, 13, 14] or duality [15], which may be also realized in QCD at finite
baryon chemical potential ($\mu_{B}$).

One of fundamental topics in dense QCD is the equation of state (EOS)
(Ref.[16] for a short review from the QCD perspective), which is a crucial
piece to understand the structure of neutron stars (NSs). Recent analyses of
NSs, with similar radii ($\simeq 12.4$ km) for $2.1$ and $1.4$ solar mass NSs
[17, 18, 19], and the nuclear physics constraints at nuclear saturation
density $n_{0}$ ($\simeq 0.16\,{\rm fm}^{-3}$), suggest that the EOS stiffens
rapidly (i.e., the pressure $P$ grows rapidly as a function of energy density
$\varepsilon$) from low baryon density to high density, of $n_{B}=4$-$7n_{0}$,
which is expected to be realized in the core of massive NSs. This stiffening
accompanies the peak of sound velocity, $c_{s}=(\partial
P/\partial\varepsilon)^{1/2}$; the $c_{s}^{2}$ is $\ll 1$ in the nuclear
domain, goes beyond the so-called conformal value $1/3$, and relaxes to $1/3$
in the relativistic limit where the quark kinetic energy dominates over the
interaction. While NS observations suggest such non-monotonic behaviors of
$c_{s}$, it is necessary to understand the mechanisms from the microscopic
physics. The sound velocity peak was first indicated in phenomenological
interpolation of hadronic and quark matter EOS [20, 21, 22, 23], discussed in
more general ground based on nuclear physics and NS observations by Ref.[24],
and further elucidated in Refs. [25, 26, 27] utilizing the quark degrees of
freedom. Recently more detailed descriptions have been attempted, see
Refs.[28, 29, 30, 31, 32, 33, 34, 35, 36, 37], but it is difficult to directly
test the scenarios. We use the isospin QCD for which lattice simulations are
available, and delineate the behavior of $c_{s}^{2}$.

Another interesting question is how the $c_{s}^{2}$ approaches the conformal
limit, $1/3$. Perturbative QCD (pQCD) [38, 39, 40, 41, 42, 43], which is
supposed to be valid at $n_{B}\gtrsim 40n_{0}$, predicts that the $c_{s}^{2}$
approaches $1/3$ from below. The domain between $n_{B}\simeq 10n_{0}$ and
$\simeq 40n_{0}$ has not been explored intensively. For this regime it is
natural to regard quarks as relevant degrees of freedom but whose properties
may be substantially renormalized by strong interaction effects [44, 45, 46];
if such interaction effects are properly absorbed into effective parameters of
quasi-quarks, the residual interactions may be treated in the same spirit as
in constituent quark models for hadron physics. If this residual corrections
are indeed smaller than the relativistic kinetic energy of quasi-particles,
the system should show the conformal behavior even before achieving weakly
correlated quark matter. How the matter reaches the conformal regime is
intensively discussed in recent works [47, 48, 49, 50, 51].

We address the non-monotonic behavior of $c_{s}^{2}$ in isospin QCD within a
renormalizable quark-meson model. The properties of the model in isospin QCD
have been analyzed in detail by Refs.[52, 53, 54]. We follow their
renormalization procedures. The advantage of using renormalizable effective
models over models with a UV cutoff (e.g., the NJL models) is that one can
temper the high density artifacts. In particular, the BCS type states have a
distorted quark occupation probability whose high momentum tail reaches very
high momenta, exceeding the UV cutoff. This is in contrast to the ideal gas
case with the occupation probability $\theta(p-p_{f})$ which discontinuously
drops to zero at the Fermi momentum $p_{f}$ before reaching the UV cutoff. In
fact, NJL studies with BCS states exhibit growing $c_{s}^{2}$ toward the high
density limit [55, 56]. In the quark-meson model such growing behaviors
disappear; the $c_{s}^{2}$ relaxes to the conformal value $1/3$, as it should.

While our model predicts $c_{s}^{2}\rightarrow 1/3$, the conformal limit is
reached from above, not from below as predicted in pQCD. The latter is due to
the density dependence induced through the running $\alpha_{s}$. In the weak
coupling limit and at high density, the only relevant scale is $\mu_{I}$ and
the $c_{s}^{2}=1/3$ follows from $\partial(P/\mu_{I}^{4})/\partial\mu_{I}=0$.
The first important corrections come from the $\Lambda_{\rm QCD}\simeq$
200-300 MeV in the running coupling constant. If we take into account
$\Lambda_{\rm QCD}$ only in this way, the $c_{s}^{2}$ is reduced from $1/3$.
Meanwhile, at the energy scale around $\sim 1\,{\rm GeV}$, it has been long
known that power corrections of $\Lambda_{\rm QCD}$, which can not be
expressed as perturbative series in $\alpha_{s}$, play important roles to
capture the qualitative features in QCD [57, 58, 59]. Parameterizing pressure
with power corrections as ($\chi_{I}$: isospin susceptibility)

$\displaystyle P(\mu_{I})=a_{0}\mu_{I}^{4}+a_{2}\mu_{I}^{2}\,,$ (1)

where $a_{2}\sim\Lambda_{\rm QCD}^{2}$, the squared sound velocity can be
expressed as

$\displaystyle
c_{s}^{2}=\frac{\,n_{I}\,}{\,\mu_{I}\chi_{I}\,}=\frac{\,2a_{0}\mu_{I}^{2}+a_{2}\,}{\,6a_{0}\mu_{I}^{2}+a_{2}\,}\,.$
(2)

For a positive $a_{2}$, the $c_{s}^{2}$ is larger than 1/3, and close to 1 if
the $a_{2}$ term dominates. In our quark-meson model the $a_{2}$ is related to
the pion condensates. We quantify the relation within our quark-meson model.
In particular we confront our model calculations with the lattice results,
finding reasonable agreement in the BEC regime but deviation at higher
density.

We briefly address the trace anomaly in dense matter which measures the
breaking of the scale invariance [47, 48, 49]. We argue that changes from the
non-perturbative to perturbative vacua add positive contributions to the trace
anomaly, while the power corrections with $a_{2}>0$ favors the negative value.
For large power corrections the trace anomaly can be negative. The sign of the
trace anomaly is very useful to characterize the non-perturbative effects in
dense quark matter.

In this paper we use nuclear saturation density in QCD, $n_{0}=0.16\,{\rm
fm}^{-3}$, as our unit for isospin density. While there is no need to address
nuclear saturation in isospin QCD, our goal is to discuss physics as a step to
understand NS EOS and $n_{0}$ is useful in this phenomenological context.

This paper is organized as follows. In Sec.II we discuss our set up for a
quark meson model. The renormalization procedures are summarized. In Sec.III
we present the renormalized thermodynamic potential and resulting EOS, as well
as the correlation between condensates and EOS. We emphasize the importance of
quark substructure which can be seen only after including quark loops. The
numerical results are confronted with the lattice data. In Sec.IV we discuss
the zero point energy in EOS which often appears as the bag constant in a
phenomenological model. In our quark-meson model this quantity can be computed
explicitly. In addition we discuss the power correction to the pQCD. The
evolution of the sound velocity at high density is presented. We also discuss
the trace anomaly as the indicator of the non-perturbative effects. Section V
is devoted to a summary.

##  II Model

The Lagrangian of the two-flavor quark-meson model is

$\displaystyle\mathcal{L}$ $\displaystyle={1\over
2}\left(\partial_{\mu}\vec{\phi}\right)^{2}-{1\over
2}m_{0}^{2}\vec{\phi}^{2}-{\lambda\over 24}(\vec{\phi}^{2})^{2}+h\sigma$
$\displaystyle-\overline{\psi}\left({\ooalign{\hfil/\hfil\crcr$\partial$}}-g(\sigma+\mathrm{i}\gamma^{5}\vec{\tau}\cdot\vec{\pi})\right)\psi\,,$
(5)

where $\psi$ is a quark field with up- and down-quark components

$\displaystyle\psi=\left(\begin{matrix}u\\\ d\end{matrix}\right).$ (6)

The $\vec{\phi}=(\sigma,\vec{\pi})$ are meson fields which correspond to the
isospin $\bf 1$ and $\bf 3$ representations. The $\tau_{i}$’s are the Pauli
matrices in flavor space.

We compute the thermodynamic potential at finite isospin density $n_{I}$,
utilizing the isospin chemical potential $\mu_{I}$ as a Lagrange multiplier.
To correctly identify the corresponding Lagrangian, we should begin with the
hamiltonian formalism. The thermodynamic potential is

$\displaystyle\Omega=H-\mu_{I}N_{I}\,,~{}~{}~{}~{}~{}~{}~{}~{}N_{I}=\int_{\bm{{x}}}n_{I}\,.$
(7)

The isospin density in terms of field variables can be identified by the
Noether theorem. Meson and quark fields transform under isospin
transformations as

$\displaystyle\pi_{a}\mapsto\exp(\mathrm{i}\theta_{i}T_{i})^{\rm
adj}_{ab}\pi_{b}\,,~{}~{}~{}~{}~{}\psi\mapsto\mathrm{e}^{\mathrm{i}\theta_{i}\tau_{i}}\psi\,,$
(8)

and corresponding conserved current can be written as

$\displaystyle
j^{\mu}_{a}=\epsilon_{abc}\pi_{b}\partial^{\mu}\pi_{c}+\delta^{\mu
0}\delta_{a3}\bar{\psi}\gamma^{0}\tau_{3}\psi\,,$ (9)

where the $\epsilon_{abc}$ is the complete anti-symmetric tensor with
$\epsilon_{123}=1$. The isospin density is now

$\displaystyle
n_{I}=j^{0}_{a=3}=\pi_{+}\mathrm{i}\partial^{0}\pi_{-}-\pi_{-}\mathrm{i}\partial^{0}\pi_{+}+\bar{\psi}\gamma^{0}\frac{\tau_{3}}{2}\psi\,.$
(10)

Writing fields collectively as $\Phi=(\vec{\phi},\psi)$, The partition
function for $\Phi$ is

$\displaystyle
Z=\int\mathcal{D}\Pi_{\Phi}\mathcal{D}\Phi\exp\left[\mathrm{i}\int_{x}\left(\dot{\Phi}\cdot\Pi_{\Phi}-\mathcal{H}+\mu_{I}n_{I}\right)\right].$
(11)

where $(\Pi_{\Phi})_{i}=\partial^{0}\Phi_{i}$ is a field conjugate to
$\Phi_{i}$. Keeping in mind that $n_{I}$ contains the conjugate fields
$\Pi_{\pi_{\pm}}=\partial^{0}\pi_{\pm}$, we integrate $\vec{\Pi}_{\phi}$ to
get

$\displaystyle
Z=\int\mathcal{D}\Phi\exp\left(\mathrm{i}\int_{x}\mathcal{L}^{B}_{\rm
dense}\right)\equiv\mathrm{e}^{\mathrm{i}\mathcal{V}_{4}V_{\rm eff}}\,,$ (12)

where $\mathcal{V}_{4}$ is the volume factor and $V_{\rm eff}$ is the
effective potential. Here the Lagrangian at finite density is

$\displaystyle\mathcal{L}_{B}(\vec{\phi}_{B})=$ $\displaystyle\,{1\over
2}\left[(\partial_{\mu}\sigma_{B})^{2}+(\partial_{\mu}\pi_{B3})^{2}\right]$
$\displaystyle+(\partial_{\mu}+2\mathrm{i}\mu_{I}\delta_{\mu}^{0})\pi_{B}^{+}\left(\partial^{\mu}-2\mathrm{i}\mu_{I}\delta^{\mu}_{0}\right)\pi_{B}^{-}$
$\displaystyle-{1\over 2}m_{0B}^{2}\vec{\phi}_{B}^{2}-{\lambda_{B}\over
24}(\vec{\phi}_{B}^{2})^{2}+h_{B}\sigma_{B}$
$\displaystyle+\overline{\psi}_{B}\left[\mathrm{i}{\ooalign{\hfil/\hfil\crcr$\partial$}}+\mu_{I}\frac{\tau_{3}}{2}\gamma^{0}-g_{B}(\sigma_{B}+\mathrm{i}\gamma^{5}\vec{\tau}\cdot\vec{\pi}_{B})\right]\psi_{B}\,.$
(15)

where we attached the subscript $B$ to emphasize the parameters and fields are
unrenormalized.

Below we construct a one-loop effective potential. First we rewrite the
Lagrangian using the renormalized parameters and fields. Following Ref.[52],
in this work we use the on-mass shell renormalization for $\sigma$ and $\pi$
fields, and decompose the renormalized and counter terms in an asymmetric way.
Explicitly,111Our parametrization differs from Ref.[52] by treatment of
$Z$-factors. For instance, Ref.[52] uses $m_{0B}^{2}=m_{0}^{2}+\delta
m_{0}^{2}$ while we write $Z_{B}m_{0B}^{2}=m_{0}^{2}+\delta m_{0}^{2}$. When
we decompose renormalized and counter terms, we always combine $Z$-factors
before decomposition.

$\displaystyle\sigma_{B}$ $\displaystyle=$ $\displaystyle
Z_{\sigma}^{1/2}\sigma\,,~{}~{}~{}~{}m_{0B}^{2}Z_{\sigma}=m_{0}^{2}+\delta
m_{0\sigma}^{2}\,,$ $\displaystyle\vec{\pi}_{B}$ $\displaystyle=$
$\displaystyle
Z_{\pi}^{1/2}\vec{\pi}\,,~{}~{}~{}~{}m_{0B}^{2}Z_{\pi}=m_{0}^{2}+\delta
m_{0\pi}^{2}\,,$ $\displaystyle\psi_{B}$ $\displaystyle=$ $\displaystyle
Z_{\psi}^{1/2}\psi\,,~{}~{}~{}~{}h_{B}Z_{\sigma}^{1/2}=h+\delta h\,,$ (16)

and

$\displaystyle\lambda_{B}Z_{\sigma}^{2}$ $\displaystyle=$
$\displaystyle\lambda+\delta\lambda_{4\sigma}\,,~{}~{}~{}\lambda_{B}Z_{\pi}^{2}=\lambda+\delta\lambda_{4\pi}\,,$
$\displaystyle\lambda_{B}Z_{\sigma}Z_{\pi}$ $\displaystyle=$
$\displaystyle\lambda+\delta\lambda_{2\sigma\pi}\,,$ $\displaystyle
g_{B}Z_{\sigma}^{1/2}Z_{\psi}$ $\displaystyle=$ $\displaystyle g+\delta
g_{\sigma}\,,~{}~{}~{}~{}~{}g_{B}Z_{\pi}^{1/2}Z_{\psi}=g+\delta g_{\pi}\,.$
(17)

The counter terms of $\sigma$ and $\vec{\pi}$ are not independent (see below)
but we write in a parallel way for book keeping purposes. The divergent
components of the counter terms respect the original $O(4)$ symmetry. We then
demand that the loop corrections to the renormalized field strength and masses
vanish at the on-mass shell. We carry out these steps in vacuum, and use the
obtained counter terms for finite density calculations.

With this renormalization scheme, the expression of the Lagrangian becomes
more complicated than the standard $O(4)$ symmetric renormalization, but the
merit is that we can discuss the vacuum condensates within the tree level
potential for the renormalized fields. We emphasize that we demand this
condition only for the vacuum renormalization procedure. At finite density the
loop corrections do modify the tree level relation.

Now the Lagrangian can be decomposed into the renormalized part and counter
terms as

$\displaystyle\mathcal{L}_{B}(\vec{\phi}_{B})=\mathcal{L}(\vec{\phi})+\mathcal{L}_{\rm
c.t.}(\vec{\phi})\,,$ (18)

where $\mathcal{L}$ is the renormalized Lagrangian where all subscripts $B$
are omitted from $\mathcal{L}_{B}$, and the counter Lagrangian is

$\displaystyle\mathcal{L}_{\rm c.t.}(\vec{\phi})$ $\displaystyle={\delta
Z_{\sigma}\over 2}(\partial_{\mu}\sigma)^{2}+{\delta Z_{\pi}\over
2}(\partial_{\mu}\pi_{3})^{2}$ $\displaystyle+\delta
Z_{\pi}(\partial_{\mu}+2\mathrm{i}\mu_{I}\delta_{\mu}^{0})\pi^{+}\left(\partial^{\mu}-2\mathrm{i}\mu_{I}\delta^{\mu}_{0}\right)\pi^{-}$
$\displaystyle-{1\over 2}\delta m_{0\sigma}^{2}\sigma^{2}-{1\over 2}\delta
m_{0\pi}^{2}\vec{\pi}^{2}+\delta h\sigma$ $\displaystyle-{1\over
24}\big{[}\delta\lambda_{4\sigma}(\sigma^{2})^{2}+\delta\lambda_{4\pi}(\vec{\pi}^{2})^{2}+2\delta\lambda_{2\sigma\pi}\sigma^{2}\vec{\pi}^{2}\big{]}$
$\displaystyle+\overline{\psi}\delta
Z_{\psi}\big{(}\mathrm{i}{\ooalign{\hfil/\hfil\crcr$\partial$}}+\mu_{I}\frac{\tau_{3}}{2}\gamma^{0}\big{)}\psi$
(21) $\displaystyle-\overline{\psi}\big{(}\delta
g_{\sigma}\sigma+\mathrm{i}\delta
g_{\pi}\gamma^{5}\vec{\tau}\cdot\vec{\pi}\big{)}\psi\,.$ (22)

The counter Lagrangian is used when we calculate loop corrections.

Actually it is often convenient to split the counter terms further into the
divergent and finite pieces. The divergent part is identified in the
$\overline{\rm MS}$ scheme. For instance,

$\displaystyle\delta\lambda_{\sigma}=\delta_{\overline{\rm
MS}}\lambda+\delta_{\rm f}\lambda_{\sigma}\,,$ (23)

and so on, and

$\displaystyle\mathcal{L}_{\rm c.t.}=\mathcal{L}_{\rm c.t.}^{\overline{\rm
MS}}+\mathcal{L}_{\rm c.t.}^{\rm f}\,.$ (24)

The $\mathcal{L}_{\rm c.t.}^{\overline{\rm MS}}$ respects the $O(4)$ symmetry
and the counter terms appear in a symmetric way for $\sigma$ and $\vec{\pi}$.

The UV finite counter terms differ for $\sigma$ and $\vec{\pi}$, but there are
relations. Starting with our definition of renormalization factors, we can
derive the relations

$\displaystyle-m_{0}^{2}\delta Z_{\sigma}+\delta
m_{0\sigma}^{2}=-m_{0}^{2}\delta Z_{\pi}+\delta m_{0\pi}^{2}\,,$ (25)

including only three independent counter terms, and

$\displaystyle\delta g_{\sigma}-\delta
g_{\pi}=\frac{\,g\,}{\,2\,}\big{(}\delta Z_{\sigma}-\delta Z_{\pi}\big{)}\,,$
(26)

including single independent counter term, and

$\displaystyle\delta\lambda_{4\sigma}-2\delta Z_{\sigma}\lambda$
$\displaystyle=$ $\displaystyle\delta\lambda_{4\pi}-2\delta Z_{\pi}\lambda$
(27) $\displaystyle=$ $\displaystyle\delta\lambda_{2\sigma\pi}-(\delta
Z_{\sigma}+\delta Z_{\pi})\lambda\,,$

including single independent counter term. In total we need five conditions to
fix five independent counter terms and the others can be deduced from the
above relations.

###  II.1 Renormalization in vacuum

We first set $\mu_{I}=0$ and fix the counter terms to define the model. In
vacuum the $\sigma$ fields condense, and we shift fields
$\sigma=\tilde{\sigma}+\langle\sigma\rangle$. The meson masses are

$\displaystyle
m_{\sigma}^{2}=m_{0}^{2}+\frac{\lambda}{\,2\,}\langle\sigma\rangle^{2}\,,~{}~{}~{}~{}m_{\pi}^{2}=m_{0}^{2}+\frac{\lambda}{\,6\,}\langle\sigma\rangle^{2}\,,$
(28)

and quarks acquire the effective quark mass

$\displaystyle M_{q}=g\langle\sigma\rangle\,.$ (29)

In the on-shell renormalization we demand that the above relations hold at any
orders of loop computations.

Below we eliminate $\langle\sigma\rangle$ in favor of $M_{q}$, and omit the
tilde of the quantum fields. We rewrite the Lagrangian

$\displaystyle\mathcal{L}=\mathcal{L}_{0}+\mathcal{L}_{2}+\mathcal{L}_{\rm
int}\,,$ (30)

whose explicit expressions are given shortly. Here the Lagrangian linear in
quantum fields is dropped off since the contributions will be cancelled when
we substitute the solutions for the gap equations. The first term
$\mathcal{L}_{0}=-V_{0}$ is

$\displaystyle V_{0}=\frac{\,m_{0}^{2}\,}{2g^{2}}M_{q}^{2}+{\lambda\over
24g^{4}}M_{q}^{4}-\frac{\,h\,}{g}M_{q}\,.$ (31)

Without the explicit breaking $h$, nonzero chiral condensates emerge only for
$m_{0\sigma}^{2}<0$.

The quantum corrections are computed with the unperturbed Lagrangian
$\mathcal{L}_{2}$ including the effective masses,

$\displaystyle\mathcal{L}_{2}$ $\displaystyle=$ $\displaystyle{1\over
2}\left[(\partial_{\mu}\sigma)^{2}+(\partial_{\mu}\vec{\pi})^{2}\right]-{m_{\sigma}^{2}\over
2}\sigma^{2}-{m_{\pi}^{2}\over 2}\vec{\pi}^{2}$ (34)
$\displaystyle+\overline{\psi}\big{(}\,\mathrm{i}{\ooalign{\hfil/\hfil\crcr$\partial$}}-M_{q}\,\big{)}\psi\,,$

and the interaction terms

$\displaystyle\mathcal{L}_{\rm int}=-{\lambda\over
24}(\vec{\phi}^{2})^{2}-{\lambda M_{q}\over
6g}\sigma\vec{\phi}^{2}-g\overline{\psi}(\sigma+\mathrm{i}\gamma^{5}\vec{\tau}\cdot\vec{\pi})\psi\,,$
(35)

which are used to compute loop corrections.

Now we consider the on-shell renormalization at one-loop. After shifting
fields $\sigma\rightarrow\sigma+M_{q}/g$, the counter terms of the quadratic
fields appears as

$\displaystyle\frac{\,1\,}{2}\bigg{[}\,\delta
m_{0\sigma}^{2}+\frac{\,\delta\lambda_{4\sigma}\,}{2}\bigg{(}\frac{M_{q}}{g}\bigg{)}^{2}\,\bigg{]}\sigma^{2}\,,$
$\displaystyle\frac{\,1\,}{2}\bigg{[}\,\delta
m_{0\pi}^{2}+\frac{\,\delta\lambda_{2\sigma\pi}\,}{6}\bigg{(}\frac{M_{q}}{g}\bigg{)}^{2}\,\bigg{]}\vec{\pi}^{2}\,,$
(36)

Using these vertices and the counter terms, the on-shell renormalization
condition at one-loop is expressed for the meson part as ($\phi=(\sigma,\pi)$)

$\displaystyle-\mathrm{i}\frac{\,\partial\Sigma_{\phi}\,}{\partial
p^{2}}\bigg{|}_{p^{2}=m_{\phi}^{2}}+\delta Z_{\phi}=0\,,$ (37)

and the mass part

$\displaystyle-\mathrm{i}\Sigma_{\sigma}(m_{\sigma}^{2})+\delta
Z_{\sigma}m_{\sigma}^{2}-\delta
m_{0\sigma}^{2}-\frac{\,\delta\lambda_{4\sigma}\,}{2}\bigg{(}\frac{M_{q}}{g}\bigg{)}^{2}=0\,,$
$\displaystyle-\mathrm{i}\Sigma_{\pi}(m_{\pi}^{2})+\delta
Z_{\pi}m_{\pi}^{2}-\delta
m_{0\pi}^{2}-\frac{\,\delta\lambda_{2\sigma\pi}\,}{6}\bigg{(}\frac{M_{q}}{g}\bigg{)}^{2}=0\,.$
(38)

It turns out that the $-\mathrm{i}\Sigma_{\phi}(m_{\phi}^{2})+\delta
Z_{\phi}m_{\phi}^{2}$ terms leave the divergent terms proportional to
$M_{q}^{2}$ which are canceled by $\delta\lambda_{4\sigma}M_{q}^{2}$ and
$\delta\lambda_{2\sigma\pi}M_{q}^{2}$ counter terms. Meanwhile the mass
counter term does not contain the divergence,

$\displaystyle\delta_{\overline{\rm MS}}m_{0\phi}^{2}=0\,.$ (39)

Since we impose the on-shell renormalization conditions on $\pi$ and $\sigma$,
we have conditions to fix four independent counter terms at this stage. For
example one can determine $\delta\lambda_{4\sigma}$ and $\delta
m_{0\sigma}^{2}$ first, and use the relations Eqs.(25) and (27) to determine
$\delta m_{0\pi}^{2}$ and $\delta\lambda_{2\sigma\pi}$. The divergent parts
are shown at the end of this section. For the finite parts, instead of
displaying each term, in the Appendix. B we give the full expression of the
effective potential in which various terms are sorted out.

The rest of counter terms is $\delta g_{\phi}$. For the quark vector and
scalar self-energies,
$\Sigma_{\psi}={\ooalign{\hfil/\hfil\crcr$p$}}\Sigma_{V}(p^{2})-\Sigma_{S}(p^{2})$,

$\displaystyle\Sigma_{V}(M_{q}^{2})+\delta Z_{\psi}=0\,,$
$\displaystyle\Sigma_{M}(M_{q}^{2})+\delta g_{\sigma}\frac{M_{q}}{g}=0\,.$
(40)

Actually, in the present model, the quark self-energies come only from meson
loops which are suppressed by $1/N_{\rm c}$. Hence in the large $N_{\rm c}$
limit one can take $\delta Z_{\psi},\delta g_{\sigma}\rightarrow 0$.

All the counter terms can be now fixed. The divergent $\overline{\rm MS}$
counter terms are

$\displaystyle\delta_{\overline{\rm MS}}Z_{\phi}=-\frac{\,4g^{2}N_{\rm
c}\,}{\,(4\pi)^{2}\epsilon\,}\,,~{}~{}~{}\delta_{\overline{\rm
MS}}m_{0\phi}^{2}=0\,,~{}~{}~{}\delta_{\overline{\rm MS}}g^{2}=0\,,$
$\displaystyle~{}~{}~{}\delta_{\overline{\rm
MS}}\lambda_{4\sigma}=\delta_{\overline{\rm
MS}}\lambda_{4\pi}=\delta_{\overline{\rm
MS}}\lambda_{2\sigma\pi}=-\frac{\,48g^{4}N_{\rm
c}\,}{\,(4\pi)^{2}\epsilon\,}\,.$ (41)

Here we postpone presenting lengthy finite terms to Appendix.

Finally we write the full expression of the effective potential. The one-loop
contribution entirely comes from the quark zero point energy and the counter
terms coupled to the condensates. The former is divergent ($N_{\rm f}=2$),

$\displaystyle V_{q}=-2N_{c}N_{\rm
f}\int_{{\bm{p}}}E_{D}({\bm{p}})\,,~{}~{}~{}E_{D}=\sqrt{{\bm{p}}^{2}+M_{q}^{2}}\,.$
(42)

We will treat these divergences by dimensional regularization $d\to
3-2\epsilon$,

$\displaystyle\int_{{\bm{p}}}=\left({e^{\gamma_{E}}\Lambda^{2}\over
4\pi}\right)^{\epsilon}\int{d^{d}p\over(2\pi)^{d}}.$ (43)

where $\Lambda$ is the renormalizing scale introduced by the $\overline{\rm
MS}$ scheme and $\gamma_{E}=0.577...$ is the Euler-Mascheroni constant. The
appearance of the scale parameter $\Lambda$ is just fictitious and can be
eliminated in favor of physical parameters such as $f_{\pi}$. The final
expression of the effective potential does not manifestly contain $\Lambda$.
After calculating the $d$-dimensional momentum integral we obtain

$\displaystyle V_{q}$ $\displaystyle=$
$\displaystyle{4N_{c}\over(4\pi)^{2}}\left({e^{\gamma_{E}}\Lambda^{2}\over
M_{q}^{2}}\right)^{\epsilon}\Gamma(-2+\epsilon)\,M_{q}^{4}$ (44)
$\displaystyle=$
$\displaystyle{2N_{c}\over(4\pi)^{2}}\bigg{(}\frac{1}{\epsilon}+\frac{3}{2}+2\ln\frac{\Lambda}{M_{q}}\bigg{)}M_{q}^{4}+O(\epsilon)\,,$

where we used

$\displaystyle\Gamma(-2+\epsilon)$ $\displaystyle={1\over
2}\left[{1\over\epsilon}-\gamma_{E}+{3\over 2}+{\cal O}(\epsilon)\right]\,.$
(45)

The divergence can be canceled with the terms $\propto M_{q}^{4}$ from the
counter Lagrangian,222To reassure the cancellation, it is useful to expand the
single particle energy in powers of condensates. The quadratic terms
correspond to the meson self-energy graphs at zero external momenta whose
divergences have been cancelled in renormalization of the $Z$-factors.
Similarly the quartic terms correspond to graphs in which four meson external
lines coupled to single quark loop. The divergence has been cancelled by
renormalization of $\lambda$. Higher order terms do not contain the UV
divergences, but include the IR divergences. The IR divergences are being
handled by resummation of graphs each of which is IR divergent but the sum is
harmless.  (no divergent counter terms $\propto M_{q}^{2}$ as we already found
$\delta_{\overline{\rm MS}}m_{0\phi}^{2}=0$)

$\displaystyle\delta_{\overline{\rm MS}}V_{q}^{\rm
c.t.}={\,\delta_{\overline{\rm MS}}\lambda_{4\sigma}\,\over
24g^{4}}M_{q}^{4}\,,$ (46)

which can be read off from $\mathcal{L}_{\rm c.t.}$ in Eq.(22) after
substituting $\sigma=\tilde{\sigma}+M_{q}/g$. Combining these, our
renormalized single particle energy density is

$\displaystyle V_{q}^{R}=V_{q}+\delta_{\overline{\rm MS}}V_{q}^{\rm c.t.}\,.$
(47)

Now our one-loop effective potential consists of the tree potential, single
quark contributions at one-loop, and the potential from the counter
Lagrangian,

$\displaystyle V_{\rm eff}=V_{0}+V_{q}^{R}+{\,\delta_{\rm
fin}m_{0\sigma}^{2}\,\over 2g^{2}}M_{q}^{2}+{\,\delta_{\rm
fin}\lambda_{4\sigma}\,\over 24g^{4}}M_{q}^{4}\,.$ (48)

The renormalized single particle energy density with the renormalization scale
$\Lambda$ is

$\displaystyle
V_{q}^{R}={2N_{c}\over(4\pi)^{2}}\bigg{(}\frac{3}{2}+2\ln\frac{\Lambda}{M_{q}}\bigg{)}M_{q}^{4}\,.$
(49)

The overall scale $\Lambda$ associated loops is related to the physical
parameter through

$\displaystyle\frac{\,\partial V_{q}^{R}\,}{\partial
M_{q}}\bigg{|}_{M_{q}=gf_{\pi}}+{\,\delta_{\rm fin}m_{0\sigma}^{2}\,\over
g^{2}}M_{q}+{\,\delta_{\rm fin}\lambda_{4\sigma}\,\over 6g^{4}}M_{q}^{3}=0\,,$
(50)

where we have used $\partial V_{0}/\partial M_{q}=0$ at $M_{q}=gf_{\pi}$. Then
$\Lambda$ can be eliminated.

###  II.2 At finite isospin density

For a large isospin chemical potential, either $\pi_{1}$ or $\pi_{2}$ can
condense while $\pi_{3}$ fields are unaffected. Without loss of generality we
assume the $\pi_{1}$ to condense. We shift fields

$\displaystyle\sigma=\tilde{\sigma}+\langle\sigma\rangle\,,~{}~{}~{}~{}~{}\pi_{1}=\tilde{\pi}_{1}+\langle\pi_{1}\rangle\,.$
(51)

Quarks acquire the chiral effective masses and pion induced gaps,

$\displaystyle
M_{q}=g\langle\sigma\rangle\,,~{}~{}~{}~{}~{}\Delta=g\langle\pi_{1}\rangle\,.$
(52)

As in the vacuum case, we eliminate $\langle\sigma\rangle$ and
$\langle\pi_{1}\rangle$ in favor of $M_{q}$ and $\Delta$, omit the tilde of
the quantum fields, and rewrite the Lagrangian

$\displaystyle\mathcal{L}=\mathcal{L}_{0}+\mathcal{L}_{2}+\mathcal{L}_{\rm
int}\,.$ (53)

The analyses go in the similar way as in the vacuum case. We do not repeat
writing the Lagrangian for the shifted variables, but merely point out
important differences from the vacuum case.

(i) The tree level potential $\mathcal{L}_{0}=-V_{0}$ takes the form

$\displaystyle V_{0}$ $\displaystyle=$
$\displaystyle\frac{\,m_{0}^{2}\,}{2g^{2}}M_{q}^{2}+\frac{\,m_{0}^{2}-4\mu_{I}^{2}\,}{2g^{2}}\Delta^{2}$
(54) $\displaystyle+{\lambda\over
24g^{4}}\big{(}M_{q}^{2}+\Delta^{2}\big{)}^{2}-\frac{\,h\,}{g}M_{q}\,.$

It is clear that a large $\mu_{I}$ reduces the energy. Eventually the pion
condensates emerge when the coefficient of the $\Delta^{2}$ term flips the
sign. At low density, the coefficient of $\Delta^{2}$ term is proportional to
$(m_{\pi}^{\rm vac})^{2}-4\mu_{I}^{2}$. Hence the pion condensates start to
form at $\mu_{I}=m_{\pi}^{\rm vac}/2$. This holds with loop corrections; at
the onset, we need only vacuum loop corrections with $\Delta=0$ which vanish
at the on-mass-shell in our renormalization scheme. After the $\Delta$ becomes
nonzero, our tree level extraction of the meson masses is no longer valid.

(ii) The quark part in the unperturbed Lagrangian $\mathcal{L}_{2}$ acquires
an extra term

$\displaystyle\mathcal{L}_{2}^{\rm
extra}=\overline{\psi}\mathrm{i}\gamma^{5}\tau_{1}\Delta\psi\,,$ (55)

with which the quark propagator becomes the BCS type propagator. The poles
exist at

$\displaystyle
E_{u}=E_{\overline{d}}=E(\mu_{I})\,,~{}~{}~{}E_{d}=E_{\overline{u}}=E(-\mu_{I})\,,$
(56)

where (see the derivation in appendix A)

$\displaystyle E(\mu_{I})=\sqrt{\left(E_{D}-\mu_{I}\right)^{2}+\Delta^{2}}\,.$
(57)

The $u$ and $\bar{d}$ quarks excite with the minimum energy of the BCS gap
$\sim\Delta$ near the Fermi surface. Meanwhile $d$ and $\bar{u}$ quarks need
large energies of $\sim M+\mu_{I}$ to get excited.

(iii) The single particle energies depend on $M_{q}$ and $\Delta$ in medium,

$\displaystyle
V_{q}=-N_{c}\int_{{\bm{p}}}\big{(}E_{u}+E_{d}+E_{\overline{u}}+E_{\overline{d}}\big{)}\,,$
(58)

whose divergent part can be extracted as

$\displaystyle V_{q}^{\rm div}=-N_{\rm
c}\int_{{\bm{p}}}\left[\sqrt{E_{D}^{2}+\Delta^{2}}+{\mu_{I}^{2}\Delta^{2}\over
2(E_{D}^{2}+\Delta^{2})^{3/2}}\right]\,,$ (59)

where we expand $V_{q}$ by $\mu_{I}$. The integral is evaluated by the
dimensional regularization as

$\displaystyle V_{q}^{\rm div.}=$ $\displaystyle\
{4N_{c}\over(4\pi)^{2}}\left({e^{\gamma_{E}}\Lambda^{2}\over
M_{q}^{2}+\Delta^{2}}\right)^{\epsilon}$
$\displaystyle\times\left[(M_{q}^{2}+\Delta^{2})^{2}\Gamma(-2+\epsilon)-2\mu_{I}^{2}\Delta^{2}\Gamma(\epsilon)\right]\,.$
(60)

The condensate dependent divergences are cancelled by the $\overline{\rm MS}$
counter terms coupled to the condensates,

$\displaystyle\delta_{\overline{\rm MS}}V_{q}^{\rm c.t.}$ $\displaystyle=$
$\displaystyle\frac{\,\delta_{\overline{\rm
MS}}Z_{\pi}\,}{\,2g^{2}\,}\big{(}-4\mu_{I}^{2}\Delta^{2}\,\big{)}$ (61)
$\displaystyle+{1\over
24g^{4}}\big{[}\delta\lambda_{4\sigma}M_{q}^{4}+\delta\lambda_{4\pi}\Delta^{4}+2\delta\lambda_{2\sigma\pi}M_{q}^{2}\Delta^{2}\big{]}_{\overline{\rm
MS}}\,,$

which can be read off from $\mathcal{L}_{\rm c.t.}$ in Eq.(22) after
substituting $\sigma=\tilde{\sigma}+M_{q}/g$ and
$\pi_{1}=\tilde{\pi}_{1}+\Delta/g$.

Combining these, our renormalized single particle energy density is

$\displaystyle V_{q}^{R}=V_{q}+\delta_{\overline{\rm MS}}V_{q}^{\rm c.t.}\,.$
(62)

In practice it is useful to use the expression

$\displaystyle V_{q}^{R}=\big{(}V_{q}-V_{q}^{\rm div}\big{)}+\big{(}V_{q}^{\rm
div}+\delta_{\overline{\rm MS}}V_{q}^{\rm c.t.}\big{)}\,.$ (63)

The first two terms are convergent momentum integral and computed numerically.
They yield the $\mu_{I}^{4}$ terms which become dominant in the high density
limit. Meanwhile the divergences in the last two terms are handled
analytically.

Summing the tree potential $V_{0}$ in Eq.(68), the renormalized single quark
loop potential $V_{q}^{R}$ in Eq.(62), and the potential made of the finite
counter terms coupled to the condensates, our one-loop effective potential is
expressed as333 The effective potential in this work slightly improves the one
in Ref.[52] in which, although not mentioned explicitly, the massless limit
$m_{\pi}\rightarrow 0$ seems to be used in evaluation of the counter terms.

$\displaystyle V_{\rm eff}(M_{q},\Delta)=V_{0}+V_{q}^{R}$
$\displaystyle~{}~{}~{}+\frac{\,\delta_{\rm
f}Z_{\pi}\,}{\,2g^{2}\,}\big{(}-4\mu_{I}^{2}\Delta^{2}\,\big{)}+{1\over\,2g^{2}\,}\big{(}\delta_{\rm
f}m_{0\sigma}^{2}M_{q}^{2}+\delta_{\rm f}m_{0\pi}^{2}\Delta^{2}\big{)}$
$\displaystyle~{}~{}~{}+{1\over
24g^{4}}\big{[}\delta\lambda_{4\sigma}M_{q}^{4}+\delta\lambda_{4\pi}\Delta^{4}+2\delta\lambda_{2\sigma\pi}M_{q}^{2}\Delta^{2}\big{]}_{\rm
f}\,.$ (64)

The expectation value $M_{q*}$ and $\Delta_{*}$ are determined by the gap
equations,

$\displaystyle\left.\frac{\partial V_{\rm eff}}{\partial
M_{q}}\right|_{M_{q*},\Delta_{*}}=0\,,~{}~{}~{}~{}\left.\frac{\partial V_{\rm
eff}}{\partial\Delta}\right|_{M_{q*},\Delta_{*}}=0\,.$ (65)

In the next section we examine the behaviors of condensates and the relation
to the thermodynamics.

##  III Equations of state

We now numerically examine the mean field EOS. We fix the model parameters to
reproduce the following quantities in vacuum:

$\displaystyle m_{\pi}$ $\displaystyle=$ $\displaystyle 140\,{\rm
MeV}\,,~{}~{}~{}~{}~{}m_{\sigma}=600\,{\rm MeV}\,,$ $\displaystyle f_{\pi}$
$\displaystyle=$ $\displaystyle 90\,{\rm
MeV}\,,~{}~{}~{}~{}~{}~{}\,M_{q}=300\,{\rm MeV}\,.$ (66)

which correspond to the following coupling constants,

$\displaystyle g\simeq 3.333\,,~{}~{}~{}~{}\lambda\simeq 126.1\,.$ (67)

Figure 1:  Chiral and pion condensates as functions of a scaled chemical
potential $\mu_{I}/m_{\pi}$.

###  III.1 Evolution of microscopic quantities

Figure 2:  Top: A meson propagator as an elementary particle. Bottom: A meson
propagator with a quark loop. Now the meson can be interpreted as a composite
particle.

First we examine how condensates evolve as functions of $\mu_{I}$. Shown in
Fig.1 are the constituent quark mass and the gap associated with the pion
condensate. There are no significant differences in $M_{q}$ between the tree-
level (dashed blue line) and one-loop (solid blue line) results. Meanwhile,
the pion condensate $\Delta$ increases linearly with chemical potential at the
tree level, whereas it converges to a finite value, $\Delta\simeq 300\,{\rm
MeV}$, at one-loop. This drastic change of behavior shows that the one-loop
correction has more contents than mere perturbative corrections.

At tree-level, the mesons are treated as elementary particles (Fig.2 (Top)).
By adding one-loop corrections, however, they no longer can be regarded as
elementary particles but must be regarded as composite particles made of
quarks and antiquarks (Fig.2 (Bottom)). As a result, the meson condensates can
behave very differently at tree level and at one-loop.

Now we closely examine the chemical potential dependence of the pion
condensates. Consider the large chemical potential limit $\mu_{I}\to\infty$.
At tree level, we suppose $M_{q}\rightarrow 0$ in the large $\mu_{I}$ limit.
Then

$\displaystyle
V_{0}~{}\rightarrow~{}-\frac{\,2\mu_{I}^{2}\,}{g^{2}}\Delta^{2}+{\lambda\over
24g^{4}}\Delta^{4}$ (68)

whose variation with respect to $\Delta$ readily leads to the scaling
$\Delta_{\rm tree}\sim\mu_{I}$. This tree level scaling may be understood as
generic; the only relevant scale at high density is $\mu_{I}$ so that
$\Delta\sim\mu_{I}$ is inevitable. Accordingly the potential satisfies the
conformal behavior $P\sim\mu_{I}^{4}$ by the $\mu_{I}^{2}\Delta^{2}$ and
$\Delta^{4}$ terms. But descriptions with large magnitudes of $\Delta$ always
leave questions on the impact of higher order terms in $\Delta$, e.g.,
$\Delta^{6}$ terms. In the present quark-meson model, such higher order terms
are discarded by demanding the renormalizability of models, but the
renormalizability should matter only when we consider loop corrections, not at
tree level discussions. We conclude that the tree level description with
$\Delta\sim\mu_{I}$ does not have solid ground and is potentially misleading.

The above qualitative feature drastically changes at one-loop and several
unwanted behaviors at tree level disappear. Mesonic fields acquire the status
of composite particles, and the size of pion condensates (or more precisely
the gap) no longer keeps growing but converges into a finite value. In the
large density limit,

$\displaystyle V_{\rm
eff}(M_{q},\Delta)~{}\rightarrow~{}V_{q}^{R}-\frac{\,2\mu_{I}^{2}\,}{g^{2}}\big{(}1+\delta_{\rm
f}Z_{\pi}\big{)}\Delta^{2}$ $\displaystyle\hskip 71.13188pt+{\delta_{\rm
f}m_{0\pi}^{2}\over\,2g^{2}\,}\Delta^{2}+{\lambda+\delta_{\rm
f}\lambda_{4\pi}\over 24g^{4}}\Delta^{4}\,.$ (69)

Through the counter terms for renormalization, the mass scale other than
$\mu_{I}$ can couple to $\Delta$. We no longer need $\Delta\sim\mu_{I}$; the
gap equation is solved by balancing $\Delta$ and the scale in the counter
terms and $V_{q}^{R}$. The dominant $\mu_{I}^{4}$ term for the pressure comes
from the single quark energy density $V_{q}^{R}$ where the kinetic energy
dominates over potentials, $E_{\rm kin}\sim\mu_{I}\gg\Delta$.

###  III.2 Equations of state

We express equations of state by the thermodynamic pressure as a function of
$\mu_{I}/m_{\pi}$. Using the effective potential at extrema with respect to
variations of $M_{q}$ and $\Delta$, the pressure is expressed as

$\displaystyle P(\mu_{I})=-V_{\rm
eff}\big{(}M_{q*}(\mu_{I}),\Delta_{*}(\mu_{I});\mu_{I}\big{)}\,.$ (70)

The isospin and energy densities are given by

$\displaystyle n_{I}=\frac{\,\partial
P\,}{\,\partial\mu_{I}\,}\,,~{}~{}~{}\varepsilon=\mu_{I}n_{I}-P\,.$ (71)

We study the sound velocity

$\displaystyle c_{s}^{2}=\frac{\,\partial
P\,}{\,\partial\varepsilon\,}=\frac{\,n_{I}\,}{\,\mu_{I}\chi_{I}\,}\,,~{}~{}~{}~{}~{}\chi_{I}=\frac{\,\partial^{2}P\,}{\,\partial\mu_{I}^{2}\,}\,.$
(72)

where $\chi_{I}$ is the isospin susceptibility.

In the following we compare our results with the lattice data in Ref.[60]. The
setup of the lattice is $N_{\rm f}=2+1$ flavors of rooted staggered quarks
with the quark masses at the physical point. The pion decays constant is
$f_{\pi}\simeq 92$-96 MeV for the lattice spacing explored (the definition of
$f_{\pi}$ differs by a factor $\sqrt{2}$ from ours). It should be noted that
their results at $T=0$ are obtained by correcting the data at small but finite
$T$ using the ChEFT.

Figure 3 shows the isospin density $n_{I}$ as a function of the isospin
chemical potential $\mu_{I}$, and Fig. 4 the pressure as a function of energy
density. As for the $n_{I}$, our results for
$(m_{\sigma},f_{\pi})=(600,90),(800,90),(600,100)$, and $(800,100)$ MeV are
all in good agreement with the lattice results in Ref.[60] up to
$\mu_{I}\lesssim m_{\pi}$. Beyond $\mu_{I}\gtrsim m_{\pi}$ or $n_{I}\gtrsim
0.5n_{0}$ the lattice data is not available. We note that the inclusion of
one-loop or, more physically, the quark substructure tempers the rapid growth
in $n_{I}$. As for the pressure at a given energy density, the quark
substructure effects enhance the pressure. We also note that the pressure is
reduced for larger $m_{\sigma}$ and $f_{\pi}$. In other words, with stronger
chiral symmetry breaking in vacuum (which increases both $m_{\sigma}$ and
$f_{\pi}$), the high density EOS after the chiral restoration becomes softer.
This point is examined in Sec.IV.1.

Figure 3:  Isospin density $n_{I}$ as a function of $\mu_{I}$ calculated from
the one-loop effective potential. The curve with a band in a small frame is
the lattice QCD result (available up to $\mu_{I}\simeq 0.9m_{\pi}$) from
Ref.[60].  Figure 4:  Pressure vs energy density from the one-loop effective
potential. The curve with a band is the lattice QCD result from Ref.[60].
Figure 5:  Squared sound velocity $c_{s}^{2}$ vs isospin density, compared
with the lattice QCD result from Ref.[60].

To further examine the variation of stiffness, we now turn to the behavior of
the sound velocity. Figure 5 shows the $c_{s}^{2}$ as a function of isospin
density $n_{I}$ normalized by the nuclear saturation density $n_{0}\sim
0.16\,{\rm fm}^{-3}$. The $c_{s}^{2}$ increases rapidly at low density, makes
a peak, and slowly relaxes to the conformal limit $1/3$ from above. This
qualitative feature seems robust and is consistent with the lattice result.
The quantitative details are different, however. The location of the
$c_{s}^{2}$ peak is near $n_{I}\sim n_{0}$ in our calculations for the
reasonable range of our parameter set for $m_{\sigma}$ and $f_{\pi}$, while
the lattice results in Ref.[60] indicate the peak at $n_{I}\sim 0.5n_{0}$. We
are not fully sure about the origin of the discrepancy. We have not found any
qualitative mechanisms to reconcile $\Delta\simeq\Lambda_{\rm QCD}$ with the
quick reduction of $c_{s}^{2}$ after making the peak. (Here we assume
$\Delta\simeq\Lambda_{\rm QCD}$ based on the lattice result for the melting
temperature of pion condensates, $T_{c}^{\rm lat}\simeq 170$ MeV, which seems
more or less constant to $\sim n_{0}$ or even higher densities, see
discussions around Eq.(90).) We rather expect that the lattice simulations
with finer lattice spacing and increasing sites would modify the EOS around
$\sim n_{0}$.

We note that the $c_{s}^{2}$ at tree level also shows the peak in the
crossover region and then the convergence to the conformal limit at high
density. As we noted in the Introduction, these tree level behaviors can be
achieved if the pressure behaves as
$P(\mu_{I})\sim\mu_{I}^{2}\Lambda^{2}+\mu_{I}^{4}$ where $\Lambda$ is some
non-perturbative scale. At low density we have $m_{\pi},m_{\sigma}\neq 0$
which plays the role of $\Lambda$, and at high density the tree level EOS is
dominated by $\mu_{I}^{2}\Delta^{2}$ and $\Delta^{4}$ terms with
$\Delta\sim\mu_{I}$. The low density behavior is qualitatively similar in
computations at tree level and at one-loop, but the mechanisms for high
density behaviors are different, as mentioned earlier. As a result, the
location of peaks, which are affected both by low and high density behaviors,
are naturally different. As claimed earlier in this section, the tree level
result which crucially depends on the scaling $\Delta\sim\mu_{I}$ is
potentially misleading from the microscopic point of view; we should not take
the tree level result at its face value.

With the above qualifications in mind, in the next section we look into more
details of our model regarding it as a model of composite particles.

###  III.3 Occupation Probability

At low density the effective degrees of freedom are pions and their internal
structure may be ignored. At higher density, the inter particle distance
becomes shorter and the quark substructure of pions becomes important.

To estimate where the quark substructure becomes important, we examine the
pion charge radius. It can be extracted from the vector form factor. The
experimental determination based on the $\pi e$ scattering and the
$e^{+}e^{-}\rightarrow\pi^{+}\pi^{-}$ process [61] yield the estimate $\langle
r^{2}\rangle_{V}=0.434(5)\,{\rm fm}^{2}$ [62], or

$\displaystyle r_{\pi}^{V}=\sqrt{\langle
r^{2}\rangle_{V}}~{}\simeq~{}0.66\,{\rm fm}\,,$ (73)

which has been well reproduced by lattice calculations [63, 64]. The typical
isospin density where pions overlap is estimated through

$\displaystyle n_{I}^{\rm overlap}=\big{(}4\pi r_{\pi}^{3}/3)^{-1}\simeq
0.83\,{\rm fm}^{-3}\simeq 5.2n_{0}\,.$ (74)

Figure 3 shows that the isospin chemical potential at $n_{I}^{\rm overlap}$ is
$\mu_{I}\simeq 3.8m_{\pi}\simeq 0.53\,{\rm GeV}$.

We note that this overlap density $n_{I}^{\rm overlap}$ is substantially
larger than the density $\sim n_{0}$ where the sound velocity peak is
developed. Hence there should be a more suitable measure to characterize the
location of sound velocity peak. Following Refs.[36, 65, 66], we interpret the
sound velocity peak as the signature of the quark saturation or the onset of
quark matter formation. The quark saturation is characterized by the
occupation probability of quark states; as density increases the low momentum
states are occupied with the probability $\simeq 1$. A newly added quark must
fill a state on top of the already occupied states.

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 A. The occupation probabilities for $u$-, $d$-, $\bar{u}$-,
$\bar{d}$-quarks are

$\displaystyle f(p)=f_{u,\bar{d}}(p)$ $\displaystyle={1\over
2}\left(1+{\mu_{I}-E_{D}\over E(\mu_{I})}\right)\,,$ (75)
$\displaystyle\bar{f}(p)=f_{\bar{u},d}(p)$ $\displaystyle={1\over
2}\left(1+{\mu_{I}+E_{D}\over E(\mu_{I})}\right)\,.$ (76)

Roughly speaking, $u$\- and $\bar{d}$-quarks occupy states up to
$\simeq\mu_{I}$ while $\bar{u}$\- and $d$-quarks are almost fully occupied as
in the Dirac sea without pion condensates.

Figure 6:  Occupation probability of $u$\- and $\bar{d}$ quarks which
corresponds to the residue of positive energy part of $\langle
u\overline{u}\rangle$. The plots are for each $0.25n_{0}$ with the gray line
from $n=0$ to $n=11n_{0}$, and the red lines are for each $2.0n_{0}$.

Shown in Fig.6 is the occupation probability $f(p)$ at various densities as a
function of quark momenta $p$. The densities we plotted are from $0$ to
$11n_{0}$ in $0.25n_{0}$ increments for gray curves and $2n_{0}$ increments
for red curves.

To relate the evolution of $f(p)$ to the stiffening of matter, it is useful to
decompose the evolution of $f(p)$ into two components. The first is the
“vertical evolution” in which $f(p)$ just increases its magnitude as $f(p)\sim
n_{I}\varphi_{\pi}^{\rm vac}(p)$ with $\varphi_{\pi}^{\rm vac}$ being the
quark distribution in a single pion (Fig.7, Left); this corresponds to the
regime where pions do not interact and quarks inside of pions are largely
unaffected. In this regime, $\varepsilon/n_{I}$ is close to a constant, and
therefore the pressure, $P=n_{I}^{2}\partial(\varepsilon/n_{I})/\partial
n_{I}$, is very small. While quarks can always contribute to the energy
density through the masses of pions, they do not directly contribute to the
pressure. The sound velocity is small in this regime. The second component is
the “horizontal evolution” in which the $f(p)$ increases in the high energy
components (Fig.7, Right). This is driven by both interactions and the Pauli
blocking effects. Here, $\varepsilon/n_{I}$ increases as in usual quark matter
and the pressure can be large. In reality with interactions, the evolution of
$f(p)$ is the mixture of these two components.

In our quark-meson model, Fig.6 suggests that, from $0$ to $\simeq n_{0}$, the
magnitude of $f(p)$ at $p=0$ grows rapidly from $0$ to $\simeq 0.9$, but at
higher density the distribution $f(p)$ develops toward the horizontal
direction. The significant turning from the vertical to horizontal evolution
occurs around $\sim n_{0}$ which roughly coincides with the location of the
$c_{s}^{2}$ peak. The low and high density limits of $c_{s}^{2}$ $\sim 0$ and
$\sim 1/3$, are easy to understood by the vertical and horizontal evolutions,
respectively.

Meanwhile, the quick growth of $c_{s}^{2}$ in the intermediate region,
$n_{I}\sim 0.3$-$1.0n_{0}$, requires us a closer look into how much the
horizontal evolution contributes. This point is elaborated in Fig.8 by taking
the ratio $f(p)/f(p)_{n_{I}=0.1n_{0}}$. The ratio should be flat during the
vertical evolution. During the horizontal evolution, the ratio is more
enhanced at higher momenta. To characterize the evolution of the Fermi
surface, it is useful to see a momentum, $p_{\rm surf}$, at which
$f(p)/\varphi_{\pi}^{\rm vac}(p)$ has the maximum (in practice we use
$f(p)/f(p)_{n_{I}=0.1n_{0}}$ regarding $f(p)_{n_{I}=0.1n_{0}}$ is within the
regime of the vertical evolution, $f(p)_{n_{I}=0.1n_{0}}\sim
0.1n_{0}\varphi_{\pi}^{\rm vac}(p)$). The evolution of $p_{\rm surf}$ can be
used to characterize transitions from the vertical to horizontal evolution;
$p_{\rm surf}$ stays almost constant up to $\simeq 0.3n_{0}$, and starts to
increase at higher density. This observation is consistent with substantial
$c_{s}^{2}\gtrsim 1/3$ at $n_{I}\gtrsim 0.3n_{0}$. This relatively early onset
of the horizontal evolution is likely due to small pion masses and
interactions among condensed pions. This point differs from QCD at finite
baryon density where baryons are heavy and $c_{s}^{2}\lesssim 0.1$ for
$n_{B}\sim 1$-$2n_{0}$ [67].

Figure 7:  Schematic figures for the evolution of the occupation probability.
Left: The “vertical” evolution; Right: The “horizontal” evolution.  Figure 8:
The ratio $f(p)/f(p)_{n_{I}=0.1n_{0}}$ for various $n_{I}$ from $0.1n_{0}$ to
$1.5n_{0}$ with $0.1n_{0}$ increments for gray curves and $0.5n_{0}$
increments for red curves. The dashed curve indicates the momentum $p_{\rm
surf}$ at which the ratio takes the maximum. Pure “vertical” evolution ceases
around $\sim 0.3n_{0}$ beyond which substantial “horizontal” evolution occurs
and increases $c_{s}^{2}$.

In our model, a quark matter with the substantial quark Fermi sea, with the
quark Fermi momenta of $\simeq\Lambda_{\rm QCD}\sim 200$-$300$ MeV or $\sim
2m_{\pi}$, is formed around $\simeq 4$-$6n_{0}$. Before reaching this domain,
the distribution $f(p)$ is not quite flat and there is little bulk domain; its
shape is largely characterized by the “surface” with the thickness
$\sim\Delta$, rather than the Fermi sea made of almost filled quark states
(Fig.7). This surface is the consequence of forming composite particles. The
above estimate of the density, $4$-$6n_{0}$, roughly coincides with the
$n_{I}^{\rm overlap}\simeq 5.2n_{0}$. In summary, the evolution of the
occupation probability can represent two characteristic scales; one is for the
quark saturation, and the other is for the overlap of composite particles. The
distinction of such two scales was emphasized in Ref.[68] which discriminate
the mode-by-mode percolation in momentum space from the conventional geometric
percolation.

##  IV Discussion

Here we address several issues not detailed in the previous sections. First we
discuss how the strength of chiral symmetry breaking in vacuum and its
restoration at high density affect EOS. For the high density domain, we
compare our results with pQCD at high density, and conjecture the importance
of the power corrections. Then we discuss the trace anomaly and the positivity
conjecture.

###  IV.1 Chiral symmetry restoration and softening

In Sec.III we have seen that larger $f_{\pi}$ and/or $m_{\sigma}$ lead to
softer EOS at high density. Here we try to explain this softening by focusing
on the chiral symmetry breaking in the vacuum and its restoration at high
density. In this context larger $f_{\pi}$ and $m_{\sigma}$ mean the stronger
chiral symmetry breaking in the QCD vacuum.

In the vacuum, the energy reduction due to the chiral symmetry breaking is
(Fig.9)

$\displaystyle B\equiv\Omega(M_{q}=m_{q})-\Omega(M_{q}=M_{q}^{\rm vac})\,,$
(77)

where the first term is the energy of the trivial vacuum while the second term
is the energy of the chiral symmetry broken vacuum. This sort of the energy
difference is often called the bag constant. Stronger breaking in the chiral
symmetry increases the size of the bag constant (Fig.10). Explicitly,
$B(m_{\sigma},f_{\pi})$ is given as

$\displaystyle B(600,90)$ $\displaystyle=$ $\displaystyle
0.389,~{}~{}~{}B(600,100)=0.483\,,$ $\displaystyle B(800,90)$ $\displaystyle=$
$\displaystyle 0.394,~{}~{}~{}B(800,100)=0.489\,,$ (78)

from which one can see that larger $f_{\pi}$ and $m_{\sigma}$ lead to a
greater $B$.

A larger bag constant softens EOS. To see this, it is useful to recall a bag
model with perturbative corrections. We note that the perturbative expansions
are performed around the trivial vacuum. Since our EOS is normalized to make
$P=0$ at $\mu=T=0$ for the non-perturbative vacuum, the perturbative
evaluation of EOS must be corrected by the non-perturbative normalization
constant. Then, the pressure and energy density are

$\displaystyle P_{\rm pert}^{\rm normalized}=P_{\rm
pert}-B,~{}~{}\varepsilon_{\rm pert}^{\rm normalized}=\varepsilon_{\rm
pert}+B.$ (79)

The bag constant associated with the chiral restoration reduces the pressure
and increases the energy density, resulting in a softer EoS. Similar
conclusions have been obtained in models with and without the $U(1)_{A}$
anomaly [69, 70].

Figure 9:  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.  Figure 10:  The energy density with
different strength of the chiral symmetry breaking. The bag constant is larger
for the stronger chiral symmetry breaking.

###  IV.2 Power corrections to pQCD at high density

Our quark-meson EOS predicts $c_{s}^{2}$ approaching $1/3$ from above as
density increases. This contradicts with the pQCD prediction in which
$c_{s}^{2}$ approaches $1/3$ from below. A possible origin of such discrepancy
would be the power corrections of $\sim\Lambda_{\rm QCD}^{2}\mu_{I}^{2}$ which
cannot be derived from perturbative computations.

In Introduction, we schematically showed how power corrections can enhance the
$c_{s}^{2}$ in Eqs.(1) and (2). The question is how large power corrections
should be to qualitatively change the perturbative behaviors of $c_{s}^{2}$.

For a given flavor $f$, the pQCD EOS up to $O(\alpha_{s}^{2})$ is given as
[71]

$\displaystyle P_{0}^{f}\\!$ $\displaystyle=$ $\displaystyle\\!\frac{N_{\rm
c}}{\,12\pi^{2}\,}\bigg{[}|\mu_{f}|u_{f}\bigg{(}\mu_{f}^{2}-\frac{\,5\,}{2}M_{q}^{2}\bigg{)}+\frac{\,3\,}{2}M_{q}^{4}\ln\frac{\,|\mu_{f}|+u_{f}\,}{M_{q}}\bigg{]}\,,$
$\displaystyle P_{1}^{f}\\!$ $\displaystyle=$
$\displaystyle\\!-\frac{\,\alpha_{s}N_{\rm
G}\,}{\,16\pi^{2}\,}\bigg{[}3\bigg{(}M_{q}^{2}\ln\frac{\,|\mu_{f}|+u_{f}\,}{M_{q}}-|\mu_{f}|u_{f}\bigg{)}^{2}-2u_{f}^{4}$
$\displaystyle+M_{q}^{2}\bigg{(}6\ln\frac{\,\Lambda_{\rm
reno}\,}{\,M_{q}\,}+4\bigg{)}\bigg{(}|\mu_{f}|u_{f}-M_{q}^{2}\ln\frac{\,|\mu_{f}|+u_{f}\,}{M_{q}}\bigg{)}\bigg{]}\,.$

In the isospin symmetric limit, $P^{u}=P^{d}$. Here, $P_{0}^{f}$ and
$P_{1}^{f}$ is zeroth and first order in $O(\alpha_{s})$, with
$u_{f}=\sqrt{\mu_{f}^{2}-M_{q}^{2}}$, $N_{\rm G}=N_{\rm c}^{2}-1$, and
$\Lambda_{\rm reno}$ being the renormalization scale. The running $\alpha_{s}$
is

$\displaystyle\alpha_{s}(\Lambda_{\rm reno})$ $\displaystyle=$
$\displaystyle\frac{\,4\pi\,}{\,\beta_{0}L\,}\bigg{[}1-2\frac{\,\beta_{1}\,}{\,\beta_{0}^{2}\,}\frac{\,\ln
L\,}{\,L\,}$ (81) $\displaystyle+$
$\displaystyle\frac{\,\beta_{1}^{2}\,}{\,\beta_{0}^{4}L^{2}\,}\big{(}\ln^{2}L-\ln
L-1+\frac{\,\beta_{2}\beta_{0}\,}{\,\beta_{1}^{2}\,}\big{)}$ $\displaystyle+$
$\displaystyle\frac{\,\beta_{1}^{3}\,}{\,\beta_{0}^{6}L^{3}\,}\big{(}-\ln^{3}L+\frac{\,5\,}{\,2\,}\ln^{2}L+2\ln
L$ $\displaystyle-$
$\displaystyle\frac{\,1\,}{\,2\,}-3\frac{\,\beta_{2}\beta_{0}\,}{\,\beta_{1}^{2}\,}\ln
L+\frac{\,\beta_{3}\beta_{0}^{2}\,}{\,2\beta_{1}^{3}\,}\big{)}\bigg{]}$

with $L=2\ln(\Lambda_{\rm reno}/\Lambda_{\overline{\rm MS}})$,
$\beta_{0}=11-2N_{\rm f}/3$, $\beta_{1}=102-19N_{\rm f}/3$,

$\displaystyle\beta_{2}\\!$ $\displaystyle=$
$\displaystyle\\!\frac{\,2857\,}{\,2\,}-\frac{\,5033\,}{\,18\,}N_{f}+\frac{\,325\,}{\,54\,}N_{f}^{2}\,,$
(82) $\displaystyle\beta_{3}\\!$ $\displaystyle=$
$\displaystyle\\!\bigg{(}\frac{\,149753\,}{\,6\,}+3564\,\xi(3)\bigg{)}$ (83)
$\displaystyle\\!-\bigg{(}\frac{\,1078361\,}{\,162\,}+\frac{\,6508\,}{\,27\,}\xi(3)\bigg{)}N_{f}$
$\displaystyle\\!+\bigg{(}\frac{\,50065\,}{\,162\,}+\frac{\,6472\,}{\,81\,}\xi(3)\bigg{)}N_{f}^{2}$
$\displaystyle\\!+\frac{\,1093\,}{\,729\,}N_{f}^{3}$

and $\Lambda_{\overline{\rm MS}}\simeq 340\,{\rm MeV}$. The central value of
$\Lambda_{\rm reno}$ is $\Lambda_{\rm reno}=\mu_{f}$, and as usual we vary
$\Lambda_{\rm reno}$ from $0.5\mu_{f}$ to $2\mu_{f}$.

Figure 11:  Running coupling from pQCD (dashed) and from the freezing coupling
(solid). For the renormalization scale $\Lambda_{\rm reno}=X\mu_{I}$, we
examine $X=0.5,1.0$, and $2.0$.

In addition to the perturbative running coupling which becomes unphysical
toward the Landau pole, we also examine the case with the freezing coupling in
the low energy limit. We divide the domain into three

$\displaystyle\alpha_{s}(Q^{2})$ $\displaystyle=$ $\displaystyle\alpha^{\rm
low}_{s}(Q^{2})\theta(t_{\rm low}-Q^{2})$ (84) $\displaystyle+$
$\displaystyle\alpha^{\rm mid}_{s}(Q^{2})\theta(t_{\rm
high}-Q^{2})\theta(Q^{2}-t_{\rm low})$ $\displaystyle+$
$\displaystyle\alpha^{\rm high}_{s}(Q^{2})\theta(Q^{2}-t_{\rm high})\,.$

where $t_{\rm low}^{1/2}=0.3\,{\rm GeV}$ and $t_{\rm high}^{1/2}=1.1$ GeV. For
the low energy limit we use the form suggested by Deur et al [72, 73]

$\displaystyle\alpha_{s}^{\rm low}=\alpha_{s}^{\rm
low}(0)\mathrm{e}^{-Q^{2}/4\kappa^{2}}\,,$ (85)

with $\alpha_{s}^{\rm low}(0)\simeq 1.22$ and $\kappa\simeq 0.51$. For the
high density we use the perturbative expression (81), and for the intermediate
region we use the interpolant

$\displaystyle\alpha_{s}^{\rm mid}(Q^{2})=\sum_{m=0}^{5}c_{m}\mu_{I}^{m}\,,$
(86)

where $c_{n}$’s are fixed by demanding the matching

$\displaystyle\frac{\,\partial^{n}\alpha_{s}^{\rm low/high}}{(\partial
Q^{2})^{n}}\bigg{|}_{Q^{2}=t_{\rm
low/high}}=\frac{\,\partial^{n}\alpha_{s}^{\rm mid}}{(\partial
Q^{2})^{n}}\bigg{|}_{Q^{2}=t_{\rm low/high}}\,.$ (87)

for $n=0,1,2$. The six boundary conditions fix the six $c_{n}$ uniquely.
Unlike in Refs.[72, 73] which needed only the continuity up to the first
derivative, in this work we use the interpolant not to generate any
discontinuities up to the second derivative, since we compute $c_{s}^{2}$.

Figure 12:  Perturbative pressure with perturbative running with the Landau
pole and infrared freezing coupling. The pressure is normalized by the
pressure in the Stefan-Boltzmann limit. Notations for the solid and dashed
lines are the same as Fig.11.

We set $Q^{2}=\Lambda_{\rm reno}^{2}$ and plot $\alpha_{s}(\Lambda_{\rm
reno})$ in Fig.11 together with the pQCD running coupling. With the IR
freezing coupling the artificial reduction of pQCD pressure is tempered and
the pressure remains positive toward the low density region (Fig.12).

Now we add power corrections which are parametrized in terms of gaps in the
pion condensed phase. The phase space factor $\sim 4\pi p_{F}^{2}\Delta$ times
the gap $\Delta$, divided by a factor $(2\pi)^{3}$, yields the naive estimate

$\displaystyle P_{\rm cond}=C\frac{\,\mu_{I}^{2}\Delta^{2}\,}{\,\pi^{2}\,}\,,$
(88)

where $C$ is a constant of $O(1)$.

In Son’s estimate [74], based on the color-magnetic long range forces, the gap
is evaluated as

$\displaystyle\Delta_{\rm color-
mag}=b|\mu_{I}|g_{s}^{-5}\mathrm{e}^{-3\pi^{2}/2g_{s}}$ (89)

with $g_{s}=g_{s}(|\mu_{I}|)$ being the running coupling constant. The gap can
be several hundreds MeV. Meanwhile, in our quark-meson model, we have found
$\Delta\simeq 300$ MeV. It is interesting to note that such $\Delta$ seems to
satisfy the BCS relation between the gap and the critical temperature,

$\displaystyle T_{c}^{\rm BCS}\simeq 0.57\,\Delta_{\rm BCS}\,,$ (90)

which implies $T_{c}\simeq 171\,{\rm MeV}$, in good agreement with the lattice
result $T_{c}^{\rm lattice}\simeq 160$-$170\,{\rm MeV}$ for the interval
$\mu_{I}\simeq 100$-$300\,{\rm MeV}$.

Figure 13:  Perturbative pressure plus power corrections divided by the
pressure in the Stefan-Boltzmann limit. The $\Delta=0,200$, and $400$ MeV. The
solid and dashed lines represent the freezing coupling and perturbative
running with the Landau pole, respectively. The $X$ for the $\Lambda_{\rm
reno}$ is fixed to $X=1$.  Figure 14:  Squared sound velocity $c_{s}^{2}$ for
pQCD+power corrections. We use the freezing coupling and vary $X$ from $0.5$
to $2$ to make bands. The solid lines in the bands are the $X=1$ case.

We set $C=1$ and examine pressure (Fig.13) and $c_{s}^{2}$ (Fig.14) for
$\Delta=0,200$, and $400$ MeV for pQCD running coupling (dashed) and freezing
coupling (solid). For $\Delta\simeq 200$ MeV, the power corrections are large
enough for $c_{s}^{2}$ to approach the conformal limit from above around $\sim
40n_{0}$. Meanwhile, at $\mu_{I}\sim 1$ GeV or $n_{I}\sim 40n_{0}$,
parametrically the power corrections in EOS are corrections of the order

$\displaystyle\bigg{(}\frac{\,\Delta\,}{\,\mu_{I}\,}\bigg{)}^{2}=0.09\bigg{(}\frac{\,\Delta\,}{\,300\,{\rm
MeV}\,}\bigg{)}^{2}\bigg{(}\frac{\,1\,{\rm
GeV}\,}{\,\mu_{I}\,}\bigg{)}^{2}\,,$ (91)

thus $\sim 10\%$ corrections. It is remarkable that even such small
corrections can change the qualitative behaviors of $c_{s}^{2}$ in the domain
where pQCD seems applicable.

###  IV.3 Trace Anomaly

Recently there has been growing interest in the trace anomaly in the context
of mechanical properties in hadrons [68, 75, 76, 77] and in neutron stars [47,
48, 49]. The latter is essentially the relation between $P$ vs $\varepsilon$
and is more fundamental than $c_{s}^{2}$ which includes only the information
of $\mathrm{d}P/\mathrm{d}\varepsilon$, not the overall magnitude of $P$. In
particular, Ref.[48] conjectured the trace anomaly to be positive. Below we
quickly mention the trace anomaly in dense matter and examine the positivity
conjecture by considering several non-perturbative effects.

The trace anomaly measures the breaking of the scale invariance and is given
by the expectation values of the operator

$\displaystyle\partial_{\mu}J^{\mu}_{D}=T^{\mu}_{\mu}=-\frac{\,\beta(g_{s})\,}{2g_{s}}G^{\mu\nu}G_{\mu\nu}+\sum_{f}m_{f}(1+\gamma_{m})\bar{q}_{f}q_{f}\,,$

where $J^{\mu}_{D}$ is the dilatation current, $\beta<0$ the QCD beta
function, and $\gamma_{m}$ the anomalous dimension of the quark mass.

For a hadronic state $|K\rangle$ with the momentum $K$, the energy momentum
tensor gives

$\displaystyle\langle K|T^{\mu\nu}(x)|K\rangle=K^{\mu}K^{\nu}/m_{H}\,,$ (93)

where the RHS is $x$-independent444 The state $|K\rangle$ with definite
momenta is a plane wave, meaning that the hadron can exist anywhere with the
probability $1/V_{\rm space}$.  and does not contain $g^{\mu\nu}$. The overall
$1/m_{H}$ factor is fixed by the condition at the rest frame,
$K_{R}^{\mu}=(m_{H},{\bf 0})$,

$\displaystyle\frac{\,\langle K_{R}|H|K_{R}\rangle\,}{\,\langle
K_{R}|K_{R}\rangle\,}=\frac{\,\langle
K_{R}|\int_{\bm{{x}}}T^{00}(x)|K_{R}\rangle\,}{\,\langle
K_{R}|K_{R}\rangle\,}=m_{H}\,.$ (94)

Thus, for a hadron at rest frame, we find

$\displaystyle\langle K_{R}|T^{\mu}_{\mu}(x)|K_{R}\rangle=\langle
K_{R}|T^{00}(x)|K_{R}\rangle=m_{H}\,,$ (95)

with vanishing spatial components, $\langle K_{R}|T^{ii}(x)|K_{R}\rangle=0$.
The trace anomaly is positive for a single hadron.

It is interesting to extend the above arguments to a many-body system. Unlike
the previous single particle case, not all particles stay at ${\bm{K}}=0$. For
instance an ideal Fermi gas leads to

$\displaystyle\langle
K_{1},\cdots|T^{ii}|K_{1},\cdots\rangle\sim\int_{{\bm{K}}}\frac{\,{\bm{K}}^{2}\,}{\,m_{H}\,}\,.$
(96)

After lowering one index, we get $\langle T^{i}_{i}\rangle<0$. Thus, the trace
anomaly in a many-body system can be negative in principle. In thermodynamic
systems, $\langle T^{\mu}_{\mu}\rangle=\varepsilon-3P$; the negative trace
anomaly means very large pressure, i.e., stiff EOS.

The trace anomaly characterizes the deviation from the relativistic or
conformal limit as expected at very high density. We first examine the impact
of the normalization in EOS. In the case of a bag model, we have

$\displaystyle\langle T^{\mu}_{\mu}\rangle_{\rm
bag}=\big{(}\varepsilon-3P\big{)}_{\rm pert}^{\rm
normalized}=\big{(}\varepsilon-3P\big{)}_{\rm pert}+4B\,.$ (97)

Changes from the non-perturbative to perturbative vacua enhances the trace
anomaly, supporting the positivity conjecture. Next we consider the impact of
power corrections using EOS similar to Eq.(1),

$\displaystyle P_{\rm with\,powers}=a_{0}\mu_{I}^{4}+a_{2}\mu_{I}^{2}\,,$ (98)

but now we include the running of coefficients $a_{0}$ and $a_{2}$ associated
with $\alpha_{s}(\mu_{I})$. The energy density can be computed as

$\displaystyle\varepsilon=3a_{0}\mu_{I}^{4}+2a_{2}\mu_{I}^{2}+\frac{\,\partial
a_{0}\,}{\,\partial\ln\mu_{I}\,}\mu_{I}^{4}+\frac{\,\partial
a_{2}\,}{\,\partial\ln\mu_{I}\,}\mu_{I}^{2}\,,$ (99)

and the trace anomaly is

$\displaystyle\langle T^{\mu}_{\mu}\rangle_{\rm
with\,powers}=-2a_{2}\mu_{I}^{2}+\frac{\,\partial
a_{0}\,}{\,\partial\ln\mu_{I}\,}\mu_{I}^{4}+\frac{\,\partial
a_{2}\,}{\,\partial\ln\mu_{I}\,}\mu_{I}^{2}\,.$ (100)

The running of $\alpha_{s}$ favors the positive trace anomaly, while the
attractive power corrections ($a_{2}>0$) favor the negative trace anomaly.

When we examine the trace anomaly, it is useful to divide it by
$3\varepsilon$,

$\displaystyle\Delta_{\rm tr}=\frac{1}{3}-\frac{P}{\varepsilon}\,.$ (101)

Shown in Fig.15 is the $\Delta_{\rm tr}$ as functions of $n_{I}$ for several
calculations, our quark meson model (QM) and pQCD results for the
renormalization scales with $X=0.5,1.0$, and $2.0$. For the pQCD, both the
perturbative (dashed) and IR freezing (solid) couplings are examined. Without
power corrections the $\Delta_{\rm tr}$ are all positive. In Fig.16, we fix
$X=1$ for these two couplings, and the dependence of $\Delta_{\rm tr}$ on the
power corrections. The $\Delta=0,200$, and $400$ MeV cases are shown. With
power corrections the $\Delta_{\rm tr}$ appears to be negative, as we
expected. Our QM model predicts $\Delta\simeq 300$ MeV and the negative trace
anomaly for wide range.

Since the pQCD corrections and bag constant favor the positive $\Delta_{\rm
tr}$, the negative $\Delta_{\rm tr}$ may be taken as an indicator of the
substantial power corrections that flip the sign.

Figure 15:  Trace anomaly $\Delta_{\rm tr}=1/3-P/\varepsilon$ as functions of
$n_{I}$ for our quark-meson model ($N_{\rm f}=2$ QM) and pQCD with the
perturbative running coupling (dashed) and IR freezing coupling (solid). We
vary $X$ from 0.5 to 2.0. The trace anomaly is all positive in pQCD but
negative for the quark meson model.  Figure 16:  Trace anomaly $\Delta_{\rm
tr}=1/3-P/\varepsilon$ as functions of $n_{I}$ for the quark-meson model, the
pQCD + powers with the perturbative running (dashed) and IR freezing coupling
(solid). The renormalization scale is fixed to $X=1$ while we examine
$\Delta=0,200,400$ MeV. For large $\Delta$ the trace anomaly can be negative.

##  V Summary

In this work we study the EOS of isospin QCD within a quark-meson model. The
model describes the BEC-BCS crossover of pion condensates. At tree level pions
look elementary, but at one-loop they acquire the status of composite
particles made of quarks and anti-quarks, resulting in behaviors qualitatively
very different from those at the tree level. The model is renormalizable and
we study its large density behaviors to study the impacts of non-perturbative
physics in the quark matter domain.

Our model exhibits the sound velocity peak with $c_{s}^{2}\simeq 0.7$ at
$n_{I}\simeq n_{0}$ substantially smaller than the density for pions to
spatially overlap $\sim 5n_{0}$. Regarding the BEC-BCS crossover as hadron-
quark crossover, we associate the sound velocity peak with the saturation of
quark states. After the bulk part of the quark Fermi surface is established,
the $c_{s}^{2}$ approaches $1/3$ as in the relativistic limit.

Our model shows that $c_{s}^{2}$ approaches 1/3 from above, mainly due to the
power corrections, $\sim\mu_{I}^{2}\Delta^{2}\sim\mu_{I}^{2}\Lambda_{\rm
QCD}^{2}$. This sort of terms is not available in pQCD calculations which
predict $c_{s}^{2}$ approaching $1/3$ from below. Which one dominates in
$c_{s}^{2}$ around $\sim 40n_{0}$ is a quantitative question. The existence of
the power corrections is related to the non-perturbative effects near the
quark Fermi surface and the structure of the QCD phase diagram. The question
is also related to the sign of the trace anomaly. The pQCD favors the positive
trace anomaly. If the trace anomaly appears to be negative, it is strong
indication of nontrivial Fermi surface structure.

The present work left several issues and should be extended to several
directions:

i) Our study should be extended to finite temperature (for recent discussions
on the quark contributions, see e.g., Refs.[78, 79]). Including thermal
effects into quark-meson models is straightforward, and the results are to be
compared with the lattice’s. Whether thermal excitations out of the quark
Fermi sea are confined or deconfined is an important issue in the context of
the quark-hadron-continuity. As for phenomenological applications to neutron
stars and heavy-ion-collisions, although several zero temperature EOS have
become available since 2012, finite temperature EOS with the continuity at the
level of excitations has not been constructed. For example, some difficulties
have been addressed for the nuclear-2SC continuity in Ref.[80]. The magnitude
of thermal corrections is much smaller than the cold matter part due to
$\sim(T/p_{F})^{2}$ suppression factors, but it can be important for NSs about
to collapse, e.g., those appearing in NS-NS mergers [81, 82, 83].

ii) The estimate of non-perturbative power corrections as well as the
normalization of EOS (bag constant) at high density should be improved.
Nowadays there has been increasing use of pQCD results to constrain the EOS at
$\simeq 5$-$40n_{0}$, with the help of general causality and thermodynamic
stability conditions (e.g., see Refs.[84, 85]). But as seen in our simple
exercise in Sec.IV.2, the power corrections of $\sim 10$% in the overall
magnitude can change the qualitative trend of quantities involving
derivatives. It should be important to see how the power corrections in
general affect the constraints at $\simeq 5$-$40n_{0}$.

iii) In this work we estimate the density for pion overlap based on the size
of pions in vacuum. But in medium pions may swell due to the quark exchange
among them. If the effective radii are larger than in the vacuum, the quark
saturation and the overlap of pions can take place at lower densities than the
estimates in this work. Changes in hadron size may occur already around
nuclear saturation density [86, 87], as indicated by the comparison of the
structure function for an isolated nucleon and nucleons in nuclei. It is
interesting to test these concepts in isospin QCD by comparing model
predictions with the lattice calculations.

###### Acknowledgements.

We thank Drs. Brandt and Endrodi for kindly providing us with their lattice
data in Ref.[60]. TK thanks Dr. Fujimoto for discussions on $c_{s}^{2}$ in
isospin QCD. The work of T.K. was supported by JSPS KAKENHI Grant No. 23K03377
and by the Graduate Program on Physics for the Universe (GPPU) at Tohoku
university.

##  Appendix A Quark propagators

We calculate the mean field quark propagator in the presence of the chiral and
pion condensates. From the propagator one can read off the excitation energy
from the pole of the propagator and the occupation probability from the
residue of the propagator.

It is convenient to introduce the projection operators for particle and
antiparticles,

$\displaystyle\Lambda_{p,a}={1\over 2}\pm{\gamma_{j}p^{j}+M_{q}\over
2E_{D}}\gamma_{0}\,,$ (102)

which satisfy

$\displaystyle\Lambda_{p}+\Lambda_{a}=1\,,~{}\Lambda_{p,a}\Lambda_{p,a}=\Lambda_{p,a}\,,~{}\Lambda_{p,a}\Lambda_{a,p}=0\,,$
(103)

as they should. The propagator of quarks can be written as

$\displaystyle{\mathrm{i}\over{\ooalign{\hfil/\hfil\crcr$p$}}+\mu_{f}\gamma_{0}-M_{q}}$
$\displaystyle=$ $\displaystyle
S_{p}(p)\gamma_{0}\Lambda_{p}+S_{a}(p)\gamma_{0}\Lambda_{a}\,,$ (106)
$\displaystyle S_{p,a}(p)$ $\displaystyle=$ $\displaystyle{\mathrm{i}\over
p_{0}+\mu_{f}\mp E_{D}}\,,$ (107)

where $S_{p,a}$ is the propagator for a particle and an antiparticle,
respectively.

The inverse of the propagator can also be separated by $\Lambda_{p,a}$ as

$\displaystyle{\ooalign{\hfil/\hfil\crcr$p$}}+\mu_{f}\gamma_{0}-M_{q}$
$\displaystyle=$ $\displaystyle(p_{0}+\mu_{f}-E_{D})\Lambda_{p}\gamma_{0}$
(111) $\displaystyle+(p_{0}+\mu_{f}+E_{D})\Lambda_{a}\gamma_{0}\,.$

The inverse of the Dirac operator is the quark propagator $S(p)$, and we can
write

$\displaystyle
S(p)^{-1}\equiv-\mathrm{i}\left(\begin{matrix}{\ooalign{\hfil/\hfil\crcr$p$}}+\mu_{u}\gamma_{0}-M_{q}&-i\gamma_{5}\Delta\\\
-i\gamma_{5}\Delta&{\ooalign{\hfil/\hfil\crcr$p$}}+\mu_{d}\gamma_{0}-M_{q}\end{matrix}\right)$
(112)

and consider its inverse. To simplify the discussion, we introduce the single-
particle propagator

$\displaystyle(G_{u,d}^{0})^{-1}=-\mathrm{i}\left({\ooalign{\hfil/\hfil\crcr$p$}}+\mu_{u,d}\gamma_{0}-M_{q}\right)\,,$
(115)

and write off-diagonal term $\Xi=\gamma_{5}\Delta$. Then our propagator must
satisfy

$\displaystyle S(p)^{-1}S(p)=\left(\begin{matrix}(G_{u}^{0})^{-1}&\Xi\\\
\Xi&(G_{d}^{0})^{-1}\end{matrix}\right)S(p)={\bf 1}.$ (116)

Introducing $(G_{u,d})^{-1}\equiv(G_{u,d}^{0})^{-1}-\Xi G_{d,u}^{0}\Xi$, the
propagator $S(p)$ can be written as follows.

$\displaystyle S(p)=\left(\begin{matrix}G_{u}&-G_{u}^{0}\Xi G_{d}\\\
-G_{d}^{0}\Xi G_{u}&G_{d}\end{matrix}\right).$ (117)

What we are interested in is the diagonal part of $S(p)$ which corresponds to
$\langle u\overline{u}\rangle$ and $\langle d\overline{d}\rangle$. Let us see
the detail of $G_{u,d}$. Its definition is

$\displaystyle(G_{u,d})^{-1}\equiv(G_{u,d}^{0})^{-1}-\Xi G_{d,u}^{0}\Xi$
$\displaystyle=-\mathrm{i}\left({\ooalign{\hfil/\hfil\crcr$p$}}+\mu_{u,d}\gamma_{0}-M_{q}\right)-\gamma_{5}\Delta{\mathrm{i}\over{\ooalign{\hfil/\hfil\crcr$p$}}+\mu_{d,u}\gamma_{0}-M_{q}}\gamma_{5}\Delta\,.$
(122)

To calculate the inverse we rewrite this formula using the projection
operators. Performing some calculations we can find

$\displaystyle\gamma_{5}\gamma_{0}\Lambda_{p,a}\gamma_{5}=-\Lambda_{a,p}\gamma_{0}.$
(124)

From the above, we obtain

$\displaystyle(G_{u,d})^{-1}$ $\displaystyle=$
$\displaystyle-\mathrm{i}\bigg{[}{\,p_{0}^{2}-(E_{D}-\mu_{u,d})^{2}-\Delta^{2}\over
p_{0}-\mu_{u,d}+E_{D}}\Lambda_{p}\gamma_{0}$ (125)
$\displaystyle~{}+\,{\,p_{0}^{2}-(E_{D}+\mu_{u,d})^{2}-\Delta^{2}\over
p_{0}-\mu_{u,d}-E_{D}}\Lambda_{a}\gamma_{0}\bigg{]}.$

Now we could separate the diagonal elements $(G_{u,d})^{-1}$ by projection
operators, and each part will not be mixed by the inverse operation.

Introducing the excitation energy

$\displaystyle\xi_{p,a}^{f}(p)=\sqrt{(E_{D}\mp\mu_{f})^{2}+\Delta^{2}},$ (126)

we obtain the propagator

$\displaystyle G_{f}$ $\displaystyle=$
$\displaystyle\mathrm{i}\bigg{[}{|u_{p}^{f}(p)|^{2}\over\,p_{0}-\xi_{p}^{f}(p)}+{|v_{p}^{f}(p)|^{2}\over\,p_{0}+\xi_{p}^{f}(p)}\bigg{]}\gamma_{0}\Lambda_{p}$
(127) $\displaystyle+$
$\displaystyle\mathrm{i}\bigg{[}{|u_{a}^{f}(p)|^{2}\over\,p_{0}-\xi_{a}^{f}(p)}+{|v_{a}^{f}(p)|^{2}\over\,p_{0}+\xi_{a}^{f}(p)}\bigg{]}\gamma_{0}\Lambda_{a}\,.$

Here we have used $\mu_{u,d}=-\mu_{d,u}$. The residues are

$\displaystyle|u_{p,a}^{f}(p)|^{2}$ $\displaystyle={1\over 2}\left(1+{\pm
E_{D}-\mu_{f}\over\xi_{p,a}^{f}}\right)$ (128)
$\displaystyle|v_{p,a}^{f}(p)|^{2}$ $\displaystyle={1\over 2}\left(1-{\pm
E_{D}-\mu_{f}\over\xi_{p,a}^{f}}\right).$ (129)

They correspond to the occupation probability and satisfy
$|u_{p}|^{2}+|v_{p}|^{2}=|u_{a}|^{2}+|v_{a}|^{2}=1$ as expected. In the main
text we use the expressions

$\displaystyle f=f_{u,\bar{d}}$ $\displaystyle={1\over
2}\left(1+{\mu_{I}-E_{D}\over
E(\mu_{I})}\right)=|v_{p}^{u}|^{2}=|u_{a}^{d}|^{2}\,,$ (130)
$\displaystyle\bar{f}=f_{\bar{u},d}$ $\displaystyle={1\over
2}\left(1+{\mu_{I}+E_{D}\over
E(\mu_{I})}\right)=|v_{a}^{u}|^{2}=|u_{p}^{d}|^{2}\,.$ (131)

##  Appendix B Full expression of the renormalized one-loop effective
potential

In the main text we express the one-loop effective potential using several
counter terms. Rewriting the counter terms using the physical parameters, the
final expression turns out to be

$\displaystyle V_{\mathrm{1-loop}}=$
$\displaystyle\frac{3}{4}m_{\pi}^{2}f_{\pi}^{2}\left[1-\frac{4{M_{q}^{\rm
vac}}^{2}N_{c}}{(4\pi)^{2}f_{\pi}^{2}}m_{\pi}^{2}F^{\prime}(m_{\pi}^{2})\right]{\Delta^{2}+M_{q}^{2}\over{M_{q}^{\rm
vac}}^{2}}$
$\displaystyle-\frac{1}{4}m_{\sigma}^{2}f_{\pi}^{2}\left[1+\frac{4{M_{q}^{\rm
vac}}^{2}N_{c}}{(4\pi)^{2}f_{\pi}^{2}}\left\\{\left(1-\frac{4{M_{q}^{\rm
vac}}^{2}}{m_{\sigma}^{2}}\right)F(m_{\sigma}^{2})+\frac{4{M_{q}^{\rm
vac}}^{2}}{m_{\sigma}^{2}}-F(m_{\pi}^{2})-m_{\pi}^{2}F^{\prime}(m_{\pi}^{2})\right\\}\right]{\Delta^{2}+M_{q}^{2}\over{M_{q}^{\rm
vac}}^{2}}$ $\displaystyle-2\mu_{I}^{2}f_{\pi}^{2}\left[1-\frac{4{M_{q}^{\rm
vac}}^{2}N_{c}}{(4\pi)^{2}f_{\pi}^{2}}\left\\{\ln\left({\Delta^{2}+M_{q}^{2}\over{M_{q}^{\rm
vac}}^{2}}\right)+F(m_{\pi}^{2})+m_{\pi}^{2}F^{\prime}(m_{\pi}^{2})\right\\}\right]{\Delta^{2}\over{M_{q}^{\rm
vac}}^{2}}$
$\displaystyle+\frac{1}{8}m_{\sigma}^{2}f_{\pi}^{2}\left[1-\frac{4{M_{q}^{\rm
vac}}^{2}N_{c}}{(4\pi)^{2}f_{\pi}^{2}}\left\\{\frac{4{M_{q}^{\rm
vac}}^{2}}{m_{\sigma}^{2}}\left(\ln\left({\Delta^{2}+M_{q}^{2}\over{M_{q}^{\rm
vac}}^{2}}\right)-\frac{3}{2}\right)\right.\right.$ $\displaystyle\hskip
85.35826pt\left.\left.-\left(1-\frac{4{M_{q}^{\rm
vac}}^{2}}{m_{\sigma}^{2}}\right)F(m_{\sigma}^{2})+F(m_{\pi}^{2})+m_{\pi}^{2}F^{\prime}(m_{\pi}^{2})\right\\}\right]{(\Delta^{2}+M_{q}^{2})^{2}\over{M_{q}^{\rm
vac}}^{4}}$
$\displaystyle-\frac{1}{8}m_{\pi}^{2}f_{\pi}^{2}\left[1-\frac{4{M_{q}^{\rm
vac}}^{2}N_{c}}{(4\pi)^{2}f_{\pi}^{2}}m_{\pi}^{2}F^{\prime}(m_{\pi}^{2})\right]{(\Delta^{2}+M_{q}^{2})^{2}\over{M_{q}^{\rm
vac}}^{4}}$ $\displaystyle-m_{\pi}^{2}f_{\pi}^{2}\left[1-{4{M_{q}^{\rm
vac}}^{2}N_{c}\over(4\pi)^{2}f_{\pi}^{2}}m_{\pi}^{2}F^{\prime}(m_{\pi}^{2})\right]{M_{q}\over{M_{q}^{\rm
vac}}}$
$\displaystyle-2N_{c}\int_{p}\left[\sqrt{\left(\sqrt{p^{2}+M_{q}^{2}}+\mu\right)^{2}+\Delta^{2}}+\sqrt{\left(\sqrt{p^{2}+M_{q}^{2}}-\mu\right)^{2}+\Delta^{2}}\right]$
$\displaystyle+4N_{c}\int_{p}\left[\sqrt{p^{2}+M_{q}^{2}+\Delta^{2}}+\frac{\mu^{2}\Delta^{2}}{2(p^{2}+M_{q}^{2}+\Delta^{2})^{3/2}}\right]\,.$
(132)

The function $F(p^{2})$ is given by

$\displaystyle F(p^{2})=2-2r\arctan\left({1\over r}\right)$ (133)
$\displaystyle p^{2}F^{\prime}(p^{2})={r^{2}+1\over r}\arctan\left({1\over
r}\right)-1$ (134)

where $r=\sqrt{4{M_{q}^{\rm vac}}^{2}/p^{2}-1}$. This parametrization suggests
that the parameters $m_{\pi}$ and $m_{\sigma}$ are restricted to
$m_{\pi},m_{\sigma}<2{M_{q}^{\rm vac}}$ to make $r$ real.

## References

  * Son and Stephanov [2001a] D. T. Son and M. A. Stephanov, QCD at finite isospin density, Phys. Rev. Lett. 86, 592 (2001a), arXiv:hep-ph/0005225 .
  * Son and Stephanov [2001b] D. T. Son and M. A. Stephanov, QCD at finite isospin density: From pion to quark - anti-quark condensation, Phys. Atom. Nucl. 64, 834 (2001b), arXiv:hep-ph/0011365 .
  * Splittorff  _et al._ [2001] K. Splittorff, D. T. Son, and M. A. Stephanov, QCD - like theories at finite baryon and isospin density, Phys. Rev. D 64, 016003 (2001), arXiv:hep-ph/0012274 .
  * Adhikari  _et al._ [2021a] P. Adhikari, J. O. Andersen, and M. A. Mojahed, Quark, pion and axial condensates in three-flavor finite isospin chiral perturbation theory, Eur. Phys. J. C 81, 449 (2021a), arXiv:2012.04339 [hep-ph] .
  * Adhikari  _et al._ [2021b] P. Adhikari, J. O. Andersen, and M. A. Mojahed, Condensates and pressure of two-flavor chiral perturbation theory at nonzero isospin and temperature, Eur. Phys. J. C 81, 173 (2021b), arXiv:2010.13655 [hep-ph] .
  * Adhikari and Andersen [2020a] P. Adhikari and J. O. Andersen, Quark and pion condensates at finite isospin density in chiral perturbation theory, Eur. Phys. J. C 80, 1028 (2020a), arXiv:2003.12567 [hep-ph] .
  * Adhikari and Andersen [2020b] P. Adhikari and J. O. Andersen, QCD at finite isospin density: chiral perturbation theory confronts lattice data, Phys. Lett. B 804, 135352 (2020b), arXiv:1909.01131 [hep-ph] .
  * Andersen  _et al._ [2019] J. O. Andersen, P. Adhikari, and P. Kneschke, Pion condensation and QCD phase diagram at finite isospin density, PoS Confinement2018, 197 (2019), arXiv:1810.00419 [hep-ph] .
  * Leggett and Zhang [2012] A. J. Leggett and S. Zhang, The bec-bcs crossover: Some history and some general observations, in  _The BCS-BEC Crossover and the Unitary Fermi Gas_, Lecture Notes in Physics, edited by W. Zwerger (2012) pp. 33–47.
  * Schrieffer [1999] J. Schrieffer,  _Theory Of Superconductivity_, Advanced Books Classics (Avalon Publishing, 1999).
  * [11] M. M. Parish, The BCS-BEC Crossover, Quantum Gas Experiments: Exploring Many-Body States. Edited by TORMA PAIVI ET AL. (World Scientific Publishing Co. Pte. Ltd, 2015).
  * Schäfer and Wilczek [1999] T. Schäfer and F. Wilczek, Continuity of quark and hadron matter, Phys. Rev. Lett. 82, 3956 (1999), arXiv:hep-ph/9811473 .
  * Hatsuda  _et al._ [2006] T. Hatsuda, M. Tachibana, N. Yamamoto, and G. Baym, New critical point induced by the axial anomaly in dense QCD, Phys. Rev. Lett. 97, 122001 (2006), arXiv:hep-ph/0605018 .
  * Baym  _et al._ [2018] G. Baym, T. Hatsuda, T. Kojo, P. D. Powell, Y. Song, and T. Takatsuka, From hadrons to quarks in neutron stars: a review, Rept. Prog. Phys. 81, 056902 (2018), arXiv:1707.04966 [astro-ph.HE] .
  * McLerran and Pisarski [2007] L. McLerran and R. D. Pisarski, Phases of cold, dense quarks at large N(c), Nucl. Phys. A 796, 83 (2007), arXiv:0706.2191 [hep-ph] .
  * Kojo [2021a] T. Kojo, QCD equations of state and speed of sound in neutron stars, AAPPS Bull. 31, 11 (2021a), arXiv:2011.10940 [nucl-th] .
  * Miller  _et al._ [2021] M. C. Miller _et al._ , The Radius of PSR J0740+6620 from NICER and XMM-Newton Data, Astrophys. J. Lett. 918, L28 (2021), arXiv:2105.06979 [astro-ph.HE] .
  * Riley  _et al._ [2021] T. E. Riley _et al._ , A NICER View of the Massive Pulsar PSR J0740+6620 Informed by Radio Timing and XMM-Newton Spectroscopy, Astrophys. J. Lett. 918, L27 (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Raaijmakers  _et al._ [2021] G. Raaijmakers, S. K. Greif, K. Hebeler, T. Hinderer, S. Nissanke, A. Schwenk, T. E. Riley, A. L. Watts, J. M. Lattimer, and W. C. G. Ho, Constraints on the dense matter equation of state and neutron star properties from NICER’s mass-radius estimate of PSR J0740+6620 and multimessenger observations, (2021), arXiv:2105.06981 [astro-ph.HE] .
  * Masuda  _et al._ [2013a] K. Masuda, T. Hatsuda, and T. Takatsuka, Hadron-Quark Crossover and Massive Hybrid Stars with Strangeness, Astrophys. J. 764, 12 (2013a), arXiv:1205.3621 [nucl-th] .
  * Masuda  _et al._ [2013b] K. Masuda, T. Hatsuda, and T. Takatsuka, Hadron–quark crossover and massive hybrid stars, PTEP 2013, 073D01 (2013b), arXiv:1212.6803 [nucl-th] .
  * Masuda  _et al._ [2016a] K. Masuda, T. Hatsuda, and T. Takatsuka, Hyperon Puzzle, Hadron-Quark Crossover and Massive Neutron Stars, Eur. Phys. J. A 52, 65 (2016a), arXiv:1508.04861 [nucl-th] .
  * Masuda  _et al._ [2016b] K. Masuda, T. Hatsuda, and T. Takatsuka, Hadron–quark crossover and hot neutron stars at birth, PTEP 2016, 021D01 (2016b), arXiv:1506.00984 [nucl-th] .
  * Bedaque and Steiner [2015] P. Bedaque and A. W. Steiner, Sound velocity bound and neutron stars, Phys. Rev. Lett. 114, 031103 (2015), arXiv:1408.5116 [nucl-th] .
  * Kojo  _et al._ [2015] T. Kojo, P. D. Powell, Y. Song, and G. Baym, Phenomenological QCD equation of state for massive neutron stars, Phys. Rev. D 91, 045003 (2015), arXiv:1412.1108 [hep-ph] .
  * Kojo [2016] T. Kojo, Phenomenological neutron star equations of state: 3-window modeling of QCD matter, Eur. Phys. J. A 52, 51 (2016), arXiv:1508.04408 [hep-ph] .
  * Fukushima and Kojo [2016] K. Fukushima and T. Kojo, The Quarkyonic Star, Astrophys. J. 817, 180 (2016), arXiv:1509.00356 [nucl-th] .
  * McLerran and Reddy [2019] L. McLerran and S. Reddy, Quarkyonic Matter and Neutron Stars, Phys. Rev. Lett. 122, 122701 (2019), arXiv:1811.12503 [nucl-th] .
  * Jeong  _et al._ [2020] K. S. Jeong, L. McLerran, and S. Sen, Dynamically generated momentum space shell structure of quarkyonic matter via an excluded volume model, Phys. Rev. C 101, 035201 (2020), arXiv:1908.04799 [nucl-th] .
  * Duarte  _et al._ [2020a] D. C. Duarte, S. Hernandez-Ortiz, and K. S. Jeong, Excluded-volume model for quarkyonic Matter: Three-flavor baryon-quark Mixture, Phys. Rev. C 102, 025203 (2020a), arXiv:2003.02362 [nucl-th] .
  * Duarte  _et al._ [2020b] D. C. Duarte, S. Hernandez-Ortiz, and K. S. Jeong, Excluded-volume model for quarkyonic matter. II. Three-flavor shell-like distribution of baryons in phase space, Phys. Rev. C 102, 065202 (2020b), arXiv:2007.08098 [nucl-th] .
  * Duarte  _et al._ [2023] D. C. Duarte, S. Hernandez-Ortiz, K. S. Jeong, and L. D. McLerran, Quarkyonic Mean Field Theory, (2023), arXiv:2302.04781 [nucl-th] .
  * Zhao and Lattimer [2020] T. Zhao and J. M. Lattimer, Quarkyonic Matter Equation of State in Beta-Equilibrium, Phys. Rev. D 102, 023021 (2020), arXiv:2004.08293 [astro-ph.HE] .
  * Cao and Liao [2020] G. Cao and J. Liao, A field theoretical model for quarkyonic matter, JHEP 10, 168, arXiv:2007.02028 [nucl-th] .
  * Margueron  _et al._ [2021] J. Margueron, H. Hansen, P. Proust, and G. Chanfray, Quarkyonic stars with isospin-flavor asymmetry, (2021), arXiv:2103.10209 [nucl-th] .
  * Kojo [2021b] T. Kojo, Stiffening of matter in quark-hadron continuity, Phys. Rev. D 104, 074005 (2021b), arXiv:2106.06687 [nucl-th] .
  * Hippert  _et al._ [2021] M. Hippert, E. S. Fraga, and J. Noronha, Insights on the peak in the speed of sound of ultradense matter, Phys. Rev. D 104, 034011 (2021), arXiv:2105.04535 [nucl-th] .
  * Freedman and McLerran [1978] B. Freedman and L. D. McLerran, Quark Star Phenomenology, Phys. Rev. D 17, 1109 (1978).
  * Freedman and McLerran [1977] B. A. Freedman and L. D. McLerran, Fermions and Gauge Vector Mesons at Finite Temperature and Density. 3. The Ground State Energy of a Relativistic Quark Gas, Phys. Rev. D 16, 1169 (1977).
  * Kurkela  _et al._ [2010] A. Kurkela, P. Romatschke, and A. Vuorinen, Cold Quark Matter, Phys. Rev. D 81, 105021 (2010), arXiv:0912.1856 [hep-ph] .
  * Annala  _et al._ [2020] E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, Evidence for quark-matter cores in massive neutron stars, Nature Phys. 16, 907 (2020), arXiv:1903.09121 [astro-ph.HE] .
  * Gorda  _et al._ [2021a] T. Gorda, A. Kurkela, R. Paatelainen, S. Säppi, and A. Vuorinen, Cold quark matter at NNNLO: soft contributions, (2021a), arXiv:2103.07427 [hep-ph] .
  * Gorda  _et al._ [2021b] T. Gorda, A. Kurkela, R. Paatelainen, S. Säppi, and A. Vuorinen, Soft interactions in cold quark matter, (2021b), arXiv:2103.05658 [hep-ph] .
  * Suenaga and Kojo [2019] D. Suenaga and T. Kojo, Gluon propagator in two-color dense QCD: Massive Yang-Mills approach at one-loop, Phys. Rev. D 100, 076017 (2019), arXiv:1905.08751 [hep-ph] .
  * Kojo and Suenaga [2021] T. Kojo and D. Suenaga, Thermal quarks and gluon propagators in two-color dense QCD, Phys. Rev. D 103, 094008 (2021), arXiv:2102.07231 [hep-ph] .
  * Fujimoto and Fukushima [2020] Y. Fujimoto and K. Fukushima, Equation of state of cold and dense QCD matter in resummed perturbation theory, (2020), arXiv:2011.10891 [hep-ph] .
  * Ma and Rho [2020] Y.-L. Ma and M. Rho, Towards the hadron–quark continuity via a topology change in compact stars, Prog. Part. Nucl. Phys. 113, 103791 (2020), arXiv:1909.05889 [nucl-th] .
  * Fujimoto  _et al._ [2022] Y. Fujimoto, K. Fukushima, L. D. McLerran, and M. Praszalowicz, Trace Anomaly as Signature of Conformality in Neutron Stars, Phys. Rev. Lett. 129, 252702 (2022), arXiv:2207.06753 [nucl-th] .
  * Marczenko  _et al._ [2022] M. Marczenko, L. McLerran, K. Redlich, and C. Sasaki, Reaching percolation and conformal limits in neutron stars, (2022), arXiv:2207.13059 [nucl-th] .
  * Blaschke  _et al._ [2023] D. Blaschke, U. Shukla, O. Ivanytskyi, and S. Liebing, Effect of color superconductivity on the mass of hybrid neutron stars in an effective model with perturbative QCD asymptotics, Phys. Rev. D 107, 063034 (2023), arXiv:2212.14856 [nucl-th] .
  * Ivanytskyi and Blaschke [2022] O. Ivanytskyi and D. B. Blaschke, Recovering the Conformal Limit of Color Superconducting Quark Matter within a Confining Density Functional Approach, Particles 5, 514 (2022), arXiv:2209.02050 [nucl-th] .
  * Adhikari  _et al._ [2017] P. Adhikari, J. O. Andersen, and P. Kneschke, On-shell parameter fixing in the quark-meson model, Phys. Rev. D 95, 036017 (2017), arXiv:1612.03668 [hep-ph] .
  * Adhikari  _et al._ [2018] P. Adhikari, J. O. Andersen, and P. Kneschke, Pion condensation and phase diagram in the Polyakov-loop quark-meson model, Phys. Rev. D 98, 074016 (2018), arXiv:1805.08599 [hep-ph] .
  * Ayala  _et al._ [2023] A. Ayala, A. Bandyopadhyay, R. L. S. Farias, L. A. Hernández, and J. L. Hernández, QCD equation of state at finite isospin density from the linear sigma model with quarks: The cold case, (2023), arXiv:2301.13633 [hep-ph] .
  * Baym  _et al._ [2019] G. Baym, S. Furusawa, T. Hatsuda, T. Kojo, and H. Togashi, New Neutron Star Equation of State with Quark-Hadron Crossover, Astrophys. J. 885, 42 (2019), arXiv:1903.08963 [astro-ph.HE] .
  * Kojo  _et al._ [2022] T. Kojo, G. Baym, and T. Hatsuda, Implications of NICER for Neutron Star Matter: The QHC21 Equation of State, Astrophys. J. 934, 46 (2022), arXiv:2111.11919 [astro-ph.HE] .
  * Novikov  _et al._ [1978] V. A. Novikov, L. B. Okun, M. A. Shifman, A. I. Vainshtein, M. B. Voloshin, and V. I. Zakharov, Charmonium and Gluons: Basic Experimental Facts and Theoretical Introduction, Phys. Rept. 41, 1 (1978).
  * Shifman  _et al._ [1979a] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, QCD and Resonance Physics. Theoretical Foundations, Nucl. Phys. B 147, 385 (1979a).
  * Shifman  _et al._ [1979b] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, QCD and Resonance Physics: Applications, Nucl. Phys. B 147, 448 (1979b).
  * Brandt  _et al._ [2022] B. B. Brandt, F. Cuteri, and G. Endrodi, Equation of state and speed of sound of isospin-asymmetric QCD on the lattice, (2022), arXiv:2212.14016 [hep-lat] .
  * Ananthanarayan  _et al._ [2017] B. Ananthanarayan, I. Caprini, and D. Das, Electromagnetic charge radius of the pion at high precision, Phys. Rev. Lett. 119, 132002 (2017), arXiv:1706.04020 [hep-ph] .
  * Tanabashi et al. [2018] M. Tanabashi et al. (Particle Data Group), Review of particle physics, Phys. Rev. D 98, 030001 (2018).
  * Koponen  _et al._ [2016] J. Koponen, F. Bursa, C. T. H. Davies, R. J. Dowdall, and G. P. Lepage, Size of the pion from full lattice QCD with physical u , d , s and c quarks, Phys. Rev. D 93, 054503 (2016), arXiv:1511.07382 [hep-lat] .
  * Wang  _et al._ [2021] G. Wang, J. Liang, T. Draper, K.-F. Liu, and Y.-B. Yang (chiQCD), Lattice Calculation of Pion Form Factor with Overlap Fermions, Phys. Rev. D 104, 074502 (2021), arXiv:2006.05431 [hep-ph] .
  * Kojo and Suenaga [2022a] T. Kojo and D. Suenaga, Peaks of sound velocity in two color dense QCD: Quark saturation effects and semishort range correlations, Phys. Rev. D 105, 076001 (2022a), arXiv:2110.02100 [hep-ph] .
  * Kojo and Suenaga [2022b] T. Kojo and D. Suenaga, Meson resonance gas in a relativistic quark model: scalar vs vector confinement and semishort range correlations, (2022b), arXiv:2208.13312 [hep-ph] .
  * Drischler  _et al._ [2022] C. Drischler, S. Han, and S. Reddy, Large and massive neutron stars: Implications for the sound speed within QCD of dense matter, Phys. Rev. C 105, 035808 (2022), arXiv:2110.14896 [nucl-th] .
  * Fukushima  _et al._ [2020] K. Fukushima, T. Kojo, and W. Weise, Hard-core deconfinement and soft-surface delocalization from nuclear to quark matter, Phys. Rev. D 102, 096017 (2020), arXiv:2008.08436 [hep-ph] .
  * Gao  _et al._ [2022] B. Gao, T. Minamikawa, T. Kojo, and M. Harada, Impacts of the U(1)A anomaly on nuclear and neutron star equation of state based on a parity doublet model, Phys. Rev. C 106, 065205 (2022), arXiv:2207.05970 [nucl-th] .
  * Minamikawa  _et al._ [2023] T. Minamikawa, B. Gao, T. Kojo, and M. Harada, Chiral restoration of nucleons in neutron star matter: studies based on a parity doublet model, (2023), arXiv:2302.00825 [nucl-th] .
  * Graf  _et al._ [2016] T. Graf, J. Schaffner-Bielich, and E. S. Fraga, Perturbative thermodynamics at nonzero isospin density for cold QCD, Phys. Rev. D 93, 085030 (2016), arXiv:1511.09457 [hep-ph] .
  * Deur  _et al._ [2015] A. Deur, S. J. Brodsky, and G. F. de Teramond, Connecting the Hadron Mass Scale to the Fundamental Mass Scale of Quantum Chromodynamics, Phys. Lett. B 750, 528 (2015), arXiv:1409.5488 [hep-ph] .
  * Deur  _et al._ [2016] A. Deur, S. J. Brodsky, and G. F. de Teramond, On the Interface between Perturbative and Nonperturbative QCD, Phys. Lett. B 757, 275 (2016), arXiv:1601.06568 [hep-ph] .
  * Son [1999] D. T. Son, Superconductivity by long range color magnetic interaction in high density quark matter, Phys. Rev. D 59, 094019 (1999), arXiv:hep-ph/9812287 .
  * Polyakov [2003] M. V. Polyakov, Generalized parton distributions and strong forces inside nucleons and nuclei, Phys. Lett. B 555, 57 (2003), arXiv:hep-ph/0210165 .
  * Fujita  _et al._ [2022] M. Fujita, Y. Hatta, S. Sugimoto, and T. Ueda, Nucleon D-term in holographic quantum chromodynamics, PTEP 2022, 093B06 (2022), arXiv:2206.06578 [hep-th] .
  * Sakai and Sasaki [2023] K. Sakai and S. Sasaki, Glueball spectroscopy in lattice QCD using gradient flow, Phys. Rev. D 107, 034510 (2023), arXiv:2211.15176 [hep-lat] .
  * Blacker  _et al._ [2023] S. Blacker, A. Bauswein, and S. Typel, Exploring thermal effects of the hadron-quark matter transition in neutron star mergers, (2023), arXiv:2304.01971 [astro-ph.HE] .
  * Peterson  _et al._ [2023] J. Peterson, P. Costa, R. Kumar, V. Dexheimer, R. Negreiros, and C. Providencia, Temperature and Strong Magnetic Field Effects in Dense Matter, (2023), arXiv:2304.02454 [nucl-th] .
  * Kojo  _et al._ [2021] T. Kojo, D. Hou, J. Okafor, and H. Togashi, Phenomenological QCD equations of state for neutron star dynamics: Nuclear-2SC continuity and evolving effective couplings, Phys. Rev. D 104, 063036 (2021), arXiv:2012.01650 [astro-ph.HE] .
  * Huang  _et al._ [2022] Y.-J. Huang, L. Baiotti, T. Kojo, K. Takami, H. Sotani, H. Togashi, T. Hatsuda, S. Nagataki, and Y.-Z. Fan, Merger and Postmerger of Binary Neutron Stars with a Quark-Hadron Crossover Equation of State, Phys. Rev. Lett. 129, 181101 (2022), arXiv:2203.04528 [astro-ph.HE] .
  * Fujimoto  _et al._ [2023] Y. Fujimoto, K. Fukushima, K. Hotokezaka, and K. Kyutoku, Gravitational Wave Signal for Quark Matter with Realistic Phase Transition, Phys. Rev. Lett. 130, 091404 (2023), arXiv:2205.03882 [astro-ph.HE] .
  * Kedia  _et al._ [2022] A. Kedia, H. I. Kim, I.-S. Suh, and G. J. Mathews, Binary neutron star mergers as a probe of quark-hadron crossover equations of state, Phys. Rev. D 106, 103027 (2022), arXiv:2203.05461 [gr-qc] .
  * Annala  _et al._ [2023] E. Annala, T. Gorda, J. Hirvonen, O. Komoltsev, A. Kurkela, J. Nättilä, and A. Vuorinen, Strongly interacting matter exhibits deconfined behavior in massive neutron stars, (2023), arXiv:2303.11356 [astro-ph.HE] .
  * Han  _et al._ [2022] M.-Z. Han, Y.-J. Huang, S.-P. Tang, and Y.-Z. Fan, Plausible presence of new state in neutron stars with masses above $0.98M_{\rm TOV}$ 10.1016/j.scib.2023.04.007 (2022), arXiv:2207.13613 [astro-ph.HE] .
  * Geesaman  _et al._ [1995] D. F. Geesaman, K. Saito, and A. W. Thomas, The nuclear EMC effect, Ann. Rev. Nucl. Part. Sci. 45, 337 (1995).
  * Saito  _et al._ [2007] K. Saito, K. Tsushima, and A. W. Thomas, Nucleon and hadron structure changes in the nuclear medium and impact on observables, Prog. Part. Nucl. Phys. 58, 1 (2007), arXiv:hep-ph/0506314 .