1 Introduction

Seminal Tolman’s work [1] about spherically symmetric and static fluid spheres driven by isotropic matter distribution is considered as a milestone regarding to the seeking of analytic solutions of Einstein field equations, which describe exciting structures such as neutron stars, white dwarfs etc. All these objects corresponding to the final stage in the life of stars provide us real laboratories (impossible to recreate on the Earth) to understand the behaviour in the strong gravitational field regime. Interestingly, some features of these kind of configurations can be obtained from direct observations. For example, by measuring the surface gravitational red-shift we can infer the mass, M, and the radius, R, of a star that constitute the so called compactness parameter, \(u\equiv M/R\), which should be interpreted as quantification of how much mass can be packaged per unit radius before the gravitational collapse of the compact object takes place. In this respect, Buchdahl [2] determined that for any isotropic matter distribution the mass–radius ratio can not exceed the upper bound \(M/R\le 4/9\) that corresponds to a maximum gravitational surface red-shift \(z_{s}=2\).

However, as pointed out by Ruderman in his pioneering work [3], celestial bodies are not necessarily made of isotropic matter but they could contain local anisotropies at least in certain very high density ranges \((\rho >10^{15}\, \mathrm{g/cm}^{3})\), where the nuclear interactions must be treated relativistically. Similarly, the existence of super-Chandrasekhar white dwarfs that allow to explain overluminous type Ia supernovae demands that the Chandrasekar limit must be violated in the presence of strong magnetic fields which may be treated as an anisotropic fluid [4,5,6,7,8,9,10]. Indeed, it is well known that a magnetic field acting on a Fermi gas produces pressure anisotropy.

In this direction, the celebrated article by Bowers and Liang [11] laid the initial basis for the study of anisotropic structures within the framework of Einstein’s general relativity (GR from now on). Regardingly, they found that the contributions coming from local anisotropies into the Tolman–Oppenheimer–Volkoff (TOV) equation [1, 12] is of Newtonian origin.

The study of Buchdahl’s ratio in different frameworks involving self-gravitating compact structures entails intriguing queries, such as: (i) How is this limit modified in the presence of additional fields? (ii) How much do corrections to Einstein–Hilbert’s action influence mass–radius ratio? (iii) Does the Buchdahl limit increase or decrease in the high dimension regime? Regardingly, these questions have been extensively studied over the last two decades leading to interesting results [13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49] and among all the mechanism behind the modification of the Buchdahl’s limit, the simplest one corresponds to introduce local anisotropies within the matter distribution [50] which lead to \(u>4/9\), namely, stellar interiors with extra packing of mass. For example, in Ref. [49] it was shown that the introduction of a Kalb–Ramond (KR) field leads to a stable compact anisotropic configuration with extra packing of mass. It should be noticed that the KR field in [49] fill the whole space and as a consequence, the compact object is surrounded not by the Schwarzschild vacuum but by a background corresponding to an exterior solution of Einstein equations sourced by a KR field. In this respect, we may wonder if it is possible to induce extra packing of mass by introducing local anisotropies into the fluid distributions, without modifying the Schwarzschild vacuum space-time. Incidentally, the gravitational decoupling (GD) introduced in [51,52,53] in the framework of GR results to be a natural manner to induce such modifications.

As it is well known, GD is a generalization of the Minimal Geometric Deformation (MGD) employed in the context of Randall–Sundrum brane-world [54,55,56,57,58,59,60,61,62,63,64,65,66,67] and as we shall review in the next section, the methodology contains two main ingredients: (i) two sources, \(\bar{T}_{\mu \nu }\) (that usually corresponds to a perfect fluid) and \(\theta _{\mu \nu }\) which only interact gravitationally and (ii) a minimal geometric deformation introduced in the \(g^{rr}\) component of the metric which allows to decouple the system into two set of equations, one for each source. It is worth mentioning that in many applications of the GD, \(\bar{T}_{\mu \nu }\) is the source of a well known isotropic interior solution, so the effect of \(\theta _{\mu \nu }\) is to introduce a local anisotropy in the system and consequently, it is said that GD leads to an extension of isotropic solutions to anisotropic domains. In this case, either the geometric deformation and the components of \(\theta _{\mu \nu }\) remain unknown and the main goal is to provide suitable extra constraints which allow to find them. In this regard, a wide range of possibilities have been proposed, among which are: (i) the so-called mimic constraint procedure [52, 68,69,70,71,72,73,74], (ii) the imposition of an adequate decoupling function f(r) meeting all the physical and mathematical requirements [75, 76] or (iii) some anisotropy mechanism [77, 78].

Given the versatility of GD by MGD, its application to deal with a variety of situations has grown considerably during the last 2 years [53, 79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104]. However, although the modification of the Buchdahl limit in the framework of GD is straightforward, surprisingly it has not been considered yet to investigate the extra packing of mass of interior stellar. For this reason, the main goal of this work is to employ the Bömer and Harko methodology in [50] to study the Buchdahl limit for anisotropic compact stars but assuming that the anisotropy is encoded in the \(\theta \)-sector induced by the GD approach.

This work is organized as follows. In Sect. 2 we present the GD through MGD scheme. In Sect. 3 the Buchdahl limit generalities and its extension to anisotropic domains are revisited. Section 4 is devoted to introduce the extra packing of mass within the arena of MGD and the constraints on the \(\alpha \) parameter and the mass–radius ratios are determined for the minimally deformed Tolman IV space-time. Finally, Sect. 5 concludes the work.

Throughout the article we shall employ the mostly negative signature \((+,-,-,-)\).

2 Anisotropic sources: gravitational decoupling by MGD

In this section we provide a short and comprehensive introduction to GD by MDG. For more details see references [51, 52], for example.

Let us consider the Einstein–Hilbert (E–H hereinafter) action describing the gravitational field coupled to a matter field, through minimal coupling-matter principle given by

$$\begin{aligned} S=S_{E-H}+S_{M}, \end{aligned}$$
(1)

where the E–H action reads

$$\begin{aligned} S_{E-H}=\frac{1}{2\kappa }\int \sqrt{-g}Rd^{4}x, \end{aligned}$$
(2)

being \(g\equiv det(g_{\mu \nu })\) the determinant of the metric tensor \(g_{\mu \nu }\), R the Ricci scalar and \(\kappa =8\pi G c^{-4}\). From now on we shall employ relativistic geometrized units where Newton’s gravitational constant G and the speed of light are taken to be the unit ı.e, \(G=c=1\). The matter sector \(S_{M}\) is given by the following general expression

$$\begin{aligned} S_{M}=\int \sqrt{-g}{\mathcal {L}}_{M}d^{4}x, \end{aligned}$$
(3)

where \({\mathcal {L}}_{M}\) stands for the Lagrangian-density matter. So, let us write \({\mathcal {L}}_{M}\) as

$$\begin{aligned} {\mathcal {L}}_{M}={\mathcal {L}}_{\bar{M}}+\alpha {\mathcal {L}}_{X}, \end{aligned}$$
(4)

and consider that the information on isotropic, anisotropic or charged fluids, among others, are encoded in the Lagrangian-density \({\mathcal {L}}_{\bar{M}}\) (throughout the text we will use barred quantities to refer the usual material content), while \({\mathcal {L}}_{X}\) encipher the new matter fields. In principle these incoming fields can be seen as corrections to general relativity [53]. So, putting all together and taking variations respect to the inverse metric tensor \(\delta g^{\mu \nu }\) in (1), we arrive at the following field equations describing the gravitational–matter interaction

$$\begin{aligned} \frac{\delta S}{\delta g^{\mu \nu }}=0 \Rightarrow G_{\mu \nu } \equiv R_{\mu \nu }-\frac{R}{2}g_{\mu \nu }=-8\pi T_{\mu \nu }, \end{aligned}$$
(5)

where \(G_{\mu \nu }\) is the Einstein’s tensor. The general expression for \(T_{\mu \nu }\) after variations is given by

$$\begin{aligned} T_{\mu \nu }=\underbrace{-2\frac{\delta {\mathcal {L}}_{\bar{M}}}{\delta g^{\mu \nu }}+g_{\mu \nu }{\mathcal {L}}_{\bar{M}}}_{\bar{T}_{\mu \nu }} +\alpha \underbrace{\left( -2\frac{\delta {\mathcal {L}}_{X}}{\delta g^{\mu \nu }}+g_{\mu \nu }{\mathcal {L}}_{X}\right) }_{\theta _{\mu \nu }}. \end{aligned}$$
(6)

So,

$$\begin{aligned} T_{\mu \nu }=\bar{T}_{\mu \nu }+\alpha \theta _{\mu \nu }, \end{aligned}$$
(7)

being \(\alpha \) a dimensionless parameter. As we are interested in the study of compact structures describing a spherically symmetric and static space-time, we supplement the field equations (5) with the most general line element, namely

$$\begin{aligned} {ds}^{2}={{e}^{\nu }}{{dt}^{2}}-{{e}^{\eta }}{{dr}^{2}}-{r}^{2}({{d\theta }^{2}} +{{sin}^{2}}\theta {{d\phi }^{2}}). \end{aligned}$$
(8)

The staticity of this space-time (8) is ensured by considering \(\nu \) and \(\eta \) as functions of the radial coordinate r only. Putting together Eqs. (5), (7) and (8) one obtains the following system of equations

$$\begin{aligned}&e^{-\eta }\left( {\frac{\eta ^{{\prime }}}{r}}-\frac{1}{r^2}\right) +\frac{1}{r^2}=8\pi \left( \bar{T}^{0}_{0}+\alpha \theta ^{0}_{0}\right) , \end{aligned}$$
(9)
$$\begin{aligned}&-e^{-\eta } \left( {\frac{\nu ^{{\prime }}}{r}}+\frac{1}{r^2}\right) +\frac{1}{r^2} =8\pi \left( \bar{T}^{1}_{1}+\alpha \theta ^{1}_{1}\right) , \end{aligned}$$
(10)
$$\begin{aligned}&-\frac{ e^{-\eta }}{4} \left( 2 \nu ^{{\prime \prime }}+\nu ^{\prime 2} +2{\frac{\nu ^{{\prime }}-\eta ^{{\prime }}}{r}}-\nu ^{\prime }\eta ^{\prime } \right) =8\pi \left( \bar{T}^{2}_{2}+\alpha \theta ^{2}_{2}\right) , \end{aligned}$$
(11)

where primes denote derivation with respect to the radial coordinate and \(\bar{T}^{2}_{2}=\bar{T}^{3}_{3}\) in accordance with the spherical symmetry. Moreover, we shall assume that the source \(\bar{T}_{\mu \nu }\) is a perfect fluid which entails \(\bar{T}^{1}_{1}=\bar{T}^{2}_{2}\). The corresponding Bianchi’s identity (conservation law) \(\nabla ^{\mu }T_{\mu \nu }=0\) associated with the system (9)–(11) reads

$$\begin{aligned}&(\bar{T}^{1}_{1})'-\frac{\nu '}{2}\left( \bar{T}^{0}_{0}-\bar{T}^{1}_{1}\right) -\frac{2}{r}\nonumber \\&+\left[ (\theta ^{1}_{1})'-\frac{\nu '}{2}(\theta ^{0}_{0} -\theta ^{1}_{1})-\frac{2}{r}(\theta ^{2}_{2}-\theta ^{1}_{1})\right] =0 \end{aligned}$$
(12)

Note that Eq. (12) represents the generalized Tolman–Oppenheimer–Volkoff (TOV) equation driven the hydrostatic equilibrium of the system. At this point one can identify an effective energy density

$$\begin{aligned} \tilde{\rho }=\bar{T}^{0}_{0}+\alpha \theta ^{0}_{0}=\bar{\rho }+\alpha \theta ^{0}_{0}, \end{aligned}$$
(13)

an effective radial pressure

$$\begin{aligned} \tilde{p}_{r}=-\bar{T}^{1}_{1}-\alpha \theta ^{1}_{1}=\bar{p}-\alpha \theta ^{1}_{1}, \end{aligned}$$
(14)

and an effective tangential pressure

$$\begin{aligned} \tilde{p}_{\perp }=-\bar{T}^{2}_{2}-\alpha \theta ^{2}_{2}=\bar{p}-\alpha \theta ^{2}_{2}. \end{aligned}$$
(15)

Note that, Eqs. (9), (10) and (11) correspond to a set of differential equations in which we have only separated the components of the matter sector. However, in the framework of MGD these equations can be successfully decoupled and at the end, the system becomes in two set of differential equations, one for each source. To this end, we need to introduce the so-called minimal geometric deformation in the \(g^{rr}\) component of the metric given by the linear map

$$\begin{aligned} e^{-\eta (r)}\mapsto e^{-\eta (r)}=\mu (r)+\alpha f(r). \end{aligned}$$
(16)

Next, plugging Eq. (16) into Eqs. (9)–(11) one arrives to

$$\begin{aligned} 8\pi \bar{\rho }= & {} \frac{1}{r^{2}}-\frac{\mu }{r^{2}}-\frac{\mu ^{\prime }}{r} \end{aligned}$$
(17)
$$\begin{aligned} 8\pi \bar{p}= & {} -\frac{1}{r^{2}}+\mu \left( \frac{1}{r^{2}} +\frac{\nu ^{\prime }}{r}\right) \end{aligned}$$
(18)
$$\begin{aligned} 8\pi \bar{p}= & {} \frac{\mu }{4}\left( 2\nu ^{\prime \prime }+\nu ^{\prime 2} +2\frac{\nu ^{\prime }}{r}\right) +\frac{\mu ^{\prime }}{4} \left( \nu ^{\prime }+\frac{2}{r}\right) , \end{aligned}$$
(19)

along with the conservation equation

$$\begin{aligned} \bar{p}^{\prime }+\frac{\nu ^{\prime }}{2}\left( \bar{\rho } +\bar{p}\right) =0, \end{aligned}$$
(20)

for the isotropic sector. Similarly, one has the following equations for the \(\theta \)-sector

$$\begin{aligned} 8\pi \theta ^{0}_{0}= & {} -\frac{f}{r^{2}}-\frac{f^{\prime }}{r} \end{aligned}$$
(21)
$$\begin{aligned} 8\pi \theta ^{1}_{1}= & {} -f\left( \frac{1}{r^{2}}+\frac{\nu ^{\prime }}{r}\right) \end{aligned}$$
(22)
$$\begin{aligned} 8\pi \theta ^{2}_{2}= & {} -\frac{f}{4}\left( 2\nu ^{\prime \prime }+\nu ^{\prime 2} +2\frac{\nu ^{\prime }}{r}\right) -\frac{f^{\prime }}{4}\left( \nu ^{\prime } +\frac{2}{r}\right) . \end{aligned}$$
(23)

The corresponding conservation equation \(\nabla ^{\nu }\theta _{\mu \nu }=0\) yields

$$\begin{aligned} \left( \theta ^{1}_{1}\right) ^{\prime }-\frac{\nu ^{\prime }}{2} \left( \theta ^{0}_{0}-\theta ^{1}_{1}\right) -\frac{2}{r} \left( \theta ^{2}_{2}-\theta ^{1}_{1}\right) =0. \end{aligned}$$
(24)

At this point some comments are in order. First, note that as can be seen from Eqs. (9)–(11) the function \(\eta (r)\) is quite involved in the full set of equations. So, the only possibility to separate \(\{\bar{\rho },\bar{p}\}\) from \(\{\theta ^{0}_{0},\theta ^{1}_{1},\theta ^{2}_{2}\}\) is performing the linear map (16). Second, as it is well-known the gravitational mass function is defined in terms of \(g^{rr}\) as follows

$$\begin{aligned} m(r)=\frac{r}{2}\left[ 1-e^{-\lambda (r)}\right] . \end{aligned}$$
(25)

Now, after replacing the linear map (16) into Eq. (25) the gravitational mass function becomes

$$\begin{aligned} m(r)=\frac{r}{2}\bigg [1-\mu (r)-\alpha f(r) \bigg ], \end{aligned}$$
(26)

from where the total gravitational mass function can be expressed as

$$\begin{aligned} m(r)=m_{0}(r)-\alpha \frac{r}{2}f(r), \end{aligned}$$
(27)

with \(m_{0}(r)\) the mass function of the isotropic system. Hence, it is clear that MGD introduces a modification in the total mass \(m_{0}(R)=M_{0}\) of the original system. Even more, if \(f(r)>0\ (f(r)<0)\) and \(\alpha <0\ (\alpha >0)\) the original total mass \(M_{0}\) is increased by a quantity \(\alpha Rf(R)/2 \) and as a consequence the mass–radius ratio is also altered, namely

$$\begin{aligned} 2u(r)=2u_{0}(r)-\alpha f(r). \end{aligned}$$
(28)

Finally, it is remarkable that although the resulting set of equations for the \(\theta \)-sector (21)–(23), correspond to the so-called quasi Einstein field equations in the sense that there is a missing \(\frac{1}{r^{2}}\), Eq. (24) corresponds to the classical hydrostatic relativistic equilibrium equation associated to \(\theta _{\mu \nu }\).

An important thing to be noted here, is that the interaction between the two sources is completely gravitational. This fact is reflected by Eqs. (20) and (24), where both sectors are individually conserved, namely \(\nabla _{\mu }\bar{T}^{\mu \nu }=0\) and \(\nabla _{\mu }\theta ^{\mu \nu }=0\).

Note that, when an isotropic solution is considered, Eqs. (17), (18) and (19) are automatically satisfied. In this respect, we can use the metric function \(\nu \) to solve for the \(\theta \)-sector which consists in four unknowns, namely \(\{\theta ^{0}_{0},\theta ^{1}_{1},\theta ^{2}_{2},f\}\) and only three Eqs. (21)–(23). Hence, in order to close the system it is necessary to prescribe some additional information. Regardingly, the so-called mimetic constraints has been broadly used [52, 70,71,72,73,74] which consist in to impose

$$\begin{aligned} \bar{p}=\theta ^{1}_{1}, \end{aligned}$$
(29)

or

$$\begin{aligned} \bar{\rho }=\theta ^{0}_{0}. \end{aligned}$$
(30)

The former means that the anisotropic behavior enters into the matter distribution through the disturbance of the isotropic pressure of the seed solution. Besides, when (29) is imposed from Eqs. (18) and (22) one gets

$$\begin{aligned} f(r)=\frac{1}{1-r\nu ^{\prime }(r)}-\mu (r). \end{aligned}$$
(31)

It is worth noticing that restriction (29) entails an interesting situation in considering the total mass contained by the fluid sphere. As said before, MGD modifies the gravitational mass function (27), and as a consequence the total mass \(m(R)=M\) contained by the object. Surprisingly, when (29) is used, the extra piece \(\alpha rf(r)/2\) vanishes at the boundary \(r=R\), which allows to conclude that this mimic constraint does not modify the total mass. Indeed, a close look reveals that Eq. (18) at \(r=R\) provides

$$\begin{aligned} R\nu ^{\prime }(R)=\frac{1}{\mu (R)}-1, \end{aligned}$$
(32)

which after being replaced into (31) evaluated at the surface, it is obtained \(f(R)=0\). From the physical point of view this result can be interpreted as a redistribution of the mass of the object.

The situation is quite different when the mimetic constraint for the density (30) is employed. In this case, after replacing (17) and (22) we obtain

$$\begin{aligned} f(r)=\mu (r)-1. \end{aligned}$$
(33)

Now, as we are dealing with a spherical mass distribution the total mass is computed as follows

$$\begin{aligned} m(R)=4\pi \int ^{R}_{0}r^{2}\rho dr=4\pi \int ^{R}_{0}r^{2} \left( \bar{\rho }+\alpha \theta ^{0}_{0}\right) dr, \end{aligned}$$
(34)

where by virtue of (30) the above expression becomes

$$\begin{aligned} m(R)=4\pi \int ^{R}_{0}r^{2}\left( 1+\alpha \right) \bar{\rho } dr =\left( 1+\alpha \right) m_{0}(R), \end{aligned}$$
(35)

from where it is evident that if \(\alpha >0\) then \(m(R)>m_{0}(R)\). As we shall see later, in this work we base our analysis of extra packing of mass in the Tolman IV solution extended by MGD using the mimic constraint for the density previously reported in [52].

3 Buchdahl’s limit: isotropic and anisotropic sources revisiting

Buchdahl’s limit states that for a spherically symmetric and static configuration describing an isotropic matter distribution, the maximum mass–radius ratio u (also known as compactness factor) is given by

$$\begin{aligned} u\equiv \frac{M}{R}\le \frac{4}{9}, \end{aligned}$$
(36)

where M and R stand for the gravitational mass contained in the sphere and the radius of the star, respectively. As can be seen from (36) one has two options (i) an equality and (ii) an inequality. The first option holds by considering a constant energy density (incompressible fluid), a degenerate metric (in the sense that \(g_{tt}(r=0)=0\)) and a non well behaved pressure (divergent pressure at the center of the compact configuration). Otherwise, the inequality requires a positive and monotonously decreasing energy density from the center to the boundary of the compact object density everywhere inside the star and a vanishing pressure at the surface of the structure. Hence, the value 4/9 is an absolute upper limit for all static fluid spheres whose density does not increase outwards. The condition (36) works in the isotropic case but as it was said before celestial bodies are not necessarily made up by perfect fluid distributions. In this respect Böhmer and Harko [50] derived the corresponding upper bound for the mass–radius relation for an anisotropic matter distribution in presence of a cosmological constant \(\varLambda \) obtaining

$$\begin{aligned} \frac{2M}{R}\le \left( 1-\frac{8\pi }{3}\lambda R^{2}\right) \left[ 1-\frac{1}{9}\frac{\left( 1-2\lambda /\langle \rho \rangle \right) ^{2}}{\left( 1-\frac{8\pi }{3}\lambda R^{2}\right) \left( 1+F\right) ^{2}}\right] . \end{aligned}$$
(37)

In the above expression \(\langle \rho \rangle \) stands for the mean energy density and F is a function proportional to the anisotropy factor \(\varDelta =p_{\bot }-p_{r}\), which is given by

$$\begin{aligned} F=2\frac{\varDelta (R)}{\langle \rho \rangle }\Bigg [\frac{\arcsin \left( \sqrt{\frac{2M\chi \left( R\right) }{R}}\right) }{\sqrt{\frac{2M\chi \left( R\right) }{R}}}-1\Bigg ], \end{aligned}$$
(38)

where \(\chi \) is given by

$$\begin{aligned} \chi (r)\equiv 1+\frac{4\pi }{3}\varLambda \frac{r^{3}}{m(r)}. \end{aligned}$$
(39)

As we are interested in studying space-time without cosmological constant from now on we will set off \(\varLambda \). Then the upper bound (37) becomes

$$\begin{aligned} 2u\le 1-\frac{1}{9 (1+F)^2}, \end{aligned}$$
(40)

and F turns

$$\begin{aligned} F=\frac{R^2 \varDelta (u)}{\frac{3 u}{8 \pi }}\left( \frac{\sin ^{-1} \left( \sqrt{2 u}\right) }{\sqrt{2 u}}-1\right) . \end{aligned}$$
(41)

At this point a couple of comments are in order. First, note that F is a positive quantity. Otherwise, the condition \(\varDelta \ge 0\) could be violated. As a consequence, non-vanishing values of F, Eq. (40) allows extra packing of mass in compact stellar structures. Second, from Eq. (41), it is straightforward to show that the bounds of the compactness parameter of the anisotropic distribution, u, is defined in the interval

$$\begin{aligned} \frac{4}{9}\le u< \frac{1}{2}. \end{aligned}$$
(42)

In the next section we will discuss in details how the local anisotropies introduced in the stellar interior by gravitational decoupling through MGD approach contributes on the maximum mass–radius ratio value allowable for relativistic anisotropic fluid spheres in the arena of GR.

4 Local anisotropy induced by MGD

In the context of MGD, the anisotropy of the total matter sector seeded by a perfect fluid, can be written as

$$\begin{aligned} \varDelta =p_{\perp }-p_{r}=\alpha (-\theta ^{2}_{2}-(-\theta ^{1}_{1})). \end{aligned}$$
(43)

Now, from Eq. (40) it is straightforward to see that the connection between the anisotropy induced by gravitational decoupling and the compactness parameter is given by

$$\begin{aligned} 2 u\le 1-\frac{1}{9 \left( 1+\alpha \frac{16 \pi R^2 (\theta ^{1}_{1}-\theta ^{2}_{2})}{3 u}\left( \frac{\sin ^{-1} \sqrt{2 u}}{\sqrt{2 u}}-1\right) \right) ^2}, \end{aligned}$$
(44)

which leads to the classic Buchdahl’s limit for the isotropic case when \(\alpha \rightarrow 0\), as expected. However, expression (44) must be considered as a formal expression because matching conditions namely, the continuity of the first and the second fundamental form of the total solution leads to non-linear equations involving \(\alpha \) and u. In this sense, obtaining an analytical expression of the bound of the compactness parameter of the total solution, u, in terms of the decoupling parameter is not possible. Instead, we can try to find the connection between u and \(\alpha \) in an alternative manner. To this end, we shall consider the MGD-extended Tolman IV solution previously reported in [52]. As it is well known, the Tolman IV solution that solve the isotropic sector parametrized by \((\nu ,\mu ,\rho ,p)\) is given by

$$\begin{aligned} e^{\nu }= & {} B^{2}\left( 1+\frac{r^{2}}{A^{2}}\right) , \end{aligned}$$
(45)
$$\begin{aligned} \mu= & {} \frac{\left( 1-\frac{r^{2}}{C^{2}}\right) \left( 1+\frac{r^{2}}{A^{2}}\right) }{1+\frac{2r^{2}}{A^{2}}}, \end{aligned}$$
(46)
$$\begin{aligned} \rho= & {} \frac{3A^{4}+A^{2}(3C^{2}+7r^{2})+2r^{2}(C^{2}+3r^{2})}{8\pi C^{2}(A^{2}+2r^{2})^{2}}, \end{aligned}$$
(47)
$$\begin{aligned} p= & {} \frac{C^{2}-A^{2}-3r^{2}}{8\pi C^{2}(A^{2}+2r^{2})}, \end{aligned}$$
(48)

where AB and C are constants. Now, using \(\nu \) given by Eq. (45) the extended solution can be obtained by imposing extra conditions to close the system (21)–(23). As it was previously commented, in this work we shall concentrate in the so-called mimic constraint for the density which after replacing (46) in (33) leads to

$$\begin{aligned} f(r)=-\frac{r^{2}}{C^{2}} \left( \frac{A^{2}+C^{2}+r^{2}}{A^{2}+2r^{2}}\right) . \end{aligned}$$
(49)

Note that, at this point the \(\theta \)-sector can be completely determined in terms of the decoupling function given by (49). In particular, a straightforward computations reveals that the anisotropy function of the extended solution can be written as [52]

$$\begin{aligned} \varDelta (r;\alpha )=\alpha \left( \theta ^{1}_{1}-\theta ^{2}_{2}\right) =\frac{\alpha r^{2}}{8\pi (A^{2}+r^{2})^{2}}. \end{aligned}$$
(50)

Now, with all of this results at hand we can obtain the source of the extended solution, namely

$$\begin{aligned} \tilde{\rho }= & {} (\alpha +1) \bar{\rho }, \end{aligned}$$
(51)
$$\begin{aligned} \tilde{p}_{r}= & {} \bar{p}-\frac{\alpha \left( \left( A^2+3 r^2\right) \left( A^2+C^2+r^2\right) \right) }{8 \pi C^2 \left( A^2+r^2\right) \left( A^2+2 r^2\right) }, \end{aligned}$$
(52)
$$\begin{aligned} \tilde{p}_{\perp }= & {} \tilde{p}_{r} +\frac{\alpha r^2}{8 \pi \left( A^2+r^2\right) ^2}. \end{aligned}$$
(53)

In order to proceed with the study of the modification of the Buchdahl limit induced by MGD we need to compare the compactness parameter of the original isotropic solution, \(u_{0}=\frac{M_{0}}{R}\), and the compactness parameter of the anisotropic solution given by \(u=\frac{M}{R}\). Incidentally, in Ref. [52] it was shown that this quantities are related by

$$\begin{aligned} 2u=2u_{0}+\alpha \frac{R^{2}}{C^{2}} \left( \frac{A^{2}+C^{2}+R^{2}}{A^{2}+2R^{2}}\right) . \end{aligned}$$
(54)

At this point a couple of comments are in order. First, note that all the quantities involved in the previous expression are positive which entails that \(u>u_{0}\). Second, it is noticeable that when the Buchdahl limit is saturated in the isotropic sector, namely \(u_{0}=\frac{4}{9}\), the anisotropic solution acquire an extra packing in the sense that

$$\begin{aligned} u\ge \frac{4}{9}. \end{aligned}$$
(55)

However, note that the above corresponds to the critical situation in which \(u_{0}\) acquires its maximum allowed value. In this sense, we may wonder how the parameters should be set to ensure extra packing of mass in more realistic situations where \(u_{0}< 4/9\). It is worth mentioning that such a set is not trivial as we shall see in what follows.

The strategy we shall employ here corresponds to express the set of arbitrary parameters \(\{A,B,C\}\) in terms of \(\{R,u,u_{0},\alpha \}\) an then we demand

$$\begin{aligned} A> & {} 0 \end{aligned}$$
(56)
$$\begin{aligned} B> & {} 0 \end{aligned}$$
(57)
$$\begin{aligned} C> & {} 0. \end{aligned}$$
(58)

for some set \(\{u,u_{0},\alpha ,R\}\) such that

$$\begin{aligned} u_{0}< & {} \frac{4}{9} \end{aligned}$$
(59)
$$\begin{aligned} \frac{4}{9}< & {} u<\frac{1}{2} \end{aligned}$$
(60)
$$\begin{aligned} \alpha> & {} 0 \end{aligned}$$
(61)
$$\begin{aligned} R> & {} 0. \end{aligned}$$
(62)

In order to do so, we shall employ the Israel–Darmois junction conditions [106, 107] which consist in to demand the continuity of the first and second fundamental forms across the boundary \(\varSigma :r=R\) of the star. Now, for the continuity of the first fundamental form it is required that the inner manifold \({\mathcal {M}}^{-}\) is joined in a smoothly way with the outer space-time \({\mathcal {M}}^{+}\) described by the vacuum Schwarzschild solution

$$\begin{aligned} ds^{2}=\left( 1-2\frac{M_{\mathrm{Sch}}}{r}\right) dt^{2} -\left( 1-2\frac{M_{\mathrm{Sch}}}{r}\right) ^{-1}dr^{2}-r^{2}d\varOmega ^{2}. \end{aligned}$$
(63)

Hence,

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

For the present situation we have

$$\begin{aligned} B^{2}\left( 1+\frac{R^{2}}{A^{2}}\right)= & {} 1-2\frac{M}{R}, \end{aligned}$$
(65)
$$\begin{aligned} \left( 1+\alpha \right) \frac{\left( 1-\frac{R^{2}}{C^{2}}\right) \left( 1+\frac{R^{2}}{A^{2}}\right) }{1+\frac{2R^{2}}{A^{2}}}-\alpha= & {} 1-2\frac{M}{R}, \end{aligned}$$
(66)

where at \(\varSigma \) the Schwarzschild mass \(M_{\mathrm{Sch}}\) coincides with the total mass M contained by the fluid sphere. The second fundamental form entails a vanishing radial pressure \(\tilde{p}_{r}\) at the surface of the compact structure. Explicitly it reads

$$\begin{aligned} \tilde{p}_{r}(R)=0. \end{aligned}$$
(67)

Now, from Eqs. (54) and (65) we get

$$\begin{aligned} B= & {} \sqrt{\frac{1-2 u_{0}-\alpha \frac{R^2}{C^2} \left( \frac{A^2+C^2+R^2}{A^2+2R^2}\right) }{1+\frac{R^2}{A^2}}} , \end{aligned}$$
(68)

and from (67)

$$\begin{aligned} C=\sqrt{\frac{(1+\alpha ) \left( A^2+R^2\right) \left( A^2+3R^2\right) }{\left( A^2+R^2\right) -\alpha \left( A^2+3 R^2\right) }}. \end{aligned}$$
(69)

Next, using (69) in (54) we obtain

$$\begin{aligned} A=\frac{\sqrt{\alpha -3u(1+\alpha )+3u_{0}(1+\alpha )}}{\sqrt{(u-u_{0})(1+\alpha )}}. \end{aligned}$$
(70)

Without loss of generality we shall set \(R=1\)Footnote 1 in what follows, so Eqs. (68) and (69) read

$$\begin{aligned} B&=\sqrt{\frac{(2 u-1) (-\alpha +3 (\alpha +1) u-3 (\alpha +1) u_{0})}{\alpha -2 (\alpha +1) u+2 (\alpha +1) u_{0}}} \end{aligned}$$
(71)
$$\begin{aligned} C&=\sqrt{\frac{\alpha (-\alpha +2 (\alpha +1) u-2 (\alpha +1) u_{0})}{(u-u_{0}) ((\alpha -1) \alpha +2 (\alpha +1) u-2 (\alpha +1) u_{0})}}. \end{aligned}$$
(72)

Now, Eqs. (56), (57), (58), constrained by (59), (60) and (61) reduce to following extra conditions

$$\begin{aligned}&3 (\alpha +1) u<\alpha +3 (\alpha +1) u_{0} \end{aligned}$$
(73)
$$\begin{aligned}&2 (\alpha +1) u<\alpha +2 (\alpha +1) u_{0} \end{aligned}$$
(74)
$$\begin{aligned}&\alpha ^2+2 (\alpha +1) u<\alpha +2 (\alpha +1) u_{0}, \end{aligned}$$
(75)

from where \(u_{0}\) and \(\alpha \) acquire extra constraints, namely

$$\begin{aligned} u_{0}> & {} 0.38889 \end{aligned}$$
(76)
$$\begin{aligned} 0.2<\alpha< & {} 0.333 \end{aligned}$$
(77)

At this point some comments are in order. First, it is worth mentioning that (76) states a lower bound for the compactness parameter of the perfect fluid solution. Even more, for values out the interval

$$\begin{aligned} 0.38889<u_{0}<\frac{4}{9}, \end{aligned}$$
(78)

the anisotropic compactness u can not surpass the Buchdahl limit for isotropic solutions and and the extra packing of mass is not possible in this case. Second, the extra packing condition allows to restrict the appropriate values for the decoupling parameter as shown in (77). Finally, it is worth mentioning that constraints given by Eqs. (76) and (77) correspond to a particular case we obtained when the system (73), (74) and (75) is reduced. However, we considered the the most simplest case in order to illustrate how the process works.

As an example we shall study the behaviour of a solution with extra packing of mass considering \(u=0.45\), \(u_{0}=0.39\) and \(\alpha =0.3\). In Fig. 1 it is shown the behavior of the density, and the effective radial \(\tilde{p}_{r}\) and tangential \(\tilde{p}_{\perp }\) pressures of the anisotropic solution. Note that all the quantities shown in the profiles are monotonously decreasing and reach their maximum value at the center of the star, as expected. Besides, the radial \(\tilde{p}_{r}\) pressure vanishes at the surface of the star.

Fig. 1
figure 1

Decreasing behavior of density and pressures of the anisotropic solution

Now we analyse the stability of the solution using the adiabatic index criteria. As it is well known, the radial and the tangential adiabatic index, \(\varGamma _{r}\) and \(\varGamma _{\perp }\)s respectively, defined as

$$\begin{aligned} \varGamma _{r}= & {} \frac{\tilde{\rho }+\tilde{p}_{r}}{\tilde{p}_{r}} \frac{d\tilde{p}_{r}}{d\tilde{\rho }} \end{aligned}$$
(79)
$$\begin{aligned} \varGamma _{\perp }= & {} \frac{\tilde{\rho }+\tilde{p}_{\perp }}{\tilde{p}_{\perp }} \frac{d\tilde{p}_{\perp }}{d\tilde{\rho }}, \end{aligned}$$
(80)

allow to connect the relativistic structure of a spherical static object and the equation of state of the interior fluid and serves as a criterion of stability. More precisely, it is said that an interior configuration is stable whenever \(\varGamma _{r}\) and \(\varGamma _{\perp }\) are grater that 4/3 for static equilibrium. Otherwise, it is said that the object is unstable and the self gravitating object undergoes a gravitational collapse.

As shown in Fig. 2, the behaviour of the adiabatic index reveals that for \(u=0.45\), we obtain a stable interior configuration from the anisotropic solution. Furthermore, note that \(\varGamma _{r}>\varGamma _{\perp }\) in accordance with the monotonously increasing of the anisotropy given that \(p_{r}\ge p_{\perp }\).

Fig. 2
figure 2

Adiabatic index showing stability for the anisotropic interior solution case (greater than 4/3)

In this sense, we conclude that MGD method not only allow to extend isotropic solutions to anisotropic domains but it can be used to obtain a physical acceptable interior with extra packing of mass. Furthermore, the extra packing condition leads to modifications in the interval of the compactness parameter of the isotropic solution. More precisely, given that \(0.38889<u_{0}<\frac{4}{9}\), not any arbitrary isotropic solution allows being extended with extra packing but only those with a high compactness comparable with observational data of neutron stars. Finally, it should be emphasized that the extra packing condition leads also to a stringent restriction in the decoupling parameter \(\alpha \) as shown in Eq. (77).

5 Concluding remarks

In this work we implemented the Gravitational Decoupling in the framework of the Minimal Geometric Deformation approach to induce extra packing of mass in stellar interiors. As anisotropic solution we considered the Tolman IV extended by using the mimetic constraint for the density. We obtained that although the induced anisotropy enters naturally in the Bömer and Harko criteria to induce extra packing, obtaining an analytical relationship between the decoupling parameters and the compactness of the star is not possible. However, we were able to set the space of parameters in an appropriate manner to obtain a well posed anisotropic solution with extra packing. We conclude that the Minimal Geometric Deformation approach allows to extend isotropic solutions to anisotropic interiors with extra packing of mass. Besides, the extra packing condition induces a lower bound on the compactness parameter of the isotropic system used as a seed and the decoupling parameter gets a restriction within a well defined interval.