# The representation and computational efficiency of the
Tolman–Oppenheimer–Volkoff equations in isotropic coordinates

Dániel Barta  HUN-REN Wigner RCP, H-1525 Budapest 114, P.O. Box 49, Hungary
<EMAIL_ADDRESS>
###### Abstract

This study aims to provide an analytical scheme for computing equilibrium
configurations of relativistic stars by solving the Tolman–Oppenheimer–Volkoff
equations directly in isotropic polar coordinates, as opposed to the commonly
applied methods of rescaling the radial profile of corresponding solutions
obtained in curvature coordinates. This study also provides evidence that the
differential equation for gravitational mass may be replaced by an algebraic
expression relating the metric potential to the energy density in the form of
the quartic equation. Nevertheless, the greater computational expense of
evaluating the algebraic equation renders its application less efficient. A
further objective of this study was to evaluate the performance of the present
computational scheme in the computational time and numerical accuracy. Our
results indicate that the computational time increases with the stiffness of
the constituent matter inside the star. Conversely, the absolute difference
between the gravitational mass obtained by the proposed method and that
computed via the use of LORENE packages initially increases rapidly with the
central energy density, but the rate of growth subsequently declines as the
maximum mass configuration is approached.

##  1 Introduction

Over a century has passed since the Einstein field equations (EFEs) for a
perfect-fluid body in a static, spherically symmetric spacetime were first
formulated in isotropic coordinates by Arthur Eddington in 1923. [1] This
formulation of the EFEs, which consists of a coupled system of second-order
ordinary differential equations (ODE) in the metric variables, also appears in
Richard Tolman’s two research papers [2, 3] in 1930, as well as Georges
Lemaître’s famous monograph [4], where he recognised in 1931 that the mass
function is interchangeable by the coefficient of the spatial metric function.
It has been realised by these early works [2, 3, 4] that two of the three non-
zero components of the EFEs are second-order ODEs in the metric variables and,
according to Ref. [7], can regarded merely as a rather elementary extension of
energy–momentum conservation. Consequently, albeit one of the ODEs is replaced
by the equation of local conservation of energy–momentum, there is still at
least one second-order equation left, making exact solutions either
analytically or numerically unattainable. [8] The issue of finding exact
solutions for these EFEs was first addressed by Tolman [5], and independently
by J. Robert Oppenheimer and George Volkoff [6] in 1939. The latter revealed
that both the use of curvature coordinates (also called Schwarzschild
coordinates) in the EFEs and introducing the function of gravitational mass
considerably simplify the ODE system [9] to a system of two first-order
differential and one algebraic equations, with the boundary condition of all
three variables.111Let us note that a recent correspondence [9, 7, 8] between
Semiz and Herrera argued that instead of “Tolman–Oppenheimer–Volkoff
equation,” the “Oppenheimer–Volkoff equation” be a more suitable name for the
solutions to the problem at hand, but this study does not wish to take side
over the appropriate nomenclature.

An extensive list of exact solutions to the EFEs had been collected in the
literature (see e.g. [10, 11] and references therein), where the solutions are
related to some specific restrictions imposed on the geometry of the spacetime
(such as spherical or axial symmetries) or a particular distribution of matter
(such as an isotropic distribution) within the spacetime. However, not all of
them can be considered as physically acceptable solutions, only those that
satisfy a certain set of conditions that serve to maintain the causality of
the spacetime structure and the energy conditions. The number of suitable
solutions can be further reduced by considering only those exact solutions
that are associated with isolated spheres of perfect fluid with boundary
conditions imposed at a finite radius. Even though some of these exact
solutions describe identical geometries in different coordinate frames, the
transformation from one frame to another, as argued by Ref. [12], is not
always well-defined.222Specifically, Ref. [13] emphasised that Schwarzschild-
like coordinates associated to an arbitrary line element can always be
transformed into isotropic polar coordinates, but the reverse process cannot
always be performed. The first exact solutions in isotropic coordinates was
discovered in 1950 by Nariai [13], where one is identical with Tolman’s
solution III333This solution is the Schwarzschild interior solution, which
also includes Tolman’s solution I (also known as Einstein’s static universe),
as a special case where the constant of integration $c$ is set equal to zero.,
another one is associated with new static and homogenous cosmological
solution, and most importantly, one represents an entirely new type of
“stellar” solution involving trigonometrical and logarithmic functions of the
radial coordinate. Further regular and well-behaving exact solutions in
isotropic coordinates were found by Refs. [14] for spherically symmetric
distributions of isotropic perfect fluid.

Despite the fact that the TOV equations in the literature are most commonly
expressed in terms of curvature coordinates, the quasi-isotropic polar
coordinates are a more commonly preferred choice to initiate dynamical
calculations that are implemented in a formalism based on the $3+1$
decomposition. The efficiency of numerical methods in these codes is greatly
improved by being performed in quasi-isotropic coordinates, where the three
spatial dimensions are treated similarly. One of the well-known examples is
the LORENE code [15] which uses a multi-domain spectral method to solve the
EFEs decomposed to a system of five quasi-linear elliptic equations. Ref. [16]
established that one may adapt maximal-slicing quasi-isotropic (MSQI)
coordinates in a stationary and circular axisymmetric spacetime for the
description of the 3+1 setting of EFEs. As a result, the solution to the EFEs
is directly obtained in quasi-isotropic coordinates but not through solving
the set of equations corresponding to the TOV equations. Another core
computational tool in relativistic astrophysics is the Einstein Toolkit [17],
which applies the modular Cactus framework, consisting of general modules
called “thorns.” The TOVSolver thorn solves the standard TOV equations
expressed in curvature coordinates, then transforms every variable into
isotropic coordinates by integrating the radius conversion formula
$\displaystyle\frac{\partial\log(r/\tilde{r})}{\partial\tilde{r}}=\frac{\tilde{r}^{1/2}-(\tilde{r}-2m)^{1/2}}{\tilde{r}(\tilde{r}-2m)^{1/2}}$,
which is subject to the boundary condition (18) in the exterior region of
spacetime. [18] In converting the solution into the variables required for a
dynamical evolution, one may assume that the metric is conformally flat, with
a conformal factor given by $\psi=\tilde{r}/r$, or equivalently, a logarithmic
conformal factor $\phi=\frac{1}{2}\log(\tilde{r}/r)$. The Einstein Toolkit
framework contains three routines that can read in publicly available data
generated by the LORENE code. Although LORENE was originally developed for
computing rapidly rotating neutron stars, it has been gradually extended to be
capable of computing various kinds of equilibrium and quasi-equilibrium
initial data sets of single or binary compact objects. Another code, COCAL was
developed by Uryū & Tsokaros [19] for the very purpose of computing such
initial data sets. The numerical method used in COCAL for computing
equilibrium configurations is based on an extension of the
Komatsu–Eriguchi–Hachisu (KEH) method [20, 21] which is less technical than
the spectral methods utilised in LORENE. Ref. [22] explains that instead of
using spectral methods in a 3+1 setting, the one-dimensional KEH solver seeks
solutions for the TOV equations in terms of curvature coordinates, then
rescales the radial profile of the solution to that of an isotropic setting by
$d\tilde{r}/dr=(\tilde{r}/r)\sqrt{1-2m(\tilde{r})/\tilde{r}}$. The Appendix B
of Ref. [22] briefly explains the transformation of the TOV equations from
curvature coordinates to isotropic coordinates, which we have reproduced using
our notations for the sake of comparison in the discussion of our results in
Sec. 6.

The present study aims, on the one hand, to present a partially uncoupled
system of first-order ODEs, derived in isotropic coordinates, which
corresponds to the TOV equations written in curvature coordinates. This
involves proving that the differential equation for gravitational mass may be
replaced by an algebraic expression relating the metric potential to the
energy density, even in the context of isotropic coordinates. On the other
hand, this study is intended to outline certain desirable attributes of
computing equilibrium configurations by directly solving the TOV equations in
isotropic coordinates rather than applying a method of rescaling the radial
profile of corresponding solutions obtained in curvature coordinates.

##  2 Centrally symmetric spacetime represented in different types of
coordinate frames

Let us consider a gravitational field that possesses central symmetry. This
aspect is expressed through the square of the line element, which is given by

$ds^{2}=A(r)dt^{2}+B(r)dr^{2}+2C(r)drdt+D(r)d\Omega^{2},$ (1)

where $A,\,B,\,C,\,D$ are all certain functions of the radial coordinate $r$
and

$d\Omega^{2}=\sin^{2}\theta d\phi^{2}+d\theta^{2}$ (2)

denotes the volume element for the 2-sphere of unit radius. One may subject
the coordinates to any transformation that does not destroy the central
symmetry of (1) over the spatial coordinates. This implies that a
transformation $r=r(\tilde{r})$ can be found where $r$ is any function of the
new coordinate $\tilde{r}$.444See the transformation laws (7) and (18) for the
interior perfect-fluid region and for the exterior vacuum region of the
Schwarzschild spacetime, respectively. One may choose the radial coordinate in
such a way that:

  1. 1.

the metric function $C(\tilde{r})$, being the coefficient of the term $2drdt$,
vanishes and the coefficient $D(\tilde{r})$ becomes equal to $\tilde{r}^{2}$;

  2. 2.

or the metric function $C(r)$ vanishes once again, but this time the metric
function $D(r)$ will be identical to $B(r)$.

Having the remaining metric-potential functions conveniently expressed in
exponential forms as

$A(\tilde{r})=e^{2\tilde{\nu}(\tilde{r})},\,B(\tilde{r})=e^{2\tilde{\lambda}(\tilde{r})}\quad\text{or}\quad
A(r)=e^{2\nu(r)},\,B(r)=e^{2\mu(r)},$ (3)

respectively, the metric of an arbitrary static spherically symmetric
spacetime can thus be expressed by spherical polar coordinates (called
Schwarzschild coordinates)
$(\tilde{t},\,\tilde{r},\,\tilde{\theta},\,\tilde{\phi})$ as

$ds^{2}=-e^{2\tilde{\nu}}d\tilde{t}^{2}+e^{2\tilde{\lambda}}d\tilde{r}^{2}+\tilde{r}^{2}(\sin^{2}\tilde{\theta}d\tilde{\phi}^{2}+d\tilde{\theta}^{2}),$
(4)

or equivalently, by isotropic polar coordinates $(t,\,r,\,\theta,\,\phi)$ as

$ds^{2}=-e^{2\nu}dt^{2}+e^{2\mu}[dr^{2}+r^{2}(\sin^{2}\theta
d\phi^{2}+d\theta^{2})].$ (5)

Setting side-by-side the line elements (4–5), and by requiring that the
angular and the radial parts to be equal, respectively, one obtains

$\tilde{r}^{2}\equiv e^{2\mu}r^{2},\quad e^{2\mu}dr^{2}\equiv
e^{2\tilde{\lambda}}d\tilde{r}^{2}.$ (6)

The earlier identity in eqs. (6), which stems from the equality of the angular
parts, defines a relation

$r=\exp\left(\int
e^{\tilde{\lambda}}\tilde{r}^{-1}d\tilde{r}\right)=\exp\left[-\operatorname{Ei}\left(-\tilde{\lambda}\tilde{r}\right)\right]$
(7)

between the radial coordinates of the two coordinate frames. In addition, the
earlier identity in eqs. (6) implies another relationship between the
derivatives:

$d\tilde{r}=e^{\mu}(r\mu^{\prime}+1)dr.$ (8)

The later identity in eqs. (6), with the substitution of $dr$ from (8), yields

$e^{2\tilde{\lambda}}=(r\mu^{\prime}+1)^{-2}.$ (9)

Figure 1: The isotropic radial coordinate $r$ as a function of the curvature
radial coordinate $\tilde{r}$ shown on logarithmic scale. In the interior
region of the spacetime, the relation (7) holds. This relation is displayed by
a set of coloured lines; each line corresponds to a configuration of
relativistic stars with different values of gravitational mass, and each line
is marked by the bullet point symbol $\bullet$ in Figure 4b, respectively. In
the exterior region, the relation is simply expressed by eq. (18) and
represented by the dashed black line. The inset figure on top shows the the
inner region $0.05\leq r/M\leq 4.50$ on linear scale.

Fig. 1 shows the relationship between the radial coordinates of the two
coordinate frames (7) for five equilibrium stellar configurations constructed
with the EoS SFHo/DD2, each with different values of gravitational mass. The
figure implies that $\tilde{r}$ is very close to being a linear function of
$r$ with an intercept of zero and a slope barely higher than one. The slope of
line is initially larger for larger gravitational mass and it decreases
slightly upon approaching the exterior region where it smoothly merges into
the parabola (18) (depicted by the dashed black line). The defining difference
between the two choices of coordinate frames lies in the interpretation of the
radial coordinate. The radial coordinate $\tilde{r}$ in (4) possesses a
natural geometric property in terms of the surface $\tilde{A}$ by being
conveniently defined as $\tilde{r}=\sqrt{\tilde{A}/4\pi}$. In contrast to
that, the radial coordinate $r$ in eq. (5) is defined so that the light cones
on constant time-slices (or hypersurfaces) $\Sigma_{t}$ with $t=t_{0}$ appear
as round spheres. The line element (5) implies that the metric restricted to
any of these hypersurfaces is

$g|_{t=t_{0},\,r=r_{0}}=e^{2\mu(r_{0})}r_{0}^{2}g_{\Omega},$ (10)

where $g_{\Omega}\equiv d\Omega^{2}$ is the Riemannian metric on the 2-sphere
of unit radius. This conveys that these nested coordinate-spheres (10) do in
fact represent geometric spheres, but the radial coordinate $r$ does not
faithfully represent radial distances within the nested coordinate-spheres. On
the other hand, angles in the constant time-slices (or hypersurfaces) are
represented without distortion by the angular coordinates $(\theta,\,\phi)$ in
isotropic charts. [23] This advantageous feature of (quasi-)isotropic
coordinates is exploited in the 3+1 formalism of general relativity. In the
3+1 formalism, the 4-dimensional spacetime is foliated by a family of
3-dimensional spacelike hypersurfaces $\Sigma_{t}$, so that the EFEs can be
decomposed into a set of four elliptic-type constraint equations and six
evolution equations. [24]

##  3 The field equations and their vacuum solution for the exterior region of
the spacetime

In accordance with definition (1), the covariant components of the metric
tensor that is associated with the line element (5) is represented by

$g_{\mu\nu}=\text{diag}(-e^{2\nu},\,e^{2\mu},\,e^{2\mu}r^{2},\,e^{2\mu}r^{2}\sin^{2}\theta).$
(11)

The 9 independent Christoffel symbols, out of the 13 non-zero components of
the metric connection, are expressed in the isotropic coordinate frame by

$\displaystyle\displaystyle\Gamma^{t}_{tr}=\nu^{\prime},\quad\Gamma^{r}_{tt}=e^{2\nu-2\mu}\nu^{\prime},\Gamma^{r}_{rr}=\mu^{\prime},\quad\Gamma^{\theta}_{r\theta}=\mu^{\prime}+\frac{1}{r},$
$\displaystyle\displaystyle\Gamma^{r}_{\theta\theta}=-r(r\mu^{\prime}+1),\quad\Gamma^{r}_{\phi\phi}=-r\sin^{2}\theta(r\mu^{\prime}+1),$
$\displaystyle\displaystyle\Gamma^{\theta}_{\phi\phi}=-\sin\theta\cos\theta,\quad\Gamma^{\phi}_{r\phi}=\mu^{\prime}+\frac{1}{r},\quad\Gamma^{\theta}_{\theta\phi}=\cot\theta,$

where the primes indicate differentiation with respect to $r$. One may, now,
continue with the straightforward computation of the covariant components of
curvature tensor $R_{\mu\nu}$ (which we do not write down in the interests of
brevity) with the purpose of constructing the Einstein field equations

$R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}=8\pi T_{\mu\nu},$ (13)

where $R$ is the Ricci scalar, $g_{\mu\nu}$, and $T_{\mu\nu}$ are the
covariant components of the metric tensor (11), and the energy–momentum tensor
[cf. eqs. (15) and (20)], respectively.

From all the above, one may now determine the non-vanishing mixed-variant
components of the EFEs in isotropic coordinates555These equations, which first
appeared in a slightly different form in Eddington’s pioneering monograph [1],
can be obtained by changing the signature from ‘Landau–Lifshitz sign
convention’ to ‘Pauli convention’ and by setting $2\nu\to\nu$ and
$2\lambda\to\mu$ in Eddington’s equations (43.5). as

$\displaystyle\displaystyle 8\pi
T^{t}_{t}=-e^{2\nu-2\mu}\frac{2r\mu^{\prime\prime}+\mu^{\prime}(r\mu^{\prime}+4)}{r},$
(14a) $\displaystyle\displaystyle 8\pi
T^{r}_{r}=\frac{\mu^{\prime}(r\mu^{\prime}+2r\nu^{\prime}+2)+2\nu^{\prime}}{r},$
(14b) $\displaystyle\displaystyle 8\pi T^{\theta}_{\theta}=8\pi
T^{\phi}_{\phi}=r\left[r(\mu^{\prime\prime}+\nu^{\prime\prime})+\mu^{\prime}+r{\nu^{\prime}}^{2}+\nu^{\prime}\right],$
(14c)

where the $(\theta,\,\theta)$-component is identical to the
$(\phi,\,\phi)$-component.

In the exterior vacuum region ($r>R$), all the mixed-variant components of the
energy–momentum tensor are identically zero, that is

$T_{\mu\nu}=0.$ (15)

Therefore, the field equations (14) can readily be integrated to yield the
Schwarzschild solution [1], expressed in isotropic polar coordinates
$(t,\,r,\,\theta,\,\phi)$ by

$ds^{2}=-\left(\frac{1-M/2r}{1+M/2r}\right)^{2}dt^{2}+\left(1+\frac{M}{2r}\right)^{4}\left[dr^{2}+r^{2}(\sin^{2}\theta
d\phi^{2}+d\theta^{2})\right],$ (16)

where $M$ denotes the total gravitational mass of the massive body. In
contrast to its much better-known counterpart, given in curvature coordinates
$(\tilde{t},\,\tilde{r},\,\tilde{\theta},\,\tilde{\phi})$ by

$ds^{2}=-\left(1-\frac{2M}{\tilde{r}}\right)dt^{2}+\left(1-\frac{2M}{\tilde{r}}\right)^{-1}d\tilde{r}^{2}+\tilde{r}^{2}(\sin^{2}\tilde{\theta}d\tilde{\phi}^{2}+d\tilde{\theta}^{2}),$
(17)

the exterior region $r>M/2$ in the isotropic coordinates corresponds to that
of $\tilde{r}>2M$ in curvature coordinates. The limit between the two regions
is represented by a black dashed line in Figures 2–3. The two systems of
coordinates are related to one another via the relations (cf. (B12–B13) in
Ref. [22])

$\tilde{r}=r\left(1+\frac{M}{2r}\right)^{2}\quad\text{and}\quad
r=\frac{\tilde{r}}{2}\left[1+\left(1-\frac{2M}{\tilde{r}}\right)^{1/2}-\frac{M}{\tilde{r}}\right],$
(18)

depicted by the dashed black line in Fig. 1, where the rest of the coordinates
are identical: $\tilde{t}=t,\,\tilde{\theta}=\theta,\,\tilde{\phi}=\phi$ as it
is explained in more detail in A. Accordingly, the value of the isotropic
radial coordinate $r$ taken at the star’s surface is expressed by

$R=\frac{\tilde{R}}{2}\left[1+\left(1-\frac{2M}{\tilde{R}}\right)^{1/2}-\frac{M}{\tilde{R}}\right]$
(19)

if the constants $R$ and $\tilde{R}$ denote the surface radius in curvature
and isotropic coordinates, respectively.

##  4 The hydrostatic-equilibrium solution for the interior region of the
spacetime

Figure 2: The radial dependence of the spherical polar metric potentials. Each
coloured line correspond to one configuration of relativistic stars,
constructed with the SFHo/DD2 EoS for different values of gravitational mass,
which were marked by the bullet point symbol $\bullet$ in Figure 4b. In the
interior region of the spacetime, the expressions $\exp(2\tilde{\nu})$ and
$\exp(-2\tilde{\lambda})$ are represented by one curve with a darker, and one
with a lighter shade of the same colour, respectively. At and beyond the
stellar surface, located at $\tilde{r}=\tilde{R}$, the Schwarzschild exterior
solution (17) (dashed black line) is valid.

In the interior region ($r<R$), let us consider the stress–energy tensor for a
perfect fluid in thermodynamic equilibrium, which can be written as

$T_{\mu\nu}=(\varepsilon+p)u_{\mu}u_{\nu}+pg_{\mu\nu},$ (20)

with $u^{\mu}$ being the components of the fluid’s four-velocity, and the
hydrostatic pressure $p$ is related to the mass–energy density $\varepsilon$
through an equation of state. The relations connecting the above quantities of
a perfect fluid with the metric of a static, spherically symmetric spacetime
(11) are directly obtained by inserting the appropriate components of the
stress–energy tensor (20) into the field equations (14), so that we obtain
three independent ODEs as

$\displaystyle\displaystyle
8\pi\varepsilon=-\frac{e^{-2\mu}}{r}\left[2r\mu^{\prime\prime}+\mu^{\prime}(r\mu^{\prime}+4)\right],$
(21a) $\displaystyle\displaystyle 8\pi
p=\frac{e^{-2\mu}}{r}\left[\mu^{\prime}(r\mu^{\prime}+2r\nu^{\prime}+2)+2\nu^{\prime}\right],$
(21b) $\displaystyle\displaystyle 8\pi
p=\frac{e^{-2\mu}}{r}\left[r(\mu^{\prime\prime}+\nu^{\prime\prime})+\mu^{\prime}+r{\nu^{\prime}}^{2}+\nu^{\prime}\right].$
(21c)

The reason for two separate expressions (21b) and (21c) connecting the
pressure with the metric lies in the fact that the metric tensor (11) is
sufficiently general to apply to any static and spherically symmetric
distribution, whereas having a perfect fluid involved, implies an equality at
each point of the fluid between the stresses in the radial and tangential
directions.[2]

Let us first consider expressing the second derivatives of the metric tensor
by their first derivatives and by thermodynamic quantities. This can be most
conveniently achieved by re-arranging the eq. (21a) in such a way that yields

$\mu^{\prime\prime}=-4\pi
e^{2\mu}\varepsilon-\frac{2\mu^{\prime}}{r}-\frac{\mu^{\prime 2}}{2}$ (22)

and by substituting this expression into eq. (21c) in order to obtain

$\nu^{\prime\prime}=4\pi
e^{2\mu}(\varepsilon+2p)+\frac{\mu^{\prime}(r\mu^{\prime}+2)+2\nu^{\prime}(r\nu^{\prime}+1)}{2r}.$
(23)

Let us use the last remaining EFE, eq. (21b), to express $\nu^{\prime}$ by

$\nu^{\prime}=\frac{8\pi
e^{2\mu}rp+\mu^{\prime}(r\mu^{\prime}+2)}{2(r\mu^{\prime}+1)},$ (24)

where $\mu^{\prime}$ is a functions of $r$ yet to be determined. On account of
the fact that eq. (21b) does not involve higher than first-order derivatives
of the metric tensor, an additional equation can be gained by a
derivation666Instead of derivation, one may alternatively compute the equation
of hydrostatic equilibrium (25) from the conservation of the fluid’s energy
and momentum which demands the covariant divergence of the density- and flux-
tensor field to vanish: $\nabla^{\mu}T_{\mu\nu}=0.$ Taking into account that
the spacetime is static implies that
$\partial_{t}\varepsilon=\partial_{t}p=0$, whereas the isotropy implies that
$\partial_{\theta}p=\partial_{\psi}p=0$. Consequently,
$\nabla_{\mu}T^{\mu}_{t}=-p^{\prime}-\frac{1}{2}(p+\varepsilon)\nu^{\prime}=0$,
as it was first derived by Tolman, Oppenheimer and Volkoff [5, 6] in the
spherical polar coordinate frame. with respect to $r$:

$p^{\prime}=-\frac{1}{2}(p+\varepsilon)\frac{8\pi
e^{2\mu}rp+\mu^{\prime}(r\mu^{\prime}+2)}{r\mu^{\prime}+1},$ (25)

where the second derivatives and $\nu^{\prime}$ were eliminated by the
substitution of eqs. (22–24).

Figure 3: The radial dependence of the isotropic metric potentials. Each
coloured line correspond to one configuration of relativistic stars,
constructed with the SFHo/DD2 EoS for different values of gravitational mass,
which were marked by the bullet point symbol $\bullet$ in Figure 4b. In the
interior region of the spacetime, the expressions $\exp(\nu)$ and
$\exp(-\mu/2)-1$ are represented by a curve with a darker and another curve
with a lighter shade of the same colour, respectively. At and beyond the
stellar surface, located at $r=R$, the Schwarzschild exterior solution (16)
(dashed black line) is valid.

One might reasonably anticipate finding an expression for $\mu^{\prime}$ which
would involve the gravitational mass $m$ by taking into account that eq. (22)
poses a direct relationship between the mass–energy density $\varepsilon$ and
the metric potential $\mu$ and its derivatives, without involving either
$\nu,\,p$ or their respective derivatives. Instead of directly using the
corresponding eq. (21a), one may as well formally integrate the
$(t,\,t)$-component of the EFE, given in curvature coordinates by
$8\pi\varepsilon=e^{-2\tilde{\lambda}}\left(\tilde{r}^{-1}\tilde{\lambda}^{\prime}-\tilde{r}^{-2}\right)+\tilde{r}^{-2}$,
with the limiting condition $\tilde{\lambda}|_{\tilde{r}=0}=0$. The
integration yields the expression

$e^{-2\tilde{\lambda}}=1-\frac{8\pi}{\tilde{r}}\int^{\tilde{r}_{1}}_{0}\varepsilon\tilde{r}^{2}d\tilde{r}$
(26)

for $\tilde{\lambda}(\tilde{r}_{1})$ at an arbitrary point
$\tilde{r}_{1}\leq\tilde{R}$ inside a spherical star of surface radius
$\tilde{R}$. In agreement with the vacuum solution (17), the metric potential
$\tilde{\lambda}$ is connected to the total-gravitational mass of the star $M$
by the relation $e^{-2\tilde{\lambda}}=1-2M/\tilde{R}$, which in turn implies
that it is rational to define the total-gravitational mass as

$M=\int^{\tilde{R}}_{0}dm=4\pi\int^{\tilde{R}}_{0}\varepsilon\tilde{r}^{2}d\tilde{r}.$
(27)

The differential gravitational mass $dm$ contained in the infinitesimal volume
of spherical shell between $r$ and $r+dr$ is given by
$dm=4\pi\varepsilon\tilde{r}^{2}d\tilde{r}=4\pi\varepsilon(e^{2\mu}r^{2})(e^{\mu}[r\mu^{\prime}+1]dr)$,
in accordance with the coordinate transformation (6).777Fig. 2 illustrates the
exponential growth of $\tilde{\lambda}$ by depicting the decrease of the
expression $\exp(-2\tilde{\lambda})$ in (26) as the gravitational mass (27)
piles up when one gradually moves from the centre towards the surface. The
higher the rate of change, the more compact stellar configuration. Upon
recognising the first bracket as the square of the term $e^{\mu}r$ and the
last bracket as the first derivative of the same term, one can considerably
simplify the previous expression and obtain

$\frac{dm}{dr}=\frac{4\pi\varepsilon}{3}\frac{d}{dr}(e^{\mu}r)^{3}.$ (28)

Taking into account that the term $e^{\mu}r$ is equivalent to the areal radius
$\tilde{r}$, the expression $d\tilde{A}\equiv d(e^{\mu}r)^{3}/dr$ evidently
conveys the meaning of a surface element of the 2-dimensional spatial surface.
Consequently, the gravitational-mass element $dm$ in (28) preserved the
physical interpretation as being the integral functional of the mass–energy
density $\varepsilon$ over a volume element $d\tilde{V}\equiv d(e^{\mu}r)^{3}$
in the 3-dimensional space.

Let us emphasise that it is misleading to assume that, in contrast to the
relationship (26) between $\lambda$ and $\varepsilon$ that is apparent in the
context of curvature coordinates, a direct relationship connecting the metric
potential $\mu$ and the energy density $\varepsilon$ cannot established for
the isotropic coordinates. We demonstrate in B that the two quantities are
linked together by a polynomial equation of fourth degree in the variable
$\exp(2\mu)$, which also involves the radical expressions of $\exp(2\mu)$ and
$\epsilon$.

###  4.1 Representation of the Tolman–Oppenheimer–Volkoff equation in
isotropic coordinates

Combining eqs. (26–27) with (9) allows us to replace $\mu^{\prime}$ by the
gravitational mass $m$ in those terms of eqs. (24–25), namely,

$\mu^{\prime}(r\mu^{\prime}+2)=\frac{2m}{e^{\mu}r^{2}},$ (29)
$r\mu^{\prime}+1=\left(1-\frac{2m}{e^{\mu}r}\right)^{1/2},$ (30)

which involved first derivatives of the metric-potential function
$\mu^{\prime}$. Eq. (24) is then re-expressed as

$\nu^{\prime}=\left(4\pi
e^{2\mu}rp+\frac{m}{e^{\mu}r^{2}}\right)\left(1-\frac{2m}{e^{\mu}r}\right)^{-1/2},$
(31)

where $\nu$ acts as an effective relativistic analogue of the Newtonian
gravitational potential. Having the derivatives eliminated from eq. (25) in a
similar manner, the condition for hydrostatic equilibrium takes on the form of

$p^{\prime}=-\frac{Gm}{e^{\mu}r^{2}}(p+\varepsilon)\left(\frac{4\pi}{c^{4}}\frac{(e^{\mu}r)^{3}p}{m}+1\right)\left(1-\frac{2G}{c^{2}}\frac{m}{e^{\mu}r^{2}}\right)^{-1/2},$
(32)

with the values of $G$ and $c$ having been restored for comparison’s sake. The
corresponding Tolman–Oppenheimer–Volkoff equation in curvature coordinates is
given by $\displaystyle
p^{\prime}=-\frac{Gm}{\tilde{r}^{2}}(p+\varepsilon)\left(\frac{4\tilde{r}^{3}p}{mc^{2}}+1\right)\left(1-\frac{2Gm}{\tilde{r}c^{2}}\right)^{-1}$.
By careful examination of eq. (32) in contrast to its better-known
counterpart, the exponent $1/2$ over the last bracket is exposed to be
different. Both expressions of the TOV equation are, therefore, of the so-
called Riccati-type, a first-order non-linear ODE that is quadratic in
isotropic pressure, irrespective of the choice of coordinate frame.

##  5 Numerically computing equilibrium sequences of relativistic stars

(a) Isotropic pressure against baryon-number density in the crust.

(b) Total gravitational mass–central energy density relations.

Figure 4: Left panel: The pressure is shown as a function of baryon-number
density / rest-mass density for the three EoSs listed in Table 1. The nuclear
saturation density $n_{0}\simeq 0.16\text{ fm}^{-3}$ is denoted by a shaded
gray line. The range of pressures at $n_{0}$ is approximately a factor of $1$
to $3$ $\text{MeV fm}^{-3}$. Right panel: The corresponding mass–radius
relation is displayed for each tabulated EoS with the orange curve
representing the sequence of relativistic stars constructed with the EoS
SFHo/DD2, whereas the blue and the green curves, respectively, denoting a
stiffer EoS (GMSR(H4)) and a softer EoS (DS(CMF)-2 Hybrid) based on different
microphysics, both of them are publicly accessible in CompOSE. The symbols
$\bullet$ mark a set of five selected configurations of relativistic stars
that were computed for different values of gravitational mass with the use of
the EoS SFHo/DD2; the maximal-mass configuration with $M=2.13M_{\odot}$, an
extremely low-mass configuration with $M=0.19M_{\odot}$, and three equidistant
points in terms of the gravitational mass with
$M=0.68M_{\odot},\,1.16M_{\odot},\,1.65M_{\odot}$, respectively. For this set
of configurations, the relationship between the two radial coordinates and the
behaviour of the metric potentials are shown in Figures 1–3.

The TOV equation (32), together with the expression for the gravitational mass
(28)888It should be pointed out that, although the differential eq. (28) can
be replaced by an algebraic relationship connecting $\mu$ and $\varepsilon$ in
the form of the quartic equation (55), the complexity of the algebraic
equation makes integrating the differential equation much more practical in
terms of computational cost. and the equations for the metric potentials
(30–31), constitutes a partially uncoupled system of first-order, non-linear
and non-autonomous ODEs [25], compiled into the following form:

$\displaystyle\displaystyle\frac{dm}{dr}=\frac{4\pi\varepsilon}{3}\frac{d}{dr}(e^{\mu}r)^{3},$
(33a)
$\displaystyle\displaystyle\frac{d\mu}{dr}=\frac{1}{r^{2}}\left[\left(1-\frac{2m}{e^{\mu}r}\right)^{1/2}-1\right],$
(33b) $\displaystyle\displaystyle\frac{d\nu}{dr}=\left(4\pi
e^{2\mu}rp-\frac{m}{e^{\mu}r^{2}}\right)\left(1-\frac{2m}{e^{\mu}r}\right)^{-1/2},$
(33c)
$\displaystyle\displaystyle\frac{dp}{dr}=-\frac{m}{e^{\mu}r^{2}}(p+\varepsilon)\left(\frac{4\pi(e^{\mu}r)^{3}p}{m}+1\right)\left(1-\frac{2m}{e^{\mu}r^{2}}\right)^{-1/2}.$
(33d)

Supplemented with an EoS that is set out in Table 1, the ODE system (33)
completely determines the stellar structure in gravitational hydrostatic
equilibrium through performing the integration from the centre of the star
with initial condition $m(0)=0$, where an arbitrary value for the central
energy-density $\varepsilon_{c}=\varepsilon(0)$, until the pressure $p(r)$
vanishes at some value $r=R$ which corresponds to the position of the star’s
surface.999 Let us remember that surface radius corresponds to different
values in curvature and isotropic coordinates. The relation (19) links the two
values associated with the star’s surface. A unique sequence of relativistic
stars in hydrostatic equilibrium corresponds to each EoS, parameterised by
$\varepsilon_{c}$. Figure 4 shows the corresponding EoS (left) and mass–radius
relations (right), respectively.

With the intention of computing configurations of relativistic stars with
considerably large gravitational mass, we used the same EoS (SFHo/DD2) in this
paper as in our earlier paper [26], where, although isotropic coordinates were
adapted in the LORENE/nrotstar code [15], outlining the computational
differences stemming from the choice of coordinate frame was not placed at the
focus of the study. This zero-temperature EoS is based on a vector-meson
extended linear $\sigma$-model that combines two relativistic mean-field (RMF)
models, the softer SFHo model developed by [27] and the stiffer DD2 model by
[28, 29]. In an effort to mitigate the risk of drawing inaccurate conclusions
about certain numerical features of our computations that may arise from the
inherent properties of a single EoS, we also measured the computational time
of equilibrium configurations built with a softer and a stiffer EoS, DS(CMF)-2
Hybrid and GMSR(H4), respectively. Both tabulated EoS data are available in
the CompOSE online repository [30]. DS(CMF)-2 Hybrid [31] is a relativistic
SU(3) chiral mean-field (CMF) model using a non-linear realisation of the
$\sigma$-model, which includes pseudo-scalar mesons as the angular parameters
for the chiral transformation. GMSR(H4) [32] corresponds to a unified EoS of
cold nuclear matter at $\beta$-equilibrium, where the effective interaction is
determined by the non-relativistic meta-model adjusted to reproduce the chiral
EFT Hamilatonian H4.

Name | Particles | $\bm{M_{\text{{max}}}}$ | $\bm{\min(n_{\text{{b}}})}$ |  $\bm{\max(n_{\text{{b}}})}$ | Ref.
---|---|---|---|---|---
GMSR(H4) | Nucleonic | $2.33M_{\odot}$ | $10^{-7}\text{ fm}^{-3}$ | 2.00 fm-3 | [32]
SFHo/DD2 | Nucleons | $2.13M_{\odot}$ | $10^{-9}\text{ fm}^{-3}$ | 1.86 fm-3 | [27]
DS(CMF)-2 Hybrid | Hybrid | $1.96M_{\odot}$ | 0.03 $\text{fm}^{-3}$ | 1.6 $\text{fm}^{-3}$ | [31]
|  |  |  |  |
Table 1: Some characteristic parameters of the two nucleonic and one hybrid
nucleon–hyperon–quark matter models used in the zero-temperature equations of
state displayed in Fig. 4a.

The numerical performance of the computational scheme proposed in this paper
is qualified in two respects. First, it is qualified by the computation time
spent on numerically solving the TOV equations (33) for each of the three
tabulated EoSs listed in Table 1. Secondly, the accuracy of the numerical
results is assessed in terms of the absolute difference in the gravitational
mass computed in the proposed scheme compared to the one computed in LORENE.
The computation time was measured as the average CPU execution time
$t_{\text{CPU}}$, averaged over 1000 runs on an equally-spaced sequence in
logarithm scale of the parameter $\varepsilon_{c}$ in the range that lies
between $100\text{ MeV/fm}^{3}$ and the central energy density of the maximum-
mass configuration, with respect to the particular EoS. Figure 5a exposes that
the length of time required to perform a computation falls in the range of
$0.7$ to $3.5$ seconds, where the evaluation, as a general rule, is faster for
stiffer EoSs. The computational time for stars with comparatively low or high
$\varepsilon_{c}$ stays constant, however, in the intermediate region, the
figure features a trend of increase in $t_{\text{CPU}}$ with increasing
$\varepsilon_{c}$ for the GMSR(H4) and the DS(CMF)-2, but decrease for
SFHo/DD2. The absolute difference in the gravitational mass given by $\Delta
M(\varepsilon_{c})\equiv
M_{\text{iso}}(\varepsilon_{c})-M_{\texttt{LORENE}}(\varepsilon_{c})$ was
measured for each EoS in the same range of central energy density as the
computation time was. Figure 5b shows that, following an initial rapid growth
in the low-energy-density region, $\Delta M$ increases steadily with
$\varepsilon_{c}$ (and therefore with $M$ as well) until the rate of growth
once again drops as the central energy density associated with the maximum
mass is approached, irrespective of the particular EoS. Another characteristic
feature of the $\Delta M$ that is apparent in the figure is that, aside from
the $\varepsilon_{c}$, it also increases with stiffness of the EoS.

(a) The computation time for evaluating equilibrium configurations, measured
as the average CPU execution time.

(b) The absolute difference in the gravitational mass between the one computed
in the proposed scheme and in LORENE.

Figure 5: The performance of computations carried out in the proposed scheme.
The computation time required for evaluating equilibrium configurations of
relativistic stars for each EoS (left panel) and the absolute difference in
the gravitational mass in $M_{\odot}$ units (right panel), are shown as
functions of the central energy density $\varepsilon_{c}$. The bullet point
symbols $\bullet$ represent the same set of configurations with varying
gravitational mass as in Fig. 4b.

##  6 Conclusions

The results of this study can be divided into an analytical part and a
computational part:

1\. By solving the EFEs set up in isotropic coordinates, we have obtained the
equation of hydrostatic equilibrium (32) which, apart from an exponent and a
minus sign, appears unchanged from the familiar form in which the TOV equation
appears in terms of curvature coordinates. Complemented by an EoS, the
partially uncoupled system of first-order ODEs (33) uniquely determines the
stellar structure of relativistic stars in hydrostatic equilibrium. We have
demonstrated in B that the differential equation (33a) for the gravitational
mass appearing in the ODE system (32) can be replaced by an algebraic
relationship connecting $\mu$ and $\varepsilon$ in the form of the quartic
equation (55), but due to the complexity of the algebraic equation, the
integration of eq. (33a) is much more practical in terms of the computational
cost.

In order to contextualise our analytical results, it is worth comparing them
with an analysis based on rescaling the radial profile, carried out by
Tsokaros, Uryū and Rezzolla [22]. In order to bring the notations of the
quoted article in line with our own, let us change the notations $r,\,\bar{r}$
used for the radial coordinates to our own $\tilde{r},\,r$. Similarly, let us
replace the notation $A,\,B,\,\alpha,\,\psi$, used for metric potentials to
$\exp(2\tilde{\nu}),\,\exp(2\tilde{\lambda}),\exp(\nu)$ and $\exp(\mu/2)$,
respectively. This article demonstrated that the derivation of new structure
equations can be avoided by transforming the ODE system of stellar structure
from Schwarzschild to isotropic coordinates by rescaling the radial coordinate
in accordance with the conversion formula (41). Eq. (41) is the result of eq.
(6) and eq. (26), which defines the relation of the Schwarzschild metric
potential $\tilde{\lambda}$ to the gravitational mass as
$e^{-2\tilde{\lambda}}=1-2m(\tilde{r})/\tilde{r}$. Without the involvement of
any relation to the gravitational mass, this transformation would revert back
to our relation (7).

Ref. [22] argues that, therefore, one could always carry out the calculations
in curvature coordinates and then rescale the radial profile of the solution
so that the stellar surface would appear in the correct position. Eq. (41)
allows one to rewrite the structure equations (B8–B11) appearing in [22] to
the form

$\displaystyle\displaystyle\frac{dm}{dr}=4\pi
e^{\mu}r^{2}\left(1-\frac{2m}{e^{\mu}r}\right)^{1/2},$ (34a)
$\displaystyle\displaystyle\frac{d\mu}{dr}=\left[\left(1-\frac{2m}{e^{\mu}r}\right)^{1/2}-1\right],$
(34b) $\displaystyle\displaystyle\frac{d\nu}{dr}=\left(4\pi
e^{2\mu}rp-\frac{m}{e^{\mu}r^{2}}\right)\left(1-\frac{2m}{e^{\mu}r}\right)^{-1/2},$
(34c)
$\displaystyle\displaystyle\frac{dp}{dr}=-\frac{m}{e^{\mu}r}(p+\varepsilon)\left(\frac{4\pi(e^{\mu}r)^{3}p}{m}+1\right)\left(1-\frac{2m}{e^{\mu}r}\right)^{-1/2},$
(34d)

that are fully consistent with eqs. (33). Here, eq. (34a) corresponds to the
expression (33a) and the eqs. (34b)–(34c) were derived with using eq. (6).

2\. A further objective of this paper was to evaluate the performance of the
computational scheme proposed here in terms of computational time and
numerical accuracy. We used the same EoS (SFHo/DD2) as in our earlier paper
[26], in which paper isotropic coordinates were implemented in the
LORENE/nrotstar code, however, it did not consider the impact of different
choices of coordinate frame on the computational cost. To avoid
unsubstantiated inferences about certain numerical aspects of our
computations, which may arise from using a single EoS, two additional EoS
models were also included in the study: a softer EoS (DS(CMF)-2 Hybrid) and a
stiffer EoS (GMSR(H4)), both of which are available in the CompOSE online
repository. Our results indicate that the computational time required to
numerically solve the TOV equations tends to increase with the stiffness of
the employed EoS. The computational time for stellar models with comparatively
low or high $\varepsilon_{c}$ remains constant, but at intermediate values of
$\varepsilon_{c}$ the results suggest a tendency for a decrease in the
computational time for the SFHo/DD2, whereas an increase for the two other EoS
models. Another aspect of evaluating the performance of the present
computational scheme was to measure the absolute difference in the
gravitational mass between the implementation of this scheme and the results
computed by LORENE code for each EoS over the same region of central energy
density as the computation time was measured. Our results indicate that after
an initial rapid growth in the low energy density region, the absolute
difference steadily increases with the central energy density until the rate
of growth decreases as the central energy density associated with the maximum
mass is reached, irrespective of the particular EoS. Another defining feature
of $\Delta M$ is that it not only increases with central energy density, but
also with stiffness of the EoS.

## Acknowledgments

We wish to acknowledge the computational resources provided by the Wigner
Scientific Computational Laboratory (WSCLAB) (formerly known as Wigner GPU
Laboratory). This research has received partial support from the National
Research, Development and Innovation Office (Hungarian abbreviation: NKFIH)
under OTKA Grant Agreement No. K138277.

## Data availability statement

The data underlying this paper will be shared on reasonable request to the
corresponding author.

##  Appendix A Labelling the coordinate frame as “isotropic” and its relation
to the Schwarzschild-coordinate frame

Provided that the spherical coordinates are converted to Cartesian coordinates
by the transformation

$(r,\,\theta,\,\phi)\mapsto(x^{1},\,x^{2},\,x^{3}):\begin{cases}x^{1}=r\sin\theta\cos\phi\\\
x^{2}=r\sin\theta\sin\phi\\\ x^{3}=r\cos\theta,\end{cases}$ (35)

the Schwarzschild metric (16), originally expressed in spherical coordinates,
assumes the form of

$ds^{2}=-\left(\frac{1-M/2r}{1+M/2r}\right)^{2}dt^{2}+\left(1-\frac{M}{2r}\right)^{4}\left[(x^{1})^{2}+(x^{2})^{2}+(x^{3})^{2}\right],$
(36)

where $r=\sqrt{(x^{1})^{2}+(x^{2})^{2}+(x^{3})^{2}}$ (e.g. [1], p. 93). Light-
like world lines correspond to the line element $ds^{2}=0$ and are, therefore,
easily expressed in the Cartesian-coordinate frame by

$\left(\frac{dx^{1}}{dt}\right)^{2}+\left(\frac{dx^{2}}{dt}\right)^{2}+\left(\frac{dx^{3}}{dt}\right)^{2}=\frac{(1-M/2r)^{2}}{(1+M/2r)^{6}},$
(37)

which implies that the coordinates $(t,\,x^{1},\,x^{2},\,x^{3})$ are
rightfully called “isotropic” on account of the components of the four-
velocity field $u^{\mu}=dx^{\mu}/dt$ that are identical in all spatial
directions. Evidently, one of the defining features of isotropic line elements
is that it exposes uniformity regardless of the orientation in space.

On account that the isotropic space is obtained from the Euclidean 3-space
$\mathbb{R}^{3}$ by substituting the usual Euclidean distance with the
isotropic distance [33], the metric (5) in isotropic coordinates has its
spacelike slices conformal to the Euclidean one, where

$d\Sigma^{2}=dr^{2}+r^{2}(\sin^{2}\theta d\phi^{2}+d\theta^{2})$ (38)

is the Euclidean metric. So let us substitute $\tilde{r}$ by $r$ and write
down our metric:

$ds^{2}=-\left(1-\frac{2M}{r}\right)dt^{2}+e^{2\mu}[dr^{2}+r^{2}(\sin^{2}\theta
d\phi^{2}+d\theta^{2})].$ (39)

By comparing this form of the metric (39) in curvature coordinates with that
of (17), and requiring both the radial and the angular parts to be equal,
respectively, we obtain eq. (6). In order to eliminate the metric function
$e^{2\mu}$, the first equation in (6) is divided by the second one, so that

$\frac{dr^{2}}{r^{2}}=\frac{d\tilde{r}^{2}}{\tilde{r}^{2}-2M\tilde{r}}.$ (40)

Finally, taking the square root of eq. (40) and integrating it will yield the
transformation (18). Rearranging eq. (40) and replacing $M$ by $m(\tilde{r})$
yields

$\frac{d\tilde{r}}{dr}=\frac{\tilde{r}}{r}\sqrt{1-\frac{2m(\tilde{r})}{\tilde{r}}},$
(41)

which acts as the relationship between the two radial coordinates inside the
star.

##  Appendix B Seeking an algebraic relationship between the metric potential
$\mu$ and the energy density $\epsilon$

The combination of the eqs. (26–27) with (9) yields an integro-differential
equation

$(r\mu^{\prime}+1)^{2}=1-\frac{8\pi}{e^{\mu}r}\int\varepsilon\cdot
e^{3\mu}r^{2}(r\mu^{\prime}+1)dr,$ (42)

then performing the differentiation on both sides results in

$e^{\mu}r^{2}\mu^{\prime}\left(2r\mu^{\prime\prime}+r{\mu^{\prime}}^{2}+3\mu^{\prime}\right)=e^{\mu}(r\mu^{\prime}+1)-8\pi\varepsilon\cdot
e^{3\mu}r^{2}(r\mu^{\prime}+1),$ (43)

where the term $\mu^{\prime\prime}$ can easily be eliminated by eq. (22). As a
result, the previous expression turns into an algebraic cubic equation in a
single variable $\mu^{\prime}$, which can be re-arranged into the standard
form

$r^{3}{\mu^{\prime}}^{3}+4r^{2}{\mu^{\prime}}^{2}+r{\mu^{\prime}}+1-8\pi\varepsilon
e^{2\mu}r^{2}=0.$ (44)

Now, restoring the ordinary form of the metric-potential function (3) and
expressing it, together with its first derivative by

$\mu(r)=\frac{1}{2}\log
B(r),\quad\mu^{\prime}(r)=\frac{B^{\prime}(r)}{2B(r)},$ (45)

allows one to eliminate its exponential form from eq. (44) and re-write it as
a general cubic equation

$a{B^{\prime}}^{3}+b{B^{\prime}}^{2}+cB^{\prime}+d=0\quad\text{with}\quad\begin{cases}a=r^{3},\\\
b=8r^{2}B,\\\ c=4r{B}^{2},\\\ d=8B(1-8\pi\varepsilon r^{2}B).\end{cases}$ (46)

After dividing by $a$ and changing the variable $B^{\prime}$ to

$t=B^{\prime}+\frac{b}{3a},$ (47)

eq. (46) can be transformed into a depressed cubic equation

$t^{3}+pt+q=0,$ (48)

where the second-degree term is missing and the coefficients of the polynomial
are real numbers, given by

$p=\frac{3ac-b^{2}}{3a^{2}},\quad q=\frac{2b^{3}-9abc+27a^{2}d}{27a^{3}}.$
(49)

Employing Vieta’s transformation $t\equiv w-\frac{p}{3w}$ and multiplying by
$w^{3}$ transforms the depressed cubic into a quadratic equation

$(w^{3})^{2}+q(w^{3})-\frac{p^{3}}{27}=0.$ (50)

(a) The discriminant $\Delta$ on log–log scale.

(b) The derivative $B^{\prime}(r)$ as the real root of the cubic equation
(46).

Figure 6: The discriminant $\Delta$ and the derivative $B^{\prime}$ are shown
on panel (a) and (b), respectively, as functions of $r/M$. Each coloured line
corresponds to one configuration of relativistic stars, constructed with the
SFHo/DD2 EoS for different values of gravitational mass, which were marked by
the bullet point symbol $\bullet$ in Figure 4b. The dashed black line
represents $B^{\prime}(r\geq R)=(M/r)\left(1+M/2r\right)$ at and beyond the
stellar surface, located at $r=R$.

By using a further substitution $u\equiv w^{3}$, the quadratic equation admits
two non-zero solutions

$u_{\pm}=-\frac{q}{2}\pm\sqrt{\Delta},$ (51)

corresponding to the sign of the discriminant $\Delta=q^{2}/4+p^{3}/27>0$. The
discriminant, shown on a log–log scale in Figure 6a, appears as an
approximately linear decreasing function of radial coordinate $r$. It has a
positive intercept, which is greater for stars with larger gravitational mass
and a constant, negative slope that does not depend on the mass. Thus, the
depressed cubic equation (48) has one real root and two complex roots. The
real root, expressed by Cardano’s formula, is the sum of the cube roots (51):

$w_{1}=\sqrt[3]{u_{+}}+\sqrt[3]{u_{-}}.$ (52)

Although the conjugate pair of complex roots don’t hold any physical
relevance, they can be obtained by multiplying the cube root by each of the
two primitive cube roots of unity101010 Let $\delta_{1},\,\delta_{2}$ and
$\delta_{3}$ denote the primitive cube roots of unity for $\delta^{3}=1$. This
can be factored into $\delta^{3}-1=(\delta-1)(\delta^{2}+\delta+1)=0$, which
has the three roots $\delta_{0}=1,\,\delta_{1}=\frac{-1+i{\sqrt{3}}}{2}$ and
$\delta_{2}=\frac{-1-i{\sqrt{3}}}{2}$.:

$w_{2}=-\frac{1}{2}(\sqrt[3]{u_{+}}+\sqrt[3]{u_{-}})+\frac{i{\sqrt{3}}}{2}(\sqrt[3]{u_{+}}-\sqrt[3]{u_{-}}),$
(53)
$w_{3}=-\frac{1}{2}(\sqrt[3]{u_{+}}+\sqrt[3]{u_{-}})+\frac{i{\sqrt{3}}}{2}(-\sqrt[3]{u_{+}}+\sqrt[3]{u_{-}}).$
(54)

Reversing the series of transformations implemented above, the real root (52)
provides us with an algebraic relationship directly between $\exp(2\mu)$ and
$\epsilon$ and, thus, renders the necessity of numerically solving the non-
linear ODE (44) obsolete. This algebraic equation of $\exp(2\mu)$ and
$\epsilon$ can be substantially simplified into a form of

$\left(\frac{237276}{\beta^{3}}+\frac{6084\sqrt[3]{2}}{\beta}+\beta^{3}+78\sqrt[3]{4}\beta+184\right)e^{4\mu}-432\pi
r^{2}\varepsilon e^{2\mu}+54=0,$ (55)

which also involves radical expressions of $\exp(2\mu)$ and $\epsilon$ that
are incorporated into the terms $\alpha$ and $\beta$:

$\displaystyle\displaystyle\alpha^{3}=-92e^{4\mu}+216\pi r^{2}\varepsilon
e^{2\mu}-27-\sqrt{\left(92e^{4\mu}-216\pi r^{2}\varepsilon
e^{2\mu}+27\right)^{2}-8788e^{8\mu}},$ (56a)
$\displaystyle\displaystyle\beta=\alpha+\sqrt[3]{\alpha^{3}-432\pi
r^{2}\varepsilon e^{2\mu}}.$ (56b)

## References

## References

  * [1] Eddington A S 1933 The Mathematical Theory of Relativity (Cambridge, UK: Cambridge University Press) ISBN 978-0-521-09165-7
  * [2] Tolman R C 1930 Phys. Rev. 35(8) 904–924
  * [3] Tolman R C and Ehrenfest P 1930 Phys. Rev. 36(12) 1791–1798
  * [4] Lemaître G 1933 Ann. Soc. Sci. Bruxelles 53(2) 51–85
  * [5] Tolman R C 1939 Phys. Rev. 55(4) 364–373
  * [6] Oppenheimer J R and Volkoff G M 1939 Phys. Rev. 55(4) 374–381
  * [7] Herrera L 2016 Tov or ov?, the whole story (Preprint 1611.06364)
  * [8] Semiz I 2017 Reply to: “tov or ov? the whole story” (Preprint 1702.06002)
  * [9] Semiz I 2016 Ov or tov? (Preprint 1610.03049)
  * [10] Delgaty M S R and Kayll L 1998 Comput. Phys. Commun. 115 395–415
  * [11] Stephani H, Kramer D, MacCallum M A H, Hoenselaers C and Herlt E 2003 Exact solutions of Einstein’s field equations Cambridge Monographs on Mathematical Physics (Cambridge, UK: Cambridge University Press) ISBN 978-0-521-46702-5, 978-0-511-05917-9
  * [12] Las Heras C and León P 2019 Eur. J. Phys. 79 990 (Preprint 1905.02380)
  * [13] Nariai H 1950 Sci. Rep. Tohoku Univ. Eighth Ser. 34 951–961
  * [14] Goldman S P 1978 Astrophys. J. 226 1079–1086
  * [15] Gourgoulhon E, Grandclément P, Marck J A, Novak J and Taniguchi K 1997 Lorene URL https://lorene.obspm.fr/
  * [16] Bonazzola S, Gourgoulhon E, Salgado M and Marck J A 1993 Astron. Astrophys. 278 421–443
  * [17] Werneck L et al. 2023 The einstein toolkit URL https://doi.org/10.5281/zenodo.7942541
  * [18] Löffler F, Faber J, Bentivegna E, Bode T, Diener P, Haas R, Hinder I, Mundim B C, Ott C D, Schnetter E, Allen G, Campanelli M and Laguna P 2012 Class. Quantum Grav. 29(11) 115001 ISSN 1361-6382
  * [19] Uryū K and Tsokaros A 2012 Phys. Rev.D 85(6) 064014 (Preprint 1108.3065)
  * [20] Komatsu H, Eriguchi Y and Hachisu I 1989 Mon. Not. R. Astron. Soc. 237 355–379
  * [21] Komatsu H, Eriguchi Y and Hachisu I 1989 Mon. Not. R. Astron. Soc. 239 153–171
  * [22] Tsokaros A, Uryū K and Rezzolla L 2015 Phys. Rev.D 91 104030 (Preprint 1502.05674)
  * [23] Gourgoulhon E 2007 3+1 formalism and bases of numerical relativity (Preprint 0703035)
  * [24] Lap-Ming L and Novak J 2006 Class. Quantum Grav. 23 4545 (Preprint 0603048)
  * [25] Martins Y X, Campos L F A, de Souza Plácido D and Biezuner R J 2019 Ann. Phys., NY 409 167929
  * [26] Kacskovics B, Barta D and Vasúth M 2023 Astron. Nachr. 344 e220109 (Preprint 2212.04885)
  * [27] Steiner A W, Hempel M and Fischer T 2013 Astrophys. J. 774 17 (Preprint 1207.2184)
  * [28] Typel S, Röpke G, Klähn T, Blaschke D and Wolter H H 2010 Phys. Rev.C 81(1) 015803 (Preprint 0908.2344)
  * [29] Hempel M and Schaffner-Bielich J 2010 Nucl. Phys.A 837 210–254 ISSN 0375-9474 (Preprint 0911.4073)
  * [30] CompOSE Core Team, Typel S, Oertel M, Klähn T, Chatterjee D, Dexheimer V, Ishizuka C, Mancini M, Novak J, Pais H, Providência C, Raduta A R, Servillat M and Tolos L 2022 Eur. J. Phys.A 58 221 (Preprint 2203.03209)
  * [31] Dexheimer V and Schramm S 2008 Astrophys. J. 683 943–948 (Preprint 0802.1999)
  * [32] Grams G, Margueron J, Somasundaram R and Reddy S 2022 Eur. J. Phys.A 58 56 (Preprint 2203.11645)
  * [33] Aydin M and Mihai I 2017 Math. Commun. 22 41–51