A general, scale-independent description of the sound speed in neutron stars

Using more than a million randomly generated equations of state that satisfy theoretical and observational constraints we construct a novel, scale-independent description of the sound speed in neutron stars where the latter is expressed in a unit-cube spanning the normalised radius, $r/R$, and the mass normalized to the maximum one, $M/M_{\rm TOV}$. From this generic representation, a number of interesting and surprising results can be deduced. In particular, we find that light (heavy) stars have stiff (soft) cores and soft (stiff) outer layers, respectively, or that the maximum of the sound speed is located at the center of light stars but moves to the outer layers for stars with $M/M_{\rm TOV}\gtrsim0.7$, reaching a constant value of $c_s^2=1/2$ as $M\to M_{\rm TOV}$. We also show that the sound speed decreases below the conformal limit $c_s^2=1/3$ at the center of stars with $M=M_{\rm TOV}$. Finally, we construct an analytic expression that accurately describes the radial dependence of the sound speed as a function of the neutron-star mass, thus providing an estimate of the maximum sound speed expected in a neutron star.


INTRODUCTION
The extreme conditions in the neutron-star interiors pose a formidable problem for the theoretic modelling of nuclear matter several times denser than the saturation density of atomic nuclei n s := 0.16 fm −3 . While effective field-theory calculations are arguably the most important tool to obtain theoretical predictions for the behavior of dense matter, the associated uncertainties become large at densities several times n s , such as those present in neutron-star cores. In addition, first-principle perturbative Quantum-Chromodynamics (QCD) calculations are only reliable at densities much larger than those realized in the neutron-star interior, but provide important consistency conditions for the modelling of matter at lower densities (Fraga et al. 2014;Annala et al. 2020;Annala et al. 2022;Gorda et al. 2022;Somasundaram et al. 2022). Our theoretical control of even basic quantities such as the equation of state (EOS) of dense nuclear and quark matter that, in the simplest case, is a relation between pressure and the energy density p(e), is therefore still very limited, often forcing the use of agnostic approaches to build the EOS of nuclear matter at neutron-star densities. At the same time, recent and upcoming observations of neutron stars and their merger events represent an unique opportunity to gain information about strongly coupled dense matter under conditions that are either difficult or impossible to create in experiments.
It is therefore important to combine our current theoretical knowledge with the available observational data to make predictions for the EOS and related quantities that determine the macroscopic properties of neutron stars. One such quantity is the (adiabatic) sound speed, which is defined as the derivative of the pressure with respect to the energy density at fixed entropy per baryon s (Rezzolla & Zanotti 2013) Among its many important properties, the sound speed provides a measure of the stiffness of matter, or in other words, it determines the amount of material pressure available to balance the gravitational pull and therefore prevent a neutron star from collapsing to a black hole. Clearly, a large speed of sound corresponds to a large stiffness, which, in turn, allows for the support of neutron stars with large radii R and large maximum masses M TOV . A similarly simple logic would suggest that the sound speed should reach its maximum value in the core of the star. There exist a number of works (Moustakidis et al. 2017;Tews et al. 2018;Margaritis et al. 2020;Kanakis-Pegios et al. 2021;Hippert et al. 2021;Altiparmak et al. 2022) addressing the question whether the sound speed as a function of density in QCD has an upper limit that is smaller than the speed of light. One natural conjecture for this bound is the value in conformally symmetric matter c 2 s = 1/3 (see, e.g., Bedaque & Steiner 2015; Alsing et al. 2018)), such as realized in QCD at asymptotically large density. However, it turns out assuming this conformal limit (c 2 s ≤ 1/3) strictly at all densities leads to a strong tension with astrophysical measurements of neutron-star masses M 2 M (Antoniadis et al. 2013;Cromartie et al. 2020;Fonseca et al. 2021), which favor stiff EOSs with c 2 s 1/3 at densities n s . In addition, various theoretic approaches, such as QCD at large isospin density (Carignano et al. 2017), two-color QCD (Hands et al. 2006), quarkyonic matter (McLerran & Reddy 2019;Margueron et al. 2021;Duarte et al. 2021), models for highdensity QCD (Pal et al. 2022;Ma & Rho 2021;Braun & Schallmo 2022;Leonhardt et al. 2020) and models based on the gauge/gravity duality (Ecker et al. 2017;Demircik et al. 2021;Kovensky et al. 2022) predict c 2 s > 1/3 at finite densities.
In other words, there now exists evidence and consensus that the sound speed at neutron-star densities exceeds the conformal limit. However, it still remains unclear what should this maximum value be and, more importantly, where in the neutron-star interior it is achieved. For instance, it is natural to expect that the maximum sound speed should always take place at the center of the star, as this is where the largest densities are achieved. As we will show below, this simple logic is valid for light stars but fails spectacularly as one considers stars near the maximum mass.
In this Letter we investigate the behavior of the sound speed in the neutron-star interior using input for the EOS from nuclear theory and perturbative QCD and impose observational constraints on neutron-star masses, radii and tidal deformabilities. Our goal is to make statements that are universal in the sense that they do not depend on a particular choice of the EOS or on any of the macroscopic scales such as the radius and maximum mass of a given EOS. To achieve such a model independence, we randomly generate more than a million EOSs that are by construction consistent with nuclear theory and perturbative QCD at low and high densities, respectively, and satisfy observational constraints of isolated and binary neutron-star merger measurements. We then generate probability density functions (PDF) for the sound speed inside non-rotating neutron stars and extract their median and 95% confidence intervals. By choosing dimensionless coordinates for the radial (r/R) and mass dependence (M/M TOV ), we obtain a novel and entirely scaleindependent description of the sound speed in the neutronstar interior.

METHODS
Our setup is similar to the one presented by Altiparmak et al. (2022) and we briefly review it here. The starting point is the construction of a large number of EOSs, which we achieve by patching together several different components. At the lowest densities, i.e., for n < 0.5 n s , we use the Baym-Pethick-Sutherland (BPS) model (Baym et al. 1971) to describe the neutron-star crust. In the range 0.5 n s < n < 1.1 n s we randomly sample polytropes to cover the entire range between the soft and stiff EOSs of Hebeler et al. (2013). At large densities, 40 n s , corresponding to a fixed baryon chemical potential of µ = 2.6 GeV, we impose the perturbative QCD result of Fraga et al. (2014) for the pressure p(X, µ) of cold quark matter in terms of a renormalization scale parameter X, which we sample uniformly in the range [1,4]. Although such high densities are not realized in neutron stars (see, e.g., Altiparmak et al. 2022), imposing p(X, µ) constrains the EOS at neutron-star densities. Finally, in the intermediate regime of densities, i.e., 1.1 n s < n 40 n s , we use the parametrization method introduced by Annala et al. (2020), which uses a continuous combination of piecewise-linear segments for the sound speed as a function of the chemical potential c 2 s (µ) as a starting point to construct thermodynamic quantities where µ i and c 2 s,i are parameters of the i-th segment in the range µ i ≤ µ ≤ µ i+1 (throughout this article we use units in which the speed of light and Newton's constant are equal to one, c = G = 1). The number density can then be expressed as where n 1 = 1.1 n s and µ 1 = µ(n 1 ) is fixed by the corresponding chemical potential of the polytrope. Integrating Eq.
(3) then gives the pressure where the integration constant p 1 is fixed by the pressure of the polytrope at n = n 1 . We integrate Eq. (4) numerically, using a fixed number of 7 segments of the form (2) for the sound speed. Following the procedure above, we construct a large number of EOSs by choosing randomly the maximum speed of sound c 2 s,max ∈ [0, 1] and uniformly sample the remaining free coefficients µ i ∈ [µ 1 , µ N +1 ], where µ N +1 = 2.6 GeV, and c 2 s,i ∈ [0, c 2 s,max ] in their respective domains. In this way, we construct 8 × 10 6 EOSs that are consistent with theoretical uncertainties in nuclear theory and perturbative QCD.
In order to impose constraints provided by the astronomical observations, we solve for each EOS the Tolman-Oppenheimer-Volkoff (TOV) equations and keep only those EOSs that are consistent with the mass measurements of J0348+0432 (Antoniadis et al. 2013) (M = 2.01 ± 0.04 M ) and of J0740+6620 (Cromartie et al. 2020;Fonseca et al. 2021) (M = 2.08 ± 0.07 M ) by rejecting those with maximum mass M TOV < 2.0 M . In addition, we impose the radius measurements by NICER of J0740+6620 (Miller et al. 2021;Riley et al. 2021) and of J0030+0451 (Riley et al. 2019;Miller et al. 2019) by rejecting EOSs with R < 10.75 km at M = 2.0 M and R < 10.8 km at M = 1.1 M , respectively. Finally, we exploit the detection of GW171817 by LIGO/Virgo to set an upper bound on the binary tidal deformabilityΛ < 720 (low-spin priors) (The LIGO Scientific Collaboration et al. 2019). Denoting respectively with M i , R i , and Λ i the masses, radii, and tidal deformabilities of the binary components, where Λ i = 2 3 k 2 (R i /M i ) 5 , i = 1, 2, and k 2 is the second tidal Love number, we compute the binary tidal deformability as Λ := 16 13 (5) For any choice of M 1,2 and R 1,2 , we then reject those EOSs withΛ > 720 for a chirp mass M chirp := (M 1 M 2 ) 3/5 (M 1 + M 2 ) −1/5 = 1.186 M and q := M 2 /M 1 > 0.73 as required for consistency with LIGO/Virgo data for GW170817 (The LIGO Scientific Collaboration et al. 2019). Imposing the observational constraints reduces our original set of 8 × 10 6 EOSs to ≈ 10 6 viable EOSs that form the basis for the results presented in the next section.
We note that, in principle, there exist also estimates for the upper bound on the maximum mass M TOV 2.16 +0.17  2021)). However, since this bound requires a number of uncertain model assumptions about the kilonova modeling of GRB170817A emitted by the merger event GW170817, we do not impose them on the results presented in the main text, but rather study their impact in the Appendix A.
3. RESULTS Figure 1 represents the synthesis and the essence of our novel scale-independent representation of the sound speed in neutron stars, which is described in a unit cube having as coordinates: the normalised radius r/R ∈ [0, 1], the normalised mass M/M TOV ∈ [0, 1], and the (normalised) sound speed squared c 2 s ∈ [0, 1] (to aid the visualisation of the data we restrict the unit cube to the most interesting region). In this unit cube, we report with the blue surface the median of the sound-speed-squared PDF, while the cyan and purple surfaces mark the 95% credibility intervals.
We first discuss the most prominent global features of the median sound speed. Obviously, close to the surface of the stars (r/R ≈ 1) the sound speed is small c 2 s ∼ O(10 −2 ) and approximately independent of the mass. This is because the underlying nuclear theory description at low density is tightly constrained and has small sound speed. Moving inside the star (r/R 0.8), the sound speed develops a nontrivial mass dependence. For M/M TOV 0.7, the sound speed changes from a monotonic to a non-monotonic function of r/R that has a single local maximum, c 2 s,max , as shown by the red dashed and solid lines in Fig. 1. Importantly, the radial location of the maximum sound speed depends on the mass (red dashed line), so that it is at the center of the stars (r/R = 0) for light stars (M/M TOV 0.7) but then moves to the outer layers of the stars (r/R 0.5 − 0.7) for heavy stars (M/M TOV 0.7). This interesting and somewhat surprising behaviour highlights how the structure of a compact star depends sensitively on its mass. As we will further discuss below, the value of this local maximum becomes independent of the mass for sufficiently large masses (M/M TOV 0.7) and attains a constant value of c 2 s,max ≈ 1/2 (see also Fig. 3). Another important feature of the sound speed is the non-monotonic behavior of the value at the center of the stars, c 2 s,center , as a function of the mass, as shown by the green solid line in Fig. 1 In Fig. 2 we show the PDFs for three different values of the dimensionless mass parameter (top panels) and for three values of the mass when expressed in solar masses (bottom panels). In all cases, the red lines represent the median of the PDFs, while the black solid lines the corresponding 95% credibility intervals and black dashed lines mark the conformal limit c 2 s = 1/3. Note how the red curves clearly show the transition from a monotonic to a non-monotonic radial dependence of the sound speed when going from light (left) to heavy stars (right). The sound speed in light stars increases relatively slow from the surface towards the center, reaching the conformal value only roughly at half its outer radius, i.e., r/R ≈ 0.5.
As a result, the outer layers of light stars are relatively soft and are therefore more sensitive to tidal disruptions than heavier stars. In stars with intermediate mass (M/M TOV ≈ 0.75) the sound speed increases more rapidly towards the center and remains at a constant value c 2 s ≈ c 2 s,max = 1/2 over a large portion (r/R 0.4) of the core region. Hence, these stars have a relatively large and stiff core compared to lighter ones. Finally, the heaviest stars (M/M TOV ≈ 1) do not have a sound speed that changes monotonically within the star, but, as mentioned above, develop a local maximum rather far from the center (r/R ≈ 0.65). Hence, heavy stars have relatively soft cores, but very stiff outer layers. The stiffening in the outer layers clearly is a consequence of imposing the two solar-mass constraint M TOV > 2 M , while the softening in the core is required to satisfy the perturbative QCD boundary conditions at large densities and, at the same time, maintain causality at all densities. A physically intuitive way of looking at this behaviour is the following: as the stellar core softens, the burden of keeping the heavy star stable against gravitational collapse has to be taken by the outer layers, that have therefore to become stiff and with large sound speeds. Hence, the appearance of a local maximum in the sound speed is a direct consequence of the inter- play between the astrophysical constraints (M TOV 2 M ) and the perturbative QCD constraints at high densities.
We note that traditional nuclear-theory EOSs that only capture confined nuclear matter typically do not take into account the perturbative QCD constraints at large densities and therefore cannot lead to the soft core that our statistical approach points out. Indeed, widely used EOSs built on the present understanding of low-density nuclear theory typically predict high-mass stars with cores that are systematically stiffer and therefore with larger sounds speed (see Fig. 2 in Appendix B).
The panels in Fig. 2 also make it very easy to appreciate that while light and intermediate-mass stars have superconformal central sound speeds (c 2 s,center > 1/3; black dashed lines), heavy stars have sub-conformal central sound speeds (c 2 s,center < 1/3), reaching a local minimum of c 2 s,center ≈ 0.3 for M/M TOV = 1.
The smooth behaviour of the medians of the sound speed shown in the panels of Fig. 2 encourages the representation of this behaviour in terms of analytic functions, namely where x := r/R and α = 0.18, β = 0.14, γ = 0.71, δ = −6.30, = −1.5, ζ = 10.00, η = 8.9 are the fitting parameters for the most interesting case, namely the median of the PDF when M/M TOV = 1. For completeness, Table 1 in Appendix C reports the values of these coefficients for different mass cuts and different credibility intervals. A few remarks on expression (6) are in order. First, while these functions are normalised in terms of the stellar radius and maximum mass -none of which are known at present - they provide useful information nevertheless. For instance, notwithstanding the limited knowledge available now, it is already possible to conclude that a neutron star with a mass close to the maximum mass of ∼ 2.2 M (see, e.g., Rezzolla et al. 2018) and an average radius of ∼ 12 km (see, e.g., Altiparmak et al. 2022) will have a maximum sound speed of c s ∼ 0.7 and at ∼ 8 km from the center. This valuable information can already be included in nuclear-theory calculations of new EOSs. Second, although a fit with seven parameters may seem excessive, they are necessary to faithfully represent the key features of the distribution, such as their values at the center (r = 0) and close to the surface of the stars (r = R), as well as the location and value of the local maximum at intermediate radii. Finally, the number of parameters could be reduced by imposing analytic conditions on the approximate behavior of the PDFs in certain limits (e.g., using the vanishing slope at r/R = 0, 1). However, we have preferred to be conservative and not to impose such constraints at the cost of a larger number of parameters. Figure 3 shows the behaviour of the median of the local maximum sound speed c 2 s,max (red line) and of the sound speed at the neutron-star center c 2 s,center (green line) as a function of the (normalized) stellar mass. The red-and greenshaded areas indicate the corresponding 95% credibility intervals, while the black dashed and dotted lines mark the values 1/2 and 1/3, respectively. In essence, Fig. 3 highlights that the maximal sound speed in stars with M/M TOV 0.7 appears in their center, because the red and green curves coincide for M/M TOV 0.7. On the other hand, for stars with M/M TOV 0.7 the two curves split: the sound speed at the stellar center decreases monotonically towards the stellar interior reaching a minimum of c 2 s,center ≈ 0.3 for the maximum-mass stars (M/M TOV = 1). Conversely, the maximum sound remains at a constant value c 2 s,max ≈ 1/2 for all masses M/M TOV 0.7. It is suggestive that this constant value is so close to half of the speed of light, but it is difficult to invoke any first-principle argument on why this should be the case.
Finally, in Fig. 4 we report the behaviour of the median of the sound speed as a function of the (normalized) stellar mass at different positions in the star, as indicated by the colorbar (these are essentially cuts at different values of r/R of the blue median surface in Fig. 1). In particular, the dark-green line represents a cut of the median at the stellar center, while the dark-blue line shows the behaviour at the stellar surface. In essence, Fig. 4 highlights how the sound speed changes when going from the surface (dark-blue line) towards the core (dark-green line) in stars as a function of their mass. At r/R ≈ 0.65 (light-blue line) the sound speed reaches in the heaviest stars its maximum value c 2 s 1/2, indicated by the black dashed line.
This shows that in the outer layers (r/R 0.65) of every star, regardless of its mass, the sound speed is a monotonically increasing function of the mass and has values c 2 s < 1/2. On the other hand, the sound speed in the core region (r/R ≈ 0), shown in dark green, is a non-monotonic function of the neutron-star mass. As a result, the maximum sound speed at the stellar center, i.e., c 2 s,center ≈ 1/2, is not attained in the lightest (M/M TOV = 0.5), nor in the heaviest stars (M/M TOV = 1), but rather at intermediate mass M/M TOV ≈ 0.75. Furthermore, the sound speed at the center of light stars (c 2 s,center ≈ 0.4) is actually larger than the corresponding value of the heaviest stars (c 2 s,center ≈ 0.3) that, as discussed above, experience a significant softening.

SUMMARY AND CONCLUSIONS
We have studied the sound speed distribution inside neutron stars using a large set of randomly generated EOSs that are consistent with nuclear theory and perturbative QCD results in their respective ranges of validity and are in agreement with astrophysical pulsar observations and gravitational wave detections from binary neutron-star mergers. Our main result is a novel and a scale-independent representation of the sound speed in a unit-cube spanning the normalised radius r/R and the mass normalized to the maximum one, M/M TOV .
This innovative way of thinking about the sound speed has allowed us to draw a number of general conclusions that increase our insight into the quantitative and qualitative behavior in neutrons stars of the sound speed in particular, and of dense nuclear matter more in general. More specifically, we find that: (i) close to the surface of the stars, the sound speed is small and approximately independent of the mass. However, moving inside the stars, the sound speed develops a non-trivial mass dependence, and for M/M TOV 0.7, it changes from a monotonic to a non-monotonic function of position r/R with a single local maximum, c 2 s,max . The radial location of the maximum sound speed depends on the mass and it is at the center of the star for light stars but then moves to the outer layers of the star for heavy stars.
(ii) the value of this local maximum becomes independent of the mass for sufficiently large masses (M/M TOV 0.7) and attains a constant value of c 2 s,max ≈ 1/2. (iii) the sound speed at the center of the stars, c 2 s,center , also exhibits a non-monotonic behavior as a function of the mass, with the maximum value of c 2 s,center being reached by stars with intermediate mass, while the minimum value is not obtained in the lightest stars but in the heaviest stars where it is c 2 s,center 0.3, thus below the conformal limit.
(iv) using the sound speed as a measure for stiffness, we find that light stars are soft in the outer layers (r/R 0.5 − 0.7) and stiff in the core, while heavy stars have a soft core and stiff outer layers. This is because the sound speed increases only slowly towards the center in light stars, but rapidly in heavy stars where it approaches a local minimum that is smaller than the conformal limit in the core.
(v) finally, we provide a simple fitting formula for the median sound speed and its confidence interval in the neutronstar interior. This information can already be included in nuclear-theory calculations of modern EOSs to constrain the behaviour of the sound speed in those regions where nucleartheory predictions have large uncertainties.
In summary, the non-trivial behaviour of the sound speed as function of the radial position inside the stars can be seen as a probe to identify changes of the matter composition. In particular, the softening of heavy neutron-star cores points to the appearance of new degrees of freedom such as hyperons or deconfined quarks.
There are a number of possible extensions to our work. First, while our method of constructing EOSs includes, in principle, cases that closely resemble a first-order phase transition, they only represent a negligible subset of our ensemble as we do not explicitly enforce them and hence their statistical weight is rather small. It would be therefore interesting to include models that, by construction, include a first-order phase transition and study their impact on our results. Second, in the current and also our previous work (Altiparmak et al. 2022), we took the so-called frequentist approach for the statistical interpretation of our results, that is, we impose the constraints with a hard-cutoff neglecting their statistical nature. An alternative approach to obtain a statistical interpretation would be a Bayesian analysis where the statistics of the observational uncertainties is taken into account. Finally, another interesting generalization of our work is its extension to study how rotation -and in particular rapid rotation -affects the properties of the sound speed in neutron stars. We plan to address many of the these issues in ongoing and future work.  Nathanail et al. 2021, for a discussion on the constraints derived from GW190814). Such a constraint is obtained using the detection of the gravitational-wave event GW170817 and the modeling of the kilonova signal of GRB170817A, together with quasi-universal relations between the maximum masses of uniformly rotating and non-rotating stars (Breu & Rezzolla 2016). Figure 1 shows the sound-speed PDF in a way that is similar to the three top panels of Fig. 2 Fig. 2, the only difference is that here an upper limit on the maximum mass MTOV < 2.16 M is imposed.
2 reveals that the overall features of the sound speed remain unchanged. The most significant difference is that the distributions are systematically narrower, simply because the stiff EOSs with large sound speeds are now penalized and rejected. This is particularly clear in the panel for heavy stars (right), where the sound speed increases less rapidly in the outer layers of the star. Interestingly, the value of the local maximum sound speed, c s,max , varies only minimally, changing from c s,max 1/2 in the absence of a maximum-mass constraint, to c s,max 0.45 when a maximum-mass constraint is imposed.

B. COMPARISON TO MICROPHYSICAL MODELS
In this Appendix we compare our sound-speed PDF with three microphysical models that have been widely used in the literature to study the properties of neutron stars and their mergers. Figure 2, in particular, compares the radial behaviour of the sound-speed PDF in maximally massive stars with the corresponding sound speeds obtained with these EOSs. The yellow line shows the sound speed of a pure nuclear-matter model, the Hempel-Schaffner-Bielich (HS) EOS with DD2 relativistic mean field interactions (Typel et al. 2010;Hempel & Schaffner-Bielich 2010). This EOS is known to be relatively stiff as can be seen from the large sound speed in the neutron-star interior. Interestingly, for r/R 0.7 the sound speed from this EOS agrees very well with the median sound speed (red line) from our approach, even though our construction does not use any input from the (HS)DD2 EOS. However, deeper inside the star, our result is very different from (HS)DD2 and gives less than half the sound speed in the neutron-star core. We should note that the (HS)DD2 EOS gives a relatively high maximum mass of M TOV ≈ 2.5 M and a tidal deformability of Λ 1.4 ≈ 690 for 1.4 M -stars, which is in tension with the upper bound Λ 1.4 580 derived from the inspiral part (low-spin prior) of GW170817 (The LIGO Scientific Collaboration et al. 2019). This means (HS)DD2 only provides a good description for dense matter n n s in light stars, but is probably too stiff to describe the properties of heavy neutron stars with cores several times denser than n s .
The green line, on the other hand, corresponds to the Banik-Hempel-Bandyopadhyay (BHBΛ) EOS (Banik et al. 2014), which is an extension of the (HS)DD2 EOS with Λ-hyperons, i.e., particles that contain strange quarks. One of the characteristic features of hyperonic degrees of freedom is that they lead to a softening of the EOS, which can be seen from the local minimum in the sound speed around r/R ≈ 0.8. The maximum mass of this EOS is relatively low M TOV ≈ 2.0 M and therefore in tension with the mass measurement M 2 M of heavy pulsars. Hence, the EOS is probably too soft to account for the properties of light stars and the outer layers of heavy stars as can be seen from the large difference between the red and green curve at r/R 0.5 in Figure 2.
Finally, the cyan line corresponds to the intermediately stiff V-QCD EOS of Demircik et al. (2021), that combines the traditional nuclear theory (HS)DD2 and the Akmal-Pandharipande-Ravenhall (APR) EOS (Akmal et al. 1998) EOSs at low densities with a string theory inspired model for QCD to describe dense baryonic and quark matter. This model gives M TOV = 2.14 M and Λ 1.4 = 511 and therefore satisfies conveniently the two-solar mass and tidal deformability constraints. Furthermore, it has been verified recently (Tootle et al. 2022) via binary neutron-star merger simulations, that this EOS is not excluded by the expected one second-long lifetime (Gill et al. 2019) of the post-merger remnant of GW170817. However, although this model passes all currently known theoretical and observational constraints, it is also not able to account for the softening of the core predicted by the red curve.
In summary, none of the three very different microphysical EOSs shown in Fig. 2 is able to describe the behavior of the sound speed predicted by our model-agnostic sampling approach. The explanation is simple and reflects the difficulties that nucleartheory calculations have in matching the high-density, perturbative QCD constraints. It is remarkable that EOSs such as those considered in Fig. 2 -that are in good agreement with all currently known observational constraints -have radial profiles of the sound speed that are quite different those predicted from our agnostic approach. This fact shows that the traditional methods of constraining EOS models only with global neutron-star properties such as the maximum mass, the radii and tidal deformabilities, might not be sufficient. Rather, additional information on the properties inside the stars is needed and can provide further nontrivial constraints. Hence, the importance of our Eq. (6) is to provide a novel constraint in terms of the radial distribution of the sound speed inside stars of a given mass. Using this constraint will ensure that the newly suggested EOSs will not only satisfy the astrophysical constraints, but are also compatible with the perturbative QCD constraints at much larger densities. We are not aware of a comparable constraint in the literature.

C. COMPREHENSIVE DESCRIPTION OF THE FITTING
In Table 1 we provide the numerical values of the fitting coefficients for the results shown in Fig. 2; these coefficients reflect the astrophysical constraint imposed here of M TOV 2 M and their values may change slightly if this constrain is increased. As briefly mentioned in the main text, some of the parameters in Eq. (6) can in principle be related by strictly imposing the approximate behavior of the distribution in certain limits. For example, the numeric results have approximately vanishing slope at r/R = 0, 1 and the value of the local maximum of the median at 0 < r/R < 1 turns out to be very close to the value 1/2 in stars with M/M TOV > 0.7. Imposing these constraints would allow to express three of the coefficients in terms of the remaining four, but at the same time make the formal expression (6) more involved. In addition, the upper and lower bounds of the credibility interval do not show the saturation behavior of the local maximum seen in the median such that it would not be possible anymore to expressed them by the same fitting law. We also remark that the sound speed loses the local maximum for M/M TOV 0.7 and becomes monotonic, allowing for a simpler fitting law with fewer parameters. This can be directly seen from the values of α provided in Table 1, which become small compared to the other parameters, meaning that the first exponential function in Eq. (6) can be neglected, reducing the number of relevant parameters to five. Table 1. Numeric values of the fitting parameters in Eq. (6) for the median sound speed and the lower and upper bound of its 95% confidence interval for three different values of M/MTOV used in Fig. 2.