Study of a Minimally Deformed Anisotropic Solution for Compact Objects with Massive Scalar Field in Brans–Dicke Gravity Admitting the Karmarkar Condition

Abstract

In the present paper, we focused on exploring the possibility of providing a new class of exact solutions for viable anisotropic stellar systems by means of the massive Brans–Dicke (BD) theory of gravity. In this respect, we used the decoupling of gravitational sources by minimal geometric deformation (MGD) ($e^{-\eta }=\Psi +\beta h$) for compact stellar objects in the realm of embedding class-one space-time to study anisotropic solutions for matter sources through the modified Einstein field equations. For this purpose, we used the ansatz for $\Psi$ relating to the prominent, well-known and well-behaved Finch–Skea model via Karmarkar condition, and the determination scheme for deformation function $h\left(r\right)$ was proposed via mimic requirement on radial pressure component: ${\theta }_1^1\left(r\right)=p_r\left(r\right)$ and matter density: ${\theta }_0^0\left(r\right)=\rho \left(r\right)$ for the anisotropic sector. Moreover, we analyzed the main physical highlights of the anisotropic celestial object by executing several physical tests for the case ${\theta }_1^1\left(r\right)=p_r\left(r\right)$. We have clearly shown how the parameters $\alpha$, $\beta$ and ${\omega }_{BD}$ introduced by massive BD gravity via the MGD approach incorporating the anisotropic profile of the matter distribution have an immense effect on many physical parameters of compact bodies such as LMC X-4, LMC X-4, Her X-1, 4U 1820-30, 4U 1608-52, SAX J1808.4–658 and many others that can be fitted.

1. Introduction

The exhaustive investigation of several sets of high-precision observational data through Cosmic Microwave Background, SuperNova type Ia, Large Scale Structure, Weak Lensing, Baryon Acoustic Oscillations and so on underpins a picture wherein the Universe is spatially flat with solely about 30% of its full energy balance in the form of luminous or dark matter types. On the other hand, the nature of the remaining 70%, living in some obscure type of perplexing “dark energy”, remains a mystery, although it may be precisely portrayed by a specific kind of anti-gravitating stress. In this regard, discovering the physics of this enigmatic dark energy, leading to the present accelerated expansion of the Universe, is a crucial objective of modern physics, specifically astrophysics and cosmology.

The intensive discussion among astrophysicists and cosmologists has revealed the possibility of different new scenarios in astrophysical and cosmological research that provide various sources for the observed acceleration and even a possibility to connect it to the associated epoch of inflationary acceleration in the extremely precocious Universe. Consequently, the general relativity theory (GRT) as a bound-together portrayal of space-time and matter gives a basic insight into the astrophysical as well as cosmological phenomena.

The careful solutions of the Einstein field equations (EFEs) represent the intrinsic characteristics of self-gravitating stellar structures. In this respect, in 1916, Schwarzschild [1] found the first spherical vacuum solution which indicates the outside geometry of perfect matter distribution. Afterward, in 1939, Tolman [2] examined smooth coordinating conditions for the resulting inside solutions with the outside one and furthermore explored some accurate solutions for perfect fluid within sight of the cosmological constant. On the other hand, in 1933, Lemaitre [3] suggested that anisotropy in the interior of heavenly structures can occur due to rotational movement, phase transition, solid core or a mixture of two fluids, and so on. At this stage, various analytical solutions are available in the literature [4,5] which are physically acceptable solutions for both isotropic and anisotropic relativistic stellar objects. Recently, several researchers have taken into account the pressure via radial and tangential components to incorporate anisotropy in the configuration of celestial bodies. Herrera and Santos [6] explored the possible parameters that prompt anisotropy in spherical stellar systems in GRT. In the same spirit, Harko and Mak [7] inferred static solutions representing the anisotropic inner side of cosmic stellar structures via the anisotropy parameter. The same authors [8] show that the anisotropy is instigated in a stellar system due to the existence of a solid core. In this respect, there are several works accessible in the literature for spherically symmetric space-time with anisotropic fluid distributions that make it possible to contrast with astrophysical observations that have been derived in [9,10,11,12,13,14,15,16,17,18,19] (and the references therein).

In order to explore salient characteristics of compact objects, viz., mechanisms of stability and hydrodynamic equilibrium, causality conditions, maximum limit of the mass–radius ratio, maximum value of the superficial red-shift, the behavior of the material content via energy conditions and so on, different strategies have been proposed. One of these strategies is gravitational decoupling through the minimal geometric deformation approach, which was originally supposed as an elective way to deform Schwarzschild space-time in sight of astrophysics braneworld. It is worth noting here that this powerful approach was originally provided in 2008 by Ovalle [20] in order to draw up such analytical stellar solutions in astrophysics as well as cosmology. In this straightforward and powerful approach, an extra source is joined in the grain distribution, relying on the prerequisite that both sources solely associate with gravitation. In this regard, the system of field equations splits into two sets with lesser degrees of freedom as contrasted to the original system when we are basing it on a geometric deformation in the radial metric ingredient. The two systems are tackled separately and their corresponding solutions are joined to acquire a solution of the full system. According to this mechanism, the authors Ovalle and Linares [21] applied the MGD approach in order to calculate an interior exact solution for an anisotropic compact celestial object; they also found that their solution is compatible with the Tolman IV solution in the background of astrophysics braneworld as well as assessing the mass consequences for the compactness of self-gravitating stellar structures. Ovalle and his colleagues [22] have also employed this approach to generate three anisotropic stellar systems via the Tolman IV solution after incorporating the anisotropy impacts in the perfect fluid distribution. Otherwise, this MGD approach has been performed to cover several other subjects such as models for neutron stars [23], anisotropic-like Tolman IV solutions [24,25], anisotropic-like Tolman VII solutions [26] as well as in alternative and modified theories of gravitation, viz., $f\left(G\right)-$gravity [27], Lovelock−gravity [28], $f(R,T)-$gravity [29,30] and Rastall−gravity [31], among others.

On the other hand, this MGD approach has been used by Cavalcanti et al. [32] to explore the gravitational lensing singularity beyond to the GRT, and by Casadio and his collaborators [33] to find the basic stability space for Bose–Einstein condensates in gravitational structure. In sight of variable tension fluid branes, it was additionally applied to examine the corrections to dark SU(N) stellar structure evident parameters [34].

The approach of gravitational decoupling by MGD is a greatly efficient procedure to construct solutions of EFEs for self-gravitating celestial bodies. Notwithstanding, the exchange of energy and momentum among the taken-into-account sources imposes some restriction on the deformation in the radial metric component, but the temporal metric component remains unchanged, which divides the EFEs. Consequently, in order to beat this weakness, the authors of [35] altered the MGD approach by incorporating deformations in both (radial/temporal) metric functions, which is valid solely in the case when the matter is absent, as well as when the Bianchi identities corresponding to self-gravitating systems filled with fluid are also non-satisfied. Later on, the author of [36] proposed a more general form of decoupling a spherical stellar structure by incorporating deformations in radial as well as temporal metric components. The major upper hand of the extended MGD approach is that it is conceivable for all space-time areas regardless of the choice of matter distribution without imposing any restrictions. In this respect, several researchers have used this extended MGD approach in the following ways: to generate anisotropic analogs of Tolman IV [37] and Krori–Barua [38] solutions, to extend black holes in $(3+1)$- and $(2+1)$-dimensional space-times [39,40,41,42], as well as to provide anisotropic solutions in GR and modified theories [27,43,44].

It is well known that various alternative or modified theories of gravity are incorporated via modifying the Einstein–Hilbert action, which is generally used to investigate the mystery of accelerated cosmic expansion as well as the truth of the dark sector i.e., dark matter and dark energy [45]. In strong-field systems, all modified gravity theories may stray broadly from GRT; however, in weak-field systems, they are in accord with GRT. It is well remarkable that the modified gravity theories in strong-field systems may prompt a correct and suitable gravity theory. The gravitational collapse phenomenon is a noticeable paradigm of a strong-field system, so its survey in modified gravity theories has pulled numerous researchers [46,47,48,49,50]. The dynamics of the collapse are modified in the case where we consider the modeling of the self-gravitational fluid distribution in generalized theories of gravity, which in turn uncovers the modification concealed in the fashioning of the Universe structure.

The Brans–Dicke (BD) theory of gravity [51,52] is the most natural suitable generalization of GRT and has been widely investigated in different astrophysical and cosmological aspects over the past decades. Their generalization is achieved by including Dirac perceptions in a scalar-tensor theory i.e., the gravitational field is interceded by the scalar field $\varphi$ as well as the curvature due to geometric part i.e., Ricci scalar. BD gravity holds a constant coupling parameter ${\omega }_{BD}$, well known as the scalar field coupling, which can be tuned to fit the perceptions [51,52,53] and also integrates a massless scalar field $\phi \left(t\right)=1/G\left(t\right)$ to analyze the involvement of the Universe. The BD coupling parameter ${\omega }_{BD}$ has been proven by a radar synchronization experiment which confirms that BD gravity is compatible with all observational tests and solar system experiments for $\left|{\omega }_{BD}\right|>4\times {10}^4$ [54,55]. On the other hand, this marvelous theory provides suitable solutions to numerous cosmic issues and also clarifies radiation-matter transition via combining observational data such as the two-degree galaxy redshift survey [56] as well as cosmic microwave background data from WMAP, VSA and CBI, which gives a broad scope of deviation from GRT [57,58,59]. We would also like to remark on recent investigations on testing BD theory and on endeavors at moderating the tensions employing specific potentials in the BD frame, see e.g., [60,61,62]. There are also tests on BD gravity with screening via scalar gravitational-wave memory, in which the author of [63] assumes that the Sun’s Vainshtein mechanism does not influence the scalar field profile generated by the distant compact celestial body and they have shown that the scalar memory effect provides a strong technique for testing BD gravity with the Vainshtein mechanism. Moreover, Sharif and his collaborators have investigated self-gravitating objects in the context of Brans–Dicke gravity [64,65,66,67].

The naturalistic systems of compact celestial bodies have widely been interpreted in scalar-tensor theories. In this regard, Ramazanoglu and Pretorius [68] investigated the allowable range taking into account the mass and scalarization of neutron stars via dispensing a mass to the scalar field. The authors of [69,70] studied how the dynamics or structure of slowly and rapidly rotating neutron stars deviate from the GRT celestial model in sight of a massive scalar field due to the most significant moment of inertia. In Ref. [71], the authors studied a certain class of scalar-tensor theories with a massive scalar field by considering a quartic self-interaction term in the potential for analyzing the comportment of static and slowly rotating neutron stars. The same authors [72] examined the impacts of a self-interacting massive scalar field on the moment of inertia and compactness of slowly rotating neutron stars by using diverse equations of state.

The coupling functions of the BD theory with a massive scalar field, which is allowed by observation, may differ significantly from those in the massless scenario, as we noted previously. Due to this occurrence, we typically come to the conclusion that anisotropic compact stellar systems with massive scalar fields could, on average, have slightly unusual topologies and properties when compared to their homologs in the massless scenario. In light of this purpose, we explore the possibility of providing exact solutions by using the decoupling of gravitational sources across MGD for viable compact stars in massive BD gravity admitting the “Karmarkar condition”.

This paper has been organized as follows: after an exhaustive introduction in Section 1, the suitable BD theory under gravitational decoupling formalism is laid down in Section 2. The basic stellar equations and the MGD approach in massive BD gravity theory are introduced in Section 3. In Section 4 and Section 5, we provide a procedure for a space-time to be embedding class-one and match the interior solution to an exterior Schwarzschild line element, respectively. The physically and mathematically adequate solutions to the modified Einstein field equations through the MGD approach with the mimic requirement on matter density and radial pressure component for anisotropy in the context of massive BD gravity are derived in Section 6. Next, in Section 7, we give the physical properties of the anisotropic stellar structure in massive BD gravity via several physical and mathematical tests. Finally, in Section 8, the conclusions are reported.

2. BD Theory under Gravitational Decoupling Formalism

In order to present the BD theory under gravitational decoupling formalism, we introduce a modified complete action in this theory that is given by

$$\begin{array}{ccc}\hfill S& =& \frac{1}{16\pi }\int \left[\mathcal{R}\phi -\frac{{\omega }_{BD}}{\phi }{\nabla }^i{\nabla }_i\phi -\mathcal{L}\left(\phi \right)\right]\sqrt{-g}{d}^{4}x\hfill \\ & +& \int {\mathcal{L}}_m\sqrt{-g}{d}^{4}x+\beta \int {\mathcal{L}}_{\theta }\sqrt{-g}{d}^{4}x.\hfill \end{array}$$

(1)

Here, $\mathcal{R}$ is the Ricci scalar while ${\mathcal{L}}_m$ and ${\mathcal{L}}_{\theta }$ denote the matter Lagrangian density and Lagrangian density for a new sector, respectively, which have been coupled by constant $\beta$, known as the coupling constant. However, $L\left(\phi \right)$ denotes a non-zero scalar field. As usual, the symbol g stands for the determinant of the metric $g_{\mu \nu }$. In this connection, the two Lagrangian density functions ${\mathcal{L}}_m$ and ${\mathcal{L}}_{\theta }$ are related to the stress-energy tensor as

$$\begin{array}{c}\hfill T_{\mu \nu }=-\frac{2}{\sqrt{-g}}\frac{\delta \left(\sqrt{-g}{\mathcal{L}}_m\right)}{\delta g^{\mu \nu }},\end{array}$$

(2)

$$\begin{array}{c}\hfill {\theta }_{\mu \nu }=-\frac{2}{\sqrt{-g}}\frac{\delta \left(\sqrt{-g}{\mathcal{L}}_{\theta }\right)}{\delta g^{\mu \nu }}.\end{array}$$

(3)

Since Lagrangian density depends solely on the metric tensor $g_{\mu \nu }$, then Equations (2) and (3) which yield the expressions for stress-energy tensors are expected to be of the form

$$\begin{array}{ccc}\hfill T_{\mu \nu }& =& g_{\mu \nu }{\mathcal{L}}_m-\frac{2\partial {\mathcal{L}}_m}{\partial g^{\mu \nu }},\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill {\theta }_{\mu \nu }& =& g_{\mu \nu }{\mathcal{L}}_{\theta }-\frac{2\partial {\mathcal{L}}_{\theta }}{\partial g^{\mu \nu }}.\hfill \end{array}$$

(5)

The modified complete action (1) for the BD gravity theory may be varied with regard to metric tensor $g^{\mu \nu }$, in order to benefit the general modified field equations as follows

$$\begin{array}{c}\hfill R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R=\frac{8\pi }{\phi }\left(T_{\mu \nu }+\beta {\theta }_{\mu \nu }\right)+\frac{1}{\phi }[{\phi }_{,\mu ;\nu }-g_{\mu \nu }\square \phi \\ \hfill +\frac{{\omega }_{BD}}{\phi }\left({\phi }_{,\mu }{\phi }_{,\nu }-\frac{g_{\mu \nu }{\phi }_{,\delta }{\phi }^{,\delta }}{2}\right)-\frac{\mathcal{L}\left(\phi \right)g_{\mu \nu }}{2}].\end{array}$$

(6)

Here, the □ represents the d’Alembert operator and it is defined as

$$\begin{array}{c}\hfill \square \phi =\frac{T^m+\Theta }{3+2{\omega }_{BD}}+\frac{1}{3+2{\omega }_{BD}}\left[\phi \frac{d\mathcal{L}\left(\phi \right)}{d\phi }-2\mathcal{L}\left(\phi \right)\right],\end{array}$$

(7)

where $T^m$ and $\Theta$ are the trace of the stress-energy tensor $T_{ij}$ and ${\theta }_{ij}$, respectively, and

$$\begin{array}{c}\hfill T^m=g^{ij}{T}_{ij},\Theta =g^{ij}{\theta }_{ij},\end{array}$$

(8)

The inner geometry of compact stellar objects is absolutely represented by a static spherically symmetric metric in Schwarzschild-like coordinates ($t,r,\theta ,\varphi$) as follows

$$\begin{array}{c}\hfill d{s}^2=e^{\xi \left(r\right)}d{t}^2-e^{\eta \left(r\right)}d{r}^2-r^2(d{\theta }^2+{sin}^2\theta d{\varphi }^2).\end{array}$$

(9)

Here, $e^{\xi \left(r\right)}$ and $e^{\eta \left(r\right)}$ represent the metric functions of an astrophysical stellar object. Moreover, all throughout this work, we use the geometrized relativistic units as {$G=c=1$} and the positive signature {$+,-,-,-$}. In this respect, we examine the physical highlights of the astrophysical stellar object with anisotropic distribution identified by the accompanying stress-energy tensor

$$T_{\mu \nu }=(\widehat{\rho }+{\widehat{p}}_t)u_{\mu }{u}_{\nu }-{\widehat{p}}_t{g}_{\mu \nu }+({\widehat{p}}_r-{\widehat{p}}_t)v_{\mu }{v}_{\nu }.$$

(10)

Here, $\widehat{\rho }$, ${\widehat{p}}_r$ and ${\widehat{p}}_t$ are the matter density, radial pressure and tangential pressure for anisotropic matter distribution, respectively, and the covariant ingredient $u_{\nu }$ is the 4-velocity with the normalization condition $u_{\mu }{u}^{\mu }=-1$ and $u_{\nu }{\nabla }^{\mu }{u}_{\mu }=0$.

3. Basic Stellar Equations and the Minimal Geometric Deformation (MGD) Approach in BD Gravity Theory

3.1. Basic Stellar Equations in BD Gravity Theory

For a physically realistic astrophysical stellar system with spherical symmetry, the combination of equations expressed in (6)–(9) and (10) provides the following differential equations:

$$\begin{array}{ccc}& & \frac{8\pi }{\phi }\left(\widehat{\rho }+\beta {\theta }_0^0\right)+\frac{1}{\phi }{f}_0^{0\left(\phi \right)}={\mathrm{e}}^{-\eta }\left(\frac{{\eta }^{\prime }}{r}-\frac{1}{r^2}\right)+\frac{1}{r^2},\hfill \end{array}$$

(11)

$$\begin{array}{ccc}& & \frac{8\pi }{\phi }\left({\widehat{p}}_r-\beta {\theta }_1^1\right)-\frac{1}{\phi }{f}_1^{1\left(\phi \right)}={\mathrm{e}}^{-\eta }\left(\frac{{\xi }^{\prime }}{r}+\frac{1}{r^2}\right)-\frac{1}{r^2},\hfill \end{array}$$

(12)

$$\begin{array}{ccc}& & \frac{8\pi }{\phi }\left({\widehat{p}}_t-\beta {\theta }_2^2\right)-\frac{1}{\phi }{f}_2^{2\left(\phi \right)}=\frac{{\mathrm{e}}^{-\eta }}{2}({\xi }^{″}+\frac{{{\xi }^{\prime }}^2}{2}+\frac{{\xi }^{\prime }-{\eta }^{\prime }}{r}\hfill \\ & & -\frac{{\xi }^{\prime }{\eta }^{\prime }}{2}),\hfill \end{array}$$

(13)

where prime indicates the derivatives with regard to the radial coordinate, r. Moreover, the scalar tensor constituents $f_0^{0\left(\phi \right)}$, $f_1^{1\left(\phi \right)}$ and $f_2^{2\left(\phi \right)}$ can be expressed in terms of $\xi$ and $\eta$, which are explicitly given by

$$\begin{array}{ccc}& & f_0^{0\left(\phi \right)}=e^{-\eta }\left[{\phi }^{″}+\left(\frac{2}{r}-\frac{{\eta }^{\prime }}{2}\right){\phi }^{\prime }+\frac{{\omega }_{BD}}{2\phi }{{\phi }^{\prime }}^2-e^{\eta }\frac{\mathcal{L}\left(\phi \right)}{2}\right],\hfill \end{array}$$

(14)

$$\begin{array}{ccc}& & f_1^{1\left(\phi \right)}=e^{-\eta }\left[\left(\frac{2}{r}+\frac{{\xi }^{\prime }}{2}\right){\phi }^{\prime }-\frac{{\omega }_{BD}}{2\phi }{{\phi }^{\prime }}^2-e^{\eta }\frac{\mathcal{L}\left(\phi \right)}{2}\right],\hfill \end{array}$$

(15)

$$\begin{array}{ccc}& & f_2^{2\left(\phi \right)}=e^{-\eta }\left[{\phi }^{″}+\left(\frac{1}{r}-\frac{{\eta }^{\prime }}{2}+\frac{{\xi }^{\prime }}{2}\right){\phi }^{\prime }+\frac{{\omega }_{BD}}{2\phi }{{\phi }^{\prime }}^2-e^{\eta }\frac{\mathcal{L}\left(\phi \right)}{2}\right].\hfill \end{array}$$

(16)

Nevertheless, from the expressions (6) and (9), we acquire,

$$\begin{array}{ccc}\hfill \square \phi & =& -e^{-\eta }\left[\left(\frac{2}{r}-\frac{{\eta }^{\prime }}{2}+\frac{{\xi }^{\prime }}{2}\right){\phi }^{\prime }\left(r\right)+{\phi }^{″}\left(r\right)\right],\hfill \end{array}$$

(17)

Taking into account Bianchi Identity, the system of equations according to BD gravity theory fulfill the accompanying preservation equation

$$\begin{array}{ccc}& & \left[-\frac{{\xi }^{\prime }}{2}(\widehat{\rho }+{\widehat{p}}_r)-{\widehat{p}}_r^{\prime }+\frac{2}{r}({\widehat{p}}_t-{\widehat{p}}_r)\right]+[-\frac{\beta {\xi }^{\prime }}{2}({\theta }_0^0-{\theta }_1^1)\hfill \\ & & +\beta {\left({\theta }_1^1\right)}^{\prime }+\frac{2\beta }{r}({\theta }_1^1-{\theta }_2^2)]+\frac{1}{8\pi }[-\frac{{\xi }^{\prime }}{2}(f_0^{0\left(\phi \right)}-f_1^{1\left(\phi \right)})\hfill \\ & & +{\left(f_1^{1\left(\phi \right)}\right)}^{\prime }+\frac{2}{r}\left(f_1^{1\left(\phi \right)}-f_2^{2\left(\phi \right)}\right)]=0.\hfill \end{array}$$

(18)

Furthermore, the mass function for the interior stellar structure can be formulated as

$$\begin{array}{c}\hfill m\left(r\right)=\frac{r(1-e^{-\eta })}{2}=4\pi \int \left(\widehat{\rho }+\beta {\theta }_0^0+\frac{f_0^{0\left(\phi \right)}}{8\pi }\right)\frac{r^2}{\phi }dr.\end{array}$$

(19)

3.2. Gravitational Decoupling by the Minimal-Geometric-Deformation (MGD) Approach

The system of equations given in (11)–(13) are a profoundly non-linear ensemble of equations. Along these lines, we shall apply the MGD method to deal with the non-linear system of Equations (11)–(13) divided into two subsystems of equations. In light of this non-linear system, let us take into account the grain anisotropic fluid solution $\{\Phi ,\Psi ,\rho ,p_r,p_t\}$ by fixing the coupling constant $\beta$ to zero, in order that the canonical metric (9) be appointed as,

$${ds}^2=e^{\Phi \left(r\right)}{dt}^2-\frac{{dr}^2}{\Psi \left(r\right)}-r^2({d\theta }^2+{sin}^2\theta {d\varphi }^2).$$

(20)

Here, $\Psi \left(r\right)=1-2\widehat{m}/r$ represents the typical term of the mass function. Currently, in order to observe the $\theta$-sector effects arising from the new sector ${\theta }_{\mu \nu }$ in the non-linear system of Equations (11)–(13), we need to include the coupling parameter $\beta$ within the anisotropic fluid distribution. In this way, so as to include the effects of the gravitational source ${\theta }_{\mu \nu }$ in an anisotropic stellar system, we take into consideration the geometric deformation functions $\Phi$ and $\Psi$, which are given as follows:

$$\begin{array}{ccc}\hfill \Phi \to \xi & =& \Phi +\beta g\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill \Psi \to e^{-\eta }& =& \Psi +\beta h,\hfill \end{array}$$

(22)

where h is the deformation linked to the constituent of the radial metrics, and g is the deformation linked to the constituent of the temporal metrics. It is worth mentioning here that the deformations expressed in (21) and (22) are absolutely radial functions; this attribute guarantees the spherical symmetry of the stellar solution. The expected MGD strategy identifies with set $h=0$ or $g=0$; for this circumstance the deformation will be carried out uniquely on the radial ingredient, surviving the temporal one unchanged, viz., $g=0$. Hence, in this regard, we obtain

$$\Psi \left(r\right)\to e^{-\eta \left(r\right)}=\Psi \left(r\right)+\beta h\left(r\right),$$

(23)

which is well known as the MGD (Minimal Geometric Deformation). Now, using this MGD approach defined in Equation (23), it is clear to observe that the non-linear system of Equations (11)–(13) divides into two systems of equations. Consequently, the first equations system is identified as follows:

$$\begin{array}{ccc}\hfill \frac{1}{\phi }\left(8\pi \widehat{\rho }+{\tilde{f}}_0^{0\left(\phi \right)}\right)& =& -\frac{{\Psi }^{\prime }}{r}-\frac{\Psi }{r{}^2}+\frac{1}{r^2},\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill \frac{1}{\phi }\left(8\pi {\widehat{p}}_r-{\tilde{f}}_1^{1\left(\phi \right)}\right)& =& \Psi \left(\frac{{\xi }^{\prime }}{r}+\frac{1}{r^2}\right)-\frac{1}{r^2},\hfill \end{array}$$

(25)

$$\begin{array}{c}\frac{1}{\phi }\left(8\pi {\widehat{p}}_t-{\tilde{f}}_2^{2\left(\phi \right)}\right)= \frac{\Psi }{4}\left({2\xi }^{″}+{{\xi }^{\prime }}^2+2\frac{{\xi }^{\prime }}{r}\right)\\ \hfill +\frac{{\Psi }^{\prime }}{4}\left({\xi }^{\prime }+\frac{2}{r}\right),\end{array}$$

(26)

it is clear that these results correspond to $\beta =0$, which implies that the fluid matter distribution is purely anisotropic. Hence, from now on, we refer to the system of Equations (24)–(26) known as the standard Einstein’s system in BD gravity theory. Moreover, we can establish straightforwardly the amounts of matter density, pressures (radial and tangential), viz., $\widehat{\rho }$, ${\widehat{p}}_r$, and ${\widehat{p}}_t$, respectively, in terms of gravitational potentials $\xi$ and $\Psi$ from expressions (24)–(26). Here, ${\tilde{f}}_0^{0\left(\phi \right)}$, ${\tilde{f}}_1^{1\left(\phi \right)}$ and ${\tilde{f}}_2^{2\left(\phi \right)}$ are explicitly indicated as

$$\begin{array}{ccc}& & {\tilde{f}}_0^{0\left(\phi \right)}=\Psi \left[{\phi }^{″}+\frac{2{\phi }^{\prime }}{r}+\frac{{\omega }_{BD}}{2\phi }{{\phi }^{\prime }}^2\right]+\frac{{\Psi }^{\prime }}{2}{\phi }^{\prime }-\frac{\mathcal{L}\left(\phi \right)}{2},\hfill \end{array}$$

(27)

$$\begin{array}{ccc}& & {\tilde{f}}_1^{1\left(\phi \right)}=\Psi \left[\left(\frac{2}{r}+\frac{{\xi }^{\prime }}{2}\right){\phi }^{\prime }-\frac{{\omega }_{BD}}{2\phi }{{\phi }^{\prime }}^2\right]-\frac{\mathcal{L}\left(\phi \right)}{2},\hfill \end{array}$$

(28)

$$\begin{array}{ccc}& & {\tilde{f}}_2^{2\left(\phi \right)}=\Psi \left[{\phi }^{″}+\left(\frac{1}{r}+\frac{{\xi }^{\prime }}{2}\right){\phi }^{\prime }+\frac{{\omega }_{BD}}{2\phi }{{\phi }^{\prime }}^2\right]+\frac{{\Psi }^{\prime }}{2}{\phi }^{\prime }\hfill \\ & & -\frac{\mathcal{L}\left(\phi \right)}{2}.\hfill \end{array}$$

(29)

The BD gravity system (24)–(26) fulfill the accompanying preservation equations,

$$\begin{array}{c}\hfill \left[-\frac{{\xi }^{\prime }}{2}(\widehat{\rho }+{\widehat{p}}_r)-{\widehat{p}}_r^{\prime }+\frac{2}{r}({\widehat{p}}_t-{\widehat{p}}_r)\right]+\frac{1}{8\pi }[-\frac{{\xi }^{\prime }}{2}({\tilde{f}}_0^{0\left(\phi \right)}\\ \hfill -{\tilde{f}}_1^{1\left(\phi \right)})+{\left({\tilde{f}}_1^{1\left(\phi \right)}\right)}^{\prime }+\frac{2}{r}\left({\tilde{f}}_1^{1\left(\phi \right)}-{\tilde{f}}_2^{2\left(\phi \right)}\right)]=0.\end{array}$$

(30)

Now, in order to establish the usual mass function, we use the definition of $\Psi =1-\frac{2\widehat{m}\left(r\right)}{r}$, which leads to the following relation:

$$\begin{array}{ccc}\hfill \widehat{m}\left(r\right)& =& {\int }_0^r\frac{1}{2\phi }\left[8\pi \widehat{\rho }\left(r\right)+{\tilde{f}}_0^{0\left(\phi \right)}\right]r^{2}dr\hfill \end{array}$$

(31)

where $\tilde{m}$ represents the mass function for the seed space-time.

Now, the second equations system, dubbed as the quasi-Einstein system in BD gravity theory, is identified in the case where we activate the coupling constant $\beta$ in the system. In this connection, we obtain

$$\begin{array}{ccc}\hfill & & \frac{1}{\phi }\left(8\pi {\theta }_0^0+{\mathcal{F}}_0^{0\phi }\right)=-\frac{h{}^{\prime }}{r}-\frac{h}{r^2},\hfill \end{array}$$

(32)

$$\begin{array}{ccc}\hfill & & \frac{1}{\phi }\left(8\pi {\theta }_1^1+{\mathcal{F}}_1^{1\phi }\right)=-h\left(\frac{{\xi }^{\prime }}{r}+\frac{1}{r^2}\right),\hfill \end{array}$$

(33)

$$\begin{array}{ccc}& & \frac{1}{\phi }\left(8\pi {\theta }_2^2+{\mathcal{F}}_2^{2\phi }\right)=-\frac{h}{4}\left({2\xi }^{″}+{{\xi }^{\prime }}^2+2\frac{{\xi }^{\prime }}{r}\right)\\ & & \hfill -\frac{h^{\prime }}{4}\left({\xi }^{\prime }+\frac{2}{r}\right).\end{array}$$

(34)

where

$$\begin{array}{ccc}& & {\mathcal{F}}_0^{0\phi }=h\left[{\phi }^{″}+\frac{2{\phi }^{\prime }}{r}+\frac{{\omega }_{BD}}{2\phi }{{\phi }^{\prime }}^2\right]+\frac{h^{\prime }}{2}{\phi }^{\prime },\hfill \end{array}$$

(35)

$$\begin{array}{ccc}& & {\mathcal{F}}_1^{1\phi }=h\left[\left(\frac{2}{r}+\frac{{\xi }^{\prime }}{2}\right){\phi }^{\prime }-\frac{{\omega }_{BD}}{2\phi }{{\phi }^{\prime }}^2\right],\hfill \end{array}$$

(36)

$$\begin{array}{ccc}& & {\mathcal{F}}_2^{2\phi }=h\left[{\phi }^{″}+\left(\frac{1}{r}+\frac{{\xi }^{\prime }}{2}\right){\phi }^{\prime }+\frac{{\omega }_{BD}}{2\phi }{{\phi }^{\prime }}^2\right]+\frac{h^{\prime }}{2}{\phi }^{\prime }.\hfill \end{array}$$

(37)

Moreover, the corresponding preservation equations related to the equations system (32)–(34) give,

$$\begin{array}{c}\hfill -\frac{{\xi }^{\prime }}{2}({\theta }_0^0-{\theta }_1^1)+\frac{d{\theta }_1^1}{dr}-\frac{2}{r}({\theta }_2^2-{\theta }_1^1)+\frac{1}{8\pi }[-\frac{{\xi }^{\prime }}{2}({\mathcal{F}}_0^{0\left(\phi \right)}\\ \hfill -{\mathcal{F}}_1^{1\left(\phi \right)})+{\left({\mathcal{F}}_1^{1\left(\phi \right)}\right)}^{\prime }+\frac{2}{r}\left({\mathcal{F}}_1^{1\left(\phi \right)}-{\mathcal{F}}_2^{2\left(\phi \right)}\right)]=0.\end{array}$$

(38)

Furthermore, we notice that the linear combination of both preservation equations given in (30) and (38) with regard to coupling constant $\beta$ provides the preservation equation for the full system, which is expressed in (18).

Actually, it is significant to note that, from now onward, we will describe the full physical parameters for the matter density, radial and tangential pressures as follows:

$$\begin{array}{ccc}\hfill \rho \left(r\right)& =& \widehat{\rho }\left(r\right)+\beta {\theta }_0^0\left(r\right),\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill p_r\left(r\right)& =& {\widehat{p}}_r\left(r\right)-\beta {\theta }_1^1\left(r\right),\hfill \end{array}$$

(40)

$$\begin{array}{ccc}\hfill p_t\left(r\right)& =& {\widehat{p}}_t\left(r\right)-\beta {\theta }_2^2\left(r\right),\hfill \end{array}$$

(41)

where $\rho$, $p_r$ and $p_t$ are physical quantities already defined above. Then again, the physical amount, namely, the anisotropy parameter ($\Delta$) can be performed by the expressions addressed in (40)–(41), read as

$$\begin{array}{c}\hfill \Delta =({\widehat{p}}_t-{\widehat{p}}_r)+\beta \left({\theta }_1^1-{\theta }_2^2\right).\end{array}$$

Presently, the internal astrophysical stellar geometry for the assisting MGD system may be specified by the accompanying metric,

$$d{s}^2=e^{\xi \left(r\right)}{dt}^2-{\left[1-\frac{2m\left(r\right)}{r}\right]}^{-1}{dr}^2-r^2({d\theta }^2+{sin}^2\theta {d\varphi }^2).$$

(42)

Here, $m\left(r\right)$ characterizes the inward mass of the MGD system, which can be read as follows:

$$m\left(r\right)=\frac{r}{2}\left[\frac{2\widehat{m}\left(r\right)}{r}-\beta h\left(r\right)\right].$$

(43)

On the other hand, it is worth mentioning here that, throughout in our study, we take the following form of the scalar function field $\mathcal{L}\left(\phi \right)$ as follows:

$$\begin{array}{c}\hfill \mathcal{L}\left(\phi \right)=\frac{1}{2}{m}_{\phi }^2{\phi }^2,\end{array}$$

(44)

where $m_{\phi }$ represents the mass of the scalar field and $\phi$ is the scalar function which is defined as $\phi =\alpha e^{n\xi \left(r\right)}$ [73] with $\alpha$ and n being two constants.

As mentioned above, our objective is to find the solution using gravitational decoupling. As we know, the MGD technique allows obtaining the solution beyond the chosen gravity. Furthermore, the MGD approach divides the original system into two parts: the first system corresponds to pure Brans–Dicke gravity while the second system is due to gravitational decoupling. We have assumed that the internal structure of the first system is made of anisotropic fluid matter distribution rather than isotropic matter distribution. Therefore, we need two extra conditions to solve the first system. Due to this reason, we are interested in the embedding class-one condition (Karmarkar condition) along with the radial metric function corresponding to the Finch–Skea model to close the first system, while the second system is solved by using the mimic constraints approach. Finally, we the solution of the original system after combining the solutions of both systems, which is the generalization of the embedding class-one solution beyond the Brans–Dicke gravity. The procedure for the embedding class-one condition is discussed in next section.

4. Procedure for a Space-Time to Be Embedding Class-One

In 1925, Eisenhart [74] presented that an implanting class-one space [$n+1$]-D space $V^{[n+1]}$ can be carried out into a [$n+2$]-D pseudo-Euclidean space $E^{[n+2]}$, which can be characterized by a [$n+1$]-D space $V^{[n+1]}$ if there remains a symmetric tensor $a_{\left[mn\right]}$ which satisfies the associated Gauss–Codazzi equations as follows:

$$\begin{array}{ccc}\hfill {\mathcal{R}}_{mnpq}& =& 2e{a}_{m[p}{a}_{q]n}\hfill \\ \hfill 0& =& a_{m[n;p]}-{\Gamma }_{\left[np\right]}^q{a}_{mq}+{\Gamma }_{m[n}^q{a}_{p]q},\hfill \end{array}$$

where ${\mathcal{R}}_{mnpq}$ and $a_{mn}$ are the Riemann curvature tensor and the coefficients of the second-order differential form, respectively, as well as $e=\pm 1$. In 1921, Kasner [75] provided that Schwarzschild’s vacuum may be submerged into 6-D pseudo-Euclidean space by a progression of coordinate transformations. This implies that the Schwarzschild outside is of class two. In the same spirit, in 1975, Gupta and Goel [76] also gave a 6-D implanted Euclidean space-time, which is specified below,

$$\begin{array}{ccc}& & z_1=k{e}^{\xi /2}cosh\left(\frac{t}{k}\right),z_2=k{e}^{\xi /2}sinh\left(\frac{t}{k}\right),z_3=f\left(r\right),\hfill \\ & & z_4=rsin\theta cos\varphi ,z_5=rsin\theta sin\varphi ,z_6=rcos\theta ,\hfill \end{array}$$

which is used to transform a generalized 4-D space-time (9) into 6-D pseudo-Euclidean space, viz.,

$$\begin{array}{ccc}\hfill d{s}^2& =& {\left(d{z}_1\right)}^2+{\left(d{z}_2\right)}^2\mp {\left(d{z}_3\right)}^2-{\left(d{z}_4\right)}^2-{\left(d{z}_5\right)}^2\hfill \\ & & -{\left(d{z}_6\right)}^2,\hfill \end{array}$$

(45)

with ${\left[f^{\prime }\left(r\right)\right]}^2=\mp \left[-\frac{(1-\Psi )}{\Psi }+\frac{1}{4}{k}^2{e}^{\xi }{{\xi }^{\prime }}^2\right]$. This also necessitates that the 4-D metric (9) may be actualized in 6-D Euclidean space; notwithstanding, there subsists a possible situation where ${\left(d{z}_3\right)}^2={\left[f^{\prime }\left(r\right)\right]}^2=0$, then the 6-D Euclidean space (45) may be reduced to 5-D pseudo-Euclidean space. This is possible only if

$$[f^{\prime }{\left(r\right)]}^2=\mp \left[-\frac{(1-\Psi )}{\Psi }+\frac{1}{4}{k}^2{e}^{\xi }{{\xi }^{\prime }}^2\right]=0,$$

(46)

or

$$\frac{1}{\Psi }=1+\frac{k^2}{4}{{\xi }^{\prime }}^2{e}^{\xi }.$$

(47)

The important condition (47) is equivalent to the embedded class-one solution (recognized as the Karmarkar condition), which fulfills a necessary condition that was originally determined by Karmarkar [77] with the help of Riemann tensor components as

$${\mathcal{R}}_{1010}{\mathcal{R}}_{2323}={\mathcal{R}}_{1212}{\mathcal{R}}_{3030}+{\mathcal{R}}_{1220}{\mathcal{R}}_{1330}.$$

(48)

In 1981, Pandey and Sharma [78] caused to notice the way that the Karmarkar condition is the only significant condition required to transform into a class-one; they found the satisfactory condition as ${\mathcal{R}}_{2323}\ne 0$. Subsequently, the basic and suitable condition to be a class-one is one which satisfies simultaneously the Karmarkar and Pandey–Sharma conditions. Thus, in terms of the gravitational potential functions, the condition expressed in (48) can be written as

$$\begin{array}{c}\hfill \frac{{\Psi }^{\prime }}{\Psi (\Psi -1)}=2({\xi }^{″}+{{\xi }^{\prime }}^2)-{{\xi }^{\prime }}^2\end{array}$$

(49)

in which, during integration, we obtain the $g_{00}$ gravitational potential function as

$$e^{\xi }={\left(A+B\int \sqrt{\frac{1}{\Psi }-1}dr\right)}^2.$$

(50)

while A and B are two integration constants. One should consistently remember that there is no class-one vacuum outside, as the static Schwarzschild’s vacuum is previously a class-two solution.

5. Exterior Space-Time: Junction Conditions

Junction conditions play a crucial and important role in the study of compact stellar structures which are bounded stellar objects by joining their interior and exterior space-time at the surface $\Sigma$ in a smooth way. In this respect, we can find all the constant parameters involved in the system. The space-time corresponding to the interior solution for field Equations (11)–(13) is expressed by Equation (42). Therefore, the exterior space-time is considered to be an empty space-time i.e., the exterior Schwarzschild solution which is given explicitly as

$$\begin{array}{ccc}\hfill d{s}_{+}^2& =& -{\left(1-\frac{2\mathcal{M}}{r}\right)}^{-1}d{r}^2-r^2(d{\theta }^2-{sin}^2\theta d{\varphi }^2)\hfill \\ & & +\left(1-\frac{2\mathcal{M}}{r}\right)d{t}^2,\hfill \end{array}$$

(51)

where $\mathcal{M}$ represents the stellar mass of the star. For analyzing the physical aspects, we need to generate a connection between the inner and outer geometries of celestial objects, which should be satisfied with the smoothness and continuity for inner and outer space-times at the hypersurface $\Sigma$ ($f=r-R=0$, where R is a radius). For this purpose, the following conditions should be met at the boundary surface $\Sigma$,

$$\begin{array}{ccc}& & {\left[{d{s}^2}_{-}\right]}_{\Sigma }={\left[{d{s}^2}_{+}\right]}_{\Sigma },{\left[{K_{ij}}_{-}\right]}_{\Sigma }={\left[{K_{ij}}_{+}\right]}_{\Sigma },\hfill \end{array}$$

(52)

$$\begin{array}{ccc}& & {[\Phi {\left(r\right)}_{-}]}_{\Sigma }={[\Phi {\left(r\right)}_{+}]}_{\Sigma },{\left[{\Phi }^{\prime }{\left(r\right)}_{-}\right]}_{\Sigma }={\left[{\Phi }^{\prime }{\left(r\right)}_{+}\right]}_{\Sigma }.\hfill \end{array}$$

(53)

1 Here, $K_{ij}$ denotes the curvature, whereas the subscripts − and + symbolize the interior and exterior solutions, respectively. Now, we use the continuity of the first fundamental form i.e., ${\left[d{s}^2\right]}_{\Sigma }$=0, and we will still find for any function $F\left(r\right)$,

$$\begin{array}{c}\hfill {\left[F\right]}_{\Sigma }\equiv F(r\to R^{+})-F(r\to R^{-})\equiv F^{+}\left(R\right)-F^{-}\left(R\right).\end{array}$$

(54)

It is worthwhile to mention that this condition supplies us $g_{rr}^{-}\left(R\right)=g_{rr}^{+}\left(R\right)$ and $g_{tt}^{-}\left(R\right)=g_{tt}^{+}\left(R\right)$. Moreover, the second fundamental form, viz., $K_{ij}$, which is equivalent to the O${}^{\prime }$Brien and Synge [79] junction condition at the hypersurface $\Sigma$, must be satisfied by the space-time (7). This condition denotes that radial fluid pressure $p_r$ vanishes at a certain radius $r=r_{\Sigma }$, which is expressed as

$$\begin{array}{c}\hfill [p_r{\left(r\right)]}_{\Sigma }=0⟹{\left[{\widehat{p}}_r-\beta {\theta }_1^1\left(r\right)\right]}_{\Sigma }=0.\end{array}$$

(55)

This requirement establishes the size of the celestial object i.e., the radius R, which implies that the material substance is restricted within the interior of the sector $0\le r\le R$, which indicates that its exterior sector does not contain any matter. Moreover, the BD scalar field $\Phi$ relating to the empty Schwarzschild solution was determined according to the procedure addressed in [80]. Now, considering ${\tau }^{-}$ and ${\tau }^{+}$ as the inner and outer space, respectively, the line element for hypersurface is defined in terms of the proper time boundary $\gamma$ as

$$\begin{array}{c}\hfill d{s}^2=d{\gamma }^2-R^2(d{\theta }^2+{sin}^2\theta d{\varphi }^2).\end{array}$$

(56)

The extrinsic curvature corresponding to the boundary $\Sigma$ is given as follows:

$$\begin{array}{c}\hfill K_{ij}^{\pm }=-m_k^{\pm }\frac{{\partial }^2{y}_{\pm }^k}{\partial n^i{n}^j}-m_k^{\pm }{\Gamma }_{ql}^k\frac{\partial y_{\pm }^q}{\partial n^i}\frac{\partial y_{\pm }^l}{\partial n^j}\end{array}$$

(57)

where $n^i$ and $m_k^{\pm }$ denote the coordinates in the boundary $\Sigma$ and the 4-velocity normal to $\Sigma$, respectively. Here, the components corresponding to the 4-velocity are defined in terms of the coordinates ($y_{\pm }^i$) of ${\tau }^{\pm }$ as

$$\begin{array}{c}\hfill m_k^{\pm }=\pm \frac{df}{d{y}^k}|g^{ql}\frac{df}{d{y}^q}\frac{df}{d{y}^l}{|}^{-1/2}\mathrm{with}{m}_k{m}^k=1.\end{array}$$

(58)

Further, the unit normal vectors for the inner and outer areas can be defined as

$$\begin{array}{ccc}& & m_k^{-}=\left[0,e^{\eta /2},0,0\right],\mathrm{and}\hfill \\ & & m_k^{-}=\left[0,{\left(1-\frac{2\mathcal{M}}{r}\right)}^{-1/2},0,0\right].\hfill \end{array}$$

(59)

Now, considering the Schwarzschild space-time (51) and using the line elements from Equations (7) and (56), we can write

$$\begin{array}{c}\hfill {\left[\frac{dt}{d\gamma }\right]}_{\Sigma }={\left[e^{-\xi /2}\right]}_{\Sigma }={\left[{\left(1-\frac{2\mathcal{M}}{r}\right)}^{-1/2}\right]}_{\Sigma }.\end{array}$$

(60)

$\mathrm{where}{\left[r\right]}_{\Sigma }=R$. The non-zero components corresponding to the curvature, $K_{ij}$, can be achieved by using expression (59) as

$$\begin{array}{ccc}& & K_{00}^{-}={\left[-\frac{{\xi }^{\prime }}{2e^{\eta /2}}\right]}_{\Sigma },K_{22}^{-}=\frac{1}{{sin}^2\theta }{K}_{33}^{-}={\left[r{e}^{-\eta /2}\right]}_{\Sigma },\hfill \\ & & K_{00}^{+}={\left[\frac{\mathcal{M}}{r^2}{\left(1-\frac{2\mathcal{M}}{r}\right)}^{-1/2}\right]}_{\Sigma },\hfill \\ & & K_{22}^{+}=\frac{1}{{sin}^2\theta }{K}_{33}^{+}={\left[r{\left(1-\frac{2\mathcal{M}}{r}\right)}^{1/2}\right]}_{\Sigma }.\hfill \end{array}$$

By using the junction condition ${\left[K_{22}^{-}\right]}_{\Sigma }={\left[K_{22}^{+}\right]}_{\Sigma }$ with ${\left[r\right]}_{\Sigma }$, we obtain the component of the radial space-time as follows:

$$\begin{array}{c}\hfill e^{-\eta \left(R\right)}=\left(1-\frac{2\mathcal{M}}{R}\right).\end{array}$$

(61)

Plugging this radial gravitational component into the matching condition ${\left[K_{00}^{-}\right]}_{\Sigma }={\left[K_{00}^{+}\right]}_{\Sigma }$ provides

$$\begin{array}{c}\hfill {\xi }^{\prime }\left(R\right)=\frac{\mathcal{M}}{R^2}{\left(1-\frac{2\mathcal{M}}{R}\right)}^{-1}.\end{array}$$

(62)

On the other hand, in principle, the $\theta$-sector may create some modification in outer space-time and matter content. As such, the explicit expression for Equation (55) can be given by

$${\widehat{p}}_r\left(R\right)-\beta {\left({\theta }_1^1\right)}^{-}\left(R\right)=-\beta {\left({\theta }_1^1\right)}^{+}\left(R\right),$$

(63)

where ${\widehat{p}}_r\left(R\right)$ refers to the radial pressure with respect to the inside matter distribution when $\beta =0$. The above condition (63) describes a second fundamental form related to the Einstein field equations expressed explicitly by Equation (6). Now, the formula expressed in (63) with the help of the value of ${\left({\theta }_1^1\right)}^{-}\left(R\right)$ for the interior space-time via Equation (33) comes out to be

$$\begin{array}{ccc}& & \frac{\beta h}{8\pi }{\left[\phi \left(\frac{{\xi }^{\prime }}{r}+\frac{1}{r^2}\right)+\left(\frac{2}{r}+\frac{{\xi }^{\prime }}{2}\right){\phi }^{\prime }-\frac{{\omega }_{BD}}{2\phi }{{\phi }^{\prime }}^2\right]}_{r=R}\hfill \\ & & +{\widehat{p}}_r\left(R\right)=-\beta {\left({\theta }_1^1\right)}^{+}\left(R\right),\hfill \end{array}$$

(64)

where ${\xi }^{\prime }\equiv {\partial }_r{\xi }^{-}$ and $\Sigma$ represent the boundary of the celestial configuration $r=R$. Presently, we employ Equations (33), (61) and (62) into Equation (64) in order to find the ${\left({\theta }_1^1\right)}^{+}\left(R\right)$ expressed in Equation (64) for the exterior space-time, leading to

$$\begin{array}{ccc}& & {\widehat{p}}_r\left(R\right)+\frac{\beta h}{8\pi }{\left[\phi \left(\frac{{\xi }^{\prime }}{r}+\frac{1}{r^2}\right)+\left(\frac{2}{r}+\frac{{\xi }^{\prime }}{2}\right){\phi }^{\prime }-\frac{{\omega }_{BD}}{2\phi }{{\phi }^{\prime }}^2\right]}_{r=R}\hfill \\ & & =\frac{\beta h_{{}_{\Sigma }}^{*}}{8\pi }[{\phi }_{{}_{\Sigma }}\left(\frac{2M}{R^2\left(R-2M\right)}+\frac{1}{R^2}\right)-\frac{{\omega }_{BD}}{2{\phi }_{{}_{\Sigma }}}{\phi }_{{}_{\Sigma }}^{\prime 2}\hfill \\ & & +\left(\frac{2}{R}+\frac{M}{R\left(R-2M\right)}\right){\phi }_{{}_{\Sigma }}^{\prime }].\hfill \end{array}$$

(65)

Here, $h_{{}_{\Sigma }}^{*}$ and ${\phi }_{{}_{\Sigma }}$ are two functions signifying the geometric deformation function and scalar function, respectively, for the outside Schwarzschild line element (51) at the boundary $r=R$ of the celestial astrophysical structure due to the new sector ${\theta }_{ij}$. It should be noted that, if the outer space-time is represented by the Schwarzschild vacuum solution given in (51), at this stage, we should put $h_{{}_{\Sigma }}^{*}=$ 0, in order to obtain the final form of Equation (65), which can be derived as follows:

$$\begin{array}{c}\hfill {\widehat{p}}_r\left(R\right)+\frac{\beta h\left(R\right)}{8\pi }[\phi \left(R\right)\left(\frac{{\xi }^{\prime }\left(R\right)}{R}+\frac{1}{R^2}\right)-\frac{{\omega }_{BD}}{2\phi \left(R\right)}{{\phi }^{\prime }}^2\left(R\right)\\ \hfill +\left(\frac{2}{R}+\frac{{\xi }^{\prime }\left(R\right)}{2}\right){\phi }^{\prime }\left(R\right)]=0,\end{array}$$

(66)

which is identical to

$$\begin{array}{c}\hfill {\widehat{p}}_r\left(R\right)-\beta {\theta }_1^1\left(R\right)=0⟹p_r\left(R\right)=0.\end{array}$$

(67)

This condition (66) or (67) establishes the size of the compact celestial body. In view of this relevant condition, we have discovered an important result that anisotropic radial pressure vanishes at the surface, which is desired for the equilibrium of a celestial object in the original Schwarzschild solution.

Finally, we want to mention that conditions (60)–(62) and (67) describe the primary and adequate conditions to establish all the constants present in the stellar system.

6. Minimally Deformed Anisotropic Solution: Mimic Constraint Approach

As the modified Einstein field equations given in (11)–(13) are divided into two subsystems according to the MGD method, the first subsystem corresponds to the field equations given in (24)–(26) with five unknown parameters $\{\Psi ,\xi ,\rho ,p_r,p_t\}$ for the seed solution in BD gravity with respect to the stress-energy tensor $T_{\mu \nu }^m$, whereas the second subsystem corresponds to the field equations expressed in (32)–(34) which examine the solution of new $\theta$-sector, which basically consists of four unknown parameters $\{h\left(r\right),{\theta }_0^0,{\theta }_1^1,{\theta }_2^2\}$. In this regard, if the first subsystem generally provides an achievable solution for anisotropic distribution, then this leads directly to the absolute evaluation of the unknown parameters of the second subsystem. In effect, for solving the system of equations via the MGD approach, we need to use the well-known embedding class-one condition to reduce the degree of freedom. Now, we take into account the ansatz for $\Psi$ related to a well-defined notable Finch–Skea model and execute a well-known embedding class-one condition (50). Hence, we obtain the accompanying seed space-time, so-called embedding class-one seed space-time, as

$$\begin{array}{ccc}\hfill d{s}^2& =& {\left(A+\frac{\sqrt{a}B}{2}{r}^2\right)}^{2}d{t}^2-\frac{d{r}^2}{(1+a{r}^2)}\hfill \\ & & -r^2(d{\theta }^2+{sin}^2\theta d{\varphi }^2).\hfill \end{array}$$

(68)

At this stage, $\Psi =1/(1+a{r}^2)$ and $\xi =ln{\left(A+\frac{\sqrt{a}B}{2}{r}^2\right)}^2$ can describe the solution of the equations system (24)–(26) in the context of the BD gravity theory. Presently, in order to close the second subsystem with the new $\theta$-sector, we necessitate the deformation function indicated by $h\left(r\right)$. Since there are diverse ways to obtain the deformation function, however, it ought to be satisfied with some fundamental prerequisites so as to be a permissible system with respect to the physical and mathematical points of view. In such a manner, a straightforward and comfortable depiction of the geometric deformation function $f\left(r\right)$ [28,81,82] ought to be satisfied with these fundamental necessities or only involve oneself in the ingredients of $\theta$-sector with a linear, polytropic or barotropic equation of state. Here, we will adopt the well-known methodology, specifically, the Mimic requirement method:

(A)

Mimic requirement on the radial pressure component of the anisotropic sector, viz., (${\theta }_1^1=p_r$),

(B)

Mimic requirement on the density of the anisotropic sector, viz., (${\theta }_0^0=\rho$).

It is mentioned that both approaches are mainly used to determine the deformation function (h). The mimic-to-pressure constraint approach i.e., ${\theta }_1^1=p_r$ involves a deformation function (h) that is free from the derivative. This approach directly determines the expression for $h\left(r\right)$ without any integration, while the mimicking of density constraint i.e., ${\theta }_1^1=p_r$ procedure provides a differential in the differential equation in the deformation function (h), whose solution depends on the seed solution. The second approach enhances the total mass of the object.

6.1. Mimic Requirement on Radial Pressure Component for Anisotropy

Here, in order to obtain a physically and mathematically adequate solution, we perform the radial pressure component expressed in Equation (25) and imitate its relationship with the anisotropic area given by Equation (33) in order to close system of Equations (32)–(34). In this scenario, we can obtain

$${\theta }_1^1\left(r\right)={\widehat{p}}_r\left(r\right).$$

(69)

This necessity derives that the stress-energy tensor for the seed solution is associated with the anisotropy in a radial direction. This implies that Equations (25) and (33) are equivalent. In this way, we determine the formula for the deformation function $h\left(r\right)$, which takes the following forms:

$$\begin{array}{c}\hfill h\left(r\right)=\frac{\Psi {{\phi }^{\prime }}^2{r}^2{\omega }_{BD}-\Psi r\phi {\phi }^{\prime }(4+{\xi }^{\prime }r)+f_{11}\left(r\right)}{2{\phi }^2(1+{\xi }^{\prime }r)+\phi {\phi }^{\prime }r(4+{\xi }^{\prime }r)-{{\phi }^{\prime }}^2{r}^2{\omega }_{BD}}.\end{array}$$

(70)

where, $f_{11}\left(r\right)=\phi r^2\mathcal{L}\left(\phi \right)-2{\phi }^2(\Psi -1+\Psi {\xi }^{\prime }r)$.

Then, the minimally deformed radial gravitational potential $e^{-\eta }$ has taken its final form, which is expressed explicitly as

$$e^{-\eta }=\Psi +\beta \frac{\Psi {{\phi }^{\prime }}^2{r}^2{\omega }_{BD}-\Psi r\phi {\phi }^{\prime }(4+{\xi }^{\prime }r)+f_{11}\left(r\right)}{2{\phi }^2(1+{\xi }^{\prime }r)+\phi {\phi }^{\prime }r(4+{\xi }^{\prime }r)-{{\phi }^{\prime }}^2{r}^2{\omega }_{BD}}.$$

(71)

Now by combining Equations (44) and (68) into Equations (70) and (71), we obtain

$$\begin{array}{ccc}\hfill h\left(r\right)& =& -\frac{r^2[{\mathcal{H}}_1\left(r\right)+{\mathcal{H}}_2\left(r\right)]}{4^{n+1}(1+a{r}^2)[-4A^2-{\mathcal{H}}_3\left(r\right)]},\hfill \end{array}$$

(72)

$$\begin{array}{ccc}\hfill e^{-\eta }& =& \frac{-r^2[{\mathcal{H}}_1\left(r\right)+{\mathcal{H}}_2\left(r\right)]}{4^{n+1}(1+a{r}^2)[{\mathcal{H}}_3\left(r\right)-4A^2]}+\frac{1}{1+a{r}^2},\hfill \end{array}$$

(73)

where,

$$\begin{array}{ccc}& & {\mathcal{H}}_1\left(r\right)=4A^2\alpha m_{\phi }^2{(2A+\sqrt{a}B{r}^2)}^{2n}-4\sqrt{a}AB[2^{3+2n}\hfill \\ & & +4^{2+n}n-\alpha m_{\phi }^2{r}^2{(2A+\sqrt{a}B{r}^2)}^{2n}]+4a^{3/2}AB{r}^2\hfill \\ & & \times [4^{1+n}+\alpha m_{\phi }^2{r}^2{(2A+\sqrt{a}B{r}^2)}^{2n}],\hfill \\ & & {\mathcal{H}}_2\left(r\right)=a^2{B}^2{r}^4[4^{1+n}+\alpha m_{\phi }^2{r}^2{(2A+\sqrt{a}B{r}^2)}^{2n}]+\hfill \\ & & a\{4A^2[4^{1+n}+\alpha m_{\phi }^2{r}^2{(2A+\sqrt{a}B{r}^2)}^{2n}]+B^2{r}^2\hfill \\ & & \times [-4^{2+n}-4^{3+n}n+\alpha m_{\phi }^2{r}^2{(2A+\sqrt{a}B{r}^2)}^{2n}\hfill \\ & & +2^{5+2n}{n}^2{\omega }_{BD}\left]\right\}.\hfill \\ & & {\mathcal{H}}_3\left(r\right)=4\sqrt{a}AB(3+4n)r^2+a{B}^2{r}^4(5+16n-8n^2{\omega }_{BD}).\hfill \end{array}$$

Now we can determine the expressions for $\theta$-sector components directly from Equations (32)–(34), namely ${\theta }_0^0$, ${\theta }_1^1$ and ${\theta }_2^2$ after plugging the deformation function $h\left(r\right)$. In this situation, we can obtain

$$\begin{array}{ccc}\hfill p_r^{\left(\mathrm{eff}\right)}& =& (1-\beta ){\widehat{p}}_r\left(r\right)\hfill \end{array}$$

(74)

$$\begin{array}{ccc}\hfill p_t^{\left(\mathrm{eff}\right)}& =& p_r^{\left(\mathrm{eff}\right)}+{\Delta }^{\mathrm{eff}},\hfill \end{array}$$

(75)

However, Equation (39) provides the explicit expression for effective matter density, viz., ${\rho }^{\left(\mathrm{eff}\right)}$ as,

$${\rho }^{\left(\mathrm{eff}\right)}=\widehat{\rho }\left(r\right)+\beta {\theta }_0^0.$$

(76)

This mimic constraint procedure suggests that the deformed gravitational potential function $e^{\eta }$ can be determined explicitly, due to the fact that the condition expressed in Equation (71) does not contain any integral term. On the other hand, the energy density $\widehat{\rho }\left(r\right)$, radial ${\widehat{p}}_r\left(r\right)$ and tangential ${\widehat{p}}_t\left(r\right)$ pressures for the seed solution are already given by Equations (24)–(26), respectively. Now, after applying the boundary conditions (60)–(62), we obtain

$$\begin{array}{ccc}& & A+\frac{\sqrt{a}B}{2}{R}^2=\sqrt{1-\frac{2\mathcal{M}}{R}},\hfill \end{array}$$

(77)

$$\begin{array}{ccc}& & \sqrt{a}BR=\frac{2\mathcal{M}}{R^2},\hfill \end{array}$$

(78)

$$\begin{array}{ccc}& & \frac{1}{1+a{R}^2}+\beta h\left(R\right)=1-\frac{2\mathcal{M}}{R}.\hfill \end{array}$$

(79)

The constant $\mathcal{M}$, which is the mass of the stellar system, will be determined numerically by using condition (67).

6.2. Mimic Requirement on Matter Density for Anisotropy

In this section, we consider another approach, namely, a mimic requirement on matter density for anisotropy, which is also an effective technique to determine the physically acceptable form of deformation function $f\left(r\right)$. In this approach, we mimic the seed matter density $\widehat{\rho }\left(r\right)$ to its temporal component ${\theta }_0^0\left(r\right)$, which is as follows

$${\theta }_0^0\left(r\right)=\widehat{\rho }\left(r\right).$$

(80)

Moreover, we acquire the first order linear differential equation in the deformation function $h\left(r\right)$ (by taking into account Equations (24)–(32) with the help of Equation (80)), which leads to

$$\begin{array}{c}\hfill \frac{dh}{dr}+F(\phi ;\Psi )h=G(\phi ;\Psi ),\end{array}$$

(81)

which provides

$$\begin{array}{c}\hfill h=e^{-\int F(\phi ;\Psi )dr}\left(\int G(\phi ;\Psi )e^{\int F(\phi ;\Psi )dr}dr+H_1\right).\end{array}$$

(82)

Here, $H_1$ is an integration constant while $F(\phi ;\Psi )$ and $G(\phi ;\Psi )$ take the following form:

$$\begin{array}{ccc}& & F(\phi ;\Psi )=\frac{{{\phi }^{\prime }}^2{r}^2{\omega }_{BD}+4r\phi {\phi }^{\prime }+2{\phi }^{″}{r}^2\phi +2{\phi }^2}{r\phi ({\phi }^{\prime }r+2\phi )},\hfill \\ & & G(\phi ;\Psi )=\frac{1}{r\phi ({\phi }^{\prime }r+2\phi )}[{{\phi }^{\prime }}^2{r}^2{\omega }_{BD}\Psi \phi +{\Psi }^{\prime }{\phi }^{\prime }{r}^2{\phi }^2+2{\phi }^3\hfill \\ & & -\mathcal{L}\left(\phi \right)r^2{\phi }^2+4\phi r\Psi {\phi }^2+2{\phi }^{″}{r}^2\Psi {\phi }^2+2\Psi r{\phi }^3+2\Psi {\phi }^3].\hfill \end{array}$$

It is worthwhile to mention here that the deformed gravitational potential $e^{-\eta }$ related to radial metric component $\Psi$ can provide via Equation (82)

$$\begin{array}{ccc}\hfill e^{-\eta }& =& \beta e^{-\int F(\phi ;\Psi )dr}\left(\int G(\phi ;\Psi )e^{\int F(\phi ;\Psi )dr}dr+H_1\right)\hfill \\ & & +\Psi .\hfill \end{array}$$

(83)

Next we can find

$$\begin{array}{ccc}\hfill p_r^{\left(\mathrm{eff}\right)}& =& {\widehat{p}}_r\left(r\right)-\beta {\theta }_1^1\hfill \end{array}$$

(84)

$$\begin{array}{ccc}\hfill p_t^{\left(\mathrm{eff}\right)}& =& {\widehat{p}}_t\left(r\right)-\beta {\theta }_2^2,\hfill \end{array}$$

(85)

also using the formula (${\rho }^{\left(\mathrm{eff}\right)}=\widehat{\rho }+\beta {\theta }_0^0$); with the assistance of condition (80), we obtain the expression of ${\rho }^{\left(\mathrm{eff}\right)}$ as

$${\rho }^{\left(\mathrm{eff}\right)}=\left(1+\beta \right)\widehat{\rho }\left(r\right).$$

(86)

Now we would like to highlight that, as we can see from Equation (83), the gravitational potential $\eta$ comprehends integral terms; in this way, it is not always possible to obtain deformed radial gravitational potential $e^{\eta }$ for a given function of $\xi$. Subsequently, the solution relating to this mimic requirement can be obtained by a numerical methodology that will be considered in future works. Therefore, in the next section, we focus and discuss the physical viability of the solutions under massive BD gravity obtained by the first approach (Section 6.1), namely: mimic requirement on radial pressure component for anisotropy.

7. Physical Properties of the Anisotropic Stellar Structure in the Massive BD Gravity Theory

In the present section, we examine the physical validity of the accomplished solutions under massive BD gravity by attempting to investigate various physical properties of the anisotropic astrophysical structure via the MGD approach in the accompanying subsections.

7.1. Thermodynamic Observable

In this subsection, we introduce the main physical amounts that identify our stellar system. These amounts are the matter density $\rho$, the radial $p_r$ and transverse $p_t$ pressures and the anisotropy parameter $\Delta$. It is worth mentioning that, in the investigation of anisotropic compact stellar configurations, all the amounts indicated above should fulfill some overall prerequisite in order to represent a well-defined astrophysical stellar interior. These overall requirements are as follows:

• Non-physical and geometrical singularities in all interior points of the stellar structure exist.
• All the physical quantities, viz., the matter density $\rho$, radial pressure $p_r$ and transverse pressure $p_t$ should have their most extreme values at the core of the stellar structure, which involves that they are monotonously decreasing functions with increasing radius towards the frontier of the spherical object.
• Concerning the anisotropy parameter $\Delta$, it should be zero at the core i.e., $\Delta =0$. This is because, at the core of the stellar structure, the radial pressure is equal to the transverse pressure ($p_r=p_t$). Then again, its values must be positive towards the boundary of the spherical object (in the case of repulsive anisotropic force).

The complete comportment of these thermodynamic observable components, viz., $\rho$, $p_r$, $p_t$ and $\Delta$ is illustrated in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11. In these graphs, we can see that all the physical amounts that characterize the stellar system satisfy the overall necessities. A few remarks concerning the physical and mathematical point of view of these quantities are given below:

• Figure 1, Figure 2 and Figure 3 show the variation of the matter density with respect to the radial coordinate for all chosen values of the parameters $\alpha$, $\beta$ and ${\omega }_{BD}$, respectively, by setting $M=1.97M_{\odot }$, $R=9.69$ km, $m_{\varphi }=0.001$, $n=-0.01$, $B=(0.0416,0.0415)/$km, $A=(0.38284,0.3857)$ and $a=(0.0163,0.016)/$km${}^2$. We can see from these plots that the matter density has its most extreme values at the center of the stellar structure and decreases monotonically towards the surface with an increasing radius, r. Additionally, it is positive everywhere inside the stellar configuration. Interestingly, the effect of the three parameters, $\alpha$, $\beta$ and ${\omega }_{BD}$, on the density energy, $\rho$, has been further demonstrated for both BD and BD+MGD scenarios, where any increase in $\alpha$ and $\beta$ shifts the energy density, $\rho$, to higher equilibrium values throughout the stellar interior for all $\alpha \in [0.3,0.9]$ and $\beta \in [0.0,1.0]$, while any increase in ${\omega }_{BD}$ shifts the energy density to lower equilibrium values throughout the stellar interior for all ${\omega }_{BD}\in [0.0,1000]$.
• Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 display the variation of the radial and transverse pressures for the same parameters used in Figure 1, Figure 2 and Figure 3. From the five graphs, we can see that these physical amounts have their most extreme values at the core of the stellar object and are monotonic decreasing functions with increasing radial coordinates, r. Moreover, the radial pressure $p_r$ disappears at the surface of the stellar object, and the transverse $p_t$ one dominates at all points. It should be noted that the equilibrium values for $p_r$ and $p_t$ increase when $\alpha \in [0.3,0.9]$ increases from $0.3$ to $0.9$, as illustrated in Figure 4 and Figure 7. In contrast, when we increase $\beta \in [0.0,1.0]$ from $0.0$ to $1.0$, the radially and tangentially outward pressures are reduced for both the BD and BD+MGD scenarios (see Figure 5). An emerging characteristic in Figure 6 and Figure 8 is that the parameter ${\omega }_{BD}$ reduces the radially outward pressure and enhances the tangential surface stresses, as the parameter ${\omega }_{BD}\in [0.0,1000]$ increases.
• Figure 9, Figure 10 and Figure 11 exhibit the variation of anisotropy parameters for all the parameters taken into account in Figure 1, Figure 2 and Figure 3, whatever the interior of the stellar structure. Its comportment is positive at all points inside the stellar structure, disappearing at the core i.e., $\Delta =0$, in impact at the core of the star $p_r=p_t$ and has a monotonic increasing function with increasing radial coordinate, r. As we clarified above, $\Delta >0$, which implies that $p_t>p_r$; this presents that the system encounters a repulsive force (outwards). This last item counteracts the gravitational gradient progressing the stability and equilibrium state as well as one having more compact and massive stellar configurations. It is interesting evidence that the three parameters, $\alpha$, $\beta$ and ${\omega }_{BD}$, for both the BD and BD+MGD scenarios enhance the turbulence throughout the stellar interior, as $\Delta$ increases for $\alpha \in [0.3,0.9]$, $\beta \in [0.0,1.0]$ and ${\omega }_{BD}\in [0.0,1000]$, respectively (see Figure 9, Figure 10 and Figure 11 for more details).

7.2. Causality and Stability Condition

All through the adjustments introduced by the anisotropies source ${\theta }_{\mu ,\nu }$ in the stellar interior, we need to examine if the irritation presented by this new source provides to the stellar system a stable equilibrium. To do this, we employ Abreu’s criteria dependent on Herrera’s cracking approach [83]. Fundamentally, this criterion establishes the stability of compact stellar structures via the speeds of propagation related to the pressure waves in the principal directions of the spherical object i.e., in radial and transverse directions. Consequently, according to this criterion, the subliminal velocities corresponding to every direction are characterized by

$$v_r^2=\frac{d{p}_r}{d\rho }\mathrm{and}{v}_t^2=\frac{d{p}_t}{d\rho }.$$

(87)

In order to obtain a genuinely allowable stellar system, the two velocities $v_r$ and $v_t$ should be limited by the velocity of light ($c=1$ in relativistic geometrized units). This reveals to us that the sound waves in the fluid do not spread at arbitrary velocities. This is known as the causality condition, which can be read as follows:

$$0\le v_r\le 1\mathrm{and}0\le v_t\le 1.$$

(88)

This causality condition is authoritative notwithstanding if the material substance of the stellar structure is isotropic or anisotropic and also has strong ramifications on the comportment of the matter distribution inside the stellar configuration. In this respect, from Figure 12, Figure 13 and Figure 14, we can observe that $v_r^2$ and $v_r^2$ fulfill the causality condition for all chosen values of the parameters $\alpha$, $\beta$ and ${\omega }_{BD}$, leading to a stable and viable stellar configuration under massive BD gravity via the MGD approach. It is intriguing to notice that $v_t$ is sensitive to changes in parameter strength, ${\omega }_{BD}$, in comparison to $v_r$. The radial variation for both remains causal: $v_r\in [0.0,1.0]$ and $v_t\in [0.0,1.0]$ throughout the stellar interior.

Consequently, according to these three graphs, we can clearly see that our stellar system is fully stable under Abreu’s criteria based on Herrera’s cracking approach.

7.3. $M-R$ and $M-I$ Diagrams with Observational Data

We analyze the $M-R$ and $M-I$ curves generated from the solution under massive BD gravity and give a pedagogical clarification of the effects presented by gravitational decoupling by means of MGD with the choices made on parameters $\alpha$ and ${\omega }_{BD}$, in order to provide a more realistic scenario and reliable stellar system. In this connection, we bound the mentioned parameters by employing genuine observational data.

From Figure 15 and Figure 16, we give the total mass in Solar mass, i.e., M in $M_{\odot }$ with the total radius R for the chosen specific values of parameters $\alpha$ and ${\omega }_{BD}$. In this respect, we have examined two cases. In the first case, we have made $\alpha$ vary according to values $0.3$, $0.5$, $0.7$, $0.9$, by fixing the BD-parameter ${\omega }_{BD}$ to 5 as shown in Figure 15. In the second case, we varied ${\omega }_{BD}$ as 0, 500, 1000, by fixing the parameter $\alpha$ to $0.5$ as illustrated in Figure 16. In both cases, we also define the parameters: $m_{\varphi }=0.001$, $n=-0.01$, $\beta =0.2$, $B=(0.0415,0.0416/$km, $A=(0.3857,0.38284)$ and $a=(0.016,0.0163)/$km${}^2$. We sum up the new insights as follows:

• Regarding the first case, Figure 15 highlights that, as the values of $\alpha$ increase, the maximal mass M and the corresponding radial coordinate R increase accordingly. Hence, this increase in the maximum mass on $M-R$ curves confirms the existence of gravitational decoupling (i.e., non-disappearing anisotropic term), which is expected. The most extreme mass in the $M-R$ curves is about $2.75M_{\odot }$, and in this manner, more massive compact stellar structures can be fitted. It also provides a direct correlation between the parameter, $\alpha$, and the maximum mass, M, and its radius, R, where any increase in $\alpha$ shifts the stable compact stellar structure to a lower M at a lower confining radius for each $\alpha \in [0.3,0.9]$. For instance, the highest maximum mass, M, is around $2.75$$M_{\odot }$, and its radius, R, is $9.35$ km with $\alpha =0.9$, while the lowest maximum mass, M, is $1.56$$M_{\odot }$, and its radius, R, is $5.40$ km with $\alpha =0.3$.
• Regarding the second case, Figure 16 shows that the most extreme value of mass progressively decrements for the incrementing values of BD-parameter ${\omega }_{BD}$. In this regard, our anisotropic astrophysical stellar system becomes more massive and transforms into more compact stellar structures. We also demonstrate the impact of the parameter ${\omega }_{BD}$ on the mass and its radius by considering values of ${\omega }_{BD}$, such that ${\omega }_{BD}\in [0.0,1000]$. We find that the maximum mass and its radius decrease as ${\omega }_{BD}$ increases, for example, where the highest maximum mass, M, is $2.5512$$M_{\odot }$ and its radius, R, is $8.76$ km for ${\omega }_{BD}=0$; the lowest maximum mass, M, is $2.484$$M_{\odot }$ and its radius, R, is greater than $8.62$ km for ${\omega }_{BD}=1000$.

Finally, for both cases, we have discovered a good adjustment for six compact stellar structures, viz., LMC X-4 (red), LMC X-4 (black), Her X-1 (blue), 4U 1820-30 (Orange), 4U 1608-52 (Cyan), SAX J1808.4-3658 (Purple) and many other compact stars that can be adjusted. On the other hand, it is clear how the two parameters $\alpha$ and ${\omega }_{BD}$ presented by massive BD gravity incorporating the anisotropic comportment of the matter distribution has a great influence on many physical parameters of the stellar structure.

On the other hand, adopting the Bejger and Haensel concept [84], one can examine the moment of inertia I connecting to a static stellar solution, which could provide a precise tool for determining the stiffness of an equation of state (EoS). This concept is given by

$$I=\frac{2}{5}\left(1+\frac{(M/R)·km}{M_{\odot }}\right)M{R}^2.$$

(89)

In Figure 17 and Figure 18, we feature the variation of the maximum moment of inertia I with respect to the total mass M for the same values of parameters taken into consideration in $M-R$ curves, respectively. It is facile to see the effect of $\alpha$ and ${\omega }_{BD}$ on $I-M$ curves from both figures for all the increasing value choices of the parameters $\alpha$ and ${\omega }_{BD}$. In this regard, we observe that the maximum moment of inertia I is always increasing from the null central value up to a specific mass value, and afterward decreases quickly. Therefore, we can conclude that the stiffness of the EoS is better through increasing the sensitivity of $I-M$ curves in both cases corresponding to the increase of the parameters $\alpha$ and ${\omega }_{BD}$. Moreover, our resulting $I-M$ curves are also equipped with observational data for the same well-known compact stellar structures fixed in $M-R$ curves, which are well fitted.

7.4. Mass–Central Density Relationship

To provide more investigations of the stability of the compact stellar system, we focus on the mass–density relation known as the static stability criterion, which is a significant thermodynamic amount. This static stability criterion was pointed out initially by Chandrashekhar [85] to describe the stability of a gaseous stellar structure according to radial perturbations. On the other hand, the development and simplification of this static stability criterion was envisaged by Harrison et al. [86] and Zeldovich and Novikov [87]. This Harrison–Zeldovich–Novikov static stability criterion infers that any solution describes static and stable stellar structure if the gravitational mass M is an increasing function of its central density ${\rho }_c$ i.e.,

$$\frac{\partial M\left({\rho }_c\right)}{\partial {\rho }_c}>0,$$

(90)

or otherwise unstable if

$$\frac{\partial M\left({\rho }_c\right)}{\partial {\rho }_c}<0,$$

(91)

under radial pulsation. From Figure 19 and Figure 20, it is easy to see that the gravitational mass is an increasing function with regard to central density, which fulfills the Harrison–Zeldovich–Novikov static stability criterion under radial pulsation. We can also observe that the stellar structures become more massive as ${\rho }_c$ is increased, and we found that the solution picks up its stability with an increase in all chosen values of different parameters, viz., $\alpha$, $\beta$ and ${\omega }_{BD}$ in both scenarios BD in solid and BD+MGD in dashed as shown in Figure 19 and Figure 20.

8. Conclusions

The considered massive BD theory is one of the most extraordinary currently attainable alternative gravity theories that have huge deviations from the GRT on the one hand and are in concurrence with all the current observations on the other. Moreover, the massive BD theory is a well-posed theory which does not suffer from intrinsic problems and it is one of the best correct generalizations of GRT. This makes the solutions provided in this article and their stellar astrophysical implications very important.

In the present paper, we study the possibility of providing a new class of exact solutions for viable anisotropic stellar systems by means of massive BD gravity. In this regard, gravitational decoupling through the MGD approach in the arena of the embedding class-one space-time has effectively been employed to explore anisotropic solutions for matter sources. Our preliminary studies confirm that the decoupling of gravitational sources by MGD for compact stellar objects in massive BD gravity admitting the Karmarkar condition is a powerful mechanism to establish a more overall stellar solution of the modified Einstein field equations. For this purpose, we have constructed modified Einstein field equations given in (11)–(13), which disintegrate into two subsystems according to the MGD approach:

• The first subsystem (standard Einstein’s system) corresponds to the field equations given in (24)–(26) with five unknown parameters, $\{\Psi ,\xi ,\rho ,p_r,p_t\}$, and portrays the seed solution in BD gravity with respect to the stress-energy tensor $T_{\mu \nu }^m$.
• The second subsystem (quasi-Einstein system) corresponds to the field equations expressed in (32)–(34) and examines the solution of a new $\theta$-sector, which comprises mostly four unknown parameters, $\{h\left(r\right),{\theta }_0^0,{\theta }_1^1,{\theta }_2^2\}$.

Hence, in order to solve the complete stellar system identified by these two subsystems, we have proposed two approaches, namely, the mimic requirement on density and pressure component for $\theta$-sector, $\{{\theta }_0^0\left(r\right)=\rho \left(r\right)$ and ${\theta }_1^1\left(r\right)=p_r\left(r\right)\}$, respectively. For a detailed analysis for the massive stellar structure, we used the ansatz for $\Psi$ relating to the notable, well-recognized and well-comported Finch–Skea model via the embedding class-one condition, as well as the deformation function indicated by $h\left(r\right)$ via mimic requirement on the radial pressure component for anisotropic sector, viz., ${\theta }_1^1\left(r\right)=p_r\left(r\right)$. It is worth mentioning that we have analyzed the stellar structures supported by an anisotropic matter distribution in the background of massive BD gravity via the MGD approach for several values of constant parameters depending on $\alpha$, $\beta$ and ${\omega }_{BD}$. The unknown parameters appearing in the anisotropic stellar system have been determined via matching the interior astrophysical stellar geometry for the system to the Schwarzschild line element that describes the vacuum in the exterior of the stellar structure. Moreover, for studying the physical validity of the accomplished solutions under massive BD gravity, we have analyzed various physical properties of the anisotropic astrophysical structure via the MGD approach, which can be summarized as follows:

• We have observed in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 that all the thermodynamic observables, in particular, matter density, radial and transverse pressures with respect to the radial coordinate r have maximum values at the center and show limited, positive and monotonic decreasing comportment gradually toward the minimum values at the boundary of the celestial bodies, while the radial pressure becomes zero at the surface, which affirms the physical suitability and agreeability of the envisaged solutions. From these plots, we affirm additionally that our celestial bodies are completely free from any physical or mathematical singularities for all various parametric values of $\alpha$, $\beta$ and ${\omega }_{BD}$.
• The graphs corresponding to the comportment of anisotropy parameter $\Delta$ versus radial coordinate r are illustrated in Figure 9, Figure 10 and Figure 11. We see that the vanished anisotropy parameter $\Delta$ at the origin obtains positively defined increases to reach its maximum value at the surface of the stellar structure. Moreover, the fact that $\Delta >0$ i.e., $p_t>p_r$ implies that the stellar system encounters a repulsive force that counteracts the gravitational gradient, progressing the stability and equilibrium state, as well as one that has more compact and massive celestial bodies.
• The anisotropic stellar model is also consistent with the causality condition which is affirmed as the components of sound velocity lie within their prescribed bounds for all chosen values of the parameters $\alpha$, $\beta$ and ${\omega }_{BD}$, leading to a stable and viable compact stellar object under massive BD gravity via the MGD approach, as shown in Figure 12, Figure 13 and Figure 14.
• Our investigation of the $M-R$ curves is very important for compact stellar objects and shows the maximal bound for the celestial bodies. In Figure 15 and Figure 16, we present the behavior of the total mass M (in Solar mass $M_{\odot }$) with respect to the total radius R for the chosen specific values of parameters $\alpha$ and ${\omega }_{BD}$, by fixing $\beta$ to $0.2$. In the current massive BD gravity model via the MGD approach, we find from Figure 15 that, due to $\alpha$, as $\alpha$ increased, the maximum mass M increased with increasing radius R, which provides us with more massive compact stellar structures. On the other hand, as observed from Figure 16, due to ${\omega }_{BD}$, as ${\omega }_{BD}$ increased, the maximum mass M decreased with increasing radius R, which provides us also a stellar system that is more compact and massive. On the other hand, we highlight the effect of $\alpha$ and ${\omega }_{BD}$, by fixing $\beta$ to $0.2$, on the behavior of the maximum moment of inertia I against the total mass M, as illustrated in Figure 17 and Figure 18. From these graphs, we can conclude that the stiffness of the EoS is better through increasing the sensitivity of $I-M$ curves in both cases corresponding to the increase of the parameters, $\alpha$ and ${\omega }_{BD}$, which implies that the $I_{max}$ is always increasing from the zero up to a $M_{max}$ value, and thereafter decreases quickly. Finally, our generating $M-R$ and $I-M$ diagrams are well fitted with observational data, viz., LMC X-4, LMC X-4, Her X-1, 4U 1820-30, 4U 1608-52, SAX J1808.4-3658 and many other compact stellar structures that can be fitted. It is clear how both parameters $\alpha$ and ${\omega }_{BD}$ presented by massive BD gravity via the MGD approach incorporating the anisotropic profile of the matter distribution has a great influence on many physical parameters of the compact stellar structures. Furthermore, the predicted radii for different compact objects are mentioned in

  Table 1. Physical parameters of the observed compact stars for constants used in Figure 15 and Figure 16.

  		Predicted Radii (km)
  Strange Stars	Observed Mass	$\mathit{\alpha }$				${\mathit{\omega }}_{\mathbf{BD}}$
  	$\left({\mathit{M}}_{\odot }\right)$	0.3	0.5	0.7	0.9	0	500	1000
  Her X-1	0.85 ± 0.15 [88]	5.58	6.83	7.76	8.52	8.11	8.096	8.066
  SAX J1808.4-3658	0.9 ± 0.3 [89]	5.63	6.91	7.86	8.64	8.24	8.21	8.18
  SMC X-1	1.04 ± 0.09 [90]	5.82	7.39	8.50	9.41	8.95	8.89	8.83
  LMC X-4	1.29 ± 0.05 [90]	5.75	7.12	8.10	8.96	8.54	8.48	8.45
  4U 1820-30	1.58 ± 0.06 [91]	-	7.51	8.77	9.77	9.27	9.19	9.11
  4U 1608-52	1.74 ± 0.14 [92]	-	7.49	8.56	9.92	8.27	8.21	8.175

  for different $\alpha$ and ${\omega }_{BD}$. It is observed that, for lower values of $\alpha$, higher mass objects are ruled out. On the other hand, when $\alpha$ increases, compactness decreases, whereas compactness increases as ${\omega }_{BD}$ increases.
• It is well recognized that, for spherically symmetric static astrophysical systems, the Harrison–Zeldovich–Novikov static stability criterion plays a crucial role under radial pulsation, which must be satisfied to ensure stability. From Figure 19 and Figure 20, we observe that this static stability criterion is well satisfied under a radial perturbation, and we can also notice that the stellar structures become more massive as ${\rho }_c$ increases. Moreover, the stellar solution regains its stability with an increase in all chosen values of different parameters, viz., $\alpha$, $\beta$ and ${\omega }_{BD}$ in both scenarios BD and BD+MGD.

Furthermore, it is mentioned that all the obtained results in this manuscript reduce to Einstein’s GR when $\phi =constant$ and ${\omega }_{BD}\to \infty$.

Finally, we would like to mention here that all these new exact solutions for viable anisotropic stellar systems by means of massive BD gravity under gravitational decoupling through the MGD approach in the realm of the embedding class-one space-time are very much able to represent and describe the anisotropic compact celestial bodies.