{'fraction_non_alphanumeric': 0.037756394586546366, 'fraction_numerical': 0.0141780928465077, 'mean_word_length': 4.066426480522405, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stiffening of matter in quark-hadron continuity

We discuss stiffening of matter in quark-hadron continuity. We introduce a model that relates quark wave functions in a baryon and the occupation probability of states for baryons and quarks in dense matter. In a dilute regime, the confined quarks contribute to the energy density through the masses of baryons, but do not directly contribute to the pressure; hence, the equations of state are very soft. This dilute regime continues until the low momentum states for quarks get saturated; this may happen even before baryons fully overlap, possibly at density slightly above the nuclear saturation density. After the saturation the pressure grows rapidly while changes in energy density are modest, producing a peak in the speed of sound. If we use baryonic descriptions for quark distributions near the Fermi surface, we reach a description similar to the quarkyonic matter model of McLerran and Reddy. With a simple adjustment of quark interactions to get the nucleon mass, our model becomes consistent with the constraints from 1.4-solar mass neutron stars, but the high density part is too soft to account for two-solar mass neutron stars. We delineate the relation between the saturation effects and short range interactions of quarks, suggesting interactions that leave low density equations of state unchanged but stiffen the high density part.

In spite of all these difficulties, the combined use of the above information and recent neutron star (NS) observations [23][24][25][26][27][28][29] allows us to get insights on the properties of dense matter (see, e.g., [30] for a short review). Recent analyses by NICER, including the radius measurements of 2.08M (M : solar mass) and 1.4M NSs, constraints from the NS merger event GW170817, and nuclear physics constraints, yielded the estimates of the radii R 2.08 R 1. 4 12.4 km [31,32]. This small variation in the radii from 1.4 to 2.08M NS suggests that the equation of state (EoS) for n B = 2 − 5n 0 should not contain substantial softening, but rather should get stiffened. This feature disfavors strong first order phase transitions in the domain n B = 2 − 5n 0 , although the weaker one is still possible.
Since the first discovery of 2M NS [25], a number of works have been devoted to the crossover description for hadron-to-quark phase transitions [33][34][35][36][37][38][39][40]. Early works [33][34][35] phenomenologically interpolated hadronic EoS at n B 2n 0 and quark matter EoS at n B 4n 0 ; the resulting EoS is consistent with the existing NS constraints, and has a novel peak structure in the speed of sound c s = (∂P/∂ε) 1/2 where P and ε are pressure and energy density, respectively. Later, such a peak structure was discussed as generic, by noting the contrast between stiffness of low and high density EoS [41]; nuclear physics calculations suggest soft low density EoS which must get stiffened rapidly to pass the 2M mass constraints [42,43]. This peak should have a mechanism specific to dense matter [44,45]; in the finite temperature crossover from a hadron resonance gas to a quark-gluon plasma, the speed of sound has a dip, instead of a peak, as one can see from the lattice simulations [46].
The microscopic mechanism for the emergence of the peak has been discussed by McLerran and Reddy (MR) [47], who used the concepts of quarkyonic matter [48]. The quarkyonic matter is a quark matter with a baryonic Fermi surface, and the excitations are confined [49][50][51][52][53][54][55][56][57][58]. In the MR model, they used a hybrid description in momentum space. A matter at low density is dominated by baryons, but as density increases, the quark Fermi sea emerges at low momenta, pushing up the baryonic states to high momenta. For suitable choices of parameters, baryons become relativistic at n B = 1.5 − 3n 0 with a substantial peak in c s . The advantage of the MR model is that relativistic baryons emerge by the quark Pauli blocking mechanism which is independent of details in nuclear forces, and hence the mechanism is qualitatively robust. Several successful descriptions of the NS have originated from this framework [59][60][61][62][63][64][65][66].
In this paper we discuss the stiffening of dense matter associated with the saturation of quark states at low momenta that may be regarded as the onset of the quark Fermi sea. The preliminary discussion was given in Ref. [67], and this paper is the fuller version. This work is basically a follow-up work of Ref. [47], but con-tains new attempts and insights which can be potentially important.
First, we describe the crossover behavior using quark degrees of freedom only, starting with the description of a single baryon, proceeding to a baryonic matter, and then to a quark matter formation. Although the descriptions are rather crude, this approach has the advantage over conventional hybrid descriptions where one uses quarks in one place and baryons (hadrons) in the other place. This removes the worries about double counting as well as the confusions associated with switching in degrees of freedom. 1 The onset density of quark matter formation is related to the size scale of a baryon, and we found that the saturation begins to occur slightly above the nuclear regime, 1 − 3n 0 , even before baryon cores overlap. The softness in baryonic matter and stiffness in quark matter are described in a unified manner. We also check how our quark descriptions are related to the MR model.
Second, we attempt to describe in-medium interactions at the level of quark descriptions. This is potentially important, considering the difficulties to predict the high density behaviors of two-, three-, and more-body baryonic forces within purely baryonic calculations. It has been known that a simple constituent quark model with one-gluon exchanges accounts for the baryon spectroscopy remarkably well [70,71]. Recent lattice calculations for baryon-baryon interactions even show that the simple constituent quark picture correctly describes the observed patterns of baryon-baryon interactions, including the hard core repulsion among nucleons as well as baryon interactions with strangeness [72][73][74]. These successes in describing semishort range correlations, at momentum scale ∼ 0.2 − 1 GeV, give us a hope to build a unified description for the properties of matter from the baryonic to quark matter regime. In this spirit, relevant interactions at ∼ 5 − 10n 0 , at lower density than for the perturbative regime, have been examined [36-38, 75, 76]. The present framework solely based on quarks is suitable for the unified treatments of interactions from low to high densities.
In Sec.II, we begin with quarks in a baryon, and in Sec.III we discuss quarks in a baryonic matter. In Sec.IV, quark matter formation and the associated stiffening are described. In Sec.V, we mention how our descriptions are related to the MR model. In Sec.VI, we discuss the spin-flavor quantum numbers of baryons and how they fill the quark spin-flavor states. In Sec.VII, we discuss interactions for matter from the dilute baryonic to dense quark regimes and examine what kind of interactions is suitable to describe the NS phenomenology. Section VIII 1 The double counting introduces serious problems into the field theoretic computations for the zero point energy in matter [68]. The UV divergences appear from both elementary and composite particles. Consistent treatments of both contributions are mandatory to cancel the divergences in medium by the vacuum subtraction of the energy [69].
is devoted to a summary. In the Appendix, we discuss the phase space density in baryonic descriptions for the quark matter domain.

II. QUARKS IN A BARYON
We consider how quark states are occupied as baryon density increases. We postulate 2 a distribution of quarks with the momentum p which belong to a baryon with the momentum P B ; the form is given by (N c : number of colors) where Λ is the scale of the order of the QCD nonperturbative scale, Λ QCD 0.2 − 0.3 GeV, related to the radius of a baryon as Λ ∼ R −1 baryon , and N is the normalization constant . The momentum distribution becomes broader for a smaller baryon radius.
We emphasize that Q in is for a quark with a particular color, flavor, and spin. When we compute baryonic quantities we have to sum contributions from all colors. The sum of N c quark momenta should be P B , and the distribution satisfies this condition at the level of the expectation value, More realistically the probability distribution of a quark must be related to those of the other N c − 1 quarks, but in this paper we consider only the averaged description for quark momentum distributions. As the quark momentum has the broad distribution, the magnitude of momentum is substantial, independent of P B , and hence quarks can be energetic compared to the baryon kinetic energy or the nuclear scale of O(1−10) MeV. In baryonic matter, this disparity will be reflected as the large energy density but small pressure.
2 For a constituent quark model with Nc = 3 and the harmonic oscillator potential [70,71], we can derive Q in within simple analytic calculations. First we compute three quark wave functions Ψ(p 1 , p 2 , p 3 ; P B ) and the corresponding probability function |Ψ(p 1 , p 2 , p 3 ; P B )| 2 . Then we integrate out two of momentum variables to get the single particle distribution Q in (p; P B ) in Eq.(1).
Below we approximate the baryon energy as the sum of the energies from N c quarks. The energy of a quark in a baryon with momentum P B is where the term linear in P B vanishes. It is important to note that the correction from finite P B is suppressed by 1/N 2 c ∼ 1/10, and the quark single particle energy is hardly affected by the baryon momentum until P B becomes very large, ∼ N c Λ.
The simplest version of our model includes only a potential localizing quarks, and the resulting baryon energy is simply For instance, in a nonrelativistic quark model and the kinetic energy of the baryon is suppressed by the large baryon mass, It is important to note that M B in our model is considerably larger than N c M q by the kinetic energy ∼ Λ of each quark. If we use usual constituent quark mass of M q ∼ 0.3 GeV for up-and down-quarks, there must be attraction of ∼ Λ to keep the picture of M B ∼ N c M q . In Sec.VII, we consider such short range interactions to modify the average single particle energy and to get the right baryon mass.

A. Occupation probability of quark states
We first write the occupation probability of quark states, f q (p; n B ), with given color, flavor, and spin as where we wrote the occupation probability of baryon states as B(P B ; n B ) for given flavor and spin 3 . The expression means that f q in dense matter is obtained by summing up the quark occupation probability from each baryonic state. Being probabilities, 0 ≤ f q ≤ 1 and 3 For the moment we treat quarks and baryons as if they have no flavor and spin species, or discuss u ↑ state in a ∆ ++ sz =3/2 = (u R ↑ , u G ↑, u B ↑) baryon. More details about spin-flavor quantum numbers will be addressed in Sec.VI. 0 ≤ B ≤ 1 must be satisfied. Note that we have assumed that the probability depends only on the size of momenta, p = |p| and P B = |P B |.
Integrating the quark momentum p with the normalization p Q in = 1, we find (9) where n R,G,B q is the quark density for a given color. (We implicitly assumed the color neutrality condition n R q = n G q = n B q . The quark number density is n q = N c n B .) In dilute regime, we can neglect interactions among baryons as they are widely separated in space. Then baryons fill the states from low momenta with the probability 1, as in an ideal gas; (N f : number of flavors) where θ(x) is the step function, P F the Fermi momentum of baryons, and 2N f the factor from spins and flavors. Let us look at how f q evolves as n B increases. In the following, we rescale momenta, (11) with which f q can be written as We note that, in a dilute regime,P F 1 or P F N c Λ, so that the domain for the integral overP B is very small, canceling the overall factor (N c Λ) 3 . AsP B in the exponent of the integrand can be regarded as small, we expand them and find Recalling N ∼ Λ −3 , the overall size of f q is ∼ (P F /Λ) 3 . As for the shape, the leading order contribution maintains the Gaussian form same as in Q in . The correction starts with the order of 1/N 2 c 1/10 for N c = 3, so the large N c should be a good approximation for the N c = 3 case, except in situations where p and P F are very large. The numerical results of Eq.(12) for various n B /n 0 are shown in Fig.1, together with the leading contribution of the 1/N c expansion. We took Λ = 0.25 GeV and N c = 3.
Here, we emphasize that we used the formula (12) for only the curves up to n B /n 0 = 1.25; for n B /n 0 > 1.25, we used another expression Eq. (20) which is explained in the next section. Here, f q for various n B /n 0 are shown for Λ = 0.25 GeV and N c = 3.
It is remarkable that the shape of f q hardly changes, while the size grows almost linearly in n B . Since f q is a probability, the growth in the size of f q must be terminated at f q = 1. We call it saturation of quark states, or more simply, quark saturation. The behavior of f q beyond n B /n 0 = 1.25 is discussed in the next section. In the large N c limit, one can derive a number of simple expressions as we may neglectP F terms in Eq. (13). The f q takes the form where the p = 0 state gives the largest f q . The p = 0 state gets saturated at n c B = N f (P c F ) 3 /3π 2 ; the corresponding baryon Fermi momentum is In Fig.2 the f q (p = 0; n B ) is shown as a function of n B for N c = 3. For a larger Λ, the baryon radius is smaller and the saturation of the p = 0 state happens at larger n B . It is important to note that, for Λ = 0.2−0.3 MeV as the reasonable scale of baryon radii, the saturation takes place at n B = 0.5 − 2n 0 ; this baryon density is within or close to the territory of conventional nuclear physics. Such density range is well below the density ∼ 4 − 7n 0 [37] where baryon cores begin to overlap. The reason why this may happen would be that quark wave functions, in general, have broader extension than the average baryon size, or the overlap of the meson cloud around a baryon take place before the cores overlap [77].

B. Equations of state in dilute baryonic matter
Next we compute the EoS using the occupation probability discussed in the previous section. For simplicity we assume the quark energy of the form The energy density is computed as The chemical potential is computed as The expression is particularly simple in the large N c limit; using Eq. (14), we find so that µ B = M B at large N c . Accordingly the pressure is small, as expected from purely baryonic matter in an ideal gas regime. This trend changes when quark states at low momenta get saturated, as we see in the next section.

IV. QUARK MATTER FORMATION
We have seen that the quark states at low momenta begin to be saturated at n c B ∼ Λ 3 . Beyond this critical density we can no longer use the ideal gas description of baryons. In fact, the EoS behaves very differently before and after the saturation; EoS are much stiffer after the saturation.
It is highly unconventional to discuss the baryon momentum distribution B after the quark Fermi sea begins to be relevant. In this situation, instead of starting with B, it is more intuitive to postulate the form of quark occupation probability f q for which we can implement the quark Pauli blocking easily. In Sec.V, we come back to the question of how B looks like for the postulated form of f q .

A. A model of saturation
We assume that the occupation probability changes smoothly just after the saturation takes place. As a trial, we postulate the form (20) where p sat is a function of n B . The first term is responsible for states occupied with the probability 1; the states to p = p sat are saturated. Meanwhile the second term is for states with p > p sat which are only partially occupied. Here we assume that the distribution f q at n c B is shifted by the occupied level, so the states up to ∼ p sat + Λ can be occupied with substantial probability. For Λ → 0, the postulated form of f q is reduced to the form for the ideal quark gas. But we keep Λ finite, assuming that quarks just after the saturation are not fully delocalized. The behavior of f q after the saturation is shown in Fig.1 for While the postulated form of f q seems a small departure from f q in a baryonic matter, it has dramatic impacts on the EoS. First, we note that the baryon densities before and after the saturation are continuous, as we postulate the continuous changes in f q . The relation between the Fermi momentum and p sat is given by Near the saturation, we keep only terms to O(p sat ). or where Λ c ∼ Λ characterizes the thickness of the distribution at n B = n c B , In Similarly ε is continuous before and after the saturation, where ε sat is the contribution from p = 0 to p sat . For a small p sat ∼ 0, we neglect ε sat and expand the integrand in the second term, where ε c is the energy density at n B = n c B . The derivative is ∂f q /∂p 2 ≤ 0, so the ε approaches ε c continuously from above for p sat → 0 + . Noting ∂p sat /∂n B π 2 /N f Λ 2 c , the chemical potential just after the saturation is which is compared to the µ B before the saturation. A number of analytic insights are obtained in the large where we wroteẼ = E/Λ. The chemical potential before the saturation is given by the baryon mass at large N c , If the saturation takes place in the relativistic regime of quarks, we may expand E q ∼ p + · · · , and find µ after where the chemical potential jumps by 0.18M B . Meanwhile, in the nonrelativistic limit, where the jump in µ B is the order of nonrelativistic corrections. We found that, while n B and ε do not contain any jumps, the derivatives do. Of course, the thermodynamics does not allow jumps in µ B and the results being presented must contain something unphysical, in spite of the seemingly reasonable form of f q in Eq. (20). But for the moment we proceed further to examine what would happen in this idealized description. Figure 3 shows µ B as a function of n B /n 0 where we chose N c = 3, M q = 0.3 GeV, and Λ = 0.25 GeV, for which M B 1.26 GeV, and the jump in µ B associated with the saturation is 0.1 GeV. It is clear that the µ B and ε/n B after the saturation grow much faster than the behavior before the saturation. For comparisons, we also show the results for the Quark-Hadron-Crossover (QHC)19-D EoS [38] as an example of EoS consistent with NS observations; for the QHC19-D, M max 2.28M , R 1.4 11.6 km, and R 2.08 11.5 km. (In Sec.VII, we make comparisons again after adjusting the baryon mass.) These jumps in µ B result in the discontinuities in the corresponding pressure P through the thermodynamic relation P = µ B n N − ε, where n B and ε are continuous but µ B are not. Now, the pressure just after the saturation is   which is much bigger than the pressure of an ideal baryon gas, P ∼ n Accordingly the squared speed of sound c 2 s = ∂P/∂ε diverges to +∞ at the saturation. The P vs ε for the setup, same as Fig.3, is shown in Fig.4.

B. Smoothing out the discontinuities
The discontinuities are the artifacts which are presumably associated with our use of the ideal baryon gas picture for baryons just before the quark saturation (Fig.  5). The validity of the ideal gas picture should break down before reaching the quark saturation point. The baryon interactions are mediated by the meson exchanges or quark exchanges, so near the quark saturation point those interactions should grow up. In realistic nuclear EoS, nuclear interactions of O(N c ) are important, and the pressure increases faster than in the ideal baryon gas.

FIG. 5.
Smoothing out the artifact: The use of ideal baryon gas description to the quark saturation point leads to a discontinuous jump in P (left panel). With baryon interactions mediated by quark exchanges, the precursory behavior should appear before reaching the quark saturation point (right panel). In the lower panels, the evolution of fq, which changes the direction from the vertical to horizontal direction at the quark saturation point, takes place more smoothly with baryon interactions. Thus, with interactions the pressure does not suddenly jump up, but there should be a precursory behavior as the system approaches the quark saturation regime.
With these considerations in mind, we introduce an ad hoc procedure to smooth out the discontinuities. We revise the previous model of quark saturation slightly, by multiplying a smearing function, where p w is a regulator which should be taken to be very small, p w Λ. For p sat p w , p sat in the step function and the Gaussian factor is negligible, so we can regard in computations of thermodynamic quantities. Technically, the factor p sat /p w 1 plays the role similar to the factor ∼ N P 3 F ∼ n B /n c B in Eq. (13). There the shape of f q was fixed but its magnitude increases until the saturation takes place. When p sat ∼ p w , the model goes back to the model in Eq.(20) moduloP F ∼ 1/N c corrections. The n B dependence of f sat q for p w = 0.01 GeV is shown in Fig.6. The qualitative behaviors are very similar to Fig.1 [the major difference comes from our neglect ofP F corrections in Eq.(33) compared to Eq. (20)].
Now, we use f sat q to examine several thermodynamic quantities. Shown in Fig.7 is µ B as a function of n B . Working with a model with the same form from low to high density, the artificial discontinuities found in the previous treatment are smoothed out, and µ B is now continuous, as it should. Accordingly, the (squared) speed of sound c 2 s is now regulated and well defined, as shown in Fig.8. The c 2 s has a peak around the density where the saturation effects become important, and it exceeds the conformal value 1/3 for a sufficiently small p w . The degree of the smearing is rather sensitive to the value of p w , and a too small p w tends to violate the causality constraint c 2 s ≤ 1. The sharpness of the peak should be dynamically determined by the interplay between baryon and quark dynamics. We choose p w = 0.01 − 0.02 GeV unless otherwise stated.
We emphasize that the existence of the peak was not driven by nuclear forces, but by the quark Pauli blocking effects. Including baryon interactions tames the peak, rather than enhancing it; stronger baryonic interactions before the quark saturation blur the peak structure, by reducing the contrast between baryonic and quark matter pressures. Taking the picture that baryon interactions are mediated by quark exchanges, baryonic matter with stronger interactions should look more similar to quark matter 4 . This viewpoint motivates a three-window modeling of dense matter from nuclear to quark matter domains [33][34][35][36][37][38][39][40], and more microscopic considerations based on mode-by-mode percolation [77].

V. BARYONS AFTER SATURATION
As we have seen, the saturation of low momentum quark states induces a rapid increase in the pressure. Such changes are difficult to imagine from purely baryonic descriptions, and we are interested in how the corresponding occupation probability of baryon states (B)  looks. Through attempts to understand B, we also discuss how to obtain expressions similar to the MR model for quarkyonic matter EoS [47].
Regarding (for a given n B ) f q and B as vectors with indices p and P B , respectively, they are related through a matrix Q in as where 0 ≤ B(P B ; n B ) ≤ 1 for any P B and n B . The equation is linear and in principle one can take inversion to determine B for a given f q . But in practice this method does not work well. Another strategy is to prepare some model B α for B with parameters α = (α 1 , α 2 , · · · ), and minimize a function "energy" functional where I cost is some functional which gives the energetic penalty if, for instance, B α violates the condition 0 ≤ the saturation takes place. The B sh convoluted with Q in leads to the quark momentum distribution f sh q around psat. For the "baryon shell" of B sh , we found h∆ ∼ 1/N 2 c for P sh ∼ NcΛ. If we fix h = 1, the thickness ∆ gets narrower as in the MR model. B(P B ; n B ) ≤ 1. For this strategy to work, we need a good guess for the form of B.
For this purpose, we first assume baryons with momenta P B Λ and discuss which form of f q should follow. We consider the form where "sh" is the abbreviation of "shell", and h is the overall size in the occupation probability (0 ≤ h ≤ 1) for states with momenta P sh − ∆(≥ 0) to P sh . Below, to simplify calculations, we assume ∆ to be small, ∆ P sh . Integrating the baryon momentum distribution, we get the contribution to f q as where we setP sh = P sh /N c Λ. For smallpP sh , and for largepP sh , where the maximum of f q appears atp ∼P sh . Below, we examine these asymptotic behaviors in some details. We can see a number of the remarkable features in the large N c limit. For P sh ∼ Λ orP sh ∼ 1/N c , we use Eq.(39) to obtain where the condition f q ≤ 1 demands h∆ ∼ Λ. For P sh ∼ N c Λ orP sh ∼ 1, we use Eq.(40) to obtain where we setP sat /p ∼ 1 as the Gaussian part has a peak atp =P sh . In this regime the condition f q ≤ 1 demands h∆ ∼ Λ/N 2 c . Now we consider the forms of h that are compatible with the scaling behaviors in Eqs. (41) and (42). For example, one can take where c 0 , c 1 , and c 2 are constants of O(1). For P sh ∼ Λ, the first term dominates. For P sh ∼ N c Λ, all these three terms can be comparable. In what follows, baryon states with large momenta P sh are compatible with the condition f q ≤ 1 only if those states are occupied with the small probability density, i.e., either h or ∆ becomes small for large P sh . In particular baryon states with momenta ∼ N c Λ are occupied with a small but nonzero probability h∆ ∼ 1/N 2 c . In the Appendix, we further explore the phase space density in baryonic descriptions.
Now we notice that the scaling f sh q ∼ e −(p−P sh ) 2 is similar to the behavior of f q (p − p sat ; n c B ) ∼ e −(p−psat) 2 in Eq. (20) for the partially occupied quark states beyond the saturated states. This suggests thatP sh p sat or P sh N c p sat , and the approximate relation Thus the shell form Eq.(37) turns out to be a good candidate for the p > p sat part of the quark distribution f q postulated in Eq. (20). The schematic picture is given in Fig.9. Substituting this form into Eq. (20), and then integrating quark momenta p, we reach a model similar to the model of McLerran and Reddy [47]. 5 With P sh = N c p sat , the baryon number and energy densities per flavor are where we have used the relations p Q in = 1 and E B (P B ) = N c p Q in (p, P B )E q (p), and neglected possible double counting in f q at p p sat which should be minor effects unless p sat is very large. We also note that 5 Our presentation here is slightly different from the paper of McLerran and Reddy [47] where the concept of the occupation probability for baryons is not used, so h = 1 from the very beginning, while they assumed the form of ∆ to be The second term plays important roles to regulate µ B and the speed of sound.

FIG. 10.
Examples of quantum numbers for baryons where ∆ ++ sz =+3/2 and Ω − sz =−3/2 are shown. The former (latter) fills the u ↑ (d ↓) states for all colors RGB, while leaving the other color-flavor-spin states empty. For a given space wave function for quarks, six baryon states completely fill color-flavor-spin states.
where the contributions from the saturated quark states are suppressed by an extra factor of 1/N c . When P sh ∼ Λ, baryons are nonrelativistic, What is remarkable is that the relativistic regime is reached already for n B ∼ Λ 3 . For P sh ∼ N c Λ, we obtain which is ∼ Λ 3 . For small P sh , the first term dominates, while the second and third terms slowly grow as P sh becomes ∼ N c Λ. This result drastically differs from naive expectation n B ∼ (N c Λ) 3 . In ordinary baryonic descriptions at n B ∼ Λ 3 , baryons approach the relativistic regime by interactions of O(N c ), but here we found that relativistic baryons can also emerge just due to the constraints on baryons as composite particles. This reproduces the remarkable findings by McLerran and Reddy [47].

VI. QUANTUM NUMBERS
We discuss the quantum numbers of baryons such as colors, flavors, and spins. The question is how baryons saturate these quantum numbers for quark states. In previous sections, we implicitly assumed that the quark states for N c colors, N f flavors, and two spins are saturated by putting N f species of baryons with two spins. In this section we explain why this description should be valid. Specifically, we consider the N c = N f = 3 cases to utilize the terminology common in hadron spectroscopy. We also ignore the mass difference among up, down and strange quarks.
The argument is simplified by assuming SU (2N f ) symmetry in which there are no energy differences associated with spins and flavors. In particular nucleons, ∆, Ω, so on, are all energetically degenerate. It is useful to focus on the following states, where quark states are totally antisymmetrized.
While there are many more S-wave baryons with spin or flavor excitations (such as baryon octet or decuplet), for given spatial wave functions we can saturate quark states using only the 2N f (= 6) states listed in Eq. (49). For example (Fig.10), by putting ∆ ++ sz=3/2 , colors RGB for up quarks with spins aligned in ↑ directions are filled at once. Therefore the number of baryon species we need is the same as the number of quark species, as we have assumed in the previous sections. The arguments are applicable to any N c .
This discussion also suggests that, once baryons with a specific spin-flavor quanta form their Fermi sea, the other baryons cannot freely enter the system due to the Pauli blocking at quark level. For instance, when ∆ ++ sz=3/2 and ∆ ++ sz=−3/2 form the large Fermi sea of up quarks, one cannot put nucleons (uud or ddu) at low energy, but must place them on top of the saturated Fermi sea. This viewpoint should be important when we describe hyperons entering dense nuclear matter. Also, this consideration should affect the treatment in the self-energy processes for protons and neutrons, as virtual states are blocked by preoccupied states.
More realistically we need to discuss the cases where protons and neutrons appear as the lowest energy states. The treatments of quantum numbers are more involved than the idealized (∆ ++ , ∆ − , Ω − ) baryon bases used here. We leave such realistic cases for our future work.

VII. INTERACTIONS
So far our discussions on EoS are entirely based on quasiparticle pictures, based on either quarks or baryons; the interactions for quarks have been taken into account only indirectly, by demanding that quark momenta are  distributed to some range as they should localize by confining effects. Now we perturb the calculations presented in the previous sections. We discuss interactions within quark descriptions.
In principle, one can construct EoS from the information of single particle propagators. The pressure at given µ B is given by where n B can be expressed as (γ 0 : the zero component of the Dirac matrices) where S q is the quark propagator and the trace runs over all quantum numbers. This relation is exact for whatever interactions; for example, in functional frameworks such as the 2PI action [78][79][80], Eq.(51) always follows from self-consistent treatments of the quark self-energy and interactions. As a result, the effects of interactions can be included into the self-energy of the propagator. If we need to include the baryon-baryon interactions, one should write the corresponding 2PI graphs and consider all possible cuts of quark propagators to generate the self-energy graph, see Fig.11 for a quark in a baryon and Fig.12 for a quark participating in baryon-baryon interactions.
In this paper, we simply assume a phenomenological parametrization for a single quark energy. We consider the form for a single particle energy, where V is the contribution from interactions that may depend on f q . Our first task is to adjust the baryon mass. As in a quark model with a phenomenological confining potential [70,71], our baryon mass is too massive due to the kinetic energies of localized quarks. To reproduce the observed baryon spectra, we need color-electric interactions for the 11.6 km, and R 2.08 11.5 km.) overall reduction in the masses for baryons, and colormagnetic interactions to get the correct mass splitting, e.g., the N -∆ and π-ρ splittings. In this work, we focus on the electric interaction but neglect the details of mass splittings. In vacuum, we consider where A indicates an antisymmetric representation in colors. We adjust V vac CE to set the baryon mass to the nucleon mass, M N 939 MeV.
In a baryon, the color-electric forces always reduce the average quark energy, as the color wave function is always antisymmetric for any combinations of quarks (e.g., Fig.11). For symmetric wave functions, the color-electric forces yield repulsive forces, and such channels are inevitable when several baryons come close together (e.g., Fig.12).
To take into account these attractive forces in the dilute regime and the repulsive forces in denser regime, we consider the parametrization (C S E ≥ 0) where S indicates symmetric channels in colors. The power β controls the impacts of partially occupied levels, and the parameter γ is chosen to either 0 or 1 to examine the effects of partially filled states. In this model, quarks with p p sat are free from the saturated levels and have the energy reduction (f q 1), as quarks in a baryon. This feature may be interpreted as the attractive correlations near the Fermi surface; since the quark states are not fully occupied, quarks can arrange their wave functions to enhance the portion of color-antisymmetric channels as in an isolated baryon. Meanwhile, quarks with p p sat have less freedom for such arrangements, and feel the overall repulsion (f q 1 for the saturated levels), due to the appearance of the repulsive channels. The repulsive energy is activated only when p sat becomes substantial. The power β determines how sharply one can distinguish the filled (f q 1) and partially filled (f q < 1) states. Now, we put these ingredients into our numerical analyses. First we set C S E = γ = 0. Shown in Fig.13 are various EoS, n B vs µ B , P vs ε, and c 2 s vs n B , together with the results of QHC-D as a guideline for EoS consistent with the NS observations [38]. We set C A E = 0.15 GeV to get M B 0.94 GeV at n B = 0. The EoS at low density behaves as those in QHC19-D, but starts to deviate from QHC19 around n B ∼ 3n 0 due to softening of our model at high density. The speed of sound has the maximum at ∼ 2.6n 0 , mainly determined by our model parameter Λ (or baryon size). We emphasize that stiffening takes place in our model solely by the saturation effects.
To get stiffer EoS at high density, we set γ = 1 to activate the effects of partially filled states. With β = 1, the low density part is much stiffer than the baryonic part in the QHC19 (which is the Togashi EoS [15]), contradicting with nuclear EoS at n B n 0 . Here, quarks lost the attractive energy too rapidly, leading to rapid growth in µ B as n B increases. One can delay the loss of attractive energies by increasing β and suppressing f β q terms for f q < 1. We found that, for β 8, our EoS gets along with the nuclear EoS for n B 1.5n 0 . Meanwhile, the high density part is considerably softer than QHC19-D.
In order to make the high density part stiffer, the simple way is to turn on C S E . We increase from C S E = 0.05 GeV to 0.20 GeV, and the corresponding results are shown in Fig.14. The repulsive forces are activated only when f q 1, or p sat is large, so leaving the low density part as before, but stiffening the high density part where p sat is substantial. The location of the maximum in the speed of sound is not very sensitive to the choice of C S E , but the height becomes larger for a greater C S E . We also note that the width of the peak is sensitive to our choice of p w , as in Fig.6.

VIII. SUMMARY
In this paper, we discussed how quark degrees of freedom stiffen EoS. In order to relate the quark dynamics for a single baryon to baryonic matter and quark matter formation, we have introduced a model that relates three relevant functions: quark momentum distribution Q in in a baryon, occupation probability of states for baryons B, and occupation probability of quark states f q . We also consider the effects of interactions at the quark level.
Below we summarize our findings: (i) In a dilute regime, the confined quarks contribute to the energy density through the mass of baryons, but do not directly contribute to the pressure; hence, the EoS are very soft (Fig.15 left). This dilute regime continues until the low momentum states for quarks are saturated (Fig.15 middle). The saturation can take place considerably before the baryons fully overlap, possibly at density close to the nuclear saturation density, n B ∼ 1 − 3n 0 . This picture is in line with the recent proposal of Soft Deconfinement as the onset of the mode-by-mode percolation [77].
(ii) After the saturation, the energy per particle, ε/n B , begins to change as in quark matter (Fig.15 right), and the pressure, P = n 2 B ∂(ε/n B )/∂n B , grows rapidly, although changes in n B and ε are modest. These features lead to a peak in speed of sound, c 2 s = ∂P/∂ε. In our model, such a peak follows by just assuming the continuity of f q before and after the saturation, while no detailed descriptions of interactions were necessary.
(iii) To take into account the constraint f q ≤ 1 after the saturation, it is easiest to directly work with the quark description. But we can also infer the baryon momentum distribution B consistent with the desired form of f q . Through such attempts, we reached a model similar to the MR model for quarkyonic matter. The resulting B differs from the pure baryonic descriptions in which baryons are treated as if elementary particles; baryons after the saturation are highly relativistic.
(iv) We do not need many baryon species to fill the quark levels. For three flavors, we need 2N f baryon states to fill the quark color-flavor-spin states for a given spatial wave function. In this respect, nucleons and hyperons should not be treated as independent when they share the same quark states [61,81].
(v) While the stiffening of matter near the saturation seems a generic trend in our modeling, a model without interactions does not lead to sufficient stiffness at high density, n B 4−5n 0 , that is required by the existence of 2M NSs. This observation is consistent with viewpoints in our previous works [38,75,82] where various shortrange interactions of p = 0.2−1 GeV were discussed. One way to stiffen the high density part is to use a model in which quarks in saturated states feel repulsions but those near the Fermi surface enjoy the attractive correlations.
Unfortunately, our discussions remain largely qualitative and the treatments of dynamics are in many senses ad hoc. There remain many things to be done, as listed below.
First, it is better to directly use a quark wave function in a constituent quark model. To get Q in in this work, we should first calculate the three-body wave function and then integrate out two momentum variables. By doing this the matter properties can be directly expressed by quantities in hadron spectroscopy.
Second, to use the framework for predictions of the NS properties, we need to include the flavor asymmetry, and, in particular, have to discuss how nucleons fill the quark states at low energy.
Third, as we use a quark model, it is desired to directly use interactions in a quark model for baryons. There have been many works to reproduce baryon properties, and the lattice QCD studies of baryon-baryon interactions support the idea that the short range part (such as the hard core repulsion among nucleons) is overall consistent with the descriptions based on quark dynamics with onegluon exchanges.
Fourth, pairing effects leading to the chiral condensates or diquark condensates should be discussed to determine the phase structures as well as EoS. In this work, we fixed the constituent quark mass as in vacuum, but it is very likely the effective mass changes with density. Finally, the model should be extended to finite temperatures. When we come to thermal excitations, we need to address whether thermal excitations are hadrons or quarks. In fact, this is a crucial step to establish the quark-hadron continuity or the quarkyonic matter scenario, because the response in EoS as well as in the transports to changes in temperature is entirely different for hadronic and quark excitations [83]. In two-color QCD, there is a hint from the lattice QCD for hadronic excitations at high density [84,85], although more data are needed to establish this idea.
In forthcoming papers we plan to give more quantitative estimates on EoS, arranging the setup for the baryon spectra and NS phenomenology.

FIG. 15.
A schematic picture for the evolution of matter from the baryonic to quark matter regime, see also Fig. 12 in Ref. [77]. In baryonic matter, the energy per particle is ε/n B const., and after the saturation, ε/n B starts to grow as in quark matter. As the pressure is given by P = n 2 B ∂(ε/n B )/∂n B , matter after the saturation has much larger pressure than baryonic matter. The speed of sound makes a peak around the density where the saturation happens. excitations in quark matter. This work is supported by NSFC Grant No. 11875144.

Appendix A: Phase space density in baryonic bases
In Sec.V, we considered the shell distribution of baryonic states after the quark saturation and have derived a constraint h∆ Λ/N 2 c . The constraint acts on the product h∆, not on each of h and ∆ separately. In the MR model, a thin shell structure with h = 1 and ∆ ∼ Λ/N 2 c is used. Here we consider another possible choice. This illustrates the peculiar structures of the baryon occupation probability after the quark saturation.
Let us discuss the increase of baryon number in two different bases, i.e., quark and baryon bases, see Fig.16. Noting the relation n B = n R,G,B q = n q /N c , we compare the change of quark Fermi sea for a specific color with that of the baryon Fermi sea.
In quark descriptions, we increase the quark Fermi momentum from p to p + δp. The density increases as where N c factors do not show up. Now we try to describe the same increase in density using purely baryonic bases. We note that the quark and baryon momenta near the Fermi surface are related as P B ∼ N c p. Similarly δP B ∼ N c δp. Then, naively one would reach which does not match with δn R,G,B q , due to the factor N 3 c . In order to avoid this contradiction, we have to conclude that the phase space for baryons is more dilute than the quark's phase space. The proper estimate should take into account the occupation probability h for baryonic states [67], In terms of the h∆ constraint discussed before, here we are thinking of the case that h ∼ 1/N 3 c and ∆ ∼ N c δp. Within this description, we consider the variation of energy density with respect to the change δp, and examine how the expression appears consistent for quark and baryon descriptions. For simplicity, we assume the case where p Λ so that quarks and baryons have the energies ∼ p and ∼ N c p, respectively.
In quark bases, we need to sum contributions from all colored quarks, δε quark ∼ N c × p 3 δp . (A4) In baryon bases, we take into account the baryon occupation probability h ∼ 1/N 3 c and get So in both bases the variation of the energy density is ∼ N c p 3 δp, as it should. This represents the dual feature of quark and baryon descriptions; the concept of occupation probability for baryonic states is the key to satisfy the consistency condition. The remarkable feature in baryonic descriptions is that while baryons occupy the high energy states with small probability, those states still make significant contributions to the baryon number and energy density. This feature is difficult to be foreseen from pure baryonic considerations; the quark saturation effects demand highly exotic descriptions for baryons.

FIG. 16.
The comparison of the occupation probabilities for quark and baryon states. We increase the quark Fermi momentum p by δp, with which the quark density for a given color (and spinflavor) increase by δn R,G,B q ∼ p 2 δp. In terms of the purely baryonic bases, the baryon number must increase by the same amount, δn B = δn R,G,B q . In order to satisfy this consistency condition, we have to conclude that the phase space in baryonic bases is more dilute than the quark's phase space.