Hadron-Quark Crossover and Massive Hybrid Stars

On the basis of the percolation picture from the hadronic phase with hyperons to the quark phase with strangeness, we construct a new equation of state (EOS) with the pressure interpolated as a function of the baryon density. The maximum mass of neutron stars can exceed $2M_{\odot}$ if the following two conditions are satisfied; (i) the crossover from the hadronic matter to the quark matter takes place at around three times the normal nuclear matter density, and (ii) the quark matter is strongly interacting in the crossover region and has stiff equation of state. This is in contrast to the conventional approach assuming the first order phase transition in which the EOS becomes always soft due to the presence of the quark matter at high density. Although the choice of the hadronic EOS does not affect the above conclusion on the maximum mass, the three-body force among nucleons and hyperons plays an essential role for the onset of the hyperon mixing and the cooling of neutron stars.

matter has been assumed to be a first-order phase transition and Gibbs phase equilibrium conditions are imposed. However, treating the point-like hadron as an independent degree of freedom is not fully justified in the transition region because all hadrons are extended objects composed of quarks and gluons. Furthermore, the system must be strongly interacting in the transition region, so that it can be described neither by an extrapolation of the hadronic EOS from the low-density side nor by an extrapolation of the quark EOS from the high-density side [24]. This is analogous to the BEC-BCS crossover realized in the many-body system of ultra-cold fermionic atoms [25]. Figure 1 illustrates the above situation in terms of the pressure as a function of baryon density (ρ). One may expect a gradual onset of quark degrees of freedom in dense matter associated with the percolation of finite-size hadrons, i.e., a smooth crossover from hadronic matter to quark matter. Such a percolation picture of hadrons has been discussed in seminal works such as Refs. [26,27]. Also, the hadron-quark continuity [28,29] and hadron-quark crossover [30,31] have been discussed in relation to the existence of color superconductivity at high density. In this paper, we show that the crossover picture can lead to a stiffening of the EOS, unlike the case of the first-order transition, if the following conditions are met: (i) the crossover takes place at relatively low density (around three times the normal nuclear matter density), and (ii) the strongly interacting quark matter has a stiff EOS. This implies that the hadron-quark crossover provides us with a novel mechanism to support massive neutron stars with quark cores. A preliminary account of our results has been reported in Ref. [33]. We note that an interpolation between hadronic matter and stiff quark matter was previously considered phenomenologically in Ref. [32].
This paper is organized as follows. In Sect. 2, the characteristic features of the hadronic EOSs (H-EOSs) to be used at low densities are summarized. In Sect. 3, we treat the strongly interacting quark matter by using a Nambu-Jona-Lasinio (NJL)-type model and derive the quark EOS (Q-EOS) to be used at high densities. In Sect. 4, we describe our interpolation procedure to obtain the EOS in the hadron-quark crossover region. In Sect. 5 and 6, numerical results and discussions are given for the bulk properties of hybrid stars, such as the M-R relationship, the maximum mass M max , and the M-ρ c (central density) relationship. We discuss how these results depend on the different choice of H-EOS and Q-EOS. A comment on the cooling of NSs with respect to the hyperon mixture inside the core is also given. Section 7 is devoted to concluding remarks.

Hadronic EOS (H-EOS)
We consider several different EOSs with Y -mixing: • TNI2, TNI3, TNI2u, and TNI3u [34,35]: TNI2 and TNI3 are obtained by the G-matrix calculation with a Reid soft-core potential for N N and a Nijmegen type-D hard-core potential for Y N and Y Y . Also, a phenomenological three-body force [36] is introduced in the form of an effective N N force to reproduce the saturation point of symmetric nuclear matter with the incompressibility κ = 250 MeV (TNI2) and κ = 300 MeV (TNI3). For TNI2u and TNI3u, the three-body interaction is introduced universally in the form of effective N N, N Y , and Y Y forces. • AV18+TBF and Paris+TBF [37]: These are obtained by the G-matrix calculation but with a different choice of potentials; AV18 and Paris potentials for N N and a Nijmegen soft-core potential for Y N and Y Y . Also, a three-body force of Urbana type is introduced in the form of an effective N N force to meet the saturation condition. • SCL3 [40]: This is based on a relativistic mean field (RMF) model with chiral SU(3) symmetry and logarithmic potential motivated by the strong coupling lattice QCD approach. The phenomenological parameters of the model are determined to reproduce the saturation condition, bulk properties of normal nuclei, and separation energies of single-and double-hypernuclei.
In Fig. 2, we plot the pressure P for Y -mixed neutron star matter with β-equilibrium and charge neutrality as a function of baryon density ρ obtained from the EOSs listed above (Paris+TBF is not shown since it is almost the same as AV18+TBF). For comparison, P for neutron star matter without 3 Table A1 of Appendix A). The colors on each line are the same as those in Fig. 2 [34,35], Paris+TBF, AV18+TBF [37][38][39], and SCL3 [40]. κ is the nuclear incompressibility and ρ th is the threshold density of hyperon-mixing with ρ 0 (= 0.17/fm 3 ) being the normal nuclear density. R and ρ c denote the radius and central density for the maximum mass (M max ) NS, respectively. The numbers in parentheses are those without hyperons. * indicates that the numbers are taken from the figures in Ref. [37]. The use of the several kinds of EOS mentioned above, from different theoretical methods (G-matrix, RMF), with various stiffness values ranging from κ ∼ 190 MeV to 300 MeV and with the variation of ρ th (Y ) (2 − 4)ρ 0 , is expected to cover the present uncertainties of the H-EOSs. For completeness, numerical values of the pressure P and the energy density ε as a function of the baryon density are tabulated in Table A1 in Appendix A.

Quark EOS (Q-EOS)
The baryon density at the central core of the NSs would be at most 10ρ 0 . Although hadrons do not keep their identities in such a high density, the chemical potentials of the quarks are about (400-500) MeV, which is not high enough for the asymptotic freedom at work. Namely, the deconfined quarks inside the NSs, even if they exist, would be strongly interacting. An analogous situation at finite temperature is expected theoretically and has recently been confirmed by the relativistic heavyion collisions at RHIC and LHC; it is now called the strongly interacting quark-gluon plasma (sQGP).
Since lattice QCD to treat the strongly interacting quark matter (sQM) at finite baryon density is unfortunately not possible due to the notorious sign problem, we adopt an effective theory of QCD, the (2 + 1)-flavor Nambu-Jona-Lasinio (NJL) model. This model is particularly useful for taking into account important phenomena such as the partial restoration of chiral symmetry at high density [42][43][44][45].
The model Lagrangian we consider is where the quark field q i (i = u, d, s) has three colors and three flavors with the current quark mass m i . The term proportional to G S is a U (3) L × U (3) R symmetric four-fermi interaction where λ a are the Gell-Mann matrices with λ 0 = √ 2/3 I. The term proportional to G D is the Kobayashi-Maskawa-'t Hooft (KMT) six-fermi interaction, which breaks U (1) A symmetry. We consider two types of vector interaction (the second line of Eq. (1)): The term proportional to g V (> 0) gives a universal repulsion among different flavors, while the one proportional to G V (> 0) gives flavor-dependent repulsion.
In the mean-field approximation, the constituent quark masses M i (i = u, d, s) are generated dynamically through the NJL interactions (G S,D ): where σ i = q i q i is the quark condensate in each flavor, and (i, j, k) corresponds to the cyclic permutation of u, d, and s. The thermodynamic potential is related to the pressure as = −T log
where n i = q † i q i is the quark number density in each flavor, and S i is the quark propagator, which can be written as where iω = (2 + 1)π T and μ eff i is an effective chemical potential [46]. There are six independent parameters in the (2+1)-flavor NJL model: the UV cutoff, , the coupling constants, G S , G D , and g V (G V ), and the quark masses, m u,d and m s . Five parameters except for g V (G V ) have been determined from hadron phenomenology. We consider three parameter sets, summarized in Table 2: HK (Hatsuda and Kunihiro), RKH (Rehberg, Klevansky, and Hufner) and LKW (Lutz, Klimt, and Weise) [42][43][44][45].
The magnitude of g V (G V ) has not been determined well: Recent studies of the Polyakov-Nambu-Jona-Lasinio model applied to the QCD phase diagram suggest that g V may be comparable to or larger than G S [47,48], so that we change its magnitude in the following range: In Sects 5 and 6, we will show our results mainly for the HK parameter set with the vector interaction of the g V type. At the end of Sect. 5, we discuss how the results change in other cases. The Q-EOS with strangeness is obtained from the above model under two conditions: (i) charge neutrality among quarks and leptons, i.e. 2 3 n u − 1 3 n d − 1 3 n s − n e − n μ = 0, and (ii) the β-equilibrium among quarks and leptons, i.e. μ d = μ s = μ u + μ e and μ e = μ μ .
In Fig. 4, the number fractions (n u,d,s,e /n tot with n tot = n u + n d + n s = 3ρ) as a function of the baryon density ρ are plotted. Also, in Fig. 5, the constituent quark masses (M i ) as a function of ρ are plotted. The HK parameter set with the g V -type interaction is used in both figures. The flavorindependent g V -type interaction leads to a pressure in Eq. (3) depending only on μ eff i . Then, the number fractions and the quark masses as a function of ρ do not depend on g V .
At low baryon densities below a threshold density ρ th 4ρ 0 , the system is composed of only u, d, and e with n d ∼ 2n u due to charge neutrality and β-equilibrium (Fig. 4). In this region, the strong interaction among quarks (mainly the G S term in the NJL model) drives the partial restoration of 6   chiral symmetry and hence a rapid decrease of the constituent masses M u,d (Fig. 5). Due to the coupling between different flavors through the G D term, the strange quark mass M s in the Dirac sea is also affected slightly.
When the baryon density exceeds ρ th , the chemical potential of the strange quark μ s becomes larger than the strange quark mass (μ s > M s ), so that the system starts to have the strangeness degree of freedom. Since the strange quark is negatively charged, electrons start to disappear from the system 7  and the d quark fraction decreases at the same time (Fig. 4). We note that the system does not have the muon, because the electron chemical potential is smaller than m μ = 106 MeV at all densities. In the high-density limit, the system approaches the flavor-symmetric u, d, s matter without leptons. Once the s-quark appears in the system, M s is also suppressed, mainly due to the G S term (Fig. 5). The strangeness threshold ρ th does not depend on g V as already mentioned, but it does depend on the NJL parameter sets in Table 2: ρ th /ρ 0 = 4.0, 3.9, and 3.0 for HK, RKH, and LKW, respectively.
In Fig. 6, we plot the pressure (P(ρ) with a normalization P(0) = 0) of the strongly interacting quark matter for the HK parameter set with different values of the vector coupling (g V /G S = 0, 1.0, 1.5 according to Eq. (5)). Due to the universal repulsion of the g V -type vector interaction, the Q-EOS becomes stiffer as g V increases. As already mentioned, the onset density of the strangeness (shown by the filled circles) does not depend on g V . We note here that the present Q-EOS has a firstorder phase transition below 2ρ 0 for g V < 0.3G S . However, it does not affect the final results of the present paper, since such a low-density region is dominated by the hadronic EOS in our hadron-quark crossover approach, to be discussed in Sect. 4.

Hadron-quark crossover
As discussed in Sect. 1, treating the point-like hadron as an independent degree of freedom loses its validity as the baryon density approaches the percolation region. In other words, the system can be described neither by an extrapolation of the hadronic EOS from the low-density side nor by an extrapolation of the quark EOS from the high-density side. Under such a situation, it does not make much sense to apply the Gibbs criterion of two phases I and II, P I (T c , μ c ) = P II (T c , μ c ), since P I and P II are not reliable in the transition region.
Since a first-principles QCD calculation at high baryon density is not available and effective models at finite baryon density with proper treatment of the confinement phenomena do not exist at present, we will consider a phenomenological "interpolation" between the H-EOS and Q-EOS as a first 8 step. Such an interpolation is certainly not unique: Here we consider the two simplest possibilities, P-interpolation and ε-interpolation, as described below.
• P-interpolation as a function of baryon density: where P H and P Q are the pressure in the hadronic matter and that in the quark matter, respectively. An interpolating function f ± similar to ours has been previously considered at finite temperature in Refs. [49][50][51]. The windowρ − ρ ρ + characterizes the crossover region in which both hadrons and quarks are strongly interacting, so that neither the pure hadronic EOS nor the pure quark EOS are reliable. The percolation picture illustrated in Fig. 1 is best implemented by the interpolation in terms of the baryon density ρ instead of the baryon chemical potential. One should not confuse Eq. (7) with the pressure in the mixed phase associated with the first-order phase transition in which f ± is considered to the volume fraction of each phase. In our crossover picture, the system is always uniform and f − ( f + ) should be interpreted as the degree of reliability of H-EOS (Q-EOS) at a given baryon density.
To calculate the energy density ε as a function of ρ in a thermodynamically consistent way, we integrate the thermodynamical relation P = ρ 2 ∂(ε/ρ)/∂ρ and obtain with g(ρ) = 2 (e X + e −X ) −2 and X = (ρ −ρ)/ . Here ε H (ε Q ) is the energy density obtained from H-EOS (Q-EOS). ε is an extra term that guarantees thermodynamic consistency. Note that the energy per baryon from the extra term ε/ρ, which receives its main contribution from the crossover region, is finite even in the high-density limit. • ε-interpolation as a function of baryon density: Other thermodynamic quantities are obtained through the thermodynamic relation and μ = (ε + P)/ρ. Here P is an extra term that guarantees thermodynamic consistency; it is a localized function in the crossover region and obeys the property P(0) = P(∞) = 0.

Interpolated EOS
In the present section we consider the case of P-interpolation. The case of ε-interpolation will be discussed in Sect. 6. We note that the crossover window in both interpolations should satisfy the following physical conditions: (i) The system is always thermodynamically stable, d P/dρ > 0, and (ii) the normal nuclear matter is well described by the H-EOS so thatρ − 2 > ρ 0 is satisfied.     P(ρ) around ρ = 0 (ρ = ∞). Therefore, naive extrapolation of H-EOS and Q-EOS beyond their applicability would miss essential physics.
In Fig. 10, we plot the interpolated EOS using TNI2u and NJL for different values of g V in a wide range of baryon density. The filled circles denote the onset of strangeness degrees of freedom, either hyperons or strange quarks. 11

Mass-radius relation
We now solve the following Tolman-Oppenheimer-Volkov (TOV) equation to obtain the M-R relationship by using EOSs with and without the hadron-quark crossover: where we have assumed spherical symmetry, with r being the radial distance from the center of the star. In Fig. 11(a), we show the M-R relationship for various H-EOSs with hyperons whose onset is denoted by the filled circles. The crosses denote the points where maximum masses are realized: In all cases, M max does not reach 2M due to the softening of EOS by the hyperon mixture.
In Fig. 11(b), we show the M-R relationship with the EOS interpolated between H-EOS and Q-EOS: For the H-EOS, we consider the same EOSs as shown in Fig. 11(a), while, for the Q-EOS, we adopt the HK parameter set with g V = G S as a typical example. The crossover window is fixed to be (ρ, ) = (3ρ 0 , ρ 0 ). Cases for different parameters in Q-EOS as well as for different window parameters are discussed in the next subsection.
The red lines in Fig. 11(b) correspond to the cases with TNI2u and TNI2, the blue lines correspond to TNI3u and TNI3, and the green lines correspond to SCL3 and AV18+TBF. The onset of strangeness and the maximum mass are denoted by filled circles and the crosses, respectively. Irrespective of the H-EOSs, the interpolated EOS can sustain a hybrid star with M max > 2M : A smooth crossover around ρ ∼ 3ρ 0 and a stiff Q-EOS due to repulsive vector interaction are two fundamental reasons for this fact. Also, we note that the radius of the hybrid star with the interpolated EOS is in the range R = (11 ± 1) km for 0.5 < M/M < 2.0, except for the SCL3 case. 1 Such a narrow window of R independent of the values of M is consistent with the phenomenological constraints on R based on recent observations of both transiently accreting and bursting sources [52,53].
In Table 3, we show the maximum mass and the associated central density of a hybrid star with the interpolated EOS with g V = G S and g V = 1.5G S . In all combinations of H-EOS and Q-EOS, M max exceeds 2M with the central density, ρ c = (5.4-6.1)ρ 0 .
Let us now turn to the internal structure of the hybrid star, in particular its strangeness content. From the location of the filled circles in Fig. 11(b), one finds that the flavor-independent universal 1 The reason that the SCL3 case is different from others can easily be seen from Fig. 2: The pressure P of SCL3 is nearly twice as large as that of the other EOSs at ρ = (1 − 2)ρ 0 . This leads to a larger R for light NSs.

12/26
Downloaded from https://academic.oup.com/ptep/article-abstract/2013/7/073D01/1571314 by guest on 26 July 2018 three-baryon repulsion in TNI2u and TNI3u increases the onset density of the strangeness inside the hybrid star. This can be seen more explicitly by plotting the radial profile of the hybrid star: The upper panels of Fig. 12 show the ρ-r relationships for 2M and 1.44M hybrid stars with TNI2 (left) and TNI2u (right). The threshold densities of the strangeness given in Table 1   double lines. In our interpolated EOSs, the above stars turn out to have almost the same radius. The lower illustrations of Fig. 12 show the cross sections of the corresponding hybrid stars.
These figures imply that, even if the mass and the radius are the same, the strangeness content of the hybrid stars can be quite different. This point is of particular interest for the cooling problem of NSs. As is well known, NSs with a Y -mixed core undergo extremely rapid cooling due to the efficient ν-emission processes called "hyperon direct URCA" (Y -Durca, e.g., → p + e − +ν e , p + e − → + ν e ) and are cooled very rapidly below the thermal X-ray detection limit. Therefore, for NSs consisting of pure hadronic components with Y , only the very light-mass NSs (M < (1.0 − 1.2)M , as in Fig. 3) can escape from Y -Durca rapid cooling. This indicates the unlikely situation that all the NSs whose T s are observed should be light-mass stars, in spite of the fact that the observed mass distribution is centered around (1.4 − 1.5)M [2]. In contrast, in the case of the hybrid star with g V = G S (1.5G S ) under consideration, NSs as heavy as up to 1.9(2.0)M can avoid this rapid cooling, allowing the T s -observed NSs to be from the light-mass to heavy-mass stars (M ≤ (1.9 − 2.0)M , as in Fig. 13). 2 It is in order here to comment on the relationship between the maximum mass and the nuclear incompressibility κ. From the properties of finite nuclei, the nuclear incompressibility κ is estimated to be (240 ± 20) MeV [54]. The interpolated EOSs with TNI2 and TNI2u are consistent with this empirical κ, and yet they can reach M max > 2M . In other words, what is important to sustain massive hybrid stars is not the value of the incompressibility, but the stiffness of the EOS at and above ∼ 3ρ 0 .

Dependence on Q-EOS
To see how the hybrid star structure changes with the stiffness of Q-EOS, we plot the M-ρ c relationship for g V /G S = 0, 1.0, 1.5 with the HK parameter set in Fig. 13(a). We take TNI2u for H-EOS 15 Table 4. The values of M max /M (ρ c /ρ 0 ) for g V /G S = 1.0, 1.5, 2.0 with (ρ, ) = (3ρ 0 , ρ 0 ) and TNI2u. The parameter sets of the NJL model, HK, RKH, and LKW, are given in Table 2.  Table 5. M max /M (ρ c /ρ 0 ) for the HK parameter set with the flavor-dependent repulsion G S . The crossover window is (ρ, ) = (3ρ 0 , ρ 0 ) and the hadronic EOS is TNI2u.
1.87 (6.6) 1.99 (6.2) 2.07 (5.8) Table 6. M max /M (ρ c /ρ 0 ) under variation of the parametersρ and , which characterize the crossover window. H-EOS and Q-EOS are obtained from TNI2u and the HK parameter set, respectively. Columns without numbers are the excluded cases corresponding toρ − 2 < ρ 0 in Sect. 4. and the same crossover window as in Fig. 11. For comparison, the M-ρ c relationship with TNI2u only is plotted by the dashed line. Figure 13(b) shows the corresponding M-R relations. As anticipated, M max increases as g V increases. In Table 4, we show how M max and ρ c depend on the choice of g V and the choice of the NJL parameter set. Although the parameter dependence is not entirely negligible, a massive hybrid star is possible for sufficiently large values of g V . Finally, we consider the flavor-dependent vector interaction proportional to G V given in Eq. (1). In the high-density limit where u, d, s quarks have equal population, u † u = d † d = s † s , the g V interaction and the G V interaction make the same contribution to the pressure in the mean-field approximation if we make the following identification: G V = 3 2 g V . Motivated by this relation, we show M max and ρ c for G V /G S = 1.5, 2.25, 3.0 in Table 5. For the density relevant to the cores of hybrid stars, the flavor SU(3) limit is not yet achieved due to the s-quark mass (see Fig. 4). Therefore, the EOS for the flavor-dependent repulsion with G V = 3 2 g V is softer than the flavor-independent repulsion with g V . This can be seen by comparing the corresponding values in Table 5 and those in Table 4. In any case, a massive hybrid star is possible for sufficiently large values of G V .

Dependence on crossover window
In Table 6, we show M max and ρ c for different choices of the crossover window parameterized bȳ ρ and . TNI2u and the HK parameter set are adopted for H-EOS and Q-EOS, respectively. As the crossover window becomes lower and/or wider in baryon density, the interpolated EOS becomes stiffer and M max becomes larger. To be compatible with the observed massive NS with M = (1.97 ± 0.04)M , the crossover needs to occur in (2 − 4)ρ 0 . 16

Sound velocity of interpolated EOS
One of the measures to quantify the stiffness of the EOS is the sound velocity v S = √ d P/dε. In Fig. 14, we plot v S for our interpolated EOS with g V /G S = 0, 1.0, 1.5 as a function of ρ. The kinks of v S at ρ 4ρ 0 are caused by the softening of the EOS by the appearance of strangeness. The enhancement of v S of the interpolated EOS relative to the pure hadronic EOS takes place just at and above the crossover window.

Stability of the hybrid star
The neutron star is gravitationally stable if the average adiabatic index¯ satisfies the inequality [55] Here = d ln P/d ln ε is the adiabatic index. Also, λG M/R, with λ being a numerical constant of order unity, is a general relativistic correction whose magnitude is much less than 1. Since of our H-EOS is about 2 at all densities and of our Q-EOS is larger than 4/3 due to the constituent quark mass and the repulsive vector interaction, Eq. (14) is always satisfied and our hybrid star is gravitationally stable.

Neutron star properties with ε-interpolation
In this section we consider an alternative interpolation procedure using the energy density ε as a function of ρ given in Eq. (10). Shown in Fig. 15 is the energy density interpolated between TNI2u for H-EOS and NJL with g V = 0.5G S for Q-EOS. The crossover window is chosen to be (ρ, ) = (3ρ 0 , ρ 0 ) and is shown by the shaded area on the horizontal axis. The pressure obtained from the interpolated energy density   using the thermodynamic relation is shown in Fig. 16. Due to the extra positive term P in Eq. (12), the full pressure is larger than P H and P Q in the crossover region with the ε-interpolation procedure. Although P is necessary for thermodynamic consistency, its physical interpretation is not clear at the moment and is left for future studies. In Figs. 17 and 18, we show P as a function of ε and the sound velocity v S as a function of ρ, respectively. Because of the effect of P, the EOS becomes 18   stiff and v S is enhanced, particularly in the crossover region. Thus, the maximum mass of the neutron star would become large even for a moderate value of g V .
In Fig. 19(a), we plot the M-ρ c relationship between TNI2u for H-EOS and NJL Q-EOS for g V /G S = 0, 0.5 with the HK parameter set. We choose the crossover window as (ρ, ) = (3ρ 0 , ρ 0 ). For comparison, the M-ρ c relationship with TNI2u only is shown by the dashed line. Figure 13(b) shows the corresponding M-R relationships. As anticipated from Fig. 18, the maximum mass is larger than the P-interpolation case for a given g V .

19/26
In Table 7, we show M max and ρ c for different H-EOSs, vector-type interactions g V , and choices of crossover window parameterized byρ and . The ε-interpolation makes the EOS stiffen more drastically than the P-interpolation. Even for (g V ,ρ) = (0, 3ρ 0 ) and (g V ,ρ) = (0.5, 5ρ 0 ), the maximum mass M max can exceed 1.97M .

Summary and concluding remarks
Recent observation of a two-solar mass NS presents a challenging problem of how to reconcile the stiff EOS suggested from the observational side with the soft EOS due to hyperon-mixing from the theoretical side. In this paper we have studied this problem on the basis of the percolation picture from hadronic matter with hyperons to quark matter with strange quarks. We have constructed an EOS by interpolation between the H-EOS at lower densities and the Q-EOS at higher densities, and found that hybrid stars could have M max ∼ 2M , compatible with observation. This conclusion is in contrast to the conventional EOS for hybrid stars derived through the Gibbs construction, in which the resultant EOS always becomes softer than the hadronic EOS and thereby leads to a smaller M max . Our qualitative conclusion is insensitive to the choice of different types of H-EOS and different types of vector interaction in Q-EOS as long as (i) the crossover between the hadronic matter and the quark matter proceeds in a relatively low-density region, (ρ = (2 − 4)ρ 0 ), and (ii) the quark matter is strongly interacting and stiff (g V /G S ∼ 1). These conditions applied to the P-interpolation procedure can be relaxed further if ε-interpolation is adopted. We found that the the sound velocity v S , which increases rapidly in the crossover window for g V /G S ≥ 1, can nicely characterize the stiffening of the interpolated EOS and associated enhancement of M max .
The idea of rapid stiffening of the EOS starting from 2ρ 0 opens a possibility that the experimental nuclear incompressibility κ = (240 ± 20) MeV at ρ ∼ ρ 0 is compatible with the existence of massive neutron stars. Also, the idea may well be checked by independent laboratory experiments with medium-energy heavy-ion collisions.
Although the M-R relationship and M max are insensitive to the existence of universal three-body repulsion, the onset density of strangeness is rather sensitive to such repulsion. If we have three-body repulsion acting universally among baryons, most of the hybrid stars with M ≤ (1.9 − 2.0)M are free from the extremely efficient hyperon direct-Urca cooling process and can avoid contradicting observations.
Finally, we remark that the crossover region may contain richer non-perturbative phases such as color superconductivity, inhomogeneous structures, and so on [1]. How these structures, as well as the associated cooling processes, affect the results of the present paper would be an interesting future problem to be examined.

Appendix. EOS tables
In this Appendix, we show the concrete values of pressure P and energy density ε as a function of baryon density x ≡ ρ/ρ 0 for H-EOSs, Q-EOSs, and interpolated EOSs.
In Table A4, EOSs obtained by the interpolation of energy density between the TNI2u H-EOS and the NJL Q EOS with the HK parameter set with g V /G S = 0, 1.0, 1.5 for (ρ, ) = (3ρ 0 , ρ 0 ) are listed. Table A2. Pressure P (MeV/fm 3 ) and energy density ε (MeV/fm 3 ) as a function of baryon density x ≡ ρ/ρ 0 for NJL Q-EOSs with the HK parameter set for g V = G S [44].  by using the least-squares method.

TNI2
TNI2u  Table A4. Pressure P (MeV/fm 3 ) and energy density ε (MeV/fm 3 ) as a function of baryon density x ≡ ρ/ρ 0 for EOSs obtained by the interpolation of energy density between the TNI2u H-EOS and the NJL Q EOS with the HK parameter set with g V /G S = 0, 0.5 for the (ρ, ) = (3ρ 0 , ρ 0 ) case.