# Adiabatic radial perturbations of relativistic stars: analytic solutions to
an old problem

Paulo Luz<EMAIL_ADDRESS>Centro de Astrofísica e Gravitação -
CENTRA, Departamento de Física, Instituto Superior Técnico - IST, Universidade
de Lisboa - UL, Av. Rovisco Pais 1, 1049-001 Lisboa, Portugal,  Departamento
de Matemática, ISCTE - Instituto Universitário de Lisboa, Portugal  Sante
Carloni<EMAIL_ADDRESS>Institute of Theoretical Physics, Faculty of
Mathematics and Physics, Charles University, Prague, V Holešovičkách 2, 180 00
Prague 8, Czech Republic,  DIME, Università di Genova, Via all’Opera Pia 15,
16145 Genova, Italy,  INFN Sezione di Genova, Via Dodecaneso 33, 16146 Genova,
Italy

###### Abstract

We present a new system of equations that fully characterizes adiabatic,
radial perturbations of perfect fluid stars within the theory of general
relativity. The properties of the system are discussed, and, provided that the
equilibrium spacetime verifies some general regularity conditions, analytical
solutions for the perturbation variables are found. As illustrative examples,
the results are applied to study perturbations of selected classical exact
spacetimes, and the first oscillation eigenfrequencies are computed.
Exploiting the new formalism, we derive an upper bound for the maximum
compactness of stable, perfect fluid stars, which is equation-of-state-
agnostic and significantly smaller than the Buchdahl bound.

##  I Introduction

Relativistic compact stellar objects are among the most complex and, at the
same time, most fascinating gravitational systems. Similarly to black holes,
these objects represent strong gravity systems. However, they are
fundamentally different in that the nature and behavior of matter play a
prominent role in their structure and evolution. This fact makes their
theoretical description particularly challenging. The complexities associated
with these objects explain, on the one hand, the predominant use of numerical
techniques and, on the other, the necessity to develop and apply perturbative
approaches to understand their properties and dynamics.

The first attempt at the description of perturbations of stars in general
relativity was made by Chandrasekhar in 1964 [1, 2]. In those works,
Chandrasekhar focused on understanding the behavior of adiabatic, radial
perturbations and developed an integro-differential equation, the so-called
Chandrasekhar radial pulsation equation, to describe this somewhat simpler
type of perturbation.

The Chandrasekhar equation has had a crucial influence on the subsequent
studies on the dynamical behavior of perturbations of self-gravitating,
massive compact objects, and various methods were developed to compute the
oscillation eigenfrequencies directly from the pulsation equation [3]. Various
reformulations of Chandrasekhar’s original equation have been proposed to ease
the numerical treatment of adiabatic, radial oscillations of self-gravitating
fluids with realistic equations of state (see e.g. [4, 5, 6, 7, 8]). However,
to our knowledge, no work has so far tackled one of the main limitations of
Chandrasekhar’s approach, i.e., the issue of gauge dependence. Indeed, the
predictions of the Chandrasekhar equation on the stability of a given matter
configuration and on the very behavior of its perturbations are intrinsically
associated with the specific coordinate system considered and, therefore,
conditioned by the choice of the gauge.

Recently, the authors of this paper derived a completely covariant and gauge-
invariant theory of perturbations for static, Locally Rotationally Symmetric
of class II (LRS II) spacetimes [9]. LRS spacetimes are characterized by local
rotational symmetry about a given spatial direction, and their nonvortical
subclass, known as LRS II in the case of perfect fluid sources, can be proven
to contain all isotropic non-rotating spacetimes suitable to describe stellar
compact objects and even some slowly rotating ones. The new perturbation
framework is based on the so-called 1+1+2 covariant formalism [10, 11]. This
formalism is based on the procedure of covariant spacetime decomposition, such
that at, each point, the spacetime is foliated by surfaces orthogonal to two
vector fields, a timelike, and a spacelike vector field. Using this covariant
decomposition, it is possible to characterize different aspects of the
geometry of LRS II spacetimes and the thermodynamics of the matter fields that
permeate them in a geometrically clear and physically meaningful way (see,
e.g., Refs. [12, 13, 14, 15, 16, 17, 17, 18, 19, 20, 21]).

The aim of this paper is to apply the covariant gauge-invariant perturbation
theory of Ref. [9] to study adiabatic, radial perturbations of static, compact
stellar objects composed of a perfect fluid. In the body of the text, we will,
however, write the perturbation equations immediately in a non-covariant way
from the classical point of view of an observer locally comoving with the
matter and use the circumferential radius to identify points within the star.
There are essentially two reasons for this choice. First, we aim to provide a
set of equations that can be applied without any specific knowledge of the
covariant formalisms. Second, the above setting is usually considered in the
standard metric-based description of this type of perturbation, providing a
clear, familiar interpretation of the quantities and the equations. We will
then present exact solutions for the perturbations of well-known models for
equilibrium spacetimes in a gauge-invariant way.

The paper is organized as follows. In Section II, we will present the
equations describing adiabatic, radial perturbations, and propose an algorithm
to obtain power series exact solutions of the system. We also give a lower
bound for the minimum eigenvalue of the system. In Section III, we apply these
methods to some classical solutions, namely, will consider the interior
Schwarzschild, Tolman IV, Kuchowicz 2-III and Heintzmann IIa (as cataloged in
[22]) spacetimes, and present the behavior of the first eigenfunctions that
characterize the perturbations in the frame of the comoving observer. In
Section IV, using results for the perturbation of the interior Schwarzschild
solution, we conjecture on a general upper bound for instability of static,
stellar compact objects composed of a perfect fluid, which is independent of
the information on the equation of state of the perturbed matter fluid. We
then draw conclusions in Section V. The paper also contains four appendices.
In Appendix A, we introduce the general definitions for the 1+1+2 covariant
quantities. In Appendix B, we give the general covariant, gauge-invariant
perturbation equations for adiabatic, radial perturbations found from the
1+1+2 formalism. In Appendix C, we show how the original Chandrasekhar radial
pulsation equation can be recovered from the new gauge-invariant equations,
demonstrating their equivalence in the considered coordinate system, and in
Appendix D we present the general intermediate matrix for the power series
solutions of the system.

Throughout the article, we will work in the geometrized unit system where
$8\pi G=c=1$, and consider the metric signature $(-+++)$.

##  II Adiabatic radial perturbations

The Einstein field equations (EFE) are a set of nonlinear partial differential
equations that are manifestly difficult to solve for general source matter
fields. To circumvent these difficulties, it is helpful to devise perturbation
schemes to linearize the field equations and study the behavior of small
deviations from equilibrium solutions. In Ref. [9], a general set of covariant
gauge-invariant equations was derived that can characterize linear
perturbations of static, spatially compact, spherically symmetric solutions of
the Einstein field equations. We will call such solutions “stars” because of
their most immediate physical application. The gauge-invariant equations for
linear perturbations were written in the language of the 1+1+2 covariant
decomposition formalism, which relates the geometry of the spacetime and the
properties of the matter fluid with the kinematical quantities that
characterize two sets of congruences: a timelike congruence and a spacelike
congruence. As mentioned above, however, we will relate all quantities
directly with the metric tensor and the matter fluid variables. Nonetheless,
since the 1+1+2 variables have intrinsic physical meaning, they are important
to interpret the results. The reader can find the basic definitions of the
1+1+2 formalism in Appendix A.

The evolution of general perturbations of stars is a compelling yet remarkably
complicated problem in gravitation theory. Here, we will focus on studying
adiabatic, radial perturbations of stars with a perfect fluid source. The
article is intended to be self-contained. Nonetheless, some general properties
will be simply stated, and we redirect the reader to Ref. [9] for technical
details.

###  II.1 The equilibrium spacetime and the perturbation variables

We will assume that the equilibrium background spacetime is static and
spherically symmetric, such that it can be characterized by a line element of
the form

$ds_{0}^{2}=-\left(g_{0}\right)_{tt}dt^{2}+\left(g_{0}\right)_{rr}dr^{2}+r^{2}d\Omega^{2}\,,$
(1)

where $d\Omega^{2}$ represents the natural line element for the unit 2-sphere,
$\left(t,r\right)$ are the standard Schwarzschild coordinates measured by an
observer at spatial infinity, and the metric components
$\left(g_{0}\right)_{tt}$ and $\left(g_{0}\right)_{tt}$ are functions of the
circumferential radius, $r$, only. We will consider the setup where two
solutions of the EFE are smoothly matched at a common timelike hypersurface,
such that the interior of the star is described by a static, spatially compact
solution with a perfect fluid source, with $r\in\left[0,r_{b}\right]$, where
$r=r_{b}$ defines the boundary of the star, and the exterior spacetime, with
$r>r_{b}$, is described by a radial branch of the vacuum Schwarzschild
solution with no event horizons. The metric coefficients
$\left(g_{0}\right)_{tt}$ and $\left(g_{0}\right)_{rr}$ are determined by the
EFE with a zero cosmological constant

$R_{\alpha\beta}-\frac{1}{2}g_{\alpha\beta}R=T_{\alpha\beta}\,,$ (2)

where $R_{\alpha\beta}:=R_{\alpha\mu\beta}{}^{\mu}$ represents the Ricci
tensor, $R:=R_{\mu}{}^{\mu}$ the Ricci scalar, and $T_{\alpha\beta}$ the
metric stress-energy tensor described by a perfect fluid source, i.e.

$T_{\alpha\beta}=\left(\mu_{0}+p_{0}\right)(u_{0})_{\alpha}(u_{0})_{\beta}+p_{0}\left(g_{0}\right)_{\alpha\beta}\,,$
(3)

where $u_{0}$ is the 4-velocity of the elements of volume of the fluid, and
$\mu_{0}$ and $p_{0}$ are, respectively, the energy density and the isotropic
pressure of the matter fluid. Here and in the following, we will use the
subscript “0” to explicitly refer to quantities of the equilibrium spacetime.

It is useful to introduce the following scalar functions

$\displaystyle\phi_{0}$
$\displaystyle=\frac{2}{r\sqrt{\left(g_{0}\right)_{tt}}}\,,$ (4)
$\displaystyle\mathcal{A}_{0}$
$\displaystyle=\frac{1}{2\left(g_{0}\right)_{tt}\sqrt{\left(g_{0}\right)_{rr}}}\frac{d\left(g_{0}\right)_{tt}}{dr}\,,$
$\displaystyle\mathcal{E}_{0}=$
$\displaystyle\frac{1}{3}\mu_{0}+p_{0}-\mathcal{A}_{0}\phi_{0}\,.$

The function $\phi_{0}$ represents the spatial expansion of the normalized
radial gradient vector field, $\mathcal{A}_{0}$ is the radial component of the
4-acceleration of the elements of volume of the fluid, and $\mathcal{E}_{0}$
is the pure radial component of the electric part of the Weyl tensor, Eq.
(53), which partially describes radial tidal forces.

To characterize the perturbed spacetime with respect to the equilibrium
background unambiguously, we have to choose variables that are identification
gauge-invariant, that is, variables that are independent of the choice of
diffeomorphism between the equilibrium and the perturbed spacetimes. Following
the Stewart-Walker Lemma [23], variables that vanish identically in the
background spacetime are gauge-invariant.

Since the equilibrium spacetime is assumed static, the proper time derivatives
of covariantly defined quantities vanish in the background and can be used to
characterize the perturbed spacetime. Indeed, to describe the perturbations,
we will consider the gauge-invariant variables

$\displaystyle\mathsf{m}$ $\displaystyle:=\dot{\mu}\,,$
$\displaystyle\mathsf{p}$ $\displaystyle:=\dot{p}\,,$
$\displaystyle\mathsf{A}$ $\displaystyle:=\dot{\mathcal{A}}\,,$
$\displaystyle\mathsf{F}$ $\displaystyle:=\dot{\phi}\,,$
$\displaystyle\mathsf{E}$ $\displaystyle:=\dot{\mathcal{E}}\,,$ (5)

where the “dot” represents the proper time derivative of an observer locally
comoving with the fluid. In addition to the variables in Eq. (5), we will also
consider two other gauge-invariant variables: the expansion scalar $\theta$
associated with the integral curves of the $u$ vector field, and the
nontrivial radial component of the shear tensor, $\Sigma$ (see Appendix A for
details). The expansion scalar, $\theta$, represents the fractional rate of
change of the sectional volume of the congruence associated with the vector
field $u$ per unit of proper time $\tau$, whereas $\Sigma$, together with
$\theta$, partially characterizes the proper time evolution of the radial
inhomogeneity of the matter fluid. Then, adiabatic, radial perturbations of
perfect fluid stars can be completely described in a frame comoving with the
fluid, by the set of gauge-invariant variables
$\left\\{\mathsf{m},\mathsf{p},\mathsf{A},\mathsf{F},\mathsf{E},\theta,\Sigma\right\\}$.

###  II.2 Harmonic decomposition

As is often the case in relativistic perturbation theory, the equilibrium
solution is found by considering highly symmetric setups. These symmetries can
be taken into account to transform the linearized system of partial
differential equations into a system of ordinary differential equations.
Indeed, in the case of a static, spherically symmetric equilibrium spacetime,
for a suitable choice of frame and at linear perturbation order, all
quantities can be can be written in terms of the Spherical Harmonics, $Y_{\ell
m}$, and the eigenfunctions of the Laplace operator in $\mathbb{R}$,
$e^{i\upsilon\tau}$, such that, a perturbed first-order scalar quantity $\chi$
can be written as

$\chi=\sum_{\upsilon}\left(\sum_{\ell=0}^{+\infty}\sum_{m=-\ell}^{\ell}\Psi_{\chi}^{\left(\upsilon,\ell\right)}Y_{\ell
m}\right)e^{i\upsilon\tau}\,,$ (6)

where $\sum_{\upsilon}$ represents either a discrete sum or an integral in
$\upsilon$, depending on the boundary conditions of the problem, $\tau$ is the
proper time measured by an observer comoving with the fluid in the background
spacetime, $\upsilon$ represent the eigenfrequencies, and
$\Psi_{\chi}^{\left(\upsilon,\ell\right)}$ are the harmonic coefficients,
which depend only on the radial coordinate $r$.

In the case of isotropic perturbations, all coefficients
$\Psi_{\chi}^{\left(\upsilon,\ell\right)}$ with $\ell\geq 1$ are identically
zero, that is, dipole and higher-order angular multipoles must be trivial,
otherwise the perturbations would induce preferred directions in the system.
Moreover, it was rigorously shown in Ref. [9] that, under certain regularity
conditions, for this type of perturbation, the eigenfrequencies are such that
$\upsilon^{2}$ are countable, real, simple, have a minimum, and are unbounded
from above. Therefore, in the study of adiabatic, isotropic perturbations a
gauge-invariant, first-order scalar quantity $\chi$ can be decomposed as

$\chi=\sum_{\upsilon^{2}=\left\\{\upsilon_{0}^{2},\upsilon_{1}^{2},...\right\\}}\Psi_{\chi}^{\left(\upsilon\right)}Y_{00}e^{i\upsilon\tau}\,,$
(7)

where the radial harmonic coefficients depend on the radial coordinate only:
$\Psi_{\chi}^{\left(\upsilon\right)}=\Psi_{\chi}^{\left(\upsilon\right)}\left(r\right)$,
and we have omitted the (trivial) dependency on $\ell$ to lighten the
notation.

Alternatively, instead of considering the proper time $\tau$, we can consider
the time coordinate $t$. In that case, given the equation (see Ref. [9])

$\upsilon\left(r\right)=\frac{\lambda}{\sqrt{\left(g_{0}\right)_{tt}}}\,,$ (8)

relating the eigenfrequencies, $\upsilon$, measured by an observer comoving
with the fluid, with the constant eigenfrequencies $\lambda$, measured by a
static observer at spatial infinity, a gauge-invariant, first-order scalar
quantity $\chi$ can be equivalently given by Eq. (7) or

$\chi=\sum_{\lambda^{2}=\left\\{\lambda_{0}^{2},\lambda_{1}^{2},...\right\\}}\Psi_{\chi}^{\left(\lambda\right)}Y_{00}e^{i\lambda
t}\,.$ (9)

###  II.3 Gauge-invariant equation of state and perturbation equations

The perturbation variables proposed previously in this section completely
describe the dynamical evolution of the adiabatically, radially perturbed
spacetime. Nonetheless, to close the system, we need to characterize the
perturbed matter fluid by providing an equation of state. As a simplifying
assumption, we will consider that the perturbed matter fluid verifies a
barotropic equation of state, such that

$p=f\left(\mu\right)\,,$ (10)

where $f$ is assumed to be twice differentiable in an open neighborhood of
$\mu_{0}$. Then, at linear order

$\mathsf{p}\approx f^{\prime}\left(\mu_{0}\right)\mathsf{m}\,,$ (11)

where prime represents derivative with respect to the function’s parameter, so
that $f^{\prime}\left(\mu_{0}\right)$ represents the square of the adiabatic
speed of sound in the perturbed fluid, to be assumed nonvanishing in the
interior of the perturbed star. Notice that the function $f$ does not have to
be equal to the equation of state of the equilibrium configuration. Indeed,
the equilibrium fluid is not even required to verify a barotropic equation of
state. This freedom is physically relevant, as it was noted, for instance, in
Refs. [6, 24] where distinct adiabatic indexes were considered for the
equilibrium and the perturbed star. Moreover, as we will see, in the case of
an interior Schwarzschild background spacetime, the choice of $f$ is
instrumental in analyzing its stability.

In Appendix B we list the covariant, nontrivial perturbation equations for the
variables
$\left\\{\mathsf{m},\mathsf{p},\mathsf{A},\mathsf{F},\mathsf{E},\theta,\Sigma\right\\}$,
that characterize adiabatic, radial perturbations. Breaking covariance, in the
Schwarzschild coordinate system $\left(t,r\right)$ and considering the
harmonic decomposition described in the previous subsection, the radial
coefficients associated with those variables verify the following system of
differential equations [9]

$\frac{d\Psi_{\mathsf{p}}^{\left(\upsilon\right)}}{dr}+\frac{4\mathcal{A}_{0}}{r\phi_{0}}\left(1+\frac{1}{3f^{\prime}\left(\mu_{0}\right)}\right)\Psi_{\mathsf{p}}^{\left(\upsilon\right)}=-\frac{2\left(\mu_{0}+p_{0}\right)}{r\phi_{0}}\left(\Psi_{\mathsf{A}}^{\left(\upsilon\right)}+\mathcal{A}_{0}\Psi_{\Sigma}^{\left(\upsilon\right)}\right)\,,$
(12)
$\frac{d\Psi_{\mathsf{A}}^{\left(\upsilon\right)}}{dr}+\left(\frac{6\mathcal{A}_{0}}{r\phi_{0}}-\frac{1}{r}\right)\Psi_{\mathsf{A}}^{\left(\upsilon\right)}=\frac{2\mathcal{E}_{0}}{r\phi_{0}\left(\mu_{0}+p_{0}\right)f^{\prime}\left(\mu_{0}\right)}\Psi_{\mathsf{p}}^{\left(\upsilon\right)}-\frac{3}{r\phi_{0}}\left(\upsilon^{2}+\mathcal{A}_{0}^{2}+\frac{1}{3}\mu_{0}-2\mathcal{E}_{0}\right)\Psi_{\Sigma}^{\left(\upsilon\right)}\,,$
(13)
$\displaystyle\frac{d\Psi_{\Sigma}^{\left(\upsilon\right)}}{dr}+\left(\frac{3}{r}-\frac{4\mathcal{A}_{0}}{3r\phi_{0}f^{\prime}\left(\mu_{0}\right)}\right)\Psi_{\Sigma}^{\left(\upsilon\right)}$
$\displaystyle=\frac{2}{3\left(\mu_{0}+p_{0}\right)f^{\prime}\left(\mu_{0}\right)}\left[\left(\frac{f^{\prime\prime}\left(\mu_{0}\right)}{f^{\prime}\left(\mu_{0}\right)}+\frac{1}{\mu_{0}+p_{0}}\right)\frac{d\mu_{0}}{dr}+\frac{1}{r\phi_{0}}\left(\frac{4}{3f^{\prime}\left(\mu_{0}\right)}+2\right)\mathcal{A}_{0}\right]\Psi_{\mathsf{p}}^{\left(\upsilon\right)}$
(14)
$\displaystyle+\frac{4}{3r\phi_{0}f^{\prime}\left(\mu_{0}\right)}\Psi_{\mathsf{A}}^{\left(\upsilon\right)}\,,$

and the constraints

$\left(\upsilon^{2}+\mathcal{A}_{0}\phi_{0}+\mathcal{A}_{0}^{2}-p_{0}\right)\left(\frac{2}{3}\Psi_{\theta}^{\left(\upsilon\right)}-\Psi_{\Sigma}^{\left(\upsilon\right)}\right)=\Psi_{\mathsf{p}}^{\left(\upsilon\right)}-\phi_{0}\Psi_{\mathsf{A}}^{\left(\upsilon\right)}\,,$
(15)

$\Psi_{\mathsf{E}}^{\left(\upsilon\right)}=\mathcal{E}_{0}\left(\frac{3}{2}\Psi_{\Sigma}^{\left(\upsilon\right)}+\frac{\Psi_{\mathsf{p}}^{\left(\upsilon\right)}}{f^{\prime}\left(\mu_{0}\right)\left(\mu_{0}+p_{0}\right)}\right)-\frac{1}{2}\left(\mu_{0}+p_{0}\right)\Psi_{\Sigma}^{\left(\upsilon\right)}\,,$
(16)

$\Psi_{\mathsf{F}}^{\left(\upsilon\right)}=\left(\frac{1}{2}\phi_{0}-\mathcal{A}_{0}\right)\left(\frac{2\Psi_{\mathsf{p}}^{\left(\upsilon\right)}}{3f^{\prime}\left(\mu_{0}\right)\left(\mu_{0}+p_{0}\right)}+\Psi_{\Sigma}^{\left(\upsilon\right)}\right)\,,$
(17)
$\Psi_{\mathsf{m}}^{\left(\upsilon\right)}=-\left(\mu_{0}+p_{0}\right)\Psi_{\theta}^{\left(\upsilon\right)}\,,$
(18)
$\Psi_{\mathsf{p}}^{\left(\upsilon\right)}=f^{\prime}\left(\mu_{0}\right)\Psi_{\mathsf{m}}^{\left(\upsilon\right)}\,,$
(19)

where Eq. (19) follows from Eq. (11). The constraint equation equation (15) is
not propagated, therefore it cannot be used to reduce the size the system of
differential equations (12)–(14). In fact, it is straightforward to show that
Eqs. (12)–(15) imply Eq. (8).

To select the physically acceptable solutions and formalize the boundary value
problem, we impose the following boundary conditions:

  1. (i)

the energy density and the pressure perturbations at the center of the star,
$r=0$, must be finite in a neighborhood of the initial instant;

  2. (ii)

the interior perturbed spacetime can be smoothly matched to an exterior vacuum
Schwarzschild spacetime at a timelike hypersurface, the boundary of the star.

From the point of view of the comoving observer, the boundary condition (ii)
implies that the pressure of the perturbed fluid is identically zero at all
times at the hypersurface and so, $\mathsf{p}$, hence
$\Psi_{\mathsf{p}}^{\left(\upsilon\right)}$, is also identically zero at the
boundary, that is

$\Psi_{\mathsf{p}}^{\left(\upsilon\right)}\left(r_{b}\right)=0\,.$ (20)

Before proceeding, we remark that, in general, the coordinate system in the
background spacetime and that of the perturbed spacetime are not necessarily
the same: any smooth mapping can be considered. Since the perturbation
variables are gauge-invariant, that choice does not affect the results. In the
particular case of isotropic perturbations, the Schwarzschild coordinate
system can always be defined since there is no gravitational wave emission,
and the spacetime is asymptotically flat, hence the time coordinate $t$ and
circumferential radius $r$ are equally defined by an observer at spatial
infinity in both spacetimes. In particular, the adoption of this coordinate
system is useful to compare our approach with the classical results of
Chandrasekhar, found from metric based perturbation theory [1, 2]. In Appendix
C we show explicitly that Chandrasekhar’s second-order radial pulsation
equation follows from the system above by relating the kinematical quantities
with the radial displacement parameter and its derivatives.

###  II.4 Analytic solutions

The system (12)–(19) with boundary conditions (i) and (ii) completely
characterizes adiabatic, radial perturbations of a star composed of a perfect
fluid. Specifying the background spacetime and the equation of state of the
perturbed fluid, numerical methods can be used to find approximate solutions.
Nonetheless, contrary to the original form of the second-order Chandrasekhar’s
pulsation equation [1, 2], or the associated first-order realizations of Refs.
[4, 6], it is possible, under rather general conditions, to find analytic
solutions for the perturbations using standard theory of systems of linear
ordinary differential equations.

To find analytic solutions for the system (12)–(19), we will further impose
that the following regularity constraints hold:

  * •

the equilibrium fluid verifies the weak energy condition;

  * •

the background spacetime is a solution of the Tolman-Oppenheimer-Volkoff
equation for real analytic, nontrivial energy density and pressure functions
for the whole range within the equilibrium star;

  * •

the square of the speed of sound of the perturbed fluid, $f^{\prime}$, is
positive and real analytic in the interior and at the boundary of the star.

Requiring real analytical background solutions is a rather strong constraint.
However, to our knowledge, all classical exact solutions for compact
astrophysical objects verify this hypothesis in some open neighborhood of the
center of the star ($r=0$). Hence, the following results are appropriate for
treating perturbations of physically relevant setups. The radius of
convergence of the power series, of course, may or may not be greater than the
radius of the equilibrium star. In the latter case, further treatment has to
be carried out to extend the solution to the boundary. Nonetheless, to prove
the following results, it suffices that the radius of convergence of the power
series around the center is nonzero.

Imposing the above conditions, the background EFE equations imply

$\displaystyle\phi_{0}$
$\displaystyle=\frac{2}{r}\sqrt{1-\frac{2M\left(r\right)}{r}}\,,$ (21)
$\displaystyle\mathcal{E}_{0}$
$\displaystyle=\frac{1}{3}\mu_{0}-\frac{2M\left(r\right)}{r^{3}}\,,$
$\displaystyle\mathcal{A}_{0}\phi_{0}$
$\displaystyle=p_{0}+\frac{2M\left(r\right)}{r^{3}}\,,$

where

$M\left(r\right):=\frac{1}{2}\int_{0}^{r}\mu_{0}x^{2}dx\,,$ (22)

is usually called the mass function. If the functions $\mu_{0}$ and $p_{0}$
verify the weak energy condition and are real analytic within the star, so are
the functions $\mathcal{A}_{0}$ and $\mathcal{E}_{0}$, and $\phi_{0}$ has a
simple pole at the center, $r=0$, but is otherwise real analytic in the
interior and boundary of the star. Therefore, in the considered setup, Eqs.
(12)–(14) form a system of ordinary differential equations with real analytic
coefficients around $r=0$, with a simple pole at $r=0$, and solutions can be
found around the singular point.

Consider the system (12)–(14) in matrix form:

$\frac{d\mathds{W}}{dr}=\left(r^{-1}\mathds{R}+\Theta\right)\mathds{W}\,,$
(23)

where

$\displaystyle\mathds{W}$
$\displaystyle=\left[\begin{array}[]{c}\Psi_{\mathsf{p}}^{\left(\lambda\right)}\\\
\Psi_{\mathsf{A}}^{\left(\lambda\right)}\\\
\Psi_{\Sigma}^{\left(\lambda\right)}\end{array}\right]\,,$
$\displaystyle\mathds{R}$ $\displaystyle=\left[\begin{array}[]{lcr}0&0&0\\\
0&1&0\\\ 0&0&-3\end{array}\right]\,,$ (24) $\displaystyle\Theta$
$\displaystyle=\frac{2}{r\phi_{0}}\left[\begin{array}[]{lcr}-2\mathcal{A}_{0}\left(1+\frac{1}{3f^{\prime}\left(\mu_{0}\right)}\right)&\quad-\left(\mu_{0}+p_{0}\right)&\quad-\left(\mu_{0}+p_{0}\right)\mathcal{A}_{0}\\\
\\\
\frac{\mathcal{E}_{0}}{f^{\prime}\left(\mu_{0}\right)\left(\mu_{0}+p_{0}\right)}&\quad-3\mathcal{A}_{0}&\quad-\frac{3}{2}\left(\upsilon^{2}+\mathcal{A}_{0}^{2}+\frac{1}{3}\mu_{0}-2\mathcal{E}_{0}\right)\\\
\\\
\frac{3f^{\prime\prime}\left(\mu_{0}\right)r\phi_{0}\partial_{r}\mu_{0}+4\mathcal{A}_{0}}{9\left(\mu_{0}+p_{0}\right)\left[f^{\prime}\left(\mu_{0}\right)\right]^{2}}+\frac{r\phi_{0}\partial_{r}\mu_{0}+2\mathcal{A}_{0}\left(\mu_{0}+p_{0}\right)}{3\left(\mu_{0}+p_{0}\right)^{2}f^{\prime}\left(\mu_{0}\right)}&\quad\frac{2}{3f^{\prime}\left(\mu_{0}\right)}&\quad\frac{2\mathcal{A}_{0}}{3f^{\prime}\left(\mu_{0}\right)}\end{array}\right]\,.$
(25)

The regularity conditions that we have imposed on the thermodynamic variables
of the equilibrium configuration, $\mu_{0}$ and $p_{0}$, and on the equation
of state of the perturbed fluid, $f$, guarantee that the matrix $\Theta$ is
real analytic at $r=0$, in particular, the fact that $\mu_{0}$ and $p_{0}$ are
differentiable functions guarantees that $r\phi_{0}$ does not vanish in the
interior of the star. From Eq. (23) we conclude that $r=0$ is a regular
singular point of the system. However, of course, this does not imply that all
solutions must be singular at the center. To select only the physically
acceptable solutions, we will impose the boundary conditions (i) and (ii).

To solve the system of ordinary differential equations, we will follow the
formalism in Ref. [25] and find solutions in a power series form. Since the
$\Theta$ matrix is assumed to be real analytic at the center of the star, it
can be expanded in a convergent power series of the form

$\Theta\left(r\right)=\sum_{n=0}^{+\infty}\Theta_{n}r^{n}\,.$ (26)

Then, the solution matrix $\mathds{W}$ can be written in a power series,
guaranteed to converge to the solution in a neighborhood of $r=0$. The radius
of convergence of the power series solution is equal, except for possibly at
$r=0$, to the radius of convergence of the power series of $\Theta$, which, of
course, depends on the considered background spacetime.

Before presenting the formal solutions, we remark that for a general static,
spherically symmetric spacetime, the general form for the perturbative
solutions can be rather complicated, however, considering the regularity
conditions imposed on $\mu_{0}$ and $p_{0}$, and assuming
$f^{\prime}\left(\mu_{0}\right)$ is not zero in a neighborhood of the center
of the star, using the Tolman-Oppenheimer-Volkoff equations it can be shown
that various entries of the coefficient-matrices
$\left\\{\Theta_{0},\Theta_{1},\Theta_{2},\Theta_{3}\right\\}$ are zero, which
greatly simplifies the general family of solutions of the system. Then, taking
into consideration the regularity of the background solution, applying the
method in Ref. [25] yields

$\displaystyle\left[\begin{array}[]{c}\Psi_{\mathsf{p}}^{\left(\lambda\right)}\\\
\Psi_{\mathsf{A}}^{\left(\lambda\right)}\\\
\Psi_{\Sigma}^{\left(\lambda\right)}\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{rcr}-1&\frac{12}{r}\left[\left(\Theta_{0}\right)_{12}\left(\Theta_{0}\right)_{23}-3\left(\Theta_{1}\right)_{13}\right]&\>0\\\
0&\>12\left(\Theta_{0}\right)_{23}\left[\left(\Theta_{1}\right)_{22}-\left(\Theta_{1}\right)_{33}\right]-18\left(\Theta_{2}\right)_{23}+4\left[\left(\Theta_{0}\right)_{23}\right]^{2}\left(\Theta_{0}\right)_{32}-\frac{12}{r^{2}}\left(\Theta_{0}\right)_{23}&r\\\
0&\frac{36}{r^{3}}&0\end{array}\right]\mathds{P}_{\mathds{W}}\left[\begin{array}[]{c}c_{1}\\\
c_{2}\\\ c_{3}\end{array}\right]\,,$ (27)

where $c_{1}$, $c_{2}$ and $c_{3}$ are integration constants, which might
depend on the constant eigenfrequencies $\lambda$, defined in Eq. (8). The
notation $\left(\Theta_{n}\right)_{ij}$ is to be interpreted as the $ij$-entry
of the $n$th order matrix coefficient of the power series expansion of
$\Theta$, Eq. (26), and $\mathds{P}_{\mathsf{W}}$ is a real analytic matrix,
such that

$\displaystyle\mathds{P}_{\mathds{W}}\left(r\right)$
$\displaystyle=\sum_{n=0}^{+\infty}\mathds{P}_{n}r^{n}\,,$ (28)
$\displaystyle\mathds{P}_{0}$ $\displaystyle=\mathds{I}_{3}\,,$
$\displaystyle\mathds{P}_{k}$
$\displaystyle=\frac{1}{k}\sum_{j=0}^{k-1}\mathds{A}_{k-1-j}\mathds{P}_{j}\,,\quad\text{for
}k\geq 1\,,$

where $\mathds{I}_{3}$ represents the $3\times 3$ identity matrix, the matrix
$\mathds{A}$ is real analytic with the same radius of convergence as $\Theta$,
and $\mathds{A}_{n}$ represents the $n$th order matrix-coefficient of its
power series expansion, that is,
$\mathds{A}\left(r\right)=\sum_{n=0}^{+\infty}\mathds{A}_{n}r^{n}$. Due to its
size, the matrix $\mathds{A}$ is presented in Appendix D.

Using $\mathds{P}_{0}=\mathds{I}_{3}$ in Eq. (27), we can directly compute the
lower order coefficients of the power series expansion of $\mathds{W}$.
Imposing the boundary condition at the center sets the coefficient $c_{2}$ to
be zero, otherwise the pressure would diverge at $r=0$ at all times. Then, we
find

$\left[\begin{array}[]{c}\Psi_{\mathsf{p}}^{\left(\lambda\right)}\\\
\Psi_{\mathsf{A}}^{\left(\lambda\right)}\\\
\Psi_{\Sigma}^{\left(\lambda\right)}\end{array}\right]=\left[\begin{array}[]{c}-c_{1}+\mathcal{O}\left(r^{2}\right)\\\
c_{3}r+\mathcal{O}\left(r^{3}\right)\\\
\mathcal{O}\left(r^{2}\right)\end{array}\right]\,.$ (29)

We see that the coefficient $c_{1}$ directly characterizes the behavior of
$\Psi_{\mathsf{p}}^{\left(\lambda\right)}$ at $r=0$ and that both
$\Psi_{\mathsf{A}}^{\left(\lambda\right)}$ and
$\Psi_{\Sigma}^{\left(\lambda\right)}$ must vanish at the center. On the other
hand, the coefficients $c_{1}$ and $c_{3}$ are not independent: considering
the regularity of the background spacetime and imposing the constraints (15),
(18) and (19) leads to

$c_{3}\stackrel{{\scriptstyle
r=0}}{{=}}-\frac{c_{1}}{3\left(\mu_{0}+p_{0}\right)f^{\prime}\left(\mu_{0}\right)}\left[\frac{\lambda^{2}}{\left(g_{0}\right)_{tt}}+\frac{1}{3}\mu_{0}+\frac{3}{2}\left(\mu_{0}+p_{0}\right)f^{\prime}\left(\mu_{0}\right)\right]\,,$
(30)

where all quantities on the right-hand side are to be evaluated at $r=0$.
Therefore, for each value of $\lambda^{2}$ there is a single arbitrary
parameter, either $c_{1}$ or $c_{3}$, to be characterized by the initial
perturbation at the center of the star. That is, the coefficients of the
Fourier transform of the initial perturbation provide the values of the
independent parameters, setting which eigenmodes are excited and the
respective initial magnitude. In what follows we will consider $c_{1}$ as the
independent parameter. Specifying the background spacetime, the equation of
state of the perturbed fluid and the values of the eigenfrequencies $\lambda$,
these results allow us to find regular analytic solutions for the
perturbations that verify the boundary conditions.

In addition to the analytical results for the eigenfunctions, in Ref. [9] it
was possible to establish lower bounds for the square of the eigenfrequencies,
$\lambda^{2}$. In particular, if the regularity conditions in the beginning of
this subsection and the boundary conditions (i) and (ii) hold, nontrivial
$\mathcal{C}^{1}$ solutions of the system (12)–(19) exist only if

$\lambda^{2}\max_{r\in\left]0,r_{b}\right[}\left(g_{0}\right)_{tt}>-\max_{r\in\left]0,r_{b}\right[}\left[\frac{\mu_{0}+p_{0}}{\phi_{0}}\left(\frac{1}{2}\phi_{0}+2\mathcal{A}_{0}\right)+\frac{\mathcal{A}_{0}^{2}}{f^{\prime}\left(\mu_{0}\right)}+\frac{r\phi_{0}\mathcal{A}_{0}}{2\left(\mu_{0}+p_{0}\right)}\frac{d\mu_{0}}{dr}\right]\,.$
(31)

This result can be useful to determine numerically the eigenfrequencies of the
system, offering a baseline to search for their values.

##  III Perturbations and fundamental eigenfrequencies of classic exact
solutions

As illustrative examples of the general results in the previous section, we
will study the properties of adiabatic, radial perturbations of some classical
solutions of perfect fluid stars within general relativity. Since the
independent parameter $c_{1}$ simply characterizes the magnitude of a specific
eigenfunction at $r=0$, without loss of generality, in this section we will
set $c_{1}=-1$.

In Table 1, considering a line element of the form (1), we present the
nontrivial metric coefficients of some classical, physically relevant
solutions of the Einstein field equations. In Table 2 we present the absolute
value of the first three eigenfrequencies for specific models for the
equilibrium background spacetime. For the interior Schwarzschild solution, at
linear order, the perturbed fluid is completely characterized by a constant
speed of sound, which must be provided as an extra parameter. For the other
models, we assume that the equation of state of the perturbed fluid is the
same as that of the equilibrium setup. For all considered models, the
eigenfrequencies take real values, therefore all equilibrium spacetimes
represent stable configurations under adiabatic, radial perturbations.

In Figures 1–4 we present the radial profile of the Fourier coefficients of
the functions $\mathsf{p}$, $\mathsf{A}$, and $\Sigma$, associated with the
eigenfrequencies presented in Table 2 for the various background spacetimes.
Figures 1–4 highlight the expected behavior for the eigenfunctions. For a
real-analytic background spacetime, in Ref. [9] it was shown that the
perturbations equations can be cast in the form of a Sturm-Liouville
eigenvalue problem. Therefore, in particular, the number of roots of the
eigenfunctions is associated with the order of the associated eigenvalue in
the sequence $\left(\lambda^{2}_{n}\right)_{n\in\mathbb{N}}$ .

Spacetime | Non-trivial metric components
---|---
Interior Schwarzschild | $\begin{aligned} \\\ \left(g_{0}\right)_{tt}=&\left(3\sqrt{1-\frac{2M}{r_{b}}}-\sqrt{1-\frac{2Mr^{2}}{r_{b}^{3}}}\right)^{2}\\\ \left(g_{0}\right)_{rr}=&\left(1-\frac{2Mr^{2}}{r_{b}^{3}}\right)^{-1}\\\ \\\ \end{aligned}$
Tolman IV | $\begin{aligned} \\\ \left(g_{0}\right)_{tt}=&B^{2}\left(\frac{r^{2}}{A^{2}}+1\right)\\\ \left(g_{0}\right)_{rr}=&\frac{\frac{2r^{2}}{A^{2}}+1}{\left(1+\frac{r^{2}}{A^{2}}\right)\left(1-\frac{r^{2}}{R^{2}}\right)}\\\ \\\ \end{aligned}$
Kuch2 III | $\begin{aligned} \\\ \,\,\left(g_{0}\right)_{tt}=&Be^{\frac{Ar^{2}}{2}}\\\ \,\,\left(g_{0}\right)_{rr}=&\left(r^{2}e^{-\frac{1}{2}Ar^{2}}\left[C-\frac{A}{2e}\text{Ei}\left(\frac{Ar^{2}}{2}+1\right)\right]+1\right)^{-1}\\\ \\\ \end{aligned}$
Heint IIa | $\begin{aligned} \\\ \left(g_{0}\right)_{tt}=&A^{2}\left(ar^{2}+1\right)^{3}\\\ \left(g_{0}\right)_{rr}=&\left(1-\frac{3ar^{2}\left[c\left(4ar^{2}+1\right)^{-\frac{1}{2}}+1\right]}{2\left(ar^{2}+1\right)}\right)^{-1}\\\ \\\ \end{aligned}$
Table 1: Metric coefficients of classical solutions of the Einstein field equations. The spacetimes are assumed to be characterized by a line element of the form of Eq. (1). We follow the naming conventions for the solutions of Ref. [22]. Spacetime | Parameters | $\left|\lambda_{0}\right|$ | $\left|\lambda_{1}\right|$ | $\left|\lambda_{2}\right|$
---|---|---|---|---
Interior Schwarzschild | $\,\,\left(M,r_{b},c_{s}^{2}\right)=\left(0.1,1,0.1\right)\,\,$ | 0.108 | 0.305 | 0.478
Tolman IV | $\left(A,B,R\right)=\left(1,1,1.5\right)$ | 1.533 | 4.281 | 6.723
Kuch2 III | $\left(A,B,C\right)=\left(5,1,-3\right)$ |  20.214 |  41.085 |  61.808
Heint IIa | $\left(a,A,C\right)=\left(1,1,1.5\right)$ | 4.004 | 10.262 | 15.939
Table 2: First absolute values of the eigenfrequencies, $\lambda$, rounded to
three decimal places, for the equilibrium solutions in Table 1 for specific
values of the spacetime parameters. In all examples, the eigenfrequencies are
real, hence all spacetimes are stable under adiabatic, radial perturbations.
We follow the naming conventions for the solutions of Ref. [22].

Except for the interior Schwarzschild solution, the results in Table 2 were
compared with the predictions of the systems in Refs. [4, 6]. Implementing a
shooting method to solve numerically each of those systems for each
equilibrium spacetime, all values for the fundamental eigenfrequencies matched
to the considered numerical accuracy. Moreover, the general solutions found in
Sec. II are regular around the center, and for the spacetime parameters in
Table 2, the solutions are exact and the radius of convergence of the power
series is greater than the radius of the star. Hence, we can explicitly
evaluate the boundary conditions at any point within the star, in particular
at the center and the surface.

Figure 1: Radial profile of the Fourier coefficients of the functions
$\mathsf{p}$, $\mathsf{A}$, and $\Sigma$, associated with the eigenfrequencies
presented in Table 2 for the interior Schwarzschild spacetime. For all
oscillation modes, it was assumed $c_{1}=-1$. Figure 2: Radial profile of the
Fourier coefficients of the functions $\mathsf{p}$, $\mathsf{A}$, and
$\Sigma$, associated with the eigenfrequencies presented in Table 2 for the
Tolman IV spacetime. For all oscillation modes, it was assumed $c_{1}=-1$.
Figure 3: Radial profile of the Fourier coefficients of the functions
$\mathsf{p}$, $\mathsf{A}$, and $\Sigma$, associated with the eigenfrequencies
presented in Table 2 for the Kuchowicz 2-III spacetime. For all oscillation
modes, it was assumed $c_{1}=-1$. Figure 4: Radial profile of the Fourier
coefficients of the functions $\mathsf{p}$, $\mathsf{A}$, and $\Sigma$,
associated with the eigenfrequencies presented in Table 2 for the Heintzmann
IIa spacetime. For all oscillation modes, it was assumed $c_{1}=-1$.

##  IV Bound on maximum compactness

The dynamical stability of an equilibrium static, self-gravitating fluid is a
crucial problem in astrophysics. Finding that a given solution of the Einstein
field equations is dynamically unstable for some type of perturbation implies
that such a solution might not be suitable to describe massive compact objects
that are expected to be empirically observable in the current universe. In
that regard, considerable effort has been devoted to finding an upper bound
for the maximum compactness of a star: the ratio between the gravitational
mass and the circumferential radius of a star, beyond which it is unstable
(cf., e.g., Refs. [26, 27] for a review for perfect and non-perfect fluids,
and Ref. [28] and references therein for the case of electrically charged
self-gravitating fluids). Nonetheless, the results in the literature are
fundamentally connected to the equation of state of the equilibrium fluid and,
to our knowledge, no universal upper bound has been found for the maximum
compactness of these types of objects. In the particular case of perfect fluid
stars, we can use the previous results and the interior Schwarzschild solution
to conjecture an equation of state agnostic upper bound for the dynamical
stability of this type of objects.

The interior Schwarzschild spacetime represents a solution of general
relativity describing the interior of a star composed of a perfect fluid with
a constant energy density that can be smoothly matched with an exterior vacuum
Schwarzschild spacetime. This solution can be characterized by a line element
of the form (1), with

$\displaystyle\left(g_{0}\right)_{tt}=$
$\displaystyle\left(3\sqrt{1-\frac{2M}{r_{b}}}-\sqrt{1-\frac{2Mr^{2}}{r_{b}^{3}}}\right)^{2}\,,$
(32) $\displaystyle\left(g_{0}\right)_{rr}=$
$\displaystyle\left(1-\frac{2Mr^{2}}{r_{b}^{3}}\right)^{-1}\,,$

where $M$ represents the ADM mass and $r_{b}$ the value of the circumferential
radius at the boundary of the star. Although the interior Schwarzschild
solution does not represent a physically reasonable configuration, this
solution is important from a conceptual point of view since it saturates the
Buchdahl bound for the maximum compactness:

$\frac{M}{r_{b}}\leq\frac{4}{9}=0.\overline{4}\,,$ (33)

hence, it can be thought of as the extreme scenario for a static, self-
gravitating perfect fluid. On this basis, this solution can be used to
conjecture the maximum compactness of a star composed of a perfect fluid,
beyond which it becomes dynamically unstable.

Considering the results of the previous section, we can study the dynamical
stability of the interior Schwarzschild solution by computing the fundamental
eigenfrequency for various values of the compactness parameter
$\frac{M}{r_{b}}$. Notice that, for an interior Schwarzschild background
spacetime with a specific value of the compactness parameter, the value of
$r_{b}$ does not affect the values of the eigenfrequencies, but only the
radial profile of the perturbation variables. Since the energy density of the
background is constant, at a linear level, the square of the speed of sound,
$f^{\prime}\left(\mu_{0}\right)$, is a constant. Then, to infer the maximum
compactness of a physical star described by the interior Schwarzschild
solution, we can impose the causality condition, and consider the extreme
scenario where $f^{\prime}\left(\mu_{0}\right)=c_{s}^{2}=1$, that is, the
square of the speed of sound of the perturbed star is the vacuum speed of
light. Applying the previous results, we find that

$M/r_{b}=0.367\Rightarrow\lambda_{0}^{2}>0\,,\qquad
M/r_{b}=0.368\Rightarrow\lambda_{0}^{2}<0\,.$ (34)

Following this reasoning, we extrapolate that a static, spherically symmetric
solution of the Einstein field equations with a perfect fluid source is
dynamically unstable if

$\frac{M}{r_{b}}\gtrsim 0.368\,.$ (35)

The accuracy of the above estimate, of course, can be increased, but this
level of accuracy is unlikely to be achieved experimentally. The result in Eq.
(35) is significantly smaller than the Buchdahl bound, Eq. (33). Moreover, it
is slightly higher, but in line with the estimates of Refs. [29, 30] for the
maximum compactness, found by considering an affine fluid model and imposing
the hypothesis that such a model yields the most compact star composed of a
perfect fluid verifying a barotropic equation of state. The analysis in this
section, however, is independent of the equation of state of the equilibrium
star and follows from simply considering the extreme case of the causality
condition for the perturbed fluid. Consequently, this upper bound for
instability is universal and it is not expected to be saturated by any perfect
fluid star solution verifying a physically meaningful equation of state.

##  V Conclusion

We have presented a system of first order differential equations to describe
general linear adiabatic, radial perturbations of spatially compact, static,
spherically symmetric solutions of general relativity with a perfect fluid
source. The new results do not rely on auxiliary test functions nor on the
introduction of a singular point at the boundary of the star. Contrary to
previous approaches, assuming some regularity conditions for the equilibrium
spacetime, the system can be solved analytically, finding general solutions
for the perturbation variables.

The results were then used to study adiabatic, radial perturbations of
classical exact solutions of the Einstein field equations, computing the first
eigenfrequencies for particular values of the spacetime parameters. For the
considered models, the equilibria proved to be stable. We have also plotted
the eigenfunctions associated with those eigenfrequencies, illustrating their
radial profile. At first glance, the plots in Figures 1–4 look remarkably
similar. However, this is a product of the very nature of the Sturm-Liouville
problem, the constraints imposed by the boundary conditions, and the
normalized radial coordinate we have used.

We have considered only a small subset of physically relevant exact solutions,
leaving aside some important solutions that have been used in the past to
study physically meaningful scenarios. As was discussed in Ref. [9], although
the thermodynamical description of the perturbed matter fluid is simpler in a
comoving frame, this leads to extra complexity to describe its dynamics such
that we end up with a three-by-three system of differential equations with a
constraint, Eqs. (12)–(15). Indeed, the formulation of the problem can be made
more computationally efficient. By changing the frame of reference, we can
significantly simplify the problem and efficiently study adiabatic, radial
perturbations of more complex background solutions. In other frames, however,
the stress-energy tensor is non-diagonal, containing fluxes and anisotropic
components. Therefore, great care has to be taken to ascertain what is meant
by adiabatic perturbations in the new frame. This will be done elsewhere.

The newly found system of equations also allowed us to conjecture,
independently of the information on the equation of state of the equilibrium
fluid, the upper bound $M/r_{b}\approx 0.368$ for the maximum compactness of a
stable, static, self-gravitating perfect fluid. The result relies on the
hypothesis that the interior Schwarzschild solution represents the extreme
case for a perfect fluid star. Hence, the maximum bound for the compactness
should be absolute, in the sense that no physically meaningful solution of
general relativity is expected to saturate it. Nonetheless, the result is a
tighter upper bound compared to the previously suggested bounds in the
literature, found by considering a specific barotropic equation of state. We
remark, however, that the reasoning for this result omits various important
issues on the stability of compact stellar objects. One example concerns the
evolutionary timescales. The eigenfrequencies
$\left\\{\lambda_{n}\right\\}_{n\in\mathbb{N}}$ can be used to determine the
oscillation period or the period of e-folds of each mode. In the case of
$\lambda_{0}\in\mathbb{R}$, the star will oscillate, and the period of
oscillation for the fundamental mode is given by
$T=2\pi/\left|\lambda_{0}\right|$. On the other hand, if $\lambda_{0}$ is
purely imaginary, the perturbations will not be bounded and the quantity
$T=2\pi/\left|i\lambda_{0}\right|$ represents the period of e-folds. Since the
fundamental mode has the shortest growth timescale, we can assess the growth
timescale of the radial instabilities in the considered idealized setup,
considering this mode. Nonetheless, the evolutionary timescale of compact
stellar objects depends on various factors omitted in our analysis, where we
have assumed the idealized scenario of an adiabatic perturbation. Indeed, the
thermodynamic evolution of the matter fields plays a predominant role in the
evolution of those objects, such that even if the star is pulsationally
unstable, the radial modes might not have time to grow and affect the
structure. Moreover, the previous discussion does not consider the issue of
quasi-stability. In principle, a dynamically unstable object, but such that
its evolution would have a characteristic time comparable to the timescale of
the universe, might be effectively considered stable. A realistic analysis of
the topic should take into consideration all these aspects. The presented
upper bound for dynamical stability should then be valid in an idealized
scenario, disregarding any other effects.

###### Acknowledgements.

PL acknowledges financial support provided under the European Union’s H2020
ERC Advanced Grant “Black holes: gravitational engines of discovery” grant
agreement no. Gravitas–101052587. Views and opinions expressed are, however,
those of the author only and do not necessarily reflect those of the European
Union or the European Research Council. Neither the European Union nor the
granting authority can be held responsible for them. PL acknowledges support
from the Villum Investigator program supported by the VILLUM Foundation (grant
no. VIL37766) and the DNRF Chair program (grant no. DNRF162) by the Danish
National Research Foundation. This project has received funding from the
European Union’s Horizon 2020 research and innovation programme under the
Marie Sklodowska-Curie grant agreement No 101007855 and No 101131233. The work
of SC has been carried out in the framework of activities of the INFN Research
Project QGSKY.

##  Appendix A The 1+1+2 decomposition

###  A.1 Projectors and the Levi-Civita volume form

Let $\left(\mathcal{M},g\right)$ be a Lorentzian manifold of dimension 4,
where $g$ represents the metric tensor, admitting in some open neighborhood
the existence of a congruence of timelike curves with tangent vector field
$u$. Let the congruence be affinely parameterized and
$u_{\alpha}u^{\alpha}=-1$. Without loss of generality, the manifold can be
foliated by 3-surfaces, $V$, pointwise orthogonal to the curves of the
congruence, such that all tensorial quantities are defined by their behavior
along the direction of $u$ in $V$. This formalism is called 1+3 spacetime
covariant decomposition. This decomposition of the spacetime manifold relies
on the existence of a pointwise projector to the cotangent space of $V$ to be
defined as

$h_{\alpha\beta}=g_{\alpha\beta}+u_{\alpha}u_{\beta}\,,$ (36)

where $g_{\alpha\beta}$ represents the components of the metric tensor in some
local coordinate system. The projector operator verifies the following
properties

$\displaystyle h_{\alpha\beta}$ $\displaystyle=h_{\beta\alpha}\,,$
$\displaystyle h_{\alpha\beta}h^{\beta\gamma}$
$\displaystyle=h_{\alpha}{}^{\gamma}\,,$ (37) $\displaystyle
h_{\alpha\beta}u^{\alpha}$ $\displaystyle=0\,,$ $\displaystyle
h_{\alpha}{}^{\alpha}$ $\displaystyle=3\,.$

The 1+1+2 decomposition builds from the 1+3 decomposition by defining a
spacelike congruence with tangent vector field $e$, such that a tensorial
quantity defined in the submanifold $V$ is characterized by its behavior along
$e$ and the 2-surfaces $W$, orthogonal to both $u$ and $e$ at each point.
Without loss of generality, let the spacelike congruence to be affinely
parameterized and $e_{\alpha}e^{\alpha}=1$. We can define a projector onto the
cotangent space of $W$ by

$N_{\alpha\beta}=h_{\alpha\beta}-e_{\alpha}e_{\beta}\,,$ (38)

verifying

$\displaystyle N_{\alpha\beta}$ $\displaystyle=N_{\beta\alpha}\,,$
$\displaystyle N_{\alpha\beta}N^{\beta\gamma}$
$\displaystyle=N_{\alpha}{}^{\gamma}\,,$ (39) $\displaystyle
N_{\alpha\beta}u^{\alpha}$ $\displaystyle=N_{\alpha\beta}e^{\alpha}=0\,,$
$\displaystyle N_{\alpha}{}^{\alpha}$ $\displaystyle=2\,.$

It is also useful to introduce the following skew-symmetric tensors derived
from the covariant Levi-Civita tensor $\varepsilon_{\alpha\beta\gamma\sigma}$,

$\displaystyle\varepsilon_{\alpha\beta\gamma}$
$\displaystyle=\varepsilon_{\alpha\beta\gamma\sigma}u^{\sigma}\,,$ (40)
$\displaystyle\varepsilon_{\alpha\beta}$
$\displaystyle=\varepsilon_{\alpha\beta\gamma}e^{\gamma}\,.$

In what follows, we will adopt the convention to indicate the symmetric and
anti-symmetric part of a tensor using parentheses and brackets, such that for
a 2-tensor $\chi$

$\displaystyle\chi_{\left(\alpha\beta\right)}$
$\displaystyle=\frac{1}{2}\left(\chi_{\alpha\beta}+\chi_{\beta\alpha}\right)\,,$
$\displaystyle\chi_{\left[\alpha\beta\right]}$
$\displaystyle=\frac{1}{2}\left(\chi_{\alpha\beta}-\chi_{\beta\alpha}\right)\,.$
(41)

###  A.2 Covariant derivatives of $u$ and $e$

Using the projector operators $h$ and $N$, the covariant derivatives of the
tangent vector fields $u$ and $e$ can be uniquely decomposed at each point in
their components along $u$, $e$ and in $W$. This decomposition yields the so-
called kinematical quantities of the congruences formed by the integral curves
of $u$ and $e$, providing a clear geometric and physical interpretation of the
behavior of the congruence.

Given a tensor quantity $\chi$, it is useful to introduce the compact notation

$\displaystyle\dot{\chi}$ $\displaystyle:=u^{\mu}\nabla_{\mu}\chi\,,$
$\displaystyle\widehat{\chi}$ $\displaystyle:=e^{\mu}\nabla_{\mu}\chi\,,$ (42)

to represent the derivative along the integral curves of the vector field $u$
and the derivative along the integral curves of the vector field $e$,
respectively.

In the 1+1+2 formalism, the covariant derivative of the tensor field $u$ can
be decomposed as

$\nabla_{\alpha}u_{\beta}=-u_{\alpha}\left(\mathcal{A}e_{\beta}+\mathcal{A}_{\beta}\right)+\frac{1}{3}h_{\alpha\beta}\theta+\sigma_{\alpha\beta}+\omega_{\alpha\beta}\,,$
(43)

where

$\displaystyle\mathcal{A}$
$\displaystyle=-u_{\mu}u^{\nu}\nabla_{\nu}e^{\mu}\,,$
$\displaystyle\mathcal{A}_{\alpha}$
$\displaystyle=N_{\alpha\mu}\dot{u}^{\mu}\,,$ (44)

and the quantities $\theta$, $\sigma_{\alpha\beta}$ and $\omega_{\alpha\beta}$
are the kinematical quantities of the congruence of the integral curves of
$u$. Namely, $\theta$ is the expansion scalar, $\sigma_{\alpha\beta}$ is the
shear tensor, and $\omega_{\alpha\beta}$ the vorticity tensor, defined as

$\displaystyle\theta$ $\displaystyle=h^{\mu\nu}\nabla_{\mu}u_{\nu}\,,$ (45)
$\displaystyle\sigma_{\alpha\beta}$
$\displaystyle=\left(\frac{h_{\alpha}{}^{\mu}h_{\beta}{}^{\nu}+h_{\alpha}{}^{\nu}h_{\beta}{}^{\mu}}{2}-\frac{1}{3}h_{\alpha\beta}h^{\mu\nu}\right)\nabla_{\mu}u_{\nu}=\Sigma_{\alpha\beta}+2\Sigma_{(\alpha}e_{\beta)}+\Sigma\left(e_{\alpha}e_{\beta}-\frac{1}{2}N_{\alpha\beta}\right)\,,$
$\displaystyle\omega_{\alpha\beta}$
$\displaystyle=\frac{1}{2}\left(h_{\alpha}{}^{\mu}h_{\beta}{}^{\nu}-h_{\alpha}{}^{\nu}h_{\beta}{}^{\mu}\right)\nabla_{\mu}u_{\nu}=\varepsilon_{\alpha\beta\mu}\left(\Omega
e^{\mu}+\Omega^{\mu}\right)\,,$

with

$\displaystyle\Sigma_{\alpha\beta}$
$\displaystyle=\left(\frac{N_{\alpha}{}^{\mu}N_{\beta}{}^{\nu}+N_{\alpha}{}^{\nu}N_{\beta}{}^{\mu}}{2}-\frac{1}{2}N_{\alpha\beta}N^{\mu\nu}\right)\sigma_{\mu\nu}\,,$
$\displaystyle\Sigma_{\alpha}$
$\displaystyle=N_{\alpha}{}^{\mu}e^{\nu}\sigma_{\mu\nu}\,,$
$\displaystyle\Sigma$ $\displaystyle=e^{\mu}e^{\nu}\sigma_{\mu\nu}\,,$ (46)

and

$\displaystyle\Omega^{\alpha}$
$\displaystyle=\frac{1}{2}N_{\gamma}{}^{\alpha}\varepsilon^{\mu\nu\gamma}\nabla_{\mu}u_{\nu}\,,$
$\displaystyle\Omega$
$\displaystyle=\frac{1}{2}\varepsilon^{\mu\nu}\nabla_{\mu}u_{\nu}\,.$ (47)

From their definitions, it is immediate to conclude that the covariantly
defined vector and 2-tensor quantities characterize the behavior of the
kinematical quantities on the surfaces $W$, whereas the scalars characterize
the behavior along $u$ or $e$.

Similarly, we can decompose the covariant derivative of $e$ along $u$, $e$ and
onto $W$ as

$\nabla_{\alpha}e_{\beta}=\frac{1}{2}N_{\alpha\beta}\phi+\zeta_{\alpha\beta}+\varepsilon_{\alpha\beta}\xi+e_{\alpha}a_{\beta}-u_{\alpha}\alpha_{\beta}-\mathcal{A}u_{\alpha}u_{\beta}+\left(\frac{1}{3}\theta+\Sigma\right)e_{\alpha}u_{\beta}+\left(\Sigma_{\alpha}-\varepsilon_{\alpha\mu}\Omega^{\mu}\right)u_{\beta}\,,$
(48)

where

$\displaystyle\phi$ $\displaystyle=N^{\mu\nu}\nabla_{\mu}e_{\nu}\,,$
$\displaystyle\zeta_{\alpha\beta}$
$\displaystyle=\left(\frac{N_{\alpha}{}^{\mu}N_{\beta}{}^{\nu}+N_{\alpha}{}^{\nu}N_{\beta}{}^{\mu}}{2}-\frac{1}{2}N_{\alpha\beta}N^{\mu\nu}\right)\nabla_{\mu}e_{\nu}\,,$
$\displaystyle\xi$
$\displaystyle=\frac{1}{2}\varepsilon^{\mu\nu}\nabla_{\mu}e_{\nu}\,,$ (49)

represent, respectively, the expansion scalar, the shear tensor, and the twist
of the congruence of the integral curves of $e$ projected on $W$, and

$\displaystyle a_{\alpha}$
$\displaystyle=e^{\mu}h_{\alpha}{}^{\nu}\nabla_{\mu}e_{\nu}\,,$
$\displaystyle\alpha_{\alpha}$
$\displaystyle=u^{\mu}h_{\alpha}{}^{\nu}\nabla_{\mu}e_{\nu}\,.$ (50)

###  A.3 Weyl and stress-energy tensors

For the Levi-Civita connection, the Riemann tensor can be defined by the Ricci
identity for an arbitrary 1-form $\chi$:

$R_{\alpha\beta\delta}{}^{\rho}\chi{}_{\rho}=\left(\nabla_{\alpha}\nabla_{\beta}-\nabla_{\beta}\nabla_{\alpha}\right)\chi_{\delta}\,.$
(51)

In the case of a manifold of dimension four, the components of the Riemann
curvature tensor, $R_{\alpha\beta\gamma\delta}$, are given by

$R_{\alpha\beta\gamma\delta}=C_{\alpha\beta\gamma\delta}+R_{\alpha\left[\gamma\right.}g_{\left.\delta\right]\beta}-R_{\beta\left[\gamma\right.}g_{\left.\delta\right]\alpha}-\frac{1}{3}R\,g_{\alpha\left[\gamma\right.}g_{\left.\delta\right]\beta}\,,$
(52)

where $C_{\alpha\beta\gamma\delta}$ represents the Weyl tensor. The Weyl
tensor plays a crucial role in relativistic gravity, describing the tidal
forces, frame-dragging, and the properties of gravitational waves. Remarkably,
the Weyl 4-tensor can be completely characterized by two 2-tensors

$\displaystyle E_{\alpha\beta}$
$\displaystyle=C_{\alpha\mu\beta\nu}u^{\mu}u^{\nu}\,,$ (53) $\displaystyle
H_{\alpha\beta}$
$\displaystyle=\frac{1}{2}\varepsilon_{\alpha}{}^{\mu\nu}C_{\mu\nu\beta\delta}u^{\delta}\,,$
(54)

respectively referred as the “electric” and “magnetic” part of the Weyl
tensor, both symmetric and traceless tensors, such that

$C_{\alpha\beta\gamma\delta}=-\varepsilon_{\alpha\beta\mu}\varepsilon_{\gamma\delta\nu}E^{\nu\mu}-2u_{\alpha}E_{\beta\left[\gamma\right.}u_{\left.\delta\right]}+2u_{\beta}E_{\alpha\left[\gamma\right.}u_{\left.\delta\right]}-2\varepsilon_{\alpha\beta\mu}H^{\mu}{}_{\left[\gamma\right.}u_{\left.\delta\right]}-2\varepsilon_{\mu\gamma\delta}H^{\mu}{}_{\left[\alpha\right.}u_{\left.\beta\right]}\,.$
(55)

In the 1+1+2 covariant formalism, the components of the Weyl tensor are
themselves decomposed as

$\displaystyle E_{\alpha\beta}$
$\displaystyle=\mathcal{E}\left(e_{\alpha}e_{\beta}-\frac{1}{2}N_{\alpha\beta}\right)+\mathcal{E}_{\alpha}e_{\beta}+e_{\alpha}\mathcal{E}_{\beta}+\mathcal{E}_{\alpha\beta}\,,$
(56) $\displaystyle H_{\alpha\beta}$
$\displaystyle=\mathcal{H}\left(e_{\alpha}e_{\beta}-\frac{1}{2}N_{\alpha\beta}\right)+\mathcal{H}_{\alpha}e_{\beta}+e_{\alpha}\mathcal{H}_{\beta}+\mathcal{H}_{\alpha\beta}\,,$

where

$\displaystyle\mathcal{E}$
$\displaystyle=E_{\mu\nu}e^{\mu}e^{\nu}=-N^{\mu\nu}E_{\mu\nu}\,,$
$\displaystyle\mathcal{H}$
$\displaystyle=e^{\mu}e^{\nu}H_{\mu\nu}=-N^{\mu\nu}H_{\mu\nu}\,,$ (57)
$\displaystyle\mathcal{E}_{\alpha}$
$\displaystyle=N_{\alpha}{}^{\mu}e^{\nu}E_{\mu\nu}=e^{\mu}N_{\alpha}{}^{\nu}E_{\mu\nu}\,,$
$\displaystyle\mathcal{H}_{\alpha}$
$\displaystyle=N_{\alpha}{}^{\mu}e^{\nu}H_{\mu\nu}=e^{\mu}N_{\alpha}{}^{\nu}H_{\mu\nu}\,,$
$\displaystyle\mathcal{E}_{\alpha\beta}$
$\displaystyle=E_{\left\\{\alpha\beta\right\\}}\,,$
$\displaystyle\mathcal{H}_{\alpha\beta}$
$\displaystyle=H_{\left\\{\alpha\beta\right\\}}\,.$

Lastly, to write the Einstein field equations in the language of the 1+1+2
formalism, we need the covariant decomposition of the metric stress-energy
tensor $T_{\alpha\beta}$. Then,

$T_{\alpha\beta}=\mu\,u_{\alpha}u_{\beta}+\left(p+\Pi\right)e_{\alpha}e_{\beta}+\left(p-\frac{1}{2}\Pi\right)N_{\alpha\beta}+2Qe_{\left(\alpha\right.}u_{\left.\beta\right)}+2Q_{\left(\alpha\right.}u_{\left.\beta\right)}+2\Pi_{\left(\alpha\right.}e_{\left.\beta\right)}+\Pi_{\alpha\beta}\,,$
(58)

with

$\displaystyle\mu$ $\displaystyle=u^{\mu}u^{\nu}T_{\mu\nu}\,,$ $\displaystyle
Q_{\alpha}$ $\displaystyle=-N_{\alpha}{}^{\mu}u^{\nu}T_{\mu\nu}\,,$ (59)
$\displaystyle p$
$\displaystyle=\frac{1}{3}\left(e^{\mu}e^{\nu}+N^{\mu\nu}\right)T_{\mu\nu}\,,$
$\displaystyle\Pi_{\alpha}$
$\displaystyle=N_{\alpha}{}^{\mu}e^{\nu}T_{\mu\nu}\,,$ $\displaystyle\Pi$
$\displaystyle=\frac{1}{3}\left(2e^{\mu}e^{\nu}-N^{\mu\nu}\right)T_{\mu\nu}\,,$
$\displaystyle\Pi_{\alpha\beta}$
$\displaystyle=\left(\frac{N_{\alpha}{}^{\mu}N_{\beta}{}^{\nu}+N_{\alpha}{}^{\nu}N_{\beta}{}^{\mu}}{2}-\frac{1}{2}N_{\alpha\beta}N^{\mu\nu}\right)T_{\mu\nu}\,.$
$\displaystyle Q$ $\displaystyle=-e^{\mu}u^{\nu}T_{\mu\nu}\,,$

The various contributions in the covariant decomposition of the stress-energy
tensor in Eq. (59) have direct physical meaning. Given an observer with
4-velocity $u$, $\mu$ represents the mass-energy density of the fluid, $p$ the
isotropic pressure, $Q$ and $Q_{\alpha}$ represent, respectively, heat and
momentum flows along $e$ and in $W$, and $\Pi$, $\Pi_{\alpha}$ and
$\Pi_{\alpha\beta}$ characterize the anisotropic pressure within the fluid.

##  Appendix B 1+1+2 scalar perturbation equations

Using the 1+1+2 formalism, briefly introduced in Appendix A, we list here the
covariant version of the equations that describe the scalar adiabatic
perturbations of a fluid distribution in a static, LRS II background. These
equations are written in terms of the variables in Eq. (5) and characterize
the perturbation from the point of view of an observer locally comoving with
the matter composing the relativistic compact stellar object. They read:

$\displaystyle\widehat{\mathsf{A}}-\ddot{\theta}$
$\displaystyle=\frac{1}{2}\left(\mathsf{m}+3\mathsf{p}\right)+\widehat{\mathcal{A}}_{0}\left(\frac{1}{3}\theta+\Sigma\right)-\left(3\mathcal{A}_{0}+\phi_{0}\right)\mathsf{A}-\mathcal{A}_{0}\mathsf{F}\,,$
(60) $\displaystyle\widehat{\mathsf{p}}$
$\displaystyle=\left(\frac{1}{3}\theta+\Sigma\right)\widehat{p}_{0}-\mathcal{A}_{0}\left(\mathsf{m}+2\mathsf{p}\right)-\left(\mu_{0}+p_{0}\right)\mathsf{A}\,,$
(61)
$\ddot{\Sigma}-\frac{2}{3}\ddot{\theta}=\frac{1}{3}\left(\mathsf{m}+3\mathsf{p}\right)-\mathsf{E}-\mathcal{A}_{0}\mathsf{F}-\phi_{0}\mathsf{A}\,,$
(62)
$\frac{2}{3}\widehat{\theta}-\widehat{\Sigma}=\frac{3}{2}\phi_{0}\Sigma\,.$
(63) $\displaystyle\mathsf{E}$
$\displaystyle=\mathcal{E}_{0}\left(\frac{3}{2}\Sigma-\theta\right)-\frac{1}{2}\left(\mu_{0}+p_{0}\right)\Sigma\,,$
(64)
$\mathsf{F}=\left(2\mathcal{A}_{0}-\phi_{0}\right)\left(\frac{1}{3}\theta-\frac{1}{2}\Sigma\right)\,,$
(65) $\mathsf{m}=-\left(\mu_{0}+p_{0}\right)\theta\,,$ (66)

Upon reorganization, using the background gravitational field equations and
applying a harmonic decomposition of the perturbation variables, the above
equations can be recast into the system in Sec. II.3.

The perturbation equations become increasingly more complicated as we relax
the constraints that we have assumed before. The general set of perturbation
equations for the variables in Eq. (5), together with additional perturbation
variables useful in other specific frames can be found in [9].

##  Appendix C  Derivation of Chandrasekhar’s radial pulsation equation

In this Appendix, we show how the system of equations for the comoving
observer in Appendix B returns the well-known Chandrasekhar radial pulsation
equation [1, 2]. We will then start by recalling the classical derivation of
the Chandrasekhar equation. Then, we will break covariance and gauge
invariance in the 1+1+2 perturbation equations to prove that they lead to the
same result.

Let us start by constructing the Chandrasekhar pulsation equation considering
perturbations of the metric directly. Consider the Schwarzschild coordinate
system $\left(t,r,\psi,\varphi\right)$ defined by an observer at spatial
infinity, such that events in both the equilibrium and the perturbed spacetime
can be identified unambiguously. Let the equilibrium spacetime be
characterized by a line element of the form

$ds^{2}=(g_{0})_{\alpha\beta}dx^{\alpha}dx^{\beta}=-e^{2\Phi_{0}(r)}dt^{2}+e^{2\Lambda_{0}(r)}dr^{2}+r^{2}d\Omega^{2}\,,$
(67)

where $d\Omega^{2}=d\psi^{2}+\sin\psi d\varphi^{2}$ represents the natural
line element of the unit 2-sphere. We will adopt the nomenclature of the body
of the text and indicate quantities in the background spacetime by a subscript
“0”.

Assuming the equilibrium spacetime to be permeated by a perfect fluid, such
that, from the point of view of an observer at rest with matter, the energy-
momentum tensor will be written as in Eq. (3)

$T_{\alpha\beta}=\left(\mu_{0}+p_{0}\right)(u_{0})_{\alpha}(u_{0})_{\beta}+p_{0}\left(g_{0}\right)_{\alpha\beta}\,,$
(68)

where

$\left(u_{0}\right)_{\alpha}dx^{\alpha}=-e^{\Phi_{0}}dt\,,$ (69)

the Einstein field equations are

$\displaystyle\Lambda^{\prime}_{0}$
$\displaystyle=\frac{1}{2}e^{2\Lambda_{0}}\mu_{0}r-\frac{e^{2\Lambda_{0}}}{2r}+\frac{1}{2r}\,,$
(70) $\displaystyle\Phi^{\prime}_{0}$
$\displaystyle=\frac{1}{2}e^{2\Lambda_{0}}p_{0}r+\frac{e^{2\Lambda_{0}}}{2r}-\frac{1}{2r}\,,$
(71) $\displaystyle\Phi^{\prime\prime}_{0}$
$\displaystyle=\frac{1}{4}e^{2\Lambda_{0}}\mu_{0}+\frac{1}{4}e^{4\Lambda_{0}}\mu_{0}+\frac{5}{4}e^{2\Lambda_{0}}p_{0}-\frac{3}{4}e^{4\Lambda_{0}}p_{0}+\frac{1}{4}e^{4\Lambda_{0}}\mu_{0}p_{0}r^{2}-\frac{1}{4}e^{4\Lambda_{0}}p_{0}^{2}r^{2}-\frac{e^{4\Lambda_{0}}}{2r^{2}}+\frac{1}{2r^{2}}\,,$
(72)

and the Bianchi identity

$p^{\prime}_{0}=-\Phi_{0}\left(\mu_{0}+p_{0}\right)\,.$ (73)

Now, radially perturbing the equilibrium configuration leads to a new metric
associated with the line element

$ds^{2}=g_{\alpha\beta}dx^{\alpha}dx^{\beta}=-e^{2\Phi(t,r)}dt^{2}+e^{2\Lambda(t,r)}dr^{2}+r^{2}d\Omega^{2}\,,$
(74)

and the energy-momentum tensor for an observer comoving with matter will be
described by

$T_{\alpha\beta}=\left(\mu+p\right)u_{\alpha}u_{\beta}+pg_{\alpha\beta}\,,$
(75)

where $u_{\alpha}$ represent, in a local coordinate system, the components of
the 4-velocity form of an observer comoving with the perturbed fluid. The
perturbation will induce a change in a given scalar variable $X$, which can be
characterized by the quantity

$\delta X=X\left(t,r\right)-X_{0}\left(r\right)\,.$ (76)

This expression constitutes the Eulerian representation of perturbations,
i.e., the description of the perturbation from the point of view of an
observer in a position with constant $r$ coordinate.

In addition to the $\delta$-variations, we need to introduce an extra
parameter to describe the radial displacement of a fluid element. A fluid
element at the coordinate $r$ in the unperturbed spacetime, is moved to a new
position $\tilde{r}$, in the perturbed spacetime, such that,

$\tilde{r}=r+\eta\left(t,r\right)\,.$ (77)

The parameter $\eta$ is usually called “radial displacement”. We need to be
especially careful in understanding the meaning of $\eta$. Even if this
function is connected to the choice of a local coordinate system, $\eta$
effectively also characterizes the mapping between the unperturbed and
perturbed spacetimes, as it is relabeling events. Therefore, it cannot be just
considered a product of a coordinate transformation, but rather a gauge. For
this reason, $\eta$ is often called “gauge parameter”. In this appendix, for
example, in our construction, we have tacitly considered the Schwarzschild
coordinate system to write the components of the perturbed metric and those of
the background metric. Therefore, we will find that the perturbations depend
directly on the gauge parameter. This, of course, is not a necessity, but it
is useful in the treatment of this type of perturbations. We stress, however,
that this approach should be used with caution. In the case of adiabatic,
radial perturbations, it is indeed possible to define the same coordinate
system in both the equilibrium and the perturbed spacetimes. Therefore, the
conclusions are not ambiguous. However, in general this is not the case, and
attributing physical observable effects to active coordinate transformations,
might lead to gauge-dependent, physically meaningless conclusions. In
addition, choosing from the start the coordinate system of the equilibrium and
the perturbed spacetimes will force the analysis to a specific gauge, which
might lead to spurious complexity of the equations.

In terms of Eulerian perturbations, at first perturbative order, that is,
disregarding quadratic and higher-order terms of the $\delta$-variations, we
have

$u_{\alpha}dx^{\alpha}=-e^{\Phi_{0}}\left(1+\delta\Phi\right)dt+e^{2\Lambda_{0}-\Phi_{0}}\eta\eta^{\ast}dr\,.\\\
$ (78)

where we have indicated the derivative with respect to $t$ with an “asterisk”.
Then the perturbed gravitational field equations associated with the metric
Eq. (74) read,

$\displaystyle\left(\delta\Lambda\right){}^{\ast}=$
$\displaystyle-\left(\Lambda^{\prime}_{0}+\Phi^{\prime}_{0}\right)\eta^{\ast}\,,$
(79) $\displaystyle\left(\delta\Lambda\right)^{\prime}=$
$\displaystyle\left(2r\Lambda^{\prime}_{0}-\frac{1}{r}\right)\delta\Lambda+\frac{1}{2}re^{2\Lambda_{0}}\delta\mu\,,$
(80) $\displaystyle\left(\delta\Phi\right)^{\prime}=$
$\displaystyle\left(2\Phi^{\prime}_{0}+\frac{1}{r}\right)\delta\Lambda+\frac{1}{2}re^{2\Lambda_{0}}\delta
p\,,$ (81) $\displaystyle\left(\delta\Phi\right)^{\prime\prime}=$
$\displaystyle
e^{2(\Lambda_{0}-\Phi_{0})}\delta\Lambda^{\ast\ast}+\frac{1}{2}e^{2\Lambda_{0}}\left(1+r\Phi^{\prime}_{0}\right)\delta\mu+\frac{1}{2}e^{4\Lambda_{0}}\left(1+r\Lambda^{\prime}_{0}-2r\Phi^{\prime}_{0}\right)\delta
p$
$\displaystyle-\frac{e^{2\Lambda_{0}}}{r^{2}}\left[2e^{2\Lambda_{0}}-r\Lambda^{\prime}_{0}\left(3+4r\Phi^{\prime}_{0}\right)+r\Phi^{\prime}_{0}\left(1+4r\Phi^{\prime}_{0}\right)\right]\delta\Lambda\,,$
(82)

and the Bianchi identities

$\displaystyle\delta\mu^{\ast}=\frac{1}{r}\left(-2\eta^{\ast}+r\eta^{\ast}\Phi^{\prime}_{0}-r\eta^{\prime\ast}\right)\left(\mu_{0}+p_{0}\right)-\eta^{\ast}\mu^{\prime}_{0}\,,$
(83) $\displaystyle\left(\delta
p\right)^{\prime}+\eta^{\ast\ast}e^{2(\Lambda_{0}-\Phi_{0})}\left(\mu_{0}+p_{0}\right)+\left(\delta\Phi\right)^{\prime}\left(\mu_{0}+p_{0}\right)+\left(\delta\mu+\delta
p\right)\Phi^{\prime}_{0}=0\,,$ (84)

where we have indicated the derivative with respect to $r$ with a “prime”.
Integrating Eq. (83) with respect to $t$ leads to

$\delta\mu=\frac{1}{r}\left(-2\eta+r\eta\Phi^{\prime}_{0}-r\eta^{\prime}\right)\left(\mu_{0}+p_{0}\right)-\eta\mu^{\prime}_{0}\,,$
(85)

where we have set the integration constant to zero as we assume $\delta\mu=0$
when $\eta=0$. As we have mentioned before, this choice of integration
constant follows from encoding the perturbations in the gauge parameter.
Consistently, in the above expressions, if $\eta=0$, we recover the background
spacetime.

Continuing, the system above is closed by providing an equation of state that
relates the pressure to the energy density. Following Chandrasekhar’s original
derivation, we will consider a barotropic equation of state, such that
$p=f\left(\mu\right)$. Using the equation of state, Eq. (85) allows us to
write $\delta p$ in terms of $\eta$ and its derivatives:

$\delta
p=\frac{\Gamma_{1}p_{0}}{r}\left(-2\eta+r\eta\Phi^{\prime}_{0}-r\eta^{\prime}\right)-\eta
p^{\prime}_{0}\,,$ (86)

where the quantity

$\Gamma_{1}=\frac{\mu_{0}+p_{0}}{p_{0}}f^{\prime}\left(\mu_{0}\right)\,,$ (87)

is the adiabatic index of the fluid. Integrating Eq. (79) with respect to $t$,
choosing the integration constant such that $\delta\Lambda=0$ when $\eta=0$,
yields

$\delta\Lambda=-\left(\Lambda^{\prime}_{0}+\Phi^{\prime}_{0}\right)\eta\,.$
(88)

Substituting the above relations and Eq. (81) in the momentum conservation
equation (84), after a laborious yet straightforward simplification, we obtain
an equation for $\eta$, the so-called Chandrasekhar radial pulsation equation.

Alternatively to the Eulerian description, we can define a Lagrangian
description of the perturbations. In this picture, an observer located at some
point with radial coordinate $r$ in the background spacetime is moved together
with matter to $\tilde{r}$ in the perturbed spacetime. We then define a
Lagrangian perturbation of a scalar quantity $X$ as

$\Delta X=X\left(t,\tilde{r}\right)-X_{0}\left(r\right)\,.$ (89)

The 1+1+2 perturbation variables we have used in this paper, correspond, by
definition, to Lagrangian perturbation variables. Indeed, since all the
quantities in Sec. II.3 are defined from the point of view of a comoving
observer, they are naturally Lagrangian.

Thus, in order to prove that the equations following from the covariant
approach imply the Chandrasekhar equation, we have to break covariance and
gauge invariance and show that the perturbation equations (60)–(63) reduce to
the perturbed field equations or to the Bianchi identities (79)–(84). In that
regard, we have to relate the Lagrangian variables
$\left\\{\mathcal{A},\phi,\theta,\Sigma\right\\}$ and
$\left\\{\mathsf{m},\mathsf{p},\mathsf{A},\mathsf{F},\mathsf{E},\theta,\Sigma\right\\}$
in terms of Eulerian variables to linear perturbation order.

We already know how to write, at first order, the 1-form associated with the
vector field $u$ in terms of the metric coefficients of the equilibrium and
the perturbed spacetime, Eq. (78). For the $e$ congruence, we have

$\displaystyle\left(e_{0}\right)_{\alpha}dx^{\alpha}$
$\displaystyle=e^{\Lambda_{0}}dr\,,$ (90) $\displaystyle
e_{\alpha}dx^{\alpha}$
$\displaystyle=-e^{\Lambda_{0}}\eta^{\ast}dt+e^{\Lambda_{0}}\left(1+\delta\Lambda\right)dr\,,$

Then, from the definitions in Appendix A and Eqs. (69), (78), (90) we obtain

$\displaystyle\mathcal{A}_{0}=$ $\displaystyle
e^{-\Lambda_{0}}\Phi^{\prime}_{0}\,,$ (91) $\displaystyle\phi_{0}=$
$\displaystyle\frac{2e^{-\Lambda_{0}}}{r}\,,$ (92) $\displaystyle\theta_{0}=$
$\displaystyle 0\,,$ (93) $\displaystyle\Sigma_{0}=$ $\displaystyle 0\,,$ (94)

and

$\displaystyle\mathcal{A}=$ $\displaystyle
e^{-\Lambda_{0}}\Phi^{\prime}_{0}+e^{-\Lambda_{0}}\delta\Phi^{\prime}-e^{-\Lambda_{0}}\Phi^{\prime}_{0}\delta\Lambda+\eta^{\ast\ast}e^{\Lambda_{0}-2\Phi_{0}}\,,$
(95) $\displaystyle\phi=$
$\displaystyle\frac{2e^{-\Lambda_{0}}}{r}\left(1-\delta\Lambda\right)\,,$ (96)
$\displaystyle\theta=$
$\displaystyle\frac{e^{-\Phi_{0}}}{r}\left(r\,\delta\Lambda^{*}+2\eta^{\ast}+r\eta^{\ast}\Lambda^{\prime}_{0}+r\eta^{\prime\ast}\right)\,,$
(97) $\displaystyle\Sigma=$
$\displaystyle\frac{2}{3}\frac{e^{-\Phi_{0}}}{r}\left(r\,\delta\Lambda^{*}+\eta^{\ast}+r\eta^{\ast}\Lambda^{\prime}_{0}+r\eta^{\prime\ast}\right)\,.$
(98)

Now, the following relation holds between the Lagrangian and Eulerian
perturbations of a scalar quantity $X$:

$\Delta X=\delta X-X^{\prime}_{0}\eta\,.$ (99)

This equation can be used to relate the Lagrangian variables
$\left\\{\mathsf{m},\mathsf{p},\mathsf{A}\right\\}$ with their Eulerian
counterparts. For the Lagrangian perturbation of the time derivative of the
pressure, $\mathsf{p}$, we find

$\mathsf{p}=\dot{p}-\mathring{p}_{0}\,,$ (100)

where we have defined

$\mathring{p}=\left(u_{0}\right)^{\alpha}\partial_{\alpha}p\,.$ (101)

Expressing the right-hand side of Eq. (100) in terms of Eulerian perturbations
yields

$\mathsf{p}=\left(\delta
p\right)^{\boldsymbol{\cdot}}+p^{\prime}_{0}\dot{\eta}=e^{-\Phi_{0}}\left[\left(\delta
p\right)^{\ast}+p^{\prime}_{0}\eta^{\ast}\right]\,.$ (102)

It is also useful to find the expression for the hat derivative of
$\mathsf{p}$ in terms of the metric perturbations:

$\widehat{\mathsf{p}}=e^{-\Lambda_{0}}\Phi^{\prime}_{0}\mathsf{p}-e^{-\left(\Phi_{0}+\Lambda_{0}\right)}\left[\left(\delta
p\right)^{\prime\ast}+p^{\prime\prime}_{0}\eta^{\ast}+p^{\prime}_{0}\eta^{\prime\ast}\right]\,.$
(103)

Similarly, we obtain the following expressions for the $\mathsf{m}$ and
$\mathsf{A}$ quantities:

$\displaystyle\mathsf{m}$
$\displaystyle=e^{-\Phi_{0}}\left[\left(\delta\mu\right)^{\ast}+\mu^{\prime}_{0}\eta^{\ast}\right]\,,$
(104) $\displaystyle\mathsf{A}$
$\displaystyle=-e^{-\Lambda_{0}-3\Phi_{0}}\left[e^{2\Phi_{0}}\Phi^{\prime}_{0}\delta\Lambda^{\ast}-e^{2\Phi_{0}}\left(\delta\Phi\right)^{\prime\ast}+e^{2\Phi_{0}}\left(\Lambda^{\prime}_{0}\Phi^{\prime}_{0}-\Phi^{\prime\prime}_{0}\right)\eta^{\ast}-e^{2\Lambda_{0}}\eta^{\ast\ast\ast}\right]\,.$
(105)

Finally, using the above results, Eqs. (4), (64) and (65), we can obtain
algebraically the expressions for $\mathcal{E}_{0}$, $\mathsf{F}$ and
$\mathsf{E}$ in terms of the coefficients of the metric in Eq. (74) and their
perturbations.

Gathering all the above results and substituting them in the perturbed
momentum conservation equation for the comoving frame, Eq. (61), yields, after
some fairly long calculations, the momentum conservation Eq. (84). Applying
the same substitutions in the remaining equations of Appendix B yield a system
of equations equivalent to combinations of the perturbed Einstein equations
(79)–(88). Thus, we conclude that in the comoving frame, the covariant, gauge-
invariant equations for adiabatic, radial perturbations are equivalent to the
ones derived from the metric, coordinate based approach.

##  Appendix D Matrix $\mathds{A}$

Consider the matrix $\Theta$ defined in Eq. (25). To shorten the expressions,
let $\Theta_{ij}$ represent the $ij$-entry of $\Theta$, and
$\left(\Theta_{n}\right)_{ij}$ be interpreted as the $ij$-entry of the $n$th
order coefficient of the power expansion of $\Theta$ at $r=0$, Eq. (26). Then,
the matrix $\mathds{A}$ in Eq. (28) is given by

$\mathds{A}=\left[\begin{array}[]{lcr}A_{11}&A_{12}&A_{13}\\\
A_{21}&A_{22}&A_{23}\\\ A_{31}&A_{32}&A_{33}\end{array}\right]\,,$

where

$A_{11}=\Theta_{11}-\frac{1}{3}r^{2}\Theta_{31}\left[\left(\Theta_{0}\right)_{12}\left(\Theta_{0}\right)_{23}-3\left(\Theta_{1}\right)_{13}\right]\,,$
(106) $\displaystyle A_{12}=$
$\displaystyle-\frac{36\Theta_{13}}{r^{3}}+\frac{12}{r^{2}}\left[\Theta_{12}\left(\Theta_{0}\right)_{23}-\left(\Theta_{0}\right)_{12}\left(\Theta_{0}\right)_{23}+3\left(\Theta_{1}\right)_{13}\right]+\frac{12}{r}\left(\Theta_{33}-\Theta_{11}\right)\left[\left(\Theta_{0}\right)_{12}\left(\Theta_{0}\right)_{23}-3\left(\Theta_{1}\right)_{13}\right]$
(107)
$\displaystyle+\frac{2}{3}r^{2}\Theta_{32}\left[\left(\Theta_{0}\right)_{12}\left(\Theta_{0}\right)_{23}-3\left(\Theta_{1}\right)_{13}\right]\left\\{6\left(\Theta_{0}\right)_{23}\left[\left(\Theta_{1}\right)_{22}-\left(\Theta_{1}\right)_{33}\right]-9\left(\Theta_{2}\right)_{23}+2\left(\Theta_{0}\right)_{32}\left[\left(\Theta_{0}\right)_{23}\right]^{2}\right\\}$
$\displaystyle-4\Theta_{32}\left(\Theta_{0}\right)_{23}\left[\left(\Theta_{0}\right)_{12}\left(\Theta_{0}\right)_{23}-3\left(\Theta_{1}\right)_{13}\right]-12\Theta_{12}\left(\Theta_{0}\right)_{23}\left[\left(\Theta_{1}\right)_{22}-\left(\Theta_{1}\right)_{33}\right]+18\Theta_{12}\left(\Theta_{2}\right)_{23}$
$\displaystyle-4\Theta_{12}\left(\Theta_{0}\right)_{32}\left[\left(\Theta_{0}\right)_{23}\right]^{2}+4r\Theta_{31}\left[\left(\Theta_{0}\right)_{12}\left(\Theta_{0}\right)_{23}-3\left(\Theta_{1}\right)_{13}\right]^{2}\,,$

$A_{13}=\frac{1}{3}r^{3}\Theta_{32}\left[\left(\Theta_{0}\right)_{12}\left(\Theta_{0}\right)_{23}-3\left(\Theta_{1}\right)_{13}\right]-r\Theta_{12}\,,$
(108)

$A_{21}=-\frac{1}{36}r^{3}\Theta_{31}\,,$ (109) $\displaystyle A_{22}$
$\displaystyle=\frac{1}{9}r^{3}\Theta_{32}\left\\{3\left(\Theta_{0}\right)_{23}\left[\left(\Theta_{1}\right)_{22}-\left(\Theta_{1}\right)_{33}\right]-\frac{9}{2}\left(\Theta_{2}\right)_{23}+\left[\left(\Theta_{0}\right)_{23}\right]^{2}\left(\Theta_{0}\right)_{32}\right\\}$
(110)
$\displaystyle+\frac{1}{3}r^{2}\Theta_{31}\left[\left(\Theta_{0}\right)_{12}\left(\Theta_{0}\right)_{23}-3\left(\Theta_{1}\right)_{13}\right]-\frac{r}{3}\left(\Theta_{0}\right)_{23}\Theta_{32}+\Theta_{33}\,,$

$A_{23}=\frac{1}{36}r^{4}\Theta_{32}\,,$ (111)

$\displaystyle A_{31}$
$\displaystyle=\frac{1}{9}\Theta_{31}\left\\{3\left(\Theta_{0}\right)_{23}\left[r^{2}\left(\Theta_{1}\right)_{22}-r^{2}\left(\Theta_{1}\right)_{33}-1\right]-\frac{9}{2}r^{2}\left(\Theta_{2}\right)_{23}+r^{2}\left(\Theta_{0}\right)_{32}\left[\left(\Theta_{0}\right)_{23}\right]^{2}\right\\}-\frac{\Theta_{21}}{r}\,,$
(112) $\displaystyle A_{32}=$
$\displaystyle-\frac{2r}{3}\Theta_{31}\left[\left(\Theta_{0}\right)_{12}\left(\Theta_{0}\right)_{23}-3\left(\Theta_{1}\right)_{13}\right]\left\\{6\left(\Theta_{0}\right)_{23}\left[\left(\Theta_{1}\right)_{22}-\left(\Theta_{1}\right)_{33}\right]-9\left(\Theta_{2}\right)_{23}+2\left[\left(\Theta_{0}\right)_{23}\right]^{2}\left(\Theta_{0}\right)_{32}\right\\}$
(113)
$\displaystyle-r^{2}\Theta_{32}\left(\Theta_{0}\right)_{23}\left[\left(\Theta_{1}\right)_{22}-\left(\Theta_{1}\right)_{33}\right]\left\\{4\left(\Theta_{0}\right)_{23}\left[\left(\Theta_{1}\right)_{22}-\left(\Theta_{1}\right)_{33}\right]-6\left(\Theta_{2}\right)_{23}+\frac{4}{3}\left[\left(\Theta_{0}\right)_{23}\right]^{2}\left(\Theta_{0}\right)_{32}\right\\}$
$\displaystyle-r^{2}\Theta_{32}\left[\left(\Theta_{0}\right)_{23}\right]^{2}\left(\Theta_{0}\right)_{32}\left\\{\frac{4}{3}\left(\Theta_{0}\right)_{23}\left[\left(\Theta_{1}\right)_{22}-\left(\Theta_{1}\right)_{33}\right]-2\left(\Theta_{2}\right)_{23}+\frac{4}{9}\left[\left(\Theta_{0}\right)_{23}\right]^{2}\left(\Theta_{0}\right)_{32}\right\\}$
$\displaystyle+\frac{36}{r^{4}}\left[\Theta_{23}-\left(\Theta_{0}\right)_{23}\right]-\frac{12}{r^{3}}\left(\Theta_{22}-\Theta_{33}\right)\left(\Theta_{0}\right)_{23}+\frac{12}{r^{2}}\Theta_{21}\left[\left(\Theta_{0}\right)_{12}\left(\Theta_{0}\right)_{23}-3\left(\Theta_{1}\right)_{13}\right]$
$\displaystyle+\frac{2}{r}\left(\Theta_{22}-\Theta_{33}\right)\left\\{6\left(\Theta_{0}\right)_{23}\left[\left(\Theta_{1}\right)_{22}-\left(\Theta_{1}\right)_{33}\right]-9\left(\Theta_{2}\right)_{23}+2\left[\left(\Theta_{0}\right)_{23}\right]^{2}\left(\Theta_{0}\right)_{32}\right\\}$
$\displaystyle+\frac{4}{3}\left(\Theta_{0}\right)_{23}\Theta_{32}\left\\{6\left(\Theta_{0}\right)_{23}\left[\left(\Theta_{1}\right)_{22}-\left(\Theta_{1}\right)_{33}\right]-9\left(\Theta_{2}\right)_{23}+2\left[\left(\Theta_{0}\right)_{23}\right]^{2}\left(\Theta_{0}\right)_{32}\right\\}$
$\displaystyle+r^{2}\Theta_{32}\left(\Theta_{2}\right)_{23}\left\\{6\left(\Theta_{0}\right)_{23}\left[\left(\Theta_{1}\right)_{22}-\left(\Theta_{1}\right)_{33}\right]-9\left(\Theta_{2}\right)_{23}+2\left[\left(\Theta_{0}\right)_{23}\right]^{2}\left(\Theta_{0}\right)_{32}\right\\}$
$\displaystyle+\frac{2}{r^{2}}\left\\{6\left(\Theta_{0}\right)_{23}\left[\left(\Theta_{1}\right)_{22}-\left(\Theta_{1}\right)_{33}\right]-9\left(\Theta_{2}\right)_{23}-2\left[\left(\Theta_{0}\right)_{23}\right]^{2}\left[\Theta_{32}-\left(\Theta_{0}\right)_{32}\right]\right\\}$
$\displaystyle+\frac{4}{r}\Theta_{31}\left(\Theta_{0}\right)_{23}\left[\left(\Theta_{0}\right)_{12}\left(\Theta_{0}\right)_{23}-3\left(\Theta_{1}\right)_{13}\right]\,,$

$A_{33}=\frac{r}{9}\Theta_{32}\left\\{3\left(\Theta_{0}\right)_{23}\left(1-r^{2}\left[\left(\Theta_{1}\right)_{22}-\left(\Theta_{1}\right)_{33}\right]\right)+\frac{9}{2}r^{2}\left(\Theta_{2}\right)_{23}-\left[\left(\Theta_{0}\right)_{23}\right]^{2}\left(\Theta_{0}\right)_{32}r^{2}\right\\}+\Theta_{22}\,.$
(114)

At first glance, it might seem that the $\mathds{A}$ matrix is singular at
$r=0$. However, after direct substitution and simplification, one can verify
that $\mathds{A}$ is real and analytic, and it has the same radius of
convergence as the matrix $\Theta$.

## References

  * [1] S. Chandrasekhar, “Dynamical instability of gaseous masses approaching the Schwarzschild limit in general relativity”, Phys. Rev. Lett. 12, 114 (1964); Erratum: Phys. Rev. Lett. 12, 437 (1964).
  * [2] S. Chandrasekhar, “The dynamical instability of gaseous masses approaching the Schwarzschild limit in general relativity”, Astrophysical Journal 140, 417 (1964); Erratum: Astrophysical Journal 140, 1342 (1964).
  * [3] J. M. Bardeen, K. S. Thorne and D. W. Meltzer, “A Catalogue of Methods for Studying the Normal Modes of Radial Pulsation of General-Relativistic Stellar Models”, Astrophysical Journal 145, 505 (1966).
  * [4] G. Chanmugam, “Radial oscillations of zero-temperature white dwarfs and neutron stars below nuclear densities”, Astrophysical Journal 217, 799 (1977).
  * [5] H. M. Väth and G. Chanmugam, “Radial oscillations of neutron stars and strange stars”, Astron. Astrophys. 260, 250 (1992).
  * [6] D. Gondek, P. Haensel, and J. L. Zdunik, “Radial pulsations and stability of protoneutron stars”, Astron. Astrophys. 325, 217 (1997).
  * [7] D. Gondek and J. L. Zdunik, “Avoided crossings in radial pulsations of neutron and strange stars”, Astron. Astrophys. 344, 117 (1999).
  * [8] K. D. Kokkotas and J. Ruoff, “Radial oscillations of relativistic stars”, Astron. Astrophys. 366, 565 (2001).
  * [9] P. Luz and S. Carloni, “Gauge invariant perturbations of static spatially compact LRS II spacetimes”, arXiv:2405.05321 [gr-qc].
  * [10] C. A. Clarkson and R. K. Barrett, “Covariant perturbations of Schwarzschild black holes”, Class. Quantum Grav. 20, 3855 (2003).
  * [11] C. A. Clarkson, “A covariant approach for perturbations of rotationally symmetric spacetimes”, Phys. Rev. D 76, 104034 (2007).
  * [12] M. Marklund et al., “Charged multifluids in general relativity”, Class. Quantum Grav. 20, 1823 (2003).
  * [13] M. Bradley et al., “Density growth in Kantowski-Sachs cosmologies with a cosmological constant”, Class. Quantum Grav. 29, 095023 (2012).
  * [14] S. Carloni and D. Vernieri, “Covariant Tolman-Oppenheimer-Volkoff equations. I. The isotropic case”, Phys. Rev. D 97, 124056 (2018).
  * [15] S. Carloni and D. Vernieri, “Covariant Tolman-Oppenheimer-Volkoff equations. II. The anisotropic case”, Phys. Rev. D 97, 124057 (2018).
  * [16] P. Luz and S. Carloni, “Static compact objects in Einstein-Cartan theory”, Phys. Rev. D 100, 084037 (2019).
  * [17] R. Törnkvist and M. Bradley, “General perfect fluid perturbations of homogeneous and orthogonal locally rotationally symmetric class II cosmologies”, Phys. Rev. D 100, 124043 (2019).
  * [18] N. F. Naidu, S. Carloni and P. Dunsby, “Two-fluid stellar objects in general relativity: The covariant formulation”, Phys. Rev. D 104, 044014 (2021).
  * [19] N. F. Naidu, S. Carloni and P. Dunsby, “Anisotropic two-fluid stellar objects in general relativity”, Phys. Rev. D 106, 124023 (2022).
  * [20] J. L. Rosa and S. Carloni, “Junction conditions for general LRS spacetimes in the $1+1+2$ covariant formalism”, arXiv:2303.12457 [gr-qc].
  * [21] M. Campbell, et al., “Reconstructing exact solutions to relativistic stars in $f(R)$ gravity”, arXiv:2403.00070 [gr-qc].
  * [22] M. S. R. Delgaty and K. Lake, “Physical acceptability of isolated, static, spherically symmetric, perfect fluid solutions of Einstein’s equations”, Comput. Phys. Commun. 115, 395 (1998).
  * [23] J. M. Stewart and M. Walker, “Perturbations of Space-Times in General Relativity”, Proc. R. Soc. Lond. A 341, 49 (1974).
  * [24] D. W. Meltzer and K. S. Thorne, “Normal Modes of Radial Pulsation of Stars at the End Point of Thermonuclear Evolution”, Astrophysical Journal 145, 514 (1966).
  * [25] E. A. Coddington and N. Levinson, _Theory of ordinary differential equations_(McGraw-Hill, New York, 1955).
  * [26] A. Alho et al., “Compactness bounds in General Relativity”, Phys. Rev. D 106, L041502 (2022).
  * [27] A. Alho et al., “Spherically symmetric elastic bodies in general relativity”, Class. Quantum Grav. 41, 073002 (2024).
  * [28] A. D. D. Masa, J. P. S. Lemos and V. T. Zanchin, “Stability of electrically charged stars, regular black holes, quasiblack holes, and quasinonblack holes”, Phys. Rev. D 107, 064053 (2023).
  * [29] J. M. Lattimer and M. Prakash, “Neutron Star Observations: Prognosis for Equation of State Constraints,” Phys. Rept. 442, 109 (2007).
  * [30] A. Urbano and H. Veermäe, “On gravitational echoes from ultracompact exotic stars,” JCAP 04, 011 (2019).