1 Introduction

The stability of general relativistic self-gravitating matter distributions has been extensively studied and reported in the literature through several techniques for many years. It is complex multivariate problem which depends on the micro-physics – bulk/shear viscosity, crust on the surface, magnetic field and so on – of the material constituents and their description through a macroscopic equation of state that characterizes the configuration (see standard texts in Relativistic Astrophysics and Neutron Stars [1,2,3,4,5] and references therein).

The studies of stability distinguish two different approaches associated to the global and local scales, where instabilities affect the structure. On the global scale it is examined through the dynamical perturbation scheme which, in the case of spherical symmetry, can be translated into how radial pulsations induce possible disruptions of a stellar body. On the other hand, local stability investigates the effects of convection and/or cracking of the material within the matter distribution. Pulsation is a global phenomenon characterized by the collective motion of the entire body while, convection/cracking occurs locally and is governed by the nearby values of the thermodynamical variables and their gradients (see an interesting review in [6]).

The dynamical instability approach studies the evolution of perturbations on the physical and geometrical variables throughout the matter distribution. In General Relativity, it arose from the seminal works of Chandrasekhar [7, 8], Tooper [9, 10] and Bardeen [11] and, a decade later, was formalized by Friedman and Schutz [12]. For an anisotropic fluid, this criterion bounds the adiabatic index as

$$\begin{aligned} \Gamma = \frac{\rho + P}{P} v^{2} \ge \frac{4}{3} \, , \end{aligned}$$
(1)

where \(\rho \) denotes the energy density, P the radial pressure, and \(v^{2}\) the radial sound speed, respectively [13,14,15,16]. Within the global stability criteria we can also identify the Harrison–Zeldovich–Novikov condition, which implies that \(\mathrm {d}M(\rho _c) /\mathrm {d} \rho _c \ge 0\), where M is the total mass of the configuration and \(\rho _c\) the central density, of the distribution [16].

The stability of a spherical star to convection implies the buoyancy principle which leads to that pressure and energy density must, decrease outwards in any hydrostatic matter configuration [17,18,19]. Finally, the cracking instability approach determines the tidal acceleration profiles generated by perturbations of the energy density and the anisotropy of pressures identifying the changes sign of the total force distribution within the system [20,21,22,23]. This approach has been applied to an anisotropic fluid with barotropic equations of state [24] and, more recently extended to take into account the perturbation of the pressure gradient in both isotropic and anisotropic matter configurations [25, 26].

In the present paper we consider convection and cracking instabilities as well as their interplay. We develop a simple criterion to identify equations of state unstable to convection, and also explore the influence of buoyancy on cracking (or overturning) of isotropic and anisotropic relativistic spheres.

This paper is organized as follows: Sect. 2 describes the general equations of the theoretical framework of General Relativity. In Section 3 we formulate the concepts of adiabatic stability while cracking for self-gravitating anisotropic matter configurations is discussed in Sect. 4. Next we present in Sect. 5 the acceptability conditions which make any model physically reasonable. The models used and a discussion of our results for isotropic and anisotropic cases are presented in the Sects. 6 and 7. Finally, in Sect. 8 we wrap-up our final remarks.

2 The field equations

Let us consider a spherically symmetric space-time whose line element is given by

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

with the regularity conditions at \(r=0: e^{2\nu (0)}= \text{ const. }, \,\, e^{-2\lambda (0)}= 1, \,\, {\nu }'(0)=\lambda '(0)=0\,. \)

At the surface of the sphere \(r=r_b\), the interior solution should match continuously the exterior Schwarzschild solution, which implies that: \(e^{2\nu (r_b)}=e^{-2\lambda (r_b)}=1-2\mu _b\), and the compactness at the surface defined as \(\mu _b=2M/r_b\).

We shall consider a distribution of matter consisting of a non-Pascalian fluid with the energy–momentum tensor:

$$\begin{aligned} T_\mu ^\nu = \text{ diag }\left[ \rho (r),-P(r),-P_\perp (r),-P_\perp (r) \right] \,, \end{aligned}$$
(3)

with energy density \(\rho (r)\), radial pressure P(r) and tangential pressure \(P_\perp (r)\) of the fluid being determined by the Einstein field equations as

$$\begin{aligned} \rho (r)= & {} \frac{ e^{-2\lambda }\left( 2 r \lambda ^{\prime }-1\right) +1 }{8\pi r^{2}}\,, \end{aligned}$$
(4)
$$\begin{aligned} P(r)= & {} \rho (r) - \frac{e^{-2\lambda }}{4\pi r^{2}}\left[ r\left( \lambda ^{\prime }-\nu ^{\prime }\right) + e^{2\lambda }-1\right] \, \qquad \text {and} \end{aligned}$$
(5)
$$\begin{aligned} P_\perp (r)= & {} -\frac{e^{-2\lambda }}{8\pi }\left[ \frac{ \lambda ^{\prime }-\nu ^{\prime }}{r}-\nu ^{\prime \prime }+\nu ^{\prime }\lambda ^{\prime }-\left( \nu ^{\prime }\right) ^2\right] \,, \end{aligned}$$
(6)

where primes denote differentiation with respect to r.

As it is well known, \(T^{\mu }_{r \; ; \mu } =0\) implies the hydrostatic equilibrium equation, or Tolman–Oppenheimer–Volkoff (TOV) equation which, if \(m(r)=\frac{r}{2} \left( 1- e^{-2\lambda } \right) \), can be written for this anisotropic fluid as

$$\begin{aligned} \frac{\mathrm {d} P}{\mathrm {d} r} +(\rho +P)\frac{m + 4 \pi r^{3}P}{r(r-2m)} -\frac{2}{r}\left( P_\perp -P \right) =0\,, \end{aligned}$$
(7)

and together with

$$\begin{aligned} \frac{\mathrm {d} m}{\mathrm {d} r}=4\pi r^2 \rho \,, \end{aligned}$$
(8)

constitute the stellar structure equations. In order to obtain the density and pressure profiles we have to provide two equations of state, \(P=P(\rho )\) and \(P_{\perp }=P_{\perp }(\rho )\), which for the present work will be assumed barotropic.

3 Adiabatic convection stability condition

The stability of a spherical star against convection can be easily understood. When a fluid element is displaced downward, if its density, \(\rho _{e}\), increases more rapidly than the surrounding density, \(\rho _{s} \), the element will sink downward and the star will be unstable. On the other hand, if the density of fluid element is less than its surroundings, it will float back and the star will be stable to convection.

Thus we can single out these three cases:

  1. 1.

    If \( \rho _{e}> \rho _{s} \), gravity will tend to push the fluid element downward further and the system will be unstable.

  2. 2.

    If \( \rho _{e} = \rho _{s} \), the system is considered neutral or metastable.

  3. 3.

    If \( \rho _{e} <\rho _{s} \), a restoring force will act on the fluid element and then the system will be stable because it tends to its original state.

Following Bondi [17], let us denote the density, \(\rho (r_p)\), of an infinitesimal fluid element at its original position \(r_p\) and displace this piece of material downward, thus:

$$\begin{aligned}&\rho (r_p) \rightarrow \rho (r_p)+\delta \rho (r)\,,\,\, \nonumber \\&\quad \text {with} \,\,\, \delta \rho (r)=\rho '(r)(-\delta r) \,\, \text {and} \,\,\, r=r_p-\delta r\,, \end{aligned}$$
(9)

where r represents the current position of the fluid element, \(r_p \) its original position and \(-\delta r\) the downward shift.

Because \(\rho '(r) <0\), then \( \delta \rho (r) \) is a positive quantity, and the density of the compressed fluid element at the new displaced position will be greater that the density at its original position \(r_p\). On the other hand, expanding the density of the environment at the displaced position we get:

$$\begin{aligned} \rho (r_p-\delta r) \approx \rho (r_p)+\rho '(r_p)(-\delta r) \, , \end{aligned}$$
(10)

The system will be stable against convection if the environment density is greater or equal than the density of the fluid element, we then have:

$$\begin{aligned} \rho (r_p)+\rho '(r_p)(-\delta r) \ge \rho (r_p)+\rho '(r)(-\delta r), \end{aligned}$$
(11)

thus \( \rho '(r_p)\le \rho '(r) \, . \)

Now, expanding \(\rho '(r)\) around \(r_p\) we get

$$\begin{aligned} \rho '(r_p)+\rho ''(r_p)\delta r \le \rho '(r_p) \,\, \Rightarrow \,\, \rho ''(r) \le 0, \end{aligned}$$
(12)

which becomes the criterion of adiabatic stability against convection. Thus, density profiles with the second derivative less or equal than zero, \(\rho ''(r) \le 0\), will be stable against adiabatic convective motions.

It is clear that parabolic density profiles, \(\rho = \alpha r^2 + \beta \), with \(\alpha \) and \(\beta \) constants, will be stable against this type of convection, because buoyancy condition must also be fulfilled at the center \(r = 0\) of the sphere: \(\rho '< 0 \,\, \Rightarrow \, \rho '_c=0,\) and \( \rho ''_c<0 \). These profiles have been implemented for the MIT Bag model through a linear equation of state, \(P = \beta (\rho -\rho _{s})\), when densities become high enough for a phase transition to quark matter to occur [3, 16].

4 Convection and cracking sources of instability

In this section we shall consider convective instabilities in the framework of cracking induced by perturbation in the density profile. Just for completeness we outline here the main concepts and equations concerning cracking for isotropic and anisotropic matter configurations, for further details, we refer interested readers to [25, 26] and references therein.

As in the previous works [25, 26], we assume that density fluctuations induce variations into all other physical variables, i.e. \(m(r), P(r), P_\perp (r)\) and their derivatives, generating a non-vanishing total radial force distribution (\(\delta \mathcal {R} \ne 0\)) within the configuration. It is important to stress we are considering local perturbations of density, that can be properly described by any function of compact support, \(\delta \rho = \delta \rho (r)\), defined in a closed interval \(\Delta r \ll r_b\), where \(r_b\) is the total radius of the configuration.

Accordingly, local density perturbations, \(\rho \rightarrow \rho + \delta \rho \), generate fluctuations in mass, radial pressure, tangential pressure and radial pressure gradient, that can be represented up to linear terms in density fluctuation as:

$$\begin{aligned} \delta P= & {} \frac{\mathrm {d}P}{\mathrm {d} \rho } \delta \rho = v^2 \delta \rho \,, \end{aligned}$$
(13)
$$\begin{aligned} \delta P_\perp= & {} \frac{\mathrm {d}P_\perp }{\mathrm {d} \rho } \delta \rho = v_\perp ^2 \delta \rho \,, \end{aligned}$$
(14)
$$\begin{aligned} \delta P'= & {} \frac{\mathrm {d} P'}{\mathrm {d} \rho } \delta \rho = \frac{\mathrm {d \ }}{\mathrm {d} \rho } \left[ \frac{\mathrm {d} P}{\mathrm {d}r} \right] \delta \rho = \frac{\mathrm {d}}{\mathrm {d}\rho } \left[ \frac{\mathrm {d}P}{\mathrm {d}\rho } \frac{\mathrm {d}\rho }{\mathrm {d} r} \right] \delta \rho \nonumber \\= & {} \frac{\mathrm {d}}{\mathrm {d}\rho } \left[ v^2 \rho ' \right] \delta \rho = \frac{1}{ \rho '} \frac{\mathrm {d}}{\mathrm {d}r} \left[ v^2 \rho ' \right] \delta \rho \nonumber \\= & {} \left[ (v^2)' + v^2\frac{\rho ''}{\rho '} \right] \delta \rho \,, \end{aligned}$$
(15)
$$\begin{aligned} \delta m= & {} \frac{\mathrm {d} m}{\mathrm {d} \rho } \delta \rho = \frac{\mathrm {d} m}{\mathrm {d} r} \left( \frac{\mathrm {d} r}{\mathrm {d} \rho } \right) \delta \rho = \frac{m'}{ \rho '} \delta \rho = \frac{4 \pi r^2 \rho }{\rho '}\delta \rho \,. \nonumber \\ \end{aligned}$$
(16)

where

$$\begin{aligned} v^{2} = \frac{\mathrm {d} P }{ \mathrm {d} \rho } \,\, \mathrm {and} \,\, v_\perp ^2 = \frac{\mathrm {d} P_\perp }{ \mathrm {d} \rho }\,, \end{aligned}$$
(17)

are the radial and tangential sound speeds, respectively.

Note that the present perturbation scenario contrasts the original one presented by Herrera and collaborators [20,21,22], where fluctuations in density and anisotropy were considered independent and simultaneous; and it is also different from a previous work [24] because there pressure gradient were not affected by the density perturbation.

Following [26], we formally expand the quantity \( \mathcal {R}\) emerging from the TOV equation as:

$$\begin{aligned} \mathcal {R} \equiv \frac{\mathrm {d} P}{\mathrm {d} r} +(\rho +P)\frac{m + 4 \pi r^{3}P}{r(r-2m)} -\frac{2}{r}\left( P_\perp -P \right) \, , \end{aligned}$$
(18)

as

$$\begin{aligned} \mathcal {R}\approx & {} \mathcal {R}_{0}(\rho , P, P_\perp , m, P^{\prime }) + \frac{\partial \mathcal {R}}{ \partial \rho } \delta \rho +\frac{\partial \mathcal {R}}{ \partial P} \delta P + \frac{\partial \mathcal {R}}{\partial P_\perp } \delta P_\perp \nonumber \\&+\frac{\partial \mathcal {R}}{ \partial m} \delta m +\frac{\partial \mathcal {R}}{ \partial P^{\prime }} \delta P' \,, \end{aligned}$$
(19)

where \(\mathcal {R}_{0}(\rho , P, P_\perp , m, P^{\prime })=0\), because initially the configuration is in equilibrium. Next, by using (13)–(16) the above Eq. (19) can be reshaped as:

$$\begin{aligned} \delta \mathcal {R}\equiv & {} \delta \underbrace{P'}_{\mathcal {R}_p} + \delta \underbrace{\left[ (\rho +P)\frac{m + 4 \pi r^{3}P}{r(r-2m)} \right] }_{\mathcal {R}_g} + \delta \underbrace{\left[ 2\frac{P}{r}- 2\frac{P_\perp }{r}\right] }_{\mathcal {R}_a} \nonumber \\= & {} \delta \mathcal {R}_p + \delta \mathcal {R}_g +\delta \mathcal {R}_a \, , \end{aligned}$$
(20)

where it is clear that density perturbations \(\delta \rho (r)\) are influencing: the distribution of reacting pressure forces \(\mathcal {R}_p\), gravity forces \(\mathcal {R}_g\) and anisotropy forces \(\mathcal {R}_a\). Depending on this effect, each perturbed distribution force can contribute in a different way to the change of sign of \(\delta \mathcal {R}\): each term can be written as

$$\begin{aligned} \delta \mathcal {R}_p= & {} \left[ \frac{P''}{\rho '} \right] \delta \rho = \left[ (v^2)' + v^2\frac{\rho ''}{\rho '}\right] \delta \rho \end{aligned}$$
(21)
$$\begin{aligned} \delta \mathcal {R}_g= & {} \left[ \frac{ \partial \mathcal {R}_g }{ \partial \rho } + \frac{ \partial \mathcal {R}_g }{ \partial P } v^2 + \frac{ \partial \mathcal {R}_g }{ \partial m } \frac{4 \pi r^2 \rho }{\rho '}\right] \delta \rho \quad \text {and} \end{aligned}$$
(22)
$$\begin{aligned} \delta \mathcal {R}_a= & {} 2\left[ \frac{ v^2 - v^2_\perp }{ r } \right] \delta \rho , \end{aligned}$$
(23)

with

$$\begin{aligned}&\frac{\partial \mathcal {R}_g}{\partial \rho } = \frac{m+ 4 \pi r^3 P }{ r(r - 2m)}\,, \,\, \frac{\partial \mathcal {R}_g}{\partial P} = \frac{ m + 4 \pi r^3 ( \rho + 2 P )}{r(r-2m)} \quad \nonumber \\&\quad \text {and} \quad \frac{\partial \mathcal {R}_g}{\partial m} = \frac{ (\rho + P)( 1 + 8\pi r^2 P) }{(r- 2 m )^2 }. \end{aligned}$$
(24)

Notice that if, as in [24], the perturbation \(\delta \rho \) is constant and does not affect the pressure gradient, we have: \(\delta \mathcal {R}_p = 0\),

$$\begin{aligned} \delta \tilde{\mathcal {R}}_g = 2\left[ \frac{m+4 \pi r^3 P }{ r(r - 2m)} +\frac{2 \pi r^2}{3} \frac{ (\rho + P)( 1 + 8\pi r^2 P) }{( 2 m - r )^2}\right] \delta \rho \,. \end{aligned}$$
(25)

Thus, only anisotropic matter distribution can present cracking instabilities because \(\delta \tilde{\mathcal {R}}_g >0\) for all r and the possible change of sign for \( \delta \mathcal {R}\) should emerge from \(\delta \mathcal {R}_a\) and the criterion against cracking is written as:

$$\begin{aligned} -1 \le v^2_\perp -v^2 \le 0 \quad \Leftrightarrow \quad 0 \ge \frac{\mathrm {d} P_\perp }{\mathrm {d} r} \ge \frac{\mathrm {d} P}{\mathrm {d} r} \, ; \end{aligned}$$
(26)

more recently this equivalence between the restriction on pressures and velocities was demonstrated in [16] and included as part of the acceptability conditions that have to be considered when building physically reasonable compact object models. We shall discuss these constraints in the next section.

5 Physical acceptability conditions

More of the spherically symmetric perfect fluid “exact solutions” of Einstein field equations found in the literature are of little physical interest because, in addition to solving the structure equations (7) and (8) for a particular set of equations of state –\(P=P(\rho )\) and \(P_{\perp }=P_{\perp }(\rho )\)– the physical and metric variables have to comply with several acceptability conditions which, over the years were recently compiled in [16] as:

  1. C1:

    Metric potentials, positive, finite and free from singularities;

  2. C2:

    Matching conditions at the surface of the star;

  3. C3:

    Decrease of interior redshift Z with the increase of r;

  4. C4:

    Positive density and pressures;

  5. C5:

    Density and pressures having a maximum at the center and decreasing monotonically outwards, with \(P_{\perp } \ge P\);

  6. C6:

    Energy conditions. Strong (SEC) \(\rho \ge P + 2P_{\perp }\) or dominant (DEC) \(\rho \ge P\) and \(\rho \ge P_{\perp }\);

  7. C7:

    Causality Conditions. \(0 \le v^2,v_{\perp }^2 \le 1\);

  8. C8:

    The adiabatic index \(\Gamma \) stability criterion as stated in Eq. (1), which is a consequence of the dynamical stability criterion;

  9. C9:

    Stability against cracking as expressed by Eq. (26);

  10. C10:

    Harrison–Zeldovich–Novikov stability condition: \(\mathrm {d}M(\rho _c)/\mathrm {d}\rho _c > 0\).

Additionally Ivanov [16], demonstrated that these conditions are not independent and can be condensed in five main inequalities:

  1. 1.

    \((m/r)' > 0\) which fulfills C1, C2 and C3 conditions;

  2. 2.

    \(0 \ge P_{\perp }' \ge P'\) accomplishing: C4, C5, C6 (SEC), C7, C9;

  3. 3.

    \(\mu = 2M/r_b \le 4/5\) executing C6 (DEC);

  4. 4.

    \(v^2 \ge 1/3\), implementing C8;

  5. 5.

    and the C10 condition: \(\mathrm {d}M(\rho _c)/\mathrm {d}\rho _c > 0\).

Clearly if we want the configuration to be stable against convection, we should add a sixth condition, i.e. the adiabatic convection stability condition \(\rho '' \le 0\), to the above mentioned set. In the next section, we shall explore its influence on the stability of isotropic and anisotropic models.

6 Isotropic and anisotropic models

In this section we select seven exact solutions – describing isotropic or anisotropic fluid spheres – which comply with the Ivanov criteria. With this selection we study the effect of convective instability and the reaction of the pressure gradient to density perturbations.

Four of seven density profiles are among most physically reasonable isotropic solutions (Tolman VII [27], Buchdahl-1 [28], Mehra [29] and Kuchowicz [30]) reported in [31]. The fifth selection corresponds to a one-parameter family of a generalized Tolman IV solution obtained in [32] and allows us to exemplify the correlation of convection stability with a decreasing profile of \((v^2)'\). Finally, we study two anisotropic solutions (Gokhroo and Mehra [33] and Sah and Chandra [34]) to illustrate the buoyancy effects in anisotropic matter configurations.

6.1 Isotropic solutions

In addition complying with the Ivanov criteria, isotropic model selection has significant physical interest in describing the interior of compact objects. The isotropic solutions shown in Table 1 are: Tolman VII and Mehra are the most frequent parabolic density profiles considered in models of stable neutron stars (see [35,36,37,38] and reference therein) while Buchdahl’s solution [28, 39] sets limits to the compactness of relativistic spheres. In [31] it is reported that the speed of sound for this solution does not decrease monotonically. We have shown that it can be attained for some particular values of the compactness \(\mu = 2M/r_b\). Finally, models of charged spheres are frequently based of Kuchowicz solution [40] and in Mehra’s solution [29] the density and the speed of sound vanish on the surface.

Table 1 Parabolic density profiles, represented by Tolman VII and Mehra, are used extensively to model stable neutron stars. In these profiles: A, R and C are constants to be determined by the boundary conditions. Rational profiles like Buchdahl sets limits to the compactness, \(\mu = 2M/r_b\) of relativistic spheres, and \(r_b\) is the boundary radius of the matter configuration. Kuchowicz solution models of charged spheres, in this case \(x=Ar^2\), the function \(F(r) = \mathrm {Ei}\left( 1 +\frac{Ar^2}{2} \right) /\mathrm {e}\) where \(F(0)= \mathrm {Ei}(1)/\mathrm {e}\) and \(\mathrm {Ei}\) is the exponential integral function, again, A and C are constants to be determined by the boundary conditions
Full size table

N-parametric isotropic Lake–Tolman IV family of solutions. In reference [32], Lake proposed an algorithm based on the choice of a single monotone seed-function, \(\mathcal {F}\), to generate all regular static spherically symmetric perfect fluid solutions of Einstein’s equations. Within this scheme we are going to consider the following family of models generated by

$$\begin{aligned} \mathcal {F}(r)=1+Cr^2 \,\, \Rightarrow \, \nu =\frac{N}{2}\ln \left[ 1+Cr^2\right] , \end{aligned}$$
(27)

where C is an arbitrary constant and N is a positive integer that produces an infinite family of analytic solutions. As we can see from the Eq. (27), different values of N recover well-known solutions: \(N = 1\) corresponds to Tolman IV solution [27]; \(N = 3\) represents Heint IIa [41] solution; \(N = 4\) and \(N=5\) are Durg IV and Durg V solutions, respectively [42]. The case \(N = 2\) is considered in [43] studying the relationship between the central barionic density and the total mass of observed neutron stars.

As can be easily guessed, the main difficulty of the method lies in how to calculate the two integrals that appearing in equation (4) of [32], but fortunately it is possible in the present case, as can seen in the appendix.

Thus, with the help of the following auxiliary function

$$\begin{aligned} \mathcal {G}(r)=1+C{r}^{2}(N+1)\,, \end{aligned}$$
(28)

we obtain the density profile for any N, as

$$\begin{aligned} \rho (r)= \frac{1}{8\pi r^{2}}\left[ \left( 1-\frac{r^2}{2}\,\frac{\Pi _1}{ \mathcal {F}^{N-2}} \right) \left( 2r\Pi _2-1 \right) +1 \right] \,, \end{aligned}$$
(29)

where:

$$\begin{aligned} \Pi _1= & {} {\frac{C\left( N-2 \right) {N}^{N-2} }{ \left( N+1 \right) ^{N-3}} \Phi }+\frac{4K}{\mathcal {G}^{\frac{2}{N+1}} }\,,\\ \Pi _2= & {} \frac{r}{2\mathcal {F}^{N-2}}\\&\left[ \frac{\left[ 1- \frac{(\mathcal {F} -1) ( N-2 )}{\mathcal {F}}\right] \Pi _1 +2(\mathcal {F} -1) \left[ {\frac{C ( N-2) ( N-3 ) }{ ( N+1 )^{N-4} \left( N+ 3 \right) N^2}\Phi }-\frac{4K}{ \mathcal {G}^{\frac{N+3}{N+1}} } \right] }{1-\frac{(\mathcal {F} -1)}{2C}\,\frac{\Pi _1}{ \mathcal {F}^{N-2}} }\right] \end{aligned}$$

and

$$\begin{aligned} \Phi = {_2\hbox {F}_1}\left( 3-N,\frac{2}{ N+1};\,{\frac{N+3}{N+1}};\, -{\frac{\mathcal {G} }{N}}\right) \,. \end{aligned}$$

Here, \({_2\hbox {F}_1}(a,b;c;d)\) is the hypergeometric function, that for certain special arguments: (a, b; c; d) automatically evaluates to exact values with K a constant.

Table 2 Florides–Stewart–Gokhroo–Mehra (FSGM) and Sah–Chandra profiles satisfy Ivanov criteria. In these profiles: \(\alpha \), K and a are constants to be determined by the boundary conditions
Full size table
Fig. 1
figure 1

Convection for isotropic models: left plates, a, c, illustrate the normalized buoyancy \(\hat{\rho }''~=~\rho ''/\rho _{c}''\) vs \(\hat{r} = r/r_b\) and right plates, b, d, display the gradient of the sound velocity \((v^2)'\) vs \({\hat{r}}\) for different values of \(\mu =2M/r_b\). a and b correspond to the Buchdahl model and c and d to Kuchowicz solution. As it can be appreciated from these plots, Buchdahl model becomes unstable to convective perturbations because \(\hat{\rho }''\) changes it sign, while Kuchowicz is stable

Full size image
Fig. 2
figure 2

Convection for isotropic N-model family: normalized buoyancy \(\hat{\rho }''\) vs \(\hat{r}\) (plates a–c–e) and gradient of the sound velocity \((v^2)'\) vs \(\hat{r}\) (plates b–d–f), for different values of N and \(\mu = 2M/r_b\). For \(N=1\), and \(N=3\) (Tolman IV [27] and Heint IIa [41], respectively) are unstable and \((v^2)'\) change sign. The \(N=6\) Model is stable to convection having \((v^2)'< 0\) for all \(\hat{r}\). It is interesting to mention, that for this family of solutions, when N increases, we obtain more convective stable models

Full size image
Fig. 3
figure 3

Convection for anisotropic models: the plate a illustrates, for different values of \(\mu = 2M/r_b\), the gradient of the sound velocity \((v^2)'\) for the FSGM model while Sah and Chandra’s model, in plates b and c, illustrate the buoyancy \({\rho }^{\prime \prime }\) and \((v^2)'\). Clearly, FSGM has a parabolic density profile and it is stable to convection. Notice in the case of \(\mu = 0.7, 0.6\) Sah and Chandra’s solution is unstable in the outer parts of the configuration and \((v^2)'\) is negative. This is because theses models do not comply Eq. (30)

Full size image
Fig. 4
figure 4

Stability against cracking, isotropic N-models: in these figures we show, the stabilizing effect of the reaction of pressure gradient to density perturbation. In these plots we present \({\delta \mathcal {R}}/{ \delta \rho }\) for two members of the Lake N-Family. As in the previous case this stabilizing effect is only present for \(N=1\), for \(N\ge 2\) no instability appears. Again, the more N increases the most stable the configuration is

Full size image

6.2 Anisotropic solutions

Local anisotropy (unequal stresses: \(P\ne P_{\perp }\)) in compact objects can be associated to different physical scenarios –such as phase transition, density inhomogeneity and electromagnetic field, just to mention a few of them – and has been considered extensively since the work of Bowers and Liang [44]. The unknown physics in the tangential equation of state, \(P_{\perp } = P_{\perp }(\rho )\) is partially compensated by using heuristic criteria: geometric, simplicity or any other assumption relating radial and tangential pressures (see [15, 45] and references therein).

Our selection of anisotropic solutions shown in the Table 2 are: Florides–Gokhroo–Mehra [33, 46, 47] and the Sah and Chandra [34].

The Florides–Gokhroo–Mehra profile was due originally to Florides [46], but also corresponds one of the different solutions considered by Stewart [47] and, more recently, by Gokhroo and Mehra [33]. The Florides–Stewart–Gokhroo–Mehra (FSGM) solution represents densities and pressures which, under particular circumstances [48], give rise to an equation of state similar to the Bethe–Börner–Sato newtonian equation of state for nuclear matter [1, 2, 49].

7 Modeling performed and discussion of some results

The first effect to be analyzed is the convective instability and its relation to the sign of the gradient of the speed of sound, \((v^2)'\), which are illustrated in Figs. 1, 2 and 3. We show this effect by plotting the normalized buoyancy \(\hat{\rho }'' =\rho ''/\rho _{c}''\) vs \(\hat{r} = r/r_b\) and the corresponding gradient \((v^2)'\) vs \(\hat{r}\). The first two figures display the convective stability for isotropic models while Fig. 3 illustrates this property for anisotropic spheres.

As we have mentioned before, parabolic density profiles – represented by Tolman VII, Mehra and Florides–Stewart–Gokhroo–Mehra in Tables 1 and 2 – are stable to convection because they have constant \(\rho ''< 0\). The stability for the other isotropic models having rational density profiles is presented in Fig. 1, where the Buchdahl model becomes unstable because \(\rho ''/\rho _c\) changes sign, while Kuchowicz is stable against convective perturbations due to \(\rho ''/\rho _c > 0\) for all r.

In Fig. 3 we illustrate the convective instability for anisotropic matter configurations. Again, models that are stable with the Ivanov criteria are revealed unstable for convection. This is the case of then Sah and Chandra model which is unstable against convective perturbation, but the Florides–Stewart–Gokhroo–Mehra solution is stable because it has a parabolic density profile.

Fig. 5
figure 5

Stability to cracking in isotropic and anisotropic models: In these plots we show, for different values of \(\mu = 2M/r_b\), the stabilizing effect of the reaction of pressure gradient to density perturbation. Plates a, b correspond to Tolman VII models; plates c, d to the Mehra models and e, f to the FSGM model. We compare \(\delta \mathcal {R}/ \delta \rho \) as a function of \(\hat{r}\) when the reaction of the pressure gradient to density perturbation are considered (plates b, d, f) and when they are not (plates a, c, e). It clear that these models are unstable to cracking when the reaction of the force distribution is not considered, i.e. \(\delta \mathcal {R}_p =0\). Notice that, contrary to what was shown in [26] Mehra models are stable when \(\delta \mathcal {R}_p\) is taken into account. This may happen because a possible mismatch in pressure and density for this solution in this reference

Full size image

In all the isotropic models analyzed we found an interesting correlation between the stability against convective perturbation and the sign of \((v^2)'\). In particular, in the regions where matter configuration has \(\rho '' < 0\) then we have \(v^2\) as a monotonous decreasing function. And this is more evident in Fig. 2 where we have ploted the buoyancy and the gradient (\(v^2\))’. We found unstable configurations for \(N=1\), Tolman IV [27] and for \(N=3\) Heint IIa [41] with \((v^2)'\) changing sign within the configuration and for \(N=6\) a convective stable configuration having \((v^2)'< 0\) for all \(\hat{r}\). This can be easily understood in those regions where \(P^{''}<0\) because

$$\begin{aligned} \rho ^{''}= & {} -\frac{(v^2)'}{(v^2)^2}P'+ \frac{P''}{v^2}\,, \,\, \text {thus, if} \,\, P^{''}< 0 \; \wedge \; (v^2)'< 0 \nonumber \\&\Rightarrow \,\, \rho ^{''} < 0\,. \end{aligned}$$
(30)

Therefore, a distribution having a monotonically decreasing concave pressure profile, i.e. \(P' < 0\) and \(P'' < 0\), will be stable against convection, if the sound velocity monotonically decreases outward from the configuration. It is interesting to see that when it happens we have

$$\begin{aligned} (v^2)' = \frac{\partial ^2 P}{\partial \rho ^2}\rho '\,, \,\, \text {thus,} \,\, \rho '< 0 \,\,\Rightarrow \,\, (v^2)' < 0 \quad \text {if } \; \frac{\partial ^2 P}{\partial \rho ^2} > 0\, , \end{aligned}$$
(31)

and this last condition is not only met by the model we have considered here, but also by any relativistic polytropic equation of \(P = K \rho ^{\gamma }\) or \(P = (\gamma -1)\rho \) and by other several more realistic numeric equation of state for ultradense matter (see [50], and references therein).

The other effect to be discussed is the stability induced by the reaction of the pressure gradient to density perturbations. This is shown in Figs. 4 and 5, where we compare the perturbation on the total force distribution when the pressure gradient is perturbed and when it is not. In these figures we plot \(\delta \mathcal {R}_p \ne 0\), and \(\delta \mathcal {R}_p = 0\), respectively. The first case, \(\delta \mathcal {R}_p \ne 0\), represents the recent cracking scheme of Gonzalez–Navarro–Nunez [25, 26] while the second one corresponds to a variation of the previous work of Abreu–Hernandez–Nunez [24]. As we have stressed above, in the first case it is assumed a non-constant density perturbation \(\delta \rho = \delta \rho (r)\), which leads to the factor \( \frac{4 \pi r^2 \rho }{\rho '} \le 0 \) in the third term in the perturbation to the gravitational force distribution \(\delta \mathcal {R}_g\) in Eq. (22), which may cause cracking instabilities.

As it is evident for the models considered, if the pressure gradient is not perturbed, i.e. \(\delta \mathcal {R}_p~=~0\), \(\delta \mathcal {R}\) may change its sign and potential cracking instabilities may appear. On the other hand, if the gradient reacts to the perturbation, \(\delta \mathcal {R}_p \ne 0\), we find that \(\delta \mathcal {R}\) does not change sign and the matter configuration becomes stable against cracking. This tendency to make models cracking-stable if the the pressure gradient reacts to the density, was previously reported incidentally in reference [26].

This induced stability can be understood if we shape the perturbation of the hydrostatic equation. Clearly, for the isotropic case, \(\mathcal {R}_a =\delta \mathcal {R}_a = 0\) and when \(\delta \mathcal {R}_p = 0\), the cracking instability emerges only from the effect of the perturbation on the gravitational force distribution, \(\delta \mathcal {R}_g\) and particularly from the third term, which is always negative. If \(\delta \mathcal {R}_p \ne 0\) and \(P'' > 0\), the reaction of the pressure gradient to the density perturbation can neutralize the effect of the negative sign of the above mentioned gravitational term.

8 Conclusions and final remarks

Stability is a key concept when considering self-gravitating stellar models, because only those in stable equilibrium are of astrophysical interest. A we have stated above, in addition to solving the structure equations (7) and (8) for a particular set of equations of state: \(P=P(\rho )\) and \(P_{\perp }=P_{\perp }(\rho )\), the emerging physical variables have to comply with the several acceptability conditions stated in Sect. 5.

The stability of a spherical star against convection has almost been forgotten in most of the stability analysis. It is a very simple criterion which implements the Archimedes principle in any hydrostatic matter configuration [17,18,19]. This criterion has proven to be interesting because several reasonable models become unstable against convection and should be included in the acceptability criteria to guarantee physically interesting models of compact objects.

In this work we have shown that:

  1. 1.

    A density profiles with its second derivative with respect to the radial marker less or equal than zero, \(\rho '' < 0\), will be stable against convective motions. This is a very simple criterion to identify potential convection instabilities within spherical matter configurations and it should be added to the above mentioned acceptability Ivanov criteria.

  2. 2.

    A decreasing concave pressure profile P(r), i.e. \(P' < 0\) and \(P'' < 0\), will be stable against convection if the radial sound velocity decreases monotonically outward in spherically matter configurations. Equivalently, if \( \frac{\partial ^2 P}{\partial \rho ^2} > 0\), then \((v^2)' < 0\) implies stability against convective motions. It should be pointed out that to illustrate this effect we have implemented a new family of exact solution by using the method propose by in Ref. [32].

  3. 3.

    From (19) we obtained the possible sources of cracking: the reaction of the pressure gradient \(\delta \mathcal {R}_p\); the perturbation to the gravitational force distribution \(\delta \mathcal {R}_g\) and the perturbation of the anisotropy \(\delta \mathcal {R}_a\). Convection may cause cracking instability only when the perturbation of the pressure gradient is considered, i.e. \(\delta \mathcal {R}_p \ne 0\), in this case cracking density perturbations affect the pressure gradient, depending on the values of \((v^2)'\) and \(v^2\frac{\rho ''}{\rho '}\).

  4. 4.

    Isotropic and anisotropic models considered can be unstable to cracking when the reaction of the pressure gradient is neglected, i.e. \(\delta \mathcal {R}_p = 0\), but if taken into account, the instabilities may vanish. Thus, there is a stabilizing effect against cracking, when the perturbation gradient is affected by density perturbation, \(\delta \mathcal {R}_p \ne 0\).

Local perturbed schemes are based on the reaction of the fluid variables to a density fluctuation that drives the system out of its equilibrium, i.e. we are exclusively considering perturbations under which the system is dynamically unstable. One way to achieve this is to assume for a barotropic fluid, that pressure gradients are not affected and we have shown that this occurs in several of the models considered. Convection contributes to the pressure gradient reaction but instead of developing or increasing cracking, it stabilizes the configuration.

Finally, it should be stressed that the results we have presented of possible instabilities for local perturbation schemes (convection and/or cracking) have to be considered as tendencies that could lead potential evolution of fluids within a relativistic matter configuration, but this should emerge from the full integration of the Einstein Equations.