1 Introduction

Einstein’s theory of general relativity (GR hereafter) was tremendously successful in describing gravitational phenomena on both astrophysical and cosmological scales. This theory of gravity leads to significant breakthroughs through the intrinsic gravitational systems such as stellar and galactic structures that make up a major part of our visible Universe. The evolution investigation of these self-gravitating systems assumes a fundamental role in uncovering different hidden aspects viz., composition, evolution, and age of the Universe. One of the main phases in the evolution of celestial structures is the process of gravitational collapse near its death which causes the formation of compact stellar objects. These compact stellar objects are the end-points in the evolution of ordinary celestial structures which prove to be ideal sources for studying the properties and composition of high dense matter. Recently, various compact stellar objects with high densities were discovered [1], which are often seen as pulsars, rotating stars with strong magnetic fields. In the present-day, our theoretical understanding of compact stellar configurations was rooted in the Fermi-Dirac statistics due to their responsibility for the high degeneracy pressure in different astrophysical issues, that prevent the stellar structures from gravitational collapse and has been pointed out by Fowler in 1926 [2]. Hereafter, Chandrasekhar [3, 4] exhibited that white dwarfs are compact stellar structures, which are upheld exclusively by a degenerate gas of electrons, to be stable if the greatest size of a stable white dwarf is around 1.4 times the mass of the Sun (approximately \(3\times 10^{30}~\hbox {kg}\)), utilizing Einstein’s special theory of relativity and the principles of quantum physics.

Nowadays, there is no comprehensive characterization of immensely dense matter in a strongly interacting system. A likely theoretical characterization of such nuclear matter in ultrahigh densities may consist not solely of nucleons and leptons, but likewise many exotic constituents in their diverse forms and stages like baryon resonances, mesons, hyperons, as well as strange quark matter.

In many respects, the compact stellar structures have been attracted a lot of interest amongst the researchers. Despite the fact that there is significantly more to be examined and investigated about the compact celestial bodies, nevertheless, they are supposed with immense masses and having small radii, which makes them ultra-high dense bodies. Exploring the exact solutions in the celestial bodies situation up being a huge evolution in gravitational physics. It was thought that the isotropic fluid was used to form celestial structures, but, recently the anisotropy fluid has been attracted much attention among researchers under several astrophysical systems such as neutron stars [5], boson stars [6], and gravastars [7] confirming that the anisotropy of spherically symmetric compact stellar structures play an important role in both phenomenon astrophysics and cosmology.

It was examined in Ref. [8] that, in spite of the spherically symmetric distribution of matter under a compact celestial object, it can be portrayed by local pressure anisotropy. A study of the more general hydrostatic equilibrium equations, considering pressure anisotropy, demonstrates that anisotropy can have a fundamental effect on the maximum equilibrium mass and gravitational surface red-shift [8,9,10]. At not very high densities, the effect of anisotropy can be analyzed in the context of Newtonian gravity [11, 12]. At higher densities i.e., when \(\rho \ge 10^{15}~\hbox {g}/\hbox {cm}^{3}\), both the GR effects and the relativistic effects become more important [7, 13,14,15,16]. It is also worthwhile to mention here that for exploring the local pressure anisotropy effect on the specific basis, it is obligatory to recognize the substantial physical reasons liable for its occurrence, such as, e.g., the existence of a strong astrophysical core [15, 16], appearance of spontaneous deformation of Fermi surfaces [17, 18], accessibility of super-fluid states with the limited orbital momentum of Cooper pairs [19,20,21,22], or limited super-fluid momentum [23, 24], or the existence of strong magnetic fields within a stellar configuration [25,26,27,28,29,30,31,32,33].

The astrophysical stellar spheres were also analyzed by employing the equation of state (EoS) in the context of anisotropic pressure [8]. The researchers were affected due to the possible merging of the dark energy paradigm for the exploration of compact stellar bodies. The exploration of the physically stable systems leads us to an analytical mechanism related to Einstein’s field equations (EFEs). One of the fundamental mechanisms consists to match the embedding class one space-time which converts a 4-dimensional variety into a higher dimensional Euclidean space. Substantial new exact stellar models to evolve in the domain of astrophysics have been obtained through this conversion of curved embedding class space-time into space-time of higher dimension. The embedding class one condition drives to a differential equation in the background of spherically symmetric space-time which linked both metric potentials viz., \(g_{rr}\) and \(g_{tt}\). This class one embedding condition is also well-known as the Karmarkar condition. This famous condition seems to be powerful and highly prestigious in the study of new stellar solutions for astrophysical systems. In fact, Schlai [34] argues that a Riemannian variety with positive defined and analytic metric can be locally and isometrically embedded as a sub-variety of a Euclidean space. The first overall embedding theorem based on the isometrics of the Riemannian variety into Euclidean space has been shown by Nash [35]. Recently, Maurya and his collaborators [36,37,38,39,40,41,42,43] were the first explorers who applying the Karmarkar condition to the anisotropic matter distributions. Other related works can be found in Refs. [44,45,46,47,48,49].

In this paper, we are exploring a new family of embedded class space-time solutions admitting field equations, Karmarkar condition, and Pandey–Sharma criterion using spherically symmetric space-time and compare with the Schwarzschild model as an exterior space-time. The physical viability of the generating stellar model for anisotropic compact stars is confirmed via performing several physical tests of the main salient features such as energy density, radial and tangential pressures, anisotropy effect, dynamical equilibrium, energy conditions (ECs), and dynamical stability, which can be compared with the abundant observational data from different compact star candidates like PSR J1416-2230 with \(1.667 \pm 0.021\) [61], PSR J1903+327 with \(1.97\pm 0.04\) [62], 4U 1820-30 with \(1.58 \pm 0.06\) [63], Cen X-3 with \(1.49 \pm 0.08\) [64].

The present paper is organized as follows: After an exhaustive introduction in Sect. 1, we briefly discuss the basic principles of field equations for anisotropic matter in Einsteinian gravity in Sect. 2. In Sect. 3 we write down the equations describing the general solutions of an anisotropic star under class one space-time in Einsteinian gravity. The Sect. 4 includes the discussion of matching conditions of spherical symmetric space-time with Schwarzschild’s model as an exterior space-time. The physical investigation was performed comprehensively in Sects. 5 and 6. Then in Sect. 7, we develop some expressions of the embedding class one solution. Finally, in Sect. 8, conclusions and astrophysical implications are reported.

2 Basic stellar equations for anisotropic matter in Einsteinian gravity

We consider the interior space-time line element for a static and spherically symmetric fluid in Schwarzschild coordinates (\(x^i=t,r,\theta , \phi \)) as

$$\begin{aligned} ds^{2}=e^{\nu }dt^{2}-e^{\lambda }dr^{2}-r^{2}\left( d\theta ^{2}+\sin ^{2}\theta d\phi ^{2}\right) , \end{aligned}$$
(1)

where \(\nu =\nu (r)\) and \(\lambda =\lambda (r)\) are the metric potentials which have functional dependence on the radial coordinate r. For our astrophysical model, the stress-energy tensor for the anisotropic matter distribution within the stellar structure is taken to be

$$\begin{aligned} T^{\nu }_{\mu }=\left( \rho +p_{t}\right) U^{\nu }U_{\mu } -p_{t}\delta ^{\nu }_{\mu }-\left( p_{t}+p_{r}\right) V^{\nu }V_{\mu }, \end{aligned}$$
(2)

where \(\rho \), \(p_{r}\) and \(p_{t}\) are the proper energy density, the radial pressure and the transverse pressure of the stellar fluid in the orthogonal direction to \(p_{r}\), respectively. Here \(U^{\nu }\) stands the four-velocity \(e^{\nu (r)/2} U^{\alpha } =\delta ^{\alpha }_{0}\), while \(V^{\alpha }\) is a unit space-like vector in the radial direction \(V^{\alpha }=e^{-\lambda /2}\delta ^{\alpha }_{1}\), which is orthogonal to \(U^{\alpha }\). By considering this last comoving stellar fluid velocity, specifically, \(e^{\nu (r)/2} U^{\alpha } =\delta ^{\alpha }_{0}\), the EFEs for the line element (1) and the stress-energy tensor (2) leads directly to the following set of independent equations

$$\begin{aligned} 8\pi {\rho }= & {} \frac{1}{r^2}-e^{-\lambda }\Bigg [\frac{1}{r^2}-\frac{\lambda ^{\prime }}{r}\Bigg ], \end{aligned}$$
(3)
$$\begin{aligned} 8\pi {p}_{r}= & {} -\frac{1}{r^2}+e^{-\lambda }\Bigg [\frac{1}{r^2}+\frac{\nu ^{\prime }}{r}\Bigg ], \end{aligned}$$
(4)
$$\begin{aligned} 8\pi {p}_{t}= & {} \frac{1}{4}e^{-\lambda }\Bigg [2\nu ^{\prime \prime }+\nu ^{\prime 2}-\lambda ^{\prime }\nu ^{\prime } +2\frac{\nu ^{\prime }-\lambda ^{\prime }}{r}\Bigg ]. \end{aligned}$$
(5)

Here the primes denotes differentiation with respect to the radial coordinate r. Furthermore, we have considered units such that the speed of light c and the constant G are set to unity, i.e., \(c=G=1\).

So, by using Eqs. (4) and (5), the anisotropy parameter \(\Delta \) which measures the anisotropy inner the celestial body can be expressed as

$$\begin{aligned} \Delta \equiv {p}_{t}-{p}_{r} = \frac{e^{-\lambda }}{8\pi }\Bigg [\frac{\nu ^{\prime \prime }}{2} -\frac{\lambda ^{\prime }\nu ^{\prime }}{4}+\frac{\nu ^{\prime 2}}{4}+\frac{\nu ^{\prime }-\lambda ^{\prime }}{2r} +\frac{e^{\lambda }-1}{r^{2}}\Bigg ]. \end{aligned}$$
(6)

3 The general solutions of an anisotropic star under class one space-time in Einsteinian gravity

3.1 Basic formulation of embedding class one

It is notable that the embedding of p-dimensional space \(V^{p}\) in a pseudo-Euclidean space \(E^{p}\) pulled in much thought as construed by the authors [50, 51]. For the situation where a p-dimensional space \(V^{p}\) can be isometrically inundated in [p + q]-dimensional space, where q is a minimum number of extra dimensions, at this stage \(V^{p}\) is considered to be q-class embedding. Usually, the static spherically symmetric line element expressed in (1) gives the 4-dimensional relativistic space-time that can be incorporated in flat space-time of class two i.e, when \(q=2\), which exhibits that it is incorporated in a 6-dimensional pseudo-Euclidean space. Then again, it ought to be noticed that one can uncover a feasible parametrization to consolidate the space-time given in (1) into a 5-dimensional pseudo-Euclidean space which prompts embedding class one [50,51,52] i.e, when \(q=1\). Moreover, in the two cases: static or non-static for a spherically symmetrical space-time to be class-one, the stellar system must be compatible with the accompanying important and appropriate conditions. In this respect, the symmetric tensor \(b_{ij}\) corresponding to a stellar system must be satisfied and determined under the related Gauss–Codazzi equations, respectively:

$$\begin{aligned}&{\mathcal {R}}_{ijhk}=\epsilon \left( b_{ih}b_{jk}-b_{ik}b_{jh}\right) , \end{aligned}$$
(7)
$$\begin{aligned}&\nabla _{h}b_{ij}-\nabla _{i}b_{hj}=0, \end{aligned}$$
(8)

where \(\epsilon =\pm 1\) is everywhere normal to the manifold is time-like \((+1)\) or space-like \((-1).\) At this stage, the Riemann constituents for the line element given (1) can be read explicitly as follows,

$$\begin{aligned}&{\mathcal {R}}_{0101}=-e^{\nu }\left( \frac{\nu ^{\prime \prime }}{2}-\frac{\lambda ^{\prime }\nu ^{\prime }}{4} +\frac{\nu ^{\prime 2}}{4}\right) ;\quad {\mathcal {R}}_{1313}=-\frac{r}{2}\lambda ^{\prime }\sin ^2\theta ;\nonumber \\&{\mathcal {R}}_{2323}=-\frac{r^{2}\sin ^{2}\theta }{e^{\lambda }}\Big (e^{\lambda }-1\Big );\quad {\mathcal {R}}_{1202}=0;\nonumber \\&{\mathcal {R}}_{0202}=-\frac{r}{2}\nu ^{\prime }e^{\nu -\lambda };\quad {\mathcal {R}}_{1303}=0;\nonumber \\&{\mathcal {R}}_{0303}=-\frac{r}{2}\nu ^{\prime }e^{\nu -\lambda } \sin ^2\theta ; \quad {\mathcal {R}}_{1212} =-\frac{r}{2}\lambda ^{\prime }. \end{aligned}$$
(9)

Subbing these Riemann constituents into Gauss equations (7) prompts directly to

$$\begin{aligned}&b_{01}b_{33}={\mathcal {R}}_{1303}=0;\quad b_{01}b_{22}={\mathcal {R}}_{1212}=0; \quad b_{00}b_{33}={\mathcal {R}}_{0303};\nonumber \\&b_{00}b_{22}={\mathcal {R}}_{0202};\quad b_{11}b_{33}={\mathcal {R}}_{1313};\quad b_{22}b_{33}={\mathcal {R}}_{2323}; \nonumber \\&b_{11}b_{22}={\mathcal {R}}_{1212};\quad b_{00}b_{11}={\mathcal {R}}_{0101}. \end{aligned}$$
(10)

These relationships expressed in (10) prompts directly to the expressions which can be written explicitly as follows,

$$\begin{aligned}&\left( b_{00}\right) ^{2}=\frac{\left( {\mathcal {R}}_{0202}\right) ^{2}}{{\mathcal {R}}_{2323}} \sin ^{2}\theta ;\quad \left( b_{11}\right) ^{2}=\frac{\left( {\mathcal {R}}_{1212}\right) ^{2}}{{\mathcal {R}}_{2323}}\sin ^{2}\theta ;\nonumber \\&\left( b_{22}\right) ^{2}=\frac{{\mathcal {R}}_{2323}}{\sin ^{2}\theta }; \quad \left( b_{33}\right) ^{2}=\sin ^{2}\theta ~{\mathcal {R}}_{2323}. \end{aligned}$$
(11)

In this respect, in order to find the relationship in terms of the Riemann constituents, we combine the last term of the expression (10) conjointly with the constituents of expression (11) which leads to

$$\begin{aligned} {\mathcal {R}}_{0202}{\mathcal {R}}_{1313}={\mathcal {R}}_{0101}{\mathcal {R}}_{2323}, \end{aligned}$$
(12)

put through to \({\mathcal {R}}_{2323}\ne 0\) well-known as Pandey–Sharma condition [53]. It should be seen that all the constituents are yielded in (11) satisfy the Codazzi equations expressed in (8). Then again, on account of an overall non-static spherically symmetric space-time, the connection between constituents for symmetric tensor \(b_{ij}\) and Riemann tensor \(R_{ijhk}\) can be read explicitly as follows,

$$\begin{aligned} b_{01}b_{22}={\mathcal {R}}_{1212} \quad \text{ and } \quad b_{00}b_{11}-\left( b_{01}\right) ^{2}={\mathcal {R}}_{0101}. \end{aligned}$$
(13)

Here \(\left( b_{01}\right) ^{2}=\sin ^{2}\theta \left( {\mathcal {R}}_{1202}\right) ^{2}/{\mathcal {R}}_{2323}\). In the present circumstance, the embedding class-one condition well-known as Karmarkar condition [51, 52] take the accompanying shape

$$\begin{aligned} {\mathcal {R}}_{0202}{\mathcal {R}}_{1313} ={\mathcal {R}}_{0101}{\mathcal {R}}_{2323}+{\mathcal {R}}_{1202}{\mathcal {R}}_{1303}. \end{aligned}$$
(14)

Albeit in our condition, the relationship (12) between Riemann constituents through the static spherically symmetric line element (1) will be comparable to expression (14). The condition expressed in (14) assumes an essential role for representing the space-time (1) to be class-one, which is also well-known as a necessary and sufficient condition. In this regards, we obtain the differential equation by embedding the Riemann constituents expressed in (14) which takes the following form,

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

with \(e^{\lambda }\ne 1\). Moreover, we find the relationship between the gravitational potentials through integrating equation (14) which lead to the following form

$$\begin{aligned} \nu (r)=2 \ln \left[ C+D\int \sqrt{e^{\lambda }-1}\,dr\right] , \end{aligned}$$
(16)

where C and D are an integration constants. Additionally, the stellar solution established by the relationship (16) is named an embedding class one solution for the line element expressed in (1). It is worth mentioning that the above approach has been used to model compact celestial bodies in different astrophysical scenarios.

Using the Karmarkar condition [51, 52] in the expression for anisotropy (6), we obtain

$$\begin{aligned} \Delta = \frac{\nu ^{\prime }}{32\pi \,e^\lambda } \bigg [\frac{2}{r} -\frac{\lambda ^{\prime }}{e^\lambda -1} \bigg ] \bigg [\frac{\nu ^{\prime } e^\nu }{2\,r\,D^2}-1 \bigg ]. \end{aligned}$$
(17)

At this stage, we ought emphasize that when \(\Delta =0\), the only limited solution simultaneously fulfilling the Karmarkar condition and pressure isotropy is the interior Schwarzschild solution. This solution endures different shortcomings including superluminal speeds under the interior of the stellar fluid. To this end we consider a solution portraying an anisotropic fluid distribution which will be taken up in the following section.

3.2 Relativistic embedding class one solution

Now, we have a heavenly system of equations comprising of 4-equations, specifically, the EFEs (3)–(6) and 5-unknowns viz., the matter content \(\{\rho , p_r, p_t\}\) and the geometry \(\{\nu , \lambda \}\). Consequently, in order to obtain the established solution of the stellar system, we have needed two specific conditions which are as per the following:

  • The first condition is focused on choosing a specific gravitational mass expression m(r) by that determining the EoS should take the particular shape: \(p=p(\rho )\).

  • The second condition is focused on building up a relationship between the gravitational potentials \(\nu \) and \(\lambda \). Subsequently, to get a well-comported stellar solution, the two gravitational potentials should meet some physical and mathematical precondition attainable that metric potentials can not be assembled arbitrarily.

In this regard, Lake [54] has carried out that the metric potential \(\nu \) should to be regular, finite, monotonic increasing function and free from any physical and mathematical singularities interior the celestial structure, which gives a physically admissible fluid celestial body solution of EFEs. In any case, \(e^{\lambda }=1+O(r^2)\) is inevitable for a physically plausible celestial configuration to be regular at the center. In this paper, we are going to characteristic a physically attainable form of the gravitational potential \(\lambda \) that satisfies the above requirements, at that point we use embedding class one condition to set up another gravitational potential \(\nu \) in order to solve completely the EFEs. Accordingly, we won’t consider any system EoS in our study. For this purpose, we suggest a new gravitational potential in order to discover a new shut structure of solutions of EFEs for anisotropic fluid celestial bodies,

$$\begin{aligned} \lambda =\ln [\alpha r^2 + (1+\beta r^4)^{2} ] - 2 \ln [(1+\beta r^4)], \end{aligned}$$
(18)

where \(\alpha \) and \(\beta \) are constants having the dimensions \([length]^{-2}\) and \([length]^{-4}\), respectively.

By combining (18) in (15) we obtain

$$\begin{aligned} \nu =2 \ln \Bigg [C+\frac{D}{2} \sqrt{\frac{\alpha }{\beta }} \frac{ \cos ({\sqrt{\beta }r^2)}}{ \sin ({\sqrt{\beta }r^2)}} \Bigg ] \end{aligned}$$
(19)

where C and D are constant parameters.

Consequently, the class one space-time reads explicitly as

$$\begin{aligned} ~ds^{2}&=\Bigg [C+\frac{D}{2} \sqrt{\frac{\alpha }{\beta }} \frac{ \cos ({\sqrt{\beta }r^2)}}{ \sin ({\sqrt{\beta }r^2)}} \Bigg ]^{2} dt^{2}\nonumber \\&\quad -\Bigg [\frac{1 + 2\beta r^4 + (\beta r^4)^{2}}{1 +\alpha r^2 + 2\beta r^4 + (\beta r^4)^{2}}\Bigg ]dr^{2}-r^{2}d\Omega ^{2}, \end{aligned}$$
(20)

where \(d\Omega ^{2}\equiv \sin ^{2}\theta ~d\phi ^{2}+d\theta ^{2}\) is the usual metric on two-sphere \({\mathcal {S}}^{2}\). On using (18) and (19), we can recast the matter content viz., density \(\rho \), radial pressure \(p_{r}\) and transverse pressure \(p_{t}\) as

$$\begin{aligned} 8\pi \rho&= \Bigg [3\alpha + \alpha ^2 r^2 -2\alpha \beta r^4 -5\alpha (\beta r^4)^{2}\Bigg ]\nonumber \\&\quad \times \Bigg [\left( 1 + \alpha r^2 +2\beta r^4 +(\beta r^4)^{2} \right) ^{2}\Bigg ]^{-1}, \end{aligned}$$
(21)
$$\begin{aligned} 8\pi p_r&= \Bigg [-2\alpha \sqrt{\beta }C+4\sqrt{\alpha \beta }D\left( 1+ \beta r^4\right) -D\alpha \frac{ \cos ({\sqrt{\beta }r^2})}{ \sin ({\sqrt{\beta }r^2)}}\Bigg ]\nonumber \\&\quad \times \Bigg [\left( 2\sqrt{\beta }C +D\sqrt{\alpha } \frac{\cos ({\sqrt{\beta }r^2)}}{ \sin ({\sqrt{\beta }r^2)}} \right) \nonumber \\&\quad \times \left( \left( 1 + \alpha r^2 +2\beta r^4 + (\beta r^4)^2 \right) ^{2} \right) \Bigg ]^{-1}, \end{aligned}$$
(22)
$$\begin{aligned} 8\pi p_t&= \Bigg [ \sqrt{\alpha }+\beta \sqrt{\alpha }r^4 \Bigg ] \nonumber \\&\quad \times \Bigg [4\sqrt{\beta } D +2\alpha \sqrt{\beta } r^2 -4\sqrt{\beta ^5}D r^8 \nonumber \\&\quad +2\sqrt{\alpha \beta }C\left( (3\beta r^4 -1)\left( 1+D\beta \frac{\cos ({\sqrt{\beta }r^2)}}{ \sin ({\sqrt{\beta }r^2)}}\right) \right) \Bigg ] \nonumber \\&\quad \times \Bigg [ \left( 2\sqrt{\beta }C +D\sqrt{\alpha } \frac{\cos {(\sqrt{\beta }r^2)}}{\sin ({\sqrt{\beta }r^2})} \right) \nonumber \\&\quad \times \left( \left( 1 + \alpha r^2 +2\beta r^4 + (\beta r^4)^2 \right) ^{2} \right) \Bigg ]^{-1}, \end{aligned}$$
(23)

Using Eqs. (22) and (23) we get the anisotropy factor \(\Delta =p_t-p_r\) as

$$\begin{aligned} \Delta&= \Bigg [ \left( \sqrt{\alpha }+\beta \sqrt{\alpha }r^4 \right) \Bigg \{4\sqrt{\beta } D +2\alpha \sqrt{\beta } r^2 -4\sqrt{\beta ^5}D r^8 + 2\sqrt{\alpha \beta }C\nonumber \\&\quad \times \left( (3\beta r^4 -1)(1+\beta D\frac{\cos ({\sqrt{\beta }r^2)}}{ \sin ({\sqrt{\beta }r^2)}} \right) \Bigg \} \nonumber \\&\quad +\Bigg \{2\alpha \sqrt{\beta }C -4\sqrt{\alpha \beta }D\left( 1- \beta r^4\right) +\alpha \beta \frac{\cos ({\sqrt{\beta }r^2)}}{ \sin ({\sqrt{\beta }r^2)}}\Bigg \}\Bigg ]\nonumber \\&\quad \times \Bigg [8\pi \left( 2\sqrt{\beta }C +\sqrt{\alpha }D \frac{\cos ({\sqrt{\beta }r^2)}}{ \sin ({\sqrt{\beta }r^2)}}\right) \nonumber \\&\quad \times \left( \left( 1 + \alpha r^2 +2\beta r^4 + (\beta r^4)^2 \right) ^{2} \right) \Bigg ]^{-1}. \end{aligned}$$
(24)

Differentiating Eqs. (21)–(23) respectively, we obtain the density, radial pressure and transverse gradients as

$$\begin{aligned} 8\pi \frac{d\rho }{ dr}&= \frac{}{}2\alpha r \Bigg [5\alpha + (\alpha ^2 +28\beta )r^2 + 6\alpha \beta r^4 + 56\beta ^2 r^5 \nonumber \\&\quad -20\beta ^2 r^6 + 17\alpha \beta ^2 r^8 - 40\beta ^3 r^9 \nonumber \\&\quad + 28\beta ^3 r^{10} - 20\beta ^4 r^{14} \Bigg ] \nonumber \\&\quad \times \Bigg [\left( 1 + \alpha r^2 +2\beta r^4 + \beta ^2 r^8 \right) ^{3}\nonumber \\&\quad \times \left( 1 + \alpha r^2 +2\beta r^4 + \beta ^2 r^8 \right) ^{3}\Bigg ]^{-1}, \end{aligned}$$
(25)
$$\begin{aligned} 8\pi \frac{dp_r}{dr}&= 2\sqrt{\alpha } r \Bigg [-4\sqrt{\alpha }\beta D^2\left( \alpha r^2 + (1 + \beta r^4)^2\right) \nonumber \\&\quad -4\sqrt{\beta } D\left( \alpha (1 - \beta r^4)+2\beta (r + \beta r^5)^2\right) \nonumber \\&\quad \times \left( 2\sqrt{\beta }C +\sqrt{\alpha }D \times \frac{\cos ({\sqrt{\beta }r^2})}{ \sin ({\sqrt{\beta }r^2})} \right) \nonumber \\&\quad +\left( \sqrt{\alpha ^3} + 4\sqrt{\alpha }\beta r^2 + 4\sqrt{\alpha } \beta ^2 r^6 \right) \nonumber \\&\quad \times \left( 2\sqrt{\beta }C +\sqrt{\alpha }D\frac{ \cos ({\sqrt{\beta }r^2)}}{ \sin ({\sqrt{\beta }r^2)}} \right) ^2\Bigg ]\nonumber \\&\quad \times \Bigg [\left( 1 + \alpha r^2 +2\beta r^4 + \beta ^2 r^8 \right) ^{2} \nonumber \\&\quad \times \left( 2\sqrt{\beta }C +\sqrt{\alpha }D \frac{\cos ({\sqrt{\beta }r^2)}}{\sin ({\sqrt{\beta }r^2)}} \right) ^2\Bigg ]^{-1}, \end{aligned}$$
(26)
$$\begin{aligned} 8\pi \frac{dp_t}{dr}&= \Bigg [4\sqrt{\alpha } r\Bigg ] \nonumber \\&\quad \times \Bigg [\sqrt{\alpha }\beta D^2\left( 2+ \alpha r^2 - 2 \beta ^4 r^8\right) \left( \alpha r^2 + (1 + \beta r^2 \right) \nonumber \\&\quad - \sqrt{\beta } D \left( 2\sqrt{\beta }C +\sqrt{\alpha }D \frac{\cos ({\sqrt{\beta }r^2})}{ \sin ({\sqrt{\beta }r^2})}\right) \nonumber \\&\quad \times \left( \alpha ^2 r^2(-1 + \beta r^4) + 4\beta r^2 (-3 + \beta r^4)(1 + \beta r^4)^3 \right. \nonumber \\&\quad \left. -\alpha (1 + \beta r^4)(3 + 13\beta ^2 r^8 ) \right) \nonumber \\&\quad + \left( 6\sqrt{\alpha }\beta r^2(1-\beta r^4) (1+\beta r^4)^2 + \sqrt{\alpha ^3}(1 + 3\beta ^2 r^8)\right) \nonumber \\&\quad \times \left( 2\sqrt{\beta }C +\sqrt{\alpha }D\frac{\cos ({\sqrt{\beta }r^2})}{\sin ({\sqrt{\beta }r^2})} \right) ^2 \Bigg ]\nonumber \\&\quad \times \Bigg [\left( 1 + \alpha r^2 +2\beta r^4 + \beta ^2 r^8 \right) ^{3} \nonumber \\&\quad \times \left( 2\sqrt{\beta }C +\sqrt{\alpha }D\frac{ \cos ({\sqrt{\beta }r^2})}{\sin ({\sqrt{\beta }r^2})} \right) ^2\Bigg ]^{-1}. \end{aligned}$$
(27)

The central density \(\rho _c\) and central pressure \(p_c\) are obtained as

$$\begin{aligned} \rho _c= & {} \rho (r=0) = \frac{3\alpha }{8\pi } >0, \end{aligned}$$
(28)
$$\begin{aligned} p_c= & {} p_r(r=0) = p_t(r=0) = \frac{+2\sqrt{\alpha }D - \alpha C}{8\pi C}>0. \end{aligned}$$
(29)

Now, as suggested by Zeldovich and Navikov [55], the inequalities \(\rho - p_r>0\), \(\rho - p_t>0\) must be satisfied in the interior of a compact stellar system which is well-known as dominant energy condition. This dominant energy condition proportionally gives \(\omega _r = p_r/\rho <1\), \(\omega _t =p_t/\rho <1\). Currently, at the center of the celestial body, the above inequalities give \(\omega _c = p_c/\rho _c<1\), which infers

$$\begin{aligned} \frac{+2 D -\sqrt{\alpha }C}{3\sqrt{\alpha }C} \le 1. \end{aligned}$$
(30)

From expressions (29) and (30) we obtain a range of ratio C/D as

$$\begin{aligned} \frac{1}{2\sqrt{\alpha }} \le \frac{C}{D} < \frac{2}{\sqrt{\alpha }}. \end{aligned}$$
(31)

4 Exterior space-time: junction conditions

In this section, we matched our interior space-time smoothly with the exterior space-time at the pressure-free boundary \(\Sigma \) (\(f=r-R=0\), where R is a radius) in order to depict the full configuration of the self-gravitating anisotropic compact celestial body. For this purpose, the outer space-time is considered to be vacuum space-time i.e., outer Schwarzschild solution which is written explicitly as,

$$\begin{aligned} ds^2 = \Bigg [1-{2M \over r}\Bigg ] dt^2 - \Bigg [1-{2M \over r}\Bigg ]^{-1}dr^2 -r^{2}d\Omega ^{2}, \end{aligned}$$
(32)

where, M is the total gravitational mass. Along these lines, matching the internal and outside space-time demands the conformity of the first and second fundamental forms. So, these joining conditions are well-known as Darmois–Isarel junction conditions [56, 57]. The first fundamental form comprises in the continuity of the gravitational potential and its derivative across the boundary \(\Sigma \), which can be read explicitly as,

$$\begin{aligned}&{[}{ds^2}_{-}]_{\Sigma } = [{ds^2}_{+}]_{\Sigma }, \end{aligned}$$
(33)
$$\begin{aligned}&e^{\lambda ^{-}}|_{r=R} = e^{\lambda ^{+}}|_{r=R}\quad \text {and} \quad e^{\nu ^{-}}|_{r=R} = e^{\nu ^{+}}|_{r=R}, \end{aligned}$$
(34)

and

$$\begin{aligned} \Big (\frac{\partial e^{\nu ^{-}}}{\partial r} \Big )_{|r=R} =\Big (\frac{\partial e^{\nu ^{+}}}{\partial r} \Big )_{|r=R}. \end{aligned}$$
(35)

Here, as usual − and \(+\) represent the inner and outer space-times, respectively. However, by employing the continuity of the first fundamental form viz., \([ds^2]_{\Sigma }=0\), we shall constantly obtain

$$\begin{aligned} {[}F]_{\Sigma } \equiv F (r\longrightarrow R^{+}) -F (r\longrightarrow R^{-})\equiv F^{+}(R)-F^{-}(R) \end{aligned}$$
(36)

for any function F(r). This condition supplies us

$$\begin{aligned} g^{-}_{rr}(R)=g^{+}_{rr}(R)\quad \text {and}\quad g^{-}_{tt}(R)=g^{+}_{tt}(R). \end{aligned}$$
(37)

Then again, the space-time (1) should fulfill the second fundamental form at the surface \(\Sigma \) which is equivalent to the O’Brien and Synge [58, 59] junction condition which is defined as

$$\begin{aligned} {[}{K_{ij}}_{-}]_{\Sigma }=[{K_{ij}}_{+}]_{\Sigma } \end{aligned}$$
(38)

while \(K_{ij}\) describes the curvature. This last condition leads directly to

$$\begin{aligned} {[}G_{ij}r^j]_{\Sigma }=0, \end{aligned}$$
(39)

where \(r^j\) symbolizes a unit radial vector. The EFEs conjointly with the condition expressed in (39) provides

$$\begin{aligned} {[}T_{ij}\, r^j]_\Sigma =0\quad \Longrightarrow \quad p_r(R)=0. \end{aligned}$$
(40)

This condition establishes the celestial body size. This is so due to the pressure diminishes as we proceed toward to the surface and the pressure at the outside of the celestial body should vanish, at that point, this will correspond to the stellar structure boundary. On the other hand, the second fundamental form proposes that the matter distribution is bound in a limited spacetime arena, as a result, the stellar object doesn’t extend indefinitely beyond \(\Sigma \). Therefore, from the first fundamental form, we can obtain

$$\begin{aligned}&\left( C+\frac{D}{2} \sqrt{\frac{\alpha }{\beta }} \frac{ \cos ({\sqrt{\beta }R^2)}}{ \sin ({\sqrt{\beta }R^2)}} \right) ^2 = 1-{2M \over R}, \end{aligned}$$
(41)
$$\begin{aligned}&\frac{1 +\alpha R^2 + 2\beta R^4 + (\beta R^4)^{2}}{1 + 2\beta R^4 + (\beta R^4)^{2}} = 1-{2M \over R}, \end{aligned}$$
(42)
$$\begin{aligned}&\frac{DR\left( 2C+\sqrt{\frac{\alpha }{\beta }} D \frac{ \cos ({\sqrt{\beta }R^2)}}{ \sin ({\sqrt{\beta }R^2)}} \right) }{1+\beta R^4} ={2M \over R^2}, \end{aligned}$$
(43)
$$\begin{aligned}&\frac{-2\alpha R+6\alpha \beta R^5}{\left( 1+\beta R^4\right) ^{3}} = \frac{2M}{R^2 \left( 1-{2M \over R^2}\right) ^{2}}, \end{aligned}$$
(44)

and from the second fundamental form i.e. \(p_r(R)=0\), we obtain

$$\begin{aligned} \alpha \sqrt{\beta }C+4\sqrt{\alpha \beta }D\left( 1+ \beta R^4\right) -\alpha D \frac{\cos ({\sqrt{\beta }R^2)}}{ \sin ({\sqrt{\beta }R^2)}} = 0,\nonumber \\ \end{aligned}$$
(45)
Table 1 Constant parameters calculated for radii and mass for some compact star candidates

but here there is a discontinuous transverse pressure. Consequently, in order to avoid this discontinuity, we establish the surface stresses i.e., the surface energy density and surface pressure at the junction boundary by employing the Darmois–Israel [56, 57] condition. The intrinsic surface energy-momentum tensor \(S_{ij}\) is defined by Lanczos equations in the explicit form:

$$\begin{aligned} S^i_j=-\frac{1}{8\pi }\left( {\kappa }^i_j-{\delta }^i_j{\kappa }^k_k\right) , \end{aligned}$$
(46)

where the term \({\kappa }^i_j\) represent the extrinsic curvatures discontinuity, \(K^i_j\) beyond the hypersurface, i.e., \({\kappa }^i_j=K^+_{ij}-K^-_{ij}\). The extrinsic curvature of this boundary \(\Sigma \) can be reads explicitly as,

$$\begin{aligned} K^{\pm }_{ij}=-m_k^{\pm }\,\frac{\partial ^2y^k_{\pm }}{\partial n^i\,n^j}-m_k^{\pm } \,\Gamma ^k_{\mu \,l}\frac{\partial y^\mu _{\pm }}{\partial n^i}\frac{\partial y^l_{\pm }}{\partial n^j}. \end{aligned}$$
(47)

Here \(n^i\) denotes the coordinates in the boundary \(\Sigma \), while \(m_k^{\pm }\) represents the 4-velocity normal to \(\Sigma \). The constituents of this 4-velocity is defined in the coordinates (\(y^\nu _{\pm }\)) of \(\tau ^{\pm }\) as

$$\begin{aligned} m_k^{\pm }=\pm \,\frac{df}{dy^{k}}\,\big |g^{\mu \,l}\frac{df}{dy^{\mu }} \frac{df}{dy^{l}}\big |^{-\frac{1}{2}} \quad \text {with}\ m_km^k=1. \end{aligned}$$
(48)

For this purpose, since the nature of our space-time is spherically symmetric, the surface energy-momentum tensor can be defined as

$$\begin{aligned} S^i_j= diag\left( -\sigma , {\mathcal {P}}\right) , \end{aligned}$$
(49)

while \(\sigma \) and \({\mathcal {P}}\) represent the surface energy density and surface pressure, respectively. The surface energy density \(\sigma \) and the surface pressure \({\mathcal {P}}\) at the junction surface \(\Sigma \) are explicitly expressed by the following formulas:

$$\begin{aligned} \sigma&=-\frac{1}{4\pi R}\Bigg [\sqrt{e^{-\lambda }}\Bigg ]^+_- =-\frac{1}{4\pi R}\sqrt{1-\frac{2M}{R}}\nonumber \\&\quad +\frac{1+\beta R^4}{4\pi R \sqrt{1+\alpha R^2+2\beta R^4 + (\beta R^4)^{2}}}, \end{aligned}$$
(50)

and

$$\begin{aligned} {\mathcal {P}}&=\frac{1}{8\pi R}\Bigg [\left\{ 1+\frac{\sqrt{\alpha }}{2} \frac{\partial \nu }{\partial R}\right\} \sqrt{e^{-\lambda }}\Bigg ]^+_- = \frac{1-\frac{M}{R}}{8\pi R \sqrt{1-\frac{2M}{R}}}\nonumber \\&\quad -\frac{2 \sqrt{\alpha \beta }DR^2+(1+\beta R^4)(2 \sqrt{\beta }C+\sqrt{\alpha }D \frac{\cos ({\sqrt{\beta }R^2)}}{\sin ({\sqrt{\beta }R^2)}}}{8\pi R \sqrt{\alpha R^2+(1+\beta R^4)(2 \sqrt{\beta }C+\sqrt{\alpha }D \frac{ \cos ({\sqrt{\beta }R^2)}}{ \sin ({\sqrt{\beta }R^2)}}}}. \end{aligned}$$
(51)
Fig. 1
figure 1

The variation of metric potentials (left panel) and energy density (right panel) with the radial coordinate r for experimental statistics of four different compact stars. For plotting these graphs, the numerical values for the constant parameters are \(\alpha =0.030390~/\hbox {km}^{2}\), \(\beta =0.000075~/\hbox {km}^{4}\), \(C=0.9985\) and \(D=0.00991\)

Fig. 2
figure 2

The variation of radial and tangential pressures (left panel) and anisotropic factor (right panel) with the radial coordinate r for experimental statistics of four different compact stars. For plotting these graphs, the numerical values for the constant parameters are \(\alpha =0.030390~/\hbox {km}^{2}\), \(\beta =0.000075~/\hbox {km}^{4}\), \(C=0.9985\) and \(D=0.00991\)

5 Physical acceptability conditions for the anisotropic stellar models

Now, we need to carry out more physical tests for the anisotropic celestial structure of embedding class one. For this purpose, we want to fix the constants \(\alpha \), \(\beta \), C and D, in order to draw the graphs of the celestial model parameters first. In this respect, from Eqs. (41)–(44) we can observe that there are four equations with six parameters, viz., \(\alpha \), \(\beta \), C, D, M and R. So to determine the numeric values of the parameters, we use a datum set for different compact celestial bodies which are studied by Rawls et al. [64] for Cen X-3, Güver et al. [63] for 4U 1820-30, Demorest et al. [61] for PSR J1416-2230, Freire et al. [62] for PSR J1903+327. Employing the observational data of these compact stars, we can solve the four equations simultaneously and obtain the values for the constants \(\alpha \), \(\beta \), C and D of each celestial body as given in Table 1. Now we will discuss the conditions which are generally well-known to be crucial for anisotropic fluid stellar systems.

5.1 Regularity

  • Metric functions at the center, \(r=0\): From the gravitational potentials expressed in (18) and (19), we found that their central expressions can be written as \(e^{\lambda (r)}|_{r=0}=1\) and \(e^{\nu (r)}|_{r=0}=C^2\) respectively. This obviously indicates that the gravitational potentials are positive and finite at the center, which validates that our stellar system is free from any physical and geometrical singularities. In addition, Fig. 1 (left panel) asserts that the stellar system has finely increasing functions.

  • Density at the center, \(r = 0\): The energy density as a function as the radial coordinate r is represented in Fig. 1 (right panel). It is easy to observe that the matter density takes a finite value at the center and decreases monotonically towards the celestial body boundary.

  • Pressure at the center, \(r = 0\): From Fig. 2 (left panel), we can show that both the radial and transverse pressures are regular at the center of the celestial structure and decrease gradually towards the stellar surface. The disappearance of the radial pressure (\(p_r\)) at such a finite value, \(r = r_\Sigma =R\) defines the boundary of the fluid spherical body. It is widely believed that if the radial pressure disappears at the stellar object boundary, the transverse pressure is not necessary. A phenomenological justification was discussed by Boonserm et al. [78]. They argued that if the transverse pressure does not disappear at the boundary of the celestial configuration, then corresponding to the Schwarzschild outside would imply that the electric field is discontinuous at the star’s surface.

  • Anisotropy parameter : The anisotropy parameter is shown in Fig. 2 (right panel), and it is easy to see that \(\Delta > 0\) at each inside point of the celestial structure. The anisotropy parameter disappears at the center of the stellar object and increases to a maximum for such a finite radius, \(r<r_{\Sigma }\) where \({\Sigma }\) represents the boundary of the stellar body. A positive incentive for \(\Delta \) (\(p_t > p_r\)) implies a repulsive force due to anisotropy. This repulsive anisotropic force is necessary to hold the compact celestial body from collapsing. It is worth mentioning here that the increase in the anisotropy parameter particularly towards the surface layers leads to greater stability in these zones (Fig. 3).

With respect to the EoS parameters, namely, \(\omega _r=p_r/\rho \) and \(\omega _t=p_t/\rho \) they must have their maximum values less than 1, signifying that the matter content is non-exotic and proves that Zeldovich’s condition is satisfied at each point interior the stellar structure. With the help of a graphical representation of EoS parameters in Fig. 4 (left panel), we have shown that our celestial model confirms that the matter content is non-exotic and Zeldovich’s condition is well-satisfies.

Furthermore, Table 2 shows the predicted values of the physical salients of different compact star candidates corresponding to central and surface density which is according to the expected ranges for a stellar configuration formed by a non-exotic matter, also the radial pressure is shown at the center of the stellar structure.

Fig. 3
figure 3

The variation of ingredients of energy density, radial pressure (left panel) and tangential pressure (right panel) with the radial coordinate r for experimental statistics of four different compact stars. For plotting these graphs, the numerical values for the constant parameters are \(\alpha =0.030390~/\hbox {km}^{2}\), \(\beta =0.000075~/\hbox {km}^{4}\), \(C=0.9985\) and \(D=0.00991\)

Fig. 4
figure 4

The variation of equations of states viz. \(p_r/\rho \) and \(p_t/\rho \) (left panel) and energy conditions (right panel) with the radial coordinate r for experimental statistics of four different compact stars. For plotting these graphs, the numerical values for the constant parameters are \(\alpha =0.030390~/\hbox {km}^{2}\), \(\beta =0.000075~/\hbox {km}^{4}\), \(C=0.9985\) and \(D=0.00991\)

Table 2 Some physical parameters calculated for radii and mass for some compact star candidates
Fig. 5
figure 5

The variation of mass function (left panel) and compactness parameter and surface red-shift (right panel) with the radial coordinate r for experimental statistics of four different compact stars. For plotting these graphs, the numerical values for the constant parameters are \(\alpha =0.030390~/\hbox {km}^{2}\), \(\beta =0.000075~/\hbox {km}^{4}\), \(C=0.9985\) and \(D=0.00991\)

5.2 Energy conditions

For the physical acceptability of our celestial model, the solution must be able to describe a physical/non-exotic fluid that should satisfy the following ECs, namely, null energy condition (NEC), weak energy condition (WEC) and strong energy condition (SEC). To fulfill these ECs, the accompanying inequalities should be hold simultaneously at each interior point of the anisotropic matter distribution:

$$\begin{aligned}&\text {NEC:} \quad \rho \ge 0, \end{aligned}$$
(52)
$$\begin{aligned}&\text {WEC:} \quad \rho +p_{t}\ge 0, \ \rho +p_{r}\ge 0, \end{aligned}$$
(53)
$$\begin{aligned}&\text {SEC:} \quad \rho +2p_{t}+p_{r}\ge 0, \end{aligned}$$
(54)
$$\begin{aligned}&\text {DEC:} \quad \rho -|p_{r}|\ge 0, \ \rho -|p_{t}|\ge 0. \end{aligned}$$
(55)

It is obvious from Fig. 4 (right panel) that all four ECs expressed in inequalities (52)–(55) are fulfilled at each interior point of the stellar body. In this way, we have concluded that our stress-energy tensor is well-behaved.

5.3 Mass function and compactness factor

Let us now turn our consideration towards the viable mass-to-radius relationship. In this way, for a static spherically symmetric perfect fluid stellar structure, Buchdahl [60] has suggested an outright limitation on the maximally suitable mass-to-radius ratio for isotropic fluid celestial configurations as \(M/R \le 4/9\) (in the units \(c=G=1\)). This fundamentally expresses that for a given radius a static isotropic fluid stellar structure can’t be arbitrarily huge. Moreover, for a more generalized formula for the mass-to-radius ratio, one refers to the work studied by Mak and Harko [66]. In this regard, for the current compact celestial model, we can obtain the mass-to-radius ratio from the relationship between \(e^{\lambda }\) and the mass function m(r), i.e.

$$\begin{aligned} e^{-\lambda }=1-\frac{2m(r)}{r}. \end{aligned}$$
(56)
Fig. 6
figure 6

The variation of gravitational red-shift function (left panel) and relativistic adiabatic index (right panel) with the radial coordinate r for experimental statistics of four different compact stars. For plotting these graphs, the numerical values for the constant parameters are \(\alpha =0.030390~/\hbox {km}^{2}\), \(\beta =0.000075~/\hbox {km}^{4}\), \(C=0.9985\) and \(D=0.00991\)

Thereby, from this expression, we get the relationship of the mass function for the compact stellar structure as follows

$$\begin{aligned} m(r)=\frac{1}{2}\alpha r^3\Bigg [1 + \alpha r^2 +2 \beta r^4 + (\beta r^4)^2\Bigg ]^{-1}. \end{aligned}$$
(57)

The compactness parameter u of the celestial body is therefore given as

$$\begin{aligned} u(r)=\alpha r^2\Bigg [1 + \alpha r^2 +2 \beta r^4 + (\beta r^4)^2\Bigg ]^{-1}. \end{aligned}$$
(58)

Now from the relationship (58), we can note that the compactness parameter u of the stellar object depends upon the mass function. The compactness parameter increases with the increase of mass, as well as their corresponding value u satisfies the Buchdahl maximal allowable mass-to-radius ratio i.e., cannot be more than 4/9. In this connection, we have shown graphically the variation of the physical quantities related to Buchdahl’s mass-to-radius ratio for isotropic fluid stellar configurations, and also mass function are plotted in Fig. 5.

5.4 Surface and central red-shift

Consequently, the surface red-shift \(Z_s\), corresponding to the compactness parameter u expressed in Eq. (58) can be determined as follows,

$$\begin{aligned} z_{s}=e^{\lambda (R)/2}-1=\frac{1-\sqrt{1-2u_s}}{\sqrt{1-2u_s}}, \end{aligned}$$
(59)

explicitly it reads

$$\begin{aligned} z_{s}=\frac{\sqrt{1 + \alpha R^2 + 2 \beta R^4 +(\beta R^4)]^2}}{\sqrt{1 - \alpha R^2 + 2 \beta R^4 + (\beta R^4)]^2}} - 1. \end{aligned}$$
(60)

Now, the gravitational red-shift of the compact celestial body is given explicitly by

$$\begin{aligned} z=e^{\nu (r)/2}-1=\Bigg [C+\frac{D}{2} \sqrt{\frac{\alpha }{\beta }} \frac{\cos ({\sqrt{\beta }r^2)}}{ \sin ({\sqrt{\beta }r^2)}}\Bigg ]^{-1} - 1, \end{aligned}$$
(61)

and the central gravitational red-shift can be obtained as

$$\begin{aligned} z_{c}=\sqrt{e^{-\nu (0)}}-1=\frac{1}{C}-1, \end{aligned}$$
(62)

which must be non-negative interior the compact astrophysical configuration, i.e., \(\frac{1}{C}-1 >0\) implies \(C<1\).

Now,

$$\begin{aligned} \frac{dz}{ dr} = -\sqrt{\alpha }C r \Bigg [1 + \beta r^4 \Bigg ]^{-1} \Bigg [C +\frac{1}{2}\sqrt{\frac{\alpha }{\beta }}D \frac{\cos ({\sqrt{\beta }r^2)}}{ \sin ({\sqrt{\beta }r^2)}} \Bigg ]^{-2},\nonumber \\ \end{aligned}$$
(63)

at \(r=0\)

$$\begin{aligned} \frac{dz}{ dr} = 0\quad \text {with}\ \frac{d^2z}{ dr^2} = -\sqrt{\alpha }\frac{D}{C^2}<0. \end{aligned}$$
(64)

More importantly, a positive anisotropy parameter \(\Delta >0\) (as in our case) does not impose an upper limit on the red-shift surface \(z_{s}\), different is the case with isotropic distributions, where the most maximum value that the red-shift surface \(z_{s}\) can hold out as \(z_{s}=4.77\) [67]. Consequently, the surface red-shift for anisotropic matter distributions is greater than its isotropic partner. Figures 5 (right panel) and 6 (left panel) are shown that both the gravitational red-shift z and surface red-shift \(z_{s}\) are positive and bounded within the stellar structure. On the other hand, it is well-known that the maximum surface red-shift can’t surpass 3.842 if an anisotropic fluid distribution is taken into account as pointed out by Ivanov [68]. In this way, from our stellar solution, we can remark that throughout the astrophysical distribution, the resulting surface red-shift is \(z_s=0.2981592319\), and therefore our celestial model is well-consistent with Ivanov’s findings.

6 Equilibrium conditions and stability for the anisotropic stellar models

For completeness, we would also like to explore the stability as well as the equilibrium conditions in the interior of the fluid sphere.

6.1 Equilibrium under three different forces

It is well known that the study of the compact celestial bodies equilibrium depends strongly on the Tolman–Oppenheimer–Volkoff (TOV) equation [69, 70]. In light of the TOV equation, we want to explore whether our current celestial model is in a stable equilibrium state under the three following forces:

  1. 1.

    The gravitational force \(F_{g}\),

  2. 2.

    The hydrostatic force \(F_{h}\),

  3. 3.

    The anisotropic repulsive force \(F_{a}\) described by the existence of a positive anisotropy factor \(\Delta \).

Fig. 7
figure 7

The variation of radial \((v_r^2)\) and transverse \((v_t^2)\) velocity of sound and stability factor \((|v_r^2 - v_t^2|)\) (left panel) and different forces in TOV equation for a static configuration (right panel) with the radial coordinate r for experimental statistics of four different compact stars. For plotting these graphs, the numerical values for the constant parameters are \(\alpha =0.030390~/\hbox {km}^{2}\), \(\beta =0.000075~/\hbox {km}^{4}\), \(C=0.9985\) and \(D=0.00991\)

As shown previously, the existence of a positive anisotropy factor counterbalances the gravitational gradient. Consequently, this deduces that the sum of three different forces becomes zero

$$\begin{aligned} +\underbrace{\frac{2}{r}\Delta }_{F_{a}} - \underbrace{\frac{dp_{r}}{dr}}_{F_{h}} -\underbrace{\frac{\nu ^{\prime }}{2}\left( \rho +p_{r}\right) }_{F_{g}}=0. \end{aligned}$$
(65)

The explicit formulas of these three forces are given as

$$\begin{aligned} F_{a}&=\frac{2}{8\pi r f_{1}(r) \left( \left( 1 + \alpha r^2 +2\beta r^4 + (\beta r^4)^2 \right) ^{2} \right) } \nonumber \\&\quad \times \Bigg [ \left( \sqrt{\alpha }+\beta \sqrt{\alpha }r^4 \right) \Bigg \{4\sqrt{\beta } D +2\alpha \sqrt{\beta } r^2 -4\sqrt{\beta ^5}D r^8 \nonumber \\&\quad + 2\sqrt{\alpha \beta }C\left( (3\beta r^4 -1)(1+\beta D \frac{ \cos ({\sqrt{\beta }r^2)}}{ \sin ({\sqrt{\beta }r^2)}}) \right) \Bigg \}\nonumber \\&\quad +\Bigg \{2\alpha \sqrt{\beta }C -4\sqrt{\alpha \beta }D \left( 1- \beta r^4\right) +\alpha \beta \frac{ \cos ({\sqrt{\beta }r^2)}}{ \sin ({\sqrt{\beta }r^2)}}\Bigg \}\Bigg ]. \end{aligned}$$
(66)
$$\begin{aligned} F_h(r)&=\frac{2\sqrt{\alpha } r}{8\pi \left( 1 + \alpha r^2 +2\beta r^4 + \beta ^2 r^8 \right) ^{2} \left( f_{1}(r) \right) ^2} \nonumber \\&\quad \times \Bigg [-4\sqrt{\alpha }\beta D^2\left( \alpha r^2 + (1 + \beta r^4)^2\right) -4\sqrt{\beta } D\nonumber \\&\quad \times \left( \alpha (1 - \beta r^4) + 2\beta (r + \beta r^5)^2\right) f_{1}(r) \nonumber \\&\quad + \left( \sqrt{\alpha ^3} + 4\sqrt{\alpha }\beta r^2 +4\sqrt{\alpha } \beta ^2 r^6 \right) \left( f_{1}(r) \right) ^2 \Bigg ], \end{aligned}$$
(67)
$$\begin{aligned} F_g(r)&= -\frac{4\alpha \sqrt{\beta }D r f_{2}(r)}{\left( f_{1}(r) \right) ^2 \left( 1 + \alpha r^2 +2\beta r^4 + (\beta r^4)^2 \right) ^{2}} , \end{aligned}$$
(68)

where,

$$\begin{aligned} f_{1}(r)&= \Bigg [2\sqrt{\beta }C +\sqrt{\alpha }D \frac{\cos ({\sqrt{\beta }r^2)}}{ \sin ({\sqrt{\beta }r^2)}} \Bigg ],\\ f_{2}(r)&=2\sqrt{\beta } \\&\quad \times \Bigg [\alpha D r^2 + D\left( 1 + 2\beta r^4 + (\beta r^4)^2 \right) +\sqrt{\alpha } C \left( 1 - 3\beta r^4 \right) \Bigg ] \\&\quad +\alpha D \left( 1 - 3\beta r^4 \right) \frac{\cos ({\sqrt{\beta }r^2)}}{ \sin ({\sqrt{\beta }r^2)}}. \end{aligned}$$

The Fig. 7 (right panel) demonstrates that gravitational force dominates both hydrostatic and anisotropic forces, and therefore the celestial model is in an equilibrium state as the gravitational force counterbalances joined impact of hydrostatic and anisotropic forces.

6.2 Causality condition and Abreu’s criterion

The speed of sound interior the compact celestial bodies can be established by employing

$$\begin{aligned} v_{r}(r)=\sqrt{\frac{dp_{r}(r)}{d\rho (r)}}, \quad v_{t}(r)=\sqrt{\frac{dp_{t}(r)}{d\rho (r)}}. \end{aligned}$$
(69)

Both the radial and transverse subliminal sound speed interior of the celestial structure ought to be less than the light speed, this phenomenon is well-known as a causality condition. As we can easily see that both these quantities indicate that our stellar model is totally stable, in light of the fact that the subliminal sound speed is less than the light speed (taking \(c=1\) in relativistic geometrized units) wherever within the celestial configuration.

We use another important criterion for the static spherically symmetric celestial structure dubbed the cracking criterion of stellar configuration suggested by Herrera [71], for which we can check whether the local anisotropic matter distribution is stable or unstable. Then, Abreu and co-workers [72] have suggested another choice for taken into account the stability of a self-gravitating anisotropic fluid sphere. More interesting, this methodology demonstrates whether the region is potentially stable where the radial sound velocity is greater than the tangential sound velocity. As per Abreu et al. [72], the region for stability compact celestial bodies can be determined via the following criterion,

$$\begin{aligned}&-1\le v^{2}_{t}-v^{2}_{r}\le 1 \nonumber \\&\quad \equiv \left\{ \begin{array}{ll} -1\le v^{2}_{t}-v^{2}_{r}\le 0 &{} \quad \text {Potentially stable } \\ 0< v^{2}_{t}-v^{2}_{r}\le 1 &{} \quad \text {Potentially unstable} \end{array}\right\} . \end{aligned}$$
(70)

It is easy to observe from Fig. 7 (left panel) that the stability factor viz., (\(v^{2}_{t}-v^{2}_{r}\)) for our stellar model satisfies the cracking criterion proposed by Herrera and Abreu and co-workers [71, 72] everywhere within the compact celestial body. Hence, we can come to the conclusion that our anisotropic stellar model is well-behaved and provides a stable configuration.

6.3 Relativistic adiabatic index

For a relativistic anisotropic stellar structure, the stability is also depends upon the adiabatic index \(\Gamma \), the ratio of two specific heats, established by the following relation as [74]

$$\begin{aligned} \Gamma =\frac{\rho +p_{r}}{p_{r}}\frac{dp_{r}}{d\rho }. \end{aligned}$$
(71)

From this point of view, the relativistic adiabatic index corresponding to the case \(\Gamma > 4/3\) leads to the condition for the stability of a Newtonian stellar configuration and for \(\Gamma =0\) being the condition for a neutral equilibrium as suggested by Bondi [75]. It is worth mentioning that this condition modifies for a relativistic isotropic stellar structure due to the regenerative impact of pressure, which delivers the stellar structure more unstable. For a relativistic anisotropic fluid stellar structure, the situation turns out to be more complicated due to the stability will rely upon the kind of anisotropy. In this respect, the stability condition via an anisotropic fluid stellar structure is given by

$$\begin{aligned} \Gamma >\frac{4}{3}\left[ 1+3\pi \frac{\rho _{0}p_{r0}}{|p^{\prime }_{r0}|}r +\frac{\left( p_{t0}-p_{r0}\right) }{|p^{\prime }_{r0}|r}\right] _{max}, \end{aligned}$$
(72)

while \(\rho _{0}\), \(p_{r0}\) and \(p_{t0}\) represent the initial energy density, initial radial pressure and initial transverse pressure, respectively, in the static equilibrium state which fulfills the TOV equation expressed in (65). Noting that the second and last terms in the square brackets represent the relativistic and anisotropic corrections, which increase the instability range of the adiabatic index. In conclusion, our outcomes have appeared in Fig. 6 (right panel), where we plot \(\Gamma \) as a function of radial coordinate, r. From this graph, the resulting \(\Gamma > 4/3 \approx 1.33\) exhibits that our celestial model is stable versus the radial adiabatic infinitesimal pulsations and increasing values of adiabatic index (\(\Gamma \)) mean the pressure development for a given increase in energy density, i.e., a stiffer EoS.

6.4 Harrison–Zeldovich–Novikov static stability criterion

In light of the above considerations, any solution describing stable astrophysical celestial structures should satisfy the stability criterion. This stability criterion determines whether the solution is static and stable under radial infinitesimal perturbations. It is worthwhile to mention here that according to the stability criterion of hydrostatic equilibrium configuration of the celestial structure, the equilibrium mass M of such structure varies with respect to its central density \(\rho _c\), which implies that the equilibrium celestial configuration with \(\partial M(\rho _c)/\partial \rho _c >0\) are stable, but those with \(\partial M(\rho _c)/\partial \rho _c <0\) are unstable [76, 77]. For our solution, the expression for equilibrium mass M against central density \(\rho _c\) is formulated via the following relation

$$\begin{aligned} M(\rho _c) = \frac{4 \pi \rho _c R^3}{3+ 8\pi \rho _c R^2 + 6\beta R^4 +3(\beta R^4)^{2}}, \end{aligned}$$
(73)

and

$$\begin{aligned} {\partial M \over \partial \rho _c} = \frac{4 \pi \rho _c R^3 \big [3+ 8\pi R^2 (\rho _c -1) +6\beta R^4 +3(\beta R^4)^{2}\Big ]}{\Big [3 + 8\pi \rho _c R^2 + 6\beta R^4 + 3(\beta R^4)^{2}\Big ]^2} >0. \end{aligned}$$
(74)

Importantly, we can easily see from Fig. 8 that the mass of the stellar structure grows with the increment of the central energy density. This outcome confirm the static stability criterion of the celestial model versus radial infinitesimal perturbations, and also we conclude that the compact stellar configuration turns out to be more massive with respect to increasing central density.

Fig. 8
figure 8

Variation of total mass with central density for experimental statistics of four different compact stars. For plotting these graphs, the numerical values for the constant parameters are \(\alpha =0.030390~/\hbox {km}^{2}\), \(\beta =0.000075~/\hbox {km}^{4}\), \(C=0.9985\) and \(D=0.00991\)

7 Herrera–Ospino–Di Prisco generators of the present embedding class one solution

The algorithm for all possible spherically symmetric static anisotropic solutions through two generating functions \(\zeta (r)\) and \(\Pi (r)\) for the EFEs has been established by Herrera and collaborators [73]. Both generators \(\zeta (r)\) and \(\Pi (r)\) are linked with the gravitational potential \(e^\nu \) and the pressure anisotropy, respectively. These both generators are formulated via the following relations,

$$\begin{aligned}&\exp \bigg [\int \bigg (2\zeta (r)-{2 \over r} \bigg ) dr \bigg ] = \exp \bigg [\nu (r)\bigg ], \end{aligned}$$
(75)
$$\begin{aligned}&\Pi (r) = 8\pi (p_r(r)-p_t(r)). \end{aligned}$$
(76)

Then, for our case, the solution of these generating functions are explicitly determined as

$$\begin{aligned} \zeta (r)&= \frac{1}{r} - \Bigg [\left( D \sqrt{\alpha } \left( 1 + \frac{\cos ({\sqrt{\beta }r^2)}}{\sin ({\sqrt{\beta }r^2)}} \right) \right) r \Bigg ] \nonumber \\&\quad \times \Bigg [ \left( C+\frac{D}{2} \sqrt{\frac{\alpha }{\beta }} \frac{\cos ({\sqrt{\beta }r^2)}}{\sin ({\sqrt{\beta }r^2)}} \right) \left( \frac{\sin ({\sqrt{\beta }r^2)}}{\cos ({\sqrt{\beta }r^2)}} \right) ^2\Bigg ]^{-1} \end{aligned}$$
(77)
$$\begin{aligned} \Pi (r)&= - \frac{1}{f_{1}(r) \left( \left( 1 + \alpha r^2 +2\beta r^4 + (\beta r^4)^2 \right) ^{2} \right) } \nonumber \\&\quad \times \Bigg [ \left( \sqrt{\alpha }+\beta \sqrt{\alpha }r^4 \right) \Bigg \{4\sqrt{\beta } D +2\alpha \sqrt{\beta } r^2 -4\sqrt{\beta ^5}D r^8 \nonumber \\&\quad + 2\sqrt{\alpha \beta }C\left( (3\beta r^4 -1)(1+\beta D \frac{\cos ({\sqrt{\beta }r^2)}}{ \sin ({\sqrt{\beta }r^2)}}) \right) \Bigg \}\nonumber \\&\quad +\Bigg \{2\alpha \sqrt{\beta }C - 4\sqrt{\alpha \beta }D \left( 1- \beta r^4\right) +\alpha \beta \frac{ \cos ({\sqrt{\beta }r^2)}}{ \sin ({\sqrt{\beta }r^2)}}\Bigg \}\Bigg ], \end{aligned}$$
(78)

where \(f_{1}(r)\) is defined above. It is worth mentioning here that the relevance of the present system from an astrophysical point of view, can be completely determined by one generating function via \(e^\nu \) and an additional ansatz in the form of embedding class one.

8 Conclusions and astrophysical implications

This work is devoted to exploring the possibility of providing a new family of exact solutions for viable anisotropic spherically symmetric systems in the background of general relativity. For this purpose, we emphasize the embedding class approach as a ground-breaking tool, for obtaining anisotropic solutions for matter sources through EFEs. To be precise, our study proceeds in the realm of class one of the embeddings, by approaches for embedding spherically symmetric static metric into the five-dimensional pseudo-Euclidean space, to determine a more general stellar solution to the EFEs. For the present analysis, we chose an appropriate forms for gravitational metric coefficient viz., \(g_{rr} \equiv \lambda (r)\) expressed in (18) that satisfies all the requirements carried out by Lake [54], and we have obtained the second gravitational metric coefficient \(g_{tt} \equiv \nu (r)\) given in (19) by means of the embedding class one approach via Karmarkar condition and Pandey–Sharma criterion. Further, we impose the boundary conditions, by choosing the Schwarzschild model i.e., vacuum space-time to describe as an exterior space-time, which helps us to find the unknown constraints. A brief qualitative analysis of all the obtained results is itemized below:

  • Foremost, it is well-known that the aspects of gravitational metric potentials play an important role in obtaining a well-behaved celestial solution. For this purpose, we have shown both gravitational metric coefficients versus radial coordinate r as illustrated in Fig. 1 (left panel), which shows that the fundamental conditions: \(e^{\lambda (r)}|_{r=0} = 1\) and \(e^{\nu (r)}|_{r=0} = C^2 \ne 0\), were surrounded by the metric potentials. Increasing monotonic attributes of gravitational metric potentials are also observed everywhere within the stellar configuration.

  • The evolution of energy density, radial and transverse pressures viz., \(\rho \), \(p_r\) and \(p_t\), can be observed in Fig. 1 (right panel) and Fig. 2 (left panel), respectively. The energy density along with radial and transverse pressures exhibit decreasing behavior i.e., behave realistically and are positive throughout the stellar configuration. The maximum value is achieved at the center while decreasing development has been observed with the increment in radii and it tends to zero towards the stellar surface. From the Figs. 1 and 2, we corroborate also that our stellar configuration is fully free from any physical or geometrical singularities.

  • The anisotropy behavior and the gradients of energy density, radial and transverse pressures according to the radial coordinate, r are represented in Figs. 2 (right panel) and 3, respectively. On the one hand, we observe that the \(p_r \ne p_t\) and \(p_t>p_r\), consequently, \(\Delta >0\), so the anisotropy is positive and implies that the celestial system encounters a repulsive force counteracts the gravitational gradient progressing the stability and equilibrium state. On the other hand, we noticed that all the derivatives of energy density, radial and transverse pressures are non-positive and show decreasing behavior i.e., \(\frac{d\rho }{ dr}<0\), \(\frac{dp_r}{ dr}<0\), \(\frac{dp_t}{ dr}<0\) . The negativity in gradients exhibits that our acquired solutions are physically agreeable.

  • We have argued the EoS parameters viz., \(\omega _r \equiv p_r/\rho \) and \(\omega _t \equiv p_t/\rho \) as well as four-types of energy bounds NEC, WEC, SEC, DEC in our study. It is obvious from Fig. 4 that both EoS parameters have their maximum values less than 1, which means that the matter content is non-exotic and proves that Zeldovich’s condition is well-satisfied throughout the stellar configuration. In the same graph, we noted that all the ECs are well-fulfilled in our case. In the same graph, we noticed that all the ECs are decreasing monotonic functions with increasing radii, which is even more satisfactory in our case.

  • We have also argued the compactness factor, mass function along with gravitational red-shift. The interesting behaviors of these physical quantities i.e., u(r), m(r), z have been shown in Figs. 5 and 6 (left panel) respectively. It can be seen that both u(r) and m(r) show increasing behavior as well as u(r) increases with the increase of m(r), and their corresponding value satisfies the Buchdahl maximal allowable mass-to-radius ratio i.e., cannot be more than 4/9. Apart from this, gravitational red-shift exhibits the decreasing attribute and the resulting surface red-shift is \(z_s = 0.30<3.482\), which is in alliance with Ivanov’s findings affirming the stability of our stellar system.

  • The equilibrium conditions via TOV equation for the present stellar configuration study have also been discussed. The balancing nature of the hydrostatic \(F_h\), gravitational \(F_g\) and anisotropic \(F_a\) forces is described in the right of Fig. 7. As all these forces i.e., \(F_h\), \(F_g\), \(F_a\) add up to zero and balance the impact of each other, demonstrates that our obtained anisotropic solutions for the compact stellar structures are stable and physically viable.

  • The stability analysis of a compact stellar configuration study is additionally a fundamental attribute. In our study, we analyzed the stability through relativistic adiabatic index, causality condition and Abreu’s criteria, and Harrison-Zeldovich-Novikov static stability criterion. The evolution of these fundamental criteria is presented in Figs. 6 (right panel), 7 (left panel) and 8 respectively. From the point of view of the causal condition and Abreu’s criterion, the stellar model of the anisotropic configuration is totally stable, due to the subliminal sound speed is less than 1 throughout the celestial configuration, as well as there is no change in sign \(v_t^2-v_r^2\) and stability factor (\(v_t^2-v_r^2\)) lies among \(-1\) and 0 for stable stellar configuration or otherwise unstable if it is between 0 and 1. Moreover, the relativistic adiabatic index \(\Gamma \) showing that all the curves involved remain in a zone \(>4/3\). The resulting \(\Gamma >4/3 \approx 1.33\) exhibits that our stellar model is stable against the radial adiabatic infinitesimal pulsations and increasing values of adiabatic index (\(\Gamma \)) mean the pressure evolution for a given increase in energy density, i.e., a stiffer EoS. Finally, the Harrison-Zeldovich-Novikov static stability criterion or \(M-\rho _c\) function plays a crucial role in ensuring the stability of spherically symmetric static stellar systems under radial pulsation has been well-satisfied. We can also notice from the data drawn in \(M-\rho _c\) curve that the stellar configurations become more massive according to increasing central density.

Finally, it would be interesting to mention here that all anisotropic spherically symmetric solutions demonstrate in the present paper satisfying obtained well-behaved celestial interiors by employing the embedding class one approach via Karmarkar condition and Pandey–Sharma criterion. It is also admitting and sharing all the physical and mathematical attributes necessary in the study of compact stars, which provide circumstantial evidence in favor of the evolution of realistic compact stars. In effect, our anisotropic stellar model supports the existence of realistic super-massive pulsars like PSR J1416-2230, PSR J1903+327, 4U 1820-30 and Cen X-3.