Implications of NICER for neutron star matter: the QHC21 equation of state

The recent NICER measurement of the radius of the neutron star PSR J0740+6620, and the inferred small variation of radii from 1.4$M_\odot$ to 2.1$M_\odot$, reveal key features of the equation of state of neutron star matter. The pressure rises rapidly in the regime of baryon density $n \sim$ 2-4 times nuclear saturation density, $n_0$ -- the region where we expect hadronic matter to be undergoing transformation into quark matter -- and the pressure in the nuclear regime is greater than predicted by microscopic many-body variational calculations of nuclear matter. To incorporate these insights into the microscopic physics from the nuclear to the quark matter regimes, we construct an equation of state, QHC21, within the framework of quark-hadron crossover (QHC). We include nuclear matter results primarily based on the state-of-the-art chiral effective field theory, but also note results of using nuclear matter variational calculations based on empirical nuclear forces. We employ explicit nuclear degrees of freedom only up to $n \sim 1.5n_0$, in order to explore the possibility of further physical degrees of freedom than nucleonic here. The resulting QHC21, which has a peak in sound velocity in $\sim 2$-$4 n_0$, is stiffer than the earlier QHC19 below 2$n_0$, predicting larger radii in substantial agreement with the NICER data.


INTRODUCTION
Recent data from the NICER observatory has yielded the important inference that the equatorial radius of the most massive neutron star known, PSR J0740+6620, with mass M/M = 2.08 ± 0.07 (Fonseca et al. 2021), is 13.7 +2.6 −1.5 km, as analyzed by Miller et al. (2021) combining NICER and XMM data for PSR J0740+6620. With nuclear physics constraints at low density, and gravitational radiation data from GW170817 added in, Miller et al.'s inferred radius drops to 12.35 ± 0.75 km for a 2.08 M neutron star. The analysis of Riley et al. (2021) using the NICER and XMM data for PSR J0740+6620 gives 12.39 +1.30 −0.98 km. These analyses lead to an updated radius of 1.4M neutron stars, 1 R 1.4 = 12.45 ± 0.65 km (Miller et al. 2021) and 12.33 +0.76 −0.81 km (Riley et al. 2021) to 12.18 +0.56 −0.79 km (Raaijmakers et al. 2021). Our aim in this paper is to delineate what the NICER measurements reveal about the microscopic properties of dense matter at baryon densities beyond nuclear saturation density n 0 , and to translate this understanding into a physically motivated new equation of state of neu-1 Here R M/M denotes the radius of a star of mass M . tron star matter, QHC21. A key observation is that the effective lack of change of the radius as the mass increases from 1.4 to 2.08 M strongly points to a very rapid growth of pressure with density in the range ∼ 2-4 n 0 , a result confirmed in the equations of state of Miller et al. (2021) used to fit the data. Such behavior is in contrast to the gradual decrease of radii of neutron stars of masses well below 1 M as the mass increases, a consequence of the pressure in nuclear matter not strongly resisting the increasing gravity. The lack of significant variation of the radii from 1.4M to ∼ 2.1M rules out substantial softening of matter (as would result from a strong first order phase transition) between 2-3n 0 and 4-5n 0 (Drischler et al. 2021a); the radical shrinkage in the radius such softening would induce is incompatible with the new NICER data.
A second key observation is that the radius of a 1.4M neutron star inferred from NICER measurements, for which the central baryon density is ∼ 2-3n 0 , indicates that matter in the nuclear regime has greater pressure than predicted earlier by microscopic many-body variational calculations of nuclear matter, as in the Akmal-Pandharipande-Ravenhall (APR) and Togashi 2 equa-tions of state (Akmal et al. 1998;Togashi et al. 2017). The radius based on the NICER analyses is ∼ 0.7-0.8 km larger than predicted using earlier equations of state based in the nuclear regime on variational calculationsfor example, in QHC19 , R 1.4 11.6 km. Since the radius of a 1.4M neutron star is essentially determined by the equation of state up to densities ∼ 2n 0 , a larger radius requires greater stiffness here. Similarly, 2.1M neutron stars were predicted in QHC19 to have radii 11.3-11.5 km. Indeed, the NICER data provide the first observational test of equation of states for neutron star matter obtained microscopically (Akmal et al. 1998;Togashi et al. 2017;Lonardoni et al. 2020;Drischler et al. 2021a).
The rapid pressure rise in the range ∼ 2-4 n 0 indicates the importance of quark degrees of freedom in this regime. Such a rapid rise does not naturally occur in descriptions of dense matter treating nucleons as elementary particles; rather, in models with realistic two-and three-nuclear forces (Akmal et al. 1998;Togashi et al. 2017), the sound velocity gently rises without stop, eventually exceeding the causal limit at densities 5n 0 . The gentle growth in sound velocity can be also seen in relativistic mean field models of hadrons (e.g., Steiner et al. (2013); Typel et al. (2010)).
Because nucleons are composed of quarks, increasing density causes the quark states at low momenta to become filled and triggers a rapid rise in pressure. This process can be understood, initially, in terms of nucleons being pushed into the relativistic regime (McLerran & Reddy 2019), or more generally in terms of higher momentum quark states forming a Fermi sea, as Kojo (2021a) and Kojo & Suenaga (2021) describe. Such duality between relativistic nucleons and quark descriptions is a useful first picture of the domain intermediate between pure nucleonic and quark matter .
The NICER analysis of Miller et al. (2021) assumes the Togashi equation of state (Togashi et al. 2017) for the crust up to baryon densities n 0.5n 0 , and parameterized equations of state -a polytrope model, a model with a parameterized speed of sound, and a Gaussian model -beyond that density. On the other hand, the analyses of Riley et al. (2021); Raaijmakers et al. (2021) use the equation of state from chiral effective field theory (χEFT, see below) and polytrope or constant speed of sound models. 3 Here we construct a new equation of state of neutron star matter, QHC21, to incorporate the messages of NICER into a microscopic, rather than phenomenological, approach which links the neutron star equations of state to the underlying physics in QCD. Microscopic nuclear forces and empirical three-body nuclear forces, computed using variational techniques for nuclear many-body problems, as in APR the equation of state; the same nuclear forces are used consistently in the nuclear liquid and crust equations of state.
3 Recent analyses similar to NICER but with nuclear constraints from χEFT (Pang et al. 2021) yield somewhat lower estimates for the radii R 1.4 = 11.94 +0.76 −0.87 km and R 2.08 = 11.96 +0.80 −0.75 km than the NICER analyses. Previously the same group (Dietrich et al. 2020) gave R 1.4 = 11.75 +0.86 −0.81 km without the data of PSR J0740+6620; the new NICER data leads to an increase in the estimate by 0.2 km, illustrating how data for 2M neutron stars constrains the properties of neutron stars of lower mass and the equation of state in the core. descriptions, once established, enable one to explore correlations among various microscopic quantities and matter composition, and reveal important observable consequences in, e.g., neutron star cooling, and neutrino emission from supernovae and neutron star mergers.
We use the chiral effective field theory (χEFT) nuclear equation of state (Lonardoni et al. 2020;Drischler et al. 2021a,b) at densities below ∼ 1.5n 0 . This approach, based essentially on an expansion in the effective momentum transfers, can provide control of the uncertainty. As discussed by Drischler et al. (2021a) (see their Fig. 1), the uncertainty in the pressure within χEFT is small at n 0 , and grows to ∼ 25% at ∼ 1.4n 0 ; the uncertainties at 2n 0 are too large to continue to use this approach at such density.
We focus on nuclear matter calculations based on microscopic forces. While phenomenological equations of state can also reproduce nuclear matter properties at densities ∼ n 0 (e.g., Serot & Walecka (1986)), the microscopic approach has the advantages that calculations can be improved both by systematically including higher order corrections and by refining the low energy constants in the chiral Lagrangian, and critically that they enable one to differentiate uncertainties in calculations of microscopic processes from uncertainties in the underlying physics, crucial in determining the importance of new beyond-nuclear physics towards high density. On the other hand, mean-field Lagrangian approaches are flexible in adding new degrees of freedom and in fitting new data, but they have difficulties in finding good rationales to truncate large powers of field variables (Weinberg 1995), which can significantly impact predicted properties of matter. Figure 1 shows the range of the pressure in nuclear equations of state vs. baryon density in χEFT approaches (labeled Lynn, Hebeler, Tews and Drischler) and variational approaches (labeled APR and Togashi). The variational approaches are significantly softer, as consequence we expect of approximations made in constructing the variational solutions. To allow for stiffer matter, we focus on the nuclear equations of state constructed by Drischler et al. (2021a) based on beyond-leading-order (N 3 LO) χEFT. Specifically we take the central value of the N 3 LO χEFT pressure at given density. This pressure, larger than in the Togashi equation of state, by ∼ 30-40% around n 0 , leads to larger radii. We use the χEFT nuclear equation of state at densities below ∼ 1.5n 0 . We also indicate results using the Togashi equation of state, thus covering the uncertainty range of microscopic nuclear calculations at lower densities.
As discussed above, the rapid rise in pressure indicated by the NICER data well motivates the consideration of the onset of higher momentum degrees of freedom around 2-3 n 0 , and a breakdown in an approach based on interactions of nucleonic degrees of freedom. Certainly at higher densities the degrees of freedom become those of strongly interacting quarks. To facilitate the transition to quark degrees of freedom, we allow a transition region between nucleonic and quark degrees of freedom within the Quark-Hadron Crossover (QHC) framework (Baym et al. 2018. We explore a greater range of equations of state above nuclear matter density that are more readily accessible by non-nucleonic descriptions, accom-  Raaijmakers et al. (2021)) and variational methods based on phenomenological nuclear forces (Akmal et al. 1998;Togashi et al. 2017). The "central" χEFT curve corresponds to the average of the data in Drischler et al. (2021a).
modating an earlier breakdown of the nuclear description, and permitting possibly stiffer interpolated equations of state, by using explicit nuclear degrees of freedom only up to n ∼ 1.5n 0 .
We use quark matter equations of state, calculated as calculated within the Nambu-Jona-Lasinio (NJL) model, now starting at a lower density, n 3.5n 0 , than ∼ 5n 0 where baryons of vacuum radius 0.8 fm begin to overlap. 4 This modification allows for the possibility that quark exchanges in baryon interactions drive earlier overlap of baryons and hence earlier onset of quark matter. In the transition region between nuclear and quark degrees of freedom we carry out a smooth and highly constrained interpolation. 5 At high densities, the NJL model provides the best available equation of state taking into account the restoration of chiral symmetry as well as color superconductivity in strongly interacting quark matter. As we see in detail below, the new equation of state well accommodates the NICER data, and indicates the extent to which the central density of PSR J0740+6620 is high enough to contain quark matter.

PHYSICS OF QHC21
We detail here the construction of QHC21. We specify the thermodynamics in terms of the pressure P as a function of baryon chemical potential µ B . Nuclear equations of state are used up to a density n L , and 4 QHC19 previously used the Togashi nuclear equation of state at n ≤ 2n 0 , and quark matter equations of state calculated within the NJL model for n 5n 0 . 5 Although theoretical calculations in the crossover region are difficult, possibly involving baryon-like and quark-like degrees of freedom at the same time (McLerran & Reddy 2019;Fukushima et al. 2020;Kojo 2021a), the physics in this regime is highly constrained by having to match the nuclear and quark matter equations of state at the boundaries of the crossover region. The constraints that the sound velocity not exceed the speed of light, and that the matter in the interpolation region be thermodynamically stable further constrain possible interpolations. The interpolating curves reveal the general trends (Masuda et al. 2013a(Masuda et al. ,b, 2016Kojo et al. 2015;Annala et al. 2020). quark equations of state above a density n U (> n L ). In the crossover between these regimes, we smoothly interpolate the pressure, the baryon density n = ∂P/∂µ B , and the susceptibility χ B = ∂ 2 P/∂µ 2 B at the matching boundaries. Thermodynamic stability imposes the requirement that the interpolated pressure leads to nonnegative χ B , and causality that the adiabatic speed of sound c s = ∂P/∂ε does not exceed the speed of light, c. Here ε = −P + µ B n is the energy density.
In QHC21 we consider primarily the χEFT equation of state in the nucleonic regime, but also note results using the Togashi equation of state there (labeled QHC21T). To describe quark matter we employ the Nambu-Jona-Lasinio (NJL) model. Although the NJL model does not include confining effects it is a good framework in which to capture dynamical effects, e.g., chiral symmetry breaking, concisely (Hatsuda & Kunihiro 1994), and in addition it reproduces quantities that are not very sensitive to confinement such as effective quark masses and low energy constants in the chiral effective Lagrangian for pions. QHC thus incorporates general features of quark matter by exploring a wide range of the NJL model parameters consistent with neutron star observables. The coupling constant G and the ultraviolet cutoff Λ of the NJL model, specified in Baym et al. (2018), are the relevant parameters describing chiral symmetry breaking. In addition, we include the short range density-density repulsion in QCD (Kunihiro 1991;Song et al. 2019), of strength measured by a coupling constant g V , and an attraction between quarks, proportional to a coupling H, which governs the color superconductivity pairing gap. 6 The quark equations of state, parametrized by given (g V , H), and the nuclear equation of state are first interpolated in a given range (n L , n U ) by using fifth order interpolating polynomials in µ B , for which the coefficients are uniquely fixed by six boundary conditions. We then determine the ranges of (g V , H) that are consistent with stability, causality, and the existence of the known high mass neutron stars. This procedure predicts M vs. R, thermodynamic quantities such as baryon density as a function of µ B , as well as more microscopic quantities such as effective quark masses and diquark paring gaps.

NEUTRON STAR STRUCTURE
Figures 2-4 display, for (n L , n U )/n 0 = (1.5, 3.5), the results of QHC21 in the (g V , H) plane for the maximum mass and the corresponding central density, as well as the radius and central density of 2.08 M and 1.4 M stars. We first determine the domain of (g V , H) that leads to equations of state satisfying the causality condition and that are consistent with the lower bound for the maximum mass, M max /M ≥ 2.08 ± 0.07; these allowed regions are shown as colored in Figs. 2-4. As a guide we choose several characteristic values in the (g V , H) plane, focussing on (g V , H)/G = (1.0, 1.50) (set A χ ), (1.1, 1.52) (set B χ ), (1.2, 1.54) (set C χ ), and (1.3, 1.56) (set D χ ). We denote these parameter sets with a subscript χ to distin- The interpolation ranges (n L , n U ) in units of n 0 , and the interaction parameter sets (g V , H) for QHC21 Aχ-Dχ, and QHC21T A T -D T . Tabulated for each parameter set are the maximum mass, Mmax, in units of M , radii R, in km, and central baryon densities, n c , in units of n 0 , of neutron stars of masses 1.4, 2.08 M , and Mmax. 12.0 12.2 12.4 12.5 11.5 11.8 11.9 12.0 R Mmax 11.7 11.5 11.4 11.2 10.9 11.1 11.1 11.3 guish them from the sets A, B, C, and D in QHC19. In addition Table 1 summarizes the parameters and results of both QHC21 and QHC21T (with parameter sets denoted by a subscript T). In Appendix B we vary the range of interpolation, and find that the domains of (g V , H) can vary by 10-20%, with no qualitative changes. The upper panel of Fig. 2 shows the maximum neutron star mass M max , in the color scale on the right, as a function of (g V , H) measured in units of the NJL coupling G. Smaller repulsion g V leads to smaller maximum masses; the requirement M max ≥ 2.08M is satisfied only for g V /G 0.84, while g V /G 0.74 leads to M max /M 2.01. The largest allowed 7 M max /M is 2.43 with (g V , H) (1.42, 1.59). The lower panel of Fig. 2 shows the central density, n c , at M max . A larger M max allows a smaller n c , consistent with the general rule that the stiffer the equation of state the larger the maximum mass and the smaller the central density at the maximum. The restriction M max /M ≥ 2.08 imposes the condition n c /n 0 6.4. At M max /M 2.40, n c /n 0 5.3. Figure 3 shows the dependence of the radius and the central density of 2.08M neutron stars on g V and H. Within the allowed band of (g V , H), the radii range from 10.8 km to 12.9 km. Radii 11 km, significantly smaller than R 1.4 12 km, appear for equations of state whose maximum mass barely reaches 2.08M , since the radius in general shrinks radically near the maximum mass.
Apropos of whether the central density of PSR J0740+6620 with a mass 2.08M is high enough to accommodate quark matter as characterized by overlaps of baryons, we see in the lower panel of Fig. 3 that the central density is 5n 0 , slightly smaller than the den-  6.8 6.8 sity where baryons of radius ∼ 0.8 fm begin to overlap, and is 3n 0 , comparable to n U = 3.5n 0 as adopted here. What is certain is that the core density is considerably larger than the density where pure hadronic calculations are applicable. The equations of state at a density slightly below n U are strongly correlated with quark matter models at densities n U .
The central density, and thus the presence of quark matter, is sensitive to how close the mass is to the maximum, M max , allowed by the equation of state. This is because the radius shrinks rapidly by ∼ 0.5-1 km in a small mass interval ∼ 0.05M below M max as seen from Fig. 5. The radius R 2.08 being comparable to R 1.4 indicates that 2.08M is still close to the lower edge of the above interval. Hence, we expect that M max is larger than 2.08M by at least 0.05M . For (g V , H) satisfying this condition, n c | 2.08 5n 0 .
The range of radii of 1.4M neutron stars, 11.8-12.7 km, is consistent with the NICER result R 1.4 = 12.4±0.6 km. The corresponding central density is 2.1-3.2n 0 . These radii are ∼ 0.3-1.2 km larger than in QHC19 (R 1.4 11.5 km, with (n L , n U )/n 0 = (2.0, 5.0)), a consequence primarily of using the stiffer χEFT nuclear equation of state. The Togashi equation of state, even if one allows more freedom for the interpolated pressure in the   density range 1.5-2.0 n 0 , leads to radii 11.6-12.2 km, on the smaller side of the NICER band. Subsequent NICER determinations of the radii of neutron stars with intermediate masses would provide a further observational test of the inferred rapid pressure rise in the equation of state. The estimate of R 1.4 can be further improved by additional gravitational wave data from neutron star mergers. The tidal deformability extracted from the GW170817 event tends to favor a smaller R 1.4 than extracted from the NICER data alone, and this trend affects the NICER inference of R 1.4 = 12.4 ± 0.6 km used in the present paper. Provided the χEFT results stay largely unchanged in the future, a larger R 1.4 ( 12.4 km) would require even more rapid stiffening than in QHC21, while a smaller one would shift the onset of rapid stiffening to higher densities.
Having described the overall trend of the masses, radii, and central densities in the parameter space, we look more closely at results for the parameter sets A χ -D χ for QHC21. The M -R relation for these parameters is shown in the upper panel of Fig. 5; the lower panel shows the corresponding results for QHC21T for the parameter sets A T -D T . The numbers next to the heavy dots indicate the central density of the star, in units of n 0 , for the given mass. As we have matched interpolating  functions with nuclear equations of state up to second derivatives at n = 1.5n 0 , the deviation of the M -R curves from what the pure, i.e., without quark matter, χEFT or Togashi equations of state would give appears only after n increases considerably beyond the matching point n = 1.5n 0 . (This trend can also be seen in the speed of sound in Fig. 7.) The central χEFT equation of state in Fig. 1, which is available only for n 2.2n 0 , predicts overall larger radii for low core densities, leading to better agreement of QHC21 with the NICER data. As we see in Fig. 6, in QHC21 the equations of state become stiffer than QHC21T and pure Togashi at ∼ 2n 0 , producing 1.4M neutron stars at a lower central density than with Togashi alone.
An important feature of the QHC21 equation of state is the behavior of the sound velocity. 8 As seen in Fig. 7, c 2 s exceeds the conformal limit c 2 s = c 2 /3 at n 2.1-2.3n 0 , a considerably lower density than 2.8n 0 with Togashi alone, and also less than the ∼ 2.6n 0 in QHC19. The peaks in c 2 s for the QHC21 set A χ -D χ or A T -D T are at 8 We note that the density dependence of the NJL couplings can change c 2 s in the quark matter region. The couplings should decrease at very large density so that the equations of state are dominated by the quark kinetic energy and c 2 s approaches 1/3, as it should in the high density limit.
2.1-3.0 n 0 , a considerably lower density than in QHC19 they are in the range 3.0-4.0n 0 .
Purely nucleonic models achieve stiff equations of state by having repulsive three-or more-body forces at short distance. Including only few-body forces, such models constrained by low density experiments lead to slow growth in c 2 s , as seen in the Togashi curve in Fig. 7. Quicker growth requires more-body forces whose dominance would invalidate using nucleons as effective degrees of freedom. The rapid growth of c 2 s found in QHC21 is unlikely to be achieved by nucleonic models, differentiating QHC from nuclear models.
The peak in the sound velocity is a novel feature of the crossover at finite baryon density, as originally pointed out phenomenologically by Masuda et al. (2013b) and recently analyzed theoretically by McLerran & Reddy (2019); Kojo (2021a); Kojo & Suenaga (2021) in terms of the quark substructure of baryons and quark Pauli blocking. In the crossover domain, quark Pauli blocking puts kinematic constraints on baryons pushing them to be relativistic. In quark descriptions, relativistic quarks start to contribute directly to the thermodynamic pressure even before the quark Fermi sea is established.
The theoretical estimates Kojo (2021a); Kojo & Suenaga (2021) suggest that quark Pauli blocking effects set in around ∼ 1-3n 0 , radically increasing the nonrelativistic pressure ∼ n 5/3 /m N (where m N is the nucleon mass) to the relativistic behavior ∼ N c n 4/3 , where N c (= 3) counts the number of quarks in a baryon. The corresponding changes in the energy density are modest; ∼ m N n in non-relativistic regime and ∼ N c n 4/3 in relativistic regime are comparable in the crossover domain. Such a rapid growth in pressure with modest change in energy density leads to rapid enhancement of c 2 s . As this transient regime comes to a close, both the pressure and energy density scale as in relativistic regime and hence c 2 s eventually relaxes to 1/3, leaving a sound velocity peak in the crossover domain.
A sound velocity peak of relativistic origin does not exist in a hot QCD plasma (Bazavov et al. 2014), where the speed of sound develops a dip instead of a peak in the crossover region. The finite temperature transition from a hadron resonance gas to a quark gluon plasma is largely driven by non-relativistic resonances which are energetically disfavored but are important due to their large entropies. Such entropic effects are absent in cold dense matter.
Finally, we compare the pressure in QHC21 with constraints deduced by the NICER analyses. Figure 8 shows the pressure in QHC21 as a function of the baryon mass density, ρ = m N n, where m N = 1.67 × 10 −24 g is the nucleon mass. Then Fig. 9 compares the pressure as a function of baryon chemical potential µ B , with the constraint "All Measurements (Gaussian Process)" in Miller et al. (2021) at a 95%CL. 9 The QHC21 equations of state are quite consistent with the NICER inferences.

DISCUSSION
Equations of state that exhibit a peak in the sound velocity occurring at the relevant intermediate densities robustly satisfy constraints from the NICER data. In the quark regime, matter can be stiffened by several physical effects (as discussed in greater detail in Appendix A). First, the vector repulsion between quarks, which allows quark matter to support neutron stars above 2M , stiffens the matter in general, and produces sound velocities above the conformal limit, c/ √ 3. Furthermore, quark pairing correlations leading to color superconductivity (Alford et al. 2005)  Within QHC21, the central density of PSR J0740+6620, 5n 0 , is well above densities where pure hadronic calculations are valid, and is entering the transition regime to strongly interacting quark matter. The cores of higher mass neutron stars, possibly to be discovered in the future, could reach the density region well beyond the peak in the sound velocity. To describe such stars quantitatively, fully microscopic calculations of matter undergoing the transition from nucleonic to quark degrees of freedom will be needed. Nonetheless the properties of quark matter strongly influence physics in the transition, possibly even slightly above saturation density.
Tables of the equations of state QHC21 for parameter sets A χ -D χ and A T -D T are available on the Com-pOSE archive at https://compose.obspm.fr/eos/236 through https://compose.obspm.fr/eos/243. In QCD at ultrahigh densities in the absence of significant length scales, the velocity approaches the conformal limit, c s = c/ √ 3, from below; this limiting value is the conformal bound. We focus in this Appendix on two questions: The first is how the equation of state determines whether the sound velocity can exceed the conformal bound. 12 In considering quark matter in neutron stars we are particularly interested in the domain where quark descriptions are natural but the matter is not dense enough for perturbative treatments (valid at 40n 0 ). While such a domain, above 5n 0 , may not show up in the cores of neutron stars, understanding c s there should give important constraints on the behavior of equations of state in the range of baryon density ∼ 2-5n 0 which is only weakly constrained by perturbative results alone. The second question we discuss is the behavior of the sound velocity in QED, a theory that is not asymptotic free. As we show, in the perturbative treatment of the massless electron gas, the sound velocity exceeds c/ √ 3. We set c = 1 here for simplicity. The sound velocity involves a microscopic interplay between the kinetic energy density, ε kin and the interaction energy, ε int . We first show schematically how the conformal bound can be violated in relativistic matter. The quark number density, n Q , is given in terms of the quark Fermi momentum, p F , for equal mass quarks, by n Q = N f p 3 F /π 2 , with N f the number of quark flavors present; the kinetic energy density for massless quarks, is thus ∼ p 4 F ∼ n 4/3 Q . Calculationally it is simplest to write the total energy density ε(n Q ) in terms of p F as where a = 3N f /4π 2 and h(p F ) = ε int /ε kin is the interaction energy relative to the kinetic energy. Then the pressure is where primes denote derivatives with respect to p F . Equations. (A1) and (A2) imply so that In the absence of h and h , the h term does not affect the sound velocity. Let us first consider h having a power-law (PL) behavior, h = bp ζ F , where b is a constant. Then, c 2 s (PL) = 1 3 1 + (4 + ζ)ζh 4 + (4 + ζ)h .
We see that a repulsive interaction, h > 0, with ζ > 0 increases the sound velocity above the conformal bound. An attractive interaction, h < 0, can also stiffen the equation of state if ζ < 0. Such density dependence is typical in correlation effects near the Fermi surface. Pairing correlations in color superconductivity (Alford et al. 2005) are in a regime with ζ ∼ 2 to within the density dependence of the gap (which is not known well in the non-perturbative regime), and thus tend to increase the sound velocity. Baryonic correlations in quarkyonic matter (McLerran & Reddy 2019) produce similar effects. For example, the repulsive vector interaction between quarks in the NJL model, ∼ g V n 2 Q , has ζ = 2, and in the three-flavor NJL model employed in the text, h NJL = (4/π 2 )g V p 2 F . Then which always exceeds 1/3 and approaches to 1 in the high density limit.
ening at low density is consistent with observational constraints (Tews et al. 2018). 11 Underlying diquark pairing is the color magnetic interaction in QCD, which plays a crucial role as well in hadron spectroscopy (De Rujula et al. 1975) and short distance baryon-baryon interactions (Oka & Yazaki 1980;Park et al. 2020).
12 Tews et al. (2018), also Greif et al. (2019), assuming the validity of nuclear calculations up to n ∼ 1.5n 0 where c 2 s is only ∼ 0.1 -0.2, point out that to have sufficient stiffening at higher density in order to reach neutron stars over two solar masses, it is necessary that c 2 s exceed 1/3 there. On the other hand, Annala et al. (2020) construct stiff equations of state by allowing c 2 s to reach ∼ 1/3 from below at n = 1.1-1.5n 0 , a velocity well above the nuclear results; such early stiffening results in stiff equations of state that do not violate the conformal limit. Such radical stiffening in the interval 1.1-1.5 n 0 could be tested in nuclear experiments.
On the other hand, in the perturbative QCD (pQCD) regime at ultrahigh densities, beyond the range one actually encounters in neutron stars, ζ = 0, and an additional physical effect comes into play, the strong coupling constant α s in dense matter runs with p F (or equivalently the chemical potential) as ∼ 1/ ln(p F /Λ QCD ) with Λ QCD the QCD scale parameter ( 340 MeV for N f = 3). To first order in α s , this dependence leads explicitly to while h is of order α 3 s . From Eq. (A4) we then find for three flavors, approaching the asymptotic value 1/3 from below. We see here that the asymptotic approach to the conformal limit is the result of a competition between the gluonic factor 33, which produces an approach from below, and the fermionic term 2N f , which produces an approach from above. In QCD the gluonic term wins, and the approach is from below. Equation (A6), valid at intermediate densities relevant to neutron star interiors, and Eq. (A9) for asymptotically high density, are smoothly joined by a density dependent vector interaction, g * V , as introduced in Song et al. (2019), Here m g 400 MeV is a non-perturbative scale proportional to Λ QCD , and the running coupling for p F < m g is frozen at α s (m g ) with g V ≡ 4πα s (m g )/27m 2 g . With this density dependent vector coupling, h = (4/π 2 )g * V p 2 F , which approaches h NJL at p F m g and h pQCD at p F m g . The Fermi momentum at which c * s 2 − 1/3 changes sign from positive to negative is 1.5Λ QCD , corresponding to a baryon density n ∼ 10n 0 .
Quantum electrodynamics, in contrast to QCD, has no "charged" gluons, so that it is not asymptotically free and the sound velocity exceeds 1/3. To see this in detail, we note that the energy density in the relativistic limit to first order in α e is so that h = α e /2π. The evolution of α e in QED is equivalent to dropping the 33, letting N f = 1 and Λ QCD → m e in α s , and multiplying by a factor 2 which arises from the algebraic difference of the interactions in QED and QCD. Then, to lowest order, we have p F ∂α e /∂p F = α 2 e /6π. The resulting sound speed, to order α 2 e , is given by which exceeds 1/3 (and grows as the Fermi momentum increases). In the low density non-relativistic limit, the sound speed of the electron gas is c 2 s (QED) = (p F /m e ) 2 /3, which is smaller than 1/3.

B. EXPLORING OTHER INTERPOLATION RANGES
While we take the interpolation range n L = 1.5n 0 and n U = 3.5n 0 for our main results in the text, we examine in this Appendix the dependence of the NJL parameters and the resulting mass-radius relations on the matching densities n L and n U . We consider n L = (1.3, 1.7)n 0 and n U = (3.5, 4.5)n 0 in the four possible combinations, continuing to use the χEFT equation of state in the nuclear regime. Figure 10 shows the allowed ranges of (g V , H), with the parameter sets A χ -D χ in the main text marked as guides. As we can see, changing the interpolation range modifies the domains of (g V , H) by only 10-20%.
Increasing the size of interpolating interval should in principle monotonically extend the window for (g V , H), as well as allow a larger mass range. For example, the results of (n L , n U ) = (1.3, 4.5)n 0 should contain the results of (1.7, 3.5)n 0 as special cases. However in our modeling where we use simple polynomials which cannot cover all allowable curves, this is not the case. We use simple structureless polynomials so as not to introduce exotic structures by hand, and yet the interpolation results in a highly nontrivial structure in the equation of state, e.g., the peak in sound velocity. Owing to our not exploring all possible interpolation curves, the maximum mass we find in a given interpolation window should be understood as a lower bound of the maximum mass. Shown in Fig. 11 are c 2 s vs. n (upper panels), M vs. R (middle panels), and M vs. n (lower panels), for the matching densities n L = 1.5n 0 and n U = (3.5, 4.0, 4.5)n 0 . For given g V = 1.1, and 1.2G, the lower (bold lines) and upper (thin lines) bounds on H are as shown in Fig. 10  There are clear correlations between the locations of the peaks in the sound velocity, the radii of 1.4M neutron stars, and the growth of the neutron star mass as a function of the central density. A peak in sound velocity at lower density indicates rapid stiffening which results in a larger radius. In the upper panels of Fig. 11, the locations of the peaks change by ∼ 1-1.5n 0 as H is varied at given g V , and the radii change by ∼ 0.3-0.7 km, as shown in the middle panels. The central value of H leads to R 1.4 12.2-12.4 km, and with the variation of H, the M -R relations remain largely consistent with the NICER data.