# On anisotropic two-fluid stellar objects in General Relativity

Nolene F. Naidu<EMAIL_ADDRESS>Department of Mathematics and
Applied Mathematics, University of Cape Town  Sante Carloni
<EMAIL_ADDRESS>DIME Sez. Metodi e Modelli Matematici, Università di
Genova,
Via All’Opera Pia 15, 16145 - Genoa, (Italy).  Peter Dunsby  Department of
Mathematics and Applied Mathematics, University of Cape Town

###### Abstract

We apply the 1+1+2 covariant semi-tetrad approach to describe a general static
and spherically symmetric relativistic stellar object which contains two
fluids with anisotropic pressure. The corresponding Tolman-Oppenheimer-Volkoff
equations are then obtained in covariant form for the anisotropic case. These
equations are used to obtain new exact solutions using direct resolution and
reconstruction techniques. Finally, we show that three of the generating
theorems known for the single fluid case can also be used to obtain two-fluid
solutions from single fluid ones.

##  I Introduction

Studying the internal structure of stars and compact objects in a relativistic
context is an important area of relativistic astrophysics. It is, however,
notoriously complex, particularly when approached using analytic methods. The
fundamental equations that describe the relativistic hydrostatic equilibrium,
the Tolman-Oppenheimer-Volkoff (TOV) equations [1, 2], though quasilinear,
present many inherent challenges in their analytical resolution. To date, only
a few physically relevant exact solutions have been found, some of them almost
a century-old [3].

The challenge is intensified by the fact that realistic compact objects have a
structure that comprises different regions with wildly different properties.
For example, neutron stars possess an inner and outer core as well as an inner
and outer crust; the compositions of which range from standard matter to
superfluid neutronium [4]. The differences in composition and conditions
between these regions are so marked that it is not advantageous to attempt to
describe them using a single equation of state. We suppose, instead, that
these compact objects comprise different fluids, each with its own equation of
state. In this picture, we expect neutron stars and other relativistic compact
objects to present a so-called “shell” structure in which every layer is
composed of different fluids or a mixture of fluids. Since we will have a
different solution of the TOV equations for each of these fluids, it becomes
crucial to join these solutions to obtain a complete solution for the entire
object. The key tool for this task is the junction conditions, i.e., the
conditions necessary to solder together different spacetimes. There are
different formulations of the junction conditions, but the most used is the
Israel formulation [5]. These conditions will also be used to match the
solutions associated with the interior of the compact object to the exterior,
typically represented by the vacuum Schwarzschild [6], or the radiating Vaidya
[7] solutions. In the following, we will use the term “complete” solution for
any solution obtained in this way, so that it describes the entire spacetime
associated with the interior of a given compact object.

Another important aspect to consider in a realistic model of compact stellar
objects is the inclusion of anisotropies, arising from the difference between
radial and transverse pressures. There are many possible reasons for the
existence of anisotropy in astrophysical objects, such as strong
electromagnetic fields and the exotic thermodynamical properties in the state
of matter [8] (and references therein). Letelier [9] also proposed that a
mixture of real gas can exhibit anisotropy.

Anisotropy adds extra degrees of freedom that make the search for analytical
anisotropic interior solutions for these objects highly non-trivial. Compact
objects which admit local anisotropies present modifications to their physical
parameters that affect the redshift, critical mass, and the stability of the
stellar objects [10]. Exotic objects, such as boson stars and generalized
classes of gravastars, are intrinsically anisotropic in pressure [11]. In
general, for a mixture of two perfect fluids in relative motion, the total
energy-momentum tensor in the comoving frame of either fluid is anisotropic in
the pressure [9].

There has been steady incremental progress on the analytical study of
anisotropic solutions of the TOV equations. Harko and Mak [12] assume a
specific value for the anisotropy parameter and give a class of physically
viable, exact solutions. Rago [13] introduced an anisotropic function to
generalize the Tolman V solution to include anisotropy, and a generating
function, forming part of the foundation for the generating theorems in [14,
15]. The work by Viaggiu [10] is a perturbative analysis applied to known
”seed” solutions such as Florides and Tolman IV. In the paper by Das et al.
[16], the anisotropic extension to Tolman IV is given. Their model is compared
to pulsar observational data and satisfies the physical viability conditions
given in [3].

The most substantial advances in this research line, however, have been
achieved numerically. There has been relevant progress with simulations and
approximations in this regard (see e.g. [17]-[22] and references therein). In
this context, analytical studies of the TOV equations are important in
addition to numerical studies for several reasons. Exact solutions can be used
to explore the full parameter space for a given metric rather than a single
set of values. They can also be utilized to test numerical codes, particularly
when they involve new approximations schemes. Furthermore, analytical
solutions can be tested using observational data to constrain the parameters
in these models [16].

In this paper, we will construct a covariant formulation of the TOV equations
describing a multifluid compact object with anisotropies. This work will be an
extension of [23], in which we developed the covariant TOV equations for
multifluid isotropic compact objects. We explore direct resolution strategies
to obtain two-fluid solutions from known single fluid ones with consideration
of the energy conditions to ensure the obtained solutions are physically
viable. In order to be considered physically viable, the behaviour of the
thermodynamical quantities in these solutions must satisfy the constraints
given for at least one set of the parameters. We extend the reconstruction
algorithms to include anisotropy and give an example of a solution,
constructed by using a combination of known single fluid solutions, the Tolman
IV and the interior Schwarzschild solutions, applied to the two-fluid model.
In Boonserm et al. [14, 15] it was shown that one can obtain new single fluid
solutions from known (seed) solutions using generating theorems. Here, we give
three examples of variations to the generating theorems, and show that they
can be applied to the two-fluid anisotropic case.

The cornerstone of this paper is the covariant 1+1+2 approach. This formalism
was pioneered by Clarkson and Barrett in 2003 [24], where it was applied to
the study of perturbations of Schwarzschild black holes. The fact that this
approach is ideally suited to the study of spherically symmetric systems and
their perturbations led Carloni to apply this framework to study single fluid
isotropic [25] and anisotropic stellar models [26], leading to physically
viable exact solutions in both cases.

The following analysis involves simplifying the structure described above by
using one or two fluids with a given equation of state. It is relatively easy
to generalise the two fluid equations to the case of the $N$ fluids. While
such generalisation would provide the opportunity to examine a wider array of
objects, it is well known that neutron stars require just two fluids to be
described at a satisfactory degree of accuracy [27]. The cases with more than
two fluids can be explored to potentially give a higher degree of accuracy to
models in the future. In order to give a truly complete and realistic
description, we would have to include other characteristics such as rotation,
phase transitions, charge, flux, interactions, etc. While the two-fluid
equations obtained can be easily generalized to include some of these effects,
we will not include such factors in this work. These aspects will be the
subject of future research.

The outline of this paper is as follows: In Sec. II the 1+1+2 semi-tetrad
equations for the two-fluid solution are presented. Section III deals with the
conditions of physical viability for a given solution of the TOV equations.
Section IV, instead, deals with the generalisation of Israel’s junction
conditions to the multifluid case. Section V contains the expressions for the
Bowers-Liang solution, an element of the class of Florides solutions, the
interior Schwarzschild solution, and the Tolman IV solution. In Sec. VI, we
explore some direct resolution strategies applied to the solutions in Sec. V.
Sec. VII contains solution reconstruction techniques. The generating theorems
of [14] are discussed and applied to a two-fluid system in Sec. VIII. A
discussion and some concluding remarks can be found in Sec. IX. The main
equations of the 1+1+2 formalism are presented in Appendix A and the $N$ fluid
generalisation of the equations in Appendix B. Appendix C contains the full
solutions for Sec. VI, while Appendix D has the full expressions for solutions
obtained in Sec. VII.

##  II The 1+1+2 equations for the two-fluid case

The 1+1+2 equations can be written easily for any number of fluids (see
Appendix B for the extension to $N$ fluids). In this section and what follows,
we will consider the application to two fluids. In [27], this is shown as
sufficient to describe systems such as neutron stars, which are one of the
primary applications of the TOV equations. We begin by defining a time-like
threading vector field $u^{a}$ that is associated to the observer’s congruence
with $u^{a}u_{a}=-1$, and a space-like vector $e^{a}$ with $e_{a}e^{a}=1$. The
$u^{a}$ and $e^{a}$ congruences allow one to foliate the spacetime in
hypersurfaces whose geometry is defined by the two projection tensors

$\begin{split}{h^{a}}_{b}&={g^{a}}_{b}+u^{a}u_{b}\hskip 5.69046pt,\hskip
5.69046pt{h^{a}}_{a}=3,\\\
{N_{a}}^{b}&={h_{a}}^{b}-e_{a}e^{b}={g_{a}}^{b}+u_{a}u^{b}-e_{a}e^{b}\hskip
2.84544pt,\hskip 2.84544pt{N^{a}}_{a}=2,\end{split}$ (1)

Here ${h^{a}}_{b}$ represents the metric of the 3-spaces orthogonal to
$u^{a}$, and ${N_{a}}^{b}$ represents the metric of the 2-spaces orthogonal to
$u^{a}$ and $e^{a}$. Any tensorial object may now be split according to the
above foliations, as described in [25]. The covariant time derivative,
orthogonally projected covariant derivative, hat-derivative and
$\delta$-derivative are given by

$\displaystyle{\dot{X}{}^{a..b}}_{c..d}$ $\displaystyle\equiv$ $\displaystyle
u^{e}{\nabla}_{e}{X^{a..b}}_{c..d},$ $\displaystyle D_{e}{X^{a..b}}_{c..d}$
$\displaystyle\equiv$
$\displaystyle{h^{a}}_{f}...{h^{b}}_{g}{h^{p}}_{c}...{h^{q}}_{d}{h^{r}}_{e}{\nabla}_{r}{X^{f..g}}_{p..q},$
$\displaystyle{\hat{X}{}_{a..b}}^{c..d}$ $\displaystyle\equiv$ $\displaystyle
e^{f}D_{f}{X_{a..b}}^{c..d},$ $\displaystyle{\delta}_{e}{X_{a..b}}^{c..d}$
$\displaystyle\equiv$
$\displaystyle{N_{a}}^{f}...{N_{b}}^{g}{N_{i}}^{c}...{N_{j}}^{d}{N_{e}}^{p}D_{p}{X_{f..g}}^{i..j},\hskip
8.5359pt$ (2)

respectively.

The full set kinematical variables associated to the covariant derivatives of
$u^{a}$ and $e^{a}$ is given in Appendix A.

The energy-momentum tensor can be decomposed as

$\displaystyle T_{ab}=$ $\displaystyle\mu
u_{a}u_{b}+(p+\Pi)e_{a}e_{b}+\left(p-\frac{1}{2}\Pi\right)N_{ab}$ (3)
$\displaystyle+2Qe_{(a}u_{b)}+2Q_{(a}u_{b)}+2\Pi_{(a}e_{b)}+\Pi_{ab},$

where the symmetrisation over the indices of a tensor is represented as
$T_{(ab)}=\frac{1}{2}(T_{ab}+T_{ba})$. The matter variables are given by

$\displaystyle\mu$ $\displaystyle=$ $\displaystyle T_{ab}u^{a}u^{b},$
$\displaystyle p$ $\displaystyle=$
$\displaystyle\frac{1}{3}T_{ab}\left(e^{a}e^{b}+N^{ab}\right),$
$\displaystyle\Pi$ $\displaystyle=$
$\displaystyle\frac{1}{3}T_{ab}\left(2e^{a}e^{b}-N^{ab}\right),$
$\displaystyle Q$ $\displaystyle=$ $\displaystyle-T_{ab}e^{a}u^{b},$
$\displaystyle Q_{a}$ $\displaystyle=$ $\displaystyle-T_{cd}{N^{c}}_{a}u^{d},$
$\displaystyle\Pi_{a}$ $\displaystyle=$ $\displaystyle
T_{cd}{N^{c}}_{a}e^{d},$ $\displaystyle\Pi_{ab}$ $\displaystyle=$
$\displaystyle T_{\\{ab\\}},$ (4)

where $\mu$ is the density, $p$ is the pressure, $Q$ and $Q_{a}$ represent the
scalar and vector parts of the heat flux, and $\Pi$ and $\Pi_{a}$ represent
the the scalar and vector components of the anisotropic pressure,
respectively. The curly brackets $\\{\\}$ denote the Projected Symmetric
Trace-Free part of a tensor with respect to $e^{a}$.

The radial and tangential pressures which appear as sources in the
gravitational field equations are given by

$\displaystyle p_{r}$ $\displaystyle=p+\Pi,$ (5a) $\displaystyle p_{\perp}$
$\displaystyle=p-\frac{1}{2}\Pi.$ (5b)

Spacetimes which have a unique preferred spatial direction at each point and
exhibit local rotational symmetry around such direction are called “LRS
spacetimes” and can be treated easily within the 1+1+2 formalism. A specific
class of non-rotating, non-twisting LRS spacetimes called LRS II spacetimes,
has the additional property that all the 1+1+2 vectors and tensors and the
quantities $\Omega,\xi,\mathcal{H}$ vanish identically. If one further chooses
to consider static, spherically symmetric spacetimes, then $\Theta,\Sigma$ and
$Q$ are zero, as well as the dot derivatives of scalar quantities. In the
following, we will deal specifically with static and spherically symmetric
LRSII spacetimes.

For these spacetimes, the relevant kinematical variables are

$\displaystyle\mathcal{A}$ $\displaystyle=e_{a}\dot{u}^{a},$ (6a)
$\displaystyle\phi$ $\displaystyle={\delta}_{a}e^{a},$ (6b)
$\displaystyle\mathcal{E}$ $\displaystyle={C}_{acbd}u^{c}u^{d}e^{a}e^{b},$
(6c)

where $C_{acbd}$ is the Weyl tensor, and the system (131) for two fluids with
zero total flux, reduces to

$\displaystyle\hat{\phi}=$
$\displaystyle-\frac{1}{2}{\phi}^{2}-\frac{2}{3}({\mu}_{1}+{\mu}_{2})-\mathcal{E}-\frac{1}{2}({\Pi}_{1}+{\Pi}_{2}),$
(7a) $\displaystyle\hat{\mathcal{E}}=$
$\displaystyle\frac{1}{3}(\hat{\mu}_{1}-\hat{\mu}_{2})-\frac{1}{2}(\hat{\Pi}_{1}+\hat{\Pi}_{2})$
$\displaystyle-\frac{3}{2}{\phi}\left(\mathcal{E}+\frac{1}{2}{\Pi}_{1}+\frac{1}{2}{\Pi}_{2}\right),$
(7b) $\displaystyle-\mathcal{A}\phi$
$\displaystyle+\frac{1}{3}({\mu}_{1}+3p_{1})+\frac{1}{3}({\mu}_{2}+3p_{2})$
$\displaystyle-\mathcal{E}+\frac{1}{2}{\Pi}_{1}+\frac{1}{2}{\Pi}_{2}=0,$ (7c)
$\displaystyle\hat{\mathcal{A}}=$
$\displaystyle-\mathcal{A}\left(\mathcal{A}+\phi\right)+\frac{1}{2}({\mu}_{1}+3p_{1})+\frac{1}{2}({\mu}_{2}+3p_{2}),$
(7d) $\displaystyle\hat{p}_{1}+\hat{\Pi}_{1}=$
$\displaystyle-{\Pi}_{1}\left(\frac{3}{2}\phi+\mathcal{A}\right)-\mathcal{A}\left({\mu}_{1}+p_{1}\right)+j_{e}^{(1,2)},$
(7e) $\displaystyle\hat{p}_{2}+\hat{\Pi}_{2}=$
$\displaystyle-{\Pi}_{2}\left(\frac{3}{2}\phi+\mathcal{A}\right)-\mathcal{A}\left({\mu}_{2}+p_{2}\right)-j_{e}^{(1,2)},$
(7f) $\displaystyle K=$
$\displaystyle\frac{1}{3}({\mu}_{1}+{\mu}_{2})-\mathcal{E}-\frac{1}{2}({\Pi}_{1}+{\Pi}_{2})+\frac{1}{4}{\phi}^{2},$
(7g) $\displaystyle\hat{Q}_{1}=$ $\displaystyle-
Q_{1}\left(\phi+2\mathcal{A}\right)+j_{u}^{(1,2)},$ (7h) $\displaystyle
Q_{2}=$ $\displaystyle-Q_{1},$ (7i)

where $j_{u}^{(1,2)}$ and $j_{e}^{(1,2)}$ represent the interaction terms
between fluid 1 and 2. Next, we introduce a useful parameter, $\rho$, such
that $\hat{X}=\phi X_{,\rho}$. In this way, the equation for the Gauss’
curvature $K$ (133) can be solved to give

$K=K_{0}^{-1}e^{-\rho}.$ (8)

Assuming $\phi\neq 0$111This assumption is might appear arbitrary. However, it
can be shown that $\phi=0$ can be related to the presence of horizons in
spacetime. As there will be no horizons in the spacetimes we will consider we
can make this assumption without loss of generality. and defining the
variables

$\displaystyle X$ $\displaystyle=\frac{{\phi}_{,\rho}}{\phi},$ $\displaystyle
Y$ $\displaystyle=\frac{\mathcal{A}}{\phi},$ $\displaystyle\mathcal{K}$
$\displaystyle=\frac{K}{{\phi}^{2}},$ (9a) $\displaystyle E$
$\displaystyle=\frac{\mathcal{E}}{{\phi}^{2}},$ $\displaystyle\mathbb{M}_{1}$
$\displaystyle=\frac{{\mu}_{1}}{{\phi}^{2}},$ $\displaystyle\mathbb{M}_{2}$
$\displaystyle=\frac{{\mu}_{2}}{{\phi}^{2}},$ (9b) $\displaystyle P_{1}$
$\displaystyle=\frac{p_{1}}{{\phi}^{2}},$ $\displaystyle P_{2}$
$\displaystyle=\frac{p_{2}}{{\phi}^{2}},$ $\displaystyle\mathbb{P}_{1}$
$\displaystyle=\frac{{\Pi}_{1}}{{\phi}^{2}},$ (9c)
$\displaystyle\mathbb{P}_{2}$ $\displaystyle=\frac{{\Pi}_{2}}{{\phi}^{2}},$
$\displaystyle\mathbb{Q}_{1}$ $\displaystyle=\frac{Q_{1}}{{\phi}^{2}},$
$\displaystyle\mathbb{Q}_{2}$ $\displaystyle=\frac{Q_{2}}{{\phi}^{2}},$ (9d)
$\displaystyle\mathbb{J}_{u}$
$\displaystyle=\frac{j_{u}^{(1,2)}}{{\phi}^{3}},$
$\displaystyle\mathbb{J}_{e}$
$\displaystyle=\frac{j_{e}^{(1,2)}}{{\phi}^{3}},$ (9e)

we can recast the system (7) as

$\displaystyle Y_{,\rho}=$
$\displaystyle-Y(X+Y+1)+\frac{1}{2}(\mathbb{M}_{1}+\mathbb{M}_{2})$
$\displaystyle+\frac{3}{2}(P_{1}+P_{2}),$ (10a)
$\displaystyle\mathcal{K}_{,\rho}=$ $\displaystyle-\mathcal{K}(1+2X),$ (10b)
$\displaystyle{P}_{1,\rho}+\mathbb{P}_{1,\rho}=$
$\displaystyle-Y(\mathbb{M}_{1}+P_{1})$
$\displaystyle-\mathbb{P}_{1}\left(2X+Y+\frac{3}{2}\right)-2XP_{1}+\mathbb{J}_{e},$
(10c) $\displaystyle{P}_{2,\rho}+\mathbb{P}_{2,\rho}=$
$\displaystyle-Y(\mathbb{M}_{2}+P_{2})$
$\displaystyle-\mathbb{P}_{2}\left(2X+Y+\frac{3}{2}\right)-2XP_{2}-\mathbb{J}_{e},$
(10d) $\displaystyle\mathbb{Q}_{1,\rho}=$
$\displaystyle-\mathbb{Q}_{1}(1+2X+2Y)+\mathbb{J}_{u},$ (10e)

with the following constraints

$\displaystyle
2(\mathbb{M}_{1}+\mathbb{M}_{2})+2(P_{1}+P_{2})+2(\mathbb{P}_{1}$
$\displaystyle+\mathbb{P}_{2})+2X$ (11) $\displaystyle-2Y+1=0,$ $\displaystyle
1-4\mathcal{K}-4(P_{1}+P_{2})+4Y-4(\mathbb{P}_{1}$
$\displaystyle+\mathbb{P}_{2})=0,$ (12) $\displaystyle
2(\mathbb{M}_{1}+\mathbb{M}_{2})+6(P_{1}+P_{2})+3(\mathbb{P}_{1}$
$\displaystyle+\mathbb{P}_{2})$ (13) $\displaystyle-6Y-6E=0.$

Rearranging (13), (11), and writing (7i) in terms of $\rho$, we can write

$\displaystyle E$ $\displaystyle=$
$\displaystyle\frac{1}{3}\left(\mathbb{M}_{1}+\mathbb{M}_{2}\right)+P_{1}+P_{2}+\frac{1}{2}\left(\mathbb{P}_{1}+\mathbb{P}_{2}\right)-Y,$
(14) $\displaystyle X$ $\displaystyle=$
$\displaystyle-\frac{1}{2}-(\mathbb{M}_{1}+\mathbb{M}_{2})-(\mathbb{P}_{1}+\mathbb{P}_{2})$
(15) $\displaystyle-(P_{1}+P_{2})+Y,$ $\displaystyle\mathbb{Q}_{1}$
$\displaystyle=$ $\displaystyle-\mathbb{Q}_{2},$ (16)

where differentiation with respect to $\rho$ is indicated by a comma and the
subscript $\rho$. The equation for the total pressure is given by

$\begin{split}{P}_{tot,\rho}+\mathbb{P}_{tot,\rho}&+P_{tot}(Y+2X)+Y\mathbb{M}_{tot}\\\
&+\mathbb{P}_{tot}\left(2X+Y+\frac{3}{2}\right)=0,\end{split}$ (17)

where, from now on the subscript “$tot$” of a given quantity represents the
total or sum of this quantity, for fluids 1 and 2. So, for example
$P_{tot}=P_{1}+P_{2}$. We can write the covariant equivalent of the TOV
equations by using the constraints to eliminate all the metric-related
variables except $\mathcal{K}$, resulting in

$\displaystyle{P}_{1,\rho}+\mathbb{P}_{1,\rho}$ $\displaystyle=$
$\displaystyle\mathbb{J}_{e}-{P_{1}}^{2}-{\mathbb{P}_{1}}^{2}+P_{1}\left[\mathbb{M}_{1}-2\mathbb{P}_{1}-3\mathcal{K}+\frac{7}{4}\right]+\mathbb{P}_{1}\left(\mathbb{M}_{1}-3\mathcal{K}+\frac{1}{4}\right)+\mathbb{M}_{1}\left(\frac{1}{4}-\mathcal{K}\right)$
(18) $\displaystyle-
P_{1}(P_{2}+\mathbb{P}_{2}-2\mathbb{M}_{2})-\mathbb{P}_{1}(P_{2}+\mathbb{P}_{2}-2\mathbb{M}_{2})-\mathbb{M}_{1}\left(P_{2}+\mathbb{P}_{2}\right),$
$\displaystyle{P}_{2,\rho}+\mathbb{P}_{2,\rho}$ $\displaystyle=$
$\displaystyle-\mathbb{J}_{e}-{P_{2}}^{2}-{\mathbb{P}_{2}}^{2}+P_{2}\left[\mathbb{M}_{2}-2\mathbb{P}_{2}-3\mathcal{K}+\frac{7}{4}\right]+\mathbb{P}_{2}\left(\mathbb{M}_{2}-3\mathcal{K}+\frac{1}{4}\right)+\mathbb{M}_{2}\left(\frac{1}{4}-\mathcal{K}\right)$
(19) $\displaystyle-
P_{2}(P_{1}+\mathbb{P}_{1}-2\mathbb{M}_{1})-\mathbb{P}_{2}(P_{1}+\mathbb{P}_{1}-2\mathbb{M}_{1})-\mathbb{M}_{2}\left(P_{1}+\mathbb{P}_{1}\right),$
$\displaystyle\mathcal{K}_{,\rho}$ $\displaystyle=$ $\displaystyle
2\mathcal{K}\left(\frac{1}{4}-\mathcal{K}+\mathbb{M}_{1}+\mathbb{M}_{2}\right),$
(20) $\displaystyle\mathbb{Q}_{1,\rho}$ $\displaystyle=$
$\displaystyle\mathbb{Q}_{1}\left[2\mathcal{K}-2(\mathbb{M}_{1}+\mathbb{M}_{2})-\frac{3}{2}\right]-\mathbb{J}_{u},$
(21) $\displaystyle{P}_{tot,\rho}+\mathbb{P}_{tot,\rho}$ $\displaystyle=$
$\displaystyle-{P_{tot}}^{2}-{\mathbb{P}_{tot}}^{2}+P_{tot}\left[\mathbb{M}_{tot}-2\mathbb{P}_{tot}-3\mathcal{K}+\frac{7}{4}\right]+\mathbb{P}_{tot}\left[\mathbb{M}_{tot}-3\mathcal{K}+\frac{1}{4}\right]+\mathbb{M}_{tot}\left(\frac{1}{4}-\mathcal{K}\right).$
(22)

The $N$ fluid generalization of the above equations can be found in Appendix
B. The above equations require an equation of state, relating the pressure and
density, in order to be complete. However, in the following we will not assume
a form for the equation of state unless it is absolutely necessary. In this
way we will be able to define several methods to solve the covariant TOV
equations (18-22).

In the following sections, we will present solutions of the above equations.
In the interest of aiding the physical interpretation of these solutions, we
present the following useful relations. For a generic metric tensor of the
form

$\displaystyle ds^{2}=$ $\displaystyle-k_{1}(x,t)dt^{2}+k_{2}(x,t)dx^{2}$
$\displaystyle+k_{3}(x,t)\left[dy^{2}+k_{4}(y)dz^{2}\right],$ (23)

with spacetime geometry represented by

$\displaystyle k_{4}(y)$
$\displaystyle=\left\\{\begin{array}[]{ll}\sin{y},&\mbox{closed geometry}\\\
y,&\mbox{flat geometry}\\\ \sinh{y},&\mbox{open geometry}\end{array}\right.,$
(27)

we can write [33]

$\begin{split}\phi=\frac{{k}_{3,x}}{k_{2}k_{3}},\quad\mathcal{A}=\frac{{k}_{1,x}}{2k_{2}k_{1}}.\end{split}$
(28)

As we will mainly use the parameter $\rho$ in our calculations, it is useful
to give the same formulae for the case in which the metric is written as

$ds^{2}=-k_{1}(\rho)dt^{2}+k_{2}(\rho)d\rho^{2}+k_{3}(\rho)d\Omega^{2},$ (29)

where

$\displaystyle k_{3}(\rho)$ $\displaystyle=K_{0}e^{\rho},$ (30) $\displaystyle
d\Omega^{2}$ $\displaystyle=d\theta^{2}+\sin^{2}\theta d\phi^{2}.$ (31)

In terms of $\rho$, we then have

$\begin{split}\phi&=\frac{1}{\sqrt{k_{2}}},\quad\mathcal{A}=\frac{k_{1,\rho}}{2k_{1}\sqrt{k_{2}}},\\\
\mathcal{K}&=\frac{k_{2}}{K_{0}e^{\rho}},\quad
Y=\frac{k_{1,\rho}}{2k_{1}},\end{split}$ (32)

where $K_{0}$ is a suitable constant. With the above relation, it is not
difficult to express the newly obtained solutions in terms of the area radius,
$r$. The relation between $r$ and $\rho$ is

$\displaystyle\rho$ $\displaystyle=2\ln\left(\frac{r}{r_{0}}\right),$ (33)
$\displaystyle r$ $\displaystyle=\sqrt{K_{0}}e^{\rho/2},$ (34)

where $r_{0}$ and $K_{0}$ are constants related by $K_{0}=r_{0}^{2}$.

The conversion from the metric coefficients in the $\rho$ coordinate to the
$r$ coordinate can be achieved by noting that $k_{1}$ and $k_{3}$ are scalars
with respect to a change of the radial parameter and that

$k_{2}(\rho)\rightarrow\frac{r^{2}}{4}k_{2}(r).$ (35)

Therefore, in terms of the area radius, we have

$ds^{2}=-k_{1}(r)dt^{2}+k_{2}(r)dr^{2}+r^{2}d\Omega^{2},$ (36)

and

$\begin{split}\phi&=\frac{2}{r\sqrt{k_{2}}},\quad\mathcal{A}=\frac{k_{1,r}}{2k_{1}\sqrt{k_{2}}},\\\
\mathcal{K}&=\frac{k_{2}}{4},\quad Y=\frac{k_{1,r}}{4k_{1}}.\end{split}$ (37)

These formulae will be useful to write the results of the TOV equations in a
form closer to the usual one in literature.

##  III Conditions for physical viability

Though many solutions to the TOV equations exist, several of these are not
physical. In this section, we provide a set of conditions that a solution must
satisfy in order to be physically relevant [3, 28, 25, 26].

We require:

  1. (i)

the sources of the Einstein equations are positive definite

$\mu\geq 0,\quad p+\Pi\geq 0,\quad p-\frac{1}{2}\Pi\geq 0,$ (38)

  2. (ii)

both fluids (individually and combined) satisfy the weak energy condition

$\mu\geq 0,\quad\mu+p+\Pi\geq 0,\quad\mu+p-\frac{1}{2}\Pi\geq 0,$ (39)

  3. (iii)

each fluid satisfies the conditions

$\mu^{\prime}<0,\quad p^{\prime}+\Pi^{\prime}<0,$ (40)

where the prime represents the derivative with respect to the area radius.
These conditions are necessary (but not sufficient) for the stability of the
solution. If $p$ is a decreasing, and $\mu$ is increasing, we have that
$\frac{dp}{d\mu}$ is negative. The requirement of a negative gradient for the
pressure and density is equivalent to the requirement that the speed of sound
is real.

  4. (iv)

for both fluids causality, is preserved and therefore the speed of sound must
never be higher than the speed of light:

$\begin{split}0<c_{s,r}^{2}&+c_{s,\perp}^{2}\leq 1,\end{split}$ (41)

where

$\begin{split}c_{s,r}^{2}&=\frac{\partial p_{r}}{\partial\mu},\\\
c_{s,\perp}^{2}&=\frac{\partial p_{\perp}}{\partial\mu},\end{split}$ (42)

  5. (v)

the 1+1+2 potentials and therefore, the metric coefficients and some of their
first derivatives are finite and positive valued at the center of the matter
distribution.

Notice that these conditions do not limit the sign of the anisotropy $\Pi$. In
fact, combining (39) we obtain

$-p\leq\Pi\leq 2p,$ (43)

which shows that the anisotropic pressure can be negative without violating
the weak energy condition.

In the case of a complete solution, we consider the additional requirement
that the total radial pressure is decreasing. Instead, the total anisotropic
pressure can have any behaviour compatible with the conditions above.

##  IV Junction conditions

An important aspect in the search for solutions is the problem of the junction
between the interior solution and the exterior vacuum spacetime, which is
characterised by the Schwarzschild metric. The procedure to determine the
junction conditions is similar to the single fluid case treated in [25, 26].
In particular, Israel’s junction conditions [5, 29] are equivalent to

$[{\cal K}]=0,\qquad[Y]=0\,.$ (44)

Using the constraint Eq. (12) in (44), we obtain

$\left[P_{1}+\mathbb{P}_{1}+P_{2}+\mathbb{P}_{2}\right]=0\,.$ (45)

This implies that in order to provide a smooth junction with the Schwarzschild
metric, the sum of the radial pressure for both fluids must be zero at the
junction. As the radial pressures are positive definite in a realistic
solution, this result implies that the radial pressure of both fluids must be
zero at the junction. A possible occurrence in a two-fluid system is that the
radial pressure of one fluid goes to zero before the other. In this case, we
will consider our solution as valid only when both pressures are positive,
regardless of what happens afterwards. Such solutions will be referred to as
representing a shell in a compact stellar object222It should be noted that the
vanishing of the radial pressure on the boundary is more of a desirable
feature than a physical necessity. In fact, there is no requirement for a
“hard boundary” in a stellar object, provided that the total mass remains
finite. Nevertheless, in the following we will always require that our
solutions present a hard boundary..

##  V Some known single fluid solutions of the TOV equations

In the sections that follow, we will consider several direct and indirect
resolution methods of the TOV equations (18). These methods rely on known
single fluid solutions. In this section, we introduce four solutions in a form
compatible with the formalism we employ: two anisotropic solutions, the
Bowers-Liang [30] and an element in the class of Florides [31] solutions, as
well as two isotropic solutions: the interior Schwarzschild solution [6], and
the Tolman IV solution [1].

###  V.1 Bowers-Liang solution

The Bowers-Liang solution, discovered in 1974 [30] is the extension of the
interior Schwarzschild solution [6] to include anisotropy. In this subsection,
we present the single fluid Bowers-Liang solution.

For a metric written in the form (36), the Bowers-Liang solution is given as

$\begin{split}k_{1}&=A_{0}\left[3^{\frac{h}{2}+1}\mathfrak{c}_{1}+\left(3-\mu_{0}r^{2}\right)^{\frac{h}{2}}\right]^{\frac{2}{h}},\\\
k_{2}&=\frac{3}{3-\mu_{0}r^{2}},\\\ k_{3}&=r^{2},\end{split}$ (46)

where $A_{0}$, ${\mathfrak{c}}_{1}$, $h$ and $\mu_{0}$ are constants. The
metric represented by (46) corresponds, via the Einstein equations, to the
following expressions for the radial pressure, tangential pressure and energy
density respectively

$\displaystyle{p_{r}}^{BL}(r)$
$\displaystyle=-\frac{\mu_{0}\left(3^{-\frac{h}{2}}\mathfrak{z}_{0}^{\frac{h}{2}}+\mathfrak{c}_{1}\right)}{3^{-\frac{h}{2}}\mathfrak{z}_{0}^{\frac{h}{2}}+3\mathfrak{c}_{1}},$
(47a) $\displaystyle{p_{\perp}}^{BL}(r)$
$\displaystyle=-\frac{\mu_{0}}{\mathfrak{z}_{0}\left(3^{\frac{h}{2}+1}\mathfrak{c}_{1}+\mathfrak{z}_{0}^{\frac{h}{2}}\right)^{2}}\left(\mathfrak{z}_{0}^{h+1}\right.$
$\displaystyle\left.+\mathfrak{c}_{1}\left\\{3^{h+1}\mathfrak{c}_{1}\mathfrak{z}_{0}+3^{1-\frac{h}{2}}\mathfrak{z}_{0}^{\frac{h}{2}}\left[12-(h+3)\mu_{0}r^{2}\right]\right\\}\right),$
(47b) $\displaystyle\mu_{BL}(r)$ $\displaystyle=\mu_{0},$ (47c)

where

$\mathfrak{z}_{0}=3-r^{2}\mu_{0}.$ (48)

In terms of the newly defined variables (9) and the parameter $\rho$, we have

$\phi_{BL}=2\sqrt{\frac{z}{3e^{\rho}K_{0}}},$ (49)

and the system (47) corresponds to

$\displaystyle\mathcal{K}_{BL}(\rho)=$ $\displaystyle\frac{3}{4z_{0}},$ (50a)
$\displaystyle P_{BL}(\rho)=$
$\displaystyle\frac{e^{\rho}K_{0}\mu_{0}}{4z_{0}^{2}\left(3^{-\frac{h}{2}}z_{0}^{\frac{h}{2}}+3\mathfrak{c}_{1}\right)^{2}}\left\\{-3^{1-h}z_{0}^{h+1}-9\mathfrak{c}_{1}^{2}z_{0}\right.$
$\displaystyle\left.+2\mathfrak{c}_{1}3^{-\frac{h}{2}}z_{0}^{\frac{h}{2}}\left[(h+5)e^{\rho}K_{0}\mu_{0}-18\right]\right\\},$
(50b) $\displaystyle\mathbb{P}_{BL}(\rho)=$
$\displaystyle-\frac{3^{-\frac{h}{2}}(h-1)e^{2\rho}K_{0}^{2}\mu_{0}^{2}\mathfrak{c}_{1}z_{0}^{\frac{h}{2}-2}}{2\left(3^{-\frac{h}{2}}z_{0}^{\frac{h}{2}}+3\mathfrak{c}_{1}\right)^{2}},$
(50c) $\displaystyle\mathbb{M}_{BL}(\rho)=$ $\displaystyle
e^{\rho}K_{0}\mu_{0}\left(\frac{3}{4z_{0}}\right),$ (50d) $\displaystyle
Y_{BL}(\rho)=$
$\displaystyle-\frac{e^{\rho}K_{0}\mu_{0}z_{0}^{\frac{h}{2}-1}}{2\left(3^{\frac{h}{2}+1}\mathfrak{c}_{1}+z_{0}^{\frac{h}{2}}\right)},$
(50e)

where

$z=3-e^{\rho}K_{0}\mu_{0}.$ (51)

###  V.2 A Florides-class solution

The Florides solution, discovered in 1974, describes a static, spherically
symmetric distribution of dust, that has vanishing radial pressure but nonzero
tangential pressure [31]. It represents the simplest element in a class of
solutions of the TOV equations for which

$\mathbb{P}=\frac{1}{6}\mathbb{M}\left(1-4\mathcal{K}\right),$ (52)

holds [26]. In this subsection, we present an example that belongs to the
class of single fluid Florides solutions characterized by a constant density.

This solution can be written in the form (36) setting

$\begin{split}k_{1}&=\frac{B_{0}\left(e^{A_{0}}+\frac{r^{2}}{r_{0}^{2}}\right)^{2}}{\sqrt{3-\mu_{0}r^{2}}},\\\
k_{2}&=\frac{3}{3-\mu_{0}r^{2}},\\\ k_{3}&=r^{2},\end{split}$ (53)

where $A_{0}$, $B_{0}$ and $\mu_{0}$ are constants. As mentioned previously,
$r_{0}$ is a constant. The metric represented by (53) corresponds, via the
Einstein equations, to the following expressions for the radial pressure,
tangential pressure and energy density respectively

$\displaystyle{p_{r}}^{F}(r)$
$\displaystyle=-\frac{4\left(\mu_{0}r^{2}-3\right)}{3\left(e^{A_{0}}r_{0}^{2}+r^{2}\right)},$
(54a) $\displaystyle{p_{\perp}}^{F}(r)$
$\displaystyle=\frac{\mu_{0}^{2}r^{2}}{4\left(3-\mu_{0}r^{2}\right)}-{p_{r}}^{F},$
(54b) $\displaystyle\mu_{F}(r)$ $\displaystyle=\mu_{0}.$ (54c)

In terms of the newly defined variables (9) and the parameter $\rho$, we have

$\phi_{F}=-2\sqrt{\frac{3-e^{\rho}K_{0}\mu_{0}}{3e^{\rho}K_{0}}},$ (55)

and the system (54) corresponds to

$\displaystyle\mathcal{K}_{F}(\rho)=$ $\displaystyle\frac{3}{4z},$ (56a)
$\displaystyle P_{F}(\rho)=$
$\displaystyle\frac{e^{2\rho}K_{0}^{2}\mu_{0}^{2}}{8z^{2}}+\frac{1}{1+e^{-\rho+\mathfrak{c}_{1}}},$
(56b) $\displaystyle\mathbb{P}_{F}(\rho)=$
$\displaystyle-\frac{e^{2\rho}K_{0}^{2}\mu_{0}^{2}}{8z^{2}},$ (56c)
$\displaystyle\mathbb{M}_{F}(\rho)=$ $\displaystyle
e^{\rho}K_{0}\mu_{0}\left(\frac{3}{4z}\right),$ (56d) $\displaystyle
Y_{F}(\rho)=$ $\displaystyle
e^{\rho}\left(\frac{K_{0}\mu_{0}}{4z}+\frac{1}{e^{\rho}+e^{\mathfrak{c}_{1}}}\right),$
(56e)

where $\mathfrak{c}_{1}$ is a constant and we have used (51).

###  V.3 The interior Schwarzschild solution

The interior Schwarzschild solution [6], the first solution for the interior
of a static spherically symmetric relativistic object can be written in the
form (36) by choosing

$\begin{split}k_{1}&=a_{0}\left(c_{1}+\sqrt{3-r^{2}\mu_{1}}\right)^{2},\\\
k_{2}&=\frac{3}{3-r^{2}\mu_{1}},\\\ k_{3}&=r^{2},\end{split}$ (57)

where $a_{0}$, $c_{1}$ and $\mu_{1}$ are constants. For the metric described
by (57), we have the following expressions for the pressure and energy
density, respectively

$\displaystyle p_{IS}(r)$
$\displaystyle=\frac{\mu_{1}\left(c_{1}\sqrt{3-\mu_{1}r^{2}}-3\mu_{1}r^{2}+9\right)}{3\mu_{1}r^{2}-9-3c_{1}\sqrt{3-\mu_{1}r^{2}}},$
(58a) $\displaystyle\mu_{IS}(r)$ $\displaystyle=\mu_{1}.$ (58b)

In terms of the parameter $\rho$, we have

$\phi_{IS}=\frac{2e^{-\frac{\rho}{2}}\sqrt{3-K_{0}\mu_{1}e^{\rho}}}{\sqrt{3}\sqrt{K_{0}}},$
(59)

and (58) corresponds to

$\displaystyle\mathcal{K}_{IS}(\rho)=$
$\displaystyle\frac{3}{12-4K_{0}\mu_{1}e^{\rho}},$ (60a) $\displaystyle
P_{IS}(\rho)=$
$\displaystyle\frac{K_{0}\mu_{1}e^{\rho}\left(c_{1}^{2}+2c_{1}\sqrt{3-K_{0}\mu_{1}e^{\rho}}+3K_{0}\mu_{1}e^{\rho}-9\right)}{4\left(K_{0}\mu_{1}e^{\rho}-3\right)\left(c_{1}^{2}+K_{0}\mu_{1}e^{\rho}-3\right)},$
(60b) $\displaystyle\mathbb{M}_{IS}(\rho)=$
$\displaystyle\frac{3\mu_{1}e^{\rho}}{12K_{0}-4\mu_{1}e^{\rho}},$ (60c)
$\displaystyle Y_{IS}(\rho)=$
$\displaystyle\frac{K_{0}\mu_{1}e^{\rho}}{-2c_{1}\sqrt{3-K_{0}\mu_{1}e^{\rho}}+2K_{0}\mu_{1}e^{\rho}-6}.$
(60d)

###  V.4 The Tolman IV solution

The Tolman IV solution presented in 1939 [1], is a well-known solution with no
singularity at $r=0$. The Tolman IV solution corresponds to the choice

$\begin{split}k_{1}&=B^{2}\left(1+\frac{r^{2}}{A^{2}}\right),\\\
k_{2}&=\frac{R^{2}(A^{2}+2r^{2})}{\left(R^{2}-r^{2}\right)\left(A^{2}+r^{2}\right)},\\\
k_{3}&=r^{2},\end{split}$ (61)

in the metric (36), where $A$, $B$ and $R$ are constants. For the metric
described by (61), the pressure and energy density are given by

$\displaystyle p_{T}(r)$
$\displaystyle=\frac{R^{2}-A^{2}-3r^{2}}{R^{2}\left(A^{2}+2r^{2}\right)},$
(62a) $\displaystyle\mu_{T}(r)$
$\displaystyle=\frac{R^{2}\left(3A^{2}+2r^{2}\right)+7A^{2}r^{2}+3A^{4}+6r^{4}}{R^{2}\left(A^{2}+2r^{2}\right)^{2}}.$
(62b)

In terms of the parameter $\rho$, we have

$\phi_{T}=-\sqrt{\frac{4(A^{2}+K_{0}e^{\rho})(R^{2}-K_{0}e^{\rho})}{K_{0}e^{\rho}R^{2}(A^{2}+2K_{0}e^{\rho})}},$
(63)

and (62) corresponds to

$\displaystyle\mathcal{K}_{T}(\rho)=$
$\displaystyle\frac{R^{2}\left(A^{2}+2K_{0}e^{\rho}\right)}{4\left(A^{2}+K_{0}e^{\rho}\right)\left(R^{2}-K_{0}e^{\rho}\right)},$
(64a) $\displaystyle P_{T}(\rho)=$
$\displaystyle\frac{K_{0}e^{\rho}\left(A^{2}+3K_{0}e^{\rho}-R^{2}\right)}{4\left(A^{2}+K_{0}e^{\rho}\right)\left(K_{0}e^{\rho}-R^{2}\right)},$
(64b) $\displaystyle\mathbb{M}_{T}(\rho)=$
$\displaystyle\frac{A^{2}\left(2A^{2}+R^{2}\right)}{4\left(A^{2}+R^{2}\right)\left(A^{2}+K_{0}e^{\rho}\right)}-\frac{A^{2}}{2\left(A^{2}+2K_{0}e^{\rho}\right)}$
$\displaystyle+\frac{R^{2}\left(3A^{2}+4R^{2}\right)}{4\left(A^{2}+R^{2}\right)\left(R^{2}-K_{0}e^{\rho}\right)}-\frac{3}{4},$
(64c) $\displaystyle Y_{T}(\rho)=$
$\displaystyle\frac{K_{0}e^{\rho}}{2\left(A^{2}+K_{0}e^{\rho}\right)}.$ (64d)

##  VI A direct resolution strategy

In this section, we present a direct method to obtain two-fluid solutions from
the single fluid ones given in Sec. V. We assume that the total isotropic and
anisotropic pressures, density, and metric coefficients are those of known
single fluid solutions given in Sec. V. Therefore, the spacetime of the
objects described by these new solutions is that of the corresponding known
single fluid solution. Using some additional conditions, we can then deduce
the thermodynamical properties of a set of two fluids, thus obtaining a
complete solution.

In the following, we will present two-fluid solutions which satisfy the
conditions in Sec. III for at least one set of parameters. Note that, for
simplicity sake, we set $K_{0}$ and $r_{0}$ to unity from this point forward.

###  VI.1 A two-fluid Bowers-Liang solution

A first example is based on the idea behind the derivation of the Bowers-Liang
solution [30]. This is a widely used approach to find single fluid solutions
[30, 8], which involves making an ansatz on the behaviour of the anisotropy
and solving for the corresponding pressure.

We choose the total anisotropic pressure $\mathbb{P}_{tot}$ as the Bowers-
Liang anisotropic pressure for a single fluid (50c):

$\mathbb{P}_{tot}=-\frac{e^{2\rho}\mu_{tot}^{2}}{2z_{1}^{2}}\left[3^{\frac{h}{2}}(h-1)\mathfrak{c}z^{\frac{h}{2}-2}\right],$
(65)

where $\mathfrak{c}$ is a generic constant, and we have used (51) and (157).
In this case, since we are extending the Bowers-Liang solution, $\mu_{1}$ and
$\mu_{2}$ are constant. We can then write, from the definition of
$\mathcal{K}$ and $\mathbb{M}_{i}$

$\displaystyle\mathbb{M}_{1}$ $\displaystyle=e^{\rho}\mu_{1}\mathcal{K},$
(66a) $\displaystyle\mathbb{M}_{2}$
$\displaystyle=e^{\rho}\mu_{2}\mathcal{K}.$ (66b)

It is easy to check that with the above choices (20) gives

$\mathcal{K}=\frac{3}{4z}=\mathcal{K}_{BL},$ (67)

and substituting these results in (22) we have

$P_{tot}=\frac{e^{\rho}\mu_{tot}}{4z^{2}z_{1}^{2}}\left[z_{2}-3z\left(3^{h+1}\mathfrak{c}^{2}+z^{h}\right)\right]=P_{BL},$
(68)

as expected, where

$z_{2}=2\times
3^{h/2}\mathfrak{c}z^{h/2}\left((h+5)e^{\rho}\mu_{tot}-18\right).$ (69)

To solve for the individual fluids we write

$\displaystyle P_{tot}=P_{1}+P_{2},$ (70a)
$\displaystyle\mathbb{P}_{tot}=\mathbb{P}_{1}+\mathbb{P}_{2}.$ (70b)

Let us now assume $\mathbb{P}_{1}=0$. Then $\mathbb{P}_{2}=\mathbb{P}_{tot}$
and we can replace $P_{2}$ in the TOV equation for fluid 1 (18) with
$P_{tot}-P_{2}$. With $\mathbb{M}_{1}$ from (66a), $\mathcal{K}$ from (67) and
$\mathbb{P}_{tot}$ from (65), we can solve the resulting TOV equation for
$P_{1}$. Subtraction from the total isotropic pressure $P$ gives $P_{2}$.

Using $\phi$ from (49) and (9), as well as (5), we can rewrite $P$ and
$\mathbb{P}$ for each fluid in terms of $r$ and derive the radial and
tangential pressures as

$\displaystyle p_{r1}=$ $\displaystyle p_{\perp
1}=-\mu_{1}+\frac{4}{3}c_{1}\mathfrak{z}_{1}^{-\frac{1}{h}},$ (71a)
$\displaystyle p_{r2}=$ $\displaystyle\frac{2\times
3^{\frac{h}{2}}\mu_{tot}\mathfrak{c}}{3^{\frac{h}{2}+1}\mathfrak{c}+\mathfrak{z}^{\frac{h}{2}}}-\frac{4}{3}c_{1}\mathfrak{z}_{1}^{-\frac{1}{h}}-\mu_{2},$
(71b) $\displaystyle p_{\perp 2}=$
$\displaystyle\frac{\mathfrak{z}\left(-3^{h+2}\mu_{tot}\mathfrak{c}^{2}-4c_{1}\mathfrak{z}_{1}^{2-\frac{1}{h}}+3\mu_{1}\mathfrak{z}_{1}^{2}\right)}{3\mathfrak{z}\mathfrak{z}_{1}^{2}}$
$\displaystyle+\frac{3^{\frac{h}{2}+1}\mu_{tot}\mathfrak{c}\mathfrak{z}^{\frac{h}{2}}\left[(h+3)\mu_{tot}r^{2}-12\right]}{3\mathfrak{z}\mathfrak{z}_{1}^{2}}$
$\displaystyle-\frac{3\mu_{tot}\mathfrak{z}^{h+1}}{3\mathfrak{z}\mathfrak{z}_{1}^{2}},$
(71c)

where

$\displaystyle\mathfrak{z}=$ $\displaystyle 3-r^{2}(\mu_{1}+\mu_{2}),$ (72a)
$\displaystyle\mathfrak{z}_{1}=$ $\displaystyle
3^{\frac{h}{2}+1}\mathfrak{c}+\mathfrak{z}^{\frac{h}{2}}.$ (72b)

The above solution is plotted in Figs. 1, 2 and 3 for three different values
of $h$. Figure 1 shows the case $h=1$ for which both fluids are isotropic. For
the same set of parameters we can plot solution (71) for $h=\frac{1}{2}$. It
is evident from Fig. 2 that in this case, the solution is not physical, as the
radial pressure of fluid 1 is negative. This result shows that not all values
of $h$ are compatible with a given set of parameters. If we take instead
$h=0.8$ (see Fig. 3), the solution is again physical and we can estimate the
effect that the presence of anisotropy has on the thermodynamical quantities.
In fact, Fig. 4 shows that the anisotropy lowers the radial pressure at the
core and raises the radial pressure at the surface, suggesting that
anisotropic compact stellar objects can have different structures with respect
to their isotropic counterparts. This is consistent with the results of [10]
already mentioned in the introduction.

Finally, as in (71) the radial pressure of fluid 1 goes to zero before that of
fluid 2, this solution can only represent a shell and will have to be joined
to another solution to describe a complete compact stellar object.

Figure 1: Radial pressure for the Bowers-Liang solution in (71) for fluid 1
(blue) and fluid 2 (green) and the tangential pressure for fluid 2 (red) vs
$r$ with $h=1$, $\mu_{1}=0.68$, $\mu_{2}=2.33$, $c_{1}=-0.45$ and
$\mathfrak{c}=-0.4$. The red and green curves overlap. Figure 2: Radial
pressure for the Bowers-Liang solution in (71) for fluid 1 (blue) and fluid 2
(green) and the tangential pressure for fluid 2 (red) vs $r$ with
$h=\frac{1}{2}$, $\mu_{1}=0.68$, $\mu_{2}=2.33$, $c_{1}=-0.45$ and
$\mathfrak{c}=-0.4$. Figure 3: Radial pressure for the Bowers-Liang solution
in (71) for fluid 1 (blue) and fluid 2 (green) and the tangential pressure for
fluid 2 (red) vs $r$ with $h=0.8$, $\mu_{1}=0.68$, $\mu_{2}=2.33$,
$c_{1}=-0.45$ and $\mathfrak{c}=-0.4$. Figure 4: Radial pressure for the
Bowers-Liang solution in (71) for fluid 2 with $h=0.8$ (solid green) and $h=1$
(dashed green) and the tangential pressure for fluid 2 (red) vs $r$ with
$\mu_{1}=0.68$, $\mu_{2}=2.33$, $c_{1}=-0.45$ and $\mathfrak{c}=-0.4$, zoomed
in for $r=0$ to $r=0.1$.

###  VI.2 Two Florides-class solutions

We utilize the single fluid Florides-class constant density solution in Sec.
V.2. For these solutions, the relation (52) always holds. For notational
convenience, we will write this relation as

$\mathbb{P}_{tot}=\frac{1}{6}\mathbb{M}_{tot}(1-4\mathcal{K}_{F}).$ (73)

Following [26], using the assumption above and defining the quantity

$\mathcal{P}_{tot}=P_{tot}+\mathbb{P}_{tot},$ (74)

related with the radial pressure, one can write the total TOV equation as

$\mathcal{P}_{tot,\rho}+\mathcal{P}_{tot}^{2}+\mathcal{P}_{tot}\left(3\mathcal{K}-\mathbb{M}_{1}-\mathbb{M}_{2}-\frac{7}{4}\right)=0,$
(75)

where $\mathcal{K}=\mathcal{K}_{F}$ is given by (56a) [or can derived by (20),
remembering that in this case $\mathbb{M}_{1}$ and $\mathbb{M}_{2}$ are given
by (66)]. The solution of Eq. (75) reads

$\mathcal{P}_{tot}=\frac{1}{e^{\mathfrak{c}-\rho}+1}.$ (76)

Now Eq. (74) can be used to determine $P_{tot}$ which reads

$P_{tot}=\frac{e^{2\rho}\mu_{tot}^{2}}{8z^{2}}+\frac{1}{e^{\mathfrak{c}-\rho}+1},$
(77)

and (73) gives

$\mathbb{P}_{tot}=-\frac{e^{2\rho}\mu_{tot}^{2}}{8z^{2}}.$ (78)

In order to find the pressures of the two fluids of the new solution we can
use the relations

$\displaystyle\mathcal{P}_{i}$ $\displaystyle=$ $\displaystyle
P_{i}+\mathbb{P}_{i},\qquad i=1,2$ (79) $\displaystyle P_{2}$ $\displaystyle=$
$\displaystyle P_{tot}-P_{1},$ (80) $\displaystyle\mathcal{P}_{2}$
$\displaystyle=$ $\displaystyle\mathcal{P}_{tot}-\mathcal{P}_{1},$ (81)

in the TOV equation for fluid 1 (18) so that it will contain only
$\mathcal{P}_{1,\rho}$, $\mathbb{P}_{1}$, $\mathcal{P}_{tot}$,
$\mathbb{M}_{1}$ and $\mathbb{M}_{2}$:

$\displaystyle\mathcal{P}_{1,\rho}=$
$\displaystyle\mathcal{P}_{1}\left[\frac{7}{4}-3\mathcal{K}+2-\mathcal{P}_{tot}+2\mathbb{M}_{tot}\right]$
(82)
$\displaystyle+\frac{1}{4}\mathbb{M}_{1}(-4\mathcal{K}-4\mathcal{P}_{tot}+1)-\frac{3}{2}\mathbb{P}_{1},$

where $\mathbb{M}_{tot}=\mathbb{M}_{1}+\mathbb{M}_{2}$ and
$\mathcal{K}=\mathcal{K}_{F}$.

We now consider two cases: The first, in which one of the fluids is isotropic
$(\mathbb{P}_{1}=0)$ and a second, in which both fluids are anisotropic. For
convenience we will suppose, without loss of generality, that fluid 1 changes
its character.

#### The case $\mathbb{P}_{1}=0$

We take $\mathbb{P}_{1}=0$, so that $\mathcal{P}_{1}={P}_{1}$ and
$\mathbb{P}_{2}=\mathbb{P}_{tot}$, and the TOV equation for the first fluid
(82) becomes

$\displaystyle\mathcal{P}_{1,\rho}=$
$\displaystyle\mathcal{P}_{1}\left[\frac{7}{4}-3\mathcal{K}+2-\mathcal{P}_{tot}+2\mathbb{M}_{tot}\right]$
(83)
$\displaystyle+\frac{1}{4}\mathbb{M}_{1}(-4\mathcal{K}-4\mathcal{P}_{tot}+1).$

We can solve this equation for $\mathcal{P}_{1}$ with (66), (67) and (76), and
obtain $\mathcal{P}_{2}$ by subtracting from the total $\mathcal{P}_{tot}$.
Since $P_{1}=\mathcal{P}_{1}$, we obtain $P_{2}$ by subtracting
$\mathcal{P}_{1}$ from $P_{tot}$. As $\mathbb{P}_{2}=\mathbb{P}_{tot}$, all
the pressures are known in terms of $\rho$ and we can utilize $\phi$ from (55)
with (9) to give $p_{1},p_{2},\Pi_{1},$ and $\Pi_{2}$ in terms of $r$.
Finally, using (5), we can write the radial and tangential pressures for each
fluid in terms of $r$ as

$\displaystyle p_{r1}=$ $\displaystyle p_{\perp
1}=-\mu_{1}+\frac{4c_{1}\mathfrak{z}^{\frac{1}{4}}}{3\left(e^{\mathfrak{c}}+r^{2}\right)},$
(84a) $\displaystyle p_{r2}=$
$\displaystyle\frac{3\mu_{1}e^{\mathfrak{c}}-r^{2}(\mu_{1}+4\mu_{2})-4c_{1}\mathfrak{z}^{\frac{1}{4}}+12}{3\left(e^{\mathfrak{c}}+r^{2}\right)},$
(84b) $\displaystyle p_{\perp 2}=$
$\displaystyle\frac{9\mu_{1}-3\mu_{2}+3\mu_{tot}+\mu_{tot}r^{2}(\mu_{2}-3\mu_{1})}{12-4\mu_{tot}r^{2}}$
$\displaystyle-\frac{4\left(\mu_{tot}r^{2}+c_{1}\mathfrak{z}^{\frac{1}{4}}-3\right)}{3\left(e^{\mathfrak{c}}+r^{2}\right)},$
(84c)

where we have used the $\mathfrak{z}_{i}$ from (72). This solution is plotted
in Fig. 5, and represents a shell as the radial pressure of fluid 2 approaches
zero before that of fluid 1.

Figure 5: Radial pressure for the constant density Florides-class solution in
(84) for fluid 1 (blue) and fluid 2 (green), and tangential pressure for fluid
2 (red) vs $r$ with $\mu_{1}=0.1$, $\mu_{2}=1.55$, $c_{1}=1.35$ and
$\mathfrak{c}=-1.85$.

#### The case $\mathbb{P}_{1}\neq 0$

When fluid 1 is anisotropic, we need an ansatz for $\mathbb{P}_{1}$ in order
to solve (82). One convenient choice is

$\mathbb{P}_{1}=\frac{1}{6}\mathbb{M}_{1}(1-4\mathcal{K}_{F}).$ (85)

In this way, equation (82) takes the form

$\mathcal{P}_{1,\rho}=\mathcal{P}_{1}\left(-3\mathcal{K}_{F}+2\mathbb{M}_{tot}-\mathcal{P}_{tot}+\frac{7}{4}\right)-\mathbb{M}_{1}\mathcal{P}_{tot}.$
(86)

As before, once we know $\mathcal{P}_{1}$ we can obtain $\mathcal{P}_{2}$ by
subtraction from $\mathcal{P}_{tot}$ and ${P}_{1}$ by subtraction of
$\mathbb{P}_{1}$ in (85) from $\mathcal{P}_{1}$. Finally $\mathcal{P}_{2}$,
$\mathbb{P}_{2}$ and $P_{2}$ are derived from the equations for the total
pressures.

Once again, we utilize $\phi$ from (55), (9) to give $p_{1},p_{2},\Pi_{1},$
and $\Pi_{2}$ in terms of $r$. With equation (5), we can write the radial and
tangential pressures for each fluid in terms of $r$ as

$\displaystyle p_{r1}=$
$\displaystyle\frac{4\mu_{1}\mathfrak{z}+4c_{2}\mu_{tot}\mathfrak{z}^{\frac{1}{4}}}{3\mu_{tot}\left(r^{2}+e^{c_{1}}\right)},$
(87a) $\displaystyle p_{\perp 1}=$
$\displaystyle\frac{\mu_{1}\mu_{tot}r^{2}}{4\mathfrak{z}}+\frac{4\left(\mu_{1}\mathfrak{z}+c_{2}\mu_{tot}\mathfrak{z}^{\frac{1}{4}}\right)}{3\mu_{tot}\left(r^{2}+e^{c_{1}}\right)},$
(87b) $\displaystyle p_{r2}=$
$\displaystyle\frac{4}{3}\left(\frac{\mathfrak{z}}{e^{\mathfrak{c}}+r^{2}}-\frac{\mu_{1}\mathfrak{z}+c_{2}\mu_{tot}\mathfrak{z}^{\frac{1}{4}}}{\mu_{tot}r^{2}+e^{c_{1}}\mu_{tot}}\right),$
(87c) $\displaystyle p_{\perp 2}=$
$\displaystyle\frac{4\mathfrak{z}}{3\left(e^{\mathfrak{c}}+r^{2}\right)}-\frac{4\mathfrak{z}\left(\frac{\mu_{1}}{\mu_{tot}}+\frac{c_{2}}{\mathfrak{z}^{3/4}}\right)}{3\left(r^{2}+e^{c_{1}}\right)}+\frac{\mu_{2}\mu_{tot}r^{2}}{4\mathfrak{z}},$
(87d)

where $c_{1}$, $c_{2}$ and $\mathfrak{c}$ are constants. This solution is
plotted in Figs. 6 and 7. We can see that for this solution, the radial
pressure of fluid 1 approaches zero before fluid 2, so this solution
represents a shell.

Figure 6: Radial pressure (blue) and tangential pressure (orange) for the
constant density Florides-class solution in (87) for fluid 1 vs $r$ with
$\mu_{1}=0.2$, $\mu_{2}=0.6$, $c_{1}=-1$, $c_{2}=-0.45$ and
$\mathfrak{c}=-0.05$. Figure 7: Radial pressure (green) and tangential
pressure (red) for the constant density Florides-class solution in (87) for
fluid 2 vs $r$ with $\mu_{1}=0.2$, $\mu_{2}=0.6$, $c_{1}=-1$, $c_{2}=-0.45$
and $\mathfrak{c}=-0.05$. We have truncated this graph to the value of $r$
where the radial pressure of fluid 1 approaches zero.

###  VI.3 Two Florides-class solutions with fractional density

We can consider a variation of the previous strategy assuming a form for the
pressure $\mathcal{P}_{tot}$ and the energy density. We choose the total
pressure variable $\mathcal{P}_{tot}$ as

$\mathcal{P}_{tot}=-\frac{3e^{\rho}\left(be^{\rho}-a\right)}{z_{3}},$ (88)

where $a$ and $b$ are constants and

$z_{3}=12a^{2}-3a(4b-1)e^{\rho}+2b(2b-1)e^{2\rho}.$ (89)

We can then derive $\mathbb{M}_{tot}$ from Eq. (75) and use this quantity to
find $\mathbb{P}_{tot}$ via (73) and solve for $\mathcal{K}$ via (20). At this
point, using (74) we can easily obtain $P_{tot}$. Utilizing this procedure, we
have

$\displaystyle P_{tot}=$ $\displaystyle\frac{e^{\rho}z_{6}}{z_{4}},$ (90a)
$\displaystyle\mathbb{P}_{tot}=$
$\displaystyle\frac{e^{2\rho}}{z_{4}}\left(4b+1\right)^{2}\left(2be^{\rho}-3a\right)z_{5},$
(90b)

where

$\displaystyle z_{4}$ $\displaystyle=24\left(a-be^{\rho}\right)z_{3}^{2},$
(91a) $\displaystyle z_{5}$
$\displaystyle=9a^{2}-7abe^{\rho}+2b^{2}e^{2\rho},$ (91b) $\displaystyle
z_{6}$ $\displaystyle=-3a^{2}b\left[8b(26b-107)+205\right]e^{2\rho}$
$\displaystyle+27a^{3}\left[8b(2b-11)+9\right]e^{\rho}+864a^{4}$
$\displaystyle+4ab^{2}\left[80(b-4)b+131\right]e^{3\rho}-4b^{3}\left[16(b-4)b\right.$
$\displaystyle\left.+37\right]e^{4\rho}.$ (91c)

We now need to give an ansatz for the energy density of the single fluids. In
fact, our choice of $\mathcal{P}_{tot}$ can only determine, via the equation
of state, the total density.

As an example, we will set $\mathbb{M}_{1}=\frac{1}{4}\mathbb{M}_{tot}$, such
that $\mathbb{M}_{2}=\frac{3}{4}\mathbb{M}_{tot}$, where

$\mathbb{M}_{tot}=-\frac{e^{\rho}z_{4}z_{5}}{6z_{3}}(4b+1),$ (92)

and

$\mathcal{K}=\frac{3\left(a-be^{\rho}\right)^{2}}{z_{3}}.$ (93)

As before, we consider two cases. The first in which one of the fluids is
isotropic $(\mathbb{P}_{1}=0)$ and a second in which both fluids are
anisotropic.

#### The case $\mathbb{P}_{1}=0$

In the case $\mathbb{P}_{1}=0$, we have $\mathbb{P}_{2}=\mathbb{P}_{tot}$ and
$\mathcal{P}_{1}=P_{1}$. With (90), (92) and (93), we can solve the resulting
TOV equation for fluid 1 (82) with (67), (76), and our ansatz on
$\mathbb{M}_{1}$. We can obtain $\mathcal{P}_{2}$ by subtracting from the
total $\mathcal{P}_{tot}$. Since $P_{1}=\mathcal{P}_{1}$, we obtain $P_{2}$ by
subtracting from $P_{tot}$. As $\mathbb{P}_{2}=\mathbb{P}_{tot}$, we proceed
to utilize $\phi$ from (55) with (9) to give $p_{1},p_{2},\Pi_{1},$ and
$\Pi_{2}$ in terms of $r$. Finally, using (5), we can write the radial and
tangential pressures for each fluid in terms of $r$ as

$\displaystyle p_{r1}=p_{\perp 1}=$
$\displaystyle\frac{1}{4a-4br^{2}}+\frac{1}{3}e^{\alpha}c_{1}\left[12a^{2}+3a(1-4b)r^{2}\right.$
$\displaystyle\left.+2b(2b-1)r^{4}\right]^{\frac{5-4b}{16b-8}},$ (94a)
$\displaystyle p_{r2}=$
$\displaystyle\frac{3}{4\left(a-br^{2}\right)}-\frac{1}{3}e^{\alpha}c_{1}\left[12a^{2}+3a(1-4b)r^{2}\right.$
$\displaystyle\left.+2b(2b-1)r^{4}\right]^{\frac{5-4b}{16b-8}},$ (94b)
$\displaystyle p_{\perp 2}=$
$\displaystyle\mathfrak{z}_{2}-\frac{1}{3}e^{\alpha}c_{1}\left[12a^{2}+3a(1-4b)r^{2}\right.$
$\displaystyle\left.+2b(2b-1)r^{4}\right]^{\frac{5-4b}{16b-8}},$ (94c)

and the energy densities of both fluids as

$\displaystyle\mu_{1}=$
$\displaystyle-\frac{(4b+1)\left(9a^{2}-7abr^{2}+2b^{2}r^{4}\right)}{48\left(a-br^{2}\right)^{3}},$
(95a) $\displaystyle\mu_{2}=$
$\displaystyle-\frac{(4b+1)\left(9a^{2}-7abr^{2}+2b^{2}r^{4}\right)}{16\left(a-br^{2}\right)^{3}},$
(95b)

where

$\displaystyle\alpha=$
$\displaystyle\frac{\sqrt{48b^{2}-24b-9}\tan^{-1}\left[\frac{a(3-12b)+4b(2b-1)r^{2}}{a\sqrt{48b^{2}-24b-9}}\right]}{8b(4b-5)+12},$
(96a) $\displaystyle\mathfrak{z}_{2}=$
$\displaystyle\left\\{-3a^{2}b\left[8b(26b-47)+109\right]r^{4}+27a^{3}\left[8b(2b-5)\right.\right.$
$\displaystyle\left.\left.+5\right]r^{2}+432a^{4}+16ab^{2}\left[5b(4b-7)+17\right]r^{6}\right.$
$\displaystyle\left.-4b^{3}\left[4b(4b-7)+19\right]r^{8}\right\\}\times\left\\{48\left(a-br^{2}\right)^{3}\right.$
$\displaystyle\left.\left[12a^{2}+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right]\right\\}^{-1}.$
(96b)

This solution, which only describes a shell, is plotted in Figs. 8 and 9.

Figure 8: Radial pressure for the Florides-class solution in (94) for fluid 1
(blue) and fluid 2 (green), tangential pressure for fluid 2 (red),
$dp_{r}/d\mu$ for fluid 1 (purple) and fluid 2 (black), and $dp_{\perp}/d\mu$
for fluid 2 (yellow) vs $r$ with $a=0.55$, $b=-1.3$, and $c_{1}=-1.3$. Figure
9: Energy density profile for the Florides-class solution in (95) for fluid 1
(gray) and fluid 2 (cyan) vs $r$ with $a=0.55$, $b=-1.3$, and $c_{1}=-1.3$.

#### The case $\mathbb{P}_{1}\neq 0$

We proceed as in the $\mathbb{P}_{1}=0$ case, and assume again (85) in (82) to
obtain the resulting TOV equation (86), which we can solve for
$\mathcal{P}_{1}$ with (66), (67), and (76). We have $\mathbb{P}_{1}$ from the
relation (85), and $P_{1}$ is given in (79). We can obtain $\mathcal{P}_{2}$,
$\mathbb{P}_{2}$ and $P_{2}$ by subtracting from the total pressures. Next, we
utilize $\phi$ from (55) with (9) to give $p_{1},p_{2},\Pi_{1},$ and $\Pi_{2}$
in terms of $r$. With equation (5), we can write the radial and tangential
pressures for each fluid in terms of $r$. The quantities $p_{r1},p_{\perp
1},p_{r2},$ and $p_{\perp 2}$ are too long to be presented here, and are given
by (152) in Appendix C.

The energy densities of both fluids are

$\displaystyle\mu_{1}=$
$\displaystyle-\frac{(4b+1)\left(9a^{2}-7abr^{2}+2b^{2}r^{4}\right)}{48\left(a-br^{2}\right)^{3}},$
(97a) $\displaystyle\mu_{2}=$
$\displaystyle-\frac{(4b+1)\left(9a^{2}-7abr^{2}+2b^{2}r^{4}\right)}{16\left(a-br^{2}\right)^{3}},$
(97b)

and are shown in Fig. 13. The solution represented by (152) is plotted in Fig.
10, 11 and 12. We observe that the radial pressure of fluid 2 approaches zero
before that of fluid 1. Thus, this solution represents a shell.

Figure 10: Radial pressure (blue) and tangential pressure (orange) for the
Florides-class solution in (152) for fluid 1 vs $r$ with $a=1.75$, $b=-1.1$,
and $c_{1}=3.35$. As the radial pressure for fluid 2 approaches zero before
fluid 1, we have truncated the graph for fluid 1. Figure 11: Radial pressure
(green) and tangential pressure (red) for the Florides-class solution in (152)
for fluid 2 vs $r$ with $a=1.75$, $b=-1.1$, and $c_{1}=3.35$. Figure 12:
$dp_{r}/d\mu$ for fluid 1 (purple) and fluid 2 (brown), and $dp_{\perp}/d\mu$
for fluid 1 (black) and fluid 2 (yellow) for the Florides-class solution in
(152) for fluid 1 vs $r$ with $a=1.75$, $b=-1.1$, and $c_{1}=3.35$. Figure 13:
Energy density profile for the Florides-class solution in (97) for fluid 1
(gray) and fluid 2 (cyan) vs $r$ with $a=1.75$, $b=-1.1$, and $c_{1}=3.35$.

###  VI.4 Two Florides-class solutions with barotropic equation of state

We follow the same procedure as in Sec. VI.3, with a different ansatz. Instead
of setting $\mathbb{M}_{1}$ and $\mathbb{M}_{2}$ as fractions of
$\mathbb{M}_{tot}$, we choose the equation of state. In particular we assume

$\mathbb{M}_{1}=\gamma\mathcal{P}_{1}.$ (98)

We now consider the usual two cases: one in which fluid 1 is isotropic
$(\mathbb{P}_{1}=0)$ and one in which both fluids are anisotropic.

#### The case $\mathbb{P}_{1}=0$

If we set the anisotropic pressure of one of the fluids to zero,
$\mathbb{P}_{2}=\mathbb{P}_{tot}$ and the equation of state (98) takes the
form

$\mathbb{M}_{1}=\gamma{P}_{1}.$ (99)

We can now solve the resulting TOV equation for fluid 1 (82) with (90), (92),
(93) and (98). Subtraction from the total pressures determines the pressures
for fluid 2.

As in the previous cases, we utilize $\phi$ from (55) with (9) to give
$p_{1},p_{2},\Pi_{1},$ and $\Pi_{2}$ in terms of $r$. Finally, using (5), we
can write the radial and tangential pressures for each fluid in terms of $r$.
Unfortunately, the full expressions for the exact solution represented by
$p_{r1},p_{\perp 1},p_{r2},$ and $p_{\perp 2}$ are too long to be reported
here. We refer the interested reader to (153) of Appendix C.

The energy densities of both fluids are

$\displaystyle\mu_{1}=$ $\displaystyle\frac{1}{3}\gamma
c_{1}\left(12a^{2}+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right)^{\mathfrak{z}_{3}}\exp\beta,$
(100a) $\displaystyle\mu_{2}=$
$\displaystyle\frac{1}{12}\left(\frac{(4b+1)\left(9a^{2}-7abr^{2}+2b^{2}r^{4}\right)}{\left(br^{2}-a\right)^{3}}\right.$
$\displaystyle\left.-4\gamma
c_{1}\left(12a^{2}+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right)^{\mathfrak{z}_{3}}\exp\beta\right),$
(100b)

where

$\mathfrak{z}_{3}=\frac{-4b+6\gamma+5}{16b-8},$ (101)

and we have used (154).

The energy densities are shown in Fig. 16. The pressures for both fluids are
plotted in Figs. 14 and 15. The radial pressures of both fluids do not
approach zero at the same value of $r$ as the radial pressure of fluid 2
approaches zero before that of fluid 1. Thus, this solution represents a
shell.

Figure 14: Radial pressure for fluid 1 (blue) and fluid 2 (green), and
tangential pressure for fluid 2 (red) for the Florides-class solution in (153)
vs $r$ with $a=0.45$, $b=-1.1$, $c_{1}=6$ and $\gamma=0.3$. Figure 15:
$dp_{r}/d\mu$ (black) and $dp_{\perp}/d\mu$ (yellow) for fluid 2 for the
Florides-class solution in (153) vs $r$ with $a=0.45$, $b=-1.1$, $c_{1}=6$ and
$\gamma=0.3$. Figure 16: Energy density profile for the Florides-class
solution in (100) for fluid 1 (gray) and fluid 2 (cyan) vs $r$ with $a=0.45$,
$b=-1.1$, $c_{1}=6$ and $\gamma=0.3$.

#### The case $\mathbb{P}_{1}\neq 0$

In the general case with both fluids being anisotropic, we follow the same
procedure as in the $\mathbb{P}_{1}=0$ case and use the equation of state
(98). The relation (85) then gives $\mathbb{P}_{1}$. Using (79), we can
calculate $P_{1}$ from $\mathcal{P}_{1}$. The corresponding quantities for
fluid 2 are calculated by subtracting from the total pressures. As usual, we
utilize $\phi$ from (55) with (9) to give $p_{1},p_{2},\Pi_{1},$ and $\Pi_{2}$
in terms of $r$. Using (5), we can write the radial and tangential pressures
for each fluid in terms of $r$. The solution is too lengthy to include here,
so we have listed it as (155) in Appendix C.

The energy densities of both fluids are

$\displaystyle\mu_{1}=$ $\displaystyle\frac{1}{3}\gamma
c_{1}\left(12a^{2}+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right)^{\mathfrak{z}_{3}}\exp\beta,$
(102a) $\displaystyle\mu_{2}=$
$\displaystyle\frac{1}{12}\left(\frac{(4b+1)\left(9a^{2}-7abr^{2}+2b^{2}r^{4}\right)}{\left(br^{2}-a\right)^{3}}\right.$
$\displaystyle\left.-4\gamma
c_{1}\left(12a^{2}+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right)^{\mathfrak{z}_{3}}\exp\beta\right),$
(102b)

where we have used (101) and (154). The energy densities are shown in Fig. 16.
The pressures for both fluids are plotted in Figs. 17, 18 and 19. The radial
pressure of fluid 2 approaches zero before fluid 1. This solution represents a
shell.

Figure 17: Radial pressure (blue) and tangential pressure (orange) for the
Florides-class solution in (155) for fluid 1 vs $r$ with $a=0.6$, $b=-0.6$,
$\mathfrak{c}_{1}=7.8$ and $\gamma=0.3$. This graph is truncated to the value
of $r$ that the radial pressure of fluid 2 approaches zero. Figure 18: Radial
pressure (green) and tangential pressure (red) for the Florides-class solution
in (155) for fluid 2 vs $r$ with $a=0.6$, $b=-0.6$, $\mathfrak{c}_{1}=7.8$ and
$\gamma=0.3$. Figure 19: $dp_{r}/d\mu$ (black) and $dp_{\perp}/d\mu$ (yellow)
for fluid 2 for the Florides-class solution in (155) for fluid 2 vs $r$ with
$a=0.6$, $b=-0.6$, $\mathfrak{c}_{1}=7.8$ and $\gamma=0.3$. Figure 20: Energy
density profile for the Florides-class solution in (102) for fluid 1 (gray)
and fluid 2 (cyan) vs $r$ with $a=0.6$, $b=-0.6$, $\mathfrak{c}_{1}=7.8$ and
$\gamma=0.3$.

##  VII Solution Reconstruction

In [23], we formulated a reconstruction algorithm for the two-fluid isotropic
TOV equations, by expanding the technique that was proposed in [25, 26] to
obtain new two-fluid solutions starting from a given metric. In this section,
we generalize this technique to include anisotropy.

We start by rearranging (10b) as

$X=-\frac{1}{2}-\frac{\mathcal{K}_{,\rho}}{2\mathcal{K}}.$ (103)

Equation (11) gives the relation

$\mathbb{M}_{tot}=-\frac{1}{4}+\mathcal{K}+\frac{\mathcal{K}_{,\rho}}{2\mathcal{K}},$
(104)

Utilizing (10a) together with (104), we obtain

$\displaystyle{P}_{tot}$ $\displaystyle=P_{1}+P_{2}$ (105)
$\displaystyle=\frac{1}{3}\left[2Y_{,\rho}+2Y^{2}+Y\right]-\frac{1}{3}\mathcal{K}-\left(\frac{2Y+1}{6}\right)\frac{\mathcal{K}_{,\rho}}{\mathcal{K}}+\frac{1}{12}.$

The total anisotropic pressure is obtained from (12) with (105), and is given
by

$\mathbb{P}_{tot}=\frac{\mathcal{K}_{,\rho}-4\mathcal{K}^{2}+\mathcal{K}-4\mathcal{K}Y+2Y\mathcal{K}_{,\rho}-4\mathcal{K}Y^{2}+4\mathcal{K}Y}{6\mathcal{K}}.$
(106)

In the isotropic case the condition $\mathbb{P}_{tot}=0$ sets the first order
relationship

$0=\frac{\mathcal{K}_{,\rho}-4\mathcal{K}^{2}+\mathcal{K}-4\mathcal{K}Y_{,\rho}+2Y\mathcal{K}_{,\rho}-4\mathcal{K}Y^{2}+4\mathcal{K}Y}{6\mathcal{K}}.$
(107)

between $\mathcal{K}$ and $Y$. Hence, whenever this relationship is violated,
the metric will be only compatible with anisotropic fluids. In seeking
physically viable solutions, it is convenient to violate the constraint (107)
using $Y$ and $\mathcal{K}$ from known physically viable isotropic solutions.
This will increase the probability of obtaining a well-behaved reconstructed
solution.

The analysis in [23] in the case of two isotropic fluids showed a degeneracy
in the equations above as a given spacetime metric might correspond to many
different combinations of fluids. However, as the conservation equations for
the individual fluids are independent equations, one can use these relations
to break this degeneracy. The same happens in the anisotropic case. Thus we
need to add either (10c) or (10d) to our reconstruction equations. Without
loss of generality, we can choose (10c) i.e. the equation for fluid 1. We then
obtain, setting the interaction to zero,

$\begin{split}{P}_{1,\rho}+\mathbb{P}_{1,\rho}=&-Y(\mathbb{M}_{1}+P_{1})-\mathbb{P}_{1}\left(2X+Y+\frac{3}{2}\right)\\\
&-2XP_{1},\end{split}$ (108)

a first order differential equation.

Depending on the properties of the source fluids, there are different possible
resolution approaches. For example, in cases where one fluid is isotropic, we
can obtain $\mathbb{P}_{tot}=\mathbb{P}_{2}$ from (106), if the metric and,
therefore, $Y$ and $\mathcal{K}$ are assigned. At this point, we can derive
$P_{tot}$ using (105). Then $P_{1}$ can be determined with an ansatz such as
the equation of state, while $P_{2}$ is derived by subtracting $P_{1}$ from
$P_{tot}$. In addition, $\mathbb{M}_{tot}$ is derived from the constraint (12)
with (104).

For the cases in which both fluids are anisotropic, additional information is
required. This can be obtained, for example, by assuming an equation of state
relating the anisotropic pressure $\mathbb{P}$ to the density $\mathbb{M}$ and
pressure $P$, or by giving a relation between $\mathbb{M}$ and $\mathcal{K}$.

Despite the sheer number of potential combinations of solutions one could use
in the reconstruction, there is no guarantee that the solution obtained will
be physically viable. Here, we will limit ourselves to illustrate one example.
In particular, we will perform the reconstruction using the interior
Schwarzschild and Tolman IV solutions.

We choose $Y$ from the Interior Schwarzschild solution and $\mathcal{K}$ from
the Tolman IV solution:

$\displaystyle Y$
$\displaystyle=\frac{e^{\rho}\mu_{1}}{-2c_{1}\sqrt{3-e^{\rho}\mu_{1}}+2e^{\rho}\mu_{1}-6},$
(109a) $\displaystyle\mathcal{K}$
$\displaystyle=-\frac{R^{2}\left(A^{2}+2e^{\rho}\right)}{4\left(A^{2}+e^{\rho}\right)\left(e^{\rho}-R^{2}\right)},$
(109b)

and we assume the equation of state

$\mathbb{M}_{1}=aP_{1},$ (110)

where $a$ is a constant. We assume fluid 1 to be isotropic, so
$\mathbb{P}_{1}=0$. Eq. (108) reduces to

${P}_{1,\rho}=-Y(\mathbb{M}_{1}+P_{1})-2XP_{1}.$ (111)

As mentioned above with $Y$ and $\mathcal{K}$ we can determine all the total
quantities. Then, using (110) we can solve (111) for $P_{1}$. In addition,
$P_{2}$ can be obtained by subtracting $P_{1}$ from $P_{tot}$ and, since
$\mathbb{P}_{1}=0$, $\mathbb{P}_{2}=\mathbb{P}_{tot}$. Finally,
$\mathbb{M}_{2}$ is obtained by subtracting $\mathbb{M}_{1}$ from
$\mathbb{M}_{tot}$. For the reconstructed quantities, we have (156) in
Appendix D.

To obtain $p_{1},p_{2},\Pi_{1},$ and $\Pi_{2}$ in terms of $r$, we use $\phi$
for the Tolman IV solution (63), with (9). Finally, using (5), we can write
the radial and tangential pressures for each fluid in terms of $r$. These
quantities are too long to be reported here, and are given in (158) of
Appendix D. The energy densities of both fluids are given in (160) of Appendix
D and are shown in Fig. 23. The pressures for both fluids are plotted in Fig.
21, while Fig. 22 shows the square of the sound speed for fluid 2 (we have
omitted fluid 1 from the plot since this quantity has a constant value of
$\frac{1}{a}$ for fluid 1). Since the radial pressures of both fluids vanish
at different values of $r$ we can conclude that this solution represents a
shell.

Figure 21: Radial pressure for the solution in (158) for fluid 1 (blue) and
fluid 2 (green), and tangential pressure for fluid 2 (red) vs $r$ with $a=3$,
$\mu_{1}=0.9$, $R=1.75$, $A=1.95$, $\mathfrak{c}_{1}=-0.2$, and $c_{1}=-3$.
Figure 22: $dp_{r}/d\mu$ (black) and $dp_{\perp}/d\mu$ (yellow) for the
solution in (158) for fluid 2 vs $r$ with $a=3$, $\mu_{1}=0.9$, $R=1.75$,
$A=1.95$, $\mathfrak{c}_{1}=-0.2$, and $c_{1}=-3$. Figure 23: Energy density
profile for the solution in (160) for fluid 1 (gray) and fluid 2 (cyan) vs $r$
with $a=3$, $\mu_{1}=0.9$, $R=1.75$, $A=1.95$, $\mathfrak{c}_{1}=-0.2$, and
$c_{1}=-3$.

##  VIII Generating theorems

Boonserm et al. developed generating/transformation theorems which “map one
perfect fluid sphere into another” [14, 15]. The term perfect fluid sphere
represents a static, isotropic, spherically symmetric fluid distribution in
curved spacetimes which, in turn, can be associated to the relativistic
stellar objects described by the TOV equations. Using generating theorems, one
can start from a known solution of the TOV equations and obtain new solutions.
Carloni et al. demonstrated that in the context of covariant TOV equations,
the generating theorems are equivalent to a linear deformation of a given
solution [25, 26]. In this section, we will show that these theorems also
exist in the case of anisotropic multifluid solutions.

Indeed one can find several different generating theorems, depending on the
quantities that can be “deformed”. For brevity, we will show explicitly some
specific cases, leaving a complete taxonomy for a future work. The selected
cases will illustrate: (i) deforming a two-fluid anisotropic solution to a new
two-fluid anisotropic solution, (ii) deforming an isotropic solution to obtain
an anisotropic one, and finally, (iii) starting with two isotropic fluids and
obtaining a solution with more than two anisotropic fluids.

For all of these examples, we will be deforming the quantity $Y$ as

$Y\rightarrow\bar{Y}+\tilde{Y},$ (112)

where, here and in the following, the bar indicates the known solution and
tilde represents the deformed quantities. We can then use the definition of
$Y$ in terms of the metric coefficient $k_{1}$ in (32), to obtain the
transformation of the metric that corresponds to this deformation:

$\displaystyle k_{1}$ $\displaystyle\rightarrow
a_{1}\bar{k}_{1}\exp\left(2\int{\tilde{Y}}d\rho\right),$ (113a) $\displaystyle
k_{2}$ $\displaystyle\rightarrow\bar{k}_{2},$ (113b) $\displaystyle k_{3}$
$\displaystyle\rightarrow\bar{k}_{3},$ (113c)

where $a_{1}$ is a constant of integration.

In what follows, we will connect the deformation of $Y$ given by (112) to
deformations in other geometrical and thermodynamical quantities. We will even
consider a case in which the deformation allows the mapping to a solution with
an additional fluid. The deformed matter variables, together with the new
metric coefficients, represent the new solution.

###  VIII.1 Case 1: Two anisotropic fluids $\rightarrow$ Two anisotropic
fluids

As a first example we will start from a solution of the TOV equations
consisting of two anisotropic fluids.

The $Y$ deformation (112) gives a new solution of Eqs. (18) if the variables
in the TOV equations are deformed in the following way

$\displaystyle{P}_{1}$ $\displaystyle\rightarrow\bar{{P}}_{1},$ (114a)
$\displaystyle{P}_{2}$ $\displaystyle\rightarrow\bar{{P}}_{2},$ (114b)
$\displaystyle\mathbb{P}_{1}$
$\displaystyle\rightarrow\bar{\mathbb{P}}_{1}+\tilde{\mathbb{P}}_{1},$ (114c)
$\displaystyle\mathbb{P}_{2}$ $\displaystyle\rightarrow\bar{\mathbb{P}}_{2},$
(114d) $\displaystyle\mathbb{M}_{1}$
$\displaystyle\rightarrow\bar{\mathbb{M}}_{1}+\tilde{\mathbb{M}}_{1},$ (114e)
$\displaystyle\mathbb{M}_{2}$
$\displaystyle\rightarrow\bar{\mathbb{M}}_{2}-\tilde{\mathbb{M}}_{1},$ (114f)
$\displaystyle{\mathcal{K}}$ $\displaystyle\rightarrow\bar{\mathcal{K}}.$
(114g)

Notice that $\mathbb{M}_{tot}$ is not deformed. Let us now prove that the
deformations can be obtained in an analytical, albeit formal way. Substituting
the deformations in the constraint (12) gives

$\tilde{Y}=\tilde{\mathbb{P}}_{1},$ (115)

i.e. the deformation in $Y$ has the same analytical form as the deformation in
the anisotropic pressure of fluid 1. Using this result, and taking into
account that Eqs. (18) hold for the barred quantities, the deformed Eq. (22)
takes the form

$\displaystyle\tilde{\mathbb{P}}_{1,\rho}$ $\displaystyle=$
$\displaystyle-\frac{1}{4}\tilde{\mathbb{P}}_{1}\left(12\bar{\mathcal{K}}-4\bar{\mathbb{M}}_{1}-4\bar{\mathbb{M}}_{2}+8\bar{P}_{1}+8\bar{\mathbb{P}}_{1}\right.$
(116)
$\displaystyle\left.+8\bar{P}_{2}+8\bar{\mathbb{P}}_{2}-1\right)-\tilde{\mathbb{P}}_{1}^{2},$

which is a Bernoulli equation for $\tilde{\mathbb{P}}_{1}$ and has a formal
general solution.

We now need to determine the energy densities of the fluids in the new
solution. Substituting the deformations in the TOV equation for fluid one (18)
and remembering that the same equation holds for the barred variables, we can
solve for $\tilde{\mathbb{M}}_{1}$ to obtain

$\displaystyle\tilde{\mathbb{M}}_{1}$ $\displaystyle=$
$\displaystyle\frac{\tilde{\mathbb{P}}_{1}}{4\left(\tilde{\mathbb{P}}_{1}+\bar{Y}\right)}\left(12\bar{\mathcal{K}}+4\bar{\mathbb{M}}_{2}+12\bar{P}_{1}+12\bar{\mathbb{P}}_{1}+16\bar{P}_{2}\right.$
(117) $\displaystyle\left.+16\bar{\mathbb{P}}_{2}-12\bar{Y}-3\right).$

Hence with the above procedure, we have started with two anisotropic fluids
and, after the deformation, we have a modified geometry and two new
anisotropic fluids. Both fluids have energy densities that are different to
those of the known solution and one of them has a modified anisotropic
pressure.

Notice that in performing the above procedure, we need to ensure that all
thermodynamical variables are positive, so we must impose that
$\mathbb{M}_{2}$ is positive as an additional condition on the solution of the
above equations. This is achieved by verifying that:

$\bar{\mathbb{M}}_{2}-\tilde{\mathbb{M}}_{1}>0.$ (118)

###  VIII.2 Case 2: Two isotropic fluids $\rightarrow$ Two anisotropic fluids

As a second example, we consider a solution comprising of two isotropic fluids
($\bar{\mathbb{P}}_{1}=0$ and $\bar{\mathbb{P}}_{2}=0$). We now deform $Y$ as
in (112) and we require that after the deformation the new solution will have
two anisotropic fluids with the same energy density. This result can be
achieved by assuming that the fluid and geometry variables are deformed as

$\displaystyle{P}_{1}$ $\displaystyle\rightarrow\bar{{P}}_{1},$ (119a)
$\displaystyle{P}_{2}$ $\displaystyle\rightarrow\bar{{P}}_{2},$ (119b)
$\displaystyle\mathbb{P}_{1}$
$\displaystyle\rightarrow\tilde{\mathbb{P}}_{1},$ (119c)
$\displaystyle\mathbb{P}_{2}$
$\displaystyle\rightarrow\tilde{\mathbb{P}}_{2},$ (119d)
$\displaystyle\mathbb{M}_{1}$ $\displaystyle\rightarrow\bar{\mathbb{M}}_{1},$
(119e) $\displaystyle\mathbb{M}_{2}$
$\displaystyle\rightarrow\bar{\mathbb{M}}_{2},$ (119f)
$\displaystyle{\mathcal{K}}$ $\displaystyle\rightarrow\bar{\mathcal{K}}.$
(119g)

Applying these deformations and taking into account that (12) holds for the
barred quantities associated to the starting solution, we obtain

$\tilde{Y}=\tilde{\mathbb{P}}_{1}+\tilde{\mathbb{P}}_{2}=\tilde{\mathbb{P}}_{tot}.$
(120)

At this point, since the barred quantities satisfy by themselves the TOV
equations, the equation for the total pressure (22) can be written as a
Bernoulli equation for the total deformed anisotropic pressure
$\tilde{\mathbb{P}}_{tot}$

$\displaystyle\tilde{\mathbb{P}}_{tot}(12\bar{\mathcal{K}}-4\bar{\mathbb{M}}_{1}-4\bar{\mathbb{M}}_{2}+8\bar{P}_{1}+8\bar{P}_{2}-1)$
$\displaystyle+4\tilde{\mathbb{P}}_{tot,\rho}+4\tilde{\mathbb{P}}_{tot}^{2}=0.$
(121)

The TOV equation for fluid 1, (18), under deformation gives the following
equation for $\tilde{\mathbb{P}}_{1}$

$\displaystyle\tilde{\mathbb{P}}_{1,\rho}$ $\displaystyle=$
$\displaystyle\tilde{\mathbb{P}}_{1}\left(2\bar{\mathbb{M}}_{1}+2\bar{\mathbb{M}}_{2}+2\bar{P}_{1}+2\bar{P}_{2}-\tilde{\mathbb{P}}_{tot}\right.$
(122)
$\displaystyle\left.-3\bar{Y}-\frac{1}{2}\right)-\tilde{\mathbb{P}}_{tot}(\bar{\mathbb{M}}_{1}+3\bar{P}_{1}).$

Then we can derive $\tilde{\mathbb{P}}_{2}$ using Eq. (120). As the Bernoulli
equations have formal general solutions we can conclude that we have obtained
a new solution with two anisotropic fluids.

###  VIII.3 Case 3: Two isotropic fluids $\rightarrow$ Three anisotropic
fluids

Our final example involves a change in both the nature and number of fluids.
We begin with two isotropic fluids and we perform the deformation

$\displaystyle{P}_{1}$ $\displaystyle\rightarrow\bar{P}_{1}+{P}_{3},$ (123a)
$\displaystyle{P}_{2}$ $\displaystyle\rightarrow\bar{P}_{2},$ (123b)
$\displaystyle\mathbb{P}_{1}$ $\displaystyle\rightarrow\mathbb{P}_{3},$ (123c)
$\displaystyle\mathbb{P}_{2}$ $\displaystyle\rightarrow-\mathbb{P}_{3},$
(123d) $\displaystyle\mathbb{M}_{1}$
$\displaystyle\rightarrow\bar{\mathbb{M}}_{1},$ (123e)
$\displaystyle\mathbb{M}_{2}$ $\displaystyle\rightarrow\bar{\mathbb{M}}_{2},$
(123f) $\displaystyle{\mathcal{K}}$
$\displaystyle\rightarrow\bar{\mathcal{K}}.$ (123g)

where now, we consider $P_{3}$ not as a change in $P_{1}$, but as the pressure
of an additional fluid present only in the new solution. In this picture, the
pressure $P_{3}$ must satisfy an equation similar to Eq. (18) for $P_{1}$ and
Eq. (19) for $P_{2}$:

$P_{3,\rho}+\mathbb{P}_{3,\rho}=P_{3}(-2X-Y)-\mathbb{P}_{3}\left(2X+Y+\frac{3}{2}\right)-\mathbb{M}_{3}Y.$
(124)

We will assume: (i) that the third fluid is anisotropic, so in the above
equation $\mathbb{P}_{3}\neq 0$, and (ii) after the deformation, fluids 1 and
2 are anisotropic as in (123). Notice that as we have seen in Sec. III, the
anisotropic term $\Pi$ can be negative, so we can consider $\mathbb{P}_{1}$ or
$\mathbb{P}_{2}$ less than zero.

Applying these deformations and solving the constraint (12) with the initial
solution, we obtain the relation

$\tilde{Y}=P_{3}.$ (125)

Applying the deformation and the result (125) to Eq. (18) with no interaction
and remembering that the barred quantities satisfy (18), we obtain the
following equation for $P_{3}$

$\begin{split}P_{3,\rho}=-&P_{3}^{2}-\bigg{(}3\bar{\mathcal{K}}-\bar{\mathbb{M}}_{1}-\bar{\mathbb{M}}_{2}+2\bar{P}_{1}+2\bar{P}_{2}-\frac{7}{4}\bigg{)},\end{split}$
(126)

which is a Bernoulli equation. Now, substituting the above results and Eq.
(15) into (124), we obtain the equation

$\displaystyle\mathbb{P}_{3,\rho}$ $\displaystyle=$
$\displaystyle\mathbb{P}_{3}\left[\frac{1}{4}\left(8\bar{\mathbb{M}}_{1}+8\bar{\mathbb{M}}_{2}+8\bar{P}_{1}+8\bar{P}_{2}-12\bar{Y}-2\right)\right.$
(127)
$\displaystyle\left.-P_{3}\right]+\frac{1}{4}P_{3}\left(12\bar{\mathcal{K}}+4\bar{\mathbb{M}}_{2}+12\bar{P}_{1}+16\bar{P}_{2}\right.$
$\displaystyle\left.-12\bar{Y}-3\right).$

Finally, substituting (126) and (127) into (22) we can solve for
$\mathbb{M}_{3}$:

$\mathbb{M}_{3}=\frac{P_{3}(\bar{\mathbb{M}}_{1}+\bar{P}_{1})}{P_{3}+\bar{Y}}.$
(128)

Hence, we have obtained three anisotropic fluids from the deformation of two
isotropic fluids. The pressure $P$ of the first fluid was deformed such that
it contained the pressure of the known solution and an additional term which
represents the pressure of the third fluid. The pressure of the second fluid
was not deformed, and the anisotropic pressures of both fluid 1 and 2 were
deformed such that they contained the anisotropic pressure of fluid 3. The
energy densities of fluids 1 and 2 were not deformed, and the energy density
of fluid 3 was calculated.

##  IX Discussion and conclusion

In this paper, we have presented the 1+1+2 semi-tetrad equations for the two-
fluid anisotropic case, along with its generalisation to $N$-fluids, and we
have used these equations to construct a covariant multifluid version of the
TOV equations. Such generalization is needed because anisotropy is a physical
property encountered in the study of most compact objects.

The general structure of these equations revealed several interesting facts.
For example, the weak energy condition showed explicitly that the anisotropic
pressure term $\Pi$ can be negative, provided that (43) holds. Allowing for
negative $\Pi$ implies that the the non radial pressure in an object is bigger
than the radial one. In addition, at the boundary, the junction conditions
require that the radial pressure vanishes, but there is no restriction on the
behaviour of the tangential pressure. Therefore, the tangential pressure does
not necessarily need to vanish at the boundary.

In order to solve the covariant TOV equations analytically, we have proposed
several resolution strategies in the case of two fluids. In particular, using
some well-known single fluid solutions we explored some direct resolution
algorithms as well as some indirect ones, namely, reconstruction methods and
generating theorems.

As a first example we derived the two-fluid generalization of the simplest
known single fluid anisotropic solution: the Bowers-Liang solution. Our
generalization has the same geometry ad the original solution, but is composed
of two different fluids. Its structure allows one to see explicitly the
changes induced on the thermodynamical quantities by the presence of
anisotropy. Specifically, the radial pressure becomes lower at the core and
higher on the surface, when compared to the isotropic case. This change in the
radial pressure implies that the matter composition at the core is different,
i.e., there is a different equation of state for this interior region. This
gives the chance to observe explicitly the structural changes [10] that
anisotropy induces in a given object.

We then applied the same strategy to other known solutions. From an element of
the Florides-class of solutions, we were able to obtain six new anisotropic
two-fluid solutions. In some of these solutions, one of the fluids was
barotropic, a feature which does not seem to be possible in the single fluid
case or when all fluids are isotropic. None of these solutions, however, could
represent an entire object as they could not be smoothly soldered to a vacuum
exterior. As we have commented, solutions such as these represent only a
specific shell in the compact object and requires matching to an additional
shell before being connected to an exterior solution, such as the
Schwarzschild [6] or Vaidya [7] solution. Such “shell structure” is expected
to arise in multifluid systems, in which the properties of matter induces
layers with different mechanical and thermodynamical properties.

A less straightforward approach, which yields interesting results, is given by
the so-called reconstruction techniques. We found that it is possible to
extend the reconstruction method proposed in [25, 26] to the case of
anisotropic multifluid objects. In particular, we found that by mixing the
metric components of different well-known isotropic solutions one can
reconstruct new multifluid anisotropic solutions with desirable features. We
presented one solution of this type, constructed from the interior
Schwarzschild and Tolman IV solutions.

Finally, we considered the extension of the generating theorems to the case of
multifluid anisotropic solutions. It turns out that one can find a great
number of such theorems that are able to connect solutions with different
types of fluids or even with a different number of fluids. We have limited
ourselves here to three specific examples. A first, in which a solution
comprising two anisotropic fluids is mapped to a new solution in which one of
the fluids has a different anisotropic pressure. A second, in which a two-
fluid anisotropic solution is obtained from a two-fluid isotropic one. Finally
a third example, in which these theorems are used to change the number and
character of fluids in a given solution. As it is well known, the generating
theorems can be used to form chains of solutions of the TOV equations. We find
here that this is even more true in the case of multifluid anisotropic
solutions. We will extend this investigation and examine some of the other
variations to these theorems in a future work.

###### Acknowledgements.

NFN would like to acknowledge funding from the CoE-MaSS and Oppenheimer
Memorial Trust.

##  APPENDIX A Additional details of the 1+1+2 covariant approach for a single
fluid

In this section we present the main aspects of the 1+1+2 covariant approach
that are not covered in Sec. II. The reader is referred to [24, 32, 33] for
more details.

The complete list of kinematical variables is given by

$\displaystyle\mathcal{A}=e_{a}\dot{u}^{a}\hskip 5.69046pt,\hskip
5.69046pt\mathcal{A}_{a}=N_{ab}\dot{u}^{b}\hskip 5.69046pt,\hskip
5.69046pt\Theta=D_{a}u^{a},$ (129a)
$\displaystyle\Omega=\frac{1}{2}\varepsilon^{abc}D_{[a}u_{b]}e_{c}\hskip
5.69046pt,\hskip
5.69046pt{\Omega}^{a}=\frac{1}{2}\epsilon^{cbd}D_{[c}u_{b]}{N_{d}}^{a},$
(129b) $\displaystyle\Sigma={\sigma}^{ab}e_{a}e_{b}\hskip 5.69046pt,\hskip
5.69046pt{\Sigma}_{a}={\sigma}_{cd}e^{d}{N_{a}}^{c},$ (129c)
$\displaystyle{\Sigma}_{ab}=\left({N^{c}}_{(a}{N_{b)}}^{d}-\frac{1}{2}N_{ab}N^{cd}\right){\sigma}_{cd},$
(129d)
$\displaystyle{\sigma}_{ab}=\left({h^{c}}_{(a}{h_{b)}}^{d}-\frac{1}{3}h_{ab}h^{cd}\right)D_{c}u_{d},$
(129e) $\displaystyle a_{b}=e^{c}D_{c}e_{b}=\hat{e}_{b}\hskip 5.69046pt,\hskip
5.69046pt\phi={\delta}_{a}e^{a},$ (129f)
$\displaystyle{\zeta}_{ab}=\left({N^{c}}_{(a}{N_{b)}}^{d}-\frac{1}{2}N_{ab}N^{cd}\right){\delta}_{c}e_{d},$
(129g) $\displaystyle\xi=\frac{1}{2}\varepsilon^{ab}{\delta}_{a}e_{b},$ (129h)
$\displaystyle\mathcal{E}={C}_{acbd}u^{c}u^{d}e^{a}e^{b},$ (129i)
$\displaystyle\mathcal{E}_{a}=C_{cdef}u^{e}u^{f}e^{c}{N^{d}}_{a}\hskip
5.69046pt,\hskip 5.69046pt\mathcal{E}_{ab}=C_{\\{acb\\}d}u^{c}u^{d},$ (129j)
$\displaystyle\mathcal{H}=\frac{1}{2}{\varepsilon}_{ade}{C^{de}}_{bc}u^{c}e^{a}e^{b},$
(129k)
$\displaystyle\mathcal{H}_{a}=\frac{1}{2}\varepsilon_{cbe}{C^{be}}_{df}u^{f}e^{c}{N^{d}}_{a},$
(129l)
$\displaystyle\mathcal{H}_{ab}=\frac{1}{2}{\varepsilon_{\\{ade}}{C^{de}}_{b\\}f}u^{f},$
(129m)

where $\varepsilon_{ab}\equiv\varepsilon_{abc}e^{c}$ and
$\varepsilon_{abc}=\eta_{dabc}u^{d}$ are the volumes of the two hypersurfaces,
and $C_{cdef}$ is the Weyl tensor. The symmetrisation over the indices of a
tensor is represented as $T_{(ab)}=\frac{1}{2}(T_{ab}+T_{ba})$, and the anti-
symmetrisation as $T_{[ab]}=\frac{1}{2}(T_{ab}-T_{ba})$. We use curly brackets
$\\{\\}$ to denote the Projected Symmetric Trace-Free part of a tensor with
respect to $e^{a}$:

$X_{\\{ab\\}}\equiv\left(N^{c}{}_{(a}N_{b)}{}^{d}-\frac{1}{2}N_{ab}N^{cd}\right)\mathds{X}_{cd}.$
(130)

According to [33], the remaining 1+1+2 scalars which fully describe the
spacetime can be divided into 3 categories:

Propagation: $\displaystyle\hat{\phi}$
$\displaystyle=-\frac{1}{2}{\phi}^{2}-\frac{2}{3}\mu-\mathcal{E}-\frac{1}{2}\Pi,$
(131b) $\displaystyle Q$ $\displaystyle=0,$ (131c)
$\displaystyle\hat{\mathcal{E}}-\frac{1}{3}\hat{\mu}+\frac{1}{2}\hat{\Pi}$
$\displaystyle=-\frac{3}{2}{\phi}\left(\mathcal{E}+\frac{1}{2}\Pi\right).$
(131d) Evolution: $\displaystyle 0$
$\displaystyle=-\mathcal{A}\phi+\frac{1}{3}(\mu+3p)-\mathcal{E}+\frac{1}{2}\Pi,$
(131e) $\displaystyle 0$ $\displaystyle=\frac{1}{2}\phi Q.$ (131f)
Propagation/evolution: $\displaystyle\hat{\mathcal{A}}$
$\displaystyle=-\mathcal{A}\left(\mathcal{A}+\phi\right)+\frac{1}{2}(\mu+3p),$
(131g) $\displaystyle\hat{Q}$
$\displaystyle=-Q\left(\phi+2\mathcal{A}\right)+j_{u},$ (131h)
$\displaystyle\hat{p}+\hat{\Pi}$
$\displaystyle=-\Pi\left(\frac{3}{2}\phi+\mathcal{A}\right)-\mathcal{A}\left(\mu+p\right)+j_{e},$
(131i) $\displaystyle K$
$\displaystyle=\frac{1}{3}\mu-\mathcal{E}-\frac{1}{2}\Pi+\frac{1}{4}{\phi}^{2},$
(131j)

where $\mu$, $p$, $Q$, and $\pi$ represent, in general, the total energy
density, pressure, total heat flux, and total anisotropic pressure of the
fluid respectively. In addition, we have included the total particle
interaction currents as $j_{u}$ and $j_{e}$ according to the definition

$j_{a}=j_{u}u_{a}+j_{e}e_{a},$ (132)

with $j_{u}$ and $j_{e}$ representing the $u_{a}$ and $e_{a}$ components,
respectively.

Finally, combining equations (131), one can show that the Gauss curvature $K$
satisfies the propagation equation

$\hat{K}=-\phi K.$ (133)

##  APPENDIX B Multifluid 1+1+2 equations

In the case of a static spherically symmetric spacetime with $N$ different
interacting fluids, the 1+1+2 equations are

$\displaystyle\hat{\phi}$
$\displaystyle=-\frac{1}{2}{\phi}^{2}-\frac{2}{3}\sum_{i=1}^{N}{\mu}_{i}-\mathcal{E}-\frac{1}{2}\sum_{i=1}^{N}{\Pi}_{i},$
(134a) $\displaystyle\sum_{i=1}^{N}Q_{i}$ $\displaystyle=0,$ (134b)
$\displaystyle\hat{\mathcal{E}}+\frac{3}{2}{\phi}\mathcal{E}$
$\displaystyle=\frac{1}{3}\sum_{i=1}^{N}{\hat{\mu}}_{i}-\frac{1}{2}\sum_{i=1}^{N}{\hat{\Pi}}_{i}-\frac{3}{2}{\phi}\frac{1}{2}\sum_{i=1}^{N}{\Pi}_{i},$
(134c) $\displaystyle\mathcal{E}+\mathcal{A}\phi$
$\displaystyle=\frac{1}{3}\sum_{i=1}^{N}({\mu}_{i}+3p_{i})+\frac{1}{2}\sum_{i=1}^{N}{\Pi}_{i},$
(134d) $\displaystyle K$
$\displaystyle=\frac{1}{4}{\phi}^{2}-\mathcal{E}+\frac{1}{3}\sum_{i=1}^{N}{\mu}_{i}-\frac{1}{2}\sum_{i=1}^{N}{\Pi}_{i},$
(134e) $\displaystyle\hat{\mathcal{A}}$
$\displaystyle=-\mathcal{A}\left(\mathcal{A}+\phi\right)+\frac{1}{2}\sum_{i=1}^{N}({\mu}_{i}+3p_{i}),$
(134f) $\displaystyle\hat{Q}_{i}$
$\displaystyle=-Q_{i}\left(\phi+2\mathcal{A}\right)+\sum_{k\neq
i}^{N}j_{u}^{(i,k)},$ (134g) $\displaystyle\hat{p}_{i}+\hat{\Pi}_{i}$
$\displaystyle=\sum_{k\neq
i}^{N}j_{e}^{(i,k)}-{\Pi}_{i}\left(\frac{3}{2}\phi+\mathcal{A}\right)-\mathcal{A}\left({\mu}_{i}+p_{i}\right),$
(134h)

where the index $i$ represents the $i$th component and $j_{e}^{(i,k)}$ is the
interaction term between the component $i$ and $k$. Notice that
$j_{u}^{(i,j)}=-j_{u}^{(j,i)}$ and $j_{e}^{(i,j)}=-j_{e}^{(j,i)}$.

Summing (134h) over $i$ we have the equation for the total energy pressure:

$\displaystyle\hat{p}_{tot}+\hat{\Pi}_{tot}$
$\displaystyle=-\sum_{i=1}^{N}\Pi_{1}\left(\frac{3}{2}\phi+\mathcal{A}\right)-\mathcal{A}\sum_{i=1}^{N}\left({\mu}_{i}+p_{i}\right).$
(135)

In terms of the variable $\rho$, defined by the relation
$K=K_{0}^{-1}e^{-\rho},$ the 1+1+2 equations become

$\displaystyle\phi{\phi}_{,\rho}$
$\displaystyle=-\frac{1}{2}{\phi}^{2}-\frac{2}{3}\sum_{i=1}^{N}{\mu}_{i}-\mathcal{E}-\frac{1}{2}\sum_{i=1}^{N}{\Pi}_{i},$
(136a) $\displaystyle\sum_{i=1}^{N}Q_{i}$ $\displaystyle=0,$ (136b)
$\displaystyle{\mathcal{E}}_{,\rho}+\frac{3}{2}\mathcal{E}$
$\displaystyle=\frac{1}{3}\sum_{i=1}^{N}{{\mu}}_{i,\rho}-\frac{1}{2}\sum_{i=1}^{N}{{\Pi}}_{i,\rho}-\frac{3}{2}{\phi}\frac{1}{2}\sum_{i=1}^{N}{\Pi}_{i},$
(136c) $\displaystyle\mathcal{E}+\mathcal{A}\phi$
$\displaystyle=\frac{1}{3}\sum_{i=1}^{N}({\mu}_{i}+3p_{i})+\frac{1}{2}\sum_{i=1}^{N}{\Pi}_{i},$
(136d) $\displaystyle K$
$\displaystyle=\frac{1}{4}{\phi}^{2}-\mathcal{E}+\frac{1}{3}\sum_{i=1}^{N}{\mu}_{i}-\frac{1}{2}\sum_{i=1}^{N}{\Pi}_{i},$
(136e) $\displaystyle\phi{\mathcal{A}}_{,\rho}$
$\displaystyle=-\mathcal{A}\left(\mathcal{A}+\phi\right)+\frac{1}{2}\sum_{i=1}^{N}({\mu}_{i}+3p_{i}),$
(136f) $\displaystyle\phi\hat{Q}_{i,\rho}$
$\displaystyle=-Q_{i}\left(\phi+2\mathcal{A}\right)+\sum_{k\neq
i}^{N}j_{u}^{(i,k)},$ (136g) $\displaystyle\hat{p}_{i}+\hat{\Pi}_{i}$
$\displaystyle=\sum_{k\neq
i}^{N}j_{e}^{(i,k)}-{\Pi}_{i}\left(\frac{3}{2}\phi+\mathcal{A}\right)-\mathcal{A}\left({\mu}_{i}+p_{i}\right).$
(136h)

Using the variables

$\displaystyle X$ $\displaystyle=\frac{{\phi}_{,\rho}}{\phi},$ $\displaystyle
Y$ $\displaystyle=\frac{\mathcal{A}}{\phi},$ $\displaystyle\mathcal{K}$
$\displaystyle=\frac{K}{{\phi}^{2}},$ (137a) $\displaystyle E$
$\displaystyle=\frac{\mathcal{E}}{{\phi}^{2}},$ $\displaystyle\mathbb{M}_{1}$
$\displaystyle=\frac{{\mu}_{1}}{{\phi}^{2}},$ $\displaystyle\mathbb{M}_{2}$
$\displaystyle=\frac{{\mu}_{2}}{{\phi}^{2}},$ (137b) $\displaystyle P_{1}$
$\displaystyle=\frac{p_{1}}{{\phi}^{2}},$ $\displaystyle P_{2}$
$\displaystyle=\frac{p_{2}}{{\phi}^{2}},$ $\displaystyle\mathbb{P}_{1}$
$\displaystyle=\frac{{\Pi}_{1}}{{\phi}^{2}},$ (137c)
$\displaystyle\mathbb{P}_{2}$ $\displaystyle=\frac{{\Pi}_{2}}{{\phi}^{2}},$
$\displaystyle\mathbb{Q}_{1}$ $\displaystyle=\frac{Q_{1}}{{\phi}^{2}},$
$\displaystyle\mathbb{Q}_{2}$ $\displaystyle=\frac{Q_{2}}{{\phi}^{2}},$ (137d)

we obtain the covariant TOV equations as

$\displaystyle{P}_{i,\rho}+\mathbb{P}_{i,\rho}$ $\displaystyle=$
$\displaystyle\sum_{i\neq
k}^{N}\mathbb{J}_{e}^{(i,k)}-{P_{i}}^{2}-{\mathbb{P}_{i}}^{2}+P_{i}\left[\mathbb{M}_{i}-2\mathbb{P}_{i}-3\mathcal{K}+\frac{7}{4}\right]+\mathbb{P}_{i}\left(\mathbb{M}_{i}-3\mathcal{K}+\frac{1}{4}\right)+\mathbb{M}_{i}\left(\frac{1}{4}-\mathcal{K}\right)$
(138) $\displaystyle-P_{i}\sum_{k\neq
i}^{N}(P_{k}+\mathbb{P}_{k}-2\mathbb{M}_{k})-\mathbb{P}_{i}\sum_{k\neq
i}^{N}(P_{k}+\mathbb{P}_{k}-2\mathbb{M}_{k})-\mathbb{M}_{i}\sum_{k\neq
i}^{N}(P_{k}+\mathbb{P}_{k}),$ $\displaystyle\mathcal{K}_{,\rho}$
$\displaystyle=$ $\displaystyle
2\mathcal{K}\left(\frac{1}{4}-\mathcal{K}+\sum_{k=1}^{N}\mathbb{M}_{k}\right),$
(139) $\displaystyle\mathbb{Q}_{i,\rho}$ $\displaystyle=$
$\displaystyle\mathbb{Q}_{i}\left[2\mathcal{K}-2\sum_{k=1}^{N}\mathbb{M}_{k}-\frac{3}{2}\right]-\sum_{i\neq
k}^{N}\mathbb{J}_{u}^{(i,k)},$ (140)

where, for notational reasons we have set

$\mathbb{J}_{u}^{(i,k)}=\frac{j_{u}^{(i,k)}}{{\phi}^{3}},\quad\mathbb{J}_{e}^{(i,k)}=\frac{j_{e}^{(i,k)}}{{\phi}^{3}}.$
(141)

In the case of two fluids the equations above reduce to

$\displaystyle\hat{\phi}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}{\phi}^{2}-\frac{2}{3}({\mu}_{1}+{\mu}_{2})-\mathcal{E}-\frac{1}{2}({\Pi}_{1}+{\Pi}_{2}),$
(142) $\displaystyle\hat{\mathcal{E}}$ $\displaystyle=$
$\displaystyle\frac{1}{3}(\hat{\mu}_{1}-\hat{\mu}_{2})-\frac{1}{2}(\hat{\Pi}_{1}+\hat{\Pi}_{2})-\frac{3}{2}{\phi}\left(\mathcal{E}+\frac{1}{2}{\Pi}_{1}+\frac{1}{2}{\Pi}_{2}\right),$
(143) $\displaystyle-\mathcal{A}\phi$ $\displaystyle+$
$\displaystyle\frac{1}{3}({\mu}_{1}+3p_{1})+\frac{1}{3}({\mu}_{2}+3p_{2})-\mathcal{E}+\frac{1}{2}{\Pi}_{1}+\frac{1}{2}{\Pi}_{2}=0,$
(144) $\displaystyle\hat{\mathcal{A}}$ $\displaystyle=$
$\displaystyle-\mathcal{A}\left(\mathcal{A}+\phi\right)+\frac{1}{2}({\mu}_{1}+3p_{1})+\frac{1}{2}({\mu}_{2}+3p_{2}),$
(145) $\displaystyle\hat{p}_{1}+\hat{\Pi}_{1}$ $\displaystyle=$
$\displaystyle-{\Pi}_{1}\left(\frac{3}{2}\phi+\mathcal{A}\right)-\mathcal{A}\left({\mu}_{1}+p_{1}\right)+j_{e},$
(146) $\displaystyle\hat{p}_{2}+\hat{\Pi}_{2}$ $\displaystyle=$
$\displaystyle-{\Pi}_{2}\left(\frac{3}{2}\phi+\mathcal{A}\right)-\mathcal{A}\left({\mu}_{2}+p_{2}\right)-j_{e},$
(147) $\displaystyle K$ $\displaystyle=$
$\displaystyle\frac{1}{3}({\mu}_{1}+{\mu}_{2})-\mathcal{E}-\frac{1}{2}({\Pi}_{1}+{\Pi}_{2})+\frac{1}{4}{\phi}^{2},$
(148) $\displaystyle\hat{Q}_{1}$ $\displaystyle=$ $\displaystyle-
Q_{1}\left(\phi+2\mathcal{A}\right)+j_{u},$ (149) $\displaystyle Q_{2}$
$\displaystyle=$ $\displaystyle-Q_{1},$ (150)

and in terms of $\rho$, we have

$\displaystyle\phi{\phi}_{,\rho}=$
$\displaystyle-\frac{1}{2}{\phi}^{2}-\frac{2}{3}({\mu}_{1}+{\mu}_{2})-\mathcal{E}-\frac{1}{2}({\Pi}_{1}+{\Pi}_{2}),$
(151a) $\displaystyle{\mathcal{E}}_{,\rho}=$
$\displaystyle\frac{1}{3}({\mu}_{1}-{\mu}_{2})-\frac{1}{2}({\Pi}_{1}+{\Pi}_{2})-\frac{3}{2}\left(\mathcal{E}+\frac{1}{2}{\Pi}_{1}+\frac{1}{2}{\Pi}_{2}\right),$
(151b) $\displaystyle-\mathcal{A}\phi$
$\displaystyle+\frac{1}{3}({\mu}_{1}+3p_{1})+\frac{1}{3}({\mu}_{2}+3p_{2})-\mathcal{E}+\frac{1}{2}{\Pi}_{1}+\frac{1}{2}{\Pi}_{2}=0,$
(151c) $\displaystyle\phi{\mathcal{A}}_{,\rho}=$
$\displaystyle-\mathcal{A}\left(\mathcal{A}+\phi\right)+\frac{1}{2}({\mu}_{1}+3p_{1})+\frac{1}{2}({\mu}_{2}+3p_{2}),$
(151d) $\displaystyle\phi{p}_{1,\rho}+\phi{\Pi}_{1,\rho}=$
$\displaystyle-{\Pi}_{1}\left(\frac{3}{2}\phi+\mathcal{A}\right)-\mathcal{A}\left({\mu}_{1}+p_{1}\right)+j_{e}^{(1,2)},$
(151e) $\displaystyle\phi{p}_{2,\rho}+\phi{\Pi}_{2,\rho}=$
$\displaystyle-{\Pi}_{2}\left(\frac{3}{2}\phi+\mathcal{A}\right)-\mathcal{A}\left({\mu}_{2}+p_{2}\right)-j_{e}^{(1,2)},$
(151f) $\displaystyle K+\mathcal{E}-\frac{1}{4}{\phi}^{2}=$
$\displaystyle\frac{1}{3}({\mu}_{1}+{\mu}_{2})-\frac{1}{2}({\Pi}_{1}+{\Pi}_{2}),$
(151g) $\displaystyle\phi{Q}_{1,\rho}=$ $\displaystyle-
Q_{1}\left(\phi+2\mathcal{A}\right)+j_{u}^{(1,2)},$ (151h) $\displaystyle
Q_{2}=$ $\displaystyle-Q_{1}.$ (151i)

The corresponding TOV equations are given in (18).

##  APPENDIX C Full solutions for Sec. VI

###  C.1 Sec. VI.3: The case $\mathbb{P}_{1}\neq 0$

$\displaystyle p_{r1}$
$\displaystyle=\frac{1}{4a-4br^{2}}+\frac{1}{3}e^{\alpha}c_{1}\left[12a^{2}+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right]^{\frac{5-4b}{16b-8}},$
(152a) $\displaystyle p_{\perp 1}$
$\displaystyle=\left\\{-3a^{2}b[8b(26b-67)+141]r^{4}+9a^{3}[24b(2b-7)+19]r^{2}+576a^{4}+4ab^{2}[40b(2b-5)+89]r^{6}\right.$
$\displaystyle\left.-4(5-4b)^{2}b^{3}r^{8}\right\\}\times\left\\{192\left(a-br^{2}\right)^{3}\left[12a^{2}+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right]\right\\}^{-1}+\frac{1}{3}e^{\alpha}c_{1}\left[12a^{2}\right.$
$\displaystyle\left.+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right]^{\frac{5-4b}{16b-8}},$
(152b) $\displaystyle p_{r2}$
$\displaystyle=\frac{3}{4\left(a-br^{2}\right)}-\frac{1}{3}e^{\alpha}c_{1}\left[12a^{2}+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right]^{\frac{5-4b}{16b-8}},$
(152c) $\displaystyle p_{\perp 2}$
$\displaystyle=\left\\{(4br+r)^{2}\left(3a-2br^{2}\right)\left(9a^{2}-7abr^{2}+2b^{2}r^{4}\right)-64e^{\alpha}c_{1}\left(a-br^{2}\right)^{3}\left[12a^{2}+3a(1-4b)r^{2}\right.\right.$
$\displaystyle\left.\left.+2b(2b-1)r^{4}\right]^{\frac{3}{8}\left(\frac{1}{2b-1}+2\right)}-6a^{2}b[8b(26b-107)+205]r^{4}+54a^{3}[8b(2b-11)+9]r^{2}+1728a^{4}\right.$
$\displaystyle\left.+8ab^{2}[80(b-4)b+131]r^{6}-8b^{3}[16(b-4)b+37]r^{8}\right\\}\times\left\\{192\left(a-br^{2}\right)^{3}\left[12a^{2}+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right]\right\\}^{-1}.$
(152d)

###  C.2 Sec. VI.4: The case $\mathbb{P}_{1}=0$

$\displaystyle p_{r1}=p_{\perp 1}$
$\displaystyle=\frac{1}{3}e^{\beta}c_{1}\left[12a^{2}+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right]^{\mathfrak{z}_{3}},$
(153a) $\displaystyle p_{r2}$
$\displaystyle=\frac{1}{a-br^{2}}-\frac{1}{3}e^{\beta}c_{1}\left[12a^{2}+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right]^{\mathfrak{z}_{3}},$
(153b) $\displaystyle p_{\perp 2}$
$\displaystyle=\left\\{r^{2}(4br+r)^{2}\left(3a-2br^{2}\right)\left(9a^{2}-7abr^{2}+2b^{2}r^{4}\right)-48e^{\beta}c_{1}r^{2}\left(a-br^{2}\right)^{3}\left[12a^{2}\right.\right.$
$\displaystyle\left.\left.+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right]^{\frac{12b+6\gamma-3}{16b-8}}-6a^{2}b[8b(26b-107)+205]r^{6}\right.$
$\displaystyle\left.+54a^{3}[8b(2b-11)+9]r^{4}+1728a^{4}r^{2}+8ab^{2}[80(b-4)b+131]r^{8}\right.$
$\displaystyle\left.-8b^{3}[16(b-4)b+37]r^{10}\right\\}\times\left\\{144r^{2}\left(a-br^{2}\right)^{3}\left[12a^{2}+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right]\right\\}^{-1},$
(153c)

where

$\beta=\frac{\sqrt{3}(4b+1)(2\gamma+1)\tan^{-1}\left[\frac{a(3-12b)+4b(2b-1)r^{2}}{a\sqrt{48b^{2}-24b-9}}\right]}{4(2b-1)\sqrt{(4b-3)(4b+1)}},$
(154)

and we have used (101) for the case $\mathbb{P}_{1}=0$.

###  C.3 Sec. VI.4: The case $\mathbb{P}_{1}\neq 0$

$\displaystyle p_{r1}$
$\displaystyle=\frac{1}{3}e^{\beta}c_{1}\left[12a^{2}+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right]^{\mathfrak{z}_{3}},$
(155a) $\displaystyle p_{\perp 1}$
$\displaystyle=\frac{1}{12}e^{\beta}c_{1}\left[12a^{2}+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right]^{\frac{-20b+6\gamma+13}{16b-8}}\left\\{48a^{2}-3ar^{2}[4b(\gamma+4)+\gamma-4]\right.$
$\displaystyle\left.+2br^{4}[4b(\gamma+2)+\gamma-4]\right\\},$ (155b)
$\displaystyle p_{r2}$
$\displaystyle=\frac{1}{a-br^{2}}-\frac{1}{3}e^{\beta}c_{1}\left[12a^{2}+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right]^{\mathfrak{z}_{3}},$
(155c) $\displaystyle p_{\perp 2}$
$\displaystyle=\left\\{-624a^{2}b^{3}r^{4}+1608a^{2}b^{2}r^{4}+432a^{3}b^{2}r^{2}-4e^{\beta}c_{1}\left(a-br^{2}\right)^{3}\left\\{48a^{2}-3ar^{2}[4b(\gamma+4)\right.\right.$
$\displaystyle\left.\left.+\gamma-4]+2br^{4}[4b(\gamma+2)+\gamma-4]\right\\}\left[12a^{2}+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right]^{\mathfrak{z}_{3}}\right.$
$\displaystyle\left.-423a^{2}br^{4}-1512a^{3}br^{2}+171a^{3}r^{2}+576a^{4}+320ab^{4}r^{6}-800ab^{3}r^{6}+356ab^{2}r^{6}-64b^{5}r^{8}\right.$
$\displaystyle\left.+160b^{4}r^{8}-100b^{3}r^{8}\right\\}\times\left\\{48\left(a-br^{2}\right)^{3}\left[12a^{2}+3a(1-4b)r^{2}+2b(2b-1)r^{4}\right]\right\\}^{-1},$
(155d)

where we have used (101).

##  APPENDIX D Full solutions for Sec. VII: Interior Schwarzschild-Tolman IV

$\displaystyle P_{1}=$
$\displaystyle\left[c_{1}e^{\rho}\left(A^{2}+2e^{\rho}\right)\left(\frac{c_{1}+z_{1}}{c_{1}}\right)^{\frac{1}{2}(-a-1)}\left(1-\frac{z_{1}}{c_{1}}\right)^{\frac{a+1}{2}}\left(z_{1}^{2}-c_{1}^{2}\right)^{\frac{1}{2}(-a-1)}\right]\times\left[\left(A^{2}+e^{\rho}\right)\left(e^{\rho}-R^{2}\right)\right]^{-1}$
(156a) $\displaystyle P_{2}=$
$\displaystyle\frac{1}{12}\left(\frac{2e^{\rho}\left(A^{4}+A^{2}\left(2e^{\rho}+R^{2}\right)+2e^{2\rho}\right)\left(-c_{1}z_{1}+2\mu_{1}e^{\rho}-3\right)}{\left(A^{2}+e^{\rho}\right)\left(A^{2}+2e^{\rho}\right)\left(e^{\rho}-R^{2}\right)\left(-c_{1}z_{1}+\mu_{1}e^{\rho}-3\right)}\frac{R^{2}\left(A^{2}+2e^{\rho}\right)}{\left(A^{2}+e^{\rho}\right)\left(e^{\rho}-R^{2}\right)}+\frac{2\mu_{1}^{2}e^{2\rho}}{z_{1}^{2}(c_{1}+z_{1})^{2}}\right.$
$\displaystyle\left.+\frac{2\mu_{1}e^{\rho}\left(c_{1}\left(\mu_{1}e^{\rho}-6\right)-6z_{1}\right)}{\left(z_{1}^{2}\right)^{3/2}(c_{1}+z_{1})^{2}}-\frac{2\mu_{1}e^{\rho}}{c_{1}z_{1}+z_{1}^{2}}+1\right)-\frac{c_{1}e^{\rho}\left(A^{2}+2e^{\rho}\right)\left(\frac{c_{1}+z_{1}}{c_{1}}\right)^{\frac{1}{2}(-a-1)}\left(1-\frac{z_{1}}{c_{1}}\right)^{\frac{a+1}{2}}\left(z_{1}^{2}-c_{1}^{2}\right)^{\frac{1}{2}(-a-1)}}{\left(A^{2}+e^{\rho}\right)\left(e^{\rho}-R^{2}\right)}$
(156b) $\displaystyle\mathbb{P}_{2}=$
$\displaystyle\left[e^{2\rho}\left(-A^{4}\mu_{1}(c_{1}+z_{1})\left(\mu_{1}R^{2}-3\right)+A^{2}\left(c_{1}^{2}\left(z_{1}^{2}\right)^{3/2}-2\mu_{1}e^{\rho}\left(2\mu_{1}R^{2}(c_{1}+z_{1})+3c_{1}\right)+3\mu_{1}R^{2}(c_{1}+z_{1})\right.\right.\right.$
$\displaystyle\left.\left.\left.+\mu_{1}^{2}e^{2\rho}(3c_{1}+2z_{1})+9(2c_{1}+z_{1})\right)+2\left(c_{1}^{2}R^{2}\left(z_{1}^{2}\right)^{3/2}-6\mu_{1}e^{\rho}R^{2}(2c_{1}+z_{1})+\mu_{1}e^{2\rho}\left(c_{1}\mu_{1}R^{2}+3(c_{1}+z_{1})\right)\right.\right.\right.$
$\displaystyle\left.\left.\left.+9R^{2}(2c_{1}+z_{1})\right)\right)\right]\times\left[6\left(z_{1}^{2}\right)^{3/2}\left(A^{2}+e^{\rho}\right)\left(A^{2}+2e^{\rho}\right)(c_{1}+z_{1})^{2}\left(e^{\rho}-R^{2}\right)\right]^{-1},$
(156c)

where

$z_{1}=\sqrt{3-\mu_{1}e^{\rho}}$ (157)

$\displaystyle p_{r1}=p_{\perp 1}$
$\displaystyle=-\frac{4\mathfrak{c}_{1}}{R^{2}}\left(\frac{c_{1}+\mathfrak{z}_{4}}{c_{1}}\right)^{-\mathfrak{z}_{5}}\left(1-\frac{\mathfrak{z}_{4}}{c_{1}}\right)^{\mathfrak{z}_{5}}\left(\mathfrak{z}_{4}^{2}-c_{1}^{2}\right)^{-\mathfrak{z}_{5}},$
(158a) $\displaystyle p_{r2}$
$\displaystyle=-\left\\{-2r^{4}\left[-A^{2}\left\\{c_{1}^{2}\mathfrak{z}_{4}^{3}+c_{1}\left[\mu_{1}r^{2}\left(3\mu_{1}r^{2}-4\mu_{1}R^{2}-6\right)+3\mu_{1}R^{2}+18\right]\right.\right.\right.$
$\displaystyle\left.\left.\left.+\mathfrak{z}_{4}\left\\{\mu_{1}\left[2\mu_{1}r^{4}+R^{2}\left(3-4\mu_{1}r^{2}\right)\right]+9\right\\}\right\\}+A^{4}\mu_{1}(c_{1}+\mathfrak{z}_{4})\left(\mu_{1}R^{2}-3\right)-2c_{1}\mu_{1}^{2}r^{4}R^{2}\right.\right.$
$\displaystyle\left.\left.-6\mu_{1}r^{2}\left[c_{1}\left(r^{2}-4R^{2}\right)+\mathfrak{z}_{4}\left(r^{2}-2R^{2}\right)\right]-2R^{2}\left[c_{1}\left(c_{1}\mathfrak{z}_{4}^{3}+18\right)+9\mathfrak{z}_{4}\right]\right]\right.$
$\displaystyle\left.-12\mathfrak{c}_{1}r^{2}\mathfrak{z}_{4}^{3}\left(A^{2}+2r^{2}\right)^{2}(c_{1}+\mathfrak{z}_{4})^{2}\left(\frac{c_{1}+\mathfrak{z}_{4}}{c_{1}}\right)^{-\mathfrak{z}_{5}}\left(1-\frac{\mathfrak{z}_{4}}{c_{1}}\right)^{\mathfrak{z}_{5}}\left(\mathfrak{z}_{4}^{2}-c_{1}^{2}\right)^{-\mathfrak{z}_{5}}+2r^{2}\mathfrak{z}_{4}^{2}(c_{1}+\mathfrak{z}_{4})\right.$
$\displaystyle\left.\left[A^{2}\left(2r^{2}+R^{2}\right)+A^{4}+2r^{4}\right]\left(c_{1}\mathfrak{z}_{4}-2\mu_{1}r^{2}+3\right)+R^{2}\mathfrak{z}_{4}^{3}\left(A^{2}+2r^{2}\right)^{2}(c_{1}+\mathfrak{z}_{4})^{2}\right.$
$\displaystyle\left.-2\mu_{1}r^{2}\mathfrak{z}_{4}^{2}\left(A^{2}+r^{2}\right)\left(A^{2}+2r^{2}\right)(c_{1}+\mathfrak{z}_{4})(r^{2}-R^{2})+\mathfrak{z}_{4}^{3}\left(A^{2}+r^{2}\right)\left(A^{2}+2r^{2}\right)(c_{1}+\mathfrak{z}_{4})^{2}(r^{2}-R^{2})\right.$
$\displaystyle\left.+2\mu_{1}r^{2}\left(A^{2}+r^{2}\right)\left(A^{2}+2r^{2}\right)(r^{2}-R^{2})\left(c_{1}\mu_{1}r^{2}-6(c_{1}+\mathfrak{z}_{4})\right)+2\mu_{1}^{2}r^{4}\mathfrak{z}_{4}\left(A^{2}+r^{2}\right)\right.$
$\displaystyle\left.\left(A^{2}+2r^{2}\right)(r^{2}-R^{2})\right\\}\times\left\\{3r^{2}R^{2}\mathfrak{z}_{4}^{3}\left(A^{2}+2r^{2}\right)^{2}(c_{1}+\mathfrak{z}_{4})^{2}\right\\}^{-1},$
(158b) $\displaystyle p_{\perp 2}$
$\displaystyle=-\left\\{r^{4}\left[-A^{2}\left\\{c_{1}^{2}\mathfrak{z}_{4}^{3}+c_{1}\left[\mu_{1}r^{2}\left(3\mu_{1}r^{2}-4\mu_{1}R^{2}-6\right)+3\mu_{1}R^{2}+18\right]\right.\right.\right.$
$\displaystyle\left.\left.\left.+\mathfrak{z}_{4}\left\\{\mu_{1}\left[2\mu_{1}r^{4}+R^{2}\left(3-4\mu_{1}r^{2}\right)\right]+9\right\\}\right\\}+A^{4}\mu_{1}(c_{1}+\mathfrak{z}_{4})\left(\mu_{1}R^{2}-3\right)-2c_{1}\mu_{1}^{2}r^{4}R^{2}\right.\right.$
$\displaystyle\left.\left.-6\mu_{1}r^{2}\left[c_{1}\left(r^{2}-4R^{2}\right)+\mathfrak{z}_{4}\left(r^{2}-2R^{2}\right)\right]-2R^{2}\left[c_{1}\left(c_{1}\mathfrak{z}_{4}^{3}+18\right)+9\mathfrak{z}_{4}\right]\right]\right.$
$\displaystyle\left.-12\mathfrak{c}_{1}r^{2}\mathfrak{z}_{4}^{3}\left(A^{2}+2r^{2}\right)^{2}(c_{1}+\mathfrak{z}_{4})^{2}\left(\frac{c_{1}+\mathfrak{z}_{4}}{c_{1}}\right)^{-\mathfrak{z}_{5}}\left(1-\frac{\mathfrak{z}_{4}}{c_{1}}\right)^{\mathfrak{z}_{5}}\left(\mathfrak{z}_{4}^{2}-c_{1}^{2}\right)^{-\mathfrak{z}_{5}}+2r^{2}\mathfrak{z}_{4}^{2}(c_{1}+\mathfrak{z}_{4})\right.$
$\displaystyle\left.\left[A^{2}\left(2r^{2}+R^{2}\right)+A^{4}+2r^{4}\right]\left(c_{1}\mathfrak{z}_{4}-2\mu_{1}r^{2}+3\right)+R^{2}\mathfrak{z}_{4}^{3}\left(A^{2}+2r^{2}\right)^{2}(c_{1}+\mathfrak{z}_{4})^{2}\right.$
$\displaystyle\left.-2\mu_{1}r^{2}\mathfrak{z}_{4}^{2}\left(A^{2}+r^{2}\right)\left(A^{2}+2r^{2}\right)(c_{1}+\mathfrak{z}_{4})(r^{2}-R^{2})+\mathfrak{z}_{4}^{3}\left(A^{2}+r^{2}\right)\left(A^{2}+2r^{2}\right)(c_{1}+\mathfrak{z}_{4})^{2}(r^{2}-R^{2})\right.$
$\displaystyle\left.+2\mu_{1}r^{2}\left(A^{2}+r^{2}\right)\left(A^{2}+2r^{2}\right)(r^{2}-R^{2})\left(c_{1}\mu_{1}r^{2}-6(c_{1}+\mathfrak{z}_{4})\right)+2\mu_{1}^{2}r^{4}\mathfrak{z}_{4}\left(A^{2}+r^{2}\right)\left(A^{2}+2r^{2}\right)\right.$
$\displaystyle\left.(r^{2}-R^{2})\right\\}\times\left\\{3r^{2}R^{2}\mathfrak{z}_{4}^{3}\left(A^{2}+2r^{2}\right)^{2}(c_{1}+\mathfrak{z}_{4})^{2}\right\\}^{-1},$
(158c)

where

$\displaystyle\mathfrak{z}_{4}=$ $\displaystyle\sqrt{3-r^{2}\mu_{1}},$ (159a)
$\displaystyle\mathfrak{z}_{5}=$ $\displaystyle\frac{a+1}{2}.$ (159b)

The energy densities of both fluids are

$\displaystyle\mu_{1}=$
$\displaystyle-\frac{4a\mathfrak{c}_{1}}{R^{2}}\left(-c_{1}^{2}-\mu_{1}r^{2}+3\right)^{\frac{1}{2}(-a-1)}\left(\frac{c_{1}+\mathfrak{z}_{4}}{c_{1}}\right)^{\frac{1}{2}(-a-1)}\left(1-\frac{\mathfrak{z}_{4}}{c_{1}}\right)^{\frac{a+1}{2}},$
(160a) $\displaystyle\mu_{2}=$
$\displaystyle-\left[-4a\mathfrak{c}_{1}\left(A^{2}+2r^{2}\right)^{2}\left(-c_{1}^{2}-\mu_{1}r^{2}+3\right)^{\frac{1}{2}(-a-1)}\left(\frac{c_{1}+\mathfrak{z}_{4}}{c_{1}}\right)^{\frac{1}{2}(-a-1)}\left(1-\frac{\mathfrak{z}_{4}}{c_{1}}\right)^{\frac{a+1}{2}}\right.$
$\displaystyle\left.-r^{2}\left(7A^{2}+2R^{2}\right)-3A^{2}\left(A^{2}+R^{2}\right)-6r^{4}\right]\times\left\\{R^{2}\left(A^{2}+2r^{2}\right)^{2}\right\\}^{-1},$
(160b)

where we have used (159).

## References

  * [1] R. C. Tolman, Phys. Rev. 55, 364 (1939).
  * [2] J. R. Oppenheimer and G. Volkoff, Phys. Rev. 55, 374 (1939).
  * [3] M. S. R. Delgaty and K. Lake, Comput. Phys. Commun. 115, 395 (1998).
  * [4] J. Schaffner-Beilich, Compact star physics (Cambridge University Press, Cambridge, England, 2020).
  * [5] W. Israel, Nuovo Cim. B 44S10, 1 (1966) Erratum: [Nuovo Cim. B 48, 463 (1967)].
  * [6] K. Schwarzschild, Sitzungsberichte der Königlich Preussischen Akademie der Wissenschaften zu Berlin, S42, 189 (1916).
  * [7] P. C. Vaidya, Proc. Indian Acad. Sci. A 33, 264 (1951).
  * [8] L. Herrera and N. O. Santos, Phys. Rept. 286, 53-130 (1997).
  * [9] P. S. Letelier, Phys. Rev. D 22, no.4, 807 (1980).
  * [10] S. Viaggiu, Int. J. Mod. Phys. D 18, 275 (2009).
  * [11] C. Cattoen, T. Faber and M. Visser, Class. Quant. Grav. 22, 4189 (2005).
  * [12] T. Harko and M. K. Mak, Annalen Phys. 11, 3-13 (2002).
  * [13] H. Rago, Astrophys Space Sci 183, 333-338 (1991).
  * [14] P. Boonserm, M. Visser and S. Weinfurtner, Phys. Rev. D 71, 124037 (2005).
  * [15] P. Boonserm, M. Visser and S. Weinfurtner, Phys. Rev. D 76, 044024 (2007).
  * [16] S. Das, B. K. Parida, K. Chakraborty and S. Ray, arXiv pre-print (2021).
  * [17] C. C. Lovekin, Front. Astron. Space Sci. 6, 77 (2019).
  * [18] C. J. Hansen, S. D. Kawaler and V. Trimble, Stellar Modeling. In: Stellar Interiors (Springer, New York, 2004).
  * [19] D. Suárez-Urango, L. A. Núñez and H. Hernández, arXiv pre-print (2021).
  * [20] F. Tello-Ortiz, M. Malaver, A. Rincón, et al., Eur. Phys. J. C 80, 371 (2020).
  * [21] M. D. Duez and Y. Zlochower, Rep. Prog. Phys. 82, 016902 (2019).
  * [22] L. Lehner and F. Pretorius, Annu. Rev. Astron. Astrophys. 52, 661-694 (2014).
  * [23] N. F. Naidu, S. Carloni and P. Dunsby, Phys. Rev. D , 044014 (2021).
  * [24] C. A. Clarkson and R. K. Barrett, Class. Quant. Grav. 20, 3855 (2003).
  * [25] S. Carloni and D. Vernieri, Phys. Rev. D 97, 124056 (2018).
  * [26] S. Carloni and D. Vernieri, Phys. Rev. D 97, 124057 (2018).
  * [27] D. Langlois, A. Sedrakian, and B. Carter, Mon. Not. R. Astron. Soc. 297, 1189 (1998).
  * [28] S. Carloni, Phys. Rev. D 90, 044023 (2014).
  * [29] C. Barrabes and W. Israel, Phys. Rev. D 43, 1129 (1991).
  * [30] R. L. Bowers and E. P. T. Liang, Astrophys. J. 188, 657 (1974).
  * [31] P. S. Florides, Proc. Roy. Soc. Lond. A 337 529 (1974).
  * [32] C. Clarkson, Phys. Rev. D 76, 104034 (2007).
  * [33] G. Betschart and C. A. Clarkson, Class. Quant. Grav. 21, 5587 (2004).