# A three-parameter characterization of neutron stars’ mass$-$radius relation
and equation of state

Dmitry D. Ofengeim1 Corresponding Author<EMAIL_ADDRESS>Peter S.
Shternin2 Tsvi Piran1 1Racah Institute of Physics, The Hebrew University,
Jerusalem 91904, Israel  2Ioffe Institute, Politekhnicheskaya 26, St.
Petersburg, 194021, Russia

###### Abstract

Numerous models of neutron star (NS) equation of state (EoS) exist based on
different superdense-matter physics approaches. Nevertheless, some NS
properties show universal (EoS-independent) relations. Here, we propose a
novel class of such universalities. Despite different physics inputs, a wide
class of realistic nucleonic, hyperonic, and hybrid EoS models can be
accurately described using only three parameters. For a given EoS, these are
the mass and radius of the maximum-mass NS (or pressure and density in its
center) and the radius of a half-maximum-mass star. With such a
parametrization, we build universal analytic expressions for mass-radius and
pressure-density relations. They form a semi-analytic mapping from the mass-
radius relation to the EoS in NS cores (the so-called inverse Oppenheimer-
Volkoff mapping). This mapping simplifies the process of inferring the EoS
from observations of NS masses and radii. Applying it to current NS
observations we set new limits on the high end of the EoS.

##  I Introduction

One of the most intriguing problems of modern physics is to unveil the mystery
of the neutron star (NS) equation of state (EoS, pressure $P$ $-$ density
$\rho$ relation) [1]. One of the few ways to study it is to measure NS masses
$M$ and radii $R$. As the star is in hydrostatic equilibrium, its EoS
determines the $M-R$ relation. Inverting this relation, we can find the EoS.
Non-rotating NS hydrostatics is described by the Tolman-Oppenheimer-Volkoff
(TOV) equations [2, 3]. The Oppenheimer-Volkoff (OV) mapping from the $P-\rho$
to $M-R$ is a bijection, i.e., there exists inverse OV mapping (IOVM) [4].
Thus unambiguous finding of EoS from $M$ and $R$ observations is, in
principle, possible. Additionally, a maximum-mass NS (MMNS) exists regardless
of the EoS model used. The MMNS characteristics, such as mass
$M_{\text{{tov}}}$, radius $R_{\text{{tov}}}$, pressure $P_{\text{{tov}}}$,
density $\rho_{\text{{tov}}}$ in the center, etc., are specific for a given
EoS.

While the true NS EoS is unique, there are hundreds of various theoretical EoS
models present in the literature [5]. Each yields its own $M-R$, $M$ — central
$\rho$ and $P$ relations and, in particular, its own $M_{\text{{tov}}}$,
$R_{\text{{tov}}}$, $\rho_{\text{{tov}}}$, and $P_{\text{{tov}}}$.
Nevertheless, it has been known for a long time that some relations between
various NS properties are universal across the EoS manifold (e.g. [6, 7, 8, 9,
10, 11, 12, 13, 14]). Some of these relations, like the I-Love-Q relation [15,
16], originate from a unifying ability of relativistic gravity [9] and are
very accurate for almost all existing models of EoS. Others are less precise
and exist due to some common features of a subclass of EoSs, which are
considered to be “realistic” (see, e.g., [17]).

In this work, we present a new class of universal relations describing the
fundamental NS features—their $M-R$ and $P-\rho$ relations (Sec. II). They are
based on the idea that $M_{\text{{tov}}}$, $R_{\text{{tov}}}$,
$P_{\text{{tov}}}$, and $\rho_{\text{{tov}}}$ set suitable universal scales
for masses, radii, pressures, and densities of NSs. Using these relations we
build an explicit semi-analytic IOVM and apply it to directly infer the high-
density EoS from NS mass-radius observations (Sec. III). The accuracy of the
obtained constraints is competitive to the results of other works (Sec. IV).

##  II Universal relations

The universal relations presented below were verified for a collection of
$169$ EoSs based on various underlying physics-motivated models of low-energy
hardronic interactions. Among those, $111$ pass the criterion
$M_{\text{{tov}}}\gtrsim 2\,$$M_{\odot}$ (based on astrophysical measurements
of NS masses [18, 19, 20]; see below). For $17$ models the maximal speed of
sound is higher than the speed of light $c$. However, we keep all these models
to obtain a more comprehensive cover of the $M-R$ and $P-\rho$ planes.
Removing these EoS models doesn’t change significantly the fits we present
below. In this collection, $100$ models exhibit nucleonic composition, $36$
have additional heavier baryons (hyperons and rarely $\Delta$-isobars), and
$33$ have hybrid composition, in which a transition to deconfined quark matter
occurs at some density. See Table S1 in Supplementary materials 111See
Supplemental Material at [URL will be inserted by publisher] which includes
Refs. [48, 49, 50, 51, 10, 52, 53, 54, 55, 24, 56, 57, 58, 59, 60, 61, 62, 63,
26, 27, 30, 31, 32, 64, 65, 66, 67, 68, 69, 38, 35, 39, 65, 70].  for details.

Figure 1: A geometrical view of the EoS manifold. The curves correspond to
different EoS models: nucleonic (blue), with hyperons and $\Delta$-isobars
(green) and hybrid (red). Except for a few outliers that correspond to
obsolete EoSs with $M_{\text{{tov}}}<2\,$$M_{\odot}$, all EoSs lie on a single
surface approximated by Eq. (1) (black grid). The histograms show
distributions of relative root-mean-square and maximal errors (the bin height
is the number of EoSs in the bin).

The $M-R$ curves can be described by the universal form:

$R/R_{\text{{tov}}}=1+\left[2\left(\sqrt{2}-1\right)r_{\text{1/2}}-a\right]\sqrt{1-M/M_{\text{{tov}}}}\\\
+\left[2\left(\sqrt{2}-1\right)r_{\text{1/2}}-2+a\sqrt{2}\right](1-M/M_{\text{{tov}}}),$
(1)

where $a=0.492$. The maximum mass $M_{\text{{tov}}}$, corresponding radius
$R_{\text{{tov}}}$, and the ratio
$r_{\text{1/2}}=R_{\text{1/2}}/R_{\text{{tov}}}$ with
$R_{\text{1/2}}=R(M_{\text{{tov}}}/2)$, are the three parameters governing the
manifold of NS $M-R$ curves. Geometrically, this means that the $M-R$
relations for different EoSs form a surface in the 3-dimensional space
$(M/M_{\text{{tov}}},R/R_{\text{{tov}}},r_{\text{1/2}})$. This is an
approximate statement (see Fig. 1). Several clear outliers correspond to
obsolete models with severe softening around phase transition points and
$M_{\text{{tov}}}<2\,$$M_{\odot}$. Moreover, this approximation is valid only
for $M\gtrsim 1$$M_{\odot}$. Below this mass, the form of $M-R$ curves is
controlled by a larger number of parameters. For the whole set of models,
errors of this approximation rarely exceed $3\%$. For the outliers, the
maximal error is $\sim 10\%$.

Figure 2: EoS-independent correlations of MMNS characteristics. Each symbol
correspond to a single EoS (nucleonic, with hyperons and isobars, hybrid). The
solid lines show the approximations (2).
$\mathcal{R}=\sqrt{P_{\text{{tov}}}/(G\rho_{\text{{tov}}}^{2})}$,
$\mathcal{M}=\rho_{\text{{tov}}}\mathcal{R}^{3}$, $\mathcal{C}=R_{\rm
TOV}^{\text{fit}}\sqrt{G\rho_{\text{{tov}}}}$, where $R_{\rm
TOV}^{\text{fit}}$ is given by Eq. (2b). The exponents $a_{i}$ and $b_{i}$,
with $i=M,R,c$, are given in Table 1. The bottom panels show relative errors
of the fits.

$M_{\text{{tov}}}$ and $R_{\text{{tov}}}$ are strongly correlated with other
parameters describing MMNS, like $\rho_{\text{{tov}}}$, $P_{\text{{tov}}}$,
and the speed of sound in the center $c_{\text{s{tov}}}$ [10, 14, 22]. Here we
present updated fits to these correlations:

$\displaystyle M_{\text{{tov}}}$
$\displaystyle=\sqrt{\frac{P_{\text{{tov}}}^{3}}{G^{3}\rho_{\text{{tov}}}^{4}}}\frac{1}{f_{M}(P_{\text{{tov}}},\rho_{\text{{tov}}})},$
(2a) $\displaystyle R_{\text{{tov}}}$
$\displaystyle=\sqrt{\frac{P_{\text{{tov}}}}{G\rho_{\text{{tov}}}^{2}}}\frac{1}{f_{R}(P_{\text{{tov}}},\rho_{\text{{tov}}})},$
(2b) $\displaystyle c_{\text{s{tov}}}$
$\displaystyle=\sqrt{\frac{P_{\text{{tov}}}}{\rho_{\text{{tov}}}}}\frac{f_{c}(P_{\text{{tov}}},\rho_{\text{{tov}}})}{f_{R}(P_{\text{{tov}}},\rho_{\text{{tov}}})}=c\
\zeta(P_{\text{{tov}}},\rho_{\text{{tov}}})\ .$ (2c) The dimensionless
functions $f_{M}$, $f_{R}$, and $f_{c}$ have the generic form
$f_{i}=c_{i}\left(\frac{P_{\text{{tov}}}}{\rho_{0}c^{2}}\right)^{a_{i}}\left(\frac{\rho_{\text{{tov}}}}{\rho_{0}}\right)^{b_{i}}+d_{i}\
.$ (2d)

The coefficients $a_{i}$, $b_{i}$, $c_{i}$, $d_{i}$ are listed in Table 1, and
$\rho_{0}=2.8\times 10^{14}\,\text{g}\,\text{cm}^{-3}$ is the nuclear
saturation density. These correlations are illustrated in Fig. 2.

Table 1:  Coefficients in the fits (2) of $M_{\text{{tov}}}-P_{\text{{tov}}},\rho_{\text{{tov}}}$, $R_{\text{{tov}}}-P_{\text{{tov}}},\rho_{\text{{tov}}}$, and $c_{\text{s{tov}}}-P_{\text{{tov}}},\rho_{\text{{tov}}}$ correlations. The last two columns show the relative root mean square and relative maximum errors of the fits. $f_{i}$ | $a_{i}$ | $b_{i}$ | $c_{i}$ | $d_{i}$ | rrms | max error
---|---|---|---|---|---|---
$f_{M}$ | $1.41$ | $-1.39$ | $5.86$ | $0.270$ | 0.013 | 0.11
$f_{R}$ | $0.744$ | $-0.839$ | $2.45$ | $0.357$ | 0.021 | 0.12
$f_{c}$ | $0.985$ | $-0.990$ | $3.02$ | $0.529$ | 0.048 | 0.22

Since the OV mapping is a bijection [4], the manifold of EoS models is also
3-parametric, at least in the density range occupying most of the volume of NS
with $M\gtrsim 1$$M_{\odot}$. Correlations (2) provide a one-to-one mapping
between the $(M_{\text{{tov}}},R_{\text{{tov}}},r_{\text{1/2}})$ and
$(\rho_{\text{{tov}}},P_{\text{{tov}}},r_{\text{1/2}})$ parameter spaces. The
latter one is suitable for constructing the universal $P-\rho$ approximation.
It is convenient to split it into a high and low density parts.

For $\rho>3\rho_{0}$, there exist a $r_{\text{1/2}}$-independent fit
$P=P_{\text{(high)}}(\rho;P_{\text{{tov}}},\rho_{\text{{tov}}})$:

$\frac{P_{\text{(high)}}}{P_{\text{{tov}}}}=\left(\frac{\rho}{\rho_{\text{{tov}}}}\right)^{\gamma(P_{\text{{tov}}},\rho_{\text{{tov}}})-a_{0}}\bigg{\\{}1+a_{0}\biggl{(}1\\\
\left.-\frac{\rho}{\rho_{\text{{tov}}}}\biggr{)}+\left[b_{0}+b_{1}\zeta(P_{\text{{tov}}},\rho_{\text{{tov}}})\right]\left(1-\frac{\rho}{\rho_{\text{{tov}}}}\right)^{p}\right\\}^{-1}\
,$ (3)

where $\zeta(P_{\text{{tov}}},\rho_{\text{{tov}}})$ is defined in Eq. (2c),
$\gamma(P_{\text{{tov}}},\rho_{\text{{tov}}})=\zeta^{2}\rho_{\text{{tov}}}c^{2}/P_{\text{{tov}}}$,
and $a_{0}=-0.392$, $b_{0}=-0.869$, $b_{1}=1.69$, $p=2.98$. Root-mean-square
errors of this approximation rarely exceed $10\%$. 222We presented a similar
fit in [22], but here we reduce the number of fitting coefficients.

Figure 3: Demonstration of the approximation (4) of $P(\rho)$ dependencies for
selected nucleonic (BSk24), hyperonic (NL3$\omega\rho$), and hybrid (QHC19-A)
EoSs. Symbols stand for exact values for given models, lines show the
approximations with $\rho_{\text{{tov}}}$, $P_{\text{{tov}}}$, and
$r_{\text{1/2}}$ taken from the corresponding EoSs. Crosses mark the
$(\rho_{\text{{tov}}},P_{\text{{tov}}})$ points. The inset shows the
distributions of relative root-mean-square and maximal errors of the fit for
the whole EoS collection (the bin height is the number of EoSs in the bin).

For $\rho<3\rho_{0}$ the dependence on $r_{1/2}$ cannot be ignored and we
built an approximation with three parameters, $\rho_{\text{{tov}}}$,
$P_{\text{{tov}}}$, and $r_{\text{1/2}}$, valid for
$\rho_{0}<\rho<\rho_{\text{{tov}}}$:

$P=P_{\text{(high)}}\exp\left\\{\left(u_{0}+u_{1}r_{\text{1/2}}+u_{2}r_{\text{1/2}}^{2}\right)\right.\\\
\left.\times\left(e^{-v\rho/\rho_{0}}-e^{-v\rho_{\text{{tov}}}/\rho_{0}}\right)\right\\},$
(4)

where $u_{0}=-61.2$, $u_{1}=81.9$, $u_{2}=-27.2$, $v=1.33$, and
$P_{\text{(high)}}$ is given by Eq. (3) (see Fig. 3). On average, the rrms
error is about $5-10\%$, it is smaller for nucleonic EoSs and larger for
hyperonic and hybrid models. Although the maximum error could reach $3-5$
times the rms, we show below that the approximation (4) together with Eqs. (1)
and (2) allows one to make valuable inferences on the EoS when applied to
observational data.

##  III Inverse OV mapping & constraints on EoS

Using these approximations one can explicitly perform IOVM. If the $M-R$
relation is known, then $M_{\text{{tov}}}$, $R_{\text{{tov}}}$, and
$r_{\text{1/2}}$ are available. From Eqs. (2a) and (2b), one derives
$P_{\text{{tov}}}$ and $\rho_{\text{{tov}}}$. Substituting these values and
$r_{\text{1/2}}$ to Eq. (4), one obtains the $P(\rho)$ curve in the range
$\rho_{0}...\rho_{\text{{tov}}}$.

In practice, instead of the full $M-R$ curve, one has access to only a few
measurements of $M$ and $R$ for a set of observed NSs. In this case, one
should fit the observations with Eq. (1) and obtain $M_{\text{{tov}}}$,
$R_{\text{{tov}}}$, and $r_{\text{1/2}}$. Specifically, one can substitute
Eqs. (2) into Eq. (1), and use the triade $\rho_{\text{{tov}}}$,
$P_{\text{{tov}}}$, and $r_{\text{1/2}}$ as fitting parameters.

The data typically strongly constrain
$R_{\text{1/2}}=r_{\text{1/2}}R_{\text{{tov}}}$ and poorly constrain the
properties of MMNS, so we add two theoretical conditions. First, we set
$c_{\text{s{tov}}}(P_{\text{{tov}}},\rho_{\text{{tov}}})<1.05\times c$ that is
the causality limit with the safety factor $1.05$ accounting for the
uncertainty of the fit (2c). Second, we impose $r_{\text{1/2}}\geqslant 1$
since we did not find an EoS model violating this condition (see Table S1 and
Fig. S1c).

Figure 4: Constraining the $M-R$ curve from NS observations. Solid-colored
contours show $90\%$ credible regions for NSs with both $M$ and $R$ measured.
Vertical bars and corresponding colored strips display $1\sigma$ mass ranges
for NSs in binaries. The thick brown line show the upper limit on
$M_{\text{{tov}}}$ adopted from the GW 170817 electromagnetic afterglow
interpretation [24]. The grey-shaded region is prohibited due to the “maximum
compactness” criterion (see, e.g., [25]). The red-shaded region is excluded
due to the causality condition applied to the correlations (2). The cross-
hatched region is formed by the $90\%$ credible intervals for radii at given
NS masses, obtained from fitting the observational data by Eq. (1). Thin
dashed and dash-dotted 90% credible contours show similar results by Annala et
al. [26, model “$c^{2}_{s,4}$”], and by Brandes et al. [27], respectively. The
thick black contour shows the 2D-credible region for $M_{\text{{tov}}}$ and
$R_{\text{{tov}}}$ for our results. Notice a statistical paradox: due to the
almost horizontal behavior of the $M-R$ curve near the MMNS point, the most
probable radii for fixed $M\sim M_{\text{{tov}}}$ are significantly greater
than $R_{\text{{tov}}}$.  Figure 5: EoS constrained from NS observations. The
cross-hatched strip is formed by $90\%$ credible intervals for pressures at
given densities, obtained from the proposed IOVM procedure. Dotted lines show
the same but with account for fit errors treated as systematic uncertainties.
The thick black contour is the $90\%$ 2D-credible area for $P_{\text{{tov}}}$
and $\rho_{\text{{tov}}}$. The orange and green bands are similar to the
cross-hatched stip but adopted from Brandes et al. [27] and Annala et al. [26,
model “$c^{2}_{s,4}$”], correspondingly. The red-shaded region is prohibited
for $P_{\text{{tov}}}$ and $\rho_{\text{{tov}}}$ by the causality condition
applied to Eq. (2c). The “pQCD” line depicts the upper limit on the EoS which
follows from the perturbative QCD matching [28]. The “$\chi$EFT” band is the
“chiral corridor” for the low-density EoS based on the chiral effective field
theory [29].

We apply this model to current $M-R$ observations, using the Bayesian
inference (see e.g. [26, 27, 30, 31, 32, 33] and Supplementary materials) to
obtain the model parameters $P_{\text{{tov}}}$, $\rho_{\text{{tov}}}$, and
$r_{1/2}$. The observations fall into two classes. The first one consists of
those which give joint information on $M$ and $R$. We use 8 such sources
(color of corresponding contours in Fig. 4 given in brackets):

  1. (i)

the NICER objects PSR J0030+0451 [34, 35, 36, 37] (blue; data from [35]) and
PSR J0740+6620 [38, 34, 32] (orange; data from [38]).

  2. (ii)

the thermally emitting isolated NS in the Cassiopeia A supernova remnant [39]
(CasA, green);

  3. (iii)

low-mass X-ray binaries which exhibit X-ray bursts during super-Eddington
accretion, namely 4U J1702–429 (red) [40], 4U 1724–307 (turquoise) and SAX
J1810.8–260 (purple) [41];

  4. (iv)

the $M-R$ data inferred from the GW 170817 signal [42] (magenta).

The second class consists of NSs with measured masses only, which limit
$M_{\text{{tov}}}$ from below. We consider 8 sources:

  1. (i)

the radio sources PSR J1614–2230 [43] and PSR J0348+0432 [18] (light blue
horizontal strips in Fig. 4) as the most robust;

  2. (ii)

the most massive ($M\gtrsim 2$$M_{\odot}$) “spiders”—binary systems containing
a NS and a low-mass companion—PSR J0952–0607, PSR J1311–3430, PSR J1810–1744,
PSR J1653–0158 [20], and PSR J2215+5135 [44] (light green strips in Fig. 4),
which are more model-dependent;

  3. (iii)

an interpretation of the electromagnetic afterglow of GW 170817 which yields
$M_{\text{{tov}}}<2.16\pm 0.16\,$$M_{\odot}$ [24] (thick brown line in Fig.
4).

The cross-hatched strip and thick black contour in Fig. 4 show the result of
fitting these data by Eq. (1). Fig. 5 shows the limits on the EoS obtained
with Eqs. (2)—(4). The one-dimensional 68% credible intervals for the key
parameters of the approximations (1)—(4) are
$M_{\text{{tov}}}=2.28^{+0.05}_{-0.06}\,$$M_{\odot}$,
$R_{\text{{tov}}}=11.2^{+0.4}_{-0.3}\,$km,
$R_{\text{1/2}}=12.0^{+0.4}_{-0.3}\,$km,
$\rho_{\text{{tov}}}=7.7^{+0.6}_{-0.3}\,\rho_{0}$, and
$P_{\text{{tov}}}=3.5^{+1.0}_{-0.4}\,\rho_{0}c^{2}$.

The credible area for the EoS matches smoothly the low-density “EoS corridor”
derived from the chiral effective field theory ($\chi$EFT) [29]. Even though
our analysis does not include $\chi$EFTs constraints as a prior (the EoS
collection contains both $\chi$EFT-consistent and $\chi$EFT-inconsistent
models), our results are consistent with the requirement to pass the “chiral
corridor” [45].

One has to keep in mind that the universal relations (1)—(4) are approximate.
Their errors could be treated as systematic uncertainties when the IOVM
procedure is performed. We accounted for these uncertainties, bootstrapping
individual errors of the approximations. Eventually, the $M-R$ area and MMNS
properties remain unaffected, and the effect for EoS is shown by the black
dotted lines in Fig. 5. Indeed, these systematics slightly widen the EoS
credible area at low densities but does not affect the high-density part.

##  IV Discussion & Conclusions

In this work, we built universal approximations describing the $M-R$ curves
and $P(\rho)$ dependencies for a wide set of EoSs. This set includes models
with various compositions (nucleonic, hyperonic, hybrid) that are based on
different approaches to the underlying physics. The key idea is that three
parameters are sufficient to describe the EoS to a good accuracy. We show that
it is convenient to choose two of those as MMNS properties (say,
$P_{\text{{tov}}}$ and $\rho_{\text{{tov}}}$) and the third as the ratio
$r_{\text{1/2}}=R_{\text{1/2}}/R_{\text{{tov}}}$.

These approximations have numerous applications in NS studies. Here we
explored the very basic one—explicit implementation 333See other explicit
implementations in, e.g., [4, 71, 72]. of IOVM and inference of the EoS from
observations. The allowed range of EoS is consistent but narrower than those
found in recent works [26, 27]. Notice that [27, 26] use an additional
theoretical condition of matching the perturbative quantum chromodynamics
(pQCD) limit. Apparently, this limit [28] coincides numerically with our
causality limit (see Fig. 5), so our results satisfy this condition as well.

Our analytic model avoids the time consuming need of repeatedly solving TOV
equations [26, 27] when attempting to infer the EoS from observed NSs masses
and radii. The agreement of our results with those of [26, 27], which are
applicable even if the real EoS differs strongly from those that are in our
collection, indicates that the true EoS is somewhere in between the existing
theoretical models. The results of [26] suggest that the true EoS should be
hybrid (see, however, [27]). While our approximations show the largest errors
for this subclass of models, they are accurate enough to make our method
useful. It may artificially smooth out kinks and bends of the real $P-\rho$
curve, but should definitely reproduce its averaged behavior.

Other applications of the universal relations presented here could be, e.g.,
simultaneous fitting of X-ray spectra of several NSs by atmospheric models
using Eq. (1) as the form of the $M-R$ relation; or using Eqs. (3) or (4) as
an analytic form of $\beta$-equilibrium EoS in hydrodynamic studies (somewhat
smoothed, however); or, together with existing constraints, Eq. (1) could be
used as a simple informed prior for any study that requires an explicit $M-R$
relation for NSs.

Finally, we speculate that the origin of the universal relations presented
here is probably some combination of the smoothing nature of the OV mapping
from the EoS to $M-R$ and some general property of low-energy nuclear
interactions, which are dominated by just three degrees of freedom. See [47]
for discussion of the sufficient number of EoS parameters, and [22] for a
detailed discussion of possible origins of universal relations similar to Eqs.
(1)—(4).

###### Acknowledgements.

We are grateful to V. F. Suleimanov and J. Nättilä for providing us the data
for the three X-ray bursters, as well as sharing us the posteriors of the
paper [26]; to L. Brandes for sharing us the contours from [27]; to E. E.
Kolomeitsev for providing us the tables for nine KVOR EoSs; to N. Jokela and
M. Järvinen for consulting us on the correct MMNS treatment for JJ EoSs. DO
and TP were supported by Advanced ERC grant MultiJets. Work of PS was
supported by the RF Ministry of Science and Higher Education (State Assignment
No. FFUG-2024-0002).

## References

  * Degenaar and Suleimanov [2018] N. Degenaar and V. F. Suleimanov, Testing the Equation of State with Electromagnetic Observations, in  _Astrophysics and Space Science Library_, Astrophysics and Space Science Library, Vol. 457, edited by L. Rezzolla, P. Pizzochero, D. I. Jones, N. Rea, I. Vidaña, and et al. (2018) p. 185.
  * Tolman [1939] R. C. Tolman, Static Solutions of Einstein’s Field Equations for Spheres of Fluid, Physical Review 55, 364 (1939).
  * Oppenheimer and Volkoff [1939] J. R. Oppenheimer and G. M. Volkoff, On Massive Neutron Cores, Physical Review 55, 374 (1939).
  * Lindblom [1992] L. Lindblom, Determining the Nuclear Equation of State from Neutron-Star Masses and Radii, Astrophys. J. 398, 569 (1992).
  * Fiorella Burgio and Fantina [2018] G. Fiorella Burgio and A. F. Fantina, Nuclear Equation of State for Compact Stars and Supernovae, in  _The Physics and Astrophysics of Neutron Stars_, Astrophysics and Space Science Library, Vol. 457, edited by L. Rezzolla, P. Pizzochero, D. I. Jones, N. Rea, and I. Vidaña (2018) p. 255, arXiv:1804.03020 [nucl-th] .
  * Andersson and Kokkotas [1998] N. Andersson and K. D. Kokkotas, Towards gravitational wave asteroseismology, Mon. Not. R. Astron. Soc. 299, 1059 (1998), arXiv:gr-qc/9711088 [gr-qc] .
  * Lattimer and Prakash [2001] J. M. Lattimer and M. Prakash, Neutron Star Structure and the Equation of State, Astrophys. J. 550, 426 (2001), arXiv:astro-ph/0002232 [astro-ph] .
  * Bejger and Haensel [2002] M. Bejger and P. Haensel, Moments of inertia for neutron and strange stars: Limits derived for the Crab pulsar, Astron. Astrophys. 396, 917 (2002), arXiv:astro-ph/0209151 [astro-ph] .
  * Yagi and Yunes [2017] K. Yagi and N. Yunes, Approximate universal relations for neutron stars and quark stars, Phys. Rep. 681, 1 (2017), arXiv:1608.02582 [gr-qc] .
  * Ofengeim [2020] D. D. Ofengeim, Universal properties of maximum-mass neutron stars: A new tool to explore superdense matter, Phys. Rev. D 101, 103029 (2020), arXiv:2005.03549 [astro-ph.HE] .
  * Sotani and Kumar [2021] H. Sotani and B. Kumar, Universal relations between the quasinormal modes of neutron star and tidal deformability, Phys. Rev. D 104, 123002 (2021).
  * Saes and Mendes [2022] J. A. Saes and R. F. P. Mendes, Equation-of-state-insensitive measure of neutron star stiffness, Phys. Rev. D 106, 043027 (2022), arXiv:2109.11571 [gr-qc] .
  * Konstantinou and Morsink [2022] A. Konstantinou and S. M. Morsink, Universal Relations for the Increase in the Mass and Radius of a Rotating Neutron Star, Astrophys. J. 934, 139 (2022), arXiv:2206.12515 [astro-ph.HE] .
  * Cai  _et al._ [2023] B.-J. Cai, B.-A. Li, and Z. Zhang, Core States of Neutron Stars from Anatomizing Their Scaled Structure Equations, Astrophys. J. 952, 147 (2023), arXiv:2306.08202 [nucl-th] .
  * Yagi and Yunes [2013] K. Yagi and N. Yunes, I-Love-Q: Unexpected Universal Relations for Neutron Stars and Quark Stars, Science 341, 365 (2013), arXiv:1302.4499 [gr-qc] .
  * Yagi and Yunes [2016] K. Yagi and N. Yunes, Binary Love relations, Classical and Quantum Gravity 33, 13LT01 (2016), arXiv:1512.02639 [gr-qc] .
  * Legred  _et al._ [2024] I. Legred, B. O. Sy-Garcia, K. Chatziioannou, and R. Essick, Assessing equation of state-independent relations for neutron stars with nonparametric models, Phys. Rev. D 109, 023020 (2024), arXiv:2310.10854 [astro-ph.HE] .
  * Antoniadis  _et al._ [2013] J. Antoniadis, P. C. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, and et al., A Massive Pulsar in a Compact Relativistic Binary, Science 340, 448 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * Fonseca  _et al._ [2021] E. Fonseca, H. T. Cromartie, T. T. Pennucci, P. S. Ray, A. Y. Kirichenko, and et al., Refined Mass and Geometric Measurements of the High-mass PSR J0740+6620, Astrophys. J. Lett. 915, L12 (2021), arXiv:2104.00880 [astro-ph.HE] .
  * Kandel and Romani [2023] D. Kandel and R. W. Romani, An Optical Study of the Black Widow Population, Astrophys. J. 942, 6 (2023).
  * Note [1] See Supplemental Material at [URL will be inserted by publisher] which includes Refs. [48, 49, 50, 51, 10, 52, 53, 54, 55, 24, 56, 57, 58, 59, 60, 61, 62, 63, 26, 27, 30, 31, 32, 64, 65, 66, 67, 68, 69, 38, 35, 39, 65, 70].
  * Ofengeim  _et al._ [2023] D. D. Ofengeim, P. S. Shternin, and T. Piran, Maximal Mass Neutron Star as a Key to Superdense Matter Physics, Astronomy Letters 49, 567 (2023), arXiv:2310.16847 [astro-ph.HE] .
  * Note [2] We presented a similar fit in [22], but here we reduce the number of fitting coefficients.
  * Rezzolla  _et al._ [2018] L. Rezzolla, E. R. Most, and L. R. Weih, Using Gravitational-wave Observations and Quasi-universal Relations to Constrain the Maximum Mass of Neutron Stars, Astrophys. J. Lett. 852, L25 (2018).
  * Lattimer and Prakash [2016] J. M. Lattimer and M. Prakash, The equation of state of hot, dense matter and neutron stars, Phys. Rep. 621, 127 (2016), arXiv:1512.07820 [astro-ph.SR] .
  * Annala  _et al._ [2023] E. Annala, T. Gorda, J. Hirvonen, O. Komoltsev, A. Kurkela, and et al., Strongly interacting matter exhibits deconfined behavior in massive neutron stars, Nat. Commun. 14, arXiv:2303.11356 (2023).
  * Brandes  _et al._ [2023a] L. Brandes, W. Weise, and N. Kaiser, Evidence against a strong first-order phase transition in neutron star cores: Impact of new data, Phys. Rev. D 108, 094014 (2023a), arXiv:2306.06218 [nucl-th] .
  * Komoltsev and Kurkela [2022] O. Komoltsev and A. Kurkela, How Perturbative QCD Constrains the Equation of State at Neutron-Star Densities, Phys. Rev. Lett. 128, 202701 (2022), arXiv:2111.05350 [nucl-th] .
  * Hebeler  _et al._ [2013] K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, Equation of State and Neutron Star Properties Constrained by Nuclear Physics and Observation, Astrophys. J. 773, 11 (2013), arXiv:1303.4662 [astro-ph.SR] .
  * Brandes  _et al._ [2023b] L. Brandes, W. Weise, and N. Kaiser, Inference of the sound speed and related properties of neutron stars, Phys. Rev. D 107, 014011 (2023b), arXiv:2208.03026 [nucl-th] .
  * Ayriyan  _et al._ [2021] A. Ayriyan, D. Blaschke, A. G. Grunfeld, D. Alvarez-Castillo, H. Grigorian, and V. Abgaryan, Bayesian analysis of multimessenger M-R data with interpolated hybrid EoS, European Physical Journal A 57, 318 (2021), arXiv:2102.13485 [astro-ph.HE] .
  * Raaijmakers  _et al._ [2021] G. Raaijmakers, S. K. Greif, K. Hebeler, T. Hinderer, S. Nissanke, and et al., Constraints on the Dense Matter Equation of State and Neutron Star Properties from NICER’s Mass-Radius Estimate of PSR J0740+6620 and Multimessenger Observations, Astrophys. J. Lett. 918, L29 (2021), arXiv:2105.06981 [astro-ph.HE] .
  * Fan  _et al._ [2024] Y.-Z. Fan, M.-Z. Han, J.-L. Jiang, D.-S. Shao, S.-P. Tang, and et al., Maximum gravitational mass MTOV=2.2 5${}_{-0.07}^{+0.08}$M⊙ inferred at about 3% precision with multimessenger data of neutron stars, Phys. Rev. D 109, 043052 (2024), arXiv:2309.12644 [astro-ph.HE] .
  * Riley  _et al._ [2019] T. E. Riley, A. L. Watts, S. Bogdanov, P. S. Ray, R. M. Ludlam, and et al., A NICER View of PSR J0030+0451: Millisecond Pulsar Parameter Estimation, Astrophys. J. Lett. 887, L21 (2019), arXiv:1912.05702 [astro-ph.HE] .
  * Miller  _et al._ [2019] M. C. Miller, F. K. Lamb, A. J. Dittmann, S. Bogdanov, Z. Arzoumanian, K. C. Gendreau, S. Guillot, A. K. Harding, W. C. G. Ho, J. M. Lattimer, R. M. Ludlam, S. Mahmoodifar, S. M. Morsink, P. S. Ray, T. E. Strohmayer, K. S. Wood, T. Enoto, R. Foster, T. Okajima, G. Prigozhin, and Y. Soong, PSR J0030+0451 Mass and Radius from NICER Data and Implications for the Properties of Neutron Star Matter, Astrophys. J. Lett. 887, L24 (2019), arXiv:1912.05705 [astro-ph.HE] .
  * Raaijmakers  _et al._ [2019] G. Raaijmakers, T. E. Riley, A. L. Watts, S. K. Greif, S. M. Morsink, and et al., A Nicer View of PSR J0030+0451: Implications for the Dense Matter Equation of State, Astrophys. J. Lett. 887, L22 (2019), arXiv:1912.05703 [astro-ph.HE] .
  * Vinciguerra  _et al._ [2024] S. Vinciguerra, T. Salmi, A. L. Watts, D. Choudhury, T. E. Riley, and et al., An Updated Mass-Radius Analysis of the 2017-2018 NICER Data Set of PSR J0030+0451, Astrophys. J. 961, 62 (2024), arXiv:2308.09469 [astro-ph.HE] .
  * Riley  _et al._ [2021] T. E. Riley, A. L. Watts, P. S. Ray, S. Bogdanov, S. Guillot, and et al., A NICER View of the Massive Pulsar PSR J0740+6620 Informed by Radio Timing and XMM-Newton Spectroscopy, Astrophys. J. Lett. 918, L27 (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Shternin  _et al._ [2023] P. S. Shternin, D. D. Ofengeim, C. O. Heinke, and W. C. G. Ho, Constraints on neutron star superfluidity from the cooling neutron star in Cassiopeia A using all Chandra ACIS-S observations, Mon. Not. R. Astron. Soc. 518, 2775 (2023).
  * Nättilä  _et al._ [2017] J. Nättilä, M. C. Miller, A. W. Steiner, J. J. E. Kajava, V. F. Suleimanov, and et al., Neutron star mass and radius measurements from atmospheric model fits to X-ray burst cooling tail spectra, Astron. Astrophys. 608, A31 (2017), arXiv:1709.09120 [astro-ph.HE] .
  * Nättilä  _et al._ [2016] J. Nättilä, A. W. Steiner, J. J. E. Kajava, V. F. Suleimanov, and J. Poutanen, Equation of state constraints for the cold dense matter inside neutron stars using the cooling tail method, Astron. Astrophys. 591, A25 (2016), arXiv:1509.06561 [astro-ph.HE] .
  * Abbott  _et al._ [2018] B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ , GW170817: Measurements of Neutron Star Radii and Equation of State, Phys. Rev. Lett. 121, 161101 (2018), arXiv:1805.11581 [gr-qc] .
  * Demorest  _et al._ [2010] P. B. Demorest, T. Pennucci, S. M. Ransom, M. S. E. Roberts, and J. W. T. Hessels, A two-solar-mass neutron star measured using Shapiro delay, Nature (London) 467, 1081 (2010).
  * Kandel and Romani [2020] D. Kandel and R. W. Romani, Atmospheric Circulation on Black Widow Companions, Astrophys. J. 892, 101 (2020).
  * Lattimer [2021] J. M. Lattimer, Neutron Stars and the Nuclear Matter Equation of State, Annual Review of Nuclear and Particle Science 71, 433 (2021).
  * Note [3] See other explicit implementations in, e.g., [4, 71, 72].
  * Lindblom [2010] L. Lindblom, Spectral representations of neutron-star equations of state, Phys. Rev. D 82, 103011 (2010), arXiv:1009.0738 [astro-ph.HE] .
  * Typel  _et al._ [2015] S. Typel, M. Oertel, and T. Klähn, CompOSE CompStar online supernova equations of state harmonising the concert of nuclear physics and astrophysics compose.obspm.fr, Physics of Particles and Nuclei 46, 633 (2015).
  * Read  _et al._ [2009] J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, Constraints on a phenomenologically parametrized neutron-star equation of state, Phys. Rev. D 79, 124032 (2009), arXiv:0812.2163 [astro-ph] .
  * Özel and Freire [2016] F. Özel and P. Freire, Masses, Radii, and the Equation of State of Neutron Stars, Ann. Rev. Astron. Astrophys. 54, 401 (2016), arXiv:1603.02698 [astro-ph.HE] .
  * Maslov  _et al._ [2016] K. A. Maslov, E. E. Kolomeitsev, and D. N. Voskresensky, Relativistic mean-field models with scaled hadron masses and couplings: Hyperons and maximum neutron star mass, Nuc. Phys. A 950, 64 (2016), arXiv:1509.02538 [astro-ph.HE] .
  * Pearson  _et al._ [2018] J. M. Pearson, N. Chamel, A. Y. Potekhin, A. F. Fantina, C. Ducoin, A. K. Dutta, and S. Goriely, Unified equations of state for cold non-accreting neutron stars with Brussels-Montreal functionals - I. Role of symmetry energy, Mon. Not. R. Astron. Soc. 481, 2994 (2018), arXiv:1903.04981 [astro-ph.HE] .
  * Järvinen [2022] M. Järvinen, Holographic modeling of nuclear matter and neutron stars, European Physical Journal C 82, 282 (2022), arXiv:2110.08281 [hep-ph] .
  * Dietrich  _et al._ [2020] T. Dietrich, M. W. Coughlin, P. T. H. Pang, M. Bulla, J. Heinzel, L. Issa, I. Tews, and S. Antier, Multimessenger constraints on the neutron-star equation of state and the Hubble constant, Science 370, 1450 (2020), arXiv:2002.11355 [astro-ph.HE] .
  * Ascenzi  _et al._ [2024] S. Ascenzi, V. Graber, and N. Rea, Neutron-star measurements in the multi-messenger Era, Astroparticle Physics 158, 102935 (2024), arXiv:2401.14930 [astro-ph.HE] .
  * Blinnikov  _et al._ [2022] S. Blinnikov, A. Yudin, N. Kramarev, and M. Potashov, Stripping Model for Short Gamma-Ray Bursts in Neutron Star Mergers, Particles 5, 198 (2022).
  * Potekhin  _et al._ [2020] A. Y. Potekhin, D. A. Zyuzin, D. G. Yakovlev, M. V. Beznogov, and Y. A. Shibanov, Thermal luminosities of cooling neutron stars, Mon. Not. R. Astron. Soc. 496, 5052 (2020).
  * Doroshenko  _et al._ [2022] V. Doroshenko, V. Suleimanov, G. Pühlhofer, and A. Santangelo, A strangely light neutron star within a supernova remnant, Nature Astronomy 6, 1444 (2022).
  * Alford and Halpern [2023] J. A. J. Alford and J. P. Halpern, Do Central Compact Objects have Carbon Atmospheres?, Astrophys. J. 944, 36 (2023), arXiv:2302.05893 [astro-ph.HE] .
  * Steiner  _et al._ [2018] A. W. Steiner, C. O. Heinke, S. Bogdanov, C. K. Li, W. C. G. Ho, A. Bahramian, and S. Han, Constraining the mass and radius of neutron stars in globular clusters, Mon. Not. R. Astron. Soc. 476, 421 (2018), arXiv:1709.05013 [astro-ph.HE] .
  * Marino  _et al._ [2018] A. Marino, N. Degenaar, T. Di Salvo, R. Wijnands, L. Burderi, and R. Iaria, On obtaining neutron star mass and radius constraints from quiescent low-mass X-ray binaries in the Galactic plane, Mon. Not. R. Astron. Soc. 479, 3634 (2018), arXiv:1806.04557 [astro-ph.HE] .
  * Lattimer [2019] J. M. Lattimer, Neutron Star Mass and Radius Measurements, Universe 5, 159 (2019).
  * Barr  _et al._ [2024] E. D. Barr, A. Dutta, P. C. C. Freire, M. Cadelano, T. Gautam, M. Kramer, C. Pallanca, S. M. Ransom, A. Ridolfi, B. W. Stappers, T. M. Tauris, V. Venkatraman Krishnan, N. Wex, M. Bailes, J. Behrend, S. Buchner, M. Burgay, W. Chen, D. J. Champion, C. H. R. Chen, A. Corongiu, M. Geyer, Y. P. Men, P. V. Padmanabh, and A. Possenti, A pulsar in a binary with a compact object in the mass gap between neutron stars and black holes, Science 383, 275 (2024), arXiv:2401.09872 [astro-ph.HE] .
  * Schweder and Hjort [2016] T. Schweder and N. L. Hjort,  _Confidence, Likelihood, Probability: Statistical Inference with Confidence Distributions_ , Cambridge Series in Statistical and Probabilistic Mathematics (Cambridge University Press, 2016).
  * Hernandez Vivanco  _et al._ [2020] F. Hernandez Vivanco, R. Smith, E. Thrane, and P. D. Lasky, A scalable random forest regressor for combining neutron-star equation of state measurements: a case study with GW170817 and GW190425, Mon. Not. R. Astron. Soc. 499, 5972 (2020).
  * Bungert and Wacker [2020] L. Bungert and P. Wacker, The lion in the attic – A resolution of the Borel–Kolmogorov paradox, arXiv e-prints , arXiv:2009.04778 (2020), arXiv:2009.04778 [math.PR] .
  * Poole and Raftery [2000] D. Poole and A. E. Raftery, Inference for Deterministic Simulation Models: The Bayesian Melding Approach, Journal of the American Statistical Association 95, 1244 (2000).
  * Steiner [2018] A. W. Steiner, Two- and Multi-dimensional Curve Fitting using Bayesian Inference, arXiv e-prints , arXiv:1802.05339 (2018), arXiv:1802.05339 [physics.data-an] .
  * Pihajoki [2017] P. Pihajoki, A geometric approach to non-linear correlations with intrinsic scatter, Mon. Not. R. Astron. Soc. 472, 3407 (2017), arXiv:1704.05466 [astro-ph.IM] .
  * Foreman-Mackey  _et al._ [2013] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, emcee: The MCMC Hammer, Publications of the Astronomical Society of the Pacific 125, 306 (2013).
  * Lindblom and Indik [2012] L. Lindblom and N. M. Indik, Spectral approach to the relativistic inverse stellar structure problem, Phys. Rev. D 86, 084003 (2012), arXiv:1207.3744 [astro-ph.HE] .
  * Soma  _et al._ [2022] S. Soma, L. Wang, S. Shi, H. Stöcker, K. Zhou, and et al., Neural network reconstruction of the dense matter equation of state from neutron star observables, J. Cosmology & Astropaticles. 2022, 071 (2022), arXiv:2201.01756 [hep-ph] .
  * Haensel  _et al._ [2007] P. Haensel, A. Y. Potekhin, and D. G. Yakovlev,  _Neutron Stars 1 : Equation of State and Structure_ , Astrophysics and Space Science Library, Vol. 326 (Springer, New York, 2007).
  * Gusakov  _et al._ [2005] M. E. Gusakov, A. D. Kaminker, D. G. Yakovlev, and O. Y. Gnedin, The cooling of Akmal-Pandharipande-Ravenhall neutron star models, Mon. Not. R. Astron. Soc. 363, 555 (2005), arXiv:astro-ph/0507560 [astro-ph] .
  * Kaminker  _et al._ [2014] A. D. Kaminker, A. A. Kaurov, A. Y. Potekhin, and D. G. Yakovlev, Thermal emission of neutron stars with internal heaters, Mon. Not. R. Astron. Soc. 442, 3484 (2014), arXiv:1406.0723 [astro-ph.HE] .
  * Potekhin  _et al._ [2013] A. Y. Potekhin, A. F. Fantina, N. Chamel, J. M. Pearson, and S. Goriely, Analytical representations of unified equations of state for neutron-star matter, Astron. Astrophys. 560, A48 (2013), arXiv:1310.0049 [astro-ph.SR] .
  * Yakovlev  _et al._ [2011] D. G. Yakovlev, W. C. G. Ho, P. S. Shternin, C. O. Heinke, and A. Y. Potekhin, Cooling rates of neutron stars and the young neutron star in the Cassiopeia A supernova remnant, Mon. Not. R. Astron. Soc. 411, 1977 (2011), arXiv:1010.1154 [astro-ph.HE] .
  * Ofengeim  _et al._ [2019] D. D. Ofengeim, M. E. Gusakov, P. Haensel, and M. Fortin, Bulk viscosity in neutron stars with hyperon cores, Phys. Rev. D 100, 103017 (2019), arXiv:1911.08407 [astro-ph.HE] .
  * Gusakov  _et al._ [2014] M. E. Gusakov, P. Haensel, and E. M. Kantor, Physics input for modelling superfluid neutron stars with hyperon cores, Mon. Not. R. Astron. Soc. 439, 318 (2014), arXiv:1401.2827 [astro-ph.HE] .

##  Supplementary A Collection of the equations of state

The EoSs for this work were selected as follows. From the CompOSE database
[48] (as of November 2022), we took the subset Cold Neutron Star EoS, keeping
only the models providing data for both the NS crust and core. Then we added
models from [49] (some time ago, this EoS collection was available by the link
given in [50]; now one can access an archived copy), nine models from [51],
and several dozen models more which are available to the authors (see [10] for
details). We discarded duplicates and tables ended at
$\rho<\rho_{\text{{tov}}}$. Ultimately, we got a sample of 169 models listed
in Table LABEL:tab:zooNuc.

The considered models are based on a variety of approaches to modeling
interactions and the microphysics of superdense matter. These include models
based on effective energy density functionals, including purely
phenomenological (PAL, PAPAL, BGN families), non-relativistic Skyrme-type
(SLy, BSk, SkI, etc.) or Gogni (D1M*) and others, numerous relativistic mean
field functionals; models derived from microscopic baryon interaction
potentials using many-body methods (APR, WFF, BBB); and also some other
models. Among the hybrid equations of state, both those in which a first-order
phase transition occurs between hadronic and quark matter (for example, CMF,
VQCD) and those in which there is a quark-hadron crossover (QHC) between these
phases are present.

We name models according to the CompOSE tradition, following the principle
original paper acronym(model name). For each EoS, tables and detailed
description (choice of microphysical parameters, stitching the core and the
crust, etc.) can be found in the references given in the “source” column. One
can find a decoding of the acronym of the original work proposed in this
particular model (typically not coinciding with the source). If neither the
source nor the references given there describe a crust model, the BSk24 crust
[52] was stitched to the core under the condition of pressure continuity.
Model names in Table LABEL:tab:zooNuc may be reasonably abbreviated with
respect to those used in the references.

Additionally, for each EoS, the mass $M_{\text{{tov}}}$, its radius
$R_{\text{{tov}}}$, pressure $P_{\text{{tov}}}$, density
$\rho_{\text{{tov}}}$, and the speed of sound $c_{\text{s{tov}}}$ at the
center of such MMNS, as well as the radius of the half-maximum-mass NS
$R_{\text{1/2}}$ are indicated in the table LABEL:tab:zooNuc. The relative
root mean square (rrms) and maximum relative fitting errors of the $M-R$ (Eq.
1) and $P(\rho)$ (Eq. 3) fits are also provided. For the latter fit, we give
the errors in two density ranges $\rho\in(3\rho_{0}...\rho_{\text{{tov}}})$
and $\rho\in(1\rho_{0}...\rho_{\text{{tov}}})$ called “all $\rho$”.

In Fig. S1 we show distributions of this sample of models over
$M_{\text{{tov}}}$ (panel a), over $c_{\text{s{tov}}}$ (b), and over the
parameter $r_{\text{1/2}}=R_{\text{1/2}}/R_{\text{{tov}}}$ (c). In Fig. S2 we
show distributions of root-mean-square and maximal relative errors of
approximations (1), (3) and (4) (panels a, b, and c, correspondingly).

Figure S1: Distributions of (a) NS maximal mass, (b) speed of sound in the
center of the star, (c) $r_{1/2}=R_{1/2}/R_{\text{{tov}}}$ ratio for the
sample of equations of state used.

Figure S2: Distributions of errors of $M-R$ approximation (a), and $P(\rho)$ approximation in the high-density domain $3\rho_{0}<\rho<\rho_{\text{{tov}}}$ (b) and for all densities $\rho_{0}<\rho<\rho_{\text{{tov}}}$ (c). Panels (a) and (c) show the same as the insets in Figs. 1 and 3 in the main text, correspondingly. Table S1:  Neutron star equations of state  Name | Source | $M_{\text{{tov}}}$ | $R_{\text{{tov}}}$ | $R_{\text{1/2}}$ | Fit $R(M)$ [%] | $\rho_{\text{{tov}}}$ | $P_{\text{{tov}}}$ | $c_{\text{s{tov}}}$ | Fit $P(\rho)$ [%]
---|---|---|---|---|---|---|---|---|---
|  |  |  |  |  |  |  |  |  | $\rho>3\rho_{0}$ | all $\rho$
|  | [$M_{\odot}$] | [km] | [km] | rrms | max | [$\rho_{0}$] | [$\rho_{0}c^{2}$] | [c] | rrms | max | rrms | max
Nucleonic models
APR(APR1) |  [49, 50]${}^{\text{a}}$ | 1.68 | 8.28 | 9.70 | 0.89 | 1.1 | 14.9 | 7.42 | 0.839 | 8.9 | 13 | 9.1 | 17
APR(APR2) |  [49, 50]${}^{\text{a}}$ | 1.81 | 8.72 | 10.51 | 0.50 | 0.62 | 13.4 | 7.63 | 1.03 | 5.9 | 10 | 7.4 | 31
APR(APR3) |  [49, 50]${}^{\text{a}}$ | 2.39 | 10.76 | 12.02 | 0.32 | 0.54 | 8.35 | 5.78 | 1.13 | 0.90 | 1.5 | 4.1 | 25
ZBL(BPAL12) |  [49, 50]${}^{\text{a}}$ | 1.45 | 9.01 | 12.69 | 1.1 | 1.3 | 14.2 | 4.29 | 0.704 | 2.1 | 5.1 | 2.5 | 6.2
EOHJBO(ENG) |  [49, 50]${}^{\text{a}}$ | 2.24 | 10.41 | 11.97 | 0.42 | 0.87 | 9.19 | 5.68 | 1.00 | 3.8 | 6.6 | 4.3 | 7.9
FP(FPS) |  [49, 50]${}^{\text{a}}$ | 1.80 | 9.28 | 11.18 | 0.20 | 0.29 | 12.1 | 5.40 | 0.870 | 3.0 | 6.0 | 3.2 | 6.4
PAL(PAL6) |  [49, 50]${}^{\text{a}}$ | 1.48 | 9.27 | 13.07 | 1.0 | 1.2 | 13.5 | 3.93 | 0.692 | 1.2 | 2.7 | 2.0 | 11
WFF(WFF1) |  [49, 50]${}^{\text{a}}$ | 2.13 | 9.41 | 10.35 | 0.68 | 0.96 | 10.8 | 8.04 | 1.18 | 5.9 | 11 | 12 | 72
WFF(WFF2) |  [49, 50]${}^{\text{a}}$ | 2.20 | 9.81 | 11.14 | 0.60 | 1.0 | 10.1 | 7.37 | 1.15 | 3.2 | 5.1 | 7.0 | 44
WFF(WFF3) |  [49, 50]${}^{\text{a}}$ | 1.84 | 9.50 | 11.06 | 0.24 | 0.36 | 11.3 | 4.86 | 0.836 | 3.0 | 5.3 | 5.6 | 20
MPA(MPA1) |  [49, 50]${}^{\text{a}}$ | 2.46 | 11.32 | 12.38 | 0.40 | 0.66 | 7.55 | 4.56 | 0.991 | 2.7 | 3.9 | 4.2 | 8.8
MS(MS1) |  [49, 50]${}^{\text{a}}$ | 2.77 | 13.33 | 14.83 | 0.71 | 1.1 | 5.53 | 2.80 | 0.892 | 2.6 | 4.9 | 6.8 | 33
MS(MS1b) |  [49, 50]${}^{\text{a}}$ | 2.78 | 13.29 | 14.52 | 0.44 | 0.74 | 5.51 | 2.80 | 0.894 | 3.0 | 6.2 | 8.2 | 31
MS(MS2) |  [49, 50]${}^{\text{a}}$ | 1.81 | 11.68 | 14.24 | 0.47 | 0.69 | 7.79 | 1.89 | 0.582 | 0.76 | 2.4 | 2.5 | 5.2
ABHT(QMC-RMF1) |  [48]${}^{\text{b}}$ | 1.95 | 10.25 | 12.02 | 0.37 | 0.47 | 9.82 | 3.96 | 0.806 | 1.5 | 2.5 | 7.0 | 24
ABHT(QMC-RMF2) |  [48]${}^{\text{b}}$ | 2.04 | 10.49 | 12.07 | 0.26 | 0.40 | 9.26 | 3.93 | 0.828 | 1.5 | 2.5 | 7.2 | 27
ABHT(QMC-RMF3) |  [48]${}^{\text{b}}$ | 2.14 | 10.71 | 12.32 | 0.66 | 0.92 | 8.81 | 4.11 | 0.862 | 2.2 | 3.3 | 6.5 | 17
ABHT(QMC-RMF4) |  [48]${}^{\text{b}}$ | 2.21 | 11.03 | 12.22 | 0.13 | 0.25 | 8.15 | 3.69 | 0.852 | 1.8 | 2.6 | 6.0 | 30
APR(APR) |  [48]${}^{\text{b}}$ | 2.19 | 9.93 | 11.37 | 0.53 | 0.75 | 9.93 | 6.87 | 1.15 | 2.0 | 3.3 | 6.6 | 42
BG(BGN1) | [73] | 2.17 | 10.91 | 13.13 | 0.14 | 0.27 | 8.66 | 4.24 | 0.939 | 2.0 | 2.7 | 2.5 | 8.0
BG(BGN2) | [73] | 2.48 | 11.74 | 13.51 | 0.11 | 0.22 | 7.21 | 4.23 | 1.04 | 2.3 | 3.1 | 4.5 | 16
BBB(BBB1) | [73] | 1.79 | 9.67 | 11.23 | 0.32 | 0.42 | 11.0 | 4.05 | 0.765 | 3.0 | 6.1 | 8.7 | 29
BBB(BHF-BBB2) |  [48]${}^{\text{b}}$ | 1.92 | 9.50 | 11.37 | 0.56 | 0.92 | 11.4 | 5.75 | 0.915 | 2.6 | 3.7 | 3.9 | 18
BL(chiral) |  [48]${}^{\text{b}}$ | 2.08 | 10.28 | 12.56 | 0.41 | 0.54 | 9.78 | 5.13 | 0.985 | 1.2 | 1.8 | 2.2 | 11
GDTB(DDH$\delta$) |  [48]${}^{\text{b}}$ | 2.15 | 11.23 | 12.49 | 0.45 | 0.65 | 7.92 | 3.13 | 0.790 | 0.73 | 1.9 | 2.1 | 5.5
GKYG(APR-I) | [74] | 1.92 | 10.31 | 12.42 | 0.31 | 0.41 | 9.85 | 3.86 | 0.801 | 0.34 | 0.91 | 3.2 | 13
GKYG(APR-II) | [74] | 1.92 | 10.27 | 12.34 | 0.29 | 0.39 | 9.91 | 3.89 | 0.802 | 0.43 | 0.85 | 3.3 | 13
GKYG(APR-III) | [74] | 1.93 | 10.38 | 12.57 | 0.35 | 0.45 | 9.75 | 3.80 | 0.799 | 0.27 | 1.1 | 3.0 | 11
GM(GM1) |  [48]${}^{\text{b}}$ | 2.38 | 12.16 | 14.23 | 0.39 | 0.63 | 6.91 | 3.07 | 0.848 | 0.90 | 1.8 | 1.6 | 6.4
GMSR(BSk14) |  [48]${}^{\text{b}}$ | 1.93 | 9.29 | 11.50 | 0.059 | 0.076 | 11.9 | 7.10 | 1.13 | 5.2 | 14 | 11 | 31
GMSR(H7) |  [48]${}^{\text{b}}$ | 2.48 | 11.23 | 12.61 | 1.4 | 2.6 | 7.64 | 5.64 | 1.38 | 10 | 16 | 10 | 17
GMSR(LNS5) |  [48]${}^{\text{b}}$ | 1.97 | 9.92 | 11.65 | 0.50 | 0.75 | 10.4 | 4.84 | 0.890 | 1.9 | 3.2 | 2.2 | 3.6
GMSR(SLy5) |  [48]${}^{\text{b}}$ | 2.10 | 10.06 | 11.93 | 0.48 | 0.74 | 10.0 | 5.67 | 1.00 | 0.90 | 1.1 | 2.4 | 13
GPPVA(DD2) |  [48]${}^{\text{b}}$ | 2.42 | 11.88 | 13.13 | 0.26 | 0.46 | 6.99 | 3.30 | 0.866 | 2.0 | 3.1 | 2.4 | 5.5
GPPVA(DDME2) |  [48]${}^{\text{b}}$ | 2.48 | 12.07 | 13.16 | 0.23 | 0.43 | 6.71 | 3.24 | 0.873 | 2.7 | 4.7 | 3.4 | 7.5
GPPVA(FSU2) |  [48]${}^{\text{b}}$ | 2.07 | 12.06 | 14.10 | 0.98 | 1.3 | 7.11 | 2.08 | 0.632 | 2.3 | 3.2 | 2.9 | 6.3
GPPVA(FSU2H) |  [48]${}^{\text{b}}$ | 2.38 | 12.36 | 13.16 | 0.38 | 0.67 | 6.36 | 2.33 | 0.703 | 6.1 | 11 | 7.1 | 11
GPPVA(FSU2R) |  [48]${}^{\text{b}}$ | 2.05 | 11.65 | 12.87 | 0.40 | 0.64 | 7.35 | 2.18 | 0.632 | 4.2 | 6.1 | 4.7 | 15
GPPVA(NL3$\omega\rho$L55) |  [48]${}^{\text{b}}$ | 2.75 | 13.00 | 13.71 | 0.23 | 0.46 | 5.66 | 2.92 | 0.894 | 4.6 | 9.6 | 7.4 | 13
GPPVA(TM1e) |  [48]${}^{\text{b}}$ | 2.12 | 11.84 | 13.07 | 0.37 | 0.61 | 7.12 | 2.24 | 0.654 | 3.9 | 5.9 | 4.9 | 19
GPPVA(TW) |  [48]${}^{\text{b}}$ | 2.08 | 10.64 | 12.48 | 0.60 | 0.81 | 9.05 | 3.95 | 0.839 | 1.3 | 2.3 | 3.9 | 10
KKPY(APR-IV) | [75] | 2.16 | 10.86 | 12.81 | 0.23 | 0.38 | 8.68 | 4.12 | 0.904 | 0.72 | 1.1 | 1.0 | 3.4
PCGS(PCSB0) |  [48]${}^{\text{b}}$ | 2.53 | 12.13 | 13.26 | 0.54 | 0.81 | 6.63 | 3.35 | 0.888 | 3.5 | 5.3 | 4.1 | 9.5
PFCPG(BSk19) | [76] | 1.86 | 9.10 | 11.08 | 0.35 | 0.47 | 12.4 | 6.76 | 1.01 | 4.3 | 8.9 | 5.5 | 14
PFCPG(BSk20) | [76] | 2.16 | 10.17 | 11.76 | 0.16 | 0.33 | 9.61 | 5.79 | 1.06 | 0.98 | 2.0 | 1.9 | 12
PFCPG(BSk21) | [76] | 2.27 | 11.04 | 12.52 | 0.23 | 0.34 | 8.17 | 4.24 | 0.963 | 0.96 | 1.5 | 1.3 | 4.1
PCPFDDG(BSk22) | [52] | 2.26 | 11.20 | 13.03 | 0.19 | 0.29 | 8.07 | 4.02 | 0.951 | 2.0 | 2.9 | 2.4 | 7.2
PCPFDDG(BSk24) | [52] | 2.28 | 11.07 | 12.51 | 0.12 | 0.19 | 8.10 | 4.15 | 0.955 | 0.74 | 1.1 | 1.2 | 3.1
PCPFDDG(BSk25) | [52] | 2.22 | 11.04 | 12.25 | 0.20 | 0.28 | 8.11 | 3.76 | 0.896 | 1.4 | 2.8 | 4.2 | 17
PCPFDDG(BSk26) | [52] | 2.17 | 10.19 | 11.79 | 0.21 | 0.38 | 9.57 | 5.76 | 1.06 | 0.86 | 1.7 | 1.9 | 12
RG(KDE0v) |  [48]${}^{\text{b}}$ | 1.96 | 9.66 | 11.71 | 0.42 | 0.54 | 11.0 | 5.82 | 0.984 | 2.0 | 3.9 | 2.2 | 6.5
RG(KDE0v1) |  [48]${}^{\text{b}}$ | 1.97 | 9.79 | 11.96 | 0.52 | 0.66 | 10.8 | 5.53 | 0.965 | 1.4 | 2.8 | 1.6 | 6.3
RG(Rs) |  [48]${}^{\text{b}}$ | 2.12 | 10.76 | 13.11 | 0.13 | 0.19 | 8.96 | 4.23 | 0.922 | 2.4 | 3.3 | 4.6 | 20
RG(Sk255) |  [48]${}^{\text{b}}$ | 2.14 | 10.85 | 13.44 | 0.46 | 0.63 | 8.88 | 4.32 | 0.937 | 2.3 | 3.0 | 2.6 | 10
RG(Sk272) |  [48]${}^{\text{b}}$ | 2.23 | 11.09 | 13.51 | 0.41 | 0.57 | 8.40 | 4.30 | 0.965 | 2.3 | 3.0 | 2.8 | 11
RG(Ska) |  [48]${}^{\text{b}}$ | 2.21 | 10.89 | 13.03 | 0.20 | 0.35 | 8.63 | 4.46 | 0.970 | 1.8 | 2.5 | 3.4 | 16
RG(Skb) |  [48]${}^{\text{b}}$ | 2.19 | 10.61 | 12.15 | 0.40 | 0.59 | 8.85 | 4.64 | 0.974 | 1.0 | 1.5 | 5.5 | 25
RG(SkI2) |  [48]${}^{\text{b}}$ | 2.16 | 11.12 | 13.64 | 0.27 | 0.41 | 8.45 | 3.89 | 0.914 | 3.5 | 4.8 | 5.0 | 18
RG(SkI3) |  [48]${}^{\text{b}}$ | 2.24 | 11.31 | 13.65 | 0.22 | 0.33 | 8.08 | 3.89 | 0.937 | 3.4 | 4.5 | 4.5 | 15
RG(SkI4) |  [48]${}^{\text{b}}$ | 2.17 | 10.67 | 12.39 | 0.16 | 0.24 | 8.86 | 4.47 | 0.949 | 0.86 | 1.3 | 2.3 | 11
RG(SkI5) |  [48]${}^{\text{b}}$ | 2.24 | 11.47 | 14.22 | 0.36 | 0.53 | 7.96 | 3.76 | 0.932 | 4.7 | 6.3 | 6.0 | 20
RG(SkI6) |  [48]${}^{\text{b}}$ | 2.19 | 10.76 | 12.51 | 0.12 | 0.18 | 8.71 | 4.42 | 0.951 | 0.94 | 1.4 | 2.1 | 10
RG(SkMp) |  [48]${}^{\text{b}}$ | 2.11 | 10.53 | 12.64 | 0.10 | 0.19 | 9.26 | 4.58 | 0.946 | 1.5 | 2.2 | 4.2 | 20
RG(SkOp) |  [48]${}^{\text{b}}$ | 1.97 | 10.13 | 12.47 | 0.31 | 0.43 | 10.2 | 4.72 | 0.906 | 0.91 | 1.4 | 2.0 | 10
RG(SLy2) |  [48]${}^{\text{b}}$ | 2.05 | 10.05 | 11.95 | 0.31 | 0.46 | 10.1 | 5.34 | 0.977 | 0.90 | 1.8 | 1.2 | 6.4
RG(SLy4) |  [48]${}^{\text{b}}$ | 2.05 | 9.99 | 11.87 | 0.30 | 0.44 | 10.2 | 5.46 | 0.985 | 1.1 | 2.1 | 1.4 | 6.9
RG(SLy9) |  [48]${}^{\text{b}}$ | 2.16 | 10.63 | 12.57 | 0.23 | 0.37 | 9.01 | 4.59 | 0.952 | 0.69 | 1.1 | 1.1 | 5.8
RG(SLy230a) |  [48]${}^{\text{b}}$ | 2.10 | 10.25 | 11.91 | 0.31 | 0.48 | 9.60 | 4.95 | 0.948 | 1.1 | 1.8 | 2.2 | 5.6
VGBCMR(D1M*) |  [48]${}^{\text{b}}$ | 2.00 | 10.21 | 11.77 | 0.22 | 0.37 | 9.72 | 4.28 | 0.855 | 1.6 | 2.7 | 4.4 | 15
XMLSLZ(DD-LZ1) |  [48]${}^{\text{b}}$ | 2.56 | 12.28 | 13.05 | 0.14 | 0.28 | 6.38 | 3.12 | 0.875 | 3.7 | 7.5 | 5.9 | 16
XMLSLZ(DDME-X) |  [48]${}^{\text{b}}$ | 2.56 | 12.35 | 13.28 | 0.15 | 0.33 | 6.36 | 3.10 | 0.876 | 3.1 | 6.2 | 4.1 | 7.1
XMLSLZ(GM1) |  [48]${}^{\text{b}}$ | 2.36 | 11.94 | 13.73 | 0.21 | 0.43 | 7.08 | 3.17 | 0.851 | 0.52 | 1.2 | 4.0 | 20
XMLSLZ(MTVTC) |  [48]${}^{\text{b}}$ | 2.02 | 10.90 | 13.32 | 0.14 | 0.26 | 8.87 | 3.46 | 0.803 | 2.0 | 2.9 | 2.5 | 9.0
XMLSLZ(NL3) |  [48]${}^{\text{b}}$ | 2.77 | 13.28 | 14.57 | 0.64 | 0.95 | 5.53 | 2.81 | 0.890 | 3.2 | 6.3 | 7.2 | 33
XMLSLZ(PK1) |  [48]${}^{\text{b}}$ | 2.31 | 12.66 | 14.40 | 0.82 | 1.2 | 6.33 | 2.14 | 0.682 | 3.1 | 4.4 | 3.4 | 12
XMLSLZ(PKDD) |  [48]${}^{\text{b}}$ | 2.33 | 11.77 | 13.69 | 0.61 | 0.85 | 7.33 | 3.30 | 0.852 | 0.31 | 1.0 | 2.5 | 7.5
XMLSLZ(TM1) |  [48]${}^{\text{b}}$ | 2.18 | 12.37 | 14.37 | 0.72 | 1.1 | 6.72 | 2.11 | 0.659 | 1.8 | 2.4 | 1.9 | 7.9
XMLSLZ(TW99) |  [48]${}^{\text{b}}$ | 2.08 | 10.62 | 12.37 | 0.49 | 0.68 | 9.06 | 3.96 | 0.839 | 1.3 | 2.1 | 3.4 | 8.8
YHSHP(PAL1_120) | [77] | 1.47 | 9.18 | 13.08 | 1.3 | 1.4 | 13.8 | 4.09 | 0.699 | 1.5 | 3.7 | 2.0 | 8.0
YHSHP(PAL1_180) | [77] | 1.74 | 9.93 | 13.06 | 0.60 | 0.71 | 11.2 | 3.97 | 0.771 | 1.1 | 1.7 | 1.8 | 8.5
YHSHP(PAL1_240) | [77] | 1.95 | 10.60 | 13.18 | 0.28 | 0.38 | 9.50 | 3.70 | 0.811 | 1.5 | 1.9 | 2.0 | 6.1
YHSHP(PAL2_120) | [77] | 1.48 | 9.72 | 14.32 | 0.32 | 0.49 | 12.6 | 3.47 | 0.677 | 2.6 | 3.8 | 2.7 | 3.9
YHSHP(PAL2_180) | [77] | 1.75 | 10.35 | 14.11 | 0.18 | 0.33 | 10.5 | 3.52 | 0.753 | 3.1 | 4.1 | 2.9 | 6.0
YHSHP(PAL2_240) | [77] | 1.97 | 10.97 | 14.06 | 0.14 | 0.27 | 9.03 | 3.35 | 0.795 | 3.2 | 4.3 | 3.1 | 4.3
YHSHP(PAL2_400_0) | [77] | 2.85 | 12.69 | 14.40 | 0.19 | 0.41 | 6.00 | 4.49 | 1.23 | 3.1 | 4.2 | 10 | 42
YHSHP(PAL3_120) | [77] | 1.42 | 8.45 | 11.56 | 1.4 | 1.6 | 15.9 | 5.14 | 0.729 | 4.8 | 9.0 | 4.8 | 9.0
YHSHP(PAL3_180) | [77] | 1.69 | 9.36 | 11.89 | 0.55 | 0.63 | 12.3 | 4.66 | 0.799 | 1.8 | 3.6 | 2.5 | 11
YHSHP(PAL3_240) | [77] | 1.91 | 10.12 | 12.22 | 0.16 | 0.26 | 10.2 | 4.18 | 0.828 | 0.50 | 0.99 | 2.4 | 12
YHSHP(PAL3_260_01) | [77] | 2.06 | 10.41 | 12.32 | 0.12 | 0.24 | 9.47 | 4.46 | 0.907 | 0.52 | 0.81 | 1.2 | 6.2
YHSHP(PAL3_300) | [77] | 1.83 | 10.51 | 12.43 | 0.16 | 0.23 | 9.47 | 3.02 | 0.709 | 0.40 | 0.66 | 4.3 | 24
YHSHP(PAL3_300_0) | [77] | 2.34 | 11.08 | 12.68 | 0.11 | 0.19 | 8.07 | 4.69 | 1.04 | 1.2 | 1.7 | 1.8 | 4.8
YHSHP(PAL3_400) | [77] | 2.83 | 12.32 | 13.31 | 0.26 | 0.70 | 6.20 | 4.92 | 1.26 | 1.1 | 1.7 | 6.3 | 18
YHSHP(PAPAL_120) | [77] | 1.44 | 8.74 | 12.24 | 1.5 | 1.7 | 15.0 | 4.70 | 0.717 | 3.3 | 6.8 | 3.5 | 6.8
YHSHP(PAPAL_180) | [77] | 1.71 | 9.59 | 12.41 | 0.67 | 0.77 | 11.9 | 4.37 | 0.786 | 1.0 | 2.1 | 1.8 | 8.6
YHSHP(PAPAL_240) | [77] | 1.93 | 10.32 | 12.64 | 0.27 | 0.36 | 9.92 | 3.98 | 0.822 | 0.61 | 0.95 | 1.9 | 8.3
MKV(KVOR) | [51] | 2.02 | 10.72 | 12.92 | 0.40 | 0.55 | 9.10 | 3.68 | 0.813 | 0.58 | 1.9 | 2.2 | 5.5
MKV(KVORcut02) | [51] | 2.27 | 12.61 | 13.74 | 0.32 | 0.47 | 6.23 | 1.99 | 0.692 | 1.8 | 5.7 | 9.1 | 22
MKV(KVORcut03) | [51] | 2.17 | 11.58 | 12.97 | 0.60 | 0.77 | 7.49 | 2.74 | 0.741 | 2.9 | 5.0 | 5.4 | 15
MKV(KVORcut04) | [51] | 2.09 | 10.97 | 12.90 | 1.2 | 1.6 | 8.54 | 3.42 | 0.780 | 2.8 | 6.4 | 5.1 | 10
MKV(MKVOR) | [51] | 2.34 | 11.33 | 12.07 | 0.37 | 0.57 | 7.52 | 3.54 | 0.846 | 5.3 | 8.0 | 10 |
${}^{\text{a}}$accessible via the archive link
https://web.archive.org/web/20200830055909/http://xtreme.as.arizona.edu/neutronstars/
${}^{\text{b}}$CompStar Online Supernova Equations of State database,
https://compose.obspm.fr ${}^{\text{c}}$See also [53]

Models with hyperons & $\Delta$-isobars G(GNH3) [49, 50]${}^{\text{a}}$ 1.96
11.39 14.62 0.49 0.81 8.49 2.85 0.747 4.2 6.9 4.7 14 LNO(H1) [49,
50]${}^{\text{a}}$ 1.55 10.97 13.88 1.5 1.9 9.16 1.95 0.541 3.4 5.4 4.3 17
LNO(H2) [49, 50]${}^{\text{a}}$ 1.67 11.51 13.87 1.1 1.5 8.02 1.71 0.564 2.8
4.4 3.5 9.8 LNO(H3) [49, 50]${}^{\text{a}}$ 1.79 11.85 14.04 0.52 0.74 7.52
1.68 0.518 1.6 2.6 3.2 9.3 LNO(H4) [49, 50]${}^{\text{a}}$ 2.03 11.76 13.69
0.77 1.1 7.51 2.30 0.653 1.8 3.3 3.7 11 LNO(H5) [49, 50]${}^{\text{a}}$ 1.73
11.31 13.38 1.3 1.8 8.19 1.92 0.596 2.9 4.4 3.7 9.0 LNO(H7) [49,
50]${}^{\text{a}}$ 1.68 10.85 13.42 1.2 1.7 9.18 2.32 0.630 3.4 5.1 3.5 8.8
RSGMT(QMC700) [50]${}^{\text{a}}$ 1.98 12.36 12.70 2.4 3.5 5.99 1.17 0.404 26
42 23 42 BG(BGN1H1) [49, 50]${}^{\text{a}}$ 1.63 9.33 13.17 7.8 12 13.6 5.48
0.876 10 48 26 80 BG(BGN1H2) [73] 1.57 9.04 13.00 9.8 14 14.6 5.95 0.897 12 53
32 100 BG(BGN2H1) [73] 1.82 9.57 13.43 6.3 11 12.4 6.36 1.01 6.8 27 34 120
BG(BGN2H2) [73] 1.74 9.12 13.42 7.0 11 13.9 7.45 1.05 10 41 45 160 DNS(CMF)
[48]${}^{\text{b}}$ 2.07 11.88 13.51 0.32 0.44 7.19 2.19 0.686 0.96 1.9 4.6 18
DS(CMF)-1 [48]${}^{\text{b}}$ 2.07 11.88 13.51 0.32 0.44 7.20 2.19 0.686 0.96
2.0 4.6 18 DS(CMF)-2 [48]${}^{\text{b}}$ 2.13 12.00 13.60 0.39 0.57 7.05 2.25
0.692 0.33 1.1 5.6 20 DS(CMF)-3 [48]${}^{\text{b}}$ 2.00 11.54 13.17 0.17 0.21
7.65 2.32 0.685 1.1 2.5 2.0 7.5 DS(CMF)-4 [48]${}^{\text{b}}$ 2.05 11.59 13.24
0.21 0.31 7.59 2.41 0.689 0.37 0.96 2.3 8.5 DS(CMF)-5 [48]${}^{\text{b}}$ 2.07
11.48 13.20 0.36 0.45 7.74 2.65 0.753 1.6 2.5 2.3 8.3 DS(CMF)-6
[48]${}^{\text{b}}$ 2.11 11.57 13.27 0.33 0.47 7.62 2.65 0.744 0.79 1.1 2.6
9.4 DS(CMF)-7 [48]${}^{\text{b}}$ 2.07 11.49 13.20 0.35 0.43 7.73 2.64 0.752
1.6 2.5 2.3 8.3 DS(CMF)-8 [48]${}^{\text{b}}$ 2.09 11.58 13.26 0.39 0.50 7.60
2.58 0.744 1.5 2.3 2.8 9.0 GHK(GM1A) [78] 1.99 11.94 13.63 1.2 2.0 7.18 2.02
0.686 4.4 6.3 6.7 15 GHK(TM1C) [78] 2.05 12.47 14.41 0.25 0.39 6.58 1.77 0.644
1.8 2.7 3.3 8.2 GHK(GM1’B) [79] 2.01 11.46 14.00 0.76 1.4 8.13 2.78 0.748 4.4
6.6 4.4 6.6 OGHF(FSU2H) [78] 1.99 11.98 13.03 0.96 1.5 6.90 1.83 0.637 2.0 5.6
10 29 OGHF(NL3$\omega\rho$) [78] 2.71 12.94 13.69 0.090 0.22 5.73 2.85 0.889
3.5 8.3 7.5 15 OPGR(GM1Y4) [48]${}^{\text{b}}$ 1.79 13.01 13.52 1.2 1.6 5.23
0.823 0.424 6.7 17 33 76 OPGR(GM1Y5) [48]${}^{\text{b}}$ 2.12 12.30 13.64 0.45
0.91 6.61 1.91 0.667 2.2 6.6 8.3 24 OPGR(GM1Y6) [48]${}^{\text{b}}$ 2.29 12.12
13.69 0.34 0.49 6.84 2.64 0.808 2.3 3.5 5.5 22 R(DD2Y$\Delta$)1.1-1.1
[48]${}^{\text{b}}$ 2.04 11.21 13.11 0.26 0.41 8.22 2.89 0.741 0.43 1.4 3.6 17
R(DD2Y$\Delta$)1.2-1.1 [48]${}^{\text{b}}$ 2.05 10.98 12.41 0.99 1.2 8.39 3.05
0.731 4.0 6.8 9.8 32 R(DD2Y$\Delta$)1.2-1.3 [48]${}^{\text{b}}$ 2.03 11.44
13.14 1.2 1.7 7.89 2.63 0.739 2.7 4.0 5.9 18 MKV(KVORHcut03$\phi$) [51] 1.98
11.76 12.98 0.79 1.1 7.17 1.93 0.630 4.4 13 6.9 14
MKV(KVORHcut03H$\phi\sigma$) [51] 2.08 11.94 12.97 1.4 2.0 6.85 1.93 0.621 8.6
16 8.8 16 MKV(MKVORH$\phi$) [51] 2.23 11.23 12.03 0.13 0.27 7.73 3.23 0.790
4.4 7.7 10 29 MKV(MKVORH$\phi\sigma$) [51] 2.30 11.43 12.05 0.60 0.89 7.35
3.09 0.782 7.2 13 12 30 Hybrid models ABPR(ALF1) [49, 50]${}^{\text{a}}$ 1.49
9.21 10.35 1.4 1.6 11.9 2.95 0.565 9.0 22 17 51 ABPR(ALF2) [49,
50]${}^{\text{a}}$ 2.09 11.96 12.86 0.58 0.87 6.94 1.95 0.553 5.7 7.4 15 45
ABPR(ALF3) [49, 50]${}^{\text{a}}$ 1.47 9.50 11.43 1.7 2.0 11.6 2.75 0.565 5.5
12 11 37 ABPR(ALF4) [49, 50]${}^{\text{a}}$ 1.94 10.88 11.58 1.6 2.3 8.28 2.33
0.506 14 21 17 35 PCL(PCL2) [49, 50]${}^{\text{a}}$ 1.48 10.18 13.44 2.1 2.6
10.8 2.54 0.599 4.6 7.6 5.5 21 BFH(QHC19-A) [48]${}^{\text{b}}$ 1.93 10.23
11.49 0.32 0.59 9.63 3.66 0.781 3.3 7.8 8.0 22 BFH(QHC19-B)
[48]${}^{\text{b}}$ 2.07 10.55 11.51 0.54 0.76 8.82 3.59 0.799 4.9 9.9 10 23
BFH(QHC19-C) [48]${}^{\text{b}}$ 2.18 10.74 11.53 1.1 1.5 8.37 3.70 0.815 7.7
14 13 25 BFH(QHC19-D) [48]${}^{\text{b}}$ 2.28 10.86 11.54 1.5 2.1 8.12 3.92
0.833 11 18 15 26 BFH(QHC18) [48]${}^{\text{b}}$ 2.04 10.38 11.62 1.5 2.0 9.22
3.88 0.804 6.9 16 13 31 DS(CMF)-1 [48]${}^{\text{b}}$ 1.97 11.19 13.49 3.3 5.8
8.54 2.98 0.784 9.7 19 12 19 DS(CMF)-2 [48]${}^{\text{b}}$ 1.96 11.14 13.55
4.1 7.2 8.70 3.08 0.787 10 21 14 24 DS(CMF)-3 [48]${}^{\text{b}}$ 1.99 11.24
13.17 1.5 2.8 8.21 2.78 0.777 6.6 14 6.2 14 DS(CMF)-4 [48]${}^{\text{b}}$ 1.98
11.24 13.23 2.0 3.8 8.27 2.82 0.778 8.2 16 8.0 16 DS(CMF)-5
[48]${}^{\text{b}}$ 2.02 11.88 13.19 0.66 1.1 6.92 1.93 0.727 5.3 9.5 4.8 11
DS(CMF)-6 [48]${}^{\text{b}}$ 2.01 11.99 13.25 0.54 1.2 6.76 1.84 0.728 7.2 13
6.4 13 DS(CMF)-7 [48]${}^{\text{b}}$ 2.02 11.89 13.19 0.66 1.1 6.92 1.92 0.727
5.5 10 4.9 11 DS(CMF)-8 [48]${}^{\text{b}}$ 2.02 11.89 13.25 0.30 1.3 6.92
1.92 0.727 6.6 13 5.8 13 JJ(VQCD)intermediate
[48]${}^{\text{b}}$${}^{\text{c}}$ 2.15 11.80 12.15 0.97 1.9 6.12 1.92 0.689
1.1 1.6 16 49 JJ(VQCD)soft [48]${}^{\text{b}}$${}^{\text{c}}$ 2.03 11.85 12.06
1.7 2.8 5.48 1.48 0.648 3.7 6.4 19 51 JJ(VQCD)stiff
[48]${}^{\text{b}}$${}^{\text{c}}$ 2.33 11.86 12.32 0.19 0.30 6.81 2.68 0.763
5.4 12 14 44 KBH(QHC21_A) [48]${}^{\text{b}}$ 2.18 11.24 12.27 0.20 0.33 7.83
3.15 0.798 3.5 12 7.9 21 KBH(QHC21_AT) [48]${}^{\text{b}}$ 2.13 10.87 11.64
0.26 0.44 8.25 3.35 0.798 4.3 12 10 26 KBH(QHC21_B) [48]${}^{\text{b}}$ 2.25
11.41 12.28 0.37 0.59 7.49 3.08 0.803 5.0 15 9.5 24 KBH(QHC21_BT)
[48]${}^{\text{b}}$ 2.20 11.07 11.65 0.14 0.23 7.86 3.25 0.803 5.7 15 12 28
KBH(QHC21_C) [48]${}^{\text{b}}$ 2.31 11.56 12.28 0.60 0.88 7.22 3.03 0.808
6.5 18 11 26 KBH(QHC21_CT) [48]${}^{\text{b}}$ 2.26 11.13 11.60 0.54 0.77 7.68
3.31 0.812 8.1 17 14 30 KBH(QHC21_D) [48]${}^{\text{b}}$ 2.36 11.69 12.27 0.82
1.2 6.99 3.01 0.813 8.2 20 13 28 KBH(QHC21_DT) [48]${}^{\text{b}}$ 2.32 11.28
11.62 0.70 0.99 7.41 3.27 0.817 9.7 20 15 32 OOS(DD2_FRG)(2)
[48]${}^{\text{b}}$ 2.05 12.53 13.08 1.4 1.9 5.92 1.36 0.551 10 31 15 36
OOS(DD2_FRG)(2+1) [48]${}^{\text{b}}$ 1.84 12.75 13.03 2.1 2.7 5.33 0.897
0.466 11 31 26 56 OOS(DD2-FRG)vec(2) [48]${}^{\text{b}}$ 2.12 12.69 13.09 1.9
2.8 5.67 1.32 0.559 17 39 17 38 OOS(DD2-FRG)vec(2+1) [48]${}^{\text{b}}$ 1.93
12.98 13.05 2.8 3.8 4.99 0.844 0.467 21 53 31 61

##  Supplementary B Data selection

Our analysis includes several $M$ and $R$ observations of individual sources.
Clearly, each data point is obtained under certain model assumptions, and as
such, it is model-dependent. Following Refs. [54, 31] and others we perform
inferences on a different subsets of the data. Ref [55] classifies the various
astrophysical measurements by the level of their model dependence. We are not
attempting to perform a similar assessment. Nevertheless, the radio timing
mass measurements and the gravitational wave data can be considered the most
robust ones. The NICER results on radii depend strongly on the complicated
emission geometry, however, they are obtained independently by two research
groups. The CasA result is based on the specific carbon atmosphere NS thermal
emission model, however being consistent with other measurements it does not
impose significant constraints. The three remaining data sets, detailed below,
are the most constraining. The ‘spiders’ data strongly constrain
$M_{\mathrm{TOV}}$ from below; at the same time, they are based on the
complicated model for optical emission of the companion, leaving room for
serious modifications. The kilonova constraints of Ref. [24] are the only ones
that limit $M_{\mathrm{TOV}}$ from above while they rely on a particular
interpretation of the collapsar evolution. See, for instance, [56] and refs.
therein for alternative scenarios. The crucial constraints on NS radii come
from the inclusion of the data on bursting sources, mainly from 4U 1702-429.
We discuss elsewhere the implication of using different combinations of data.

We do not consider other thermally emitting isolated NSs [57] since they
either yield much more significant uncertainties of $M$ and $R$, or the
inferred mass range is beyond the range of validity of our approximations (as
for XMMU J1732–344 [58] with $M=0.8\pm 0.2\,$$M_{\odot}$; see also [59]). We
also do not consider observations of LMXBs in quiescence, which can provide
mass and radii information, since the results are currently inferior to other
NS radii constraints [60, 61]. Many other NS mass measurements exist in binary
systems (see, e.g., [62]), but they are less constraining than those we use
due to lower masses. Recently, a binary system of a pulsar and a compact
object of mass $2.35^{+0.20}_{-0.18}$ was discovered [63]. The nature of a
compact component is unknown; it can be either a black hole or a NS. In the
latter case, this mass measurement is consistent with the data based on the
spider pulsars, and we checked that it does not provide an additional
significant constraint. Because of that, and since the source can be a black
hole, we don’t include it in our fit.

##  Supplementary C Methods

We infer constraints on the model parameters using the Bayesian approach,
e.g., [26, 27, 30, 31, 32]. The construction is standard, and we briefly
describe the key points here. We are interested in the posterior probability
$\mathcal{P}(\bm{\theta}|\mathrm{data};\,\mathcal{M})$ of the model parameters
$\theta$ given the observed data and the model $\mathcal{M}$ which in our case
is the $M-R$ curve parameterization described in the main text. According to
the Bayes rule, the posterior probability of the parameters given data can be
expressed as

$\mathcal{P}(\bm{\theta}|\mathrm{data};\,\mathcal{M})=\frac{\mathcal{L}(\mathrm{data}|\bm{\theta};\
\mathcal{M})\pi(\bm{\theta}|\mathcal{M})}{p(\mathrm{data}|\mathcal{M})},$ (5)

where vector of model parameters assuming $N$ pairs of mass and radii
measurements is

$\bm{\theta}=\left\\{M_{i},\,R_{i}\right\\}_{i=1}^{N}\bigcup\bm{\theta}_{\mathrm{EoS}},\quad\bm{\theta}_{\mathrm{EoS}}=(P_{\text{{tov}}},\,\rho_{\text{{tov}}},\,r_{1/2}).$
(6)

The inference on individual $\\{M_{i},\,R_{i}\\}$ is updated given the
structural information contained in the model $\mathcal{M}$ therefore, they
are treated as the model parameters as well as the EoS parameters
$\bm{\theta}_{\mathrm{EoS}}$. In Eq. (LABEL:eq:Bayes_general),
$\mathcal{L}(\mathrm{data}|\bm{\theta};\ \mathcal{M})$ is the observed data
likelihood, $\pi(\bm{\theta}|\mathcal{M})$ is the prior probability density of
$\bm{\theta}$ (given the model) and $p(\mathrm{data}|\mathcal{M})$ is a
normalization factor (or the evidence – the prior probability of data given
the model). The likelihood of the data for each source depends only on its
mass and radius. Thus one can express $\mathcal{L}(\mathrm{data}|\bm{\theta};\
\mathcal{M})=\prod_{i=1}^{N}\mathcal{L}_{i}(\mathrm{data_{i}}|M_{i},\,R_{i})$.
Then Eq. (LABEL:eq:Bayes_general) becomes

$\mathcal{P}(\bm{\theta}|\mathrm{data};\,\mathcal{M})=\frac{\pi(\bm{\theta}_{\mathrm{EoS}}|\mathcal{M})}{p(\mathrm{data}|\mathcal{M})}\prod_{i}\mathcal{L}_{i}(\mathrm{data}_{i}|M_{i},\,R_{i})\mathcal{P}(M_{i},\,R_{i}|\bm{\theta}_{\mathrm{EoS}},\,\mathcal{M}),$
(7)

where $\mathcal{P}(M_{i},\,R_{i}|\bm{\theta}_{\mathrm{EoS}},\,\mathcal{M})$ is
the probability of a particular source having mass $M_{i}$ and radius $R_{i}$
given the EoS parameter vector $\bm{\theta}_{\mathrm{EoS}}$, and
$\pi(\bm{\theta}_{\mathrm{EoS}}|\mathcal{M})$ is the prior probability of the
EoS parameters. Both terms under the product sign require specific
consideration. The likelihood for the individual source can be constructed
from scratch based on the data generation statistical model, basically
repeating the analysis of the original works on the source, which is usually
unfeasible in practice. Alternatively, the likelihood needs to be deduced from
the published analysis results which are usually available in the form of the
posterior probabilities $\mathcal{P}(M_{i},R_{i}|\mathrm{data}_{i})$ if
Bayesian approach is used or in the form of confidence regions or confidence
distributions [64] if the frequentist approach is adopted. In rare cases, the
likelihood function itself is reported, like in the case of the toast
likelihood for the GW170817 [65], see below. When Bayesian approach was used
for the $i$th source in the reference study, the likelihood can be expressed
via the Bayes rule

$\mathcal{L}_{i}(\mathrm{data_{i}}|M_{i},\,R_{i})=\frac{p(\mathrm{data}_{i})\mathcal{P}(M_{i},\,R_{i}|\mathrm{data})}{\pi_{\mathrm{ref}}(M_{i},R_{i})}.$
(8)

Therefore, either the prior $\pi_{\mathrm{ref}}(M_{i},\,R_{i})$ of the
reference study needs to be explicitly specified, or the likelihood function
can be approximated by the posterior probability itself provided the prior is
sufficiently flat in the region of support of the final posterior probability
(so it can be treated as a constant during marginalization).

The conditional probability
$\mathcal{P}(M_{i},\,R_{i}|\bm{\theta}_{\mathrm{EoS}})$ is also a subtle term
since it involves conditioning over a null set. Depending on a choice of a
parameterization, the relationship between $M_{i},\,R_{i}$ and
$\bm{\theta}_{\mathrm{Eos}}$ can be expressed in different ways, leading, in
principle, to different inferences. For instance, specification of the model
as the relation $R=R_{i}(M_{i},\,\bm{\theta}_{\mathrm{EoS}})$ supplemented
with prior distributions on $\bm{\theta}_{\mathrm{EoS}}$ and $M_{i}$ induces a
specific prior distribution on $R_{i}$ which is not expected to be uniform or
in any other ways non-informative alike. Adopting a different parameterization
of the $M-R$ curve would result in different $M_{i},\,R_{i}$ prior
distribution given the structural model $\mathcal{M}$ and, consequently, in
different posterior distributions. This difficulty is a manifestation of the
Borel-Kolmogorov paradox, e.g., [66]. One of the strategies to overcome it is
in some way to combine, or ‘meld’ the induced priors with the prior
distributions on $M_{i},\,R_{i}$ specified before the one is informed about
the model $\mathcal{M}$ [67] putting in a sense a selection of a specific
parametrization in Bayesian context. In the astrophysical context, the
dependence of the inference on low-dimensional relations embedded in higher-
dimension space on the parameterization of this space and a related problem of
choice of priors on the parameters is discussed, e.g., in [68, 69].

Having said this we nevertheless follow the previous studies and formulate the
$M-R$ relationship in a form $R=R_{i}(M_{i},\,\bm{\theta}_{\mathrm{EoS}})$
with uniform prior on $M_{i}$ in the range
$M_{\mathrm{min}}...M_{\text{{tov}}}(\bm{\theta}_{\mathrm{EoS}})$, where
$M_{\mathrm{min}}=1M_{\odot}$ in our case. In a sense, a specific choice of
the parametrization can be regarded as a part of the structure model. It is
believed that the data likelihood mainly determines the posteriors; therefore,
variations in a choice of the parametrization would not produce much harm. An
alternative would be to select some physics-motivated prior on $M_{i}$, e.g.,
equate it to the observed NS mass distribution. Given the observed NS
distribution and assuming that the structural relationship $\mathcal{M}$ is
the correct one, we then a priori expect a specific NS radii distribution.
However, the sample of the NSs with measured masses as well as the sample of
NSs for which $M-R$ measurements exist are not fair samples from the NS
population due to unquantified selection effects. In the absence of a better
solution, the conditional probability for mass and radius of $i$th source
given EoS parameters is then

$\mathcal{P}(M_{i},\,R_{i}|\bm{\theta}_{\mathrm{EoS}})=\mathcal{P}(R_{i}|M_{i},\,\bm{\theta}_{\mathrm{EoS}})\mathcal{P}(M_{i}|\bm{\theta}_{\mathrm{EoS}})=\delta\left(R_{i}-R(M_{i},\bm{\theta}_{\mathrm{EoS}})\right)\frac{\Theta(M_{\text{{tov}}}(\bm{\theta}_{\mathrm{EoS}})-M_{i})\Theta(M_{i}-M_{\mathrm{min}})}{M_{\text{{tov}}}(\bm{\theta}_{\mathrm{EoS}})-M_{\mathrm{min}}},$
(9)

with $\Theta(x)$ for the Heaviside step function. Inserting Eqs.
(LABEL:eq:Likei) and (LABEL:eq:Peos) in Eq. (LABEL:eq:Post_fin) we can infer
the posterior probability densities on the parameters $\bm{\theta}$. For NICER
and CasA data, the samples from posterior distributions on $M_{i}$ and $R_{i}$
are available [38, 35, 39] and we use the kernel density estimation (KDE) for
$\mathcal{P}(M_{i},R_{i}|\mathrm{data}_{i})$ in Eq. (LABEL:eq:Likei). For the
GW data we also use 4D KDE of the
$M_{\mathrm{GW}1},\,R_{\mathrm{GW}1},\,M_{\mathrm{GW}2},\,R_{\mathrm{GW}2}$
samples from Ref. [32]. In the latter case it is also possible to use directly
the marginalized (overall nuisance parameters of GW wavefront modeling)
likelihood
$L_{\mathrm{GW}}(\mathrm{data}|M_{\mathrm{GW}1},\,R_{\mathrm{GW}1},\,M_{\mathrm{GW}2},\,R_{\mathrm{GW}2})$
provided by Ref. [65]. We checked, and both approaches give similar results.

When only the mass $M_{i}$ of a given source is measured, the sole information
one gets is that $M_{\text{{tov}}}>M_{i}$. Then Eq. (LABEL:eq:Peos) can be
marginalized over $R_{i}$ giving

$\mathcal{P}(M_{i},|\bm{\theta}_{\mathrm{EoS}})=\mathcal{P}(M_{i}|M_{\text{{tov}}})=\frac{\Theta(M_{\text{{tov}}}(\bm{\theta}_{\mathrm{EoS}})-M_{i})\Theta(M_{i}-M_{\mathrm{min}})}{M_{\text{{tov}}}(\bm{\theta}_{\mathrm{EoS}})-M_{\mathrm{min}}}.$
(10)

The likelihood of such observation, marginalized over $M_{i}$, is then

$\mathcal{L}(\mathrm{data}_{i}|\bm{\theta}_{\mathrm{EoS}})=\mathcal{L}(\mathrm{data_{i}}|M_{\mathrm{TOV}})=\int\mathrm{d}M_{i}\,\mathcal{L}(\mathrm{data}_{i}|M_{i})\mathcal{P}(M_{i}|M_{\text{{tov}}})\propto\mathrm{CDF}_{i}(M_{\text{{tov}}})$
(11)

if we approximate likelihood $\mathcal{L}(\mathrm{data}_{i}|M_{i})$ as the
posterior distribution of $M_{i}$ as above. Here
$\mathrm{CDF}_{i}(M_{\text{{tov}}})$ is the cumulative posterior distribution
for $M_{i}$, i.e., the posterior probability that $M_{i}<M_{\text{{tov}}}$
given the data. When the inference is summarized in the form of a frequentists
confidence interval, we assume the normal (Gaussian) likelihood [64], which
leads to the same expression with $\mathrm{CDF}_{i}(M_{\text{{tov}}})$ being
the normal distribution cumulative distribution function (namely, the error
function). Likewise, for the kilonova constraint, where $M_{\text{{tov}}}$ is
limited from above, the complementary cumulative distribution function is used

$\mathcal{P}_{\mathrm{kilonova}}(\mathrm{data}|M_{\mathrm{TOV}})\propto
1-\mathrm{CDF}_{\mathrm{kilonova}}(M_{\text{{tov}}}).$ (12)

Apart from the causality condition (see Sec. III in the main text) and priors
on $M_{i}$ described above, we use the following weak non-informative priors
on the rest of model parameters: $1<r_{1/2}<2$,
$0.1<P_{\text{{tov}}}/\rho_{0}<100$, and
$0.1<\rho_{\text{{tov}}}/\rho_{0}<100$. Additionally, we set the conservative
limit $M_{\text{{tov}}}<3\,$$M_{\odot}$. To sample from posteriors and obtain
credible regions for the model parameters we employ affine-invariant Markov
Chain Monte Carlo (MCMC) sampler emcee [70]. Given the number of model
parameters, 32 walkers were used. The convergence of the MCMC chain was
controlled by the autocorrelation time following the emcee manual
recommendations. The chain is considered to be converged when the
autocorrelation time estimate is stable and the number of MCMC steps is more
than 50 times the autocorrelation time. A typical number of steps in the
converged chain was about 20000, and the first half of the chain was burned
in. All credible intervals and regions given in the tables and shown in the
plots correspond to the highest posterior density intervals.

$P_{\mathrm{TOV}}/(\rho_{0}c^{2})$ | $\rho_{\mathrm{TOV}}/\rho_{0}$ | $r_{1/2}$ | $M_{\mathrm{TOV}}$ | $R_{\mathrm{TOV}}$ | $R_{1/2}$
---|---|---|---|---|---
|  |  | $[M_{\odot}]$ | [km] | [km]
$3.5^{+1.0}_{-0.4}$ | $7.7^{+0.6}_{-0.3}$ | $1\dots 1.10$ | $2.28^{+0.05}_{-0.06}$ | $11.19^{+0.39}_{-0.34}$ | $11.97^{+0.37}_{-0.31}$

The model parameters are summarized in Table LABEL:tab:allmod_params and
marginalized 1D and 2D posterior distributions are presented in Fis.
LABEL:fig:tri1. We show the fit parameters
$P_{\mathrm{TOV}}\,\rho_{\mathrm{TOV}},\,r_{1/2}$ as well as the derived
parameters $M_{\mathrm{TOV}},\,R_{\mathrm{TOV}},\,R_{1/2}$. The uncertainties
in Table LABEL:tab:allmod_params correspond to 68% highest posterior density
intervals for all parameters except $r_{1/2}$. Since for the latter the lower
bound of such interval is 1 in our case (see Fig. LABEL:fig:tri1), we show the
range which corresponds to the 68% quantile of the marginalized $r_{1/2}$
posterior distribution.

##  Supplementary D Individual sources

| GW170817-1 | GW170817-2 | 0740 | CasA | 0030 | 1702 | 1724 | 1810
---|---|---|---|---|---|---|---|---
$M$ [$M_{\odot}$] | $1.42^{+0.11}_{-0.05}$ | $1.30^{+0.04}_{-0.09}$ | $2.05^{+0.08}_{-0.05}$ | $1.58^{+0.13}_{-0.12}$ | $1.37^{+0.06}_{-0.11}$ | $1.74^{+0.17}_{-0.12}$ | $1.42^{+0.33}_{-0.25}$ | $1.39^{+0.37}_{-0.24}$
$R$ [km] | $12.20^{+0.21}_{-0.33}$ | $12.25^{+0.25}_{-0.25}$ | $12.20^{+0.21}_{-0.33}$ | $12.25^{+0.25}_{-0.25}$ | $12.13^{+0.32}_{-0.30}$ | $12.30^{+0.22}_{-0.25}$ | $12.22^{+0.28}_{-0.34}$ | $12.22^{+0.29}_{-0.33}$

Bayesian procedure allows to update inferences on masses and radii of
individual sources, given they all have the same EoS in their interiors. The
corresponding 68% credible intervals are given in Table
LABEL:tab:ind_masses_radii. We also illustrate this in Fig.
LABEL:fig:individuals. In each plot, thin lines correspond to the initial
measurements, while thick lines are the updated results. One can notice the
shrinkage of the credible intervals on $R$ for each source (the least
pronounced for 4U J1702–429).