1 Introduction

Spherical shear-free self-gravitating systems of perfect fluids are perhaps the most well studied systems in the historical quest for exact solutions of the Einstein field equations, and as rightly pointed out by [1], is rich in rediscoveries. All the known solutions are presented in [2,3,4,5], while in all later works these solutions are always contained as special cases. However, it is not particularly clear how these solutions physically exist due to the fact that in the Einstein field equations the properties of shear, inhomogeneity and tidal effects are intertwined. Generally, there are three different approaches to find exact solutions of the field equations for shear-free and spherical perfect fluid distributions. The first approach is via an adhoc ansatz for one of the metric functions, while the second one is to look for symmetries including Lie point symmetries, contact symmetries and Noether symmetries in the field equations. The third approach is more general and rigorous, where the solutions have Painlevé properties [4]. As a matter of fact, all the known solutions belong to this latter class.

The quest to relate the three features of the shear, the Weyl tensor and inhomogeneity has contributed to a rich tapestry of research over the last 40 years. One of the first attempts to connect these attributes is the Penrose conjecture which assumes that the Weyl tensor [6], or some functions of it [7,8,9], gives rise to a gravitational arrow of time. This assumption is informed by the evolution of inhomogeneity through the effects of tidal forces on the gravitating fluid and was ultimately met with some doubt [10,11,12]. In particular, in [10] two arrows of time, one arising from radiation processes and the other being the common definition associated with the Weyl and Ricci tensors, were observed to point in opposite senses, making the assumption unsatisfactory.

Thereafter investigations on the interplay of these attributes for a general fluid distribution with various symmetry configurations were published in several papers by Herrera and co-authors. These configurations include spherical symmetry [13,14,15], cylindrical symmetry [16], and static [17] and dynamic [18] axial symmetry. In [13], dissipation and pressure isotropy were discussed in detail relative to the effect it has on the relationship between the Weyl tensor and inhomogeneity. A general dissipative fluid, geodesic fluids, and non-dissipative fluids with either isotropic or anisotropic pressure were analyzed in [13]. In particular we are interested in non-dissipative isotropic fluids for this paper. From [13] we learned that for locally isotropic perfect fluids, the vanishing of the Weyl tensor and density inhomogeneity are equivalent and either one of them would imply the shear-free condition. This is not the case for locally anisotropic fluids, both dissipative and non-dissipative. It was further shown in [15] that a scalar function, arising from the trace-free part of the double dual of the Riemann tensor, is responsible for the inhomogeneity in non-dissipative fluids [15] and is a combination of the Weyl tensor (scalar part) and pressure anisotropy. These results are summarized succinctly in [19].

Along with dissipation and pressure anisotropy, it was shown in [14] that electric charge also affects the relationship between the Weyl tensor and inhomogeneity. For the cylindrically symmetric case [16], the analysis uncovered that the shear-free condition is also related to the magnetic part of the Weyl tensor but the cause of inhomogeneity was indecipherable. Whereas in the axially symmetric case [17], it was possible to identify the inhomogeneity factor which are scalars and scalar functions relating to the electric part of the Weyl tensor and pressure anisotropy.

Contributing to more links between the three attributes, in [20] the evolution of the shear-free condition was found to be closely related to the absence of dissipation and pressure anisotropy. Additionally in [21], the three features and dissipative flux are shown to generate pressure anisotropy even if isotropy is initially assumed. All arguments considered, the aforementioned Penrose conjecture relating inhomogeneity solely to the Weyl tensor clearly does not hold. As demonstrated in the literature, there exists a profound connection between the shear, the Weyl tensor, pressure isotropy and dissipative fluxes. In this paper, yet another interesting connection is added to the list which involves the equation of state. Our focussed interest lies in determining the physical reason why a self-gravitating, shear-free, inhomogeneous system with tidal effects exists at all and that is what we present.

In this paper, we investigate this well studied system through a more geometrical perspective. In fact, we generalize the spacetime geometry to locally rotationally symmetric (LRS-II) spacetimes [22], of which spherical symmetry is a subclass. By covariantly decomposing the spacetime using the fluid flow congruence and a preferred spherical congruence which is guaranteed by LRS-II symmetry, we classify all the possible solutions of the field equations in terms of the covariantly defined fluid acceleration and the fluid expansion.

Shear-free matter distributions have been used to model both static and radiating stars, justifying its assumption. Spherical symmetry simplifies the relativistic equations even further. Spherically symmetric shear-free solutions have been used to model many physical applications leading to many classes of solutions [23,24,25,26,27,28,29,30] where it is demonstrated that zero shear restricts the possible choices of equations of state [25, 26]. Recent applications of this model relate to gravitational collapse [31,32,33,34]; symmetry analysis [35, 36] and methods [37, 38]; and particular solutions [39,40,41].

Interestingly, our analysis sheds new light on the class of solutions where fluid acceleration and expansion are strictly non-vanishing. It can easily be seen that the Weyl curvature is also non-vanishing for this class. This class is interesting as it is dynamical and necessarily inhomogeneous but still shear-free. It is well known that spacetime shear and the electric part of the Weyl form a feedback loop with each other. The electric part of the Weyl is a source term for the shear evolution equation while the shear is the source term for the evolution of the electric part of the Weyl.

Therefore, constraining one of the positive feedback loop quantities to vanish identically, while the other is strictly non-vanishing, would imply a stringent constraint on other geometrical and thermodynamical quantities in the spacetime.

From our analysis, we transparently show that indeed such a stringent constraint exists on the possible equations of state that give rise to this class. We find the governing ODE, which is highly nonlinear, the solution to which gives the possible classes of equations of state that may give rise to this class of solutions. We present a numerical solution for this ODE and clearly show that the usual equations of state that are linear, or a finite combination of power law solutions, will not generate the given class.

The paper is organized as follows: In the next two sections we discuss the covariant semitetrad decomposition of the spacetime manifold, and the field equations written in terms of geometrical and thermodynamical variables that emerge due to the decomposition. In the subsequent two sections we discuss the special case of shear-free spacetimes, and show how all possible solutions can be classified in terms of the fluid acceleration and the fluid expansion. The next section is dedicated to the interesting class of inhomogeneous and shear-free spacetimes, and the differential equation that governs the equation of state of the perfect fluid. Finally, we conclude our results in the last section.

2 1+1+2 formalism

The 1+1+2 covariant approach, first introduced by Greenberg in [42], further reviewed by van Elst and Ellis [22] and subsequently expanded by Clarkson and Barrett [43, 44] is an extension of the 1+3 formalism [45, 46]. The latter is one of the most widely used tetrad approaches which formulates the equations of general relativity, using a timelike congruence, as first order differential equations unlike the coordinate approach involving second order partial derivatives of the metric functions.

In the 1+3 covariant approach the essential ‘time’ coordinate is separated from the 3-‘space’ so that the metric tensor takes the form

$$\begin{aligned} g_{ab}=h_{ab}-u_au_b. \end{aligned}$$
(1)

Two important derivatives are defined. The covariant time derivative along the observers’ worldlines, denoted by ‘ \({^{\cdot }}\) ’, is defined using the vector \({u^{a}}\), as

$$\begin{aligned} \dot{Z}^{a \ldots b}{}_{c \ldots d} = u^{e}\nabla _{e} Z^{a \ldots b}{}_{c \ldots d}, \end{aligned}$$
(2)

for any tensor \({Z^{a\ldots b}{}_{c\ldots d}}\). The fully orthogonally projected covariant spatial derivative, denoted by ‘ D ’, is defined using the spatial projection tensor \({h_{ab}}\), as

$$\begin{aligned} D_{e} Z^{a\ldots b}{}_{c\ldots d} = h^r{}_{e} h^p{}_{c}\ldots h^q{}_{d} h^a{}_{f}\ldots h^b{}_{g}\nabla _{r} Z^{f\ldots g}{}_{p\ldots q}, \end{aligned}$$
(3)

with total projection on all the free indices. In the context of locally rotationally symmetric class II (LRS-II) spacetimes [22], the geometrical quantities defined for the timelike congruence are the expansion scalar \((\Theta = D_{a} u^{a} )\), the acceleration 3-vector \((\dot{u}^a = u^{b}\nabla _{b}u^a )\) and the shear 3-tensor \([\sigma _{ab} = \left( h^{c}{}_{(a} h^{d}{}_{b)} - \frac{1}{3}h_{ab} h^{cd}\right) D_{c} u_{d}]\). The timelike congruence uniquely defines the electric part of the Weyl tensor \((E_{ab} = C_{acbd} u^{c} u^{d} = E_{<ab>})\) and the magnetic part vanishes identically. The energy momentum tensor is also decomposed according to the timelike congruence to produce the energy density \((\mu = T_{ab} u^{a} u^{b})\), the isotropic pressure \((p = \frac{1}{3} h_{ab} T^{ab})\), the heat flux 3-vector \((q_a = q_{<a>} = - h^{c}{}_{a} T_{cd} u^{d})\) and the anisotropic stress 3-tensor \((\pi _{ab} = T_{cd} h^c{}_{<a} h^d{}_{b>})\). Angle brackets denote orthogonal projections of covariant time derivatives along \({u^{a}}\) as well as represent the projected, symmetric and trace-free part of tensors as follows

$$\begin{aligned} v_{<a>} = h^{b}{}_{a}\dot{v}_{b}, \quad \quad Z_{<ab>} = \left( h^{c}{}_{(a} h^{d}{}_{b)} - \frac{1}{3}h_{ab} h^{cd}\right) Z_{cd}. \end{aligned}$$
(4)

All these quantities have a direct geometrical and physical meaning and are described by tensorial quantities that remain valid in all coordinate systems which is the principal advantage of using spacetime decomposition.

A further decomposition of the 1+3 quantities gives rise to the useful 1+1+2 set of variables featured in [43]. The splitting of the LRS-II spacetime is performed with respect to the timelike unit vector \(u^a\) (\(u^au_a = -1\)) as well as the isolation of a preferred spatial direction \(e^a (e^ae_a = 1)\) chosen orthogonal to \(u^a (u^a e_a = 0)\). The metric tensor (1) is decomposed further into

$$\begin{aligned} g_{ab}=-u_a u_b+e_a e_b+N_{ab}, \end{aligned}$$
(5)

where \(N_{ab}\) is the 2-dimensional metric on the spherical 2-shell. We introduce two new derivatives for any tensor \({\Psi _{a\ldots b}{}^{c\ldots d}}\):

$$\begin{aligned} \hat{\Psi }_{a\ldots b}{}^{c\ldots d}\equiv & {} e^{f} D_{f} \Psi _{a\ldots b}{}^{c\ldots d}, \end{aligned}$$
(6)
$$\begin{aligned} \delta _{f}\Psi _{a\ldots b}{}^{c\ldots d}\equiv & {} N_{f}{}^{j} N_{a}{}^{l}\ldots N_{b}{}^{g} N_{h}{}^{c}\ldots N_{i}{}^{d} D_{j}\Psi _{l\ldots g}{}^{h\ldots i}, \end{aligned}$$
(7)

defined by the spatial congruence \({e^{a}}\). The hat-derivative (6) is the spatial derivative along the \({e^{a}}\) vector field in the surfaces orthogonal to \({u^{a}}\), and the delta-derivative (7) is the projected spatial derivative onto the 2-sheet, with projection on every free index. In a similar way, there are geometrical quantities generated for the preferred spatial congruence with the only non-vanishing quantity for LRS-II spacetimes being the 2-sheet volume expansion \((\phi = \delta _a e^a)\). Using the preferred spatial congruence we can then extract a set of covariant scalars that completely govern the dynamics of the system in the following manner

$$\begin{aligned} \mathcal {A} = \dot{u}^a e_a, \quad \Sigma = \sigma _{ab}e^a e^b, \quad \mathcal {E} = E_{ab}e^a e^b, \quad Q = q^a e_a, \quad \Pi = \pi _{ab}e^a e^b. \end{aligned}$$
(8)

Evidently the 1+1+2 decomposition method is well suited to spacetimes that have a preferred spatial direction such as LRS-II spacetimes where \(u^a\) and \(e^a\) are hypersurface orthogonal. These spacetimes have the inherent property that there exists a unique preferred spacial direction at each point that creates a local axis of symmetry. Hence all the physics and geometry of the spacetime are described by well defined kinematic and dynamic scalar variables that generate the field equations. The 1+1+2 covariant method has generated new results in studies involving LRS spacetimes [47], the Kerr spacetime [48], the Vaidya spacetime [49] and a general spacetime admitting conformal symmetry [50].

3 Field equations for LRS-II spacetimes

As pointed out by Goswami and Ellis [51], the only variables in the 1+1+2 formalism which characterize the kinematics are

$$\begin{aligned} \mathcal {D} = \left\{ \mathcal {A},\Theta ,\phi ,\Sigma ,\mathcal {E},\mu ,p,\Pi ,Q\right\} . \end{aligned}$$
(9)

For these variables the propagation, evolution and constraint equations can be derived as

$$\begin{aligned} \hat{\phi }= & {} -\frac{1}{2}\phi ^2+\left( \frac{1}{3}\Theta +\Sigma \right) \left( \frac{2}{3}\Theta - \Sigma \right) -\frac{2}{3}\mu -\frac{1}{2}\Pi -\mathcal {E}, \end{aligned}$$
(10)
$$\begin{aligned} \hat{\Sigma }-\frac{2}{3}\hat{\Theta }= & {} -\frac{2}{3}\phi \Sigma -Q, \end{aligned}$$
(11)
$$\begin{aligned} \hat{\mathcal {E}}-\frac{1}{3}\hat{\mu }+\frac{1}{2}\hat{\Pi }= & {} -\frac{3}{2}\phi \left( \mathcal {E} +\frac{1}{2}\Pi \right) +\left( \frac{1}{2}\Sigma -\frac{1}{3}\Theta \right) Q, \end{aligned}$$
(12)
$$\begin{aligned} \dot{\phi }= & {} -\left( \Sigma -\frac{2}{3}\Theta \right) \left( \mathcal {A}-\frac{1}{2}\phi \right) +Q, \end{aligned}$$
(13)
$$\begin{aligned} \dot{\Sigma }-\frac{2}{3}\dot{\Theta }= & {} -\mathcal {A}\phi +2\left( \frac{1}{3}\Theta -\frac{1}{2}\Sigma \right) ^2 -\mathcal {E} +\frac{1}{3}(\mu +3p)+\frac{1}{2}\Pi , \end{aligned}$$
(14)
$$\begin{aligned} \dot{\mathcal {E}}-\frac{1}{3}\dot{\mu }+\frac{1}{2}\dot{\Pi }= & {} \left( \frac{3}{2}\Sigma -\Theta \right) \mathcal {E}+\frac{1}{4}\left( \Sigma -\frac{2}{3}\Theta \right) \Pi +\frac{1}{2}\phi Q \nonumber \\{} & {} -\frac{1}{2}(\mu +p)\left( \Sigma -\frac{2}{3}\Theta \right) , \end{aligned}$$
(15)
$$\begin{aligned} \dot{\mu }+\hat{Q}= & {} -\Theta (\mu +p) -\frac{3}{2}\Sigma \Pi -(\phi +2\mathcal {A})Q, \end{aligned}$$
(16)
$$\begin{aligned} \dot{Q}+\hat{p}+\hat{\Pi }= & {} -\left( \frac{3}{2}\phi +\mathcal {A}\right) \Pi -(\mu +p)\mathcal {A} -\left( \frac{4}{3}\Theta +\Sigma \right) Q, \end{aligned}$$
(17)
$$\begin{aligned} \hat{\mathcal {A}}-\dot{\Theta }= & {} -(\mathcal {A}+\phi )\mathcal {A}+\frac{1}{3}\Theta ^2 +\frac{3}{2}\Sigma ^2+\frac{1}{2}\left( \mu +3p\right) , \end{aligned}$$
(18)

according to [44]. The field equations are not a closed set of equations because there are no explicit equations for \(\dot{\mathcal {A}}\) (evolution equation for \(\mathcal {A}\)) and \(\hat{\Theta }\) (propagation equation for \(\Theta \)). Hence an equation of state, governing the thermodynamical quantities, is needed of the general form \(F(\mu ,p,\Pi ,Q)=0\). This form becomes simplified for a shear-free perfect fluid model where there are no shear stresses or heat flux \((\Pi = Q =0\)). Hence our equation of state takes the form

$$\begin{aligned} p=p(\mu ), \end{aligned}$$
(19)

where the isotropic pressure p is a function of \(\mu \), the effective energy density. At this stage, we need to point out that the derivative operators ‘\(\delta ^a\)’, ‘\({^\cdot {}}\)’ and ‘ \(\hat{}\) ’ do not commute and give rise to an interesting result later. Instead, according to [44], the commutation relation when acting on some scalar \(\beta \) for our imposed LRS-II spacetime conditions is given by

$$\begin{aligned} \hat{\dot{\beta }}-\dot{\hat{\beta }}=-\mathcal {A}\dot{\beta }+ \left( \frac{1}{3}\Theta +\Sigma \right) \hat{\beta }, \end{aligned}$$
(20)

which will be used in our calculations in the subsequent sections.

4 Shear-free perfect fluid LRS-II equations

Thus far and later in this paper, we consider a shear-free perfect fluid LRS-II spacetime with an equation of state (19). Perfect fluids were studied in the context of the 1+3 decomposition method [22] and for the 1+1+2 decomposition method [44]. Such a fluid appears to be a good description of the observed universe on a large scale. The absence of shear stresses, viscosity and heat conduction is a great advantage as the relativistic equations become simpler. For this shear-free perfect fluid LRS-II model we have

$$\begin{aligned} \Sigma =Q=\Pi =0. \end{aligned}$$
(21)

Hence the above general LRS-II equations (10)–(18) simplify to

$$\begin{aligned} \hat{\phi }= & {} -\frac{1}{2}\phi ^2+\frac{2}{9}\Theta ^2-\frac{2}{3}\mu -\mathcal {E}, \end{aligned}$$
(22)
$$\begin{aligned} \hat{\Theta }= & {} 0, \end{aligned}$$
(23)
$$\begin{aligned} \hat{\mathcal {E}}-\frac{1}{3}\hat{\mu }= & {} -\frac{3}{2}\phi \mathcal {E}, \end{aligned}$$
(24)
$$\begin{aligned} \dot{\phi }= & {} \frac{2}{3}\Theta \mathcal {A}-\frac{1}{3}\Theta \phi , \end{aligned}$$
(25)
$$\begin{aligned} -\frac{2}{3}\dot{\Theta }= & {} -\mathcal {A}\phi +\frac{2}{9}\Theta ^2+\frac{1}{3}(\mu +3p)-\mathcal {E}, \end{aligned}$$
(26)
$$\begin{aligned} \dot{\mathcal {E}}-\frac{1}{3}\dot{\mu }= & {} -\Theta \mathcal {E}+\frac{1}{3}\Theta (\mu +p), \end{aligned}$$
(27)
$$\begin{aligned} \dot{\mu }= & {} -\Theta (\mu +p), \end{aligned}$$
(28)
$$\begin{aligned} \hat{p}= & {} -\mathcal {A}(\mu +p), \end{aligned}$$
(29)
$$\begin{aligned} \hat{\mathcal {A}}-\dot{\Theta }= & {} -(\mathcal {A}+\phi )\mathcal {A}+\frac{1}{3}\Theta ^2+\frac{1}{2}(\mu +3p). \end{aligned}$$
(30)

We immediately note that due to our imposed conditions (21) we obtain an explicit expression for \(\hat{\Theta }\) in (23); however, we still do not have an explicit expression for \(\dot{\mathcal {A}}\). At this point, we make \(\mathcal {A}\) the central focus of our study and use it to classify shear-free perfect fluid LRS-II solutions in the next section.

5 Classes of shear-free perfect fluid LRS-II solutions

For the equations outlined in the preceding section, we obtain various classes of shear-free perfect fluid LRS-II solutions characterized by the acceleration scalar \(\mathcal {A}\). First, it is pertinent to consider the simplest case where acceleration is absent.

5.1 LRS-II spacetime with no acceleration: \(\mathcal {A} = 0\)

In this case, we immediately see from Eq. (29) that the pressure must be homogeneous with \(\hat{p}=0\). Now comparing Eqs. (26) and (30), implies \( \mathcal {E}=0\). Equation (24) then implies \(\hat{\mu }=0\), and hence this scenario describes a dynamic and homogeneous matter distribution. Since we have \(\mathcal {A} = \mathcal {E} = 0\) and the matter is homogeneous, the spacetime has to strictly be Friedmann–Lemaître–Robertson–Walker (FLRW). It is also important to note that from the general LRS-II field equations (10)–(18), we cannot have a perfect fluid solution where \(\Theta = \mathcal {A}=0\), so we do not consider that solution branch.

5.2 LRS-II spacetime with acceleration: \(\mathcal {A} \ne 0\)

From the Eqs. (22)-(30), we noted earlier that we do not have an explicit equation for \(\dot{\mathcal {A}}\). However with the imposed equation of state (19), we can obtain this crucial evolution equation so that we have a closed set of equations, and we now describe this process in detail.

Considering the definition of the isentropic speed of sound \(c_s^2 = \left( \partial p/ \partial \mu \right) _{s = \text {constant}}\), we note that we can write

$$\begin{aligned} \hat{p}= & {} c_s^2 \hat{\mu } = \frac{\partial p}{\partial \mu }\,\hat{\mu }, \end{aligned}$$
(31)
$$\begin{aligned} \dot{p}= & {} c_s^2\dot{\mu } =\frac{\partial p}{\partial \mu }\,\dot{\mu }. \end{aligned}$$
(32)

Therefore, from Eq. (29),

$$\begin{aligned} \mathcal {A}= & {} -\frac{\hat{p}}{(\mu +p)},\nonumber \\ \implies \dot{\mathcal {A}}= & {} -\frac{\dot{\hat{p}}}{(\mu +p)}+\mathcal {A}\Theta \left( 1+\frac{\partial p}{\partial \mu }\right) , \end{aligned}$$
(33)

using (28) and (32) for simplification. Obtaining an expression for \(\dot{\hat{p}}\) from the commutation relation (20) and substituting into Eq. (33) we get

$$\begin{aligned} \dot{\mathcal {A}}=\mathcal {A}\Theta \left( \frac{\partial p}{\partial \mu }-\frac{1}{3}-\frac{(\mu +p)}{c_s^2}\frac{\partial ^2 p}{\partial \mu ^2}\right) =\mathcal {A}\Theta \mathcal {F}, \end{aligned}$$
(34)

after simplification using (32) and where we have set

$$\begin{aligned} \mathcal {F} \equiv \mathcal {F}(\mu ) = \left( \frac{\partial p}{\partial \mu }-\frac{1}{3}-\frac{(\mu +p)}{c_s^2}\frac{\partial ^2 p}{\partial \mu ^2}\right) . \end{aligned}$$
(35)

We now investigate how the obtained acceleration evolution equation (34) is compatible with the system from the integrability condition; to achieve this we make use of the commutation relation (20). According to the left hand side of Eq. (20), we need to obtain \(\hat{\dot{\mathcal {A}}}\) and \(\dot{\hat{\mathcal {A}}}\). Using (34), we have

$$\begin{aligned} \hat{\dot{\mathcal {A}}}= & {} \hat{\mathcal {A}}\Theta \mathcal {F}+\mathcal {A}\hat{\Theta }\mathcal {F}+\mathcal {A}\Theta \hat{\mathcal {F}}\nonumber \\= & {} \Theta \mathcal {F}\left( -\mathcal {A}^2+\frac{1}{2}\mathcal {A}\phi +\frac{3}{2}\mathcal {E}\right) -\mathcal {A}^2\Theta \frac{(\mu +p)}{c_s^2}\mathcal {F}', \end{aligned}$$
(36)

where \(\mathcal {F}'=\frac{\partial \mathcal {F}(\mu )}{\partial \mu }\). Using (26) in the Raychaudhuri equation (30) we get

$$\begin{aligned} \hat{\mathcal {A}}= & {} -\mathcal {A}^2+\frac{1}{2}\mathcal {A}\phi +\frac{3}{2}\mathcal {E}, \nonumber \\ \implies \dot{\hat{\mathcal {A}}}= & {} -2\mathcal {A}^2\Theta \mathcal {F}+\frac{1}{3}\mathcal {A}^2\Theta -\frac{1}{6}\mathcal {A}\Theta \phi +\frac{1}{2}\mathcal {A}\Theta \phi \mathcal {F}-\frac{3}{2}\Theta \mathcal {E}. \end{aligned}$$
(37)

Then subtracting (37) from (36), as required by the commutation relation (20), yields

$$\begin{aligned} \hat{\dot{\mathcal {A}}}-\dot{\hat{\mathcal {A}}}= & {} \mathcal {A}^2\Theta \mathcal {F}+\frac{3}{2}\Theta \mathcal {E}\mathcal {F} -\mathcal {A}^2\Theta \frac{(\mu +p)}{c_s^2}\mathcal {F}' -\frac{1}{3}\mathcal {A}^2\Theta +\frac{1}{6}\mathcal {A}\Theta \phi +\frac{3}{2}\Theta \mathcal {E}.\nonumber \\ \end{aligned}$$
(38)

However according to the commutation relation (20), we expected to obtain

$$\begin{aligned} \hat{\dot{\mathcal {A}}}-\dot{\hat{\mathcal {A}}}= & {} -\mathcal {A}\dot{\mathcal {A}}+\frac{1}{3}\Theta \hat{\mathcal {A}} \nonumber \\= & {} -\mathcal {A}^2\Theta \mathcal {F}-\frac{1}{3}\mathcal {A}^2\Theta +\frac{1}{6}\mathcal {A}\Theta \phi +\frac{1}{2}\Theta \mathcal {E}. \end{aligned}$$
(39)

We immediately notice the interesting discrepancy between (38) and (39). They do not match and this needs to be resolved for the existence of solutions. The condition for \(\dot{\mathcal {A}}\) to be compatible with the system is obtained by setting the difference between (38) and (39) to zero as follows

$$\begin{aligned} 2\mathcal {A}^2\Theta \mathcal {F}+\frac{3}{2}\Theta \mathcal {E}\mathcal {F}-\mathcal {A}^2\Theta \frac{(\mu +p)}{c_s^2}\mathcal {F}'+\Theta \mathcal {E} = 0, \end{aligned}$$
(40)

which simplifies to

$$\begin{aligned} \Theta \left[ \mathcal {A}^2\left( 2\mathcal {F}-\frac{(\mu +p)}{c_s^2}\mathcal {F}'\right) +\mathcal {E}\left( \frac{3}{2}\mathcal {F}+1\right) \right] = 0. \end{aligned}$$
(41)

Equation (41) implies that either \(\Theta = 0\) or the square-bracketed term is zero.

5.2.1 Case 1: static fluid \(\left( \Theta = 0\right) \)

Since we have \(\Sigma = Q =\Pi =\Theta = 0\) then the spacetime is necessarily static. This means that any reasonable equation of state of the form \(p=p(\mu )\) will solve the system.

5.2.2 Case 2: nonstatic fluid \(\left( \Theta \ne 0\right) \)

In this case \(\Theta \) is strictly not equal to zero. Hence the square-bracketed term in (41) has to be zero, that is

$$\begin{aligned} \left[ \mathcal {A}^2\left( 2\mathcal {F}-\frac{(\mu +p)}{c_s^2}\mathcal {F}'\right) +\mathcal {E}\left( \frac{3}{2}\mathcal {F}+1\right) \right] =0. \end{aligned}$$
(42)

Two possible subcases arise.

  • Subcase 1: If both \(\mathcal {A}=\mathcal {E}=0\) then the spacetime is necessarily Friedmann-Lemaître-Robertson-Walker (FLRW). This type reduces to the case contained in section 5.1.

  • Subcase 2: If \(\mathcal {A}\) and \(\mathcal {E}\) are well defined and both nonzero then from (41) we have

    $$\begin{aligned}{} & {} \mathcal {A}^2\left( 2\mathcal {F}-\frac{(\mu +p)}{c_s^2}\mathcal {F}'\right) +\mathcal {E}\left( \frac{3}{2}\mathcal {F}+1\right) = 0, \nonumber \\{} & {} \quad \implies \frac{\mathcal {E}}{\mathcal {A}^2}= \frac{\frac{(\mu +p)}{c_s^2}\mathcal {F}'-2\mathcal {F}}{\frac{3}{2}\mathcal {F}+1}. \end{aligned}$$
    (43)

Subcase 2 gives us a new constraint equation (43) that needs to be consistently satisfied by the equation of state. Now for this constraint equation to time-evolve consistently, its time derivative must also be zero, and evolving (41) once yields

$$\begin{aligned}{} & {} \Theta \left[ \mathcal {A}^2\left( 2\mathcal {F}\mathcal {G}-(\mu +p)\mathcal {G}'\right) - \mathcal {E}\left( \frac{3}{2}\mathcal {F}+1+\frac{3}{2}(\mu +p)\mathcal {F}'\right) \right] =0, \end{aligned}$$
(44)

where we have set

$$\begin{aligned} \mathcal {G} \equiv \mathcal {G}(\mu ) = 2\mathcal {F}-\frac{(\mu +p)}{c_s^2}\mathcal {F}', \end{aligned}$$
(45)

and \(\mathcal {G}'=\frac{\partial \mathcal {G}(\mu )}{\partial \mu }\). As we shall see in the next section, any further time evolution of the above constraint is not required as the constraints (43) and (44) give us the governing differential equation for the equation of state of the matter that gives rise to this given subcase. We summarize the findings in this section in a diagram depicting the various classes of spherically symmetric shear-free perfect fluid solutions characterized by the acceleration scalar \(\mathcal {A}\) in Fig. 1.

Fig. 1
figure 1

Diagram illustrating the various configurations of the property of acceleration and the spacetime classification

Full size image

6 Dynamic and inhomogeneous fluid

In the event that \(\Theta \ne 0\) and \(\mathcal {A} \ne 0\), setting the square-bracketed term in Eq. (44) to zero yields

$$\begin{aligned} \frac{\mathcal {E}}{\mathcal {A}^2}=\frac{2\mathcal {F}\mathcal {G}-(\mu +p)\mathcal {G}'}{\frac{3}{2}\mathcal {F}+1+\frac{3}{2}(\mu +p)\mathcal {F}'}. \end{aligned}$$
(46)

For consistency, Eqs. (43) and (46) need to be equivalent and this condition produces a master differential equation for the equation of state given by

$$\begin{aligned} \frac{\frac{(\mu +p)}{c_s^2}\mathcal {F}'-2\mathcal {F}}{\frac{3}{2}\mathcal {F}+1}= \frac{2\mathcal {F}\mathcal {G}-(\mu +p)\mathcal {G}'}{\frac{3}{2}\mathcal {F}+1+\frac{3}{2}(\mu +p)\mathcal {F}'}, \end{aligned}$$

which simplifies to

$$\begin{aligned}{} & {} \left( \frac{(\mu +p)}{c_s^2}\mathcal {F}'-2\mathcal {F}\right) \left( \frac{3}{2}\mathcal {F}+\frac{3}{2}(\mu +p)\mathcal {F}' + 1 \right) \nonumber \\{} & {} -\left( \frac{3}{2}\mathcal {F}+1\right) \left( 2\mathcal {F}\mathcal {G}-(\mu +p)\mathcal {G}'\right) = 0. \end{aligned}$$
(47)

Equation (47) imposes a constraint on the equation of state. Finding an equation of state which is a solution to (47) will then necessarily make the two constraints (43) and (44) identically compatible. Since these two constraints are identically satisfied by the equation of state (which is the solution of the above equation), any further time evolution of these constraints will also be identically satisfied and the system of equations is well posed. From Eq. (46), both \(\Theta \) and \(\mathcal {A}\) are nonzero hence \(\hat{\mu }\) and \(\hat{p}\) are also nonzero and so these solutions are necessarily dynamic and inhomogeneous but shear-free.

The resulting differential equation (47) is not autonomous as it stands. However, it can be transformed into that type via a change of variables. Letting \(X=(\mu +p)\) yields the autonomous ordinary differential equation given by

$$\begin{aligned}{} & {} -\frac{3}{2} \left( X'-1\right) \left( XX'' -\left( X'\right) ^{2}+\frac{5}{3}X'-\frac{2}{3}\right) X^{3} X'''' \nonumber \\{} & {} \quad + \frac{3}{2}XX'''\left\{ -\left( X -2 X' +2\right) \left( XX''\right) ^{2} \right. \nonumber \\{} & {} \quad +\left( X'-1\right) \left[ -4 \left( X'\right) ^{2} +\left( 4 X +\frac{16}{3}\right) X' +X^{2}-\frac{11}{3}X-\frac{4}{3}\right] XX'' \nonumber \\{} & {} \quad \left. -4\left( X'-1\right) ^{2} \left( X'-\frac{2}{3}\right) \left( -\frac{1}{2}\left( X'\right) ^{2}+\left( X+\frac{1}{2}\right) X' -\frac{7}{6}X\right) \right\} \nonumber \\{} & {} \quad -\frac{3}{2} \left( -13 X' +X+9\right) \left( XX''\right) ^{3} -6\frac{\left( XX''\right) ^{4}}{X'-1} \nonumber \\{} & {} \quad -3 \left( X'-1\right) \left( 8 \left( X'\right) ^{2}+\left( X -12\right) X' -\frac{3}{2} X+\frac{35}{9}\right) \left( XX''\right) ^{2} \nonumber \\{} & {} \quad +3 \left( X'-1\right) ^{2} \left( 6 \left( X'\right) ^{3}+\left( X-\frac{95}{6}\right) \left( X'\right) ^{2}+\left( -\frac{7}{3}X + \frac{233}{18}\right) X' \right. \nonumber \\{} & {} \quad \left. +\frac{4}{3} X - \frac{29}{9}\right) XX'' -6 \left( X'-1\right) ^{4} \left( X'-\frac{2}{3}\right) \left( X'-\frac{4}{3}\right) \left( X'-\frac{5}{6}\right) =0.\nonumber \\ \end{aligned}$$
(48)

This equation is not satisfied by either linear or power laws of the form \(p=\kappa \mu \) or \(p=\kappa \mu ^\gamma ,(\gamma >2)\) respectively. This infers that it has a complicated solution and so carrying out a numerical analysis on the mathematical software Maple [52] yields the solution given in Fig. 2.

Fig. 2
figure 2

Numerical solution for the autonomous system obtained from the mathematical software Maple

Full size image

Now in order to conduct a dynamical analysis, we investigate the stationary points of the differential equation (48). We first create a system of four differential equations by letting

$$\begin{aligned} X'= & {} Y, \nonumber \\ Y'= & {} Z, \nonumber \\ Z'= & {} W, \nonumber \\ W'= & {} f(X,Y,Z,W). \end{aligned}$$
(49)

Further using (48), \(X'''' \equiv W'\), the function f in (49) can be written explicitly giving

$$\begin{aligned} W'= & {} \left\{ -\frac{2}{3} \left( -\frac{3}{2}X W \left( - \left( X^2 Z -2XYZ +2 \right) ^{2} + \left( Y \left( -4 Y ^{2}+ \left( 4X +\frac{16}{3}\right) Y \right. \right. \right. \right. \right. \nonumber \\{} & {} \quad \left. \left. + X ^{2}-\frac{11}{3}X -\frac{4}{3}\right) -1\right) X Z -4 \left( Y \left( Y -\frac{2}{3} \right) \left( -\frac{1}{2} Y ^{2}\right. \right. \nonumber \\{} & {} \quad \left. +\frac{1}{6}\left( 6X +3 \right) Y -\frac{7}{6}X -1\right) ^{2} +\frac{3}{2} \left( -13XYZ +X^2 Z +9 \right) ^{3} \nonumber \\{} & {} \quad +6\frac{ \left( XZ \right) ^{4}}{Y -1}+3 \left( Y \left( 8 Y ^{2}+ \left( X -12 \right) Y -\frac{3}{2}X +{\frac{35}{9}} \right) XZ -1 \right) ^{2} \nonumber \\{} & {} \quad -3 \left( Y \left( \frac{1}{18} \left( 18 Y ^{2}-42Y +24 \right) X +6 Y^{3} -{\frac{95}{6}} Y ^{2}+\frac{233}{18}Y -\frac{29}{9} \right) \right. \nonumber \\{} & {} \left. \left. \left. \left. -1 \right) ^{2}X Z +6 \left( Y \left( Y -\frac{2}{3}\right) \left( Y -\frac{4}{3} \right) \left( Y -\frac{5}{6}\right) -1\right) ^{4}\right) \right) \right\} \Bigg / \nonumber \\{} & {} \left\{ \left( Y \left( X Z - Y ^{2}+\frac{5}{3}Y -\frac{2}{3} \right) -1 \right) X^3\right\} . \end{aligned}$$
(50)

At the stationary points \(X'=Y'=Z'=0\) which imply that \(Y=Z=W=0\). Substituting \(Y=Z=W=0\) into (50) and simplifying yields

$$\begin{aligned} W'&= \frac{40}{9X^3}. \end{aligned}$$
(51)

Clearly \(W'=0\) only if \(X = (\mu +p) \rightarrow \infty \). This means that the stationary points for the autonomous system exist at infinity indicating the presence of a singularity there.

7 Discussion

In this paper we asked and answered a very important question: what are the physical factors that can make a self-gravitating system inhomogeneous with tidal effects but shear-free? We investigated the behaviour of the self-gravitating system of a shear-free perfect fluid with LRS-II symmetry through a geometrical perspective. As demonstrated in the literature [19,20,21], there exists a deep connection between the properties of the shear, tidal effects, pressure isotropy and dissipative fluxes. The new feature that emerged transparently through our geometrical analysis is: the role of the equation of state of a shear-free perfect fluid in facilitating an inhomogeneous dynamical evolution with tidal deformations. Through this geometrical analysis we enrich the tapestry of research in this avenue further. We also note that here we are considering exact solutions which are unstable with respect to anisotropic perturbations. Although there are previous works involving the connections between shear, anisotropy and tidal effects in spherical symmetry, this paper directly connects them to the equation of state of a perfect fluid.

We have provided a comprehensive answer to the question posed in the title of our paper, in the most generic covariant way, to extract the highly nonlinear differential equation that governs the equation of state. The novel results that emerged from our analysis can be summarized as follows:

  • We can completely characterize all spherically symmetric shear-free perfect fluid spacetimes under a single classification that depends on the covariantly defined fluid acceleration and expansion.

  • One of the most important classes of shear-free spherical symmetry is the inhomogeneous dynamical class. We pinpoint the exact constraint, given by a differential equation (48), governing the equation of state. It is interesting to note that the autonomous ODE (48) obtained places a restriction solely on the equation of state. It does not place a restriction on any of the kinematical variables.

  • Since the governing ODE is highly nonlinear in nature, we easily showed that usual equations of state used in astrophysical settings (for example, a linear equation of state or a finite combination of power law equations of state) will not satisfy this equation. Hence all the solutions belonging to this class are extremely special in terms of the matter content. We also presented a numerical solution for this ODE.