Evidence for quark-matter cores in massive neutron stars

The theory governing the strong nuclear force, Quantum Chromodynamics, predicts that at sufficiently high energy densities hadronic nuclear matter undergoes a deconfinement transition to a new phase of quarks and gluons. Although this has been observed in ultrarelativistic heavy-ion collisions, it is currently an open question whether quark matter exists inside neutron stars. By combining astrophysical observations and theoretical ab-initio calculations in a model-independent way, we find that the inferred properties of matter in the cores of neutron stars with mass corresponding to 1.4 solar masses are compatible with nuclear model calculations. However, the matter in the interior of maximally massive, stable neutron stars exhibits characteristics of the deconfined phase, which we interpret as evidence for the presence of quark-matter cores. For the heaviest reliably observed neutron stars with masses of about two solar masses, the presence of quark matter is found to be linked to the behaviour of the speed of sound c_s in strongly interacting matter. If the conformal bound (c_s)^2<1/3 is not strongly violated, massive neutron stars are predicted to have sizable quark-matter cores. This finding has important implications for the phenomenology of neutron stars, and affects the dynamics of neutron star mergers with at least one sufficiently massive participant.

nuclear model calculations. However, the matter in the interior of maximally massive, stable neutron stars exhibits characteristics of the deconfined phase, which we interpret as evidence for the presence of quark-matter cores. For the heaviest reliably observed neutron stars 5,6 with mass M ≈ 2M ⊙ , the presence of quark matter is found to be linked to the behaviour of the speed of sound c s in strongly interacting matter. If the conformal bound c 2 s ≤ 1/3 7 is not strongly violated, massive neutron stars are predicted to have sizable quark-matter cores.
This finding has important implications for the phenomenology of neutron stars, and affects the dynamics of neutron star mergers with at least one sufficiently massive participant.
Observations of neutron stars (NSs) inform us about the properties of matter inside their cores in an indirect way. To translate them to statements about NS matter requires modeling strongly interacting matter all the way from the crust to the highest densities reached inside the stars. The lack of accurate first-principles predictions at densities beyond the nuclear matter saturation (baryon number) density n 0 ≈ 0.16/fm 3 has so far prevented the determination of the phase of matter inside NS cores, and it is unlikely that the question will be answered based on gravitational wave (GW) data alone, at least in the near future 8 . Nevertheless, recent observations are beginning to offer empirical constraints so strong that a model-independent approach to the problem has become feasible.
The equation of state (EoS) of NS matter-the relation p(ǫ) between the pressure and energy density of beta-equilibrated matter interacting under quantum chromodynamics (QCD) at temperature T = 0-is known in two opposing limits. From the well-studied NS crust region 9 to the density n CET ≡ 1.1n 0 , where matter resides in the hadronic-matter phase, modern nucleartheory machinery, such as chiral effective field theory (CET), provides the EoS to good precision, currently better than ±24% 10, 11 . In the opposite limit of very high densities, perturbative-QCD (pQCD) techniques, rooted in high-energy particle phenomenology and built on deconfined quark and gluon degrees of freedom 12,13 , become accurate, providing the quark-matter EoS to the same accuracy at densities n 40n 0 ≡ n pQCD .
In the above two limits, QCD matter is known to exhibit markedly different properties. Highdensity quark matter is approximately scale invariant, or conformal, whereas in hadronic matter the number of degrees of freedom is much smaller and, additionally, scale invariance is violated by the breaking of chiral symmetry. These qualitative differences are reflected in the values taken by different physical quantities. The speed of sound takes the constant value c 2 s = 1/3 in exactly conformal matter, and slowly approaches this number from below in high-density quark matter 12 .
By contrast, in hadronic matter the quantity varies considerably: below saturation density, CET calculations indicate c 2 s ≪ 1/3, while at higher densities most hadronic models predict max(c 2 s ) 0.5. The polytropic index γ ≡ d(ln p)/d(ln ǫ) on the other hand obtains the value γ = 1 in conformal matter, while both CET calculations and hadronic models generically predict γ ≈ 2.5 around and above saturation density. Finally, the number of degrees of freedom is reflected in the pressure normalized by that of free quark matter (the Fermi-Dirac, or FD, limit), p/p FD 12 Fig. 1, which roughly coincides with the energy density inside free nucleons and to the location of the deconfinement transition at high temperatures 22, 23 . At even higher densities, the interpolated EoSs approach the pQCD predictions for c 2 s , p/p FD , and γ 12 .
On each of the EoS lines, the density reached in the centers of 1.4M ⊙ NSs is marked with solid blue (interpolated EoSs) or empty cyan (nuclear EoSs) diamonds. The significant overlap between the two distributions shows that the material properties inside 1.4M ⊙ stars are consistent with a description in terms of hadronic degrees of freedom. The same is, however, not true for the majority of maximally massive stars (filled red and empty magenta circles), for which the two families of points are clearly separated in In order to make the above observations somewhat more quantitative, it would clearly be valuable to establish a connection between the physical phase of QCD matter and its EoS. To this end, we note that our results from Fig. 2 and the distinct values the polytropic index obtains in nuclear and quark matter calculations both suggest using the values of γ as a good approximate criterion. Given that γ = 1.75 is both the average between its pQCD and CET limits and very close to the minimal value the quantity obtains in viable hadronic models (see Fig. 2 and our discussion in the Methods section), we are led to choose the following criterion for separating hadronic from quark matter: given an interpolated EoS, the smallest density from which γ < 1.75 continuously to asymptotic densities is identified with the onset of quark matter. We emphasize, however, that this is only an approximate rule to guide our analysis, and not a robust or rigorous result.
As a first application of the above criterion, for NSs with M = 1.4M ⊙ we find that the central polytropic index always satisfies γ 2, implying that the stars are composed of hadronic matter as expected. In contrast, maximally massive stars have γ values much closer to unity, indicating that they typically contain quark matter. In Fig. 3, we display the sizes of quark cores in the latter NSs. The core has a significant extent, M core > 0.25M ⊙ , for all those EoSs that satisfy c 2 s < 0.5.
However, for extreme EoSs in which the speed of sound almost reaches that of light, the core may be significantly smaller or even absent. Note that this is consistent with the fuller picture given by If the maximal value of c 2 s exceeds 0.7, we find a small class of EoSs where even maximally massive stars do not contain quark cores according to our criterion. We find that each of these EoSs exhibits an interval in ǫ where γ < 0.5, which destabilizes the star. This corresponds to a rapid change in the EoS, and is practically indistinguishable from a first-order phase transition. The minimal latent heat (i.e. the extent of the interval with γ < 0.5) required for the destabilization is (∆ǫ) lat > 130 MeV/fm 3 , corresponding to a relative discontinuity (∆ǫ) lat /ǫ > 0.2 at the beginning of the transition. We thus find that in order for all stable NSs to be composed of hadronic matter alone, the EoS must both significantly violate the conformal limit and feature a sufficiently strong phase transition. Finally, two-solar-mass stars contain a quark core for all EoSs that satisfy c 2 s < 0.4, irrespective of the properties of the phase transition; for subconformal EoSs, featuring c 2 s < 1/3 at all densities, the radius of the core R ≈ 6.5 km is roughly half of the entire star's radius. By contrast, if the EoS supports substantially higher maximal masses M max > 2.25M ⊙ , quark cores are absent in 2M ⊙ stars (see Extended Data Fig. 4).
In conclusion, our model-independent analysis has demonstrated that the existence of quark cores in massive NSs should be considered the standard scenario, not an exotic alternative: for all stars to be made up of hadronic matter, the EoS of dense QCD matter must be truly extreme. This view is also consistent with recent NS radius measurements, which are compatible with the larger radii predicted by the less extreme EoSs (see discussion under Methods, and Extended Data Fig. 5).
Note, however, that our analysis does not preclude the possibility of massive, purely hadronic stars even with less extreme EoSs, as quark cores may appear only at very high masses, even beyond 2M ⊙ . Nuclear matter EoSs that predict purely hadronic 2M ⊙ stars (see, e.g., 24-28 ) may therefore be compatible with our results until very high densities.
The existence of quark cores in at least some NSs, and that the nucleation of quark matter begins so close to the maximum-mass limit, may have very interesting observable consequences.
In NS mergers, the core may lead to shock waves reflecting from the quark-hadron interface inside hypermassive NSs. This may be particularly amplified if the conformal limit is strongly violated in hadronic matter, leading to large differences in the speeds of sound between the phases. In addition, the onset of the transition may give rise to dissipation during the merger in the form of a large effective bulk viscosity that may lead to an enhanced damping of the ringdown. Importantly, both of these have the potential to lead to observable effects in GW signals from NS mergers and the associated EM counterparts.
Finally, our results are systematically improvable with more observations. For example, there are several candidates for NSs with very large masses (see e.g. ref. 29 ). If even one of these stars turns out to have a mass significantly larger than 2M ⊙ , this would impose strong new constraints on the EoS and e.g. imply that the conformal bound must be broken. Similarly, with many binary-NS merger observations currently recorded by LIGO/Virgo, the current limits on the tidal deformability will inevitably become tighter, enabling additional improvements to our analysis. With these advances and the road map laid out in our work, further significant progress in understanding the nature of ultra-dense matter inside NSs can be expected in the near future.

Methods
New speed-of-sound interpolation The starting point of the new interpolation method is to consider the squared speed of sound c 2 s as a function of the baryon chemical potential µ, and use this quantity to construct all other thermodynamic functions, in particular the pressure p(µ). In practice, the speed of sound is first integrated from the CET matching point n CET = 1.1n 0 to higher densities to give the baryon density where µ CET is the baryon chemical potential corresponding to the density n CET , that is, n CET ≡ n(µ CET ). This result is then further integrated to arrive at the pressure, where p CET ≡ p(µ CET ).
The above relations must be solved numerically in general, but in the following simple case that we have implemented in our analysis, they may be dealt with analytically. Namely, we first take sequence of N p pairs with µ 1 = µ CET , µ Np = 2.6 GeV, and µ i−1 < µ i < µ i+1 for all other i. We then construct a c 2 s curve as a piecewise-linear function connecting these points; that is, for each i = 1, . . . , N p − 1, At the matching points µ 1 and µ Np , we require p and c 2 s to match the corresponding values given by the CET and pQCD EoSs, respectively. In addition, we take n to be continuous at each matching point, but note that our construction allows for EoSs that mimic discontinuous first order transitions arbitrarily closely.
For a given N p , we have N p − 2 independent matching chemical potentials µ i and N p − 2 independent speed-of-sound points c 2 s i , from which one of both is determined through matching to the high-density EoS, leaving 2N p − 6 parameters for given low-and high-density EoSs. If we instead write this in terms of the number of interpolating segments N ≡ N p −1, the result becomes 2N − 4. This is one free parameter fewer than the number of free parameters needed to define a polytropic EoS composed of the same number of segments 16 .
The above procedure is used to construct individual EoSs by choosing N = 3, 4, 5, and then randomly picking values for the matching points µ i , speeds of sound c i , and the pQCD parameter X pQCD 30 . The parameter values are taken from uniform distributions µ i ∈ [µ CET , 2.6 GeV], 1,4], in addition to which we choose roughly the same number of the "hard" or "soft" nuclear EoSs of ref. 14 . Finally, we vary the extreme EoSs in the ǫ, p plane within each c 2 s band plotted in our paper, to ensure that we satisfactorily probe the size of these regions. This leads to the ensemble studied above, which consists of approximately 570,000 individual EoSs.
Roughly 160,000 of these fulfill the astrophysical constraints described in the main text, while approx. 70,000 of the allowed EoSs contain at least one first order phase transition. We have carefully made sure that these ensemble sizes are sufficiently high, so that our results are stable with respect to increasing the number of EoSs.
Finally, we note that while the interpolation method described above is genuinely new, a number of related articles have recently appeared, in which the NS matter EoS has been constructed starting from the speed of sound 31-36 . While most of these works introduce a nontrivial ansatz function for the quantity, thus being more restrictive than our approach, in ref. 31 the speed of sound is allowed to behave in a more general way. The main difference between the EoS bands constructed in this reference and our current work originates from our high-density pQCD constraint, which effectively forces the EoS to be softer at high densities.

Comparison of different interpolations
To quantify the potential bias introduced into our results by the selection of the speed-of-sound interpolation method, we compare our EoS and MR ensembles to ones obtained with the following two schemes: 1. A piecewise polytropic interpolation of the pressure as a function of baryon density, in terms of Chebyshev polynomials.
Both of these interpolation methods have been abundantly discussed in the literature 14  Again, since the spectral method leads to smoother interpolations, it is natural that it does not allow these rapidly changing EoSs.
Another difference between the interpolation schemes is that the polytropic interpolation does not allow for as massive stars as the other two. We attribute this to the fact that in order to achieve very large maximal masses, the EoS needs to stay very stiff, c s ≈ 1, throughout an extensive density window, which is difficult to realize with polytropic interpolation functions. This difference between different interpolations is somewhat ameliorated when upper limits are placed on the tidal deformability.
Polytropic index and its relation to the phase of QCD matter As stated in the main text, our criterion for identifying the phase of QCD matter in NS cores is based on analyzing the behavior of the polytropic index γ ≡ d(ln p)/d(ln ǫ), i.e., the slope of the EoS in Fig. 1 and Extended Data Figs. 1, panel a; 2; and 3. Here, we comment on the physics behind this statement, and explain our choice to identify the presence of quark matter using as the quantitative criterion γ < 1.75 continuously up to asymptotic densities.
Matter that exhibits exact conformal symmetry, i.e., does not possess intrinsic mass scales, is characterized by γ = 1 independent of the strength of the coupling. This is so because in the absence of any dimensionful parameters, the energy density and pressure must to be proportional to each other, leading to γ = 1. This symmetry can also be shown to lead to a speed of sound squared c 2 s = 1/3 for the system.
In low-and moderate-density QCD matter, it is well known that the ground state does not exhibit the approximate chiral symmetry of the underlying Lagrangian (see e.g. ref. 39 for details).
This spontaneous breaking of the symmetry leads to the emergence of the fundamental scales of nuclear matter, such as hadron masses, and scale-dependent interactions. These mass scales lead to a highly non-conformal behavior for the EoS, which is reflected in the polytropic index taking large values, typically in excess of 2, in viable models of high-density nuclear matter.
Collections of γ values that different nuclear physics models predict are available through the related adiabatic index Γ = ǫ+p ǫ γ, tabulated e.g. in Table III of ref. 40 and plotted in figure 5.9 of ref. 41  In high-density quark matter, on the other hand, the underlying approximate chiral symmetry of QCD is restored, and, as a result, the system exhibits approximate conformal symmetry. Minor violations of conformal behavior arise from the masses of the up, down, and strange quarks, which are, however, very small compared to the nucleon masses. Moreover, the interactions between quarks and gluons also lead to a mild breaking of conformal symmetry in the quark-matter phase, manifesting as a logarithmic dependence of the strong coupling constant on the baryon chemical potential. To a good accuracy, however, quark matter always behaves as an ultrarelativistic gas of interacting quarks and gluons, which becomes even more pronounced in the controlled, perturbative high-density region of the QCD phase diagram, where the polytropic index γ quickly approaches unity.
In the most nontrivial density range near the deconfinement transition, QCD matter evolves from the highly nonconformal hadronic behavior to one characteristic of quark matter. This transition may take place either as a discontinuous jump in energy density, in which case the value γ = 1.75 may never be reached, or in a smooth crossover manner, whereby a crisp phase identification may not always be feasible (then there may even exist an overlap region where the system can be described both in terms of nuclear and quark degrees of freedom). In either case, our results indicate that the cores of typical pulsars and maximal-mass NSs very likely do not reside here, but rather safely belong to the nuclear and quark matter regimes, respectively. This is reflected in the fact that our qualitative conclusions are not sensitive to the exact choice of the critical polytropic index as long as it resides between the hadronic and quark matter regimes. Indeed, only our detailed numerical conclusions would be somewhat modified, should we vary the number γ = 1.75 in a moderate way. shows that the EoSs that make up the boundaries of the our EoS band must exhibit both very large speeds of sounds as well as rapid changes in material properties.

Analysis of EoS smoothness
For our analysis of the central values of γ in stars of different masses, we have used ∆ln ǫ > 0.5, which has a minor effect on the global characteristics of the EoS family. In particular, this cut has a minor effect on the mass-radius relation for stars above 1.4M ⊙ , shifting the extremes of the allowed radius for a given mass by about 0.3 km at most. Moreover, we note that for completeness we have also allowed the first inflection point to approach n CET without limit, finding that all the results presented in the main text remain unchanged.
Comparison with recent mass and radius constraints Our ensemble of EoSs can also be transformed to the M-R plane in order to compare its behaviour to recent radius observations of NSs. In Extended Data Fig. 5, we overlay a few representative X-ray MR-measurements on top of our family of MR-curves that are obtained from the EoS ensemble of Fig. 1 Data Availability Source data are available for this paper. For three tabulated sample EoSs and the boundary regions for the p(ǫ) (Fig. 1) and MR regions (Extended Data Fig. 5), please consult the Supplemental files and Source Data files. The authors will in addition be happy to provide more EoS tables upon request.
Code availability The code used to construct the interpolated EoSs is available from the authors upon request.
Author contributions The authors conducted the research together and participated in the preparation and revision of the manuscript. The speed-of-sound interpolation was implemented by T.G. and the spectral interpolation by E.A., both based on an earlier code written by A.K. For comparison, the black lines stand for the different hadronic EoSs we have obtained from refs. 9,20,21 . Finally, the light blue regions correspond to the CET and pQCD EoSs of 12,14 , and the rough location of the deconfinement transition in hot quark-gluon plasma, ǫ QGP , is indicated for illustrative purposes.      the time-evolving X-ray burst spectra (corresponding to the low-mass X-ray binary