Sound velocity peak and conformality in isospin QCD

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.


I. INTRODUCTION
The quantum chromodynamics (QCD) with a large isospin chemical potential (µ 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., Ref. [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,16], which may be also realized in QCD at finite baryon chemical potential (µ B ).
One of fundamental topics in dense QCD is the equation of state (EOS) (Ref.[17] 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 (≃ 12.4 km) for 2.1 and 1.4 solar mass NSs [18][19][20], and the nuclear physics constraints at nuclear saturation density n 0 (≃ 0.16 fm −3 ), suggest that the EOS stiffens rapidly (i.e., the pressure P grows rapidly as a function of energy density ε) 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 = (∂P/∂ε) 1/2 ; the c 2 s is ≪ 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 [21][22][23][24], discussed in more general ground based on nuclear physics and NS observations by Ref. [25], and further elucidated in Refs.[26][27][28] utilizing the quark degrees of freedom.Recently more detailed descriptions have been attempted, see Refs.[29][30][31][32][33][34][35][36][37][38], 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 2 s .
Another interesting question is how the c 2 s approaches the conformal limit, 1/3.Perturbative QCD (pQCD) [39][40][41][42][43][44], which is supposed to be valid at n B ≳ 40n 0 , predicts that the c 2 s approaches 1/3 from below.The domain between n B ≃ 10n 0 and ≃ 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 [45][46][47]; if such interaction effects are properly absorbed into effective parameters of quasi-quarks, it is possible that 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 [42,[48][49][50][51][52][53].
We address the non-monotonic behavior of c 2 s in isospin QCD within a renormalizable quark-meson model.The properties of the model in isospin QCD have been analyzed in detail by Refs.[54][55][56]. 1 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 θ(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 2 s toward the high density limit [57,58].In the quark-meson model such growing behaviors disappear; the c 2 s relaxes to the conformal value 1/3, as it should.
While our model predicts c 2 s → 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 α s .In the weak coupling limit and at high density, the only relevant scale is µ I and the c 2 s = 1/3 follows from ∂(P/µ 4 I )/∂µ I = 0.The first important corrections to the conformal limit come from the Λ QCD ≃ 200-300 MeV in the running coupling constant.If we take into account Λ QCD only in this way, the c 2 s is reduced from 1/3.Meanwhile, at the energy scale around ∼ 1 GeV, it has been long known that power corrections of Λ QCD , which can not be expressed as perturbative series in α s , play important roles to capture the qualitative features in QCD [59][60][61].Parametrizing pressure with power corrections as (χ I : isospin susceptibility) where a 2 ∼ Λ 2 QCD , the squared sound velocity can be expressed as For a positive a 2 , the c 2 s 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.
We briefly address the trace anomaly in dense matter which measures the breaking of the scale invariance [48][49][50].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.In this respect the sign of the trace anomaly is very useful to characterize the non-perturbative effects in dense quark matter.
For quantitative aspects, we confront our model calculations with the lattice results from two groups [62,63].Ref. [62] have more focus on the BEC regime while Ref. [63] covers more global nature up to the pQCD regime.Both groups agree in the BEC regime, and our model results are consistent with the lattice results.At high density, our model captures the overall trend of Ref. [63], especially the sound velocity peak and negative trace anomaly.
Since we are not sure about the convergence of loop expansion, as supplement studies we perform several parametric studies of EOS to examine several qualitative effects which we believe to be important.They are used to delineate the results of the quark-meson model in Sec.IV.
In this paper we use nuclear saturation density in QCD, n 0 = 0.16 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 an indicator of the nonperturbative effects.Section V is devoted to a summary.

II. MODEL
The Lagrangian of the two-flavor quark-meson model is where ψ is a quark field with up-and down-quark components, ψ = (u, d) T .The ⃗ ϕ = (σ, ⃗ π) are meson fields which correspond to the isospin 1 and 3 representations.The τ 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 µ I = µ u = −µ d as a Lagrange multiplier.To correctly identify the corresponding Lagrangian, we should begin with the hamiltonian formalism.The thermodynamic potential is The isospin density in terms of field variables can be identified by the Noether theorem.Meson and quark fields transform under isospin transformations as and corresponding conserved current can be written as where the ϵ abc is the complete anti-symmetric tensor with ϵ 123 = 1.The isospin density is now Writing fields collectively as Φ = ( ⃗ ϕ, ψ), the partition function for Φ is Here the Lagrangian at finite density is where we attached the subscript B to emphasize that the parameters and fields are unrenormalized.
Below we construct a one-loop effective potential within the leading order of the 1/N c expansion.In this approximation, meson loop effects on quarks are neglected, while quark loop effects play crucial roles in renormalizing meson parameters as well as the amplitude of meson condensates.This quark substructure affects the density evolution of meson condensates and hence the EOS.
First we rewrite the Lagrangian using the renormalized parameters and fields.We begin with the O(4) symmetric scheme and later relate those renormalized parameters to those in the on-shell scheme.The bare parameters are written with O(4) symmetric renormalized fields and couplings as We also define δZ i = Z i − 1 for i = ϕ, ψ, g, • • • and so on.
The Zi represents the radiative corrections without those for the external lines.In our model, the loop corrections to the quark self-energies and quark-meson vertices appear only through meson-loops and hence Meanwhile, the meson self-energies and tadpole contain quark loops of O(N c ) which are combined with g 2 ∼ 1/N c vertices to yield and hence one must keep the corrections.It is useful to note that the relation in the large N c limit.The first relation means that the dynamically generated quark mass and gap are renormalization group (RG) invariant.The second relation tells that we need to study the meson propagators to describe the running of g 2 .Now the Lagrangian is decomposed into the renormalized part and counter terms as Here L is the renormalized Lagrangian where all subscripts B are omitted from L B and couplings are replaced with the renormalized couplings.The counter terms necessary in the large N c limit is The counter Lagrangian is used when we calculate loop corrections.
We construct the effective potential Γ(ϕ 0 ) with the MS normalization of fields.The effective potentials defined at different renormalization schemes are related as while the kinetic terms are always normalized to 1. Actually, it is more convenient to work with a Γ(gϕ 0 ) in which gϕ 0 is the RG invariant in the large N c limit.We specify M q = gσ 0 and ∆ = g(π 1 ) 0 as variables for the effective potential and then Γ R (M q , ∆) = Γ R ′ (M q , ∆), i.e., we need only to take into account changes in (g, λ, , • • • ) when we change the renormalization conditions.
A. Parameter fixing with vacuum quantities

Renormalization of effective potenital
Now we fix the counter terms by renormalizing physical parameters in vacuum.The simplest scheme to obtain the renormalized effective potential is the MS scheme.The effective potential takes the form The V q is the one-loop contributions from the quark energy where we write M q = gσ 0 .We will treat these divergences by dimensional regularization where Λ is the renormalizing scale introduced by the MS scheme and γ E = 0.577... is the Euler-Mascheroni constant.The quark energy now reads There is no 1/ϵ pole in the linear and quadratic terms, while the 1/ϵ in the quartic term is cancelled by δ λ, The effective potential in vacuum now reads We demand the effective potential to be RG invariant, i.e., the effective potential does not change by replacement of Λ → Λ ′ .This must be valid for a given M q , so each coefficient of M q must be invariant.The invariance of M q and M 2 q terms requires (m/g, h/g) do not run in the large N c limit, as consistent with Eq. (11).Meanwhile the M 4 q terms contain the ln Λ 2 factor so that or The running of g 2 is obtained from the analyses of field normalizations of ϕ, thanks to the relation g 2 (Λ) = g 2 B /Z ϕ (Λ).

Renormalization of meson propagators
Now we consider the renormalization conditions for mesons to fix the Z ϕ .We write M 0 as a solution to minimize the effective potential and use it to compute the meson self-energies.We demand that the σ and π have the pole at p 2 = m 2 σ and p 2 = m 2 π , respectively, The self-energies include the quark one-loop contributions Σ q σ,π and counter terms The quark loop Σ q σ,π is UV divergent.The counter term δ λ automatically cancels the UV divergences coupled to M 0 .The δZ ϕ is arranged to cancel m σ -and m π -dependent UV divergences in Σ q σ,π , with which where G σ and G π are functions of p 2 and M 2 0 , It is useful to note G σ (4M 2 0 ) = 0 and G π (0) = 0. Later we also make use of F (0) = 0 and F (4M 2 0 ) = −2.We note that the parameter Λ manifestly appears in the self-energies but the pole positions should be RG invariant.This demands Eqs.( 14) and ( 24) lead to The effective potential and the pole locations with running parameters are RG invariant, so below we choose Λ = M 0 to get rid of the ln Λ terms.
Finally we also mention how the MS and on-shell renormalization schemes are related.Here we discuss only Z OS σ as we will fix f π by ⟨σ OS ⟩ = f π .We note where the residue of ⟨ϕϕ⟩ OS is normalized to 1. Thus In the parameter range of our interest, we find Z ϕ < Z OS σ .For instance, for theories with m σ = 2M 0 , the inequality is verified by noting F (4M 2 0 ) = −2.

Parameter fixing
We fix the values of parameters in our model.To evaluate the effective potential, we need to fix four parameters (m, g, λ, h) at Λ = M 0 .Our input is (f π , M 0 , m σ , m π ).
First we fix the value of g.We note that M 0 is RG invariant (in the large N c limit), We can fix g OS = M 0 /f π while g can be fixed by the relation For typical parameter set M 0 ∼ 300 MeV and f π ∼ 90 MeV, g OS ∼ 3.3 which is large.In the MS scheme g 2 is smaller and the expansion of g 2 is slightly better in systematics.
Having g fixed, we can determine m 2 and λ from the pole conditions for m σ and m π , To get analytic insights, it is again useful to consider the case m σ = 2M 0 and m π = 0. Then In this limit, it is clear that m 2 , which drives the σ condensation at tree level, becomes more negative than the tree level counterpart by radiative corrections.This limit also suggests that λ is typically large, of O(10-100).
Finally we fix h.Using the parameters defined at Λ = M 0 , the effective potential takes the form The gap equation at M q = M 0 fixes the value of h, Using the condition for the pion pole, one can rewrite it as where h OS = m 2 π f π , the standard expression in the chiral EFT.

B. At finite isospin density
For a large isospin chemical potential, either π 1 or π 2 can condense while π 3 fields are unaffected.Without loss of generality we assume the π 1 to condense.The quark part in the unperturbed Lagrangian acquires an extra term with which the quark propagator becomes the BCS type propagator.The poles exist at where (see the derivation in appendix A) The u and d quark excitations cost at least the energies of the BCS gap ∼ ∆.Meanwhile d and ū quarks need large energies of ∼ M + µ I to get excited.The effective potential in the MS scheme is We note that the µ I dependent term contains a UV divergent counter term which is necessary to cancel a µ I dependent UV divergence from V q .
The single particle energies depend on M q and ∆ in medium, To get analytic insights, we split where the upper script specifies the power of µ I , whose computations can be carried out with the dimensional regularization, We have extracted up to µ 2 I terms as they contain the UV divergences, while V R q is UV finite and contains terms which scale as µ 4 I and vanish at µ I = 0.At large µ I , V R q dominates over the other terms as far as ∆ and M q do not grow as ∼ µ I .We numerically evaluate V R q and found that The effective potential in the MS scheme now reads The effective potential rewritten with hadronic parameters is shown in Appendix.B and we use it to evaluate M q and ∆, as was done in Refs.[54,55].The expectation value M q * and ∆ * are determined by the gap equations, 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.Unless otherwise stated, we fix the model parameters to satisfy the following vacuum parameters 2 : which correspond to the following on-shell coupling constants, g OS ≃ 3.33 and λ OS ≃ 126. 3  The large couplings in the present one-loop analyses are worrisome.Meanwhile it has been known that constituent quark type models with couplings of O(1) work remarkably well without rigorous justifications.In this 2 Here we have used the sigma mass as the renormalization condition but in reality the sigma or f 0 (500) state has a broad width.This width has been studied and confirmed in the linear sigma model, which is very similar to this model, by Ref. [64] considering the σ → ππ scattering process.In our study at large Nc, quark loops enter only condensed mesons and counter terms for mesons, but do not affect mesonic fluctuations or meson excitations, and hence the impacts of meson width are not addressed. 3In the MS scheme, the couplings are smaller.The details depend on the choice of mσ which is more uncertain than the other input parameters.For mσ = 2M 0 , there is simple relation from which g ≃ 2.03.This reduces the value of λ by a factor ∼ (g/g OS ) 2 ≃ 0.37.The value of λ becomes even smaller for mπ → mσ.
work we simply hope that the similar situation holds in our studies.We also note that, in the case of the nucleon-meson model, whose structure is very similar to the quark-meson model, the tree and one-loop results are qualitatively different, but the difference between oneloop results and the functional renormalization group results are quantitative one, the order of ∼ 30% [65].Thus we expect our one-loop results to be useful to gain some qualitative insights into the overall trend of isospin QCD.With this qualification in mind, we proceed to the examination of the EOS.For comparison to the lattice data in Ref. [63], later we also examine the m π = 170 MeV case with (m σ , f π , M 0 ) kept the same as the m π = 140 MeV case.First we examine how condensates evolve as functions of µ I .Shown in Fig. 1 are the constituent quark mass and the gap associated with the pion condensate.For M q , there are no significant differences between the tree-level (dashed blue line) and one-loop (solid blue line) results.Meanwhile, the pion condensate ∆ at tree level increases linearly with µ I , whereas, at one-loop, it converges to a finite value, ∆ ≃ 300 MeV.This drastic change of behavior indicates that the one-loop correction has more physical contents than mere perturbative corrections.

A. Evolution of microscopic quantities
At tree-level, the Lagrangian makes no reference to quarks so that the mesons are treated as elementary particles (Fig. 2 (Top)).By adding quark loops, however, they no longer can be regarded as purely elementary particles.If we regard mesons as fundamental, quark loops are regarded as corrections to the meson dynamics.But if we regard quark descriptions as more fundamental, mesons are intermediate states appearing in the quarkantiquark scattering processes (Fig. 2 (Bottom)).
In this study, we keep only the leading N c contributions and hence the quark substructure effects on meson fluctuations are not reflected in EOS (as meson loops are suppressed by 1/N c ).However, the quark substructure effects do affect condensed mesons by tempering the amplitudes.The quark loops change the structure of the present theory and in this sense it may not be appropriate to call quark loops as corrections; rather they should be regarded as leading order contributions.Now we examine how quark loops qualitatively change the behavior of ∆.To address this question we work with the MS expression for the moment as it takes a concise form.We consider a large µ I and assume M q ≪ ∆.At tree level, the effective potential behaves as and the solution of the gap equation is found by balancing µ 2 ∆ 2 and ∆ 4 terms.Hence ∆ tree ∼ µ I inevitably follows.
Note also that λ > 0, like the hard core repulsion, plays an essential role to stop the growth of ∆.
Including quark loops, however, the coefficient of µ 2 I term acquires the ln(∆/M 0 ) 2 term which, before ∆ 4 becomes dominant, can stop the growth of ∆.At large µ I and assuming ∆ ≪ µ I , where V R q ∼ µ 4 I weakly depends upon ∆.Then the gap equation is determined by the coefficient of The solution is µ I -independent, 4 4 The g 2 -dependence in this expression may be confusing and here we give supplemental comments.At large g ≫ 1, ∆ 2 → M 2 0 e −1 which looks smaller than M 2 0 .This reduction is fictitious; if we hold (m, λ, h) fixed and increase g, M 0 ∼ gfπ increases and ∆ also increases.On the other hand, at small g ≪ 1, apparently ∆ becomes larger, but our assumption of µ I ≫ ∆ and hence our estimate is violated.In this situation, the size of ∆ is primarily determined by the tree level relation as the hadronic and quark sectors decouple for g → 0.
For our parameter set, (4π) 2 /4g 2 N c ∼ 1 and the exponent is small; we find ∆ ≃ M 0 as shown in Fig. 1.
Substituting the solution into Eq.(54), the 1/g 2 and the logarithmic terms cancel, leaving the −µ 2 I ∆ 2 * term.As a result the pressure P (µ where P R q = −V R q (M * q , ∆ * ).The µ 2 I dependence can be interpreted as the Fermi surface effects with the phase space ∼ 4πp 2 F with p F being the quark Fermi momentum p F ∼ µ I .
In this expression for large µ I , hadronic parameters disappear.The hadronic parameters m 2 and λ are neglected because they appear as m 2 ∆ 2 and λ∆ 4 terms much smaller than µ 2 I ∆ 2 and µ 4 I , while g 2 is absorbed into the expression of ∆ * .The resulting expression can be most naturally understood in terms of quarks with non-perturbative effects near the Fermi surface whose strength depends upon the hadron physics.

B. Equations of state
Starting with the thermodyamic pressure P (µ q , ∆ * ; µ I ), the isospin and energy densities are given by We study the sound velocity where χ I is the isospin susceptibility.
In the following we compare our results with the lattice data in Refs.[62] and [63].The setup of the former is N f = 2 + 1 flavors of rooted staggered quarks with the quark masses at the physical point.The pion decay constant is f π ≃ 92-96 MeV for the lattice spacing explored (the definition of f π differs by a factor √ 2 from ours and we have corrected it).It should be noted that their results at T = 0 are obtained by correcting the data at small but finite T using the ChEFT.Beyond µ I ≳ m π or n I ≳ 0.5n 0 the lattice data is not available in Ref. [62].Meanwhile, the lattice data in Ref. [63] using m π ≃ 170 MeV and a different formalism is more suitable to explore high density region up to µ I ∼ 7.5m π ≃ 1.3 GeV (our definition of µ I is a half of that in Ref. [63], taken into account in our figures).
Figure 3 shows the isospin density n I as a function of the isospin chemical potential µ I , for the global (upper panel) and low density (lower panel) behaviors.We use (m σ , f π ) = (600, 90) MeV and consider m π = 140 and 170 MeV for comparison with the lattice data of Refs.[62] [62] (available up to µI ≃ 0.9mπ with mπ ≃ 140 MeV) and bands are from Abbott et al. [63] with mπ ≃ 170 MeV.(The definition of µI in Ref. [63] differs from ours by a factor 2 and this is taken into account in our figures.) and [63].The upper panel in Fig. 3 is specialized for the examination of global features to high density.As expected from the qualitative difference in the behavior of condensates, the tree and one-loop results become very different toward high density.In purely hadronic descriptions, we found P ∼ λ∆ 4 tree ∼ λµ 4 I whose asymptotic behavior is controlled by the hadronic coupling λ, the strength of "hard core repulsion" between mesons.This scaling behaviors are changed by quark loops, with which the scaling P ∼ cµ 4  I is controlled by the phase space factor for quarks, rather than parameters for hadronic interactions.
In the lower panel of Fig. 3, we make more detailed comparison at low density.We note that the ChEFT results of Ref. [7] including chiral loop corrections agree well with the lattice results of Ref. [62] (see Fig. 2 in Ref. [7] for various comparisons), while our model results slightly underestimate n I for a larger µ I .We note that the chiral loops in the ChEFT and quark-loops cover different types of quantum fluctuations.In Fig. 4, we also show the relation between P and µ I for the low density and global behaviors.
Shown in Fig. 6 is the pressure as a function of energy density.For the density range to ∼ 10n 0 , the pressure with quark loops is larger than that that in the tree level by 10-20%.This means that the quark substructure effects enhance the pressure from the purely hadronic one.On the other hand, toward high density the difference in P vs ε becomes much smaller than in n I vs µ I .Such degeneracy is reached when P enters the µ 4 I scaling regime; for whatever coefficients of the µ 4 I term, P ≃ ε/3 is achieved when µ 4 I terms dominate.We also note that the pressure is reduced for larger m σ and f π , as shown in Fig. 5.In other words, with stronger chiral symmetry breaking in vacuum (which increases both m σ and f π ), the high density EOS after the chiral restoration becomes softer.This point is examined in Sec.IV A. To further examine the variation of stiffness, we now turn to the behavior of the sound velocity.Figure 7 shows the c 2 s as a function of isospin chemical potential and also as a function of isospin density.The c 2 s 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 results in Refs.[62] and [63].However, the quantitative agreement beyond the BEC regime depends on the lattice results.The location of the c 2 s peak is near µ I ≃ 1.2m π for n I ≃ 5n 0 in our calculations for the reasonable range of our parameter set for m σ and f π .The lattice results in Ref. [62] indicate the peak at µ I ≃ 0.8m π or n I ∼ 0.5n 0 , lower than our model results.Meanwhile our results agree better with the results of Ref. [63], although our sound velocity peak is located at density slightly lower than found in the lattice simulations, 6-7n 0 .We are not sure about the origin of the discrepancy between results of Refs.[62] and [63] as they seem to contain different systematic errors.But after performing several parametric studies as given in Sec.IV, we could not find any qualitative mechanisms to reconcile ∆ ≃ Λ QCD with the quick reduction of c 2 s after making the peak as seen in Ref. [62].Here we assume ∆ ≃ Λ QCD based on the lattice result for the melting temperature of pion condensates, T lat c ≃ 170 MeV, which seems more or less constant to ∼ n 0 or even higher densities, see discussions around Eq. (78).For this reason, in the beyond-BEC regime, the results of Ref. [63] seem more natural to us than those of Ref. [62] whose simulations are more optimized for the low density region.
We note that the c 2 s 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 behaviors can be acheived by P (µ I ) ∼ c 2 µ 2 I Λ 2 + c 4 µ 4 I with Λ being some scale.Purely hadronic models may achieve this condition, but this by itself does not mean that the EOS is described correctly, as we have mentioned in discussion of P vs ε and n I vs µ I .In our standpoint, the tree level results, which crucially depend on the scaling ∆ tree ∼ µ I , becomes potentially misleading at high density.
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.

C. 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 refer to the pion charge radius.It can be extracted from the vector form factor.The experimental determination based on the πe scattering and the e + e − → π + π − process [66] yield the estimate ⟨r 2 ⟩ V = 0.434(5) fm 2 [67], or which has been well reproduced by lattice calculations [68,69].The typical isospin density where pions overlap is estimated through5 ) Figure 3 shows that the isospin chemical potential at We note that this overlap density n overlap I ∼ 10n 0 is substantially larger than the density 2-3n 0 where the tree and one-loop results begin to differ substantially, and the density ∼ 5n 0 where c 2 s develops a peak.This would indicate that the quark substructure of hadrons become important before hadrons overlap.In this respect, there should be a more suitable measure to characterize the location of sound velocity peak.One of possible explanations is the quark saturation [37,70,71].As density increases, quark states at low momentum are inevitably occupied and then a newly added quark must fill a state on top of the already occupied states.Quark states at large momenta are the source of large pressure.
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-, ū-, d-quarks are Roughly speaking, u-and d-quarks occupy states up to ≃ µ I while ū-and d-quarks are almost fully occupied as in the Dirac sea without pion condensates.Shown in Fig. 9 is the occupation probability f (p) at various densities as a function of quark momenta p.The densities we plotted are from 0 to 10.0n 0 in 0.2n 0 increments for gray curves and 1.0n 0 increments for red curves.The blue dots, where f (p) takes the half value of f (p = 0), are the measure of typical momentum at the Fermi surface.
For later convenience we define the quark distribution in a single pion as φ vac π (p) ≡ lim n I →0 f (p)/n I .It turns out that φ vac π (p) is approximated well by a simple monopole Ansatz or Breit-Wigner form with a ≃ 87.6 MeV.This suggests that, in our quark- meson model, the pion wavefunction in the coordinate space has the exponentially-decaying form.
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) ≃ n I φ vac π (p) (Fig. 10, Left); this corresponds to the regime where pions do not interact and quarks inside of pions are largely unaffected.In this regime, ε/n I is close to a constant, and therefore the pressure, P = n 2 I ∂(ε/n I )/∂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. 10, Right).This is driven by both interactions and the Pauli blocking effects.Here, ε/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. 9 suggests that, from 0 to ∼ 2n 0 , the magnitude of f (p) at p = 0 grows rapidly from 0 to ≃ 0.6, but at higher density the distribution f (p) develops toward the horizontal direction.If we treated pions as if elementary and non-ineteracting particles, the f (p = 0) would violate the Pauli principle around ≃ 2n 0 .The c 2 s peak is located around ≃ 5n 0 where f (p = 0) ≃ 0.9.Beyond this density the horizontal evolution dominates over the vertical evolution and c 2 s relaxes toward 1/3 as in a relativistic quark gas.We note that, the quark substructure effects are already significant at 1-2n 0 and develops a peak in c 2 s at ∼ 5n 0 , at density substantially smaller than the naive estimate of the pion overlap, n overlap I ∼ 10n 0 .This suggests that 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. [72] which discriminates the mode-by-mode percolation in momentum space from the conventional FIG. 10.Schematic figures for the evolution of the occupation probability.Left: The "vertical" evolution; Right: The "horizontal" evolution.
Finally we comment on some difference between nucleonic matter and pionic matter in isospin QCD.In nuclear matter the evolution of c 2 s is much slower than in isospin QCD, c 2 s ≲ 0.1 for n B ∼ 1-2n 0 [73].Nucleons are much heavier than pions and c 2 s is naturally small because of the non-relativistic regime.Two-and three-nucleon repulsions increase c 2 s , but their effects basically enter as powers of ∼ n 2 B and ∼ n 3 B whose growth are rather slow and c 2 s goes beyond 1/3 only at n B ≳ 2-3n 0 .This aspect differs from pionic matter in isospin QCD where pions can be relativistic already at ≃ n 0 and c 2 s ≥ 1/3 is achieved already at n I ∼ 2n 0 .

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.

A. Chiral symmetry restoration and softening
In Sec.III we have seen that larger f π and/or m σ 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 π and m σ mean the stronger chiral symmetry breaking in the QCD vacuum.
In the vacuum, the energy reduction due to the chiral symmetry breaking is (Fig. 11) where the first term is the energy of the trivial vacuum while the second one 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. 12).
A larger bag constant softens EOS at high density.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 µ = 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 The bag constant associated with the chiral restoration reduces the pressure and increases the energy density, resulting in a softer EoS at high density where the chiral symmetry is restored.Similar conclusions have been obtained in models with and without the U (1) A anomaly [74,75].
FIG. 11.The energy density as a function of the chiral effective mass Mq.After chiral restoration the minimum energy is realized at Mq = 0.In the broken phase the minimum is realized at Mq ̸ = 0 and the energy is smaller than that of Mq = 0.This gap in the zero-point energy density is the bag constant.

B. Power corrections to pQCD at high density
Our quark-meson EOS predicts c 2 s approaching 1/3 from above as density increases.This contradicts with the pQCD prediction in which c 2 s approaches 1/3 from below.A possible origin of such discrepancy would be the power corrections of ∼ Λ 2 QCD µ 2 I which cannot be derived from perturbative computations.
In Introduction, we schematically showed how power corrections can enhance the c 2 s in Eqs.( 1) and ( 2).The question is how large power corrections should be to qualitatively change the perturbative behaviors of c 2 s .For a given flavor f , the pQCD EOS up to O(α 2 s ) is given as [76] (we use the current quark mass, m u,d ≃ 5 MeV and µ I = µ u = −µ d in the present work) In the isospin symmetric limit, P u = P d .Here, P f 0 and P f 1 is the zeroth and first order in O(α s ), with and Λ reno being the renormalization scale.The running α s is FIG.13.Running coupling from pQCD (dashed) and from the freezing coupling (solid).For the renormalization scale Λreno = XµI , we examine X = 1, 2, and 4.
with L = 2 ln(Λ reno /Λ MS ), and Λ MS ≃ 340 MeV.The central value of Λ reno is Λ reno = µ f , and as usual we vary Λ reno from µ f to 4µ f .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 where t 1/2 low = 0.3 GeV and t 1/2 high = 1.1 GeV.For the low energy limit we use the form suggested by Deur et al [77,78] with α low s (0) ≃ 1.22 and κ ≃ 0.51.For the high density we use the perturbative expression (69), and for the intermediate region we use the interpolant 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. 13.
where c n 's are fixed by demanding the matching for n = 0, 1, 2. The six boundary conditions fix the six c n uniquely.Unlike in Refs.[77,78] 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 2 s .We set Q 2 = Λ 2 reno and plot α s (Λ reno ) in Fig. 13 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. 14).Now we add power corrections which are parametrized in terms of gaps in the pion condensed phase.The phase space factor ∼ 4πp 2 F ∆ times the gap ∆, divided by a factor (2π) 3 , yields the naive estimate where C is a constant of O(1).For our quark meson model C ≃ N c /2, see Eq. ( 57).In Son's estimate [79], based on the color-magnetic long range forces, the gap is evaluated as  We set C = 1 and examine pressure (Fig. 15) and c 2 s (Fig. 16) for ∆ = 0, 200, and 300 MeV for pQCD running coupling (dashed) and freezing coupling (solid).For ∆ ≃ 200 MeV, the power corrections are large enough for c 2 s to approach the conformal limit from above around ∼ 40n 0 .Meanwhile, at µ I ∼ 1 GeV or n I ∼ 40n 0 , parametrically the power corrections in EOS are corrections of the order thus ∼ 10% corrections.It is remarkable that even such small corrections can change the qualitative behaviors of c 2 s in the domain where pQCD seems applicable.

C. Trace Anomaly
Recently there has been growing interest in the trace anomaly in the context of mechanical properties in hadrons [72,[80][81][82] and in neutron stars [48][49][50].The latter is essentially the relation between P vs ε and is more fundamental than c 2 s which includes only the information of dP/dε, not the overall magnitude of P .In particular, Ref. [49] 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 where J µ D is the dilatation current, β < 0 the QCD beta function, and γ m the anomalous dimension of the quark mass.
For a hadronic state |K⟩ with the momentum K, the energy momentum tensor gives where the RHS is x-independent6 and does not contain g µν .The overall 1/m H factor is fixed by the condition at the rest frame, where we divide by ⟨K R |K R ⟩ to cancel the volume factor in the numerator.Thus, for a hadron at rest frame, we find with vanishing spatial components, ⟨K R |T ii (x)|K R ⟩ 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 K = 0.For instance an ideal Fermi gas leads to After lowering one index, we get ⟨T i i ⟩ < 0. Thus, the trace anomaly in a many-body system can be negative in principle.In thermodynamic systems, ⟨T µ µ ⟩ = ε − 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 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), but now we include the running of coefficients a 0 and a 2 associated with α s (µ I ).The energy density can be computed as and the trace anomaly is The running of α 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ε, which should not be confused with the BCS gap ∆.Shown in Fig. 17 is the ∆ tr as functions of n I for several calculations, our quark meson model (QM) and pQCD results for the renormalization scales with X = 1, 2, and 4. For the pQCD, both the perturbative (dashed) and IR freezing (solid) couplings are examined.Without power corrections the ∆ tr are all positive.In Fig. 18, we fix X = 2 for these two couplings, and the dependence of ∆ tr on the power corrections.The ∆ = 0, 200, 300, and 400 MeV cases are shown.With power corrections the ∆ tr appears to be negative, as we expected.Our QM model predicts ∆ ≃ 300 MeV and the negative trace anomaly for wide range.Finally we make a comparison between our QM model results and the lattice results in Ref. [63], as shown in Fig. 19.The QM model seems to capture the overall trend of the lattice data.Since the pQCD corrections and bag constant favor the positive ∆ tr , the negative ∆ tr may be taken as an indicator of the substantial power corrections.Trace anomaly ∆tr = 1/3 − P/ε as functions of µI /mπ for the quark-meson model with mπ = 140, 170 MeV and lattice results of Ref. [63] with mπ ≃ 170 MeV for different lattice spacing and volume.

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, tempering meson fields compared to the tree level amplitudes.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 going beyond the conformal limit c 2 s = 1/3 at ∼ 2n 0 , and making a peak at ∼ 5n 0 , at densities substantially smaller than the density for pions to spatially overlap, ∼ 10n 0 .The quark occupation probability at p = 0, f (p = 0), is ∼ 0.6 at 2n 0 and ∼ 0.9 at 5n 0 .The sound velocity peak is located around 5n 0 where the quark states around p = 0 are almost fully saturated, and it makes sense to associate the sound velocity peak with the saturation of quark states.After the bulk part of the quark Fermi surface is established, the c 2 s approaches 1/3 as in the relativistic limit.
Our model shows that c 2 s approaches 1/3 from above, mainly due to the power corrections, ∼ µ 2 I ∆ 2 ∼ µ 2 I Λ 2 QCD .This sort of terms is not available in pQCD calculations which predict c 2 s approaching 1/3 from below.Which one, perturbative or power corrections, dominates in c 2 s around ∼ 40n 0 is a quantitative question.The existence of the power corrections is related to the nonperturbative 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.Lattice results of Ref. [63] seem to support the negative trace anomaly in the domain between the BEC and the pQCD domains.Since the presence of non-perturbative effects in quark matter is a fundamental question, further clarifications by several lattice calculations with different systematics are highly desired to establish the findings in Ref. [63].
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.[83,84]).Including thermal effects into quarkmeson 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 quarkhadron-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. [85].The magnitude of thermal corrections is much smaller than the cold matter part due to ∼ (T /p F ) 2 suppression factors, but it can be important for NSs about to collapse, e.g., those appearing in NS-NS mergers [86][87][88].
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 ≃ 5-40n 0 , with the help of general causality and thermodynamic stability conditions (e.g., see Refs.[53,89]).But as seen in our simple exercise in Sec.IV B, the power corrections of ∼ 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 ≃ 5-40n 0 .The present one-loop analyses of quark-meson models should be also improved, using e.g., the functional renormalization group to include quark and meson fluctuations [90].
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 [91,92], 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. (B1) The function F (p 2 ) is given by where r = 4M 0 2 /p 2 − 1.This parametrization suggests that the parameters m π and m σ are restricted to m π , m σ < 2M 0 to make r real.
2  and b 2 ≃ 0. Hence, the µ 2I components of the effective potential are well-saturated by the V(2)q .
FIG. 1. Chiral and pion condensates as functions of a scaled chemical potential µI /mπ.

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

FIG. 3 .
FIG. 3.Isospin density nI as a function of µI , for the global (upper) and low density (lower) behaviors.The data points at low density are from Brandt et al.[62] (available up to µI ≃ 0.9mπ with mπ ≃ 140 MeV) and bands are from Abbott et al.[63] with mπ ≃ 170 MeV.(The definition of µI in Ref.[63] differs from ours by a factor 2 and this is taken into account in our figures.)

FIG. 4 .FIG. 5 .
FIG. 4. Pressure P as a function of µI for the low density (upper panel) and global (lower panel) behaviors.

FIG. 6 .
FIG.6.Pressure vs energy density from the tree and oneloop effective potential, for the low density(upper panel) and global (lower panel) behaviors.

FIG. 8 .
FIG. 8. (Upper) The occupation probability of u-and d quarks which corresponds to the residue of positive energy part of ⟨uu⟩.The densities we plotted are from 0.2 to 10.0n0 in 0.2n0 increments for gray curves and 1.0n0 increments for red curves.The blue dots are the locations of the surface of the distribution where f (p) has the half maximum.(Lower) The evolution of the occupation probability for the p = 0 state, f (p = 0).

FIG. 9 .
FIG. 9.The fit of the dipole function in Eq.(64) to the φ vac π

FIG. 12 .
FIG.12.The energy density with different strength of the chiral symmetry breaking.The bag constant is larger for the stronger chiral symmetry breaking.
µ u = -µ d [GeV] FIG. 14.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.13.

FIG. 15 . 2 . 2 µ
FIG.15.Perturbative pressure plus power corrections divided by the pressure in the Stefan-Boltzmann limit.The ∆ = 0, 200, and 300 MeV.The solid and dashed lines represent the freezing coupling and perturbative running with the Landau pole, respectively.The X for the Λreno is fixed to X = 2.

FIG. 17 .FIG. 18 .
FIG. 17. Trace anomaly ∆tr = 1/3 − P/ε as functions of nI for our quark-meson model (N f = 2 QM) and pQCD with the perturbative running coupling (dashed) and IR freezing coupling (solid).We vary X from 1 to 4. The trace anomaly is all positive in pQCD but negative for the quark meson model.