New perspectives on the TOV equilibrium from a dual null approach

When one describes spherical stars in equilibrium it is well known that the matter distribution must satisfy the Tolman–Oppenheimer–Volkoff (TOV) equation [1, 2]. The TOV equation appear as the relativistic counterpart of the classical condition for hydrostatic equilibrium. It gives a first approximation to describe virtually any body in the sky which is large enough such that its dynamics is dominated by gravity, and stationary enough to enable us to assume it to be in static equilibrium, such as planets and stars. Due to their evident relevance, many solutions for this type of configurations have been found [3–15], and more specifically, different formalisms [16–20] and solutions generating techniques have been developed [21, 22]. Extensions of the TOV equation have also been investigated in the framework of modified gravity theories [23–29]. Yet a unified characterisation of the underlying features of the TOV equation has attracted little attention, and this is what concerns us in the present work. As it happens with the generalised perception of Birkhoff theorem [30] as being restricted to the spherically symmetric case, which is a misled assumption as shown for instance by Bona [31], and more recently discussed in [32], a similar idea is very much spread regarding the TOV equation. In the present work we aim at showing that this is a restrictive view, and that a generalised TOV equation also characterises the equilibrium of models endowed with other symmetries besides the spherical.

Focusing on GR fluid dynamics (for a review of the pioneering work on this viewpoint see for instance, and references therein [33–35]), we can trace it back to Hawking and Penrose's singularity theorems [36], and it eventually enables one to tackle 'small scale' problems with the same tools as those applicable to 'large scale' ones. In this context, an approach using the properties of light cones is most likely to reveal the structures of spacetime both at small and large scales.

The dual null formalism⁶ offers a description of the spacetime based on the properties of the optical flow. The latter is characterised by two linearly independent null congruences which are orthogonal to some codimension-two foliation of the spacetime. This approach shares some of the convenience of the choice of dual null coordinates, but it has the significant advantage of being a coordinate free formalism. Furthermore it also reveals, by construction, the causal structure of the spacetime in a natural way. It has been originally introduced to study general relativistic problems associated with the behaviour of dynamical black holes [40, 41], but it is also most convenient to analyse other diverse questions. For instance, it has been considered in connection with the definition of energy in more general geometries [42], with the gravitational collapse of fluids [43], and even with the definition of generalised horizons in modified gravity [44]. The dual null formalism has been useful to explicit the 'linear' behaviour of gravity for sources that satisfy the hypotheses of the Birkhoff theorem [32], as well.

Here, we apply the dual null formalism to analyse, in a unified way, the spacetimes which admit a codimension-two foliation with constant curvature leaves. This comprises the spherical, planar and hyperbolic symmetries, sourced by a perfect fluid. Aside from the Killing vector fields that are tangent to those surfaces of symmetry, we assume the existence of an additional symmetry generated by a Killing vector field orthogonal to those surfaces at each event. A particular case of this setup, where the symmetry is spherical and the Killing vector is timelike, corresponds to the spherically symmetric perfect fluid in hydrostatic equilibrium, which leads us to the well-known TOV equation. As we will show this celebrated equation arises most naturally in the dual-null framework which, moreover, allows us to generalise it for the planar and hyperbolic cases. To the best of our knowledge, this generalisation of the geometry underlying the TOV equation, stepping beyond spherical symmetry, has never been seen before, and leads to consequences which are far from trivial.

Since in planar and hyperbolic geometries the spatial hypersurfaces are open, this extension requires the novel introduction of a mass-energy 'parameter'. In fact one needs to promote a generalisation of the Misner-Sharp/Hawking-Hayward definitions of the mass-energy distribution, which overcomes the problem of the divergence of the latter quantities due to the natural threading with infinite surfaces.

Finally and unexpectedly, the third novelty of our explorations of the planar and hyperbolic geometries stems from the physically significant appearance of violations of the weak energy conditions in order to maintain hydrostatic equilibrium. This takes the form of negative mass, physically translating repulsive curvature effects, which suggest a link to repulsive source models, as those proposed to mimic dark energy, generate bouncing universes, or support classical wormholes [45, 46].

When the metric is characterised by a spacelike Killing vector, we have spatially homogeneous spacetimes that, as we will show, correspond to some of the Bianchi spacetimes, or Kantowski–Sachs, as expected. The hydrostatic equilibrium on those spacetimes is only possible when their source is a cosmological constant in the non flat cases, also implying the violation of energy conditions.

We proceed as follows: in section 2 we start by giving a short introduction to the dual null formalism, using from the onset the symmetries assumed in our class of problems to simplify expressions. In section 3 we then prove a proposition that states that the Killing vector two-expansion always vanishes. In order to interpret the geometrical quantities that appear, and to establish their underlying physical content, we discuss, in section 4, the mass-energy definition in these spacetimes. From the property found in section 3, we show that one derives either the equation of hydrostatic equilibrium when the Killing is timelike, in section 5.1, or the evolution equation when the Killing is spacelike, in section 5.2. Finally, in section 6 we look for planar and hyperbolic symmetric analogs of Schwarzschild interior solution, by assuming that the fluid is uncompressible and solving the unified TOV equations.

We consider metrics that have a codimension-two maximally symmetric foliation, and can be written as

$$\begin{align}\hfill \mathrm{d}{s}^{2}& ={N}_{ab}\mathrm{d}{x}^{a}\mathrm{d}{x}^{b}+{Y}^{2}\left({x}^{c}\right)\left(\mathrm{d}{\theta }^{2}+{S}_{{\epsilon}}^{2}\mathrm{d}{\phi }^{2}\right),\hfill \end{align} \tag{ 1 }$$

where

$$\begin{equation*}{S}_{{\epsilon}}=\left\{\begin{array}{c}\hfill \mathrm{sin} \theta ,\;\text{for}\;{\epsilon}=1\\ \hfill \quad 1,\;\text{for}\;{\epsilon}=0\\ \hfill \mathrm{sinh}\;\theta ,\;\text{for}\;{\epsilon}=-1\end{array}\right.,\end{equation*}$$

and where we divide the tangent space $\mathcal{T}$ at each event in two orthogonal subspaces $\mathcal{T}=\mathcal{N}\oplus \mathcal{S}$. Here $\mathcal{S}$ is the subspace generated by the orbits of (θ, ϕ) and $\mathcal{N}$, the subspace of $\mathcal{T}$ orthogonal to $\mathcal{S}$. The x^a coordinates are chosen orthogonal to $\mathcal{S}$, which gives the metric in the warped sum form of equation (1).

We denote s_ab = Y²γ_ab the induced metric in each leaf of the foliation where Y(x^c) is the warp factor. Evidently, ${\gamma }_{ab}{:=}{\delta }_{a}^{\theta }{\delta }_{b}^{\theta }+{S}_{{\epsilon}}^{2}{\delta }_{a}^{\phi }{\delta }_{b}^{\phi }$ has constant curvature and does not depend on the coordinates x^a which identify each leaf ${{\Sigma}}_{{x}^{c}}$, defined as the locus spanned by the orbits of θ and ϕ for fixed x^c. We define an orthonormal two dimensional basis (n^a, e^a) for $\mathcal{N}$, whose induced metric is N_ab, according to equation (1). This basis satisfies

$$\begin{equation}-{n}^{a}{n}_{a}={e}^{a}{e}_{a}=1,\quad {n}^{a}{e}_{a}={n}^{a}{s}_{ab}={e}^{a}{s}_{ab}=0.\end{equation} \tag{ 2 }$$

We may also define a dual null basis for the same subspace from n^a and e^a by

$$\begin{align}\hfill {k}^{a}& =\frac{1}{2}\left({n}^{a}+{e}^{a}\right),\quad {l}^{a}=\frac{1}{2}\left({n}^{a}-{e}^{a}\right),\quad {n}^{a}={k}^{a}+{l}^{a},\quad {e}^{a}={k}^{a}-{l}^{a},\hfill \end{align} \tag{ 3 }$$

which satisfies

$$\begin{equation}{k}^{a}{k}_{a}={l}^{a}{l}_{a}=0,\quad {k}^{a}{l}_{a}=-\frac{1}{2}.\end{equation} \tag{ 4 }$$

The metric g_ab can be written as

$$\begin{equation}{g}_{ab}=\frac{2}{{k}^{c}{l}_{c}}{k}_{\left(\right.a}{l}_{b\left.\right)}+{s}_{ab}.\end{equation} \tag{ 5 }$$

We associate the null expansion for each null vector as follows

$$\begin{equation}{{\Theta}}_{k}=\frac{1}{2}{s}^{ab}{\mathcal{L}}_{k}{s}_{ab}=\frac{1}{2}{Y}^{-2}{\gamma }^{ab}{\mathcal{L}}_{k}{Y}^{2}{\gamma }_{ab}=\frac{2}{Y}{k}^{a}{\partial }_{a}Y.\end{equation} \tag{ 6 }$$

We may extend the definition of null expansion to timelike and spacelike vectors in $\mathcal{N}$, calling it the two-expansion, since it measures the rate of variation of area, as in the null case. We may define the mean curvature form ${\mathcal{K}}_{a}={\partial }_{a}\;\mathrm{ln}\;{Y}^{2}$, such that, we obtain for the two-expansion Θ_(u) of any vector u^a in $\mathcal{N}$

$$\begin{equation}{{\Theta}}_{\left(u\right)}={u}^{a}{\mathcal{K}}_{a}.\end{equation} \tag{ 7 }$$

We describe our spacetimes by means of the behaviour of the null expansion, casting the Einstein equations, G_ab = 8πT_ab in terms of expansions, i.e., by writing the Raychaudhuri equations [36, 40, 47, 48]. In the latter equations T_ab is, as usual, the energy–momentum tensor T_ab, and we express it as a general fluid under a 1 + 1 + 2 decomposition along n^a and e^a [49] thus reading

$$\begin{equation}{T}_{ab}=\left(\rho +P\right){n}_{a}{n}_{b}+P{g}_{ab}+2{q}_{\left(\right.a}{n}_{b\left.\right)}+{{\Pi}}_{ab}.\end{equation} \tag{ 8 }$$

In the latter decomposition ρ is the energy density, and P is the isotropic pressure, both measured by an observer moving with four-velocity n^a, ${q}^{a}=q{e}^{a}+{\mathfrak{q}}^{a}$ is the heat flow vector, decomposed into its scalar part along e^a and a two-vector ${\mathfrak{q}}^{a}$ on the maximally symmetric surfaces, and, finally, Π_ab is the anisotropic stress tensor. Π_ab is also decomposed as ${{\Pi}}_{ab}={\Pi}{P}_{ab}+2{{\Pi}}_{\left(\right.a}{e}_{b\left.\right)}+{{\Pi}}_{{\left\langle ab\right\rangle }_{s}}$ into a scalar part along the flow-orthogonal symmetric tensor P_ab = s_ab − 2e_ae_b, which is traceless, into a vector part semi-orthogonal to the maximally symmetric surface with inner vector Π^a, and into a fully embedded tensor ${{\Pi}}_{{\left\langle ab\right\rangle }_{s}}$.

In terms of the metric (5), the energy–momentum tensor can be expressed as

$$\begin{align}\hfill {T}_{ab}=& \left[\rho +P+2\left(q-{\Pi}\right)\right]{k}_{a}{k}_{b}+\left[\rho +P-2\left(q+{\Pi}\right)\right]{l}_{a}{l}_{b}\hfill \\ \hfill & +2{k}_{\left(\right.a}{l}_{b\left.\right)}\left[\rho -P+2{\Pi}\right]+\left[P-{\Pi}\right]{s}_{ab}+2{\mathfrak{q}}_{\left(\right.a}{\left[k+l\right]}_{b\left.\right)}\hfill \\ \hfill & +2{{\Pi}}_{\left(\right.a}{\left[k-l\right]}_{b\left.\right)}+{{\Pi}}_{{\left\langle ab\right\rangle }_{s}},\hfill \end{align} \tag{ 9 }$$

and thus its projected components yield

$$\begin{equation}{T}_{ab}{k}^{a}{k}^{b}=\frac{1}{4}\left[\rho +P-2\left(q+{\Pi}\right)\right],\end{equation} \tag{ 10 }$$

$$\begin{equation}{T}_{ab}{k}^{a}{l}^{b}=\frac{1}{4}\left[\rho -P+2{\Pi}\right],\end{equation} \tag{ 11 }$$

$$\begin{equation}{T}_{ab}{l}^{a}{l}^{b}=\frac{1}{4}\left[\rho +P+2\left(q-{\Pi}\right)\right].\end{equation} \tag{ 12 }$$

The Raychaudhuri and constraint equations then read

$$\begin{equation}{\mathcal{L}}_{k}{{\Theta}}_{\left(k\right)}={\nu }_{k}{{\Theta}}_{\left(k\right)}-\frac{{{\Theta}}_{\left(k\right)}^{2}}{2}-8\pi {T}_{ab}{k}^{a}{k}^{b},\end{equation} \tag{ 13a }$$

$$\begin{equation}{\mathcal{L}}_{l}{{\Theta}}_{\left(l\right)}={\nu }_{l}{{\Theta}}_{\left(l\right)}-\frac{{{\Theta}}_{\left(l\right)}^{2}}{2}-8\pi {T}_{ab}{l}^{a}{l}^{b},\end{equation} \tag{ 13b }$$

$$\begin{align}\hfill {\mathcal{L}}_{k}{{\Theta}}_{\left(l\right)}+{\mathcal{L}}_{l}{{\Theta}}_{\left(k\right)}=-{{\Theta}}_{\left(l\right)}{\nu }_{k}-{{\Theta}}_{\left(k\right)}{\nu }_{l}-2\;{{\Theta}}_{\left(k\right)}\;{{\Theta}}_{\left(l\right)}+{\epsilon}\frac{2\;{k}^{a}{l}_{a}}{{Y}^{2}}+16\pi \;{T}_{ab}{k}^{a}{l}^{b},\end{align} \tag{ 13c }$$

where we included the inaffinities ν_k and ν_l, defined as

$$\begin{equation}{\nu }_{k}=\frac{1}{{k}^{c}{l}_{c}}{l}^{b}{k}^{a}{\nabla }_{a}{k}_{b}\quad \quad {\nu }_{l}=\frac{1}{{k}^{c}{l}_{c}}{k}^{b}{l}^{a}{\nabla }_{a}{l}_{b}.\end{equation} \tag{ 14 }$$

In this work we adapt our vector basis to a fluid source, such that n^a gives its flow. Therefore, it will be useful to relate our quantities to this flow. By construction, the flow n^a is always orthogonal to the surfaces of symmetry and will be characterised by two quantities

$$\begin{equation}\mathcal{A}={e}^{a}{\dot {n}}_{a}={e}^{a}{n}^{b}{\nabla }_{b}{n}_{a},\quad \quad \mathcal{B}={e}^{a}{n}_{a}^{\prime }={e}^{a}{e}^{b}{\nabla }_{b}{n}_{a}.\end{equation} \tag{ 15 }$$

The scalar $\mathcal{A}$ gives us the acceleration of the flow, a positive sign meaning that the acceleration is outwards in the spherical, compact case. The scalar $\mathcal{B}$ gives the change of direction of n^a as we travel along e^a. It is the e–e component of the extrinsic curvature K_ab of the three-space orthogonal to this flow, since

$$\begin{equation}{K}_{ab}=\frac{1}{2}{\mathcal{L}}_{n}{h}_{ab}\,,\end{equation} \tag{ 16 }$$

where h_ab = g_ab + n_an_b. We may also write

$$\begin{equation}{h}_{ab}={e}_{a}{e}_{b}+{Y}^{2}{\gamma }_{ab}\,,\end{equation} \tag{ 17 }$$

which gives

$$\begin{equation}{K}_{ab}=\mathcal{B}\;{e}_{a}{e}_{b}+\frac{{{\Theta}}_{\left(n\right)}}{2}\;{Y}^{2}{\gamma }_{ab}\,.\end{equation} \tag{ 18 }$$

The trace of equation (18) gives us the flow of the volumetric expansion ${{\Theta}}_{3}={\nabla }_{a}{n}^{a}={K}_{a}^{a}$ as

$$\begin{equation}{{\Theta}}_{3}=\mathcal{B}+{{\Theta}}_{\left(n\right)}\end{equation} \tag{ 19 }$$

In order to relate our quantities with the flow scalars, we compute the shear scalar σ, by taking the traceless part of K_ab. We obtain

$$\begin{equation}\sigma =\frac{{{\Theta}}_{\left(n\right)}}{6}-\frac{\mathcal{B}}{3},\end{equation} \tag{ 20 }$$

which implies

$$\begin{equation}\frac{{{\Theta}}_{3}}{3}+\sigma =\frac{{{\Theta}}_{\left(n\right)}}{2},\end{equation} \tag{ 21 }$$

in agreement with the result obtained in reference [43].

Using the inaffinities of the null basis vectors, $\mathcal{A}$ and $\mathcal{B}$ can be expressed as

$$\begin{equation}\mathcal{A}={\nu }_{k}-{\nu }_{l},\quad \quad \mathcal{B}={\nu }_{k}+{\nu }_{l}.\end{equation} \tag{ 22 }$$

We now assume that our metric has a Killing vector orthogonal to maximally symmetric surfaces. Our symmetry requirements imply that it commutes with the symmetry generators on the foliation. We denote this hypersurface orthogonal Killing vector field χ^a. It satisfies the Killing equation,

$$\begin{equation}{\mathcal{L}}_{\chi }{g}_{ab}=0.\end{equation} \tag{ 23 }$$

Proposition 3.1. If a spacetime is described by a metric of the form (1) and admits an orthogonal Killing vector ${\chi }^{a}\in \mathcal{N}$, then Θ_χ = 0.

Proof. We may write, from equation (6), ${{\Theta}}_{\chi }=\frac{1}{2}{s}^{ab}{\mathcal{L}}_{\chi }{s}_{ab}$$=\frac{1}{2}{Y}^{-2}{\gamma }^{ab}{\mathcal{L}}_{\chi }{Y}^{2}{\gamma }_{ab}$, and

$$\begin{equation}{g}_{ab}={N}_{ab}+{Y}^{2}{\gamma }_{ab}.\end{equation} \tag{ 24 }$$

Then

$$\begin{align}\hfill 0& ={Y}^{-2}{\gamma }^{ab}{\mathcal{L}}_{\chi }{g}_{ab}={Y}^{-2}{\gamma }^{ab}{\mathcal{L}}_{\chi }{N}_{ab}+2{{\Theta}}_{\chi }=-{Y}^{-2}{N}^{ab}{\mathcal{L}}_{\chi }{\gamma }_{ab}+2{{\Theta}}_{\chi }.\hfill \end{align} \tag{ 25 }$$

However

$$\begin{equation}{\mathcal{L}}_{\chi }{\gamma }_{ab}=0,\end{equation} \tag{ 26 }$$

since χ^a does not admit components in $\mathcal{S}$ and γ_ab does not depend on coordinates along $\mathcal{N}$. Therefore, equation (25) implies that Θ_χ = 0. □

Consequently, if there is an extra symmetry with orbits orthogonal to those of the maximally symmetric leaves of the foliation, the two-expansion of its generator vanishes. This also implies that if dY is spacelike, then χ_a is timelike and vice-versa. If dY is null, the Killing vector will also be null.

In order to properly interpret our spacetimes, we have to understand their mass-energy content. There is a widely known mass-energy definition suitable to the spherically symmetric case, namely the Misner–Sharp mass-energy [50, 51], defined regardless of asymptotic assumptions. However, as we also intend to analyse nonspherical spacetimes in this work, we are lead to a more general mass-energy definition such as the Hawking–Hayward's one (hereafter HH)[42, 52]. The HH mass-energy gives the mass-energy content inside a closed compact surface in terms of an integral over that surface, in a manner similar to the Gauss law in Newtonian gravity. This quasilocal mass-energy has been explored in different contexts, such as seen in references [53–55] .

In the case where the ${{\Sigma}}_{{x}^{c}}$ are spheres, which are compact, we can compute the Hawking–Hayward mass-energy enclosed by Σ (we drop the x^c index for short) as

$$\begin{align}\hfill & \quad \quad \quad {M}_{{\Sigma}}=\frac{1}{8\pi }\sqrt{\frac{A}{16\pi }}{\int }_{{\Sigma}}\left[\mathcal{R}-\frac{1}{{k}^{a}{l}_{a}}\left({{\Theta}}_{\left(k\right)}{{\Theta}}_{\left(l\right)}-\frac{1}{2}{\sigma }_{\left(k\right)\;ab}{\sigma }^{\left(l\right)\;ab}-2{\omega }_{a}{\omega }^{a}\right)\right]\mathrm{d}{\Sigma}\hfill \end{align} \tag{ 27 }$$

where $\mathcal{R}$ is the two-dimensional Ricci scalar and A is the area of Σ, σ_(k)ab and σ^(l)ab are the two-dimensional shear tensors along Σ, associated with the k and l congruences, respectively, and ω^a is the twist vector given by the projection on Σ of the commutator of the null basis vectors. We have included a factor of $-\frac{1}{{k}^{a}{l}_{a}}$ in the optical scalars part of the mass-energy, compared with the formula present in reference [42], in order to take account of our different normalization of the null normals.

Our symmetry assumptions imply that the only non-vanishing optical scalar on Σ is the null expansion. Therefore, the Hawking–Hayward mass-energy is reduced to

$$\begin{equation}{M}_{{\Sigma}}=\frac{1}{8\pi }\sqrt{\frac{A}{16\pi }}{\int }_{{\Sigma}}\left[\mathcal{R}-\frac{1}{{k}^{a}{l}_{a}}{{\Theta}}_{\left(k\right)}{{\Theta}}_{\left(l\right)}\right]\;\mathrm{d}{\Sigma}\end{equation} \tag{ 28 }$$

Since we assume that Σ is maximally symmetric, we have $\mathcal{R}=\frac{2{\epsilon}}{{Y}^{2}}$. We also have

$$\begin{align}\hfill & \hfill {{\Theta}}_{\left(k\right)}{{\Theta}}_{\left(l\right)}={k}^{a}{\partial }_{a}\;\mathrm{ln}\;{Y}^{2}\;{l}^{b}{\partial }_{b}\;\mathrm{ln}\;{Y}^{2}={k}^{a}{l}^{b}{\partial }_{a}\;\mathrm{ln}\;{Y}^{2}{\partial }_{b}\;\mathrm{ln}\;{Y}^{2}\\ \hfill & ={k}^{\left(\right.a}{l}^{b\left.\right)}{\partial }_{a}\;\mathrm{ln}\;{Y}^{2}{\partial }_{b}\;\mathrm{ln}\;{Y}^{2}=\frac{{k}^{c}{l}_{c}}{2}{g}^{ab}{\partial }_{a}\;\mathrm{ln}\;{Y}^{2}{\partial }_{b}\;\mathrm{ln}\;{Y}^{2}\hfill \\ \hfill & \hfill \quad \quad \quad =\frac{1}{2}{k}^{c}{l}_{c}{{\Vert}\mathrm{d}\;\mathrm{ln}\;{Y}^{2}{\Vert}}^{2}={k}^{c}{l}_{c}\frac{2}{{Y}^{2}}{{\Vert}\mathrm{d}Y{\Vert}}^{2}{\Rightarrow}\frac{{{\Theta}}_{\left(k\right)}{{\Theta}}_{\left(l\right)}}{{k}^{c}{l}_{c}}=\frac{2}{{Y}^{2}}{{\Vert}\mathrm{d}Y{\Vert}}^{2}\,,\end{align} \tag{ 29 }$$

where we used equation (5) in the fourth step.

For the spherical case = 1 and A = 4πY², we obtain the known interpretation of ||dY|| in terms of the Misner–Sharp mass-energy, which coincides with the Hawking–Hayward one

$$\begin{equation}{M}_{{\Sigma}}=\frac{Y}{2}\left(1-{\Vert}\mathrm{d}Y{{\Vert}}^{2}\right){\Leftrightarrow}{\Vert}\mathrm{d}Y{{\Vert}}^{2}=1-\frac{2M}{Y}.\end{equation} \tag{ 30 }$$

In the planar and hyperbolic cases ( = 0 and = −1, respectively), the Hawking–Hayward mass is not conveniently defined for the integration domain set by our preferred foliation, as it requires a closed compact surface.

In this work, we aim to treat all three symmetry types in the same manner. Therefore, we need to find a mass-energy definition which might be equivalent to the HH mass-energy, but suitable to deal with non compact domains in order to take advantage of the planar or hyperbolic symmetry. We can make such an extension of the HH mass-energy, as long as their boundary correspond to a pair Σ of symmetric two-surfaces of symmetry corresponding to the same warp factor Y. Of course, those domains are infinite and have an infinite mass-energy content in general. However, as they are homogeneous along the surfaces of symmetry, we can successfully adapt the HH mass-energy definition in order to obtain a finite mass-energy parameter for those cases. They then describe an infinite mass-energy distribution, homogeneous along the surfaces of symmetry, with a finite density.

We proceed by first making the replacement

$$\begin{equation}\frac{1}{8\pi }\sqrt{\frac{A}{16\pi }}\to \frac{Y}{4\pi \kappa },\end{equation} \tag{ 31 }$$

in order to keep its dimensionality, and eliminating the explicit dependence on the area of Σ. Evidently, by setting κ = 4 we recover the Hawking–Hayward mass-energy in the spherical case. This step is justified by the fact that originally this factor was introduced to correct the dimensionality of the mass-energy, and to make it match the Arnowitt–Deser–Misner (ADM) mass [56], where both are well defined. Since our symmetric spacetimes allow an 'areal scalar' as the warp factor Y, we can replace $\sqrt{A}$ by Y as the quantity with dimension of length associated to each surface of symmetry.

We then define the quasi-local mass-energy parameter μ(Y) by

$$\begin{equation}{M}_{{\Sigma}}=\frac{\mu \left(Y\right)}{4\pi }\int {S}_{{\epsilon}}^{2}\left(\theta \right)\mathrm{d}\theta \mathrm{d}\phi ,\end{equation} \tag{ 32 }$$

and we write

$$\begin{align}\hfill & \frac{Y}{4\pi \kappa }\left[\mathcal{R}-\frac{{{\Theta}}_{\left(k\right)}{{\Theta}}_{\left(l\right)}}{{k}^{c}{l}_{c}}\right]{\int }_{{\Sigma}}\mathrm{d}{\Sigma}=\frac{Y}{4\pi \kappa }\left[2{\epsilon}-2{\Vert}\mathrm{d}Y{{\Vert}}^{2}\right]\int {S}_{{\epsilon}}^{2}\left({\Theta}\right)\mathrm{d}\theta \mathrm{d}\phi .\hfill \end{align} \tag{ 33 }$$

We equate equations (32) and (33) and eliminate the improper area integral on both sides

$$\begin{equation}{\Vert}\mathrm{d}Y{{\Vert}}^{2}={\epsilon}-\frac{\kappa \mu \left(Y\right)}{2Y}.\end{equation} \tag{ 34 }$$

An alternative route to equation (34) can be obtained by computing the HH mass-energy in a finite domain, symmetric with respect to the central plane or wire, Y = 0, and taking the limit where the domain tends to be the whole surface. The finite integration domain consist of the union of

• (a)  
  a subset of the Σ_Y, that we denote Γ_r, bounded by a circle γ_r of radius r on the (θ, ϕ) coordinate plane and
• (b)  
  a compact surface given by the surfaces Δ_r defined by γ_r transported along Y orbits.

It forms a closed surface, corresponding to a part of a cylinder bounded by Y = constant surfaces in the space of coordinates (Y, θ, ϕ). Therefore, the HH mass-energy enclosed by those surfaces will by finite, and given by

$$\begin{equation}{M}_{r}=\frac{1}{8\pi }\sqrt{\frac{{A}_{r}}{16\pi }}\left({\int }_{{{\Gamma}}_{r}}\left(\dots \;\right){S}_{{\epsilon}}^{2}\mathrm{d}\theta \mathrm{d}\phi +{\int }_{{{\Delta}}_{r}}\left(\dots \;\right)\mathrm{d}{\Delta}\right)\end{equation} \tag{ 35 }$$

where (...) replaces the integrand of equation (27). In the limit r → ∞, the first integral in equation (35) scales as r² while the second one scales as r. This means that, in the limit r → ∞, and repeating the replacement in equation (31) we obtain

$$\begin{equation}\frac{{M}_{r}}{{A}_{r}}\to \frac{\mu \left(Y\right)}{4\pi {Y}^{2}}\end{equation} \tag{ 36 }$$

Equation (34) coincides with the known mass function (see equation (15.7a) in [57]) which appears as we integrate the Einstein equations of specific spacetimes with metrics of the form (1) for planar and hyperbolic symmetries. Hereafter, we will consider equation (34) with the choice κ = 4 as the mass-energy definition.

5.1. Timelike Killing vector

We assume χ^aχ_a < 0. In this case, the spacetime is static, and n^a ∼ χ^a. Therefore, from proposition 3.1, Θ_(n) = 0 everywhere, and dY is spacelike, since it is orthogonal to n^a. If Θ_(n) vanishes everywhere, this means that the fluid has no radial velocity, therefore we are dealing with a static fluid with a flow parallel to the Killing vector field.

In order to characterise its static equilibrium, we need to compute the derivative of the flow two-expansion along the flow itself:

$$\begin{equation}{\mathcal{L}}_{n}{{\Theta}}_{\left(n\right)}=0,\end{equation} \tag{ 37 }$$

since Θ_(n) = 0 everywhere.

We may write ${\mathcal{L}}_{n}{{\Theta}}_{\left(n\right)}$ in terms of the null expansions as

$$\begin{equation}{\mathcal{L}}_{n}{{\Theta}}_{\left(n\right)}={\mathcal{L}}_{k}{{\Theta}}_{\left(k\right)}+{\mathcal{L}}_{l}{{\Theta}}_{\left(l\right)}+{\mathcal{L}}_{l}{{\Theta}}_{\left(k\right)}+{\mathcal{L}}_{k}{{\Theta}}_{\left(l\right)}.\end{equation} \tag{ 38 }$$

Substituting the equations (13a), (13b), and (13c), we obtain

$$\begin{align}\hfill {\mathcal{L}}_{n}{{\Theta}}_{\left(n\right)}=-\frac{{\left({{\Theta}}_{\left(k\right)}+{{\Theta}}_{\left(l\right)}\right)}^{2}}{2}-{{\Theta}}_{\left(k\right)}{{\Theta}}_{\left(l\right)}-\frac{{\epsilon}}{{Y}^{2}}-8\pi {T}_{ab}{e}^{a}{e}^{b}+\mathcal{A}\left({{\Theta}}_{\left(k\right)}-{{\Theta}}_{\left(l\right)}\right).\end{align} \tag{ 39 }$$

Recall that we are assuming Θ_(n) = Θ_(k) + Θ_(l) = 0, and that Θ_(k) −Θ_(l) = Θ_(e), using equations (7) and (3). We identify here Θ_(k)Θ_(l) as the mass term, since it equals $\frac{2}{{Y}^{2}}{\Vert}\mathrm{d}Y{{\Vert}}^{2}=\frac{2}{{Y}^{2}}\left({\epsilon}-\frac{2\mu \left(Y\right)}{Y}\right)$.

Taking the source to be a perfect fluid, then the energy momentum tensor (8) reduces to

$$\begin{equation}{T}_{ab}=\rho {n}_{a}{n}_{b}+P\left({e}_{a}{e}_{b}+{s}_{ab}\right).\end{equation} \tag{ 40 }$$

Contracting the conservation of the energy–momentum tensor with e_b (Euler equation in [58]) we get

$$\begin{align}\hfill {e}_{b}{\nabla }_{a}{T}^{ab}& =\left(\rho +P\right){\dot {n}}^{b}{e}_{b}+{e}^{a}{\nabla }_{a}P=0{\Rightarrow}\mathcal{A}=-\frac{{e}^{a}{\partial }_{a}P}{\rho +P}.\hfill \end{align} \tag{ 41 }$$

Since Θ_(n) = 0, this implies that e^a is proportional to ∂_Y, and as e^a is normalized, we have ${e}_{a}=\frac{1}{{\Vert}\mathrm{d}Y{\Vert}}{\partial }_{a}Y$. Imposing e^ae_a = 1 we obtain

$$\begin{equation}{e}^{a}={\Vert}\mathrm{d}Y{\Vert}{\left({\partial }_{Y}\right)}^{a},\end{equation} \tag{ 42 }$$

which gives us

$$\begin{equation}\mathcal{A}{{\Theta}}_{e}=-{\Vert}\mathrm{d}Y{{\Vert}}^{2}\;\frac{2}{Y}\frac{{\partial }_{Y}P}{\rho +P}.\end{equation} \tag{ 43 }$$

Therefore, replacing ||dY||² by its meaning in terms of mass, ${\mathcal{L}}_{n}{{\Theta}}_{\left(n\right)}=0$ corresponds to

$$\begin{equation}\left({\epsilon}-\frac{2\mu \left(Y\right)}{Y}\right)\frac{2}{Y}\frac{{\partial }_{Y}P}{\rho +P}=-2\frac{\mu \left(Y\right)}{{Y}^{3}}-8\pi P,\end{equation} \tag{ 44 }$$

or, alternatively,

$$\begin{equation}\frac{{\partial }_{Y}P}{\rho +P}=-\left(\frac{\mu \left(Y\right)}{{Y}^{2}}+4\pi PY\right){\left({\epsilon}-\frac{2\mu \left(Y\right)}{Y}\right)}^{-1},\end{equation} \tag{ 45 }$$

which is what we call the unified TOV equation. It reduces to the well-known TOV equation for spherically symmetric spacetimes when ɛ = 1, and it corresponds to the equation of hydrostatic equilibrium for planar and hyperbolic geometries, in the cases where ɛ = 0 and ɛ = −1, respectively. This underlines the fact that the TOV equation is a hydrostatic equilibrium equation, and not an equation of state, as it is erroneously stated sometimes.

In order to determine μ(Y) we consider the ${\mathcal{L}}_{e}{{\Theta}}_{\left(e\right)}$ Raychaudhuri equation

$$\begin{align}\hfill {\mathcal{L}}_{e}{{\Theta}}_{\left(e\right)}& =\mathcal{B}{{\Theta}}_{\left(n\right)}-\frac{{{\Theta}}_{\left(e\right)}^{2}}{2}+\frac{1}{4}\left({{\Theta}}_{\left(n\right)}^{2}-{{\Theta}}_{\left(e\right)}^{2}\right)+\frac{{\epsilon}}{{Y}^{2}}-8\pi {T}_{ab}{n}^{a}{n}^{b},\hfill \end{align} \tag{ 46 }$$

which, by using Θ_(n) = 0, and equation (42) leads us to

$$\begin{equation}{\Vert}\mathrm{d}Y{\Vert}{\partial }_{Y}\left(\frac{2}{Y}{\Vert}\mathrm{d}Y{\Vert}\right)=-\frac{3}{{Y}^{2}}{\Vert}\mathrm{d}Y{{\Vert}}^{2}+\frac{{\epsilon}}{{Y}^{2}}-8\pi \rho .\end{equation} \tag{ 47 }$$

Substituting equation (34) into equation (47), we obtain

$$\begin{equation}{\partial }_{Y}\mu =4\pi \rho {Y}^{2},\end{equation} \tag{ 48 }$$

which looks like the mass-energy equation of spherical symmetry. Here, it should be interpreted as the mass-energy equation in the spherical case, and as a mass-energy parameter equation in the planar and hyperbolic cases. Furthermore, equations (48) and (34) imply that if the weak energy condition (WEC) [59] holds, only the spherically symmetric case admits static regular solutions. Indeed, as those solutions require ||dY||² > 0, that implies μ < 0 for ⩽ 0 and, as in regular spacetimes,

$$\begin{equation}\mu \left(Y\right)=4\pi {\int }_{0}^{Y}\rho \left(y\right){y}^{2}\mathrm{d}y,\end{equation} \tag{ 49 }$$

this imposes ρ < 0.

With equation (48), the last requirement to solve equation (45) is the equation of state of the fluid, f(ρ, P) = 0 which should come from specific physical considerations.

5.2. Spacelike Killing vector

In the spacelike Killing vector case, dY is timelike, the flow n_a is orthogonal to the Killing vector, and the unitary base vector e^a is parallel to it. This imposes no constraint on the sign of µ according to equation (34), and thus there is no need to violate energy conditions in order to consider these solutions, thoroughly studied in cosmology [60].

One dynamical equation is given by ${\mathcal{L}}_{e}{{\Theta}}_{\left(e\right)}=0$, which according to equation (46) gives:

$$\begin{equation}\mathcal{B}{{\Theta}}_{\left(n\right)}-\frac{{{\Theta}}_{\left(e\right)}^{2}}{2}+\frac{1}{4}\left({{\Theta}}_{\left(n\right)}^{2}-{{\Theta}}_{\left(e\right)}^{2}\right)+\frac{{\epsilon}}{{Y}^{2}}-8\pi \rho =0.\end{equation} \tag{ 50 }$$

From proposition 3.1, we have Θ_(e) = 0. Replacing equation (20) in equation (46), we obtain

$$\begin{equation}\frac{3}{4}{{\Theta}}_{\left(n\right)}^{2}-3\sigma {{\Theta}}_{\left(n\right)}+\frac{{\epsilon}}{{Y}^{2}}=8\pi \rho ,\end{equation} \tag{ 51 }$$

and, using equation (21), we can express this equation (51) in terms of the volume expansion Θ₃ obtaining

$$\begin{equation}\frac{{{\Theta}}_{3}^{2}}{3}-3{\sigma }^{2}=8\pi \rho -\frac{{\epsilon}}{{Y}^{2}},\end{equation} \tag{ 52 }$$

which corresponds to the generalised Friedmann constraint equation for the evolution of a homogeneous and anisotropic universe. In the case σ = 0, we may identify Θ₃ = 3H, and we recover the usual Friedmann equation for the flat ( = 0) and open ( = 1) spatially isotropic universes. Notice though that σ = 0 also yields anisotropic, cosmological solutions when the matter content is not a perfect fluid [61].

The ${\mathcal{L}}_{n}{{\Theta}}_{\left(n\right)}$ Raychaudhuri equation gives the evolution of Θ_(n). Using equation (39)

$$\begin{align}\hfill {\mathcal{L}}_{n}{{\Theta}}_{\left(n\right)}& =-\frac{{{\Theta}}_{\left(n\right)}^{2}}{2}-{{\Theta}}_{\left(k\right)}{{\Theta}}_{\left(l\right)}-\frac{{\epsilon}}{{Y}^{2}}\;-8\pi {T}_{ab}{e}^{a}{e}^{b}{\Rightarrow}\hfill \\ \hfill {\mathcal{L}}_{n}{{\Theta}}_{\left(n\right)}& =-\frac{3}{4}{{\Theta}}_{\left(n\right)}^{2}-\frac{{\epsilon}}{{Y}^{2}}-8\pi {T}_{ab}{e}^{a}{e}^{b}.\hfill \end{align} \tag{ 53 }$$

Subtracting equation (51) from equation (53), we obtain

$$\begin{equation}{\mathcal{L}}_{n}{{\Theta}}_{\left(n\right)}=-3\sigma {{\Theta}}_{\left(n\right)}-8\pi \left(\rho +P\right),\end{equation} \tag{ 54 }$$

which, together with an equation of state relating ρ and P closes our system. By adding half of the equation (54) with one third of equation (52), we obtain:

$$\begin{align}\hfill {\mathcal{L}}_{n}\left(\frac{{{\Theta}}_{3}}{3}\right)& +{\left(\frac{{{\Theta}}_{3}}{3}\right)}^{2}=-2\sigma \left(\frac{{{\Theta}}_{3}}{3}+\sigma \right)-\frac{{\epsilon}}{3{Y}^{2}}-\frac{4\pi }{3}\left(\rho +3P\right).\hfill \end{align} \tag{ 55 }$$

Those homogeneous and anisotropic spacetimes belong to a subclass of Bianchi models [62, 63], with the case = 0 corresponding to Bianchi type I universes, = −1 corresponding to the Bianchi type III models, and = 1 to the Kantowski–Sachs spacetimes [64].

An extensive classification and evolution analysis of this family of spacetimes is certainly worthy of interest, and a great deal of work has already been carried out in the literature in this connection [60]. However, in this work we are mainly focused on the hydrostatic equilibrium situations. Thus, our interest will be directed to understanding whether it is possible to find a correspondence between the TOV equation of static equilibrium, and some condition applying to the spatially homogeneous models.

Imposing staticity amounts in the present case to have Θ₃ = 0, ${\mathcal{L}}_{n}{{\Theta}}_{3}=0$, and σ = 0 in equations (52) and (55). Reconciling the reduced equations simply requires

$$\begin{equation}\rho +P=0,\end{equation} \tag{ 56 }$$

in the = ±1 cases, and has no realisation when = 0, as ρ = 0 from equation (52). Hence we conclude that the TOV condition interpreted as a cornerstone of stability yields the well-known equation of state characterising a cosmological constant in the non flat cases. In hindsight, one could have anticipated this result, which emerges here in a self-consistent way. Moreover we see that strong energy condition (SEC) is violated for both cases = ±1, whilst the WEC is additionally violated for the = −1 case, as follows from equation (52).

Using our unified TOV equation, equation (45), we may look for static perfect fluid solutions for all three symmetries considered here. By choosing a timelike coordinate T along the flow, making n_a = −α(Y)dT, and the warp factor Y, we obtain the following line element in the (T, Y) coordinates:

$$\begin{equation}\mathrm{d}{s}^{2}=-{\alpha }^{2}\left(Y\right)\mathrm{d}T+\frac{\mathrm{d}{Y}^{2}}{{\epsilon}-\frac{2\mu \left(Y\right)}{Y}}+{Y}^{2}\mathrm{d}{{\Omega}}_{{\epsilon}},\end{equation} \tag{ 57 }$$

where $\mathrm{d}{{\Omega}}_{{\epsilon}}=\left(\mathrm{d}{\theta }^{2}+{S}_{{\epsilon}}^{2}\mathrm{d}{\phi }^{2}\right)$ and the functions α and μ will be given by solving Einstein equations, i.e., equations (45) and (34).

Here, we will apply our unified treatment to find the analogs of Schwarzschild interior solution, that is, we will use the equation of state of an incompressible fluid ρ = ρ₀ constant. It is important to note that, as we have discussed in section 5.1, the static solutions with ≠ 1 violate the WEC, therefore we should take ρ₀ < 0 in those cases.

Equation (41) implies

$$\begin{equation}\frac{{\alpha }^{\prime }}{\alpha }=-\frac{{P}^{\prime }}{\rho +P}{\Rightarrow}\alpha =\frac{{c}_{0}}{{\rho }_{0}+P},\end{equation} \tag{ 58 }$$

where c₀ is an integration constant that can be set by rescaling the time coordinate and the prime denotes Y differentiation.

Equation (48) gives us

$$\begin{equation}\mu \left(Y\right)=\frac{4\pi {\rho }_{0}{Y}^{3}}{3},\end{equation} \tag{ 59 }$$

which we replace in equation (45) to obtain

$$\begin{equation}P\left(Y\right)={\rho }_{0}\left(\frac{2\sqrt{\vert {\epsilon}-\frac{{Y}_{\mathrm{s}}}{{Y}_{\mathrm{g}}}\vert }}{3\sqrt{\vert {\epsilon}-\frac{{Y}_{\mathrm{s}}}{{Y}_{\mathrm{g}}}\vert }-\sqrt{\vert {\epsilon}-\frac{{Y}_{\mathrm{s}}{Y}^{2}}{{Y}_{\mathrm{s}}^{3}}\vert }}-1\right).\end{equation} \tag{ 60 }$$

where Y_g is the analog of the radius of the object and is the least positive number that satisfy P(Y_g) = 0, ${Y}_{\mathrm{s}}=\frac{8\pi {\rho }_{0}{Y}_{\mathrm{g}}^{3}}{3}$ is the analog of the Schwarzschild radius, although it cannot be interpreted as a location since it will be a negative number. This gives

$$\begin{equation}\alpha =\frac{1}{2}\left(3\sqrt{\left\vert {\epsilon}-\frac{{Y}_{\mathrm{s}}}{{Y}_{\mathrm{g}}}\right\vert }-\sqrt{\left\vert {\epsilon}-\frac{{Y}_{\mathrm{s}}{Y}^{2}}{{Y}_{\mathrm{g}}^{3}}\right\vert }\right)\end{equation} \tag{ 61 }$$

which has a similar form to the interior Schwarzschild solution, where we only change the sign of the mass-energy parameter and change the value of in the formula. Of course the physical properties are very distinct, since the solutions violate the WEC.

In figure 1 we compare the pressure for the three cases. From the slope of the curves, we notice that only the hyperbolic case presents P' > 0, compensating the repulsive gravity force in this setup. This is the opposite of the more familiar situations presented in the spherical and planar cases, where gravity is attractive, with P' < 0 sustaining the weight of the configuration. We can also see that the planar case admits a positive pressure for 0 < Y < Y_g. That means that, as long as mass-energy is negative, we may have static plane configurations over a finite Y interval. On the other hand, the hyperbolic solution only admits positive pressure for Y > Y_g, so there is no analog of the Schwarzschild interior solution for this foliation, although it can be interpreted as an exterior fluid solution to an internal void. It can thus be matched to a hyperbolic vacuum solution for Y < Y_g, found as a particular case in reference [32]:

$$\begin{equation}\mathrm{d}{s}^{2}=-\left(\frac{2\;m}{Y}-1\right)\mathrm{d}{t}^{2}+\frac{\mathrm{d}{Y}^{2}}{\frac{2\;m}{Y}-1}+{Y}^{2}\left(\mathrm{d}{\theta }^{2}+{\mathrm{sinh}}^{2}\;\theta \mathrm{d}{\phi }^{2}\right),\end{equation} \tag{ 62 }$$

where the parameter m = |μ|. The peculiarities of the hyperbolic solutions with regard to the energy conditions are also found in one of the coauthors' work [46].

Zoom In Zoom Out Reset image size

The equation of state consisting of a negative energy density with a positive pressure might be achieved by some kind of phantom field, but a Lagrangian description of the fluid is beyond the scope of this work. However, our simple incompressible model, with constant energy density, but varying pressure, is reminiscent of a constant time surface of a McVittie or Shah-Vaidya spacetime, which admits a Lagrangian scalar field as source [65, 66].This suggests the possibility that there exist field models in the literature which can source the solutions presented in this paper.

We notice that the planar solutions may also represent a subclass of cylindrical solutions (see references [67, 68]) if we select one coordinate along the plane to be periodic. Thus our planar solution may be interpreted as a static cylinder of fluid with boundary given by Y = Y_g. At this surface, the solution must be matched with a vacuum solution.

Actually, all fluid configurations found can be matched with the corresponding static solutions presented in reference [32] which arise from applying Birkhoff theorem for external fluid sources that satisfy the hypotheses of the theorem. In those cases there is matter content present outside, the most common examples being an electromagnetic field, and a cosmological constant. Therefore the matching surface will correspond to a surface where P(Y) matches the pressure of the exterior solution, in a manner similar to the way in which an uncompressible charged sphere is matched to a Reissner–Nordström solution in [69].

In this paper we analysed spacetimes with a two-dimensional maximally symmetric foliation sourced by a perfect fluid. We proved that in those cases, if there is a Killing vector orthogonal to the leaves, its two-expansion vanishes, which allows us to simplify our dynamical equations in terms of the two-expansion of a unitary vector orthogonal to the Killing field.

When the Killing vector χ^a is timelike, we find that the flow lines must be tangent to χ^a, and as this is true at all times, the equations describe a hydrostatic equilibrium, governed by a (generalised) TOV equation. When the Killing vector is spacelike, we have instead a spatially homogeneous dynamical spacetime. The result is a subclass of Bianchi universes, with only one shear degree of freedom. The corresponding equation gives the evolution of expansion and shear scalars.

We have discussed the geometric meaning of the mass-energy in such spacetimes, and our procedure matches the traditional mass parameter found in those cases by usual methods of integration of Einstein equations. Our approach relates the mass parameter to the geometrically defined quasi-local mass-energies of Misner–Sharp and Hawking–Hayward by slightly changing its definition in order to apply it to our infinite mass-energy distributions. This innovation is in itself a step towards addressing the open issue of defining mass/energy in gravitation and cosmology, cf the recent works of [70], and others [71, 72] on this subject.

Using these concepts we could recover the physical interpretation of the geometrical quantities appearing in the equilibrium/evolution equations, translate the dual null formalism to the more usual relativistic fluid dynamics framework, and show that the TOV equation arises as a particular case of those equations. Henceforth the generalisations of the TOV equation appear automatically by just setting = 0 or −1 accordingly.

From this treatment it emerges the fact that the only static fluid solutions that satisfy the WEC are the spherical ones, as the other two cases require a negative energy density.

In what regards the spatially homogeneous spacetimes, the hydrostatic equilibrium condition also implies a violation of SEC for the non planar solutions, constraining the equation of state for the perfect fluid to be that of a cosmological constant.

In order to illustrate the analogy between the planar, hyperbolic and spherical cases we studied the static solutions for an incompressible fluid. We found that, besides the known case of spherically symmetric spacetimes, we can obtain a static interior fluid configuration only in the case of planar symmetric spacetimes. In the hyperbolic case, the static configuration is an exterior solution that can surround an inner vacuum region.

Our unified way to describe three classes of spacetimes foliated with codimension-two leaves of constant curvature leads the way to further generalisations, as those spacetimes are of interest in many domains, from braneworld models to AdS/CFT duality. The adaptation of our formalism to N-dimensional spacetimes is straightforward.

The introduction of the mass parameter to generalise Hawking–Hayward's mass also may have impact on other studies of compact objects in open backgrounds. In addition, we discovered that some of the models require sources that violate energy conditions. Some popular models for modified gravity theories that aim to explain large scale phenomena, such as cosmological inflation and dark energy, also violate energy conditions [73, 74], and there are several arguments in the literature suggesting those conditions should be abandoned as a criterion of viability [75].

One of the consequences of matter sources with such equations of state is that solutions with very different and intriguing properties arise, even considering the simplifying assumptions we made in order to obtain analytic results. Seeking a Lagrangian model in the modified gravity literature that could generate our fluid solutions appears as an interesting continuation of our findings.

AM wishes to thank the hospitality of Instituto de Astrofísica e Ciências do Espaço (IA), at the FCUL in Lisbon, where a part of this work was carried out. The work of MLeD has been supported by Lanzhou University starting fund and by the Fundamental Research Funds for the Central Universities (Grant No.lzujbky-2019-25). The work of JPM was supported by FCT/MCTES through national funds (PIDDAC) by this grant UID/FIS/04434/2019, and by the project PTDC/FIS-OUT/29048/2017. The authors wish to thank Xu Yumeng for helpful discussions.