# I – Love – $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $: Approximately
universal relations for the average neutron star stiffness

Jayana A. Saes<EMAIL_ADDRESS>Illinois Center for Advanced Studies of
the Universe, Department of Physics, University of Illinois at Urbana-
Champaign, Urbana, IL 61801, USA  Raissa F. P. Mendes<EMAIL_ADDRESS>Instituto de Física, Universidade Federal Fluminense, Niterói, Rio de Janeiro,
24210-346, Brazil  CBPF - Centro Brasileiro de Pesquisas Físicas, Rio de
Janeiro, 22290-180, Brazil  Nicolás Yunes<EMAIL_ADDRESS>Illinois Center
for Advanced Studies of the Universe, Department of Physics, University of
Illinois at Urbana-Champaign, Urbana, IL 61801, USA

(September 3, 2024)

###### Abstract

The accurate observations of neutron stars have deepened our knowledge of both
general relativity and the properties of nuclear physics at large densities.
Relating observations to the microphysics that govern these stars can
sometimes be aided by approximate universal relations. One such relation
connects the ratio of the central pressure to the central energy density and
the compactness of the star, and it has been found to be insensitive to
realistic models for the equation of state to a $\sim 10\%$ level. In this
paper, we clarify the meaning of the microscopic quantity appearing in this
relation, which is reinterpreted as the average of the speed of sound squared
in the interior of a star, $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $.
The physical origin of the quasi-universality of the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ -C$ relation is then investigated. Making
use of post-Minkowskian expansions, we find it to be linked to the Newtonian
limit of the structure equations, as well as to the fact that the equations of
state that describe NSs are relatively stiff. The same post-Minkowskian
approach is also applied to the relations between $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $, the moment of inertia, and the tidal
deformability of a neutron star, arriving at similar conclusions.

##  I Introduction

Neutron stars (NSs) are highly compact astrophysical objects, in which matter
exists at extremely dense levels, with their maximum density possibly reaching
up to $\sim$10 times the saturation density $\rho_{sat}=2.7\times
10^{14}$g/cm3. This unique environment renders a NS an ideal laboratory for
fundamental physics, encompassing tests of general relativity (GR) and
investigations into the equation of state (EoS) of supra-nuclear matter. To
realize this potential, accurate measurements of NS observables are needed.
Recent gravitational wave (GW) detections of binary NS mergers [1, 2], radio
observations of pulsars [3, 4, 5], and pulse-profiling of X-ray data from the
Neutron-Star-Interior Composition Explorer (NICER) [6, 7, 8, 9], are beginning
to provide the first accurate measurements of NS radii, masses, and tidal
deformabilities. Next-generation GW observatories such as the Einstein
Telescope [10] and the Cosmic Explorer [11], as well as future X-ray missions,
will allow an increase in abundance and precision of such measurements.

An important tool for the interpretation of NS observational data and its
linkage to the microphysics governing the interior of these stars lies in
approximate universal relations, i.e., EoS-insensitive relations between NS
properties. One example is the I-Love-Q relations [12], which connect (a
dimensionless version of) the moment of inertia, tidal deformability, and
quadrupole moment to each other, and can be used to break degeneracies in GW
parameter estimation. Another example is the binary Love relations [13, 14],
which connect (a dimensionless version of) the tidal deformabilities of two
NSs in a binary to each other and can enable the measurement of individual
tidal deformabilities from GW data. For a review of approximate universal
relations, we refer the reader to Ref. [15].

Approximate universal relations between macroscopic observables and
hydrodynamic quantities are also known to exist. One example is the relation
between the pressure at densities $(1$–$2)\,\rho_{sat}$ and the radius and
tidal deformability of a NS with a mass of $(1$–$1.4)M_{\odot}$ [16, 17],
which has been used to constrain the EoS space at these densities using
observational data [18, 19]. In a recent study [20], a novel approximate
universal relation of this type was established. This relation links the ratio
between the central pressure and the central energy density, denoted as
$p_{c}/\epsilon_{c}$, to the compactness $C$, dimensionless moment of inertia
$\bar{I}$, and dimensionless tidal deformability $\Lambda$. This relation was
shown to hold for a large range of realistic and parameterized EoSs commonly
used to describe NSs, which we illustrate in Fig. 1 for a set of tabulated
EoSs. Through these relations, a measure of $C$, $\Lambda$ or $\bar{I}$ would
therefore yield an inference of $p_{c}/\epsilon_{c}$.

Figure 1: Relation between the average speed of sound squared (or,
equivalently, the ratio of central pressure to central energy density,
$\langle c_{s}^{2}\rangle\leavevmode\nobreak\ =p_{c}/\epsilon_{c}$) and the
compactness ($C$), tidal deformability ($\Lambda$), and dimensionless moment
of inertia ($\bar{I}$) of NSs obeying a set of realistic EoSs. In particular,
KDE0v [21], KDE0v1 [21], Rs [22], SK255 [21], SK272 [21], SKI2-6 [23], SkMp
[24], SkOp [23], SLY2 [25], SLY230a [25], and SLY9 [25] EoSs [26, 27], with
tables obtained from the CompOSE repository. Curves for an incompressible
fluid and the Tolman VII solution are also shown.

The perhaps surprising insensitivity of this relation to the choice of EoS,
despite it connecting hydrodynamic quantities ($p_{c}/\epsilon_{c}$) to
macroscopic observables, raises questions about the underlying reasons for its
existence. Previous research, e.g. [12, 28, 29, 30], has been dedicated to
identifying the causes of the emergence of insensitive relations in NSs,
particularly regarding the I-Love-Q one. In these studies, several aspects of
the approximate universal relation were analyzed, such as analytic
approximations for the observables, the contributions of different segments of
the NS interior and EoS to the approximate universality, and which assumptions
are more relevant to the insensitivity, among others. Inspired by the
strategies adopted in these works, we here study the physical origins of the
quasi-universality of the $p_{c}/\epsilon_{c}$ — $C/\Lambda/\bar{I}$ relation
presented in [20].

Figure 2: Speed of sound as a function of the energy density for two EoSs with
varying stiffness, KDE0v [27, 26, 31] and SLY2 [32]. The filled circles
indicate the central energy density and the speed of sound squared that
corresponds to a NS with compactness $C=0.2$ for each EoS.

We find that the cause of the quasi-universality is closely tied to the fact
that $p_{c}/\epsilon_{c}$ is nothing but the average speed of sound squared
inside a NS, an interpretation that we establish formally in this paper.
Figure 2 shows the speed of sound squared $c_{s}^{2}$ as a function of the
energy density $\epsilon$ for two EoSs. The filled circles in this figure
indicate the central energy density and the maximum speed of sound inside a
star with $C=0.2$. The speed of sound squared, averaged over the energy
density, is then nothing but the “area under the curve” (i.e. the shaded
regions) divided by the central energy density (because the energy density at
the surface is zero). This average, which we denote $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $, has a natural interpretation as a
measure of the average stiffness inside the NS.

The approximate universality between $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $and $C$ is then connected to the fact
that two NSs with different EoSs but the same compactness will typically have
different central energy densities. As we can see from Fig. 2, as one
considers an EoS with a smaller stiffness, one naturally decreases the rate at
which the speed of sound increases as a function of energy density in the NS
interior. This, then, implies that, to obtain the same compactness as a NS
with a stiffer EoS, the central energy density must be higher. This naturally
increases the area under the $c_{s}^{2}$-$\epsilon$ curve, in such a way that
the ratio of this new area to the new central energy density is roughly
unchanged.

To go beyond the qualitative explanation presented above, we study these
approximately universal relations analytically, both in the Newtonian limit
and in post-Minkowskian expansions. In the Newtonian limit, we find that the
same relation, $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ =C/2$, applies
for various EoSs (i.e. an incompressible fluid, an $n=1$ polytrope, and a
Tolman VII density profile [33]), which, despite their simplicity, are
relevant to the description of NSs. We interpolate between these models by
considering two one-parameter families of EoSs, namely, a family of polytropes
and a family of generalized Tolman VII profiles, the latter of which we
introduce in this paper. We show that, in the Newtonian limit, the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$ relation is remarkably flat over
the range of model parameters most relevant to the description of NSs, with an
EoS-dependence of only $\sim 3\%$ and $\sim 6\%$ for polytropic EoSs and the
generalized Tolman VII model, respectively (cf. Fig. 4 in Sec. III.2.1).

We confirm and extend the above conclusions by considering NSs in the post-
Minkowskian expansion by considering solutions as a power series in the NS
compactness. We find that the $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $–
$C$ relation becomes increasingly sensitive to the EoS as one includes higher
powers of compactness (i.e. higher order GR contributions). The EoS-dependence
increases to $\sim 9\%$ for both the polytrope family and the generalized
Tolman VII family, when considering up to third-order contributions in
compactness for all stable NSs with these EoSs (cf. Figs. 5 and 6 in Sec.
III.3.1). We assess to which point these semi-analytic results for simple,
one-parameter families of EoSs can be informative of realistic EoSs, and
conclude that this is the case as long as the star is sufficiently compact
($C\gtrsim 0.1$).

With an understanding of the $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $–
$C$ relation under our belts, we then move on to consider the
$I$-Love-$\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $relation. We carry out
a similar study in the Newtonian limit and through post-Minkowskian
expansions. Unlike in $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$
case, we find that the $I$-Love-$\langle c_{s}^{2}\rangle\leavevmode\nobreak\
$relation presents the same degree of universality in the Newtonian limit (cf.
Figs. 9 and 10 in Sec. IV.1.2) and in the presence of increasing GR
contributions (cf. Figs. 11 and 12 in Sec. IV.2). For example, we find that
the sensitivity of the $I$-Love-$\langle c_{s}^{2}\rangle\leavevmode\nobreak\
$relation to variations of the generalized Tolman VII family remains at the
$\sim 10-11\%$ level to all orders in compactness.

The post-Minkowskian expansions performed here not only enable one to access
the influence of relativistic effects in the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C/\bar{I}/\Lambda$ relations, but it
also allows for the construction of analytic (e.g. Padé) approximants that can
more efficiently (i.e., with fewer parameters) describe these relations (see
e.g. [34]). For instance, for the realistic EoSs shown in Fig. 1, the Padé-
inspired fitting function

$\langle c_{s}^{2}\rangle\leavevmode\nobreak\ =C\frac{1+p_{1}C}{q_{0}+q_{1}C}$
(1)

with parameters $(p_{1},q_{0},q_{1})\approx(-0.73,1.84,-4.62)$ provides a
simple description that captures well the overall behavior of the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$ relation for all stable NSs.

The remainder of this paper presents the details that support the conclusions
summarized above and it is organized as follows. In Sec. II, we elucidate the
physical significance of $p_{c}/\epsilon_{c}$ as the average speed of sound
squared. In Secs. III and IV, we delve into the origins of the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$ relation and $I$-Love-$\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $relations, correspondingly, following a
post-Minkowksian approach. Section V summarizes our conclusions.

##  II Re-interpreting $p_{c}/\epsilon_{c}$ as a measure of the average NS
stiffness

Let us begin by providing a new interpretation for the ratio
$p_{c}/\epsilon_{c}$ between the central pressure and central energy density.
This ratio is important because a quasi-universal relation between
$p_{c}/\epsilon_{c}$ and NS observables, such as the NS compactness $C$, was
established in [20]. As argued in that reference, this ratio can be physically
interpreted as a measure of the mean stiffness of nuclear matter inside a NS,
indicative of how rapidly the pressure increases with the energy density111The
ratio $p_{c}/\epsilon_{c}$ is also indicative of the conformality at the
center of the star via the trace anomaly $\Delta_{\rm
conf}=1/3-p_{c}/\epsilon_{c}$ [35].. Given this, let us use the fact that the
energy density is zero when the pressure is null, i.e., $\epsilon(p=0)=0$, so
that then

$\displaystyle\frac{p_{c}}{\epsilon_{c}}=\frac{1}{\epsilon_{c}}\int_{0}^{\epsilon_{c}}\frac{dp}{d\epsilon}d\epsilon=\frac{1}{c^{2}\epsilon_{c}}\int_{0}^{\epsilon_{c}}c_{s}^{2}d\epsilon=:\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \,.$ (2)

The quantity $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $is the
dimensionless speed of sound squared averaged over the range of energy
densities present inside the NS. This expression reaffirms the interpretation
provided in [20] because the speed of sound is a typical indicator of the
stiffness of an EoS, and thus, its average would represent the mean stiffness.

A visual representation of this average is shown in Fig. 2, where we present
the function $c_{s}^{2}(\epsilon)$ for two representative EoSs. The average
$\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $associated with a specific
compactness $C$ is computed by calculating the area under the curve of
$c_{s}^{2}(\epsilon)$ up to the central energy density corresponding to the
star with compactness $C$, and then, dividing this area by the central energy
density. For two different EoSs to yield the same compactness, the softer EoS,
characterized by a lower speed of sound, must have a higher central energy
density. This accounts for the fact that, even though different EoSs predict
different speeds of sound in the center of the star (with well-defined
compactness) and different central energy densities, they yield a similar
$\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $.

The physical interpretation of $\langle c_{s}^{2}\rangle\leavevmode\nobreak\
$as a stiffness measure is visually represented in Fig. 3, which shows various
EoSs and NS properties for three representative EoSs with different stiffness.
The top panels describe the EoS itself (i.e. $p(\epsilon)$ on the left and
$c_{s}^{2}(\epsilon)$ on the right), while the bottom panels describe
properties of equilibrium sequences of NSs for each EoS (i.e. the M–R curve on
the left and the $\langle c_{s}^{2}\rangle\leavevmode\nobreak\
$–$\epsilon_{c}$ curve on the right). The stiffer the EoS, the more rapidly
$p$ increases with $\epsilon$, and thus, the larger the slope of the
$p(\epsilon)$ curve and the larger the speed of sound. Therefore, the top left
and right panels tell us that the SK272 is the stiffest EoS considered,
followed by SLY2, and then by SkOp. A stiffer EoS typically implies a larger
pressure force to balance gravity, leading to NSs with larger radii and mass,
as one can see in the bottom, left panel. On the other hand, a stiffer EoS
typically implies a larger averaged speed of sound squared, as one can see on
the bottom, right panel.

Figure 3: Curves for three realistic EoSs, SkOp [26, 27, 23], SLY2 [32] and
SK272 [27, 21, 26] which are ordered with increasing stiffness. The legend and
colors of the top left panel are preserved for all plots. Top left panel:
Pressure versus energy density curves. Top right panel: Speed of sound versus
energy density for the three EoSs. Bottom left panel: Mass-radius curves found
by integrating the structure equations for NSs varying the central energy
density up to the one that corresponds to maximum mass. Bottom right panel:
$\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $versus central energy density
for each of the stars presented in the bottom left panel. All of the curves
are cut off at the values of energy density that generate the star with
maximum mass.

##  III $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$ relation

With this clearer physical interpretation in mind, we now proceed to the
investigation of the origins of the quasi-universal relation between $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $and the NS compactness, $C\equiv
GM/Rc^{2}$. The latter and $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $are
known to display a quasi-universal relation, first presented in [20], with an
EoS-dependence of $\mathcal{O}(10\%)$ for realistic EoSs. To gain a deeper
understanding of the underlying causes behind this quasi-universality, in this
section, we first explore its presence within analytic solutions for the TOV
equation [36, 33]. Namely, we study the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$ relation using the incompressible,
Tolman VII and Buchdahl solutions, and then analyze the Newtonian limit of
these solutions. Using the Newtonian equations of structure, we follow this up
with a study of the relation for 2 one-parameter families of EoSs: polytropes
and a generalized version of the Tolman VII EoS. Next, we conduct a post-
Minkowskian expansion within these solutions. Finally, we contrast these
simple models with realistic EoSs and discuss in which regime the results
obtained for the former can be informative of the latter.

###  III.1 Exact Solutions in GR

Exploring analytic solutions in GR not only helps us elucidate the origin of
the EoS-insensitive relation but also inspires improved analytic formulae that
serve as better approximations of such a relation for realistic EoSs. While
several known analytic solutions to the TOV equations, not all of them are
relevant and accurate descriptions for NSs. Here, as in [16], we shall present
the $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$ relation for three
exact solutions of the TOV equation: the incompressible fluid, a Tolman VII
fluid, and a Buchdahl fluid.

For concreteness, and to further establish notation, let us review the TOV
equations in GR. The relativistic generalization of the equation of
hydrostatic equilibrium is

$\displaystyle\frac{dp}{dr}=-\frac{G}{c^{2}}(\epsilon+p)\frac{m+4\pi
r^{3}p/c^{2}}{r\left(r-\frac{2Gm}{rc^{2}}\right)}\,,$ (3)

where $m$ is the enclosed mass, $p$ is the internal pressure, $\epsilon$ is
the energy density, and $r$ is the radial coordinate. The enclosed mass
satisfies the equation:

$\displaystyle\frac{dm}{dr}=4\pi r^{2}\frac{\epsilon}{c^{2}}\,.$ (4)

These two equations form a closed system once an EoS is prescribed and their
solution requires boundary conditions. As usual, we impose regularity at the
center of the star, enforcing $m(r=0)=0$, while $p(r=0)=p_{c}$, where $p_{c}$
is the central pressure. With these conditions, the radius of the star $R$ is
defined as the radial coordinate at which the pressure vanishes, $p(r=R)=0$,
while the stellar mass is the enclosed mass evaluated at this radial
coordinate $M=m(r=R)$.

####  III.1.1 Incompressible fluid

We first examine the incompressible case, which describes a fluid with
constant density. The EoS of this fluid is given by

$\displaystyle\epsilon(p)$ $\displaystyle=\epsilon_{c}=\text{constant},$ (5)

and it results in the following analytic solutions to the TOV and the enclosed
mass equations

$\displaystyle p(r)$
$\displaystyle=\frac{\epsilon_{c}\left(R\sqrt{1-2C}-\sqrt{R^{2}-2Cr^{2}}\right)}{\sqrt{R^{2}-2Cr^{2}}-3R\sqrt{1-2C}},\quad
m(r)=\frac{4}{3}\pi r^{3}\frac{\epsilon_{c}}{c^{2}},$ (6)

where the central energy density $\epsilon_{c}$ can be expressed in terms of
the stellar mass and radius by

$\displaystyle\epsilon_{c}=\frac{3Mc^{2}}{4\pi R^{3}}.$ (7)

Evaluating $p(r=0)$, we obtain

$\displaystyle p_{c}=\epsilon_{c}\frac{1-\sqrt{1-2C}}{3\sqrt{1-2C}-1},$ (8)

which leads to the simple analytic relation222Although an incompressible fluid
does not satisfy the condition $\epsilon(p=0)=0$ required for Eq. (2) to hold,
and has a diverging speed of sound, throughout the paper we adhere to the
definition $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ \equiv
p_{c}/\epsilon_{c}$.

$\displaystyle\langle c_{s}^{2}\rangle\leavevmode\nobreak\
=\frac{p_{c}}{\epsilon_{c}}$
$\displaystyle=\frac{1-\sqrt{1-2C}}{3\sqrt{1-2C}-1}.$ (9)

This can be expanded for small values of compactness to find

$\displaystyle\langle c_{s}^{2}\rangle\leavevmode\nobreak\
\approx\frac{C}{2}+C^{2}+\mathcal{O}(C^{3}).$ (10)

The above solutions present some interesting, albeit well-known behavior.
First, the averaged speed of sound diverges when $C\to 4/9$, which corresponds
to the Buchdahl limit, i.e. an equilibrium sequence of incompressible
solutions to the TOV equations ceases to exist at this value of compactness.
Second, the average speed of sound squared exceeds unity when $C>3/8$. For
extremely non-relativistic stars, the average speed of sound tends to zero,
and, for weakly-gravitating stars, it scales with compactness. The radius of
convergence of the post-Minkowskian expansion of this series is simply the
Buchdahl limit, as one can confirm by inspection or through the ratio test.

####  III.1.2 Tolman VII fluid

The Tolman VII solution, introduced in [33], is characterized by the density
profile

$\displaystyle\epsilon(r)$
$\displaystyle=\epsilon_{c}\left[1-\left(\frac{r}{R}\right)^{2}\right].$ (11)

To simplify the resulting analytic expressions presented in [34], we introduce
the dimensionless variable $x=r/R$. The mass aspect function and pressure can
be expressed as follows:

$\displaystyle m(x)$
$\displaystyle=\frac{4\pi\epsilon_{c}}{R^{3}}\left(\frac{x^{3}}{3}-\frac{x^{5}}{5}\right),$
(12) $\displaystyle p(x)$ $\displaystyle=\frac{c^{4}}{4\pi
GR^{2}}\left[\sqrt{3Ce^{-\lambda}}\tan\phi-\frac{C}{2}(5-3x^{2})\right],$ (13)

where $\epsilon_{c}$ can be expressed in terms of $R$ and the stellar mass
$M$, or, alternatively, the compactness $C$, as

$\displaystyle\epsilon_{c}=\frac{15Mc^{2}}{8\pi R^{3}}=\frac{15c^{4}C}{8\pi
GR^{2}}.$ (14)

The functions $\lambda(x)$ and $\phi(x)$ are given by

$\displaystyle e^{-\lambda(x)}$ $\displaystyle=$ $\displaystyle
1-Cx^{2}\left(5-3x^{2}\right),$ (15) $\displaystyle\phi(x)$ $\displaystyle=$
$\displaystyle
k-\frac{1}{2}\log\left(x^{2}-\frac{5}{6}+\sqrt{\frac{e^{-\lambda}}{3C}}\right),$
(16)

where

$\displaystyle
k=\arctan\sqrt{\frac{C}{3(1-2C)}}+\frac{1}{2}\log\left(\frac{1}{6}+\sqrt{\frac{1-2C}{3C}}\right).$
(17)

At the center of the star ($x=0$), we have that $\lambda=0$ and

$\displaystyle
p_{c}=\frac{\epsilon_{c}}{15}\left(\frac{2\sqrt{3}\tan\phi_{c}}{\sqrt{C}}-5\right),$
(18)

where $\phi_{c}=\phi(0)=k-\frac{1}{2}\log\left(1/\sqrt{3C}-5/6\right)$.

All of the above allows us to easily calculate the averaged speed of sound
squared analytically, namely

$\displaystyle\langle c_{s}^{2}\rangle\leavevmode\nobreak\
=\frac{1}{15}\left(\frac{2\sqrt{3}\tan\phi_{c}}{\sqrt{C}}-5\right).$ (19)

This expression can be Taylor expanded in small $C$ to obtain

$\displaystyle\langle c_{s}^{2}\rangle\leavevmode\nobreak\
\approx\frac{C}{2}+\frac{133C^{2}}{120}+\mathcal{O}(C^{3}).$ (20)

There are a few features of Eq. (19) and its expansion that are noteworthy.
First, this result shows that the central pressure, and consequently, the
averaged squared speed of sound becomes infinite when the compactness reaches
the point where $\phi_{c}(C)=\pi/2$, which corresponds to $C\approx 0.386$.
This value of the compactness is lower than the Buchdahl limit of
incompressible case fluids, and it defines the end of the equilibrium branch
of the Tolman VII solution. However, before diverging, $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $reaches and subsequently exceeds unity
at $C\approx 0.335$, which sets the limit beyond which the fluid becomes
superluminal, on average. The Taylor series for this analytic relation
converges for $C\lesssim 0.386$, as one can check with the ratio test, which
coincides with the value beyond which the central pressure diverges.

####  III.1.3 Buchdahl fluid

The Buchdahl EoS [37] is given by

$\displaystyle\epsilon(p)=12(p_{*}p)^{1/2}-5p,$ (21)

where $p_{*}$ is a constant. This equation of state reduces, in the Newtonian
limit, to a $n=1$ polytrope, which we consider below (see Sec. III.2.1). With
this EoS, one can find an analytic solution to the TOV equation, which, as
presented in [16], is

$\displaystyle p(r)=\frac{c^{4}}{8\pi G}A^{2}u^{2}(1-2C)(1-C+u)^{-2},$ (22)

where

$\displaystyle r$ $\displaystyle=r^{\prime}(1-C+u)(1-2C)^{-1},$ (23)
$\displaystyle u$ $\displaystyle=C(Ar^{\prime})^{-1}\sin Ar^{\prime},$ (24)
$\displaystyle A^{2}$ $\displaystyle=288\pi p_{*}Gc^{-4}(1-2C)^{-1}.$ (25)

At the center of the star, $r=r^{\prime}=0$ and $u=C$; thus, the central
pressure can be expressed as

$\displaystyle p_{c}=36p_{*}C^{2},$ (26)

and, from Eq. (21), the central energy density is

$\displaystyle\epsilon_{c}=72p_{*}C\left(1-\frac{5C}{2}\right).$ (27)

With all of the above, we can now compute the averaged speed of sound square
as

$\displaystyle\langle c_{s}^{2}\rangle\leavevmode\nobreak\ =\frac{C}{2-5C},$
(28)

and for small values of compactness:

$\displaystyle\left<c_{s}^{2}\right>\approx\frac{C}{2}+\frac{5C^{2}}{4}+\mathcal{O}(C^{3}).$
(29)

Examining these expressions we first note that the central pressure, given by
Eq. (26), is regular for all real values of the compactness, meaning that,
unlike the previous solutions, this solution does not diverge for any value of
$C$. We can observe, on the other hand, that the central energy density goes
to zero when $C=2/5$, which will be the limiting factor that marks the end of
the equilibrium sequence for the Buchdahl solution. This also means that at
that value of $C$, $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $diverges.
Before that value, when $C=1/3$, the average squared speed of sound becomes
greater than one, and the fluid becomes superluminal, on average. The
expansion for low compactness values, Eq. (29), has a convergence radius of
$C=2/5$, as expected from the analytic expression, and as can also be found
via the ratio test. As a final remark, notice that, for the Buchdahl solution,
the relation between $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $and $C$,
given in Eq. (28), is exactly the ratio of two linear polynomials, which
succinctly captures both the $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ \to
C/2$ limit as $C\to 0$ and the pole at $C=2/5$. This same structure, but
suitably generalized, as in Eq. (1), can also provide useful fitting formulae
for the case of realistic EoSs.

####  III.1.4 Comparing solutions

Let us now compare the different $\langle c_{s}^{2}\rangle\leavevmode\nobreak\
$– $C$ analytic relations we found above. As one can see from the left panel
of Fig. 1, for small compactness, all curves approach each other. More
quantitatively, the Taylor expansions of the analytic expressions shown above
imply that $\langle c_{s}^{2}\rangle\leavevmode\nobreak\
=C/2+{\cal{O}}(C^{2})$, and the coefficient of the $C^{2}$ term is of order
unity. This observation leads us to speculate that the emergence of
universality in this limit may be attributed to the predominance of Newtonian
effects over relativistic effects, at least for EoSs close to those considered
here. To further explore this point, in the next subsection we study the
$\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$ relation in the Newtonian
limit for two families of one-parameter EoSs, which interpolate between the
models considered in this subsection. We will later perform a post-Minkowskian
expansion to elucidate the role of GR in the quasi-universality.

###  III.2 One-Parameter Families of EoSs:
Newtonian Limit

Given the common Newtonian limit of the averaged speed of sound squared,
regardless of the EoS considered above, it is of interest to study the
Newtonian limit of the TOV equations directly for other EoSs. By “Newtonian
limit,” we here mean the leading-order expansion of the equations of structure
when $c\to\infty$, or equivalently, when $G\to 0$. Below, we present the
$\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$ relation for a polytropic
EoS and a generalized version of the Tolman VII solution. Notably, both EoSs
possess analytic solutions, at least for specific values of the parameters
that characterize them.

Let us then begin by presenting the Newtonian equations of structure. Using
that in the Newtonian limit $\epsilon=\rho c^{2}$, where $\rho$ is the rest
mass density, Eqs. (3) and (4) reduce to

$\displaystyle\frac{dp}{dr}=-\rho\frac{m}{r^{2}},\qquad\frac{dm}{dr}=\rho\frac{Gm}{r^{2}}.$
(30)

Just as in the relativistic case, these equations must be closed by choosing
an EoS and prescribing boundary conditions.

####  III.2.1 Polytropic EoSs

A polytropic EoS is defined as

$\displaystyle p=K\,\rho^{1+\frac{1}{n}},$ (31)

where $K$ is the polytropic constant, and $n$ is the polytropic index. For
Newtonian polytropes, it is customary to define the dimensionless variables

$\displaystyle\xi\equiv r/r_{0},\qquad\mu^{(N,P)}\equiv m/m_{0}^{(N,P)},$ (32)

where

$\displaystyle r_{0}^{2}\equiv\frac{(n+1)p_{c}}{4\pi G\rho_{c}^{2}},\qquad
m_{0}^{(N,P)}=4\pi\rho_{c}r_{0}^{3},$ (33)

are given in terms of the central rest mass density $\rho_{c}=\rho(r=0)$ and
the superscript $(N,P)$ stands for “Newtonian polytrope.” The stellar radius
is then given by $R=r_{0}\xi_{1}$ and is determined by finding the value of
the dimensionless variable $\xi$ at which $p(\xi=\xi_{1})=0$. Similarly, the
stellar mass is given by $M=m_{0}\mu_{1}^{(N,P)}$, where
$\mu_{1}^{(N,P)}=\mu^{(N,P)}(\xi=\xi_{1})$.

With this in hand, we can now compute the averaged speed of sound squared.
From Eq. (33) it follows that

$\displaystyle\frac{Gm_{0}}{r_{0}c^{2}}=(n+1)\frac{p_{c}}{\rho_{c}c^{2}}=(n+1)\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ ,$ (34)

so that the stellar compactness is linearly related to $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $via

$\displaystyle\langle c_{s}^{2}\rangle\leavevmode\nobreak\
=C\left[(n+1)\frac{\mu^{(N,P)}_{1}}{\xi_{1}}\right]^{-1}.$ (35)

In this expression, it is evident that both $(n+1)$ and
$\mu_{1}^{(N,P)}/\xi_{1}$ are highly dependent of the EoS parameter $n$.
Consequently, any universality must arise from the combination of both terms.

The polytropic index values that suitably approximate the behavior of
realistic EoSs for NSs are limited to $n\in[0,1]$. For the boundaries, $n=0$
(the incompressible limit) and $n=1$ (the Newtonian limit of a Buchdahl
fluid), the equations of hydrostatic equilibrium can be analytically solved to
find

$\displaystyle n=0\longrightarrow$ $\displaystyle\quad
p(\xi)=p_{c}\left(1-\frac{\xi^{2}}{6}\right),\quad\mu^{(N,P)}(\xi)=\frac{\xi^{3}}{3},$
$\displaystyle\quad\xi_{1}=\sqrt{6},\quad\mu_{1}^{(N,P)}=2\sqrt{6}.$ (36)
$\displaystyle n=1\longrightarrow$ $\displaystyle\quad
p(\xi)=p_{c}\left(\frac{\sin\xi}{\xi}\right)^{2},\quad\mu^{(N,P)}(\xi)=\sin\xi-\xi\cos\xi,$
$\displaystyle\quad\xi_{1}=\pi,\quad\mu_{1}^{(N,P)}=\pi.$ (37)

Thus, one sees that in both cases

$\displaystyle\left<c_{s}^{2}\right>=\frac{C}{2}\qquad(n=0\;{\textrm{or}}\;1).$
(38)

In fact, by solving the Lane-Emden equation (i.e. the Newtonian version of the
TOV equation in the reduced variables defined above) numerically for values of
$n\in[0,1]$, we find that the maximum fractional difference in $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \\!/C$ with respect to a fiducial EoS
with $n=0.5$ (the EoS median) is $\sim 3\%$. This can be seen in Fig. 4, where
we show the ratio $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ \\!/C$ for the
extended interval $n\in[0,2]$. This ratio is notably flat for $n\in[0,1]$,
with a stationary point at $n\approx 0.454$, but rapidly increases for higher
values of $n$, eventually diverging as $n\to 5$.

Figure 4: Ratio between the average speed of sound squared and the compactness
for 2 one-parameter families of solutions to the Newtonian equations of
structure, the polytropic EoS (in pink) and the generalized Tolman VII
solution (in dashed cyan). We use the polytropic index $n$ to parameterize the
former, and for the latter, we use the quantity $2/k$, which reduces to the
incompressible case when $2/k=0$ and to the Tolman VII case when $2/k=1$.

####  III.2.2 Generalized Tolman VII EoS

Let us now consider fixing the density profile through a generalized Tolman
VII EoS, which we define as

$\displaystyle\rho(r)=\rho_{c}\left[1-\left(\frac{r}{R}\right)^{k}\right]\,.$
(39)

The constant $k$ must here satisfy $k>1$ to ensure that $d\rho/dr|_{r=0}=0$.
When $k=2$, the above profile reduces to the Newtonian version of the Tolman
VII profile shown in Eq. (11), which is thought to be a reliable approximation
for the density profile inside NSs. As $k\to\infty$, we recover the flat
profile characteristic of an incompressible fluid (see Eq. (5)). Therefore,
Eq. (39) provides an interesting way of “interpolating” between the Tolman VII
and the incompressible fluid profiles, which allows for the analytic analysis
to be developed further.

Let us define the dimensionless variables

$\displaystyle x\equiv r/R,\qquad\mu^{(N,T)}\equiv\frac{m}{m_{0}^{(N,T)}},$
(40)

where $R$ is the stellar radius, and the superscript $(N,T)$ indicates the
Newtonian generalized Tolman VII solution, where we set

$\displaystyle m_{0}^{(N,T)}\equiv 4\pi\rho_{c}R^{3}.$ (41)

Then, for arbitrary values of $k$, the equation of hydrostatic equilibrium and
the equation for the enclosed mass can be solved analytically to find

$p(x)=p_{c}\left[1-\frac{\beta}{\langle c_{s}^{2}\rangle\leavevmode\nobreak\
}\left(\frac{x^{2}}{6}-\frac{x^{k+2}(k+6)}{3(k+2)(k+3)}+\frac{x^{2k+2}}{2(k+1)(k+3)}\right)\right],$
(42) $\mu^{(N,T)}(x)=\frac{x^{3}}{3}-\frac{x^{k+3}}{k+3},$ (43)

where anticipating the notation we will use in our post-Minkowskian
calculations, we have defined $\beta\equiv Gm_{0}/Rc^{2}$.

With this in hand, we can now compute the averaged speed of sound squared.
From the fact that $p(x=1)=0$, and considering that

$\displaystyle
C=\frac{GM}{Rc^{2}}=\beta\mu^{(N,T)}(1)=\beta\left(\frac{1}{3}-\frac{1}{k+3}\right),$
(44)

one finds the simple relation

$\displaystyle\left<c_{s}^{2}\right>=\left(2+\frac{1}{k}-\frac{3}{k+4}\right)^{-1}C,$
(45)

from where we see that

$\displaystyle\left<c_{s}^{2}\right>=\frac{C}{2}\quad\textrm{for }k=2\textrm{
and }k\to\infty.$ (46)

For $k\in[2,\infty)$ in Eq. (45), the maximum fractional difference in
$\langle c_{s}^{2}\rangle\leavevmode\nobreak\ \\!/C$ with respect to the EoS
median is $\sim 6\%$, which characterizes the dependence of the relation with
respect to $k$ in such interval. Here and below, when referring to the
generalized Tolman VII model, the median is defined using a uniform
probability distribution in the variable $q=2/k\in[0,1]$, i.e., it refers to a
fiducial model with $q=0.5$ ($k=4$). The ratio $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \\!/C$ is shown in Fig. 4 for the full
range $k\in(1,\infty)$. Similarly to the polytropic case, $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \\!/C$ is relatively flat over the range
$k\in[2,\infty)$, displaying a stationary point between the incompressible and
the Tolman VII solutions, when $k=2(1+\sqrt{3})\approx 5.46$, or equivalently,
$q\approx 0.37$.

###  III.3 One-Parameter Families of EoSs:
Post-Minkowskian Expansion

In the previous subsection, we observed that the emergence of quasi-
universality is rooted in the Newtonian limit. Furthermore, upon examining the
exact solutions in GR, we noted that the contribution of higher-order terms
becomes increasingly dependent on the EoS. To more precisely quantify the
impact of GR on this dependence, we perform a post-Minkowskian expansion for
the aforementioned one-parameter families of EoSs, which allows us to
systematically study and analyze such effects semi-analytically.

####  III.3.1 Polytropic EoSs

The relativistic generalization of the polytropic EoS family is

$\displaystyle p(\rho)=K\rho^{1+\frac{1}{n}},\quad\text{with
}\quad\epsilon=\rho c^{2}+np,$ (47)

and recall, from Eq. (31), that $K$ is a constant, $\rho$ is the rest mass
density and $n$ is the polytropic index. We will once more use the
dimensionless quantities defined in Eq. (32), $\xi$ and $\mu^{(P)}$ here,
where the superscript $(P)$ indicates that we are referring to the polytropic
EoS. Additionally, we use the definitions $\bar{p}\equiv p/p_{c}$,
$\bar{\rho}\equiv\rho/\rho_{c}$,
$\bar{\epsilon}\equiv\epsilon/(\rho_{c}c^{2})$, where the subscript $c$ stands
for the central value of a particular variable. The relativistic structure
equations can then be written as

$\displaystyle\frac{d\mu^{(P)}}{d\xi}=\xi^{2}\bar{\epsilon},\qquad\frac{d\bar{p}}{d\xi}=-(n+1)\frac{(\bar{\epsilon}+\sigma\bar{p})(\mu^{(P)}+\sigma\xi^{3}\bar{p})}{\xi^{2}(1-2(n+1)\sigma\mu^{(P)}/\xi)}.$
(48)

The dimensionless parameter $\sigma\equiv p_{c}/(\rho_{c}c^{2})$, is related
to the average speed of sound squared via

$\langle c_{s}^{2}\rangle\leavevmode\nobreak\
=\frac{p_{c}}{\rho_{c}c^{2}+np_{c}}=\frac{\sigma}{1+n\sigma}.$ (49)

In the Newtonian limit, we have that $\sigma\ll 1$, and therefore $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \sim\sigma\ll 1$ which then implies that
all three quantities ($\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $,
$\sigma$, and $C$) are natural expansion parameters in the post-Minkowskian
approximation.

By considering $\sigma$ as the small parameter, we can search for power series
solutions of the form

$\mu^{(P)}(\xi;\sigma)=\sum_{i=0}^{N}\mu_{(i)}^{(P)}(\xi)\sigma^{i}+\mathcal{O}(\sigma^{N+1}),$
(50)
$\bar{p}(\xi;\sigma)=\sum_{i=0}^{N}p_{(i)}(\xi)\sigma^{i}+\mathcal{O}(\sigma^{N+1}).$
(51)

Additionally, the dimensionless radius of the star and dimensionless mass of
the star can also be expanded in the form

$\xi_{1}(\sigma)=\sum_{i=0}^{N}\xi_{1,i}\sigma^{i},\qquad\mu_{1}^{(P)}(\sigma)=\sum_{i=0}^{N}\mu_{1,i}^{(P)}\sigma^{i},$
(52)

and can be obtained by solving the resulting equations perturbatively (see
Appendix A for details).

The averaged speed of sound squared as a function of the compactness can now
be found as follows. The compactness can be found by evaluating

$C=(n+1)\sigma\frac{\mu_{1}^{(P)}(\sigma)}{\xi_{1}(\sigma)}$ (53)

as a power series in $\sigma$. This series can be rewritten in terms of
$\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $, and can subsequently be
inverted into a power series in $C$. For example, for an $n=1$ polytrope, the
next-to-leading-order term in the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$ relation can be obtained
analytically, and the subsequent terms can be obtained numerically, namely

$\displaystyle\langle c_{s}^{2}\rangle\leavevmode\nobreak\ \approx$
$\displaystyle\frac{C}{2}+\frac{C^{2}}{16}\left[8+9\text{Ci}(\pi)-9\text{Ci}(3\pi)+9\ln
3\right]+2.97C^{3}+\mathcal{O}(C^{4})$ $\displaystyle\approx$
$\displaystyle\frac{C}{2}+1.15C^{2}+2.97C^{3}+\mathcal{O}(C^{4}).$ (54)

Observe that, similarly to the analytic solutions, the coefficient of the
$C^{2}$ term is of order unity.

The post-Minkowskian expansion of the averaged speed of sound with $n\in[0,1]$
is shown in Fig. 5, which demonstrates the effects of higher-order terms on
the $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$ relation. A
significant observation derived from this plot is that incorporating the
$C^{2}$ and $C^{3}$ terms results in an increased sensitivity to the EoS in
the $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $ – $C$ relation. For a
fiducial compactness of $C=0.2$, the maximal fractional deviation with respect
to the EoS median, for polytropes with $n\in[0,1]$, raises from $\sim 3\%$ at
$\mathcal{O}(C)$ (Newtonian limit) to $\sim 7\%$ at $\mathcal{O}(C^{2})$,
$\sim 9\%$ at $\mathcal{O}(C^{3})$, and $\sim 17\%$ in full GR. This
observation is consistent with the notion that the insensitivity to the EoS
(for $n\in[0,1]$), which stems from the Newtonian limit, deteriorates (i.e.
the $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$ relation becomes more
sensitive to the EoS) as GR contributions are included.

Figure 5: The variation of $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $in
relation to the post-Minkowskian expansion for the polytropic equation of
state, considering values of $n\in[0,1]$ along the bottom axis and
representing different compactness levels with distinct colors. The expansion
analysis starts from the first-order term and progressively incorporates
higher-order contributions, up to the third order, proceeding from left to
right. The black dotted line indicates the numerical solution in full GR for a
NS with $C=0.2$ for varying values of $n$.

####  III.3.2 Generalized Tolman VII

Let us now consider the case of the generalized Tolman VII density profile,
which in the relativistic limit becomes

$\epsilon(r)=\epsilon_{c}\left[1-\left(\frac{r}{R}\right)^{k}\right],$ (55)

where $k>1$. This family of solutions interpolates between the Tolman VII
solution, for $k=2$, and the incompressible case, that once again corresponds
to $k\to\infty$. Taking the Newtonian limit, $\epsilon=\rho c^{2}$, one
recovers the results presented in Sec. III.2.2. As done before, we define the
dimensionless quantities in Eq. (40), but replacing Eq. (41) by

$m_{0}^{(T)}\equiv 4\pi\epsilon_{c}R^{3}.$ (56)

We use the same definition of $\beta$, with $\bar{p}=p/p_{c}$ and
$\bar{\epsilon}=\epsilon/\epsilon_{c}$. The equation for the enclosed mass can
be directly integrated and results in

$\mu^{(T)}(x)=\frac{x^{3}}{3}-\frac{x^{k+3}}{k+3},$ (57)

where the superscript ${(T)}$ indicates the generalized Tolman VII solution
and where once again the compactness and the parameter $\beta$ are linearly
related by

$C=\frac{GM}{Rc^{2}}=\beta\mu^{(T)}(1)=\beta\left(\frac{1}{3}-\frac{1}{k+3}\right).$
(58)

The dimensionless pressure then obeys the differential equation

$\frac{d\bar{p}}{dx}=-\frac{\beta}{\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ }\frac{(\bar{\epsilon}+\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \bar{p})(\mu^{(T)}+\langle
c_{s}^{2}\rangle\leavevmode\nobreak\
x^{3}\bar{p})}{x^{2}(1-2\beta\mu^{(T)}/x)}.$ (59)

Let us now choose the stellar compactness as our small post-Minkowskian
expansion parameter to obtain power series approximations to the solution for
the pressure.

$\bar{p}(x;C)=\sum_{i=0}^{N}p_{i}(x)C^{i},\qquad\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ (C)=\sum_{i=0}^{N}d_{i+1}C^{i+1},$ (60)

where we write the relevant boundary condition as $\bar{p}(0;C)=1$, while the
requirement that $\bar{p}(1;C)=0$ determines the coefficients in the expansion
of $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ (C)$. For the generalized
Tolman VII models, the structure equations can be solved analytically, order
by order, for any $k$ and up to high values of $N$. Here we show such analytic
solution up to order $C^{3}$ for generic values of $k$

$\displaystyle\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $
$\displaystyle=\frac{1+2\left(\frac{2}{k}\right)}{\left(1+\left(\frac{2}{k}\right)\right)\left(2+\left(\frac{2}{k}\right)\right)}C$
$\displaystyle+\frac{1+\frac{17}{3}\left(\frac{2}{k}\right)+\frac{239}{24}\left(\frac{2}{k}\right)^{2}+\frac{133}{24}\left(\frac{2}{k}\right)^{3}}{3\left(1+\left(\frac{2}{k}\right)\right)^{2}\left(3+8\left(\frac{2}{k}\right)+4\left(\frac{2}{k}\right)^{2}\right)}C^{2}$
$\displaystyle+\frac{17}{2}\frac{1+\frac{33}{4}\left(\frac{2}{k}\right)+\frac{3165}{136}\left(\frac{2}{k}\right)^{2}+\frac{458}{17}\left(\frac{2}{k}\right)^{3}+\frac{1497}{136}\left(\frac{2}{k}\right)^{4}}{\left(1+\left(\frac{2}{k}\right)\right)^{3}\left(4+15\left(\frac{2}{k}\right)+9\left(\frac{2}{k}\right)^{2}\right)}C^{3}$
$\displaystyle+\mathcal{O}(C^{4}).$ (61)

The influence of higher-order post-Minkowskian terms (i.e. higher order in $C$
terms) on the EoS-dependence of the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \\!-C$ relation is shown in Fig. 6.
Introducing the $C^{2}$ and $C^{3}$ terms leads to an increased variability in
the $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ \\!-C$ relation. This
observation again, as in the polytropic case, aligns with the notion that the
insensitivity to the EoS (for generalized Tolman VII models), which originates
in the Newtonian limit, diminishes as GR contributions are added. For a
fiducial compactness of $C=0.2$, the maximum fractional deviation of $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $with respect to the EoS median, for
$q=2/k\in[0,1]$, increases from $\sim 6\%$ at $\mathcal{O}(C)$ (Newtonian
limit) to $\sim 8\%$ at $\mathcal{O}(C^{2})$, $\sim 9\%$ at
$\mathcal{O}(C^{3})$, and $\sim 13\%$ in full GR.

Figure 6: The variation of $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $in
relation to the post-Minkowskian expansion for the generalized Tolman VII
solution, considering values of $k\in[2,\infty)$ along the bottom axis and
representing different compactness levels with distinct colors. The expansion
analysis starts from the first-order term and progressively incorporates
higher-order contributions, up to the third order, proceeding from left to
right. The black dotted line indicates the numerical solution in full GR for a
NS with $C=0.2$ for varying values of $k$.

###  III.4 Comparison to Realistic EoSs

To conclude this section, let us assess how useful the results obtained above
(for simple, one-parameter families of EoSs) can be for understanding the
$\langle c_{s}^{2}\rangle\leavevmode\nobreak\ \\!-C$ relation when using
realistic EoSs. To address this question, we consider a set of 15 tabulated
EoSs, which is the same as that employed in Fig. 1. In particular, we consider
the KDE0v [21], KDE0v1 [21], Rs [22], SK255 [21], SK272 [21], SKI2-6 [23],
SkMp [24], SkOp [23], SLY2 [25], SLY230a [25], and SLY9 [25] EoSs [26, 27],
with tables obtained from the CompOSE repository.

Realistic EoSs encode the rich behavior of nuclear matter from the outer crust
to the inner core. The adiabatic index, $\Gamma=d(\log p)/d(\log\rho)$
predicted by these EoSs typically spans a large interval, starting from $\sim
1.33$ at the outer crust, decreasing after neutron drip, and then increasing
considerably (typically above 3) after the crust-core interface [32].
Naturally, a single polytrope is not capable of capturing the full behavior of
the nuclear EoS, nor is the density profile fully described by a simple
analytic function as in Eq. (55). Still, it will be useful to construct
effective polytropic and generalized Tolman VII models to NSs described by our
set of realistic EoSs, in order to verify to which extent (or in which domain)
our results for the former models can be used to build intuition for the
latter.

From the many possible definitions of an effective polytropic index of an EoS
inside a NS, we select the simplest one, consisting of

$n_{\text{eff}}=\frac{1}{\Gamma(\rho_{c})–1}\,,$ (62)

and thus, $n_{\text{eff}}$ is fully determined by the adiabatic index at the
center of the star. Additionally, we define an effective generalized-Tolman-
VII exponent, $k_{\text{eff}}$, by fitting an expression of the form of Eq.
(55) to the energy-density profile of a given NS. To assess the adequacy of
these effective descriptions, we introduce the following error measure,

$E=\sqrt{\left(\frac{\Delta\langle c_{s}^{2}\rangle\leavevmode\nobreak\
}{\langle c_{s}^{2}\rangle\leavevmode\nobreak\ }\right)^{2}+\left(\frac{\Delta
C}{C}\right)^{2},}$ (63)

where $\Delta$ denotes the difference between the quantities predicted by the
effective model (the polytropic or generalized Tolman VII mimickers) and the
original tabulated EoSs, which is then divided by the quantity predicted by
the original EoS.

Figure 7: Effective indices and error norm as a function of the stellar
compactness. The top panel shows the effective polytropic index
$q_{\text{eff}}=n_{\text{eff}}$ (solid lines) and the effective generalized
Tolman VII exponent $q_{\text{eff}}=k_{\text{eff}}$ (dashed lines) for fits to
NSs described by a set of 15 tabulated EoSs. The bottom panel shows the error
measure, as defined in Eq. (63).

Figure 7 shows the effective indices
$q_{\text{eff}}=\\{n_{\text{eff}},2/k_{\text{eff}}\\}$ as a function of the
stellar compactness, for our set of tabulated EoSs, while the bottom panel
shows the error $E$ as defined in Eq. (63). In the range of compactness
considered here, $n_{\text{eff}}\in(0,1)$ for a polytropic description of the
EoS, but it fails to reasonably describe the NS properties for $C\lesssim
0.1$, as the error $E$ rapidly increases. Physically, this is due to the fact
that the effective polytrope underestimates the stellar radius for low
compactness NSs, for which the crust occupies a larger portion of the stellar
volume. This can be visualized in Fig. 8, where the energy density profile is
represented for very low ($C=0.05$), medium ($C=0.2$), and high-compactness
($C=0.28$) NSs.

Analogously, an effective generalized Tolman VII model also fails to describe
low compactness stars, with $k_{\text{eff}}$ eventually dropping below 1,
which is physically unacceptable. By construction, these models preserve the
radius of the star but may fail to reasonably predict its mass (cf. Fig. 8),
giving rise to large errors for small values of $C$. For $C\gtrsim 0.1$,
$E\lesssim 0.1$ and our results for these simple one-parameter families of
EoSs can be extended to interpret the case of realistic EoSs.

Figure 8: Energy density profiles for a low (left panel), medium (middle
panel), and a high-compactness (right panel) NS. In each panel, the radial
profile is shown for the tabulated (SLY9) EoS, the effective polytrope, and
the effective generalized Tolman VII model, as defined in the main text.

##  IV I – Love – $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $

Alongside the near EoS-insensitive relation between $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $and $C$, [20] presented a similar
relation between the average speed of sound squared and both the dimensionless
moment of inertia $\bar{I}$ and the dimensionless tidal deformability
$\Lambda$, which is insensitive to the EoS at the $\sim 10\%$ level for
realistic EoSs. Having explored the origins of the quasi-universality of the
$\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$ relation, and arguing
that, for relevant models that describe the NS interior, this quasi-
universality emerges from the Newtonian limit and subsequently decays with
increasing contributions of GR, we will now investigate the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $ – $\bar{I}/\Lambda$ relations. Due to
the lack of analytic solutions in GR for both these quantities, we start by
analyzing the relations in the Newtonian regime and proceed to study their
corresponding post-Minkowskian corrections.

Before we begin, let us define the dimensionless quantities that shall be used
throughout this section. The dimensionless moment of inertia and the
dimensionless tidal deformability are respectively defined as

$\displaystyle\bar{I}\equiv\frac{c^{4}I}{G^{2}M^{3}},\qquad\Lambda=\frac{2k_{2}}{3C^{5}},$
(64)

where $I$ is the moment of inertia and $k_{2}$ is the $l=2$ tidal Love number.

###  IV.1 Newtonian Limit

To delve into the question of the role of the Newtonian limit in the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $ – $\bar{I}/\Lambda$ relations, we study
them in this limit and compare any universality, to that present in the
$\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$ relation, or lack
thereof. The moment of inertia in the Newtonian regime, which we shall name
$I^{(N)}$, is defined as

$I^{(N)}=\frac{8\pi}{3}\int_{0}^{R}r^{4}\rho(r)dr.$ (65)

The tidal deformability is found through the perturbation of a non-rotating
equilibrium state. If a static quadrupolar tidal field perturbs such state,
the Eulerian change $h_{2}(r)$ to the gravitational potential can be found
through the Newtonian limit of Eq. (15) in [38]:

$\frac{d^{2}h_{2}}{dr^{2}}+\frac{2}{r}\frac{dh_{2}}{dr}-\left[\frac{6}{r^{2}}-4\pi
G\rho\frac{d\rho}{dp}\right]h_{2}=0,$ (66)

and noticing that regularity at $r=0$ implies that $h_{2}(r)\sim r^{2}$ as
$r\to 0$. Defining $y\equiv Rh_{2}^{\prime}(R)/h_{2}(R)$, the Newtonian $l=2$
Love number can be computed from

$k_{2}^{(N)}=\frac{1}{2}\left(\frac{2-y}{y+3}\right),$ (67)

which is then related to the dimensionless tidal deformability via Eq. (64).

####  IV.1.1 Polytropic EoSs

Here, we return to polytropic EoSs, defined in the Newtonian limit in Eq.
(31). For the analytic solutions presented for polytropes in Sec. III.2, when
$n=0$ and $n=1$, we obtain analytic solutions for the moment of inertia

$\displaystyle n=0\longrightarrow$
$\displaystyle\quad\bar{I}^{(N,P)}=\frac{2}{5}\frac{1}{C^{2}},$ (68)
$\displaystyle n=1\longrightarrow$
$\displaystyle\quad\bar{I}^{(N,P)}=\left(\frac{2}{3}-\frac{4}{\pi^{2}}\right)\frac{1}{C^{2}}\,,$
(69)

which can then be related to the averaged speed of sound squared via

$\displaystyle n=0\longrightarrow$ $\displaystyle\quad\langle
c_{s}^{2}\rangle\leavevmode\nobreak\
=\frac{1}{\sqrt{10}}\left(\frac{1}{\bar{I}^{(N,P)}}\right)^{1/2},$ (70)
$\displaystyle n=1\longrightarrow$ $\displaystyle\quad\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \approx
0.26\left(\frac{1}{\bar{I}^{(N,P)}}\right)^{1/2}.$ (71)

Unlike the case of the $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$
relation, there is no common limit of the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $– $\bar{I}$ relation between the
incompressible model and the $n=1$ polytrope. For polytropes with $n\in[0,1]$,
the quantity $\left(\bar{I}^{(N,P)}\right)^{1/2}\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $ displays a fractional difference with
respect to the EoS median of (at most) $\sim 14\%$. This quantity is shown in
Fig. 9 as a function of $n\in[0,1]$.

Such EoS dependence is also present in the relation between the dimensionless
tidal deformability and the average speed of sound squared. When $n=0$ and
$n=1$ we find that

$\displaystyle n=0\longrightarrow$
$\displaystyle\quad\Lambda^{(N,P)}=\frac{1}{2C^{5}}\,,$ (72) $\displaystyle
n=1\longrightarrow$
$\displaystyle\quad\Lambda^{(N,P)}=\frac{15-\pi^{2}}{3\pi^{2}}\frac{1}{C^{5}}\,,$
(73)

and thus

$\displaystyle n=0\longrightarrow$ $\displaystyle\quad\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \approx
0.44\left(\Lambda^{(N,P)}\right)^{-1/5}$ (74) $\displaystyle
n=1\longrightarrow$ $\displaystyle\quad\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \approx
0.35\left(\Lambda^{(N,P)}\right)^{-1/5}.$ (75)

For polytropes with $n\in[0,1]$, the quantity $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \left(\Lambda^{(N,T)}\right)^{1/5}$
displays a fractional difference with respect to the EoS median of (at most)
$\sim 14\%$. This quantity is shown in Fig. 9 as a function of $n\in[0,1]$.

####  IV.1.2 Generalized Tolman VII

In accordance with the methodology employed in the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$ case, we now examine the
$\bar{I}/\Lambda$ – $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $ relations
using the generalized Tolman VII profile. The relevant EoS expression and
equilibrium equations pertaining to this solution were already discussed in
Sec. III.2.2.

By substituting the density profile, $\rho(r)$, presented in Eq. (39), into
the integral shown in Eq. (65), we find that

$\bar{I}^{(N,T)}=\frac{8\pi}{3}\\!\\!\int_{0}^{R}\\!\\!\\!r^{4}\rho_{c}\left[1-\left(\frac{r}{R}\right)^{k}\right]dr=\frac{8\pi}{3}R^{5}\rho_{c}\left(\frac{1}{5}-\frac{1}{k+5}\right),$
(76)

where the superscript $(N,T)$ indicates we are in the Newtonian regime and
using the generalized Tolman VII solution, and recall that $\rho_{c}$ stands
for the central rest-mass density. By referring to Eq. (44), we can establish
the relation $4\pi\rho_{c}R^{2}={3C(k+3)}/{k}$. Consequently, substituting
this relation into the aforementioned result yields

$\displaystyle\bar{I}^{(N,T)}$
$\displaystyle=\frac{2}{C^{2}}\left(1+\frac{3}{k}\right)\left(\frac{1}{5}-\frac{1}{k+5}\right)$
$\displaystyle=\frac{k^{2}(k+3)(k+4)^{2}}{10(k+5)\left(k^{2}+3k+2\right)^{2}}\frac{1}{\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \\!\\!^{2}},$ (77)

where we have used the $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$
relation for the generalized Tolman VII solution in the Newtonian limit. With
this expression, we obtain the following limits when $k=2$ and
$k\rightarrow\infty$

$\displaystyle k=2\longrightarrow$
$\displaystyle\quad\bar{I}^{(N,T)}=\frac{1}{14\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \\!\\!^{2}}\longrightarrow\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \approx
0.27\left(\bar{I}^{(N,T)}\right)^{-\frac{1}{2}}$ (78) $\displaystyle
k\rightarrow\infty\longrightarrow$
$\displaystyle\quad\bar{I}^{(N,T)}=\frac{1}{10\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \\!\\!^{2}}\longrightarrow\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \approx
0.32\left(\bar{I}^{(N,T)}\right)^{-\frac{1}{2}}.$ (79)

As expected, when $k\rightarrow\infty$, which corresponds to the
incompressible case, one obtains the same limit as that found using the
polytropic EoS when $n=0$. Also, similarly to that case, there is no common
limit between the incompressible case and the Tolman VII profile. For
$q=2/k\in[0,1]$, the quantity $\langle c_{s}^{2}\rangle\leavevmode\nobreak\
\left(\bar{I}^{(N,T)}\right)^{1/2}$ displays a maximum fractional difference
with respect to the EoS median of $\sim 10\%$. This quantity is shown in Fig.
9 as a function of $q$.

Figure 9: $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ \bar{I}^{1/2}$ as a
function of $q=n$ for a polytropic EoS model (dot-dashed pink) and $2/k$ for a
generalized Tolman VII model (dashed cyan). Observe that for $q\in[0,1]$, the
quantity $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ \bar{I}^{1/2}$ varies
little, thus exhibiting an approximate EoS universality.

For the tidal deformability, we find that

$\displaystyle k=2\longrightarrow$ $\displaystyle\quad\Lambda^{(N,T)}\approx
0.22\frac{1}{C^{5}}\longrightarrow\langle c_{s}^{2}\rangle\leavevmode\nobreak\
\approx 0.37\left(\Lambda^{(N,T)}\right)^{-\frac{1}{5}}$ (80) $\displaystyle
k\rightarrow\infty\longrightarrow$
$\displaystyle\quad\Lambda^{(N,T)}=\frac{1}{2C^{5}}\longrightarrow\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \approx
0.44\left(\Lambda^{(N,T)}\right)^{-\frac{1}{5}}.$ (81)

For $q=2/k\in[0,1]$, the quantity $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \left(\Lambda^{(N,T}\right)^{1/5}$
displays a maximum fractional difference with respect to the EoS median of
$\sim 10\%$. This quantity is also shown in Fig. 10 as a function of $q$.

Figure 10: $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ \Lambda^{1/5}$ for
two types of solutions to the Newtonian equations of structure, the polytropic
EoS (in pink) and the generalized Tolman VII solution (in dashed cyan). We use
the polytropic index $n$ for the former, varying it in the interval that is
relevant for NSs, and for the latter, we use the quantity $2/k$ that varies
between the incompressible case, when $2/k=0$, and Tolman VII, when $2/k=1$.

####  IV.1.3 Comparing solutions

The $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ \\!-\bar{I}/\Lambda$
relations present a stronger dependence on the EoS than the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \\!-C$ relation, at least in the
Newtonian regime. The variation in the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \\!-\bar{I}/\Lambda$ relations is
depicted in Figs. 9 and 10, respectively, with an EoS-sensitivity of $\sim
14\%$ for a polytropic EoS with $n\in[0,1]$ and $\sim 10\%$ for the
generalized Tolman VII solution with $q=2/k\in[0,1]$. To put these numbers in
perspective, we can compare them to the EoS-sensitivity of the (Newtonian)
$\langle c_{s}^{2}\rangle\leavevmode\nobreak\ \\!-C$ relation, $3\%$ for
polytropes and $6\%$ for the generalized Tolman VII solution.

The parallels between the $\langle c_{s}^{2}\rangle\leavevmode\nobreak\
\\!-\bar{I}$ and $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ \\!-\Lambda$
relations are noteworthy, although unsurprising, as these observables are well
known for having a strong universal relation both in GR and in the Newtonian
limit [12]. Given the known quasi-universal relations $I$-Love-$C$ in GR [39,
40], one might wonder whether the influence of incorporating GR corrections,
as observed in the $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ \\!-C$
relation, will similarly manifest in the moment of inertia and tidal
deformability. We tackle this question next.

###  IV.2 Post-Minkowskian Expansion

To assess whether increasing GR contributions will lead to a reduction in the
universality for both $\langle c_{s}^{2}\rangle\leavevmode\nobreak\
\\!-\bar{I}/\Lambda$ relations, we perform a post-Minkowskian expansion. This
enables us to derive GR contributions in an order-by-order manner,
facilitating a more precise assessment of the impact of each contribution. In
this section, this will be done for the generalized Tolman VII models alone,
defined in Eq. (55).

The relativistic moment of inertia can be derived under the context of slow
rotation, as originally outlined by Hartle [41]:

$I=\frac{8\pi
R^{5}}{3}\int_{0}^{1}dx\,x^{4}\frac{(\epsilon+p)}{c^{2}}e^{(\lambda-\nu)/2}\bar{\omega}.$
(82)

Here, $g_{tt}=-e^{\nu}$ and $g_{rr}=e^{\lambda}$ are metric coefficients
describing the unperturbed (i.e., non-rotating) spacetime, while
$g_{t\phi}=-\Omega(1-\bar{\omega})r^{2}\sin^{2}\theta$ describes the
perturbation due to rotation (to first order in the angular velocity
$\Omega$). The function $\bar{\omega}(x)$, with $x\equiv r/R$, obeys

$\frac{1}{x^{4}}\frac{d}{dx}\left(\,x^{4}j\frac{d\bar{\omega}}{dx}\right)+\frac{4}{x}\frac{dj}{dx}\bar{\omega}=0,$
(83)

with $j\equiv e^{-(\nu+\lambda)/2}$, along with the boundary conditions
$d\bar{\omega}/dx|_{x=0}=0$ and $\lim_{x\to\infty}\bar{\omega}(x)=1$.

The function $\lambda(x)$, for the generalized Tolman VII models, can be
directly derived from Eq. (57) as $\lambda(x)=-\log(1-\beta\mu/x)$. For the
metric potentials $\nu$ and $\bar{\omega}$, we employ a power series
approximation analogous to Eq. (60):

$\nu(x;C)=\sum_{i=0}^{N}\nu_{i}(x)C^{i},\qquad\bar{\omega}(x;C)=\sum_{i=0}^{N}\bar{\omega}_{i}(x)C^{i}.$
(84)

We can now solve Eq. (83) and the corresponding equation for $\nu$,

$\frac{d\nu}{dx}=-\frac{2\langle c_{s}^{2}\rangle\leavevmode\nobreak\
}{\bar{\epsilon}+\langle c_{s}^{2}\rangle\leavevmode\nobreak\
\bar{p}}\frac{d\bar{p}}{dx}\,,$ (85)

analytically, order by order in $C$, which allows the relativistic moment of
inertia (82) to be computed in successive post-Minkowskian approximations. In
the Newtonian limit (i.e. when $N=0$), we recover Eq. (IV.1.2) for the
dimensionless quantity $\bar{I}$.

The relation $\bar{I}-\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $ can be
found using Eq. (III.3.2). As our interest lies in the inverse relation,
$\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $as a function of $\bar{I}$, we
perform a series expansion of the moment of inertia corresponding to an
expansion of $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $about 0. Since in
the Newtonian limit we have that $\langle c_{s}^{2}\rangle\leavevmode\nobreak\
\propto\bar{I}^{-1/2}$, we expand around $\bar{I}$ approaching infinity, or,
equivalently, $\bar{I}^{-1/2}\rightarrow 0$. The resulting series converges
for $\bar{I}^{-1/2}\lesssim 0.49$, or $\bar{I}\gtrsim 4.2$ (with the radius of
convergence being roughly independent of $k$). The effects of this expansion
in powers of $\bar{I}^{-1/2}$ on the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $– $\bar{I}$ relation is depicted in Fig.
11, where we show values of $\bar{I}$ for which $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ <1$ in GR. Additionally, we show the full
GR curve for a NS with $\bar{I}=10$, which agrees with inferences obtained in
[42]. Observe that, as we consider higher-order post-Minkowskian expansions,
$\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $ approaches the full GR value.
We expect that the correct solution is obtained when $N\rightarrow\infty$, for
values of $\bar{I}^{-1/2}$ within the radius of convergence of the series.

Importantly, we observe that the increasing post-Minkowskian orders seem to
have minimal effects on the EoS-dependence of the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $– $\bar{I}$ relation. For instance, for
the fiducial dimensionless moment of inertia $\bar{I}=10$, and considering
$q=2/k\in[0,1]$, the quantity $\langle c_{s}^{2}\rangle\leavevmode\nobreak\
$displays a maximum fractional difference with respect to the EoS median of
$10\%$ to $11\%$ to all orders. Thus, the degree of universality observed in
the Newtonian limit remains consistent as we move to higher post-Minkowskian
orders. This is to be compared to the effect of post-Minkowskian corrections
to the $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $– $C$ relation, which
tends to deteriorate the EoS quasi-universality, as discussed earlier.

Figure 11: The variation of $\langle c_{s}^{2}\rangle\leavevmode\nobreak\
\\!-\bar{I}$ in relation to the post-Minkowskian expansion for the generalized
Tolman VII solution, considering values of $k\in[2,\infty)$ along the bottom
axis and representing different moment of inertia levels with distinct colors.
The expansion analysis starts from the first-order term and progressively
incorporates higher-order terms, up to the third order, proceeding from left
to right. The black dotted line indicates the numerical solution in full GR
for a NS with $\bar{I}=10$ for varying values of $k$. The minimum value for
the moment of inertia is determined by the point where $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $ becomes greater than 1 for the GR
solution.

Now let us consider the tidal deformability. This quantity is obtained from
the generalization of Eq. (66) into

$\displaystyle\frac{d^{2}h_{2}}{dr^{2}}+\frac{dh_{2}}{dr}\left\\{\frac{2}{r}+e^{\lambda}\left[\frac{2Gm}{r^{2}c^{2}}-\frac{4\pi
Gr}{c^{4}}(\epsilon-p)\right]\right\\}$
$\displaystyle+h_{2}\left[-\frac{6e^{\lambda}}{r^{2}}+\frac{4\pi
Ge^{\lambda}}{c^{4}}\left(5\epsilon+9p+\frac{\epsilon+p}{dp/d\epsilon}\right)-\left(\frac{d\nu}{dr}\right)^{2}\right]=0,$
(86)

where $h_{2}(r)=-e^{-\nu}\delta g_{tt}$ represents the perturbation to the
$(t,t)$ component of the metric, in terms of which all other metric potentials
can be obtained. The relativistic $l=2$ Love number can be computed from Eq.
(50) of Ref. [43], recalling that $y\equiv Rh_{2}^{\prime}(R)/h_{2}(R)$.
Furthermore, we again define the dimensionless tidal deformability through Eq.
(64).

Despite $h_{2}(r)$ having an analytic solution in terms of hypergeometric
functions for the generalized Tolman VII profiles in the Newtonian limit,
exact analytic solutions could not be found to higher post-Minkowskian orders.
Still, to access the impact of post-Minkowskian corrections, we introduce the
power series approximation

$h_{2}(x;C)=\sum_{i=0}^{N}h_{2,i}(x)C^{i}$ (87)

where $x=r/R$ and $h_{2,i}(x)$ is calculated numerically.

The relation between the dimensionless tidal deformability and $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $can be found by rewriting $C$ as a
function of $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $using Eq.
(III.3.2). As we are interested in $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $as a function of $\Lambda$, we perform a
series expansion of the tidal deformability in powers of $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $, as done previously for the moment of
inertia. Since in the Newtonian limit we have that $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \propto\Lambda^{-1/5}$, we expand around
$\Lambda$ approaching infinity, or, equivalently, $\Lambda^{-1/5}\rightarrow
0$. This expansion, up to the third order, is displayed in Fig. 12, where we
show values of $\Lambda$ for which $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \\!<1$ in GR. Additionally, we show the
numerical solution for a NS with $\Lambda=100$, which agrees with measurements
for a usual NS [2]. Observe again that as we include higher and higher post-
Minkowskian orders, $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $ approaches
the full GR result. We expect that the correct solution is obtained when
$N\rightarrow\infty$, for values of $\bar{\Lambda}^{-1/5}$ within the radius
of convergence of the series solution.

Importantly, we observe that, for the values for which the series is
convergent, the increasing post-Minkowskian orders seem again to have minimal
effects on the EoS-dependence of the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \\!-\Lambda$ relation. For the fiducial
tidal deformability $\Lambda=100$, and considering $q=2/k\in[0,1]$, the
quantity $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $displays a maximum
fractional difference with respect to the EoS median of $10\%$ to $11\%$ to
all orders. Thus, the degree of universality observed in the Newtonian limit
remains consistent as we move to higher post-Minkowskian orders.

Figure 12: The variation of $\langle c_{s}^{2}\rangle\leavevmode\nobreak\
\\!-\Lambda$ in relation to the post-Minkowskian expansion for the generalized
Tolman VII solution, considering values of $k\in[2,\infty)$ along the bottom
axis and representing different tidal deformability levels with distinct
colors. The expansion analysis starts from the first-order term and
progressively incorporates higher-order terms, up to the third order,
proceeding from left to right. The black dotted lines indicate the numerical
solution in full GR for a NS with $\Lambda=100$ for varying values of $k$.

##  V Conclusions

We have here explored the underlying physics responsible for the approximately
universal relations between the ratio of the central pressure to the central
energy density, $p_{c}/\epsilon_{c}$, and dimensionless astrophysical
observables, namely, the compactness $C$, the tidal deformability $\Lambda$
and the moment of inertia $\bar{I}$ of NSs [20]. The quantity
$p_{c}/\epsilon_{c}$ was reinterpreted as the squared speed of sound averaged
over the range of energy densities present inside the star, $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \\!$. The averaged speed of sound squared
has a natural interpretation as a measure of the mean stiffness of the EoS up
to the central energy density of a given NS.

We first studied the $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ \\!-C$
relation using analytic solutions of the TOV equation for an incompressible
fluid, a Tolman VII profile, and the Buchdahl fluid. We found that, albeit
quite different from each other, these EoSs led to $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \\!-C$ relations with the same Newtonian
limit, $\langle c_{s}^{2}\rangle\leavevmode\nobreak\ =C/2$. This prompted
further study of the Newtonian limit and its relativistic corrections, which
we performed for 2 one-parameter families of EoSs that interpolate between the
(Newtonian limit of the) models considered previously: polytropic EoSs and a
generalized Tolman VII profile. For the range of parameters most relevant for
the description of NSs, the $\langle c_{s}^{2}\rangle\leavevmode\nobreak\
\\!-C$ relation was shown to be remarkably flat (with respect to EoS
variation) in the Newtonian limit. This insensitivity to the EoS then
deteriorated as relativistic corrections were included. For instance, for
polytropic EoSs with $n\in[0,1]$, the maximal fractional difference of the
$\langle c_{s}^{2}\rangle\leavevmode\nobreak\ \\!-C$ relation with respect to
the EoS median increased from $3\%$ at $\mathcal{O}(C)$ (Newtonian limit), to
$9\%$ at $\mathcal{O}(C^{3})$ and 17% in full GR.

When contrasting realistic EoSs to the one-parameter families considered
previously, we found that the latter provides reasonable approximants to the
former, as long as the compactness of the NS is sufficiently large ($C\gtrsim
0.1$). Thus, the EoS insensitivity of the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \\!-C$ relation, for realistic EoSs, can
be linked to EoSs being relatively stiff throughout a large portion of the NS
for sufficiently large values of the stellar compactness.

In contrast to the case of the $\langle c_{s}^{2}\rangle\leavevmode\nobreak\
\\!-C$ relations, we found that the approximate universality of the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \\!-\bar{I}/\Lambda$ relations is not
enhanced in the Newtonian limit. Indeed, for the generalized Tolman VII
profile with $k\in[2,\infty)$, we showed that the EoS dependence remains at
the $10-11\%$ level to all orders in compactness. This is also the case when
using realistic EoSs, as long as the NS is sufficiently compact.

Our work suggests a few avenues for future work. One possibility would be to
study potential approximately universal relations between $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $and the quadrupole and higher-order
multipole moments of rotating NSs. This could be done by considering the slow-
rotation expansion [41], and constructing equilibrium sequences of stars with
fixed rotational periods. Such a study would also benefit from the post-
Minkowskian analysis we laid out in this work.

###### Acknowledgements.

N. Y. acknowledges support from NSF Award PHY-220765. R. M. acknowledges
support from the National Council for Scientific and Technological Development
(CNPq), and by the Carlos Chagas Filho Research Support Foundation (FAPERJ).
We thank Jaki Noronha-Hostler for helpful discussions about this paper.

##  Appendix A Post-Minkowskian expansion of $C-\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $ for polytropic EoSs

Differently from the generalized Tolman VII density profile, for which the
post-Minkowskian expansion of the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \\!-C$ relation can be found
analytically, in the case of a generic polytropic EoS, the $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ \\!-C$ relation must be obtained
numerically. In this Appendix, we sketch the numerical procedure and discuss
certain subtleties that render the straightforward computation of high post-
Minkowskian coefficients of that expansion numerically challenging.

The compactness $C$ is related to $\langle
c_{s}^{2}\rangle\leavevmode\nobreak\ $by

$C=(n+1)\sigma\frac{\mu_{1}^{(P)}(\sigma)}{\xi_{1}(\sigma)},$ (88)

where $\sigma=(1/\langle c_{s}^{2}\rangle\leavevmode\nobreak\ -n)^{-1}$ [cf.
Eq. (49)], which can be expanded to any desired order. Here,
$\mu_{1}^{(P)}(\sigma)=\mu^{(P)}(\xi_{1}(\sigma))$, and the dimensionless
stellar radius, $\xi=\xi_{1}$, is determined by

$\bar{p}(\xi_{1}(\sigma);\sigma)=0.$ (89)

Inserting the Taylor expansion of Eq. (52) for $\xi_{1}(\sigma)$ in the above
equation, we can solve the latter order by order in a post-Minkowskian
expansion. For example, at Newtonian order, one obtains

$p_{(0)}(\xi_{1,0})=0,$ (90)

which implicitly defines $\xi_{1,0}$. The coefficient $\xi_{1,i}$ with $i\geq
1$ then depends on the respective pressure coefficient $p_{(i)}$, as well as
on derivatives of the lower-order pressure coefficients ($p_{(j)}$, $j<i$),
evaluated at $\xi_{1,0}$. Explicitly, for $i=1,2$, one has

$\displaystyle\xi_{1,1}$
$\displaystyle=-\frac{p_{(1)}(\xi_{1,0})}{p^{\prime}_{(0)}(\xi_{1,0})},$ (91)
$\displaystyle\xi_{1,2}$
$\displaystyle=-\frac{p_{(2)}(\xi_{1,0})+\xi_{1,1}p^{\prime}_{(1)}(\xi_{1,0})+\frac{1}{2}\xi_{1,1}^{2}p^{\prime\prime}_{(0)}(\xi_{1,0})}{p^{\prime}_{(0)}(\xi_{1,0})}.$
(92)

The numerical procedure then consists of the following steps:

  * (i)

numerically solving the Newtonian equations for $p_{(0)}$ and
$\mu_{(0)}^{(P)}$, and extracting the value of $\xi_{1,0}$ such that Eq. (90)
holds, as well as $\mu_{1,0}^{(P)}$;

  * (ii)

numerically solving the equations for $p_{(i)}$, $\mu_{(i)}^{(P)}$, for $i\geq
1$ up to the desired order in the fixed domain $\xi\in[0,\xi_{1,0}]$, and
computing the coefficients $\xi_{1,i}$ and $\mu_{1,i}^{(P)}$;

  * (iii)

computing the $C-\langle c_{s}^{2}\rangle\leavevmode\nobreak\ $ relation from
the Taylor-expanded version of Eq. (88).

Now, it can be shown that, around $\xi=\xi_{1,0}$, the $i$-th pressure
coefficient has the leading behavior $p_{(i)}\propto(\xi-\xi_{1,0})^{n-i+1}$.
In the range of interest $0\leq n\leq 1$, Eq. (91) presents a $0/0$
indeterminacy, and must be computed by taking the limit of that expression as
$\xi\to\xi_{1,0}$. This limit is expected to be finite and is easily handled
numerically. However, in the same range $0\leq n\leq 1$, Eq. (92) displays an
apparent divergence in the limit $\xi\to\xi_{1,0}$, since the numerator goes
as $(\xi-\xi_{1,0})^{n-1}$ and the denominator goes as $(\xi-\xi_{1,0})^{n}$.
A careful analysis of the differential equation for $p_{(2)}$ shows, however,
that the coefficient of the dominant term in the numerator of Eq. (92), i.e.,
$\propto(\xi-\xi_{1,0})^{n-1}$, vanishes exactly, so that the true leading-
order behavior of both the numerator and denominator is $(\xi-\xi_{1,0})^{n}$,
yielding a finite value for the ratio in the limit $\xi\to\xi_{1,0}$. A
similar situation is expected for higher-order coefficients ($\xi_{1,i}$,
$i>2$). However, numerically the cancellation of diverging terms described
above is not guaranteed due to the inevitable presence of numerical errors,
and a higher-precision calculation may be required.

## References

  * Abbott  _et al._ [2017] B. P. Abbott _et al._ (LIGO Scientific, Virgo), GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2018] B. P. Abbott _et al._ (LIGO Scientific, Virgo), GW170817: Measurements of neutron star radii and equation of state, Phys. Rev. Lett. 121, 161101 (2018), arXiv:1805.11581 [gr-qc] .
  * Kramer  _et al._ [2006] M. Kramer _et al._ , Tests of general relativity from timing the double pulsar, Science 314, 97 (2006), arXiv:astro-ph/0609417 .
  * Weisberg  _et al._ [2010] J. M. Weisberg, D. J. Nice, and J. H. Taylor, Timing Measurements of the Relativistic Binary Pulsar PSR B1913+16, Astrophys. J. 722, 1030 (2010), arXiv:1011.0718 [astro-ph.GA] .
  * Fonseca  _et al._ [2014] E. Fonseca, I. H. Stairs, and S. E. Thorsett, A Comprehensive Study of Relativistic Gravity using PSR B1534+12, Astrophys. J. 787, 82 (2014), arXiv:1402.4836 [astro-ph.HE] .
  * Miller  _et al._ [2019] M. C. Miller _et al._ , 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] .
  * Riley  _et al._ [2019] T. E. Riley _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._ [2021] M. C. Miller _et al._ , The Radius of PSR J0740+6620 from NICER and XMM-Newton Data, Astrophys. J. Lett. 918, L28 (2021), arXiv:2105.06979 [astro-ph.HE] .
  * Riley  _et al._ [2021] T. E. Riley _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] .
  * Maggiore  _et al._ [2020] M. Maggiore _et al._ , Science Case for the Einstein Telescope, JCAP 03, 050, arXiv:1912.02622 [astro-ph.CO] .
  * Evans  _et al._ [2021] M. Evans _et al._ , A Horizon Study for Cosmic Explorer: Science, Observatories, and Community (2021), arXiv:2109.09882 [astro-ph.IM] .
  * Yagi and Yunes [2013] K. Yagi and N. Yunes, I-Love-Q Relations in Neutron Stars and their Applications to Astrophysics, Gravitational Waves and Fundamental Physics, Phys. Rev. D 88, 023009 (2013), arXiv:1303.1528 [gr-qc] .
  * Yagi and Yunes [2016] K. Yagi and N. Yunes, Binary Love Relations, Class. Quant. Grav. 33, 13LT01 (2016), arXiv:1512.02639 [gr-qc] .
  * Yagi and Yunes [2017a] K. Yagi and N. Yunes, Approximate Universal Relations among Tidal Parameters for Neutron Star Binaries, Class. Quant. Grav. 34, 015006 (2017a), arXiv:1608.06187 [gr-qc] .
  * Yagi and Yunes [2017b] K. Yagi and N. Yunes, Approximate Universal Relations for Neutron Stars and Quark Stars, Phys. Rept. 681, 1 (2017b), arXiv:1608.02582 [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 .
  * Lim and Holt [2018] Y. Lim and J. W. Holt, Neutron star tidal deformabilities constrained by nuclear theory and experiment, Phys. Rev. Lett. 121, 062701 (2018), arXiv:1803.02803 [nucl-th] .
  * Lattimer and Steiner [2014] J. M. Lattimer and A. W. Steiner, Constraints on the symmetry energy using the mass-radius relation of neutron stars, Eur. Phys. J. A 50, 40 (2014), arXiv:1403.1186 [nucl-th] .
  * Drischler  _et al._ [2021] C. Drischler, J. W. Holt, and C. Wellenhofer, Chiral Effective Field Theory and the High-Density Nuclear Equation of State, Ann. Rev. Nucl. Part. Sci. 71, 403 (2021), arXiv:2101.01709 [nucl-th] .
  * 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] .
  * Agrawal  _et al._ [2003] B. K. Agrawal, S. Shlomo, and V. K. Au, Nuclear matter incompressibility coefficient in relativistic and nonrelativistic microscopic models, Phys. Rev. C 68, 031304 (2003), arXiv:nucl-th/0308042 .
  * Friedrich and Reinhard [1986] J. Friedrich and P. G. Reinhard, Skyrme-force parametrization: Least-squares fit to nuclear ground-state properties, Phys. Rev. C 33, 335 (1986).
  * Reinhard  _et al._ [1999] P. G. Reinhard, D. J. Dean, W. Nazarewicz, J. Dobaczewski, J. A. Maruhn, and M. R. Strayer, Shape coexistence and the effective nucleon-nucleon interaction, Phys. Rev. C 60, 014316 (1999), arXiv:nucl-th/9903037 .
  * Bennour  _et al._ [1989] L. Bennour, P.-H. Heenen, P. Bonche, J. Dobaczewski, and H. Flocard, Charge distributions of Pb-208, Pb-206, and Tl-205 and the mean-field approximation, Phys. Rev. C 40, 2834 (1989).
  * Chabanat  _et al._ [1997] E. Chabanat, J. Meyer, P. Bonche, R. Schaeffer, and P. Haensel, A Skyrme parametrization from subnuclear to neutron star densities, Nucl. Phys. A 627, 710 (1997).
  * Gulminelli and Raduta [2015] F. Gulminelli and A. R. Raduta, Unified treatment of subsaturation stellar matter at zero and finite temperature, Phys. Rev. C 92, 055803 (2015), arXiv:1504.04493 [nucl-th] .
  * Danielewicz and Lee [2009] P. Danielewicz and J. Lee, Symmetry Energy I: Semi-Infinite Matter, Nucl. Phys. A 818, 36 (2009), arXiv:0807.3743 [nucl-th] .
  * Yagi  _et al._ [2014] K. Yagi, L. C. Stein, G. Pappas, N. Yunes, and T. A. Apostolatos, Why I-Love-Q: Explaining why universality emerges in compact objects, Phys. Rev. D 90, 063010 (2014), arXiv:1406.7587 [gr-qc] .
  * Sham  _et al._ [2015] Y. H. Sham, T. K. Chan, L. M. Lin, and P. T. Leung, Unveiling the universality of I-Love-Q relations, Astrophys. J. 798, 121 (2015), arXiv:1410.8271 [gr-qc] .
  * Chan  _et al._ [2016] T. K. Chan, A. P. O. Chan, and P. T. Leung, Universality and stationarity of the I-Love relation for self-bound stars, Phys. Rev. D 93, 024033 (2016), arXiv:1511.08566 [gr-qc] .
  * Agrawal  _et al._ [2005] B. K. Agrawal, S. Shlomo, and V. K. Au, Determination of the parameters of a Skyrme type effective interaction using the simulated annealing approach, Phys. Rev. C 72, 014310 (2005), arXiv:nucl-th/0505071 .
  * Douchin and Haensel [2001] F. Douchin and P. Haensel, A unified equation of state of dense matter and neutron star structure, Astron. Astrophys. 380, 151 (2001), arXiv:astro-ph/0111092 .
  * Tolman [1939] R. C. Tolman, Static solutions of Einstein’s field equations for spheres of fluid, Phys. Rev. 55, 364 (1939).
  * Jiang and Yagi [2020] N. Jiang and K. Yagi, Analytic I-Love-C relations for realistic neutron stars, Phys. Rev. D 101, 124006 (2020), arXiv:2003.10498 [gr-qc] .
  * Fujimoto  _et al._ [2022] Y. Fujimoto, K. Fukushima, L. D. McLerran, and M. Praszalowicz, Trace Anomaly as Signature of Conformality in Neutron Stars, Phys. Rev. Lett. 129, 252702 (2022), arXiv:2207.06753 [nucl-th] .
  * Oppenheimer and Volkoff [1939] J. R. Oppenheimer and G. M. Volkoff, On massive neutron cores, Phys. Rev. 55, 374 (1939).
  * Buchdahl [1967] H. A. Buchdahl, General-Relativistic Fluid Spheres. III. a Static Gaseous Model, Astrophys. J. 147, 310 (1967).
  * Hinderer [2008] T. Hinderer, Tidal Love numbers of neutron stars, Astrophys. J. 677, 1216 (2008), arXiv:0711.2420 [astro-ph] .
  * Hinderer  _et al._ [2010] T. Hinderer, B. D. Lackey, R. N. Lang, and J. S. Read, Tidal deformability of neutron stars with realistic equations of state and their gravitational wave signatures in binary inspiral, Phys. Rev. D 81, 123016 (2010), arXiv:0911.3535 [astro-ph.HE] .
  * Postnikov  _et al._ [2010] S. Postnikov, M. Prakash, and J. M. Lattimer, Tidal Love Numbers of Neutron and Self-Bound Quark Stars, Phys. Rev. D 82, 024016 (2010), arXiv:1004.5098 [astro-ph.SR] .
  * Hartle [1967] J. B. Hartle, Slowly rotating relativistic stars. 1. Equations of structure, Astrophys. J. 150, 1005 (1967).
  * Silva  _et al._ [2021] H. O. Silva, A. M. Holgado, A. Cárdenas-Avendaño, and N. Yunes, Astrophysical and theoretical physics implications from multimessenger neutron star observations, Phys. Rev. Lett. 126, 181101 (2021), arXiv:2004.01253 [gr-qc] .
  * Damour and Nagar [2009] T. Damour and A. Nagar, Relativistic tidal properties of neutron stars, Phys. Rev. D 80, 084035 (2009), arXiv:0906.0096 [gr-qc] .