# Algorithms for Generating All Static Spherically Symmetric (An)isotropic
Fluid Solutions of Einstein’s Equations

M M Akbar111E-mail<EMAIL_ADDRESS>& R Solanki222E-mail:
<EMAIL_ADDRESS>

∗Department of Mathematical Sciences
†Department of Physics
University of Texas at Dallas
Richardson, TX 75080, USA

###### Abstract

We study the Einstein equations of the static spherically symmetric
anisotropic fluid system in curvature coordinates to find algorithms that
generate all solutions and all solutions that are regular at the center. All
possible combinations of input functions from the set of four functions that
characterize the anisotropic system are considered and all equivalent
conditions for central regularity are determined (for both isotropic and
anisotropic systems). We provide the first regularity analysis of the known
algorithm that uses the potential function and anisotropy as inputs. The one-
parameter family of solutions produced by this algorithm and by its isotropic
counterpart are regular at the center for the usual initial conditions on the
input functions; these conditions by themselves are inadequate for geometric
regularity, but the Einstein equations fill the gap. By slightly reformulating
the anisotropic algorithm, all regular solutions can be obtained in a 1–1
relation with the input. This allows us to interpret the parameter for the
isotropic and anisotropic systems as the central density. For three other
choices of input function pairs (any two of the potential function, density,
or radial pressure), a remarkably straightforward algorithm follows, which is
very efficient in generating regular anisotropic solutions. This is because
the equivalency of the three pairs in this algorithm arises precisely from the
same algebraic relation that made the different equivalent sets of regularity
conditions possible. In addition, the choice of functions makes this algorithm
very suitable for finding particular solutions that admit other desirable
physical properties; we construct three examples. This algorithm does not
admit an isotropic limit although all isotropic solutions are produced as part
of the anisotropic system. The remaining two choices of input function pairs
(anisotropy with the radial pressure or density) lead to the old barriers one
encounters in the isotropic system: Riccati and Abel equations. However, with
any solution generated by the new and existing algorithms, one can now
construct the general solution of the corresponding Riccati equation to obtain
a one-parameter family of geometries for each input solution. We discuss the
regularity of the resulting solutions.

##  1 Introduction

The most general spherically symmetric static-fluid configurations possible in
general relativity have anisotropic pressure. Einstein’s equations for this
system reduce to a system of two equations in four unknown functions: density,
radial pressure, (one) tangential pressure, and the potential function. The
difference between the two pressures is the anisotropy function, and when this
vanishes identically, one obtains the system of isotropic fluid
configurations. Isotropic and anisotropic systems, thus, require one and two
input functions, respectively, to produce a solution. One can also use an
equation of state as an input. Isotropic fluid spheres (or perfect fluid
spheres, as they are more commonly known) have been studied extensively since
the early days of relativity and continue to be of physical and mathematical
interest [27, 14, 10]. There is a large corpus of exact solutions of the
isotropic system, which were obtained by employing various ad hoc techniques.
However, only a very few of them can be considered physical [10] and, thus, to
obtain physically reasonable solutions, one generally uses numerical
integration. Despite Einstein’s cluster solution in the late 30s [12],
anisotropic fluid spheres have received much less attention. In recent times,
they have attracted fresh interest, as anisotropy has been found to be
relevant in several physical situations [19], especially in modeling static
galactic dark matter halos (see, for example, [3]) as well as compact stars,
both of which are currently attracting a considerable amount of research
interest [4, 11, 29].

More recently, a number of generation algorithms have been developed that can
obtain all static isotropic fluid spheres by quadrature from a single
function—the generating function—in various coordinates [2, 20, 25, 23]. This
approach is comparable to the well-studied static axisymmetric vacuum system
for which one can generate all solutions from an arbitrary (harmonic)
generating function (see, for example, [27]). An algorithm gives considerable
insight into the structure of the solution space and can be useful in
constructing new solutions, exact or otherwise, of particular interest. There
is also an algorithm for the anisotropic system in canonical coordinates [18],
which extends the perfect fluid algorithm in [20]. The two inputs used in this
algorithm are the potential function (equivalently, the $g_{tt}$ component of
the metric) and the anisotropy. Another algorithm uses two inputs that are
functions of fluid variables to generate anisotropic solutions in canonical
coordinates. A less general algorithm that generates a class of anisotropic
solutions from isotropic Newtonian solutions is known [22]. Generating
functions are not _a priori_ required to be quantities that are readily
physically interpretable and they could be complicated expressions in terms of
the physical quantities of the system (see, for example, [15]). In any case,
even when they can produce all solutions, it is not easy to determine
beforehand what inputs will give a desired output solution (like solutions
that are regular at the center and have other properties of physical
interest). Addressing such issues will be central to this paper.

In this paper, we will systematically study all possible generation algorithms
in canonical coordinates using the basic functions of the (an)isotropic system
(as in the algorithms of [20] and [18]). We will not consider possibilities
using other functions as inputs (as in [15, 23]). A generating algorithm that
finds all solutions, or solutions of a certain type, of the system via
quadrature is not in general not the best way to find _particular_ solutions
and these are generally considered separate questions. However, using the
basic variables of the system, they can be simultaneously approached since
physical conditions in canonical coordinates are given in terms of the basic
variables and their derivatives. As we will see in this paper, regularity can
be considered within a new algorithm which generates all regular solutions in
canonical coordinates. This will streamline the search for particular
solutions that satisfy other physical properties. We have three main goals:
(i) to identify _all_ possible algorithms using these basic variables such
that each algorithm generates all solutions of the system via quadrature, (ii)
to generate all (and only) regular solutions using each algorithm found via
(i),333Note that generating solutions that are regular at the center (i.e.,
scalar polynomial invariants constructed from that Riemann tensor that are
finite) is an altogether different problem from checking whether an existing
solution is regular. There are no discussions in the literature on the
regularity of the algorithm in [18]. We will address this as well as
regularity in all other algorithms we find.  and (iii) to use (i) and (ii) to
find particular examples that satisfy other physical conditions.

This paper is arranged as follows. In Section 2, we introduce the Einstein
equations for the system, and in Section 3 we examine the geometric conditions
of regularity in the light of Einstein’s equations. In Section 4, we study the
existing algorithm(s) and reformulate it in terms of the tangential pressure
function and show how both can produce solutions that are regular at the
center. We then present a new algorithm that generates all solutions from
three other pairs of functions; this algorithm makes generating regular
anisotropic solutions very straightforward. Then, in Section 5, we reconsider
the Riccati equations—which on their own could not be integrated—in light of
the algorithms. In Section 6, we find solutions that satisfy other physical
conditions, staring from the new algorithm. We conclude in Section 7.

##  2 Spherically Symmetric Static Fluid Systems

We will be working in the curvature/Schwarzschild coordinates [24]:

$ds^{2}=-e^{2\Phi(r)}dt^{2}+e^{2\Psi(r)}dr^{2}+r^{2}(d\theta^{2}+\sin^{2}{\theta}d\phi^{2}),$
(2.1)

for which the most general fluid energy–momentum tensor is (see, for example,
[13])

$T^{t}_{t}=-\rho(r),\quad T^{r}_{r}=p(r),\quad
T^{\theta}_{\theta}=T^{\phi}_{\phi}=P(r),$ (2.2)

where $\rho(r)$ is the energy density, $p(r)$ is the radial pressure, and
$P(r)$ is the tangential pressure. The anisotropy is then measured by the
anisotropic pressure function: $\chi(r)\equiv{(P(r)-p(r))}/{3}$. This is often
used instead of $P(r)$. Integrating $G^{t}_{t}=8\pi T^{t}_{t}$ from the
origin, we get the “mass” function:

$m(r)=\int^{r}_{0}4\pi\rho(x)x^{2}dx,\quad\frac{dm}{dr}=4\pi\rho(r)r^{2}.$
(2.3)

One can also consider a “core” for $r\leq r_{c}\neq 0$ (which could be made of
another fluid solution, a vacuum Schwarzschild solution, or something else):

$m(r)=\int^{r}_{r_{c}}4\pi\rho(x)x^{2}dx+m_{c},\quad\frac{dm}{dr}=4\pi\rho(r)r^{2},$
(2.4)

where $m_{c}$ is the “core mass.” In either case, we get the following
“simplified” line element:

$ds^{2}=-e^{2\Phi(r)}dt^{2}+\frac{dr^{2}}{1-2m(r)/r}+r^{2}(d\theta^{2}+\sin^{2}{\theta}d\phi^{2}).$
(2.5)

For the remaining two Einstein equations, it is customary to use the
conservation equation, $\nabla_{\mu}T^{\mu}_{r}=0$, and $G^{r}_{r}=8\pi
T^{r}_{r}$ to obtain

$P=\frac{r}{2}\left(\frac{dp}{dr}+(\rho+p)\frac{d\Phi}{dr}\right)+p$ (2.6)

and

$\frac{d\Phi}{dr}-\frac{m+4\pi r^{3}p}{r(r-2m)}=0,$ (2.7)

respectively. The algebraic implication of this equation for the anisotropic
system will become important later in this paper.

Substituting (2.7) into (2.6) gives the generalized Tolman–Oppenheimer–Volkoff
(TOV) equation:

$\frac{dp}{dr}+\left[{m+4\pi\rho r^{3}\over
r(r-2m)}+\frac{2}{r}\right]p+\left(\frac{4\pi
r^{2}}{r-2m}\right)p^{2}=\frac{2P}{r}-\frac{\rho m}{r(r-2m)},$ (2.8)

which takes a slightly tidier form with $\chi(r)$:

$\frac{dp}{dr}+\left[{m+4\pi\rho r^{3}\over r(r-2m)}\right]p+\left(\frac{4\pi
r^{2}}{r-2m}\right)p^{2}=\frac{6\chi}{r}-\frac{\rho m}{r(r-2m)}.$ (2.9)

In the following, we will often refer to (2.7) and (2.8) as “the system” and
will switch between (2.8) and (2.9) freely. The perfect or isotropic system is
obtained by setting $p(r)=P(r)$ or $\chi(r)=0$ in the generalized TOV equation
and will be referred to as the “isotropic system.”

##  3 Central Regularity

The simplified line element (2.5) has reduced the number of independent
ordinary differential equations to two, equations (2.6) and (2.7) or equations
(2.7) and (2.8), in four unknowns, and thus two inputs are needed to integrate
the system (one for the isotropic case). In stellar models with isotropic
fluids, for example, the equation of state $p=p(\rho)$ acts as the functional
input needed. One takes a finite initial value $\rho(0)$ and integrates (2.4)
and (2.9) (with $\chi=0$) simultaneously until a radius $r=r_{b}$ is reached
where $p(r)$ vanishes (see, for example, [16, 24]).444One also generally
requires the solution to be surrounded by a Schwarzschild exterior solution
(i.e., $T^{\mu}_{\nu}=0\;\text{for}\;r>r_{b}$). Continuity of the metric
components requires that $e^{2\Phi(r_{b})}=\left(1-{2M}/{r_{b}}\right)$, where
$M=m(r_{b})$ is the mass of the Schwarzschild solution. Moreover, the
differentiability of the metric components requires the continuity of $p(r)$
and $\rho(r)$. However, physically, for hydrostatic equilibrium, the radial
pressure is required to be continuous across the boundary, i.e., $p(r_{b})=0$,
but $\rho(r_{b})$ can be nonzero (see e.g., [6]). Any initial finite value of
$\rho(0)$, which returns a finite central pressure via the equation of state
$p(0)=p(\rho(0))$, ensures regularity at the center.

We will be interested in generating all solutions via solution-generating
algorithms — infinite fluid spheres as well as those that terminate at finite
radii, as in stellar models. We will also pay special attention to generating
solutions that are regular at the center. As we will discuss in detail later,
in any solution generation algorithm the generating functions are the input
functions which are freely specifiable and which determine the remaining
functions via quadratures. For the resulting solutions to be regular at the
center one has to ensure that: (i) the input/generating functions have all the
properties consistent with the geometric condition of regularity and (ii)
their initial values also ensure that the initial values of the remaining
variables are consistent with regularity. If these conditions cannot be
fulfilled one will not be able to specialize to regular solutions, although
one may be able to generate _all_ solutions. With these in mind, we take a
closer look at the regularity conditions below.

As usual in the literature, by regularity, we mean the finiteness of scalar
polynomial invariants constructed from the Riemann tensor. This has been
discussed by a number of authors for both isotropic and anisotropic fluid
spheres (see, for example, [7, 21, 10, 22]). For spherically symmetric static
spacetimes it is sufficient that the Kretschmann scalar
$\kappa=R_{\mu\nu}^{\;\;\;\;\rho\sigma}R_{\rho\sigma}^{\;\;\;\;\mu\nu}$ is
finite for scalar polynomial invariants constructed from the Riemann tensor
(see, for example, [7]).

For such a spacetime in static coordinates, like (2.1) and (2.5), the
algebraically independent Riemann tensor is pairwise diagonal with the
following nonzero components:
$R_{01}^{\;\;\;01},R_{02}^{\;\;\;02}=R_{03}^{\;\;\;03},R_{12}^{\;\;\;12}=R_{13}^{\;\;\;13},R_{23}^{\;\;\;23}$.555For
a line element in stationary coordinates where $g_{01}\neq 0$, for example,
$R_{12}^{\;\;\;02}\neq 0$. Thus, the Kretschmann scalar becomes a sum of four
squares:

$\kappa=4K_{1}^{2}+8K_{2}^{2}+8K_{3}^{2}+4K_{4}^{2},$ (3.1)

where $K_{1}=-R_{01}^{\;\;\;01}$, $K_{2}=-R_{02}^{\;\;\;02}$,
$K_{3}=-R_{13}^{\;\;\;13}$ and $K_{4}=-R_{23}^{\;\;\;23}$. This is finite iff
each square term is finite.666Recall, for the Schwarzschild exterior, that
$R_{1212}=-M/(r-2M)$ i.e., singular at $r=2M$ but
$R_{12}^{\;\;\;12}=-M/r^{3}$.. Thus, the finiteness of $\kappa$ is sufficient
to ensure the regularity of these scalars [21].

In the simplified line element (2.5), one geometric and one fluid variables
appear, $\Phi(r)$ and $\rho(r)$, and777The tensor calculations were performed
with Maple using the GRTensor III package.

$\kappa=4\left[\left(1-{2m\over r}\right)(\Phi^{\prime\prime}+\Phi^{\prime
2})+\Phi^{\prime}\left(\frac{m}{r^{2}}-4\pi\rho r\right)\right]^{2}\\\
+8\left[\left(1-\frac{2m}{r}\right)\frac{\Phi^{\prime}}{r}\right]^{2}+8\left({m\over
r^{3}}-4\pi\rho\right)^{2}+4\left(\frac{2m}{r^{3}}\right)^{2}.$ (3.2)

As $r\rightarrow 0$, ${m}/{r^{3}}\rightarrow 4\pi\rho/3$ from the definition
of $m(r)$, equation (2.4). The last two terms in $\kappa$ near $r=0$ are
separately proportional to $\rho$, so one needs $\rho(0)$ to be finite. For
the second term to be finite at the center, one needs $\Phi^{\prime}(0)=0$ and
$\Phi^{\prime\prime}(0)$ finite simultaneously, which make the first term
finite without any extra conditions. Thus, the necessary and sufficient
conditions that $\kappa$ is finite at $r=0$ are

$\displaystyle\rho(0)\text{ finite},$ (3.3)
$\displaystyle\Phi^{\prime}(0)=0\text{ and }\Phi^{\prime\prime}(0)\text{
finite}.$ (3.4)

Under the above conditions, the first and second Weyl invariants, given by
$(C(r))^{2}/6$ and $(C(r))^{3}/36$, respectively, where

$\displaystyle
C(r)=\left(1-\frac{2m}{r}\right)(\Phi^{\prime\prime}+\Phi^{\prime
2})+4\pi\rho-\frac{3m}{r^{3}}-\frac{\Phi^{\prime}}{r}\left(4\pi\rho
r^{2}-{3m\over r}+1\right),$ (3.5)

vanish identically as $r\rightarrow 0$. This shows that the metric is
conformally flat at the origin under these conditions [10, 21]. These
conditions are geometric in that none of the Einstein equations except the
$G^{t}_{t}$ component has been used, which amounts to merely renaming
$\Psi(r)$. These apply to both isotropic and anisotropic fluids spheres.

#### Equivalent Conditions of Regularity

We now explore the consequence of these geometric conditions on the other
variables of the system. Applying them to the $G^{r}_{r}$ and
$G^{\theta}_{\theta}$/$G^{\phi}_{\phi}$ equations,888 $G^{r}_{r}\equiv
2\left(1-\frac{2m}{r}\right)\frac{\Phi^{\prime}(r)}{r}-\frac{2m}{r^{3}}=8\pi
p(r)$ (3.6) and
$G^{\theta}_{\theta}=G^{\phi}_{\phi}=\left(1-\frac{2m}{r}\right)\left(\Phi^{\prime\prime}+\frac{\Phi^{\prime}}{r}+\Phi^{\prime
2}\right)+\left(1+r\Phi^{\prime}\right)\left(\frac{m}{r^{3}}-4\pi\rho\right)=8\pi
P(r).$ (3.7) we find

$\displaystyle p(0)=\frac{1}{4\pi}\Phi^{\prime\prime}(0)-\frac{1}{3}\rho(0),$
(3.8) $\displaystyle
P(0)=\frac{1}{4\pi}\Phi^{\prime\prime}(0)-\frac{1}{3}\rho(0).$ (3.9)

The first relation between nonzero finite quantities is known [20] and they
will be useful for our analysis. In terms of anisotropy one can write them as

$\displaystyle p(0)=\frac{1}{4\pi}\Phi^{\prime\prime}(0)-\frac{1}{3}\rho(0),$
(3.10) $\displaystyle\chi(0)=0.$ (3.11)

Note that any two of $\rho(0)$, $p(0)$, and $\Phi^{\prime\prime}(0)$ will fix
the initial value of the third quantity via (3.8) as well as imply that
$\Phi^{\prime}(0)=0$ and $P(0)=p(0)$ (or $\chi(0)=0$) as a consequence of the
Einstein equations. In the isotropic case, $P(0)=p(0)$ holds trivially and
finiteness of $\rho(0)$ ensures finiteness of $p(0)$ via an equation of state
as in stellar models. Thus, the following combinations are equivalent for
regularity:

$\displaystyle\rho(0)\,\,\,\mathrm{finite},\,\,\,\Phi^{\prime}(0)=0,\,\,\,\Phi^{\prime\prime}(0)\,\,\,\mathrm{finite}$
(3.14) $\displaystyle
p(0)\,\,\,\mathrm{finite},\,\,\,\Phi^{\prime}(0)=0,\,\,\,\Phi^{\prime\prime}(0)\,\,\,\mathrm{finite}$
$\displaystyle\rho(0)\,\,\,\mathrm{finite},\,\,\,p(0)\,\,\,\mathrm{finite}$

Note that none of the above uses $P(0)$. Combinations of $P(0)$ with other
functions like

$\displaystyle
P(0)\,\,\,\mathrm{finite},\,\,\,\Phi^{\prime}(0)=0,\,\,\,\Phi^{\prime\prime}(0)\,\,\,\mathrm{finite}$
(3.15)
$\displaystyle\rho(0)\,\,\,\mathrm{finite},\,\,\,P(0)\,\,\,\mathrm{finite}$
(3.16)

are not sufficient to reproduce the geometric conditions of regularity. One,
thus, needs a supplemental condition to be able to use them. As we will see
when generating regular solutions using $(\Phi^{\prime}(r),P(r))$, such a
condition arises quite fortuitously from the considerations of the
differential equation. We will also see that

$\displaystyle\chi(0)=0,\,\,\,\Phi^{\prime}(0)=0,\,\,\,\Phi^{\prime\prime}(0)\,\,\,\mathrm{finite}$
(3.17)

and

$\displaystyle\Phi^{\prime}(0)=0,\,\,\,\Phi^{\prime\prime}(0)\,\,\,\mathrm{finite}$
(3.18)

can serve as sufficient conditions for regularity for the anisotropic and
isotropic systems, respectively. This is because Einstein’s equations will
produce the right behavior for $m(r)$, which ensures that $\kappa$ is finite
at the origin.

#### TOV and Momentum Conservation with Regularity

The only quantity that has not appeared so far in the discussion is
$p^{\prime}(r)$. To see the effect of regularity on $p^{\prime}(r)$, it is
natural to turn to the generalized TOV equation.999The TOV equation can be
turned into a differential equation in $m(r)$ and $\Phi(r)$, as we will see
later. For the anisotropic case, it is easy to see that the TOV applies only
for $r>0$. However, we argue that, even for the isotropic case, one needs to
look at the original momentum conservation equation (2.6) where
$p^{\prime}(r)$ appears only in $rp^{\prime}(r)$, which is clearly valid at
$r=0$. The conditions of regularity then imply that $p^{\prime}(r)$ will
satisfy the following condition:

$\lim_{r\to 0}rp^{\prime}(r)=0\,\,\,\,(\text{isotropy and anisotropy})$ (3.19)

under the above regularity conditions. This condition can be fulfilled by any
finite $p^{\prime}(0)$ (including zero), which would have been the only
conclusion had we used the TOV equation. However, (3.19) clearly shows $p(r)$
does not have to be differentiable at $r=0$ while $p(0)$ is finite. Such
possibilities should be included if one is interested in generating all
regular solutions.

Further corroboration of condition (3.19) comes from casting the Kretschmann
scalar in terms of $p(r)$ using (2.7):

$\kappa=4\left(4\pi(\rho+p)-\frac{2m}{r^{3}}+rp^{\prime}(r)+\frac{(m+4\pi
r^{3}p)(p+\rho)}{r-2m}\right)^{2}\\\ +8\left({m\over r^{3}}+4\pi
p\right)^{2}+8\left({m\over
r^{3}}-4\pi\rho\right)^{2}+4\left(\frac{2m}{r^{3}}\right)^{2}.$ (3.20)

The last two terms are the same as before and the second term implies that
$p(0)$ should be finite and that $\lim_{r\to 0}rp^{\prime}(r)$ is _finite_.
However, for conformal flatness at the center, the vanishing of

$\displaystyle C(r)=4\pi\left(rp^{\prime}(r)+\frac{(m+4\pi
r^{3}p)(p+\rho)}{r-2m}+2\rho-\frac{6m}{4\pi r^{3}}\right)$ (3.21)

requires $\lim_{r\to 0}rp^{\prime}(r)=0$. Also note that, in either of the
expressions above, $p^{\prime}(r)$ does not appear except in $rp^{\prime}(r)$.
That the finiteness of $\kappa$ in terms of $\rho(r)$ and $p(r)$ does not
guarantee conformal flatness at the center, unlike when $\rho(r)$ and
$\Phi^{\prime}(r)$ were used, shows that something can get lost in
translation.

To summarize, in choosing $p(r)$ as a generating function, it has to be chosen
such that

$\displaystyle p(0)\,\,\,\mathrm{finite},\,\,\,\lim_{r\to 0}rp^{\prime}(r)=0$
(3.22)

to ensure all possible regular solutions. Other variables in (3.14)–(3.16) do
not need any such qualification except that they are continuous and
differentiable as shown. It is easy to see that, given the three quantities in
the algebraic equation (3.8), no one function alone can ensure regularity at
the center. This makes regularity in the anisotropic systems, where two inputs
are needed, much easier to handle.

_A Simple Nontrivial Example:_ We have discussed above how it is natural to
conclude that $p^{\prime}(0)$ is finite if one uses the TOV equation. This
narrative is reinforced by the fact that $p^{\prime}(0)=0$ is often preferred
physically. However, it is easy to see that $\lim_{r\to 0}rp^{\prime}(r)=0$
(with finite $p(0)$) can be satisfied by functions like the following:

$p(r)=c-r^{1/n}e^{r},$ (3.23)

where $n>1$ and $c>0$ is a constant. It is monotonically decreasing and
$\lim_{r\to 0^{+}}p^{\prime}(r)=\mathrm{DNE}$, but $\lim_{r\to
0^{+}}rp^{\prime}(r)=0$ and $p(0)=c$. It is easy to construct similar
examples. Such solutions should be included in the consideration of all
regular solutions for both the isotropic and anisotropic systems.

##  4 Generating Functions and Solutions

Although for any pair of input functions one can find a solution of the system
(2.7) and (2.8), for a pair of generating functions the system admits a
general form solution and the metric can be expressed formally in terms of
integrals involving these functions only. One then has what is called a
generation algorithm.101010Whether those integrations can be performed exactly
to obtain an exact solution is considered an _a posteriori_ issue and
considered individually. In most cases, this will require numerical
integration, as mentioned in the introduction.

As mentioned in the introduction, generating functions do not have to be from
the set of basic variables, although this is the case in the existing
algorithms in canonical coordinates, which we will revisit below [20, 18]. We
will consider all six combinations of two inputs that are possible out of
$\rho(r)$, $p(r)$, $P(r)$, and $\Phi(r)$. First, note that the following sets
are equivalent:

  1. 1.

Linear combinations of $\rho(r)$, $p(r)$, and $P(r)$ are equivalent. This just
involves rewriting the Einstein equations using corresponding linear
combinations. This makes $P(r)$ and $\chi(r)$ equivalent inputs when used
alongside $p(r)$, which we have seen in (2.8) and (2.9).111111If $p(r)$ is not
used, one may be preferred over the other, as we have seen in the regularity
conditions (3.16). However, if one wants to take advantage of the simplified
line element (2.5), one must keep $\rho(r)$ as it is, as we do.

  2. 2.

Specifying a function is equivalent to specifying another function if they are
in a 1-1 algebraic correspondence. This obvious fact will help identify
equivalent functions with the basic functions of the system.

  3. 3.

For complete spheres, $r\geq 0$, $\rho(r)$ and $m(r)$ are in 1-1
correspondence thanks to (2.4). Either is equivalent to specifying $e^{\pm
2\Psi(r)}$. In the presence of a core, $\rho(r)$ determines $m(r)$ up to an
arbitrary constant, the unspecified mass of the core. This point will be
elaborated in Sections 4.2 as we will revisit the existing algorithms.

  4. 4.

In both (2.6) and (2.7), only $\Phi^{\prime}(r)$ appears. The additive
constant resulting from the integration of $\Phi(r)$ becomes a multiplicative
constant in the line element and, as such, can be absorbed into the time
coordinate of (2.5). Thus, we can treat $\Phi^{\prime}(r)$ as the basic
function instead of $\Phi(r)$.

###  4.1 Algorithm for Generating All Isotropic Solutions

The isotropic system needs one input to generate a solution. This could be an
equation of state, one of $\Phi(r)$, $\rho(r)$, and $p(r)$, or any suitable
combination of these. Exact solutions of this system have been studied
extensively and were mostly obtained using _ad hoc_ assumptions on $\rho(r)$,
$e^{2\Phi(r)}$ (equivalently, $\Phi^{\prime}(r)$), and $e^{-2\Psi(r)}$
(equivalently, $m(r)$). Out of the three basic functions, only one can serve
as a generating function, as we will see below:

  * •

$\rho(r)$: If one uses $\rho(r)$ as the input, the TOV equation (2.8) is a
(first-order) Riccati equation in $p(r)$ of the form
$p^{\prime}+f(r)p+g(r)p^{2}=h(r)$. By replacing $p(r)$ in (2.8) and using
(2.6), on the other hand one gets a Riccati equation in $\Phi^{\prime}(r)$. In
either case, the basic problem is that the general Riccati is not solvable by
quadratures. Thus, one cannot create an algorithm to generate all solutions
using $\rho(r)$ as the input.

  * •

$p(r)$: With $p(r)$ as the input, the TOV equation is an Abel equation of the
second kind in $m(r)$ (see, for example, [32]). This is a more difficult
equation to solve than the Riccati equation, even with special choices of
coefficients, which perhaps explains why there is not a single example of
$p(r)$ being used as the input to obtain an exact solution for the isotropic
system [27].

  * •

$\Phi^{\prime}(r)$: With $\Phi^{\prime}(r)$ as input, the TOV equation turns
into a first-order linear differential equation in $m(r)$, which, as Weynman
noticed in 1949, “can always be solved by quadratures” [31] and thus, one can
have an algorithm that generates all solutions [20]. The details of this
algorithm are contained in the corresponding anisotropic algorithm discussed
below. This has, thus, been a popular route for obtaining an exact solution,
too [27]. However, it is not easy to know beforehand what choice of $\Phi(r)$
would produce the desired properties of $\rho(r)$ and $p(r)$ or a particular
equation of state. To generate solutions that are regular at the center,
$\Phi(r)$ must have the property that $\Phi^{\prime}(0)=0$ and
$\Phi^{\prime\prime}(0)$ must be finite, as a necessary condition, as we
discussed above in Section 3 and in [20]. We will address this below and see
how Einstein’s equations will fill the gap for regularity.

###  4.2 Generating Regular Anisotropic Solutions from the Existing Algorithm

For the two simultaneous inputs needed for the anisotropic system, if one
combines anisotropy $\chi(r)$ (or, equivalently, $P(r)$) with $\rho(r)$,
$p(r)$, or $\Phi(r)$, the qualitative behaviors remain the same as above. For
$(\rho(r),P(r))$ and $(\Phi(r),P(r))$, the generalized TOV is still a Riccati
equation and an Abel equation of the second kind in $m(r)$, respectively, with
no gain in their general solvability. The only generating pair is
$(\Phi(r),P(r))$, which a natural extension of the isotropic case, and has
been worked out using $\chi(r)$ [18]. Below we will rederive this and look for
specialized algorithms that generate all regular solutions.

#### Existing Algorithm for (An)isotropic Solutions: A Second Look

The generalized TOV can be turned into a first-order differential equation in
$m(r)$ for either $(\Phi(r),P(r))$ or $(\Phi(r),\chi(r))$ as input. One can
solve (2.7) algebraically for $p(r)$ and substitute it in (2.8) or (2.9) with
$m^{\prime}(r)=4\pi\rho r^{2}$ to get121212One can also use
$G^{\theta}_{\theta}$ and $G^{r}_{r}-G^{\theta}_{\theta}$ equations for this.:

$\frac{dm(r)}{dr}+a(r)m(r)=b(r).$ (4.1)

Unlike the Riccati or the Abel equation, this can be solved exactly for
arbitrary forms of its variable coefficients:

$m(r)=\frac{\int b(r)e^{\int a(r)dr}dr+C}{e^{\int a(r)dr}}.$ (4.2)

The coefficients for $(\Phi(r),P(r))$ are

$\displaystyle a(r)$ $\displaystyle=\left(\frac{2r^{2}(\Phi^{\prime
2}+\Phi^{\prime\prime})+r\Phi^{\prime}-1}{r(r\Phi^{\prime}+1)}\right),$ (4.3)
$\displaystyle b(r)$
$\displaystyle=\frac{r(\Phi^{\prime}+r\Phi^{\prime\prime}+r\Phi^{\prime
2}-8\pi rP)}{(r\Phi^{\prime}+1)}.$ (4.4)

Thus, there is a one-parameter family of solutions $m(r)$ for any input
$(\Phi(r),P(r)$), which will give a parameter-dependent $p(r)$ algebraically
via (2.7). The resulting one-parameter family of geometries found by
substituting (4.2) and $\Phi(r)$ in the line element (2.5) will differ
nontrivially in its $g_{rr}$ component for different values of $C$. The
isotropic case is recovered by setting $P(r)=p(r)$; in this case, one would
have eliminated at the previous step the $P(r)$ appearing in $b(r)$. If one
uses $\chi(r)$ instead of $P(r)$, the coefficients of (4.1) are then as
follows:

$\displaystyle a(r)$
$\displaystyle=\left(\frac{2r(\Phi^{\prime\prime}+\Phi^{\prime
2})}{(r\Phi^{\prime}+1)}-\frac{3}{r}\right),$ (4.5) $\displaystyle b(r)$
$\displaystyle=\frac{r(r\Phi^{\prime\prime}+r\Phi^{\prime
2}-\Phi^{\prime}-24\pi r\chi)}{(r\Phi^{\prime}+1)}.$ (4.6)

The algorithm presented in [18] uses the inputs $(z(r),\Pi(r))$:

$\displaystyle z(r)$ $\displaystyle=\Phi^{\prime}+1/r,$ (4.7)
$\displaystyle\Pi$ $\displaystyle=8\pi(p-P),$ (4.8)

which, from the discussion at the beginning of this section, are equivalent to
the $(\Phi^{\prime}(r),P(r))$ or $(\Phi^{\prime}(r),\chi(r))$ that we used
above. It also used $y(r)=1-2m(r)/r$ instead of $m(r)$. This is
$e^{-2\Psi(r)}$. With these, it is easy to check that (4.1) is precisely
equation (8) of [18]:

$y^{\prime}+y\left(\frac{2z^{\prime}}{z}+2z-\frac{6}{r}+\frac{4}{zr^{2}}\right)=-\frac{2}{z}\left(\Pi+\frac{1}{r^{2}}\right).$
(4.9)

We do not reproduce the final one-parameter $y(r)$ or the metric expressed in
terms of $z$ and $\Pi$ [18]. Although completely equivalent, using $m(r)$
instead of $y(r)$ makes it slightly easier to interpret the free parameter in
the general solution.

#### Regular Solutions using $\Phi(r)$ and $P(r)$

There is no discussion in the literature on the role of the parameter that
appears in the solution of $y(r)$ in [18] or in $m(r)$ in [20] for the
isotropic case. We will now identify the parameter by looking at (4.1) from
the point of view of an initial-value problem first and discuss its subsequent
role as $r\to 0$. Assuming, there is a core, i.e., $r=r_{c}>0$, one needs to
solve the following initial-value problem for $r\geq r_{c}$:

$\frac{dm(r)}{dr}+a(r)m(r)=b(r),\quad m(r_{c})=m_{c},$ (4.10)

where $m_{c}$ is the mass of the core. This has the unique solution:

$m(r)=\frac{\int_{r{{}_{c}}}^{r}b(x)e^{\int_{r{{}_{c}}}^{x}a(s)ds}dx+m_{c}}{e^{\int_{r{{}_{c}}}^{r}a(s)ds}},$
(4.11)

which, on comparison with (4.2), shows that $C=m_{c}$ precisely. Thus, the
parameter $C$ in our formulation above, and in that of [20], represents the
mass of the (nonzero) core.131313Note that this does not follow from the look
of (4.1) since $C$ is in its numerator. (We will soon see that for complete
spheres it will assume a different meaning.) One can further qualify $C$ on
physical grounds for nonzero cores. For example, for a Schwarzschild core, one
must have:

$m_{c}\leq\frac{r_{c}}{2}.$ (4.12)

Thus, contrary to the impression given by the general solution (4.2) of
$m(r)$, not all possible values of the parameter $C$ are admissible.

For complete fluid spheres with $r\geq 0$, unfortunately, one cannot set
$m(0)=0$ to get $C=0$ as $r=0$ is a singular point of the differential
equation ($\lim_{r\to 0^{+}}a(r)$ diverges). However, $r=0$ is a regular
singular point since $\lim_{r\to 0^{+}}ra(r)$ is finite (see, for example,
[9]). To analyze the solution near it, we consider the dominant terms in
$a(r)$ and $b(r)$ as $r\to 0$:

$\displaystyle a(r)$ $\displaystyle=-\frac{1}{r},$ (4.13) $\displaystyle b(r)$
$\displaystyle=r^{2}\left[2\Phi^{\prime\prime}(0)-8\pi P(0)\right],$ (4.14)

where we have used $\lim_{r\to 0^{+}}\Phi(r)/r=\Phi^{\prime\prime}(0)$. This
yields the following general solution:

$m(r)=\left[\Phi^{\prime\prime}(0)-4\pi P(0)\right]r^{3}+C_{1}r.$ (4.15)

Imposing the initial value $m(0)=0$ will not allow us to determine the
arbitrary constant $C_{1}$, illustrating the breakdown of the uniqueness
theorem at the singular point $r=0$. This linear term will prevent $\lim_{r\to
0^{+}}m(r)/r^{3}$ from being finite and will give a non-regular solution for
any $C_{1}\neq 0$.141414The condition $m(0)=0$ is a necessary but not
sufficient condition of regularity; it avoids a conical singularity [30] and
one needs $m(r)/r^{3}$ to be finite as we discussed earlier. This is expected
from our discussions on regularity in Section 3 since $\Phi^{\prime}(0)=0$,
$\Phi^{\prime\prime}(0)$ finite, and $P(0)$ finite are not sufficient for
regularity.

The cubic term in (4.15) comes from the first term of (4.2). The linear term,
on the other hand, comes from the second term of (4.2), $Ce^{-\int a(r)dr}$,
which gives the difference between any two solutions of (4.1). It is, thus,
easy to see that $C_{1}=C$ in (4.15). So, if we can find a condition that
requires $C_{1}=0$ in (4.15), we can set $C=0$ in general.

Fortunately, this can be argued from the general solution (4.2) as a limiting
process as follows. After integrating the diverging term $1/r$, (4.2) can be
rewritten as

$\frac{m}{r^{3}}=\frac{C}{r^{2}}\,e^{-A(r)}+\frac{e^{-A(r)}}{r^{2}}\left(\int\frac{b(r)}{r}e^{A(r)}dr\right)$
(4.16)

where

$A(r)=\int\alpha(r)dr=2\int\frac{\Phi^{\prime}+r(\Phi^{\prime\prime}+\Phi^{\prime
2})}{(r\Phi^{\prime}+1)}dr.$ (4.17)

It is easy to see as $r\to 0$, both $\alpha(r)$ and $b(r)/r$ go to zero. As
$r\to 0$ the first term of the right hand side diverges since the integrand
$\alpha(r)\to 0$. However, the other part is finite:

$\lim_{r\to
0}\left(\int\frac{b(r)}{r}e^{A(r)}dr\right)/r^{2}=\Phi^{\prime\prime}(0)-4\pi
P(0).$ (4.18)

Since $\lim_{r\to 0^{+}}m/r^{3}=4\pi\rho(0)/3$, one get from (4.16)

$\frac{4\pi}{3}\rho(0)=\Phi^{\prime\prime}(0)-4\pi P(0).$ (4.19)

This is the regularity condition (3.9). We, thus, have the following theorem.

Theorem 4.1: For a differentiable function $\Phi^{\prime}(r)$ and continuous
function $P(r)$ defined on the common interval $[0,r_{b})$, $0<r_{b}<\infty$,
such that $\Phi^{\prime}(0)=0$ and $\Phi^{\prime\prime}(0)$ is finite, and
$P(0)$ is finite, the input $(\Phi^{\prime}(r),P(r)$) generates a unique
anisotropic fluid sphere that is regular at the center with mass given by

$m(r)=re^{-\int_{0}^{r}\alpha(s)ds}\int_{0}^{r}\frac{b(x)}{x}e^{\int_{0}^{x}\alpha(s)ds}dx,$
(4.20)

and central density given by

$\rho(0)=\frac{3}{4\pi}\Phi^{\prime\prime}(0)-3P(0).$

#### Regular Solutions using $\Phi(r)$ and $\chi(r)$

The above (i.e., the vanishing of $C_{1}$ and, hence, $C$) can also be argued
by invoking the vanishing of anisotropy at the center in regular solutions,
i.e., $\chi(0)=0$. Using anisotropy, i.e., (4.5) and (4.6) for $a(r)$ and
$b(r)$, and keeping the dominant terms as before one gets:

$\displaystyle a(r)$ $\displaystyle=-\frac{3}{r},$ (4.21) $\displaystyle b(r)$
$\displaystyle=r^{4}\Phi^{\prime\prime}(0),$ (4.22)

which gives

$m(r)=Cr^{3}+\frac{1}{2}\Phi^{\prime\prime}(0)r^{5}.$ (4.23)

Here, opposite to (4.15), the first and second terms come, respectively, from
the second and first terms of the corresponding general solution
(4.2).151515One would get an identical result using (4.9). Other terms in the
expansions of $\Phi^{\prime}(r)$ and $\chi(r)$ would produce terms with powers
higher than three, as one can easily check.

Comparing (4.15) with (4.23) clearly shows that we have to set $C_{1}=0$ in
(4.15), and hence, $C=0$ in (4.2), to achieve isotropy and regularity
simultaneously at the center. This also shows that in (4.23)
$C=\Phi^{\prime\prime}(0)-4\pi P(0)$, which could not be achieved using only
anisotropy. The resulting $m(r)$ and $p(r)$ will be parameter-free and so the
metric (2.5) will be unique for every appropriate choice of
$(\Phi^{\prime}(r),P(r)$).

The general solution (4.2) with $\chi(r)$ can also be used for $r=0$ provided
one imposes the conditions of regularity as $r\to 0^{+}$. This can be seen by
rewriting the general solution as

$\frac{m}{r^{3}}=e^{-A(r)}C+e^{-A(r)}\left(\int\frac{b(r)}{r^{3}}e^{A(r)}dr\right)$
(4.24)

where

$A(r)=\int\alpha(r)dr=\int\frac{2r(\Phi^{\prime\prime}+\Phi^{\prime
2})}{(r\Phi^{\prime}+1)}dr.$ (4.25)

It is easy to check that $\lim_{r\to 0^{+}}\alpha(r)=0$ and $\lim_{r\to
0^{+}}b(r)/r^{3}$ is finite if, for $\chi(0)=0$, $\chi^{\prime}(0)$ and
$\Phi^{\prime\prime\prime}(0)$ are finite. However, these extra conditions do
not appear in the Einstein equations for the system. So, the right way to
recover the solution at $r=0$ from the general form (4.24) above (which was
obtained for $r\neq 0$) is, first, to impose the conditions of regularity:
$\chi(0)=0$ and $\lim_{r\to 0^{+}}\Phi(r)/r=\Phi^{\prime\prime}(0)$. After
integration, the second term in (4.24) will then precisely reproduce the
$\Phi^{\prime\prime}(0)r^{5}/2$ term of (4.23). With this, (4.24) is a regular
one-parameter family of solutions.161616For the isotropic system, this freedom
with the parameter is what was used in Theorem (P2) in [5] to generate new
solutions (with $C=4\pi(\delta\rho)/3$ where $\delta\rho$ is shift in central
density).

Theorem 4.2: For two differentiable functions $\Phi^{\prime}(r)$ and $\chi(r)$
defined on the common interval $[0,r_{b})$, $0<r_{b}<\infty$, such that
$\Phi^{\prime}(0)=\chi(0)=0$ and $\Phi^{\prime\prime}(0)$ is finite, each
input $(\Phi^{\prime}(r),\chi(r)$) generates a one-parameter family of
anisotropic fluid spheres that are regular at the center.

Thus, the conditions $\Phi^{\prime}(0)=\chi(0)=0$ and finite
$\Phi^{\prime\prime}(0)$ are sufficient conditions for regularity for
anisotropic fluid spheres. Note that these conditions alone do not make
$\kappa$ finite. However, they prove sufficient for the system because the
Einstein equations lead to $m(r)$ having the right behavior under these
conditions, making $\kappa$ finite. We, thus, have the following for the
isotropic system:

Corollary 4.3: For a differentiable function $\Phi^{\prime}(r)$ defined on the
interval $[0,r_{b})$, $0<r_{b}<\infty$, such that $\Phi^{\prime}(0)=0$ and
$\Phi^{\prime\prime}(0)$ is finite, the input function $\Phi^{\prime}(r)$
generates a one-parameter family of isotropic fluid spheres that are regular
at the center.

###  4.3 Three New Choices and an Algorithm

We now consider the remaining three combinations of the four variables of the
anisotropic system: $(\rho(r)$, $p(r)$), $(\rho(r),\Phi(r))$, and
$(p(r),\Phi(r))$. The line element (2.5) can be formally integrated and be
expressed it in terms of $\rho(r)$ and $p(r)$171717This is obtained by
rewriting $\Phi(r)$:  $\displaystyle\frac{d\Phi}{dr}$
$\displaystyle=\frac{m(r)+4\pi r^{3}p}{r(r-2m(r))}$ (4.26)
$\displaystyle=\frac{m(r)}{r(r-2m(r))}+\frac{4\pi r^{3}p(r)}{r(r-2m(r))}$
(4.27)
$\displaystyle=\frac{1}{2}\left[\frac{1}{r-2m(r)}-\frac{1}{r}\right]+\frac{4\pi
r^{2}p(r)}{r-2m(r)}$ (4.28)
$\displaystyle=\frac{1}{2}\left[\frac{1-2m^{\prime}(r)}{r-2m(r)}-\frac{1}{r}\right]+\frac{4\pi
r^{2}p(r)}{r-2m(r)}+\frac{m^{\prime}(r)}{r-2m(r)}.$ (4.29) :

$ds^{2}=-\left(1-\frac{2m(r)}{r}\right){e^{8\pi\int\frac{\left(\rho(r)+p(r)\right)r^{2}}{r-2m(r)}\,dr}}dt^{2}+\frac{dr^{2}}{1-\frac{2m(r)}{r}}+r^{2}(d\theta^{2}+\sin^{2}{\theta}d\phi^{2}).$
(4.30)

This metric requires two integrations to generate any particular solution, but
it is now abundantly clear that one can generate _all_ solutions by specifying
$(\rho(r),p(r))$. However, what about the other two pairs and regular
solutions? The answer can be summarized as follows.

Lemma 4.4: For the anisotropic system of fluids in canonical coordinates, any
pair of the input functions $(\rho(r)$, $p(r)$), $(\rho(r),\Phi(r))$, or
$(p(r),\Phi(r))$:

  1. (a)

are equivalent,

  2. (b)

generate all solutions,

  3. (c)

determine the geometry (2.5) uniquely upon two integrations,

  4. (d)

determine $P(r)$ uniquely.

Proof: We first note that the algebraic relation (2.6) is linear in
$\Phi^{\prime}(r)$, $m(r)$, and $p(r)$. Thus, any two determine the third one
uniquely.181818Explicitly:  $\displaystyle\Phi^{\prime}(r)$
$\displaystyle=\frac{m(r)+4\pi r^{3}p}{r(r-2m(r))},$ (4.31) $\displaystyle
p(r)$ $\displaystyle=\frac{r(r-2m(r))\Phi^{\prime}(r)-m(r)}{4\pi r^{3}},$
(4.32) $\displaystyle m(r)$ $\displaystyle=\frac{r^{2}\Phi^{\prime}(r)-4\pi
r^{3}p(r)}{1+2r\Phi^{\prime}(r)}.$ (4.33) We have already noted that
$\Phi^{\prime}(r)$ and $\Phi(r)$ are equivalent and that $\rho(r)$ and $m(r)$
are equivalent if we have the origin at $r=0$ (with zero cental mass) or a
core $r\leq r_{c}$ with a known mass. Thus, $(m(r)$, $p(r)$),
$(\rho(r),\Phi(r))$, $(p(r),\Phi(r))$, and $(\rho(r)$, $p(r)$) are equivalent
inputs. In the absence of any such boundary condition, purely at the
differential equation level, specifying $\rho(r)$ will return a family of
functions $m(r)$ with an arbitrary additive constant (parameter) $C$. The
equivalence will continue to hold, since ($\rho(r),p(r))$ and
($\rho(r),\Phi(r))$ will then, respectively, determine $\Phi(r)$ and $p(r)$
via the same algebraic relation above, up to the same parameter. On the other
hand, using $(p(r),\Phi(r))$ as input will determine $m(r)$ and hence,
$\rho(r)$ uniquely.

Given that $(\rho(r)$, $p(r)$), $(\rho(r),\Phi(r))$, and $(p(r),\Phi(r))$ are
equivalent inputs, it is sufficient to work with the pair $(\rho(r)$, $p(r)$)
to see that it generates all solutions. $\rho(r)$ determines $m(r)$ and hence,
$g_{rr}$ and the pair $(\rho(r)$, $p(r)$) fixes $\Phi(r)$ and hence, $g_{tt}$
uniquely and determines $P(r)$ uniquely algebraically via (2.8). Thus, for any
of $(\rho(r)$, $p(r)$), $(\rho(r),\Phi(r))$, or $(p(r),\Phi(r))$, one requires
two integrations to find $m(r)$ and $\Phi(r)$. For fluid solutions starting
from $r=0$ or surrounding a central core of known mass, every choice of input
pair will generate a unique sphere. On the other hand, if one does not
consider a boundary value imposed by the center, as we discussed above, one
gets a one-parameter family.

It is important to appreciate that nothing trivial is happening here. The line
element (4.30) is equally valid for the isotropic case, but there $\rho(r)$
and $p(r)$ cannot be specified simultaneously. For the same reason, the
algebraic equivalence of the three pairs loses its significance (except at the
level of regularity) in the isotropic case, since one can specify only one
input. On the other hand, specifying $\rho(r)$ or $p(r)$ leads to Abel or
Riccati equations, which are very difficult to solve in general. Here, we can
choose $\rho(r)$ and $p(r)$ simultaneously and independently, which makes it
immensely simpler to satisfy the regularity conditions. This is a clear
advantage of the anisotropic system over the isotropic system.

The isotropic limit of any anisotropic solution with the same potential
function can be found by applying the algorithm of [20]. For a given potential
function $\Phi(r)$, if $m_{A}(r)$ and $m_{I}(r)$ are the mass functions for
the anisotropic and isotropic solutions, respectively, then by introducing a
parameter $\alpha$, one can write:

$m_{A}(r)=m_{I}(r)+\alpha F(r).$ (4.34)

Since the field equations (3.6) and (3.7) are linear (as an algebraic equation
and as a differential equation) in the mass function, the radial and
tangential pressures have linear relations with the isotropic pressure. (See
[8] or substitute (4.34) into (3.6) and (3.7) to obtain $p(r)$ and $P(r)$ in
terms of the isotropic pressure and $F(r)$.)

#### Regular Solutions via the New Algorithm

Combining with the discussions on regularity in Section 3, it is very clear
how one can generate all regular anisotropic solutions using the above
algorithm in addition to Theorem 4.1

Theorem 4.5: All anisotropic fluid solutions of Einstein equations in
canonical coordinates that are regular at the center can be generated using
any of the following:

  1. (a)

$(\rho(r)$, $p(r)$) with $\rho(0)$ finite, $p(0)$ finite, and $\lim_{r\to
0}rp^{\prime}(r)=0$,

  2. (b)

$(\rho(r),\Phi(r))$ with $\rho(0)$ finite, $\Phi^{\prime}(0)=0$, and
$\Phi^{\prime\prime}(0)$ finite,

  3. (c)

$(p(r),\Phi(r))$ with $p(0)$ finite, $\lim_{r\to 0}rp^{\prime}(r)=0$,
$\Phi^{\prime}(0)=0$, and $\Phi^{\prime\prime}(0)$ finite.

This algorithm cannot be specialized to produce only isotropic solutions.
However, it does generate all isotropic solutions as the special case when the
pair $(\rho(r),p(r))$ satisfies the TOV equation of the isotropic system.

##  5 Harvesting the Riccati

An interesting feature of a Riccati equation is that if a particular solution
is known, one can find the general solution by adding to the particular
solution the general solution of an associated Bernoulli equation.191919There
is a nice discussion in the appendix of [5] on the possible forms of general
solutions of the Riccati equation when one, two, or three particular solutions
are known; in all cases there is one single constant of integration, as
expected. There are a number of transformations that convert a general Riccati
into a homogeneous linear second-order ordinary differential equations whose
general solutions have two arbitrary constants; however, all of these map back
to a one-parameter solution of the Riccati equation. Finding a particular
solution, thus, is as difficult as finding the general solution, unless the
particular solution is exported from elsewhere. The two algorithms can do just
that.

### Riccati in $p(r)$

Any “seed” solution $(\rho_{i}(r)$,
$p_{i}(r),\Phi^{\prime}_{i}(r),\chi_{i}(r))$ or, equivalently, $(\rho_{i}(r)$,
$p_{i}(r),\Phi^{\prime}_{i}(r),P_{i}(r))$ generated by either algorithm can be
seen as a particular solution of equation (2.8) and (2.8). They both are
Riccati equations in $p(r)$:

$\frac{dp}{dr}+f(r)p+g(r)p^{2}=h(r),$ (5.1)

where the coefficients $f(r)$, $g(r)$, and $h(r)$ are determined by
$\rho_{i}(r)$ and $\chi_{i}(r)$ as in (2.9) or by $\rho_{i}(r)$ and $P_{i}(r)$
as in (2.8). The general solution of (2.9) is

$p_{N}(r)=p_{i}(r)+\frac{Ce^{-\int\left(2g(r)p_{i}(r)+f(r)\right)dr}}{1+C\int
g(r)e^{-\int\left(2g(r)p_{i}(r)+f(r)\right)dr}dr}.$ (5.2)

It is easy to check that as $r\rightarrow 0$, both $f(r)\rightarrow 0$ and
$g(r)\rightarrow 0$ under the same conditions of regularity discussed in
Section 3. Thus,

$p_{N}(0)=p_{i}(0)+C.$ (5.3)

However, since anisotropy has been kept the same, and thus still vanishes at
the center, in the new solution the two pressures increase or decrease by the
same amount $C$ and that the new solution is regular for any $C$. Thus one
gets a one-parameter family of regular solutions for any regular seed
solution.

On the other hand, the general solution of (2.8), where $f(r)$, $g(r)$, and
$h(r)$ are determined by $\rho_{i}(r)$ and $P_{i}(r)$, the general solution is

$p(r)=p_{i}(r)+\frac{Ce^{-\int\left(2g(r)p_{i}(r)+f_{0}(r)\right)dr}}{r^{2}\left(1+C\int{g(r)\over
r^{2}}e^{-\int\left(2g(r)p_{i}(r)+f_{0}(r)\right)dr}dr\right)},$ (5.4)

where $f_{0}(r)=f(r)-2/r$. It is easy to see that unless the second term
vanishes, there is no way one can get $p(r)=P_{i}(r)$. So there will be
nonzero anisotropy at the center and the one-parameter solution is not
regular.

### Riccati in $\Phi^{\prime}(r)$

As we mentioned earlier, the generalized TOV can also be seen as a Riccati
equation in $\Phi^{\prime}(r)$,

$\Phi^{\prime\prime}(r)+\left[\frac{r-m_{i}-4\pi
r^{3}\rho_{i}(r)}{r(r-2m_{i})}\right]\Phi^{\prime}(r)+\Phi^{\prime
2}(r)=\frac{8\pi r^{3}P_{i}(r)+4\pi r^{3}\rho_{i}(r)-m_{i}}{r^{2}(r-2m_{i})}.$
(5.5)

One can thus generate another one-parameter solution using a regular seed
solution and it is not difficult to see that

$\Phi^{\prime}_{N}(0)=\Phi_{i}^{\prime}(0)$ (5.6)

and thus the resulting one-parameter solutions can be regular. However, since
this uses the same input $(\rho(r),P(r))$ as the Riccati equation in $p(r)$,
all solution of the latter can be combined with $\rho(r)$ algebraically via
(2.7) to give the same $\Phi^{\prime}(r)$ and hence one will get the same one-
parameter family of solutions. We, thus, have the following theorem:

Theorem 5.1: For any (an)isotropic fluid solution of the Einstein equations in
canonical coordinates that is regular at the center, one can generate a one-
parameter family of regular (an)isotropic solutions by solving the
(generalized) TOV equation as Riccati equations in $p(r)$ or
$\Phi^{\prime}(r)$.

This would apply to all solutions including those generated using equations of
state. The resulting one-parameter solutions, however, will not in general
obey the equation of state of the seed.

##  6 Physical Solutions: Two New Classes of Anisotropic Solutions

In addition to central regularity and matching with the Schwarzschild metric
at the boundary, the following properties are often sought in a physical fluid
sphere (see, for example, [10, 17]):

  1. 1.

$\rho(r)$, $p(r)$ and $P(r)$ are positive and monotonically decreasing.

  2. 2.

The solution is regular at the center.

  3. 3.

$p(r)$ vanishes at the boundary.202020$P(r_{b})$ is not necessarily required
to vanish.

  4. 4.

Fluid variables satisfy the energy conditions: $\rho+p+2P\geq 0$ and $\rho\geq
p+2P$.

  5. 5.

Speed of sound is less than speed of light, i.e., $0\leq dp/d\rho\leq 1$ and
$0\leq dP/d\rho\leq 1$.212121Actually, the quantities, $dp/d\rho$ and
$dP/d\rho$, don’t necessarily represent sound speed; additional significant
assumptions (e.g., barotropic equation of state or adiabatic star) are
required to enforce the relationship, $0\leq dp/d\rho\leq 1$ and $0\leq
dP/d\rho\leq 1$ (see, for example, section $3.8$ of [25]).

Apart from these conditions, it is required that the metric component $g_{rr}$
is positive (i.e., $r>2m$). However, this is a consequence of conditions (1)
and (2) (see, for example, [1], where it was shown, through regularity
analysis of the Riccati equation (2.8), that for positive density,
differentiable tangential pressure and finite isotropic central pressure,
$r>2m$).

### Example I

Consider the following density and radial pressure as two input functions:

$\rho(r)=\alpha\left(r_{b}^{2}-r^{2}\right),\quad
p(r)=\beta\left(r_{b}^{2}-r^{2}\right)^{2},\quad\alpha,\beta>0,\quad
r\in[0,r_{b}],$ (6.1)

where $r_{b}>0$ is a boundary. Thus, the radial pressure and density follow
the relationship, $p\propto\rho^{2}$, and vanish simultaneously at the
boundary. Furthermore,

$\rho^{\prime}(r)=-2\alpha\,r<0,\quad
p^{\prime}(r)=-4\beta\,r\left(r_{b}^{2}-r^{2}\right)<0.$ (6.2)

From condition (3.14), the solution is regular the center. The mass function
is given by

$m(r)=\left(\frac{4\pi\alpha\,r_{b}^{2}}{3}\right)r^{3}-\left(\frac{4\pi\alpha}{5}\right)r^{5}.$
(6.3)

Thus, the mass function is same as in the Tolman-VII solution 222222This
solution is one of $16$ isotropic solutions (out of $127$) which satisfy all
the physical conditions listed above [10]. [28]. It turns out that the
tangential pressure vanishes at the boundary as well and is given by

$P(r)=\frac{\left(r^{2}-r_{b}^{2}\right)\left[9\beta\,r^{2}\left(5+4\pi\alpha
r^{4}\right)-\beta\,r_{b}^{2}\left(15+68\pi\alpha
r^{4}\right)+2\pi\alpha^{2}r^{2}\left(3r^{2}-5r_{b}^{2}\right)+30\pi\beta^{2}\,r^{2}\left(r^{2}-r_{b}^{2}\right)^{3}\right]}{15+8\pi\alpha\,r^{2}\left(3r^{2}-5r_{b}^{2}\right)}.$
(6.4)

Furthermore, the $g_{tt}$ component of metric, from (4.30) reduces to

$g_{tt}=-\left(1-\frac{2m}{r}\right)^{-\xi}\exp\left[\frac{5\beta}{2\alpha}\left(r^{2}-\frac{5}{6}r_{b}^{2}\right)+\frac{5}{12\alpha\delta}\left(\frac{\beta
r_{b}^{4}}{6}+\alpha
r_{b}^{2}-6\beta\delta^{2}\right)\tan^{-1}\left(\frac{6r^{2}-5r_{b}^{2}}{6\delta}\right)\right]$
(6.5)

where

$\xi=\frac{1}{4}+\frac{5\beta
r_{b}^{2}}{12\alpha}\quad\text{and}\quad\delta^{2}=\frac{5}{8\pi\alpha}-\frac{25}{36}\;r_{b}^{4}\,>0\quad\text{i.e.,}\quad\alpha<\frac{9}{10\pi
r_{b}^{4}}.$ (6.6)

Figure 1:  Here, $r_{b}=1$, $\alpha=0.15$, and $\beta=0.018$. Thus, the
Schwarzschild mass is $M=0.25$ and radius is $R=3.9\,M$. As can be seen, this
solution satisfy all the physical conditions.

Figure 2: If the relationship between the radial pressure and density is
changed to $p\propto\rho^{3}$, while $\rho(r)$ is still the same, then the new
solution still satisfy all the physical conditions. Here, $r_{b}=1$,
$\alpha=0.15$ and $\beta=0.018$.

### Example II

Consider the following potential function and density profile as two input
functions:232323 Parametrizing this way yields $\Phi^{\prime}(r)>0$, recall
equation (2.7).

$\Phi(r)=\gamma^{2}r^{2},\quad\rho(r)=\rho_{0}-\left(\rho_{0}-\rho_{b}\right)\left(\frac{r}{r_{b}}\right)^{2},\quad
0\leq\rho_{b}<\rho_{0},\quad r\in[0,r_{b}].$ (6.7)

In other words, the $g_{tt}$ component of the metric is Gaussian. Moreover,

$\Phi^{\prime}(0)=0,\quad\Phi^{\prime\prime}(r)=2\gamma^{2},\quad\rho(0)=\rho_{0},\quad\rho(r_{b})=\rho_{b}.$
(6.8)

Therefore, from (3.14), the solution is regular at the center. Since

$\rho^{\prime}(r)=-\frac{2r}{r_{b}}(\rho_{0}-\rho_{b})\leq 0,$ (6.9)

density is monotonically decreasing and $\rho(r)\geq\rho(r_{b})\geq 0$. We
will apply the boundary condition $p(r_{b})=0$ to determine $\rho_{b}$ in
terms of $\rho_{0}$ and $r_{b}$. The second boundary condition
$\exp[{2\Phi(r_{b})}]=1-(2m(r_{b})/r_{b})$ can be satisfied by rescaling the
time coordinate with appropriate constant. In other words, the two parameters
one can choose in this approach are $\rho_{0}$ and $r_{b}$. The mass function
reduces to

$m(r)=\frac{4\pi\rho_{0}r^{3}}{3}-\frac{4\pi
r^{5}}{5}\left(\frac{\rho_{0}-\rho_{b}}{r_{b}^{2}}\right)=\frac{r^{3}}{a^{2}}-\frac{r^{5}}{b^{2}}.$
(6.10)

Once again, the mass function is the same as that in the Tolman-VII solution
(but the density is non-zero at the boundary) where

$a^{2}=\frac{3}{4\pi\rho_{0}},\quad
b^{2}=\frac{5r_{b}^{2}}{4\pi(\rho_{0}-\rho_{b})}.$ (6.11)

From the field equations

$\displaystyle p$ $\displaystyle=$
$\displaystyle\left(\frac{8\gamma^{2}}{b^{2}}\right)r^{4}+\left({2\over
b^{2}}-{8\gamma^{2}\over
a^{2}}\right)r^{2}+\left(4\gamma^{2}-\frac{2}{a^{2}}\right)$ (6.12)
$\displaystyle P$ $\displaystyle=$
$\displaystyle\left(\frac{8\gamma^{4}}{b^{2}}\right)r^{6}+\left({16\gamma^{2}\over
b^{2}}-{8\gamma^{4}\over a^{2}}\right)r^{4}+\left({4\over
b^{2}}+4\gamma^{4}-{12\gamma^{2}\over
a^{2}}\right)r^{2}+\left(4\gamma^{2}-\frac{2}{a^{2}}\right).$ (6.13)

We take $\gamma=1$ for simplicity. Thus, $p(0)=P(0)>0$ if

$a^{2}>\frac{1}{2}\quad\text{or}\quad\rho_{0}<\frac{3}{2\pi}.$ (6.14)

Since the radial pressure is of degree four in $r$ and is symmetric under the
transformation $r\to-r$, equation $p(r)=0$ has a maximum of two roots for
$r\geq 0$. If $r_{b}$ is the smaller root, $\rho_{b}$ can be found in terms of
$\rho_{0}$ and $r_{b}$. One way to find specific solution would be to choose
$\rho_{0}$ first, and then choose $r_{b}$ such that $0\leq\rho_{b}<\rho_{0}$.
Here, we directly choose values of $a^{2}$ and $b^{2}$ such that all physical
conditions, except (21), are satisfied.

Figure 3:  Here, $\gamma^{2}=1$, $a^{2}=0.519$ and $b^{2}=0.325$. The boundary
is at $r_{b}=0.129$ and the total mass is $M=0.004$. Therefore, the boundary
in terms of total mass is given by $r_{b}=(32.269)M$. All the physical
conditions, except (21), are satisfied.

##  7 Conclusion

We studied all possible combinations of the basic functions in curvature
coordinates to produce all solutions of the (an)isotropic fluid system via
quadrature and paid particular attention to sub-algorithms that can generate
all (and only) regular solutions. This first required us to revisit and take a
closer look at the geometric conditions of central regularity under the
Einstein equations of the system and work out different equivalent conditions
of regularity. The conditions of regularity, in general, are three conditions
on two functions, with the only exception in the case of isotropy when two
conditions on the potential function are found to be sufficient. We also found
that one can allow $p(r)$ to be non-differentiable at $r=0$ as long as
$\lim_{r\to 0}rp^{\prime}(r)=0$ to generate regular solution; one should
include such functions to generate all (regular) solutions when using $p(r)$
as one of the inputs.

There is no discussion in the literature on the regularity of the existing
algorithm for the anisotropic system [18]. We, therefore, revisited this first
using a slightly different set of variables, $\Phi^{\prime}(r)$ and $P(r)$. We
noted that the parameter arising in the general solution of this algorithm
could be interpreted as the core mass in the presence of a nonzero core,
which, therefore, cannot be negative or unbounded on physical grounds. For
solutions starting from $r=0$, i.e., complete spheres (possibly terminating at
a finite radius), the same analysis gives multiple solutions with zero central
mass, and they are regular if and only if the parameter is set to zero. Thus,
all anisotropic spheres are produced in a 1-1 fashion using $\Phi^{\prime}(r)$
and $P(r)$ as generating functions. When $\Phi^{\prime}(r)$ and $\chi(r)$ are
used as input functions, as in [18], the parameter in the general solution
reverses its role and does not have to be zero; one, thus, have a one-
parameter family of solutions that are regular under the usual initial
conditions on $\Phi^{\prime}(r)$ and $\chi(r)$ for regularity. In addition, in
comparison with the algorithm using $\Phi^{\prime}(r)$ and $P(r)$, we can now
interpret the parameter as the central density. This also means that for pure
isotropy, we get a one-parameter family of regular solutions, and the
parameter, the central density, remains free; this analysis thus complements
the analysis for the isotropic system in [20].

We found that the three other combinations—$(\rho(r)$, $p(r)$),
$(\rho(r),\Phi^{\prime}(r))$, and $(p(r),\Phi^{\prime}(r))$—lead to another
algorithm for the anisotropic system. The equivalency of the three pairs in
this algorithm arises precisely because of the same algebraic relation between
these three functions that made the different equivalent sets of regularity
conditions possible.242424We could not find a formal statement in the
literature to this effect, but that $(\rho(r),p(r))$ can serve as generating
functions for all anisotropic solutions is amply clear if one reads the
generalized TOV equation as the defining algebraic equation for $P(r)$ in
terms of $\rho(r)$ and $p(r)$.This fact, for example, was noted in [5] and
this was used in [22] to generate a certain class of anisotropic solutions
(using the Newtonian relationship $p^{\prime}=-m\rho/r^{2}$). The line element
(4.30) is very suitable for this algorithm. We noted that in regular solutions
finiteness of $p(0)$ has to be be supplemented by $\lim_{r\to
0}rp^{\prime}(r)=0$. This extra condition comes from our analysis of (2.6) and
is reconfirmed by the tensorial conditions of regularity. As for the other two
choices—$(\rho(r),\Phi^{\prime}(r))$ and $(p(r),\Phi^{\prime}(r))$—we find
that they are algebraically equivalent to $(\rho(r)$, $p(r)$) and can be used
with their accompanying regularity conditions to generate all regular
solutions. It is quite conceivable that one would convert any such pair to
$(\rho(r)$, $p(r)$) and use (4.30). Note that $\rho(r)$ and $p(r)$ are output
variables in the existing algorithms of isotropic and anisotropic systems and
as such the line element (4.30) is especially significant only in this new
algorithm for the anisotropic system. The conditions of regularity, as well as
any energy condition, become a simple matter of choosing the appropriate
algebraic forms of two generating functions in this algorithm. These exhaust
all possible combinations of the basic functions of the anisotropic system.
These results can be used with equations of state with little change.

An interesting feature of a Riccati equation is that if a particular solution
is known one can work out the general solution. We have now three
algorithms—with $(\rho(r)$, $p(r)$), $(\Phi^{\prime}(r),P(r))$, and
$(\Phi^{\prime}(r),\chi(r))$—that generate all solutions and each of them is
capable of generating all regular solutions. Any given solution generated by
them can be seen as the particular solution of the Riccati equation (in $p(r)$
or in $\Phi^{\prime}(r)$) and one gets two distinct one-parameter family of
solutions. We found that one of them necessarily results in solutions that are
regular at the center, for both isotropy and anisotropy, if the particular
solution used is regular.

Finally, following the new algorithm, we used $(\rho(r),\Phi(r))$ and
$(\rho(r),p(r))$ to generate examples of regular solutions that satisfy other
physical conditions in addition to the condition of central regularity. In
another work, we will discuss various maps between and within isotropic and
anisotropic solutions and see how they can facilitate the process of finding
solutions that satisfy all physical conditions.

## Acknowledgements

We thank Pankaj Joshi, Kayll Lake, and Malcolm MacCallum for useful
communications. Special thanks to Charlie Brewer for useful comments on the
manuscript. RS acknowledges the support of a Julia Williams Van Ness Merit
Scholarship.

## References

  * [1] T. W. Baumgarte and A. D. Rendall, “Regularity of spherically symmetric static solutions of the Einstein equations,” Class. Quantum Grav. 10 327 (1993)
  * [2] S. Berger, R. Hojman and J. Santamarina, “General Exact Solutions of Einstein Equations for Static Perfect Fluids With Spherical Symmetry,” J. Math. Phys. 28, 2949 (1987). doi:10.1063/1.527697
  * [3] C. G. Boehmer and T. Harko, “On Einstein Clusters as Galactic Dark Matter Halos,” Mon. Not. Roy. Astron. Soc. 379, 393 (2007)
  * [4] C. G. Boehmer and T. Harko, “Bounds on the Basic Physical Parameters for Anisotropic Compact General Relativistic Objects,” Class. Quant. Grav. 23, 6479 (2006)
  * [5] P. Boonserm, M. Visser and S. Weinfurtner, “Solution Generating Theorems for the TOV Equation,” Phys. Rev. D 76, 044024 (2007)
  * [6] R. L. Bowers and E. P. T. Liang, “Anisotropic Spheres in General Relativity,” Astrophys. J. 188, 657-665 (1974)
  * [7] K. A. Bronnikov and S. G. Rubin (2012), _Black Holes, Cosmology and Extra Dimensions_ (World Scientific, Singapore).
  * [8] R. Casadio, E. Contreras, J. Ovalle, A. Sotomayor and Z. Stuchlick, “Isotropization and Change of Complexity by Gravitational Decoupling,” Eur. Phys. J. C 79, no.10, 826 (2019)
  * [9] E. Coddington and N. Levinson (1955) _Theory of Ordinary Differential Equations_ (McGraw-Hill, New York).
  * [10] 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)
  * [11] K. Dev and M. Gleiser, “Anisotropic Stars: Exact Solutions,” Gen. Rel. Grav. 34, 1793 (2002)
  * [12] A. Einstein, “On a Stationary System with Spherical Symmetry Consisting of Many Gravitating Masses,” Annals Math. 40, 922 (1939).
  * [13] G. F. R. Ellis, “Relativistic Cosmology,” Gen. Rel. Grav. 41, 581 (2009)
  * [14] M. Finch, “The Painleve-Gambier Equation and the Relativistic Static Fluid Sphere,” Ph.D. Thesis (1987)
  * [15] G. Fodor, “Generating spherically symmetric static perfect fluid solutions,” [arXiv:gr-qc/0011040 [gr-qc]].
  * [16] J. B. Hartle (2003), _Gravity: An Introduction to Einstein’s General Relativity_ , (Addison-Wesley, San Francisco).
  * [17] T. Harko and M. K. Mak, “Anisotropic relativistic stellar models,” Annalen Phys. 11, 3-13 (2002)
  * [18] L. Herrera, J. Ospino and A. Di Prisco, “All Static Spherically Symmetric Anisotropic Solutions of Einstein’s Equations,” Phys. Rev. D 77, 027502 (2008)
  * [19] L. Herrera and N. O. Santos, “Local anisotropy in self-gravitating systems,” Phys. Rept. 286, 53-130 (1997)
  * [20] K. Lake, “All Static Spherically Symmetric Perfect Fluid Solutions of Einstein’s Equations,” Phys. Rev. D 67, 104015 (2003)
  * [21] K. Lake and P. Musgrave, “The Regularity of Static Spherically Cylindrically and Plane Symmetric Spacetimes at the Origin,” Gen. Relat. Gravit. 26, 917-925 (1994)
  * [22] K. Lake, “Generating Static Spherically Symmetric Anisotropic Solutions of Einstein’s Equations from Isotropic Newtonian Solutions,” Phys. Rev. D 80, 064039 (2009)
  * [23] D. Martin and M. Visser, “Algorithmic construction of static perfect fluid spheres,” Phys. Rev. D 69, 104028 (2004)
  * [24] C. W. Misner, K. S. Thorne and J. A. Wheeler, _Gravitation_ , San Francisco (1973)
  * [25] S. Rahman and M. Visser, “Space-time Geometry of Static Fluid Spheres,” Class. Quant. Grav. 19, 935-952 (2002) doi:10.1088/0264-9381/19/5/307 [arXiv:gr-qc/0103065 [gr-qc]].
  * [26] H. Rago, “Anisotropic Spheres in General Relativity,” Astrophys Space Sci 183, 333-338 (1991).
  * [27] H. Stephani, D. Kramer, M. A. H. MacCallum, C. Hoenselaers and E. Herlt (2003). _Exact Solutions of Einstein’s Field Equations_ (Cambridge University Press, UK).
  * [28] R. C. Tolman, “Static solutions of Einstein’s field equations for spheres of fluid,” Phys. Rev. 55, 364-373 (1939)
  * [29] S. Viaggiu, “Modelling Usual and Unusual Anisotropic Spheres,” Int. J. Mod. Phys. D 18, 275 (2009)
  * [30] R. M. Wald (1984) _General Relativity_ (University of Chicago Press, Chicago).
  * [31] M. Wyman, “Radially Symmetric Distributions of Matter,” Phys. Rev. 75, 1930-1936 (1949)
  * [32] D. Zwillinger (1998), _Handbook of Differential Equations_ , 3rd ed. (Academic Press, San Diego).