Quark-matter cores in neutron stars

Eemeli Annala, Tyler Gorda, Aleksi Kurkela, Joonas Nättilä, and Aleksi Vuorinen Department of Physics and Helsinki Institute of Physics, P.O. Box 64, FI-00014 University of Helsinki, Finland Department of Physics, University of Virginia, Charlottesville, Virginia 22904-4714, USA Theoretical Physics Department, CERN, Geneva, Switzerland and Faculty of Science and Technology, University of Stavanger, 4036 Stavanger, Norway Nordita, KTH Royal Institute of Technology and Stockholm University, SE-10691 Stockholm, Sweden


I. INTRODUCTION
The theory governing the strong nuclear force, Quantum Chromodynamics (QCD), predicts that at sufficiently high energy densities nuclear matter undergoes a deconfinement transition to a new phase composed of quarks and gluons. An experimental program with ultrarelativistic heavy-ion collisions, carried out at the SPS, RHIC and LHC colliders, has led to the conclusion that high-temperature quark-gluon plasma (QGP) has indeed been successfully produced in laboratory conditions [1]. According to lattice Monte-Carlo simulations, at small baryon densities the crossover transition from the hadronic phase to the QGP takes place at an energy density of approximately 500 MeV/fm 3 [2,3].
A complementary system where deconfined matter may exist lies within the cores of neutron stars (NSs)the densest astrophysical objects in our Universe-where gravity compresses matter to energy densities comparable to those reached in heavy-ion collisions [4,5]. However, whether the conditions in the centers of NSs are extreme enough to lead to the formation of a core consisting of a new state of deconfined matter, cold and dense quark matter (QM), remains an open question.
The past decade has witnessed remarkable advances in NS observations: thanks to the discovery of extremely massive NSs [6,7], qualitative improvements in the understanding of systematic uncertainties in x-ray radius measurements [8][9][10][11][12][13][14], as well as the famous LIGO/Virgo detection of gravitational waves (GWs) originating from the NS-NS merger GW170817 [15], robust and increasingly stringent constraints have been placed on the EoS of NS matter [16][17][18][19][20][21][22][23][24][25][26][27]. What has thus far hindered firm conclusions about the presence of QM inside NSs is the lack of accurate first-principles predictions for the properties of QCD matter at high baryon densities. The numerical lattice-simulation techniques that form the theoretical foundation of the description of hot QGP fail in the cold, baryon-rich conditions of NSs due to the infamous sign problem [28]. While steady progress has occurred both in the theoretical description of moderate-density nuclear matter (see, e.g., [29,30]) and ultrahigh-density QM [31][32][33], no reliable results exist in the crucial regime between approximately one and ten nuclear saturation densities. Combined theoretical and experimental studies of hot QGP have taught that whether matter is in a partonic or hadronic phase is ultimately determined by the thermodynamic and dynamical properties of the system that reflect the underlying degrees of freedom. Lattice simulations have demonstrated that at high temperatures the hadronic and partonic phases are separated by a smooth crossover transition [2,3]. However, despite the lack of a discontinuous transition, the two phases are clearly iden-tifiable by, e.g., the qualitative behavior of the Equation of State (EoS). At low temperatures, the EoS of dilute QCD matter follows the predictions of the hadronresonance-gas model [34,35], while at high temperatures it resembles the nearly conformal EoS of resummed perturbative QCD (pQCD), building on quark and gluon degrees of freedom [36][37][38].
For the high densities and low temperatures realized inside NSs, there is no strict reason for why the nature of the deconfinement transition should qualitatively change. While it is naturally possible that a line of first-order transitions begins at a tricritical point, as predicted by a number of model calculations (see, e.g., the recent review [39]), the latent heat of such a transition would in any case be limited by its non-observation at smaller baryon densities [40]. In either case, the presence of QM inside the cores of NSs does not rely on a first-order phase transition separating it from hadronic matter (HM), but is ultimately determined by whether or not there is a rapid qualitative change in the properties of dense QCD matter, and whether or not the central densities of NSs reach this (pseudo)critical density. Concretely, the following scenarios exhaust the possibilities: 1. The deconfinement transition takes place so slowly that it is altogether impossible to delineate at which density HM turns to QM.
2. The transition proceeds as a first order one or a rapid crossover at an energy density c , but the central densities in the most massive stable NSs do not exceed it, so no QM is present inside any star.
3. The central densities of at least some stable NSs lie clearly above the transition region, so that QM undeniably exists in them and the stars become hybrid stars.
Note that QM may also be found inside so-called Twin Stars [41], which we will not discuss in this work.
In the present paper, we combine the most recent theoretical and observational inputs to maximally constrain the material properties of NS matter, and use this information to discriminate between the three scenarios listed above. By applying a set of novel and versatile interpolation functions that allow for complex structures in the EoS, yet respect the known low-and high-density limits of the quantity, we discover a range of possible outcomes that Nature can choose from. Interestingly, we find that Scenario 1 is not present in the set of viable EoSs, and moreover that the central densities of maximally massive NSs always reach c , reducing the question about the presence of QM inside their cores to the properties of the deconfinement transition. A closer inspection of the individual EoSs further reveals that very specific and extreme conditions need to be met in order for all stars be composed of HM alone. In particular, should the maximal speed of sound in the physical EoS be moderate, the QM cores of massive NSs are typically sizable; see Fig. 1.

II. METHODOLOGY
The EoS-or energy density as a function of pressure, (p)-of beta-equilibrated strongly interacting matter at non-zero baryon number density n and zero temperature T is known reliably in two opposing limits. From the wellstudied NS crust region [42] to densities slightly above the nuclear saturation density n 0 = 0.16/fm 3 , nuclear theory methods offer a robust way of determining the properties of nuclear matter. As described in [30], the uncertainty in the state-of-the-art Chiral Effective Theory (CET) EoS can be estimated by varying a cutoff scale Λ, which indicates an uncertainty of ±24% in the pressure at a baryon density n CET ≡ 1.1n 0 , growing thereafter rapidly with increasing density. In the opposite limit of very high n, the asymptotic freedom of QCD guarantees that the strong coupling constant α s ( 1/4 ) ∝ 1/ log( ) is small enough to admit a weak coupling expansion of the QM EoS, which has so far been calculated to O(α 3 s log 2 α s ) [33]. The systematic uncertainty in the pQCD EoS is quantified via a variation of the renormalization scaleΛ by a factor 2 around its fiducial value, indicating convergence similar to the CET one at baryon chemical potentials µ B 2.6 GeV, or densities n 40n 0 ≡ n pQCD .
In the intermediate-density range n CET < n < n pQCD , no reliable theoretical calculations are available to offer guidance about the behavior of the EoS. At the same time, it is now well-established that knowledge of the lowand high-density limits of the quantity restricts its behavior at all densities-a fact that is most easily demonstrated by introducing a set of basis functions to interpolate the EoS through the intermediate-density region [16,17,43,44]. While many sets of basis functions exist that allow complex structures in the EoS, any particular choice may introduce a bias in the results [45]. To account for this, we have used multiple interpolation methods and compared the corresponding results: in terms of Chebyshev polynomials [46], 3. A piecewise-linear interpolation of the squared speed of sound, c 2 s (µ B ) = dp d . Of these schemes, the first two have been abundantly discussed in the literature [16-18, 43, 44, 46], but the third one is new and will be introduced in detail in Appendix A (see also [47][48][49] for related approaches).
Our first new result concerns the comparison of the EoS and mass-radius (MR) bands composed using the above basis functions, displayed in Figs. 8 and 9 of Appendix B. As will be discussed there in detail, we find the results obtained to be mostly in good agreement with each other, after setting the numbers of free parameters in the different interpolation schemes roughly equal and imposing observational constraints from NS properties (see also Sec. III.A). This illustrates that they are not subject to a significant bias arising from the choice of basis functions, and a posteriori strengthens the conclusions made in previous works [16-18, 43, 44, 46]. As the three interpolations agree, in the following, we choose to use the speed-of-sound interpolation. We note that the added benefit of this method is that it allows one to keep track of the stiffness of the EoS in a natural way.

III. CONSTRAINING THE NS-MATTER EOS
The next two sections are devoted to a detailed analysis of our ensemble of NS-matter EoSs, constructed with the speed-of-sound method. As detailed in Appendix A, the approximately 570.000 EoSs are built from randomly generated functions c 2 s (µ B ), containing up to 5 linear intervals, whereafter we vary the outlier EoSs to make sure that the boundaries of the EoS band are stable. Note that while we do not add discontinuous first-order transitions to our EoSs by hand, our interpolation functions allow crossover transitions that may be arbitrarily strong, thus closely mimicking discontinuous phase transitions and mixed phase constructions [50].

A. Properties of the EoS band
In Fig. 2, we display our ensemble of NS-matter EoSs obtained with the speed-of-sound interpolation method. In deriving the result, we have required that the EoSs support a 1.97M NS [6,7] and that the tidal deformability Λ for a 1.4M star satisfy 70 < Λ(1.4M ) < 580, consistent with the LIGO/Virgo bound from the GW170817 observation [18]. As noted earlier (see, e.g., [16,44]), the two-solar-mass constraint forces the EoS to be relatively stiff at low densities, which is reflected in the rapid rise of the interpolation functions for the pressure as a function of energy density. At the same time, the constraint on Λ(1.4M ) sets an upper limit for the stiffness, constraining the EoS band in a complementary direction.
While the astrophysical observations significantly constrain the behavior of the EoS in the intermediate-density region, and the new band is more restrictive than, e.g., that of [16], the range of allowed EoSs still remains relatively wide. A partial reason for this is the high versatility of our interpolation method, which allows for very complex structures and extreme states of matter, some of which are unlikely to appear in Nature. Instead of imposing a theoretical bias and restricting the set of EoSs by hand, we have chosen to classify the functions based on their extremeness as quantified by the maximum value that the speed of sound reaches and the level of fine structure that each EoS contains.
In Fig. 2, the speed-of-sound classification is performed following a coloring scheme where EoSs corresponding to a lower maximal value of c 2 s are drawn on top of the higher ones. While we are not aware of a proven theo- with the speed-of-sound interpolation method introduced in this paper. The color coding refers to the maximal value that c 2 s reaches at any density, while the black lines denote the extrapolations of the low-and high-density theoretical bands to higher/lower densities [33,56]. The rough location of the deconfinement transition in hot QGP is indicated as QGP.
rem that would exclude speeds of sound exceeding the conformal value c 2 s = 1/3 (see, however, [51] for an attempt in this direction), we note that the bound appears to be a very nontrivial one to break. In hot QGP, nonperturbative lattice simulations have shown that the speed of sound remains subconformal [52], and in QCD matter at asymptotically high energy density the quantity is known to approach the conformal limit from below [31]. In holographic calculations the bound has been violated, but only in finely tuned constructions that do not directly correspond to quantum field theories realized in Nature [53,54]. As discussed in [55], having c 2 s > 1/3 furthermore corresponds to matter in which the number of degrees of freedom decreases as a function of energy density, which strongly goes against the partonic picture of hadrons arising from QCD. Based on these considerations, we conclude that there is a strong theoretical reason to expect that the speed of sound never exceeds the conformal value by a sizable amount in QCD matter. As seen from Fig. 2, excluding those EoSs for which the conformal limit is strongly violated, say c 2 s > 0.6, would lead to significantly tighter limits for the allowed EoSs.
Another way in which some of the EoSs generated by the speed-of-sound interpolation method are extreme is that the interpolation functions allow for very quick changes in the material properties of the medium in arbitrarily small density windows. While such versatility is in principle a desirable feature of the interpolator, these structures are clearly not very likely to appear in Nature. To quantify the level of local structure in our EoSs, we classify them according to the smallest (logarithmic) energy density interval where structures appear. In practice, this is implemented by demanding that the energy densities at two successive inflection points i and i+1 where the speed of sound changes its behavior, satisfy Note that imposing this constraint does not exclude, e.g., first order phase transitions or rapid crossovers. Fig. 3 demonstrates the effect of imposing various lower limits for the parameter ∆ln on the band of allowed EoSs. We observe that placing small limits on it (∆ln 1) affects the band mainly around the matching points where the EoS is best known but does not have a significant effect at intermediate densities. However, somewhat larger values (∆ln 1) begin to significantly constrain the band at all densities. This shows that the EoSs that make up the boundaries of the band must have both very large speeds of sounds as well as rapid changes in material properties.

B. Transition to quark matter
Returning to the full EoS ensemble of Fig. 2, we next inspect the overlaid extrapolations of the low-and highdensity EoS bands, taken from [33,56]. While the reliability of the low-density EoS becomes questionable beyond the density of 1.1n 0 , it is remarkable how the interpolations accurately follow the trend set by the theoretical calculation, corresponding to γ ≡ d(log p)/d(log ) ≈ 2.5, until considerably larger densities. Shortly thereafter, the rapid rise of energy density slows down, forming a visible kink in the EoS band at an energy density of c ≈ 400 − 700 MeV/fm 3 , corresponding roughly to the energy density inside free nucleons [57] and to where the chiral phase transition occurs at high temperatures [2,3]. The nascent onset of the kink has been observed in several works with and without the high-density pQCD constraint (see, e.g., [18,43]), and is essentially caused by the requirement that the EoS stay subluminal.
A ballpark estimate for the location of the kink can be obtained as follows. Below the kink, the EoS roughly follows a trend set by p ∼ 0.01 0 ( / 0 ) γ , with 0 ≈ 100 MeV/fm 3 and γ ≈ 2.5. When γ is larger than one, the ratio of p/ grows rapidly as a function of energy density. The speed of sound is correspondingly given by where the logarithmic derivative (i.e. γ) corresponds to the slope in Fig. 2. To maintain c 2 s < 1, this parameter must begin to change at the latest when p/ ∼ 1/γ, which corresponds to ≈ 1 GeV/fm 3 .
At even higher densities, the pressure begins to rise more slowly as a function of energy density, so that it may fulfill the high-density pQCD constraint. The requirement of reaching the pQCD EoS restricts the slope both from below and above. Again, we observe that the interpolations above the kink accurately follow the trend set by a naive extrapolation of the upper limit of the pQCD band with γ ≈ 1 [74]. It is worth stressing that the divergence of the lower limit of the uncertainty band need not signal a qualitative change in the properties of the matter, but merely indicates that the naive perturbative series may not be the most efficient way to organize the weak coupling calculation.
The rapid change in the material properties of the system, as indicated by Fig. 2, together with the theoretical interpretations of the two parts of the EoS band, strongly suggest the following picture: below the kink, NS matter is described by γ 2.5 clearly corresponding to HM, while above it, we have γ ≈ 1, corresponding to QM. The fact that we observe a clear separation of the two phases moreover rules out Scenario 1 of Sec. I, leaving as the key question whether the densities reached inside physical NSs are high enough to form a sizeable QM core.  figure, we observe that the central densities of maximally massive stars appear to all lie to the right of the kink, suggesting that at least small amounts of QM can be expected to be found inside them. At the same time, there are a number of EoSs for which the two-solar-mass and in particular 1.44M stars do not lie above the kink, so that the fate of QM inside them is less clear. In addition, it should be noted that while all the individual EoSs exhibit a transition to the QM phase, the location of the (pseudo)critical density varies between them, and may in some cases be somewhat displaced from the center of the kink. For this reason, a case-by-case analysis of the EoS family is clearly needed. In order to inspect the EoSs case by case, we next investigate the ranges of polytropic indices γ found at the centers of NSs with different masses, recalling their distinct values in the HM and QM phases. For this analysis, we need to exclude from our full ensemble a small number of EoSs that contain very sharp local structures, leading to rapidly varying γ values at the centers of NSs that do not reflect the overall trend of the EoS in question. A sufficient cut having the desired effect is ∆ln > 0.5, which according to Fig. 3 has a minor effect on the global characteristics of the EoS family [75].
The result of the polytropic index analysis is displayed in Fig. 5, where we reproduce the ranges of polytropic indices γ found at the centers of M max and 1.44M NSs for different maximal values of c 2 s . For NSs with M = 1.44M , we always find the central polytropic index to satisfy γ 2, which clearly corresponds to the HM phase. On the other hand, for the maximally massive stars we typically find γ values slightly above unity, indicating that the matter is in the QM phase. Fig. 1 displays the size of the quark core, which we define as the continuous region at the center of the NS where γ remains below 1.75 (denoted by the dashed vertical line in Fig. 5). The core has a significant extent, M core > 0.25M , for all those EoSs that satisfy c 2 s < 0.5. However, for EoSs that strongly violate the conformal limit, the core may be significantly smaller or even absent.
If the maximal value of c 2 s exceeds 0.7 (or 0.5 for ∆ln = 0), we find a small class of EoSs which do not lead to QM cores even for maximally massive stars. These EoSs correspond to that part of the M max cloud in Fig. 5 that extends to the right of the dashed vertical line. Inspecting this set of EoSs further, we find that they all exhibit a first order phase transition, which we define as an interval in where γ < 0.5, where the pressure is approximately flat as function of energy density. Further analysis confirms that in these cases it is indeed the phase transition itself that destabilizes the star. To study how large latent heats are required for the destabi- lization, we inspect this set of EoSs with different values of ∆ln > 0.05, thereby making sure that the size of the latent heat is not limited by the smoothing procedure. This analysis shows that a phase transition with (∆ ) lat > 130 MeV/fm 3 , or (∆ ) lat / > 0.2 at the beginning of the transition, can prevent the formation of a quark core for EoSs for which max(c 2 s ) > 0.5. Finally, we find that two-solar-mass stars contain a quark core for all EoSs that satisfy c 2 s < 0.4 (as well as many with c 2 s > 0.4). The respective sizes of these cores are displayed in Fig. 6 together with the maximal masses obtained with the same EoSs. We observe that for subconformal EoSs, where the maximal masses are close to 2M , the two-solar-mass stars contain large quark cores of R ≈ 6.5 km. We note that for these EoSs, the MR measurements of the core are nearly identical to those within the maximal mass NSs in Fig. 1. On the contrary, for those EoSs that lead to substantially higher maximal masses, M max > 2.25M , the quark cores are absent in 2M stars, indicating that the formation of a soft core quickly leads to a destabilization of the star even in the absence of a strong phase transition.

V. CONCLUSIONS
Although increasingly precise constraints have been placed on the EoS of neutron-star matter [16,17,43,44], the microscopic composition of QCD matter deep inside NS cores has so far been addressed only within the context of specific phenomenological models. In the present paper, we have shown that current astrophysical and theoretical constraints are starting to be restrictive enough so that this question can be addressed in a more robust, model-independent way. In particular, we have demonstrated that the NS-matter EoS has a clear two- phase structure: at low densities, it is characterized by a hadronic polytropic index γ ≡ d(log p)/d(log ) 2.5, while at high densities we have γ ≈ 1, corresponding to nearly conformal quark matter. The transition between the two phases happens at an energy density of 400 − 700 MeV/fm 3 , comparable to that where quarkgluon plasma is generated in heavy-ion collisions. We have observed that these results are moreover independent of the choice of basis functions used for interpolating the EoS through the intermediate-density regime.
We have studied the constrained EoS as a function of the maximum speed of sound it attains, motivated by the fact that the conformal bound c 2 s ≤ 1/3 appears to be nearly universally respected in physical systems. We have identified a large number of EoSs that are not only consistent with this bound but also all other observational and theoretical constraints. We note that while this conclusion is seemingly ostensibly different from that of [55], in their analysis the authors of this work too find a number subconformal EoSs that lead to 2M stars.
We also note that the EoSs with non-extreme speeds of sounds are in addition in good agreement with the most recent simultaneous NS mass-radius measurements. This can be seen from Fig. 7, where we compare the MRrelation stemming from our EoSs to the most recent measurements corresponding to NSs in the low-mass x-ray binary systems 4U 1702−429, 4U 1724−307, and SAX J1810.8−2609, obtained with the x-ray-burst cooling-tail method [12,13]. We emphasize that this data was not used to constrain our ensemble of EoSs. In addition, we note that low speeds of sound are consistent with bounds from SSS17a and GRB170817a, the EM counterparts of GW170817 [63][64][65][66][67][68][69][70], suggesting Λ > 300 [58] and M max < 2.16M [59][60][61][62].
An important finding of ours is that the cores of stars with different masses have strikingly different properties. On the one hand, we find that, e.g., typical binary pulsars with M ≈ 1.4M do not reach central energy densities  Fig. 2, together with three recent simultaneous MR-measurements [12,13]. The dark regions on the left and top correspond to regions excluded by constraints based on the EM counterpart of GW170817, Λ > 300 [58] and Mmax < 2.16M [59][60][61][62]. The dashed red curve denotes the masses below which there are no NSs containing QM cores.
high enough for a QM core to form irrespective of the maximal speed of sound attained; indeed, below the thick red dashed curve in Fig. 7, no NSs contain QM cores of any kind. At the same time, maximally massive stable stars contain (typically large) quark cores unless the EoS is truly extreme with c 2 s > 0.7 and a phase transition strong enough to destabilize the star.
Finally, we find that if the maximal mass of NSs is smaller than 2.25M , the two most massive NSs known to date may contain very large QM cores up to R core ≈ 7 km; in particular, if c 2 s < 0.4, then the 2M NSs contain at least a 3 km quark core. That the subconformal EoS predicts a large QM core to be present in the known J1614−2230 and J0348+0432 NSs may open up a phenomenological way of answering an open fundamental problem in QCD concerning whether the speed of sound exceeds the conformal bound: if there is no quark core inside these two stars, then we know that the bound has been violated in QCD matter.
The existence of massive quark cores in at least some physical NSs-or that the nucleation of QM begins so close to the maximum mass limit-may have interesting observable consequences. In NS mergers, currently under intense observational and theoretical scrutiny [71], the core may lead to shock waves reflecting from the QM-HM interface inside hypermassive NSs. This may be particularly amplified, if the conformal limit is strongly violated in the HM phase, leading to large differences in the speeds of sound between the two phases. In addition, the onset of the transition may give rise to increased dissipation in the form of a large effective bulk viscosity that may lead to an enhanced damping of the ringdown [72]. Importantly, both of these have the potential to lead to observable effects in NS merger GW signals and the associated kilonova explosions and gamma-ray bursts.

Acknowledgments
We would like to thank Ingo Tews for helpful discussions and collaboration in early stages of this work. In addition, we acknowledge useful discussions with Niko Jokela, Simonetta Liuti, Anton Rebhan, Sanjay Reddy, and Kent Yagi. The work of EA, TG, and AV has been supported by the European Research Council, grant no. 725369. In addition, EA gratefully acknowledges support from the Finnish Cultural Foundation and TG from the U.S. Department of Energy Grant No. de-sc0007984.

Appendix A: 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 µ B , and use this quantity to construct all other thermodynamic functions, in particular the pressure p(µ B ). In practice, the speed of sound is first integrated from the CET matching point n CET = 1.1n 0 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 a sequence of N p pairs with µ 1 = µ CET , µ Np = 2.6 GeV, and µ i ∈ (µ 1 , µ Np ) for all other i. We then construct a c 2 s curve as a piecewiselinear function connecting these points; that is, for each i = 1, . . . , N p − 1, and for µ ∈ [µ i , µ i+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 arbitrarily closely mimic discontinuous first order transitions.
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.
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 (see [73] for the definition).
The parameter values are taken from uniform distributions µ i ∈ [µ CET , 2.6 GeV], c 2 s i ∈ (0, 1), X pQCD ∈ [1,4], in addition to which we choose roughly the same number of the "hard" or "soft" nuclear EoSs of [43]. We in addition 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 in our paper, which consists of approximately 570.000 individual EoSs. Roughly 160.000 EoSs fulfill the astrophysical constraints described in the main text of which approx. 70.000 EoSs contain at least one first order phase transition, defined as in Sec. IV above.
Appendix B: Comparison of the different interpolations Fig. 8 contains EoS bands corresponding to each of the three interpolation methods studied, with the astrophysical constraints of Sec. III implemented. To make the EoS families comparable to each other, we not only make sure that the ensembles are of roughly similar size but in addition choose the numbers of free parameters in the EoSs approximately equal. Following the approach of [16], we allow up to 4 independent segments in the piecewise polytropic interpolation (amounting to 5 free parameters), while for the spectral interpolation proposed by Lindblom [46], we use Chebyshev polynomials of degree 5 (4 free parameters). Finally, for the speed of sound interpolation, we use up to 5 independent segments (6 free parameters) in this comparison. In each case, we randomly generate large ensembles of interpolation functions, ensure that the resulting EoSs are causal and thermodynamically consistent, and in the end discard those EoSs that are in disagreement with the observational constraints discussed in Section III. Again, we add no explicit first order transitions to the EoSs, but allow continuous transitions that are arbitrarily strong, thus closely mimicking discontinuous phase transitions.
Our conclusion from the EoS plot is that the speedof-sound and polytropic interpolations produce nearly identical results, while the spectral interpolation leads to a somewhat more constrained band. The latter fact is interesting to contrast with our findings concerning the smoothing of the speed-of-sound EoS family in Fig. 3, where the shape of the EoS band obtained for large values of ∆ln can be seen to be very close to the outcome of the spectral interpolation. This fact should not come as a surprise, considering that the spectral method does not build on piecewise-defined interpolating functions, so that the resulting EoSs are smooth by construction. Figure 9 displays the families of MR curves obtained by integrating the Tolman-Oppenheimer-Volkov (TOV) equations using the above three ensembles of randomly generated EoSs. Our first observation is that, by and large, the three EoS families lead to similar results. The minimal and maximal radii for a fixed mass agree well between the different interpolations, with the spectral interpolation occupying a slightly more restricted area for stars with M < 1.2M .
The agreement between different interpolations also persists as a function of tidal deformability. Constraining the tidal deformability of a 1.4M NS according to 70 < Λ(1.4M ) < 580 [18], we see that the different interpolations still give similar maximal radii as functions of the NS mass as long as M 1.4M . In particular, the maximal radii at M = 1.4 are in excellent quantitative agreement between the different interpolation methods, as was to be expected from the previously observed tight correlation between NS radii and tidal deformabilities [16]. Considering stars with smaller masses, we observe that the speed-of-sound and polytropic interpolations allow EoSs that are extremely hard at low densities, leading to large radii R ≈ 14 km for M ≈ M , but rapidly soften at larger densities, such that for M = 1.4M the radii are smaller and consistent with the upper limits for the tidal deformability. 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 s interpolation, with the cyan color standing for EoSs that fail to reproduce 2M stars, green for those that give 70 < Λ(1.4M ) < 580, and red for the excluded region. The black dotted lines finally correspond to the 4-trope polytropes, and the orange dashed lines to the spectral interpolation method. 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. Of these three EoSs, I and II utilize the "hard" EoS of [43] until n CET , while III follows the "soft" one. We note that these three EoSs also respect the multimessenger bounds of Sec. III.
In the auxiliary material, we give each of these EoSs in a tabulated format, where the baryon number density, pressure, energy density, and speed of sound squared are all reproduced for a variety of densities from n = 0 to well above the central densities reached in maximal-mass  stars. Should the reader be interested in obtaining further example EoSs, we will be happy to provide them upon request. Finally, the corresponding EoS and MR curves are shown in Figs. 10 and 11 with our entire re-spective ensembles in the background, and the internal structures of maximal mass stars built with these EoSs in Figs. 12 and 13.