1 Introduction

The final outcome of gravitational collapse has been the focus of attention since Laplace and Michell first conceived of the idea of a black or invisible star. One of the early attempts to determine the result of continued gravitational collapse of a homogeneous dust sphere was carried out by Oppenheimer and Snyder in 1939 [1]. The resulting singularity remains hidden behind the trapping horizon allowing us to conclude that the final fate of collapsing homogeneous dust cloud leads to a Schwarzschild black hole. Although highly simplified the Oppenheimer–Snyder collapse model sparked an interest in seeking more general collapse scenarios [2,3,4,5]. The Cosmic Censorship Conjecture hypothesizes that any reasonable matter distribution undergoing gravitational collapse leads to the formation of a black hole i.e., singularity remains hidden behind the horizon at all times [6]. There have been a number of counterexamples to the Cosmic Censorship Conjecture with the discovery of naked singularities as possible end-states of gravitational collapse [7,8,9]. A natural question that arises from these investigations is how the initial static configuration (the state of the stellar fluid just before the onset of collapse) affects the outcome of gravitational collapse. To this end there have been various approaches in modeling dissipative collapse starting from an initial static configuration [10,11,12]. It has been shown that pressure anisotropy, shear, inclusion of charge, dimensionality of spacetime and equation of state of the initially static core affects the subsequent collapse. In a recent investigation, Naidu and Govender [13] showed that two initially static stellar models with the same masses and radii but different pressure profiles undergoing collapse lead to very different temperature profiles, particularly during the latter stages of their evolution. Finding exact solutions of the Einstein field equations describing bounded matter distributions are important in understanding the subsequent gravitational collapse of these objects.

The Einstein field equations describing localized bodies is a system of highly nonlinear partial differential equations which are difficult in general. In seeking solutions to these equations various novel ideas ranging from an ad-hoc specification of the gravitational potentials, imposing an equation of state, prescribing the behavior of the density, pressure or anisotropy profiles ab initio and specifying the spacetime symmetry have been utilized. It is the very nature of the Einstein field equations which connects the curvature of spacetime to the matter content which allows one to either specify the geometry or the matter distribution to determine the behavior of the other. The first successful attempt at modeling the interior of a spherically symmetric star was carried out by Schwarzschild in 1916, in which he considered a matter distribution with uniform density. The Schwarzschild solution is conformally flat and is characterized by isotropic pressure. Conformal flatness implies vanishing of the Weyl tensor which equates to the vanishing of tidal forces. The study of matter at ultra-high densities of the order of \(10^{15} \mathrm{g~cm}^{-3}\) indicate that the transverse and radial stresses within the stellar fluid may not be equal. Local anisotropy may drastically affect the stability of self–gravitating systems as was shown by Chan et al. [14]. Various scenarios have been proposed to incorporate local anisotropy in stellar models some of which are pion condensation (Hartle et al. [15]), neutrino trapping at high densities [16] and different types of phase transitions [17]. The relaxation of the pressure isotropy condition has led to an explosion of exact solutions of the Einstein field equations describing compact objects. Based on fundamental particle interactions the standard linear equation of state has been extended to include the bag constant. This equation of state has been used extensively to model compact objects with anisotropic pressure profiles as well as a non-vanishing electromagnetic field in the stellar interior. These models are well behaved and were shown to mimic neutron stars, pulsars, and strange star candidates. The quadratic equation of state has also been successfully used to model stellar interiors of compact objects such as Her X-1, RXJ 1856-37, SAX J1808.4-3658(SS1) and SAX J1808.4-3658(SS2). Utilizing curvature coordinates Herrera and Barreto derived an algorithm to generate relativistic polytropes with anisotropic pressures [18]. Motivated by the existence of dark energy, Lobo et al. hypothesized the existence of dark stars with an equation of state of the form \(p = \alpha \rho \) where \(-1< \alpha < -1/3\) [19]. A more exotic form of matter distribution is the so-called Chaplygin gas and the generalized Chaplygin gas which reduce to the linear equation of state in the appropriate limit. Stable dark stars are remnants of gravitational collapse which are formed as a result of the repulsive nature of dark energy. The repulsion is sufficiently strong to halt collapse leading to the formation of stable stars free of any singularity [20, 21].

Higher order gravity theories have been fruitful in producing models of compact stellar objects. Various authors have shown that modifications to 4-D classical Einstein gravity feature in the thermodynamical properties of the stellar fluid [22,23,24]. The braneworld scenario provides a natural mechanism for the existence of anisotropic pressure within the stellar fluid [26]. In addition, it was shown that in the Randall–Sundrum II type braneworld, the exterior spacetime of spherical star is filled with radiative-type stresses induced by five-dimensional graviton effects and is not necessarily the vacuum Schwarzschild solution as in the 4-D case [27]. Recently, Dadhich et al. [28] have shown that within the framework of pure Lovelock gravity there cannot exist self-gravitating bounded distributions \(d =2N + 1\) dimensions. This is to say that there is no finite radius for which the pressure vanishes [28]. The transition from classical 4-D gravity to higher-dimensional gravity theories has sparked immense interest in studying phenomenological processes which reside in extra dimensions. One of the main proponents of these investigations is Dadhich and his collaborators who proved the universality of the Schwarzschild constant density sphere, i.e., it was shown that this solution carries over to Einstein–Gauss–Bonnet gravity and Lovelock gravity [25, 29].

It is widely believed that the four fundamental interactions in Nature were once a manifestation of a single, unified force. Furthermore, the dimensionality of spacetime could have evolved in such a manner so as to reveal four dimensions which we observe today. Kaluzua–Klein theories have shown that the electromagnetic interaction manifests naturally in five-dimensional spacetime. These observations of physical phenomena transcending the dimensionality of spacetime have generated widespread interest in embedding our standard four-dimensional spacetime into higher-dimensional spacetimes [30]. It is well known that any pseudo-Riemannian manifold, \((\mathcal{V}_n)^{-}\) with dimensionality n may be locally embedded into a pseudo-Euclidean space, \((\mathcal{V}_m)^{+}\) of dimension \(m = \frac{n(n + 1)}{2}\). It follows that the embedding class of \((\mathcal{V}_n)^{-} \le m - n = \frac{n(n - 1)}{2}\). For the relativistic 4-D spacetime \((\mathcal{V}_4)^{-}\), the embedding class is 6. A recent and popular approach in deriving exact solutions of the Einstein field equations describing compact stars is to make use of the Karmarkar condition [31,32,33,34,35,36,37,38,39]. The necessary and sufficient condition for a spherically symmetric spacetime to be of embedding class I was first derived by Karmarkar [41]. It is a mathematical simplification which reduces the problem of obtaining exact solutions to a single-generating function. The approach is to choose one of the gravitational potentials on physical grounds and to then integrate the Karmarkar condition to fully specify the gravitational behavior of the model. In this paper we utilize the Karmarkar condition to derive solutions which describe compact objects in general relativity. We subject our solutions to rigorous physical tests which ensure that they do describe physically observable objects in the universe.

2 Einstein field equations for anisotropic fluid distributions

The interior of the super-dense star is assumed to be described by the line element

$$\begin{aligned} \mathrm{d}s^2 = -e^{\nu (r)}\mathrm{d}t^2 + e^{\lambda (r)}\mathrm{d}r^2 + r^2(\mathrm{d}\theta ^2 + \sin ^2{\theta }\mathrm{d}\phi ^2) \end{aligned}$$
(1)

where the gravitational potentials \(\nu (r)\) and \(\lambda (r)\) are yet to be specified. The Einstein field equations describing an anisotropic fluid distribution are given as (in the unit \(G=c=1\))

$$\begin{aligned} -8\pi T^\mu _\xi = {\mathscr {R}}^\mu _\xi -{1\over 2}{\mathscr {R}}~g^\mu _\xi \end{aligned}$$
(2)

where

$$\begin{aligned} T^\mu _\xi= & {} \rho v^\mu v_\xi + p_r \chi _\xi \chi ^\mu + p_t(v^\mu v_\xi -\chi _\xi \chi ^\mu -g^\mu _\xi ) ~ \end{aligned}$$
(3)

is the energy-momentum tensor, \({\mathscr {R}}^\mu _\xi \) is the Ricci tensor, \({\mathscr {R}}\) represents the scalar curvature, \(p_r\) and \(p_t\) denote radial and transverse pressures, respectively, \(\rho \) is the density of the fluid distribution, \(v^\mu \) the four velocity, and \(\chi ^\mu \) is the unit space-like vector in the radial direction.

The Einstein field equations (2) for the line element (1) are

$$\begin{aligned} 8\pi \rho (r)= & {} \frac{1 - e^{-\lambda }}{r^2} + \frac{\lambda 'e^{-\lambda }}{r} \end{aligned}$$
(4)
$$\begin{aligned} 8\pi p_r(r)= & {} \frac{\nu ' e^{-\lambda }}{r} - \frac{1 - e^{-\lambda }}{r^2} , \end{aligned}$$
(5)
$$\begin{aligned} 8\pi p_t(r)= & {} \frac{e^{-\lambda }}{4}\left( 2\nu '' + {\nu '}^2 - \nu '\lambda ' + \frac{2\nu '}{r}-\frac{2\lambda '}{r}\right) , \end{aligned}$$
(6)

where primes denote differentiation with respect to the radial coordinate r. In generating the above field equations we have utilized geometrized units where the coupling constant and the speed of light are taken to be unity. Using Eqs. (5) and (6) we obtain the anisotropic parameter

$$\begin{aligned} \Delta (r)= & {} 8\pi (p_t-p_r) \nonumber \\= & {} e^{-\lambda }\left[ {\nu '' \over 2}-{\lambda ' \nu ' \over 4}+{\nu '^2 \over 4}-{\nu '+\lambda ' \over 2r}+{e^\lambda -1 \over r^2}\right] , \end{aligned}$$
(7)

which vanishes in the case of an isotropic pressure.

Eisenhart [40] has mentioned that, for any Riemannian space to be class I, a necessary and sufficient condition is that there exists a second-order symmetric tensor \(b_{\mu \alpha }\) satisfying the following equations:

$$\begin{aligned} {\mathscr {R}}_{\mu \nu \alpha \beta }= & {} \varepsilon (b_{\mu \alpha }b_{\nu \beta }-b_{\mu \beta }b_{\nu \alpha }) \end{aligned}$$
(8)
$$\begin{aligned} 0= & {} b_{\mu \nu ;\alpha }-b_{\mu \alpha ;\nu } \end{aligned}$$
(9)

where \(\varepsilon = \pm 1\) (\(+\) when the normal to the manifold is space-like or − when the normal to the manifold is time-like) and ‘(;)’ represents covariant differentiation.

For the line element (1), the non-zero components of the Riemann curvature tensor are given by

$$\begin{aligned} {\mathscr {R}}_{1414}= & {} -e^\nu \left( {\nu '' \over 2}+{\nu '^2 \over 4}-{\lambda ' \nu ' \over 4}\right) \end{aligned}$$
(10)
$$\begin{aligned} {\mathscr {R}}_{2323}= & {} -e^\lambda r^2 \sin ^2\theta ~ (e^\lambda -1) \end{aligned}$$
(11)
$$\begin{aligned} {\mathscr {R}}_{1334}= & {} {\mathscr {R}}_{1224}~\sin ^2\theta = 0 \end{aligned}$$
(12)
$$\begin{aligned} {\mathscr {R}}_{1212}= & {} {1\over 2}r \lambda '\end{aligned}$$
(13)
$$\begin{aligned} {\mathscr {R}}_{3434}= & {} -{1\over 2}r\sin ^2\theta ~ \nu ' e^{\nu -\lambda } \end{aligned}$$
(14)

The non-zero components of the tensor \(b_{\mu \alpha }\) corresponding to (1) are \(b_{11},~b_{22},~b_{33}, ~b_{44}\) and \(b_{14}\) with \(b_{33}=b_{22}\sin ^2 \theta \). With these components, (8) reduces to

$$\begin{aligned} {\mathscr {R}}_{1414}={{\mathscr {R}}_{1212}{\mathscr {R}}_{3434}+{\mathscr {R}}_{1224}{\mathscr {R}}_{1334} \over {\mathscr {R}}_{2323}}, \end{aligned}$$
(15)

which is known as the Karmarkar condition [41] in the literature.

Using (10)–(14) in (15) leads to the following differential equation:

$$\begin{aligned} {\lambda ' e^\lambda \over e^\lambda -1} = {2\nu '' \over \nu '}+\nu ', \end{aligned}$$
(16)

which can easily be integrated to give a relationship between \(\nu (r)\) and \(\lambda (r)\):

$$\begin{aligned} e^\lambda = 1 + \frac{K{\nu '}^2e^\nu }{4} \end{aligned}$$
(17)

where K is a constant of integration.

By using (17) we can rewrite (7) as

$$\begin{aligned} \Delta (r) = {\nu ' \over 4e^\lambda }\left[ {2\over r}-{\lambda ' \over e^\lambda -1}\right] ~\left[ {\nu ' e^\nu \over 2rB^2}-1\right] ,~~~~B={1\over \sqrt{K}}. \end{aligned}$$
(18)

However, Pandey and Sharma [42] argued that satisfying the Karmarkar condition alone is insufficient for a spherically symmetric spacetime to be class I. As an example, they presented the following spacetime:

$$\begin{aligned} \mathrm{d}s^2=-e^\nu \mathrm{d}t^2+\mathrm{d}r^2+r^2(\mathrm{d}\theta ^2+\sin ^2 \theta \mathrm{d}\phi ^2), \end{aligned}$$
(19)

which does satisfy (15). This spacetime (19) has \(e^\lambda =1\) implying \({\mathscr {R}}_{2323}= 0\) from (11). \(e^\lambda =1\) or \({\mathscr {R}}_{2323}= 0\) also implies (19) is spatially flat.

Now the non-zero components of curvature tensor for (19) are \({\mathscr {R}}_{1414},~{\mathscr {R}}_{2424}\) and \({\mathscr {R}}_{3434}\) only. Using these components, (8) implies inconsistent equations:

$$\begin{aligned} b_{22}b_{33}=0,~b_{24}b_{33}=0,~b_{22}b_{44}-b_{24}^2 \ne 0,~b_{33}b_{44}\ne 0.\nonumber \\ \end{aligned}$$
(20)

Therefore, the spacetime given in (19) does satisfy the Karmarkar condition but fails to satisfy (8) i.e. (19) is not a class I spacetime due to \({\mathscr {R}}_{2323}= 0\). Hence, any symmetric spacetimes are called class I if they satisfy the Karmarkar condition and the Pandey–Sharma condition (\({\mathscr {R}}_{2323}\ne 0\)) simultaneously. The condition \({\mathscr {R}}_{2323}= 0\) or equivalently \(e^\lambda = 1\) gives the spacetime (19), which in fact describes a perfect fluid sphere with zero density [42]. Hence, in order to describe a perfect fluid with non-vanishing density we require \({\mathscr {R}}_{2323} \ne 0\). It is also important to note that all the spherically symmetric spacetimes are in general of class II unless they simultaneously satisfy the Karmarkar and Pandey–Sharma conditions.

3 Isotropic Class I solutions

For isotropy in pressure, the anisotropy factor \(\Delta = 0\). Assuming that \(\nu '(r) \ne 0\), we will get from (18) either

$$\begin{aligned} \left[ {2\over r}-{\lambda ' \over e^\lambda -1}\right]= & {} 0~~~~\text{ or } \end{aligned}$$
(21)
$$\begin{aligned} \left[ {\nu ' e^\nu \over 2rB^2}-1\right]= & {} 0 \end{aligned}$$
(22)

or both. The first condition (21) leads to Schwarzschild’s constant density model [43] and the second condition (22) leads to Kohler–Chao solution [44].

3.1 Schwarzschild interior solution

Integration of (21) yields

$$\begin{aligned} e^{-\lambda }=1-cr^2. \end{aligned}$$
(23)

Using (23) in (17), we obtain

$$\begin{aligned} e^\nu =\Big (A-{B \over \sqrt{c}} \sqrt{1-cr^2}\Big )^2. \end{aligned}$$
(24)

The above solution is the well-known interior Schwarzschild model, which describes an incompressible, static sphere with uniform density. For completeness we present the physical quantities of this solution as determined from (4) and (5),

$$\begin{aligned} \rho (r)= & {} {3c \over 8\pi },\end{aligned}$$
(25)
$$\begin{aligned} P(r)= & {} {c \over 8\pi } \left( \frac{2 B \sqrt{1-c r^2}}{A \sqrt{c}-B \sqrt{1-c r^2}}-1\right) ,\end{aligned}$$
(26)
$$\begin{aligned} {P(r) \over \rho (r)}= & {} \frac{1}{3} \left( \frac{2 B \sqrt{1-c r^2}}{A \sqrt{c}-B \sqrt{1-c r^2}}-1\right) . \end{aligned}$$
(27)

The interior Schwarzschild solution has been extensively studied by various authors including Schwarzschild himself [43]. This solution serves as a toy model for self-gravitating bounded configurations. One of its main shortcomings is the fact that it leads to an infinite speed of sound within the interior of the sphere.

3.2 Kohler–Chao solution: a cosmological solution

Integrating (22) we obtain

$$\begin{aligned} e^{\nu }=A+B r^2 \end{aligned}$$
(28)

and using (28) in (17) yields

$$\begin{aligned} e^{\lambda }={A+2Br^2 \over A+Br^2}, \end{aligned}$$
(29)

which is the Kohler–Chao–Tikekar solution [44, 46].

The corresponding expressions for density, pressure, and equation of state parameter can be written as

$$\begin{aligned} 8\pi \rho (r)= & {} {B(3A+2Br^2) \over (A+2Br^2)^2},\end{aligned}$$
(30)
$$\begin{aligned} 8\pi P(r)= & {} {B \over A+2Br^2} ,\end{aligned}$$
(31)
$$\begin{aligned} {P(r) \over \rho (r)}= & {} \frac{A+2 B r^2}{3 A+2 B r^2}, \end{aligned}$$
(32)

with \(B > 0\). However, we can see clearly from (31) that the pressure at the surface of any configuration cannot be zero for a finite boundary unless the boundary itself is infinite. This property of an infinite boundary does have the property of a cosmological solution. The same discussion is also given in Maurya et al. [45].

Maurya et al. [45] comment on the charged isotropic solutions of embedding class I, i.e., the Schwarzschild interior and Kohler–Chao solutions. If the charge vanishes in these two solutions, then the remaining neutral counterpart will only be either the Schwarzschild interior solution or the Kohler–Chao solution, otherwise either the charge cannot be zero or the surviving spacetime metric will become flat.

It is well known that an isotropic spherically symmetric conformally flat metric is necessarily to be a class I solution. However, does the converse hold good, i.e., is a spherically symmetric class I solution representing an isotropic fluid sphere necessarily conformally flat? This was resolved by Tikekar [46] concluding that “it is not necessary a spherically symmetric class I solution representing isotropic fluid sphere to be conformally flat” and he gave an example by rediscovering the Kohler–Chao solution. Here what we want to stress is that the conformally flat solution, i.e., the Schwarzschild interior solution is the only isotropic class I solution that can represent a bounded stellar configuration. However, the conformally non-flat solution, i.e., the Kohler–Chao solution, cannot describe a finite bounded configuration, although it can qualify as a cosmological solution.

It is well known that all the non-vanishing components of the Weyl tensor are proportional to

$$\begin{aligned} W= & {} \frac{r^3e^{-\lambda }}{6}\left[ \frac{e^\lambda }{r^2} - \frac{1}{r^2} + \frac{\nu ' \lambda '}{4} - \frac{\nu '^2}{4} - \frac{\nu ''}{2} +\frac{\nu '-\lambda '}{2r}\right] .\nonumber \\ \end{aligned}$$
(33)

Here the Schwarzschild interior solution yields a vanishing Weyl tensor (\(W=0\)) postulating that it is a conformally flat space. However, the Kohler–Chao solution yields a non-vanishing Weyl tensor where

$$\begin{aligned} W = \frac{2B^2r^5}{3(A + 2Br^2)^2}, \end{aligned}$$
(34)

implying that the Kohler–Chao solution is not conformally flat. In general, the Karmarkar condition and the pressure isotropy do not imply conformal flatness. However, the converse is true: conformally flat, perfect fluid spheres obey the Karmarkar condition.

4 Generating a new family of embedding class I models

We now seek relativistic stellar models which satisfy the Karmarkar condition. In the light of our findings in the previous section we relax the condition of pressure isotropy. This implies that the radial and tangential stresses are unequal throughout the fluid distribution. It is well known that pressure anisotropy plays an important role during dissipative collapse. In a recent study by Govender et al. [47] it has been shown that the dynamics of a collapsing core is closely related to the radial pressure and energy density of the stellar fluid. By assuming a linear equation of state for the initial static configuration, of the form \(p_r = \alpha \rho - \beta \), where \(\alpha \) and \(\beta \) are constants, they demonstrated that the subsequent collapse is sensitive to the interplay between the radial pressure and energy density. They also demonstrated that the equation of state parameter, \(\alpha \), influences the behavior of the temperature profile of the collapsing body.

We now proceed to obtain a family of solutions which describe anisotropic matter configurations obeying the Karmarkar condition. In order to completely specify the gravitational behavior of our model we assume

$$\begin{aligned} e^\lambda = a r^2 \sin ^2\left( b r^2+c\right) +1 \end{aligned}$$
(35)

where a, b, and c are constants which are determined from the boundary conditions. The sinusoidal behavior of the gravitational potential has been widely used in various contexts in both cosmology and astrophysics. Dadhich and Raychaudhuri demonstrated that it was possible to obtain an oscillating cosmological model without Big Bang singularity [48]. An interesting feature of this model is that it allows for the prediction of blue-shifts without violating the basic postulates of general relativity. In modeling a dissipative gravitational collapse of a spherically symmetric star in which the Weyl stresses vanish, Maharaj and Govender [49] showed that the solution of the boundary condition admits oscillatory solutions. The extension from 4-D to 5-D gravity of the Finch and Skea stellar model leads to sinusoidal behavior of the gravitational potentials [22].

Using the metric potential (35) in (17), we get

$$\begin{aligned} e^\nu = \left[ A-\frac{\sqrt{a} B }{2 b}~\cos \left( b r^2+c\right) \right] ^2. \end{aligned}$$
(36)

Using (35) and (36), we can rewrite the expression of density, \(p_r\), \(\Delta \), and \(p_t\) as

$$\begin{aligned} 8\pi \rho (r)= & {} {a \over \left[ a r^2 \sin ^2\left( b r^2+c\right) +1\right] ^2} ~\bigg [a r^2 \sin ^4\left( b r^2+c\right) \nonumber \\&+3 \sin ^2\left( b r^2+c\right) +2 b r^2 \sin \left\{ 2 \left( b r^2+c\right) \right\} \bigg ] ,\nonumber \\ \end{aligned}$$
(37)
$$\begin{aligned} 8\pi p_r(r)= & {} \frac{\sqrt{a} ~\sin \left( b r^2+c\right) }{2 \left[ a r^2 \sin ^2\left( b r^2+c\right) +1\right] } \end{aligned}$$
(38)
$$\begin{aligned}& \times \frac{4 \sqrt{a} A b \sin \left( b r^2+c\right) -a B \sin \left\{ 2 \left( b r^2+c\right) \right\} -8 b B}{\sqrt{a} B \cos \left( b r^2+c\right) -2 A b} ,\nonumber \\\end{aligned}$$
(39)
$$\begin{aligned} \Delta (r)= & {} \frac{r ~\csc ^4\left( b r^2+c\right) }{4 \left[ a r^2+\csc ^2\left( b r^2+c\right) \right] ^2} \nonumber \\&\times \frac{a \cos \left\{ 2 \left( b r^2+c\right) \right\} -a+4 b \cot \left( b r^2+c\right) }{2 A b-\sqrt{a} B \cos \left( b r^2+c\right) } \nonumber \\&\times \bigg [2 a A b r \cos \left\{ 2 \left( b r^2+c\right) \right\} \nonumber \\&-\,2 a A b r +4 \sqrt{a} b B r \nonumber \\&\times \sin \left( b r^2+c\right) + a^{3/2} B r \sin \left( b r^2+c\right) \nonumber \\&\times \sin \left\{ 2 \left( b r^2+c\right) \right\} \bigg ],\end{aligned}$$
(40)
$$\begin{aligned} 8\pi p_t(r)= & {} 8\pi p_r(r)+\Delta (r). \end{aligned}$$
(41)

5 Properties of the new model

The central values of pressure and density are given by

$$\begin{aligned} 8\pi p_{rc}= & {} 8\pi p_{tc} \nonumber \\= & {} \frac{\sqrt{a}~ \sin c \big (4 \sqrt{a} A b \sin c-a B \sin (2 c)-8 b B\big )}{2 \big (\sqrt{a} B \cos c-2 A b\big )},\nonumber \\\end{aligned}$$
(42)
$$\begin{aligned} 8\pi \rho _{c}= & {} 3 a \sin ^2c. \end{aligned}$$
(43)

To satisfy Zeldovich’s condition at the interior, \(p_{rc}/\rho _c\) at the center must be \(\le 1\). Therefore,

$$\begin{aligned} \frac{4 \sqrt{a} A b \sin c-a B \sin (2 c)-8 b B}{3 \sqrt{a} \sin c\big [2 \sqrt{a} B \cos c-4 A b\big ]} \le 1. \end{aligned}$$
(44)

On using (42) and (44) we generate a constraint on B / A given by

$$\begin{aligned} {8b+a \sin (2c) \over 4\sqrt{a} Ab \sin c} <{A \over B} \le {a \sin (2c) +2b \over 4\sqrt{a} b \sin c}. \end{aligned}$$
(45)

6 Matching of physical boundary conditions

The exterior spacetime of our static model is the vacuum Schwarzschild solution given by

$$\begin{aligned} \mathrm{d}s^2= & {} \left( 1-{2M\over r}\right) \mathrm{d}t^2-\left( 1-{2M\over r}\right) ^{-1}\mathrm{d}r^2 \nonumber \\&-r^2(\mathrm{d}\theta ^2+\sin ^2 \theta \mathrm{d}\phi ^2). \end{aligned}$$
(46)

By matching the first and second fundamental forms the interior solution (1) and exterior solution (46) at the boundary \(r=R\) (Darmois–Israel junction conditions) we obtain

$$\begin{aligned}&e^{\nu _b} = 1-{2M \over R} = \left[ A-\frac{\sqrt{a} B }{2 b}~\cos \left( b R^2+c\right) \right] ^2,\end{aligned}$$
(47)
$$\begin{aligned}&e^{-\lambda _b} = 1-{2M \over R} = \Big [a R^2 \sin ^2\left( b R^2+c\right) +1\Big ]^{-1} ,\end{aligned}$$
(48)
$$\begin{aligned}&p_r(R) = 0 . \end{aligned}$$
(49)

Using the boundary condition (47)–(49), we get

$$\begin{aligned} B= & {} \frac{4 \sqrt{a} A b \sin \left( b R^2+c\right) }{a \sin \left\{ 2 \left( b R^2+c\right) \right\} +8 b} ,\end{aligned}$$
(50)
$$\begin{aligned} A= & {} \frac{\sqrt{1-2 M/R} \Big [a \sin \left\{ 2 \left( b R^2+c\right) \right\} +8 b\Big ]}{a \sin \left\{ 2 \left( b R^2+c\right) \right\} -a \sin \left( 2 b R^2+2 c\right) +8 b},\end{aligned}$$
(51)
$$\begin{aligned} a= & {} \bigg [\frac{1}{1-2 M/R}-1\bigg ]~{1 \over R^2 \sin ^2\left( b R^2+c\right) } , \end{aligned}$$
(52)

and we have chosen b,  c,  M and R as free parameters and the rest of the constants a,  A and B are determined from Eqs. (50)–(52).

The gravitational red-shift of the stellar system is given by

$$\begin{aligned} Z(r)= & {} \left[ A-\frac{\sqrt{a} B }{2 b} ~\cos \left( b r^2+c\right) \right] ^{-1}-1. \end{aligned}$$
(53)

The mass–radius relation and the compactness parameter of the solution can be determined using the equation given by

$$\begin{aligned} m(r)= & {} 4\pi \int _0^r \rho r^2 \mathrm{d}r = \frac{a r^3 \sin ^2\left( b r^2+c\right) }{a r^2+2-a r^2 \cos \left( 2 b r^2+2 c\right) },\nonumber \\\end{aligned}$$
(54)
$$\begin{aligned} u(r)= & {} {2m(r) \over r} = \frac{2a r^2 \sin ^2\left( b r^2+c\right) }{a r^2+2-a r^2 \cos \left( 2 b r^2+2 c\right) }. \end{aligned}$$
(55)

7 Equilibrium and stability conditions

7.1 Condition for equilibrium

For a stellar system in equilibrium under different forces, the generalized Tolman–Oppenheimer–Volkoff (TOV) equation must be satisfied [50], i.e.,

$$\begin{aligned} {2\Delta \over r} = {\mathrm{d}p_r \over \mathrm{d}r}+{M_g(\rho +p_r) \over r^2}~e^{(\lambda -\nu )/2} \end{aligned}$$
(56)

where \(M_g(r)\) is the effective gravitational mass contained within a sphere of radius r and is defined by the Tolman–Whittaker formula viz.,

$$\begin{aligned} M_g(r)= & {} 4 \pi \int _0^r \big (T^t_t-T^r_r-T^\theta _\theta -T^\phi _\phi \big ) r^2 ~e^{(\nu +\lambda )/2}\mathrm{d}r. \end{aligned}$$
(57)

For Eqs. (4)–(6), Eq. (57) reduces to

$$\begin{aligned} M_g(r)= & {} {1\over 2}~r^2 \nu ' ~e^{(\nu -\lambda )/2}. \end{aligned}$$
(58)

Equation (56) can be written in terms of balanced force equation due to the anisotropy (\(F_a\)), gravity (\(F_g\)) and hydrostatic force (\(F_h\)), i.e.,

$$\begin{aligned} F_g+F_h+F_a=0. \end{aligned}$$
(59)

Here

$$\begin{aligned} F_g= & {} -{M_g(\rho +p_r) \over r^2}~e^{(\lambda -\nu )/2},\end{aligned}$$
(60)
$$\begin{aligned} F_h= & {} -{\mathrm{d}p_r \over \mathrm{d}r},\end{aligned}$$
(61)
$$\begin{aligned} F_a= & {} {2\Delta \over r}. \end{aligned}$$
(62)

The TOV equation (59) can be represented graphically showing the interplay amongst \(F_g\), \(F_h\), and \(F_a\) required to bring about equilibrium as evidenced in Fig. 1.

7.2 Relativistic adiabatic index and stability

For a relativistic anisotropic sphere the stability is related to the adiabatic index \(\Gamma \), the ratio of two specific heats, defined by [14]

$$\begin{aligned} \Gamma =\frac{\rho +p_r}{p_r}\frac{\mathrm{d}p_r}{\mathrm{d}\rho }. \end{aligned}$$
(63)

Now \(\Gamma >4/3\) gives the condition for the stability of a Newtonian sphere and \(\Gamma =4/3\) is the condition for a neutral equilibrium proposed by [51]. This condition changes for a relativistic isotropic sphere due to the regenerative effect of pressure, which renders the sphere more unstable. For an anisotropic general relativistic sphere the situation becomes more complicated, because the stability will depend on the type of anisotropy. For an anisotropic relativistic sphere the stability condition is given by [14],

$$\begin{aligned} \Gamma >\frac{4}{3}+\left[ \frac{4}{3}\frac{(p_{t0}-p_{r0})}{|p_{r0}^\prime |r}+\frac{1}{2}\kappa \frac{\rho _0p_{r0}}{|p_{r0}^\prime |}r\right] , \end{aligned}$$
(64)

where \(p_{r0}\), \(p_{t0}\), and \(\rho _0\) are the initial radial, tangential, and energy density in static equilibrium satisfying (56). The first and last terms inside the square brackets represent the anisotropic and relativistic corrections, respectively, and both quantities are positive, which increases the unstable range of \(\Gamma \) [14, 52].

Fig. 1
figure 1

Balancing of different forces in TOV equation for static configurations of 4U1538-52, LMC X-4, and PSR J1614-2230 are plotted with radial coordinate r

Full size image
Fig. 2
figure 2

Variation of \(v_r^2\) and \(v_t^2\) with radial coordinate r for 4U1538-52, LMC X-4, and PSR J1614-2230 with their respective parameters given in Table 1

Full size image
Fig. 3
figure 3

Variation of \(\mathrm{d}\rho /\mathrm{d}r,~\mathrm{d}p_r/\mathrm{d}r\) and \(\mathrm{d}p_t/\mathrm{d}r\) (km\(^{-1}\)) with radial coordinate r for 4U1538-52, LMC X-4 and PSR J1614-2230 with their respective parameters given in Table 1

Full size image
Fig. 4
figure 4

Variation of anisotropy (km\(^{-2}\)) with radial coordinate r for 4U1538-52, LMC X-4, and PSR J1614-2230 with their respective parameters given in Table 1

Full size image
Fig. 5
figure 5

Variation of stability factor \(v_t^2-v_r^2\) with radial coordinate r for 4U1538-52, LMC X-4, and PSR J1614-2230 with their respective parameters given in Table 1

Full size image
Fig. 6
figure 6

Variation of mass with central density \(8\pi \rho _c~(0-2.68 \times 10^{17}~\mathrm{g}/\mathrm{cm}^3)\) for \(R=7{-}10~\mathrm{km}\)

Full size image
Fig. 7
figure 7

Variation of metric potentials with radial coordinate r for 4U1538-52, LMC X-4, and PSR J1614-2230 with their respective parameters given in Table 1

Full size image
Fig. 8
figure 8

Variation of interior pressures (km\(^{-2}\)) with radial coordinate r for 4U1538-52, LMC X-4, and PSR J1614-2230 with their respective parameters given in Table 1

Full size image
Fig. 9
figure 9

Variation of density (km\(^{-2}\)) with radial coordinate r for 4U1538-52, LMC X-4, and PSR J1614-2230 with their respective parameters given in Table 1

Full size image

7.3 Causality and stability condition

The radial and tangential speeds of sound of our compact star model are given by

$$\begin{aligned} v_r^2=\frac{\mathrm{d}p_r}{\mathrm{d}\rho }={\mathrm{d}p_r/\mathrm{d}r \over \mathrm{d}\rho /\mathrm{d}r},~~~v_t ^2=\frac{\mathrm{d}p_t}{\mathrm{d}\rho }={\mathrm{d}p_t/\mathrm{d}r \over \mathrm{d}\rho /\mathrm{d}r}. \end{aligned}$$
(65)

The profiles of \(v_r^2\) and \(v_t^2\) are given in Fig. 2, which indicates that both the radial and transverse velocity satisfy the causality conditions, i.e., both \(v_r^2,\,v_t^2\) are less than 1 and are monotonic decreasing functions of r.

The stability of anisotropic stars under the radial perturbations is studied by using the concept of [53] known as Hererra’s “cracking” method. Using the concept of cracking, [54] showed that the region of the anisotropic fluid sphere where \(-1\le v_{t}^{2}-v_{r}^{2}\le 0\) is potentially stable but the region where \(0<v_{t}^{2}-v_{r}^{2}\le 1\) is potentially unstable. We have

$$\begin{aligned} \frac{\mathrm{d}p_t}{\mathrm{d}\rho }= & {} \frac{\mathrm{d}p_r}{\mathrm{d}\rho }+\frac{\mathrm{d}\Delta }{\mathrm{d}\rho } = \frac{\mathrm{d}p_r}{\mathrm{d}\rho }+\frac{\mathrm{d}\Delta /\mathrm{d}r}{\mathrm{d}\rho /\mathrm{d}r} \nonumber \\ i.e.,~~~ v_t^2-v_r^2= & {} \frac{\mathrm{d}p_r}{\mathrm{d}\rho }+\frac{\mathrm{d}\Delta /\mathrm{d}r}{\mathrm{d}\rho /\mathrm{d}r}. \end{aligned}$$
(66)

In order to maintain \(-1\le v_{t}^{2}-v_{r}^{2}\le 0\) throughout the fluid distribution it is required that \(\mathrm{d}\Delta /\mathrm{d}r>0\) (from (66)) as we have \(\mathrm{d}\rho /\mathrm{d}r<0\) (see Fig. 3), i.e., it is required that \(\Delta \) is an increasing function of r, which is already satisfied by our model (see Fig. 4). With the help of a graphical representation we have also shown that \(v_t^2-v_r^2<0\) in Fig. 5 everywhere inside the fluid sphere, which renders our model stable.

Fig. 10
figure 10

Variation of pressure to density ratios with radial coordinate r for 4U1538-52, LMC X-4, and PSR J1614-2230 with their respective parameters given in Table 1

Full size image
Fig. 11
figure 11

Variation of red-shift with radial coordinate r for 4U1538-52, LMC X-4, and PSR J1614-2230 with their respective parameters given in Table 1

Full size image
Fig. 12
figure 12

Variation of relativistic adiabatic index with radial coordinate r for 4U1538-52, LMC X-4, and PSR J1614-2230 with their respective parameters given in Table 1

Full size image
Fig. 13
figure 13

Variation of \(\rho -p_r,~\rho -p_t\) and \(\rho -p-r-2p_t\) (km\(^{-2}\)) with radial coordinate r for 4U1538-52, LMC X-4, and PSR J1614-2230 with their respective parameters given in Table 1

Full size image
Fig. 14
figure 14

Variation of interior mass with radial coordinate r for 4U1538-52, LMC X-4, and PSR J1614-2230 with their respective parameters given in Table 1

Full size image
Fig. 15
figure 15

Variation of compactness parameter with radial coordinate r for 4U1538-52, LMC X-4, and PSR J1614-2230 with their respective parameters given in Table 1

Full size image

7.4 Harrison–Zeldovich–Novikov static stability criterion

The stability analysis adopted by [55, 56], among other treatments, requires the determination of eigen-frequencies of all the fundamental modes. However, [56] and [57] provide a simpler formalism to study the stability of the stellar model. They have assumed that the adiabatic index of a pulsating star is the same as in slowly deformed matter. This leads to a stable configuration only if the mass of the star is increasing with central density, i.e., \(\mathrm{d}M/\mathrm{d}\rho _c > 0\) and unstable if \(\mathrm{d}M/\mathrm{d}\rho _c \le 0\).

In our model, the mass as a function of the central density can be written as

$$\begin{aligned} M = \frac{8\pi \rho _c R^3 \sin ^2(c+bR^2)/3\sin ^2c}{8\pi \rho _c R^2 /3\sin ^2c -8\pi \rho _c R^2 \cos (2c+2bR^2)/3\sin ^2c+2}, \end{aligned}$$
Table 1 Parameters of three well-known compact stars that give masses and radii compatible with observational data.
Full size table

which gives us (for a given radius)

$$\begin{aligned} {\mathrm{d} M \over \mathrm{d} \rho _c} = \frac{12 \pi R^3 \sin ^2c ~~\sin ^2(b R^2+c)}{\left[ 8 \pi R^2 \rho _c \sin ^2c ~~\sin ^2(b R^2+c)+3\right] ^2} > 0. \end{aligned}$$
(67)

Figure 6 shows that our models are stable according to the static stability criterion. It is interesting to note that the stability of our configurations is enhanced with increasing radii and plateaus after attaining a maximum value for the respective central matter densities. Wherever the curve starts leveling off, it implies that \(\mathrm{d}M/\mathrm{d}\rho _c = 0\), indicating that the configuration is rendered unstable.

8 Discussion of results

Graphical analyses of the physical parameters \(\big (e^{-\lambda },~p_r,\) \(p_t,~\rho ,\) \(~p_r/\rho ,~p_t/\rho ,~v_r^2,~v_t^2,~Z\big )\) show that they are finite at the center and monotonically decreasing outward (Figs. 2, 7, 8, 9, 10, and 11). Figures 4, 7, and 12 show that \(e^\nu \), the anisotropy parameter, \(\Delta \), and \(\Gamma \) are increasing radially outward.

The null energy condition \(\big (\rho -p_i \ge 0\big )\), dominant energy condition \(\big (\rho -p_i \ge 0\), \(\rho \ge 0\big )\) and strong energy condition \(\big (\rho -p_i \ge 0\), \(\rho -p_r-2p_t \ge 0\big )\) are simultaneously satisfied by our solution (Fig. 13). The solution can also represent static and stable stellar configurations as the stability factor \(v_t^2-v_r^2\) lies between the limits \(-1\) and 0 (Fig. 5). For a non-collapsing stellar configuration, the adiabatic index must also be greater than 4/3 for positive values of anisotropy, which can be seen from (Fig. 12). Furthermore, the gravitational force \(F_g\) in the configuration is balanced by the combined effect of hydrostatic force \(F_h\) and anisotropic force \(F_a\) (Fig. 1) and thus the solution satisfies the TOV equation, Eq. (56). The mass and the compactness parameter also monotonically increase from the center to the surface of the star and the compactness parameter is also within the Buchdahl limit, i.e., \(u \le 8/9\) (Figs. 14 and 15). The negative values of the gradients of density and pressures signify that the density and pressures are decreasing radially outward (Fig. 3).

The well-behaved nature of the solution depends on the parameter c for a particular star. For 4U 1538-52 the solution behaves well for \(0.1574 \le c \le 0.46\) and for the values of a,  b,  A,  B,  M,  R given in Table 1, which corresponds to \(1\ge v_{r0}^2 \ge 0.13\), \(0.91 \ge v_{t0}^2 \ge 0.04\) and \(17.8 \ge \Gamma _0 \ge 3.8\). For LMC X-4 the solution behaves well for \(0.1235 \le c \le 0.35\) and for the values of a,  b,  A,  B,  M,  R given in Table 1, which yield \(0.99\ge v_{r0}^2 \ge 0.15\), \(0.91 \ge v_{t0}^2 \ge 0.06\) and \(15.35 \ge \Gamma _0 \ge 3.35\). Finally for PSR J1614-2230 the solution behaves well for \(0.05 \le c \le 0.13\) along with the values of a,  b,  A,  B,  M,  R given in Table 1, corresponding to \(1\ge v_{r0}^2 \ge 0.21\), \(0.94 \ge v_{t0}^2 \ge 0.13\) and \(8.34 \ge \Gamma _0 \ge 2.34\). Hence, we can conclude that smaller values of c lead to a stiffer equation of state and vice versa. The calculated masses and radii of the present stars are well fitted with those provided by Gangopadhyay et al. [58].