1 Introduction

Recent statistical analyses of astrophysical and cosmological datasets have once again confirmed the concordance \(\varLambda \)CDM model [1,2,3,4]. Despite its successes, the model shows some shortcomings. On one side, the fundamental nature of the two most important energy density components, namely Dark Energy (DE) and Dark Matter (DM), is still unknown [5,6,7,8]. Many candidates have been proposed without being able to solve the puzzle [9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25]. On the other hand, it is well known that General Relativity (GR) is not a Quantum Theory of gravity and it cannot provide a description of the Universe at the quantum scales needed to solve the fine-tuning of initial conditions [26,27,28]. As a consequence, many modified theories of the gravity have been proposed to solve the puzzle [29,30,31,32,33,34,35,36,37]. Nevertheless, having alternative explanations demands to test the modified gravity models and other basic tenets of the \(\varLambda \)CDM cosmology at all scales, both in the strong and weak field regime. In particular, let us remember that the constraints at the Solar System scale must be matched by any theory of gravity under consideration [38,39,40].

Here, we will compute the weak field limit of energy–momentum-squared-gravity (EMSG), recently introduced by [41, 42], to study the collapse of a self-gravitating system. The main idea behind EMSG is to resolve the Big Bang singularity in a non-quantum description. It is important to mention that, GR inherently leads to the singularity in the early universe. On the other hand, as already mentioned, in the early universe, i.e., at the Planck scale, the quantum gravity effects play an important role. Therefore GR predicts space-time singularity in a physical situation in which its viability is seriously doubted. EMSG’s action functional is obtained by adding scalar terms proportional to \(T_{\mu \nu }T^{\mu \nu }\) (where \(T_{\mu \nu }\) is the energy–momentum tensor) to the Einstein-Hilbert action, and it leads to interesting cosmological behaviours This kind of corrections, naturally induce squared contributing terms like \(\rho ^2\), \(p^2\) and \(\rho p\) to the Friedman equations governing the background cosmological evolution. Where \(\rho \) and p are the energy density and pressure of the cosmic fluid. As a consequence, there are bouncing cosmological solutions in this model, and the cosmic scale factor cannot be smaller than a minimal length scale. In other words, there is a finite maximum energy density. This directly means that EMSG can prevent the Big Bang singularity in a completely non-quantum way. More importantly, EMSG does not alter the cosmological evolution. Its only effect is to resolve the singularity (for more details we refer the reader to [41]). However, it is should be stressed that EMSG’s effects can appear also in the stellar configurations. For example, it is shown in [43] that EMSG can lead to more massive neutron stars than in GR. This fact is satisfactory in the sense that there are difficulties in GR for explaining the internal structure of massive neutron stars, especially their high mass, using ordinary equations of state (for more details see [44] and [45]).

The study of the collapse of a self-gravitating system is somehow the first test to do to probe any modified theory of gravity. Indeed, Jeans Instability for a spherically symmetric self-gravitating systems causes the collapse of a gas cloud under the gravitational force giving rise to the formation of self-gravitating structure such as stars and galaxies among the others [46]. For stability, the cloud must be in hydrostatic equilibrium, and this physical condition holds only on certain scales determined by the so-called Jeans length, \(\lambda _J^2 = \frac{c_s^2}{2 G\rho }\) where \(c_s\) is the sound speed, G is the Newton’s gravitational coupling constant and \(\rho \) is the matter density. All perturbations having wavelengths larger than it are unstable. On the contrary, smaller wavelengths are stable. Since such a scale is strongly dependent by the underlying theory of gravity, it has been used to probe several modified theories of gravity [47,48,49,50].

Besides the stability criteria for spherically symmetric perturbations, [51] investigated the stability condition of all local axisymmetric perturbations introducing the dimensionless parameter \(Q=\frac{c_s \kappa }{\pi G \varSigma }\), where \(\kappa \) is the epicyclic frequency and \(\varSigma \) is the surface density of the system. Thus, any cloud or disk is stable if the condition \(Q>1\) holds. As it was for Jeans instability, it has been shown that also the Toomre’s criterion can be used to check the validity of several modified theories of gravity [52,53,54]. Generalizing both criteria for the local stability in the framework of EMSG could provide a very remarkable tool to describe the dynamics of self-gravitating system such as the collapse of spherical clouds, the collapse of massive star into Black Hole and/or the accretion disks around a massive object, leading to new results that could potentially be used to retain/rule out the theory.

The paper is organized as follows: in Sect. 2 we briefly introduce the EMSG and derive its field equations. In Sect. 3 we perform the weak field limit of EMSG. In particular, we write down the modified Poisson’s equation. In Sect. 4, we give the modified continuity and Euler equation for EMSG. In Sects. 5 and 6, we compute the Jeans’s length and the Toomre parameter for EMSG, respectively. In both cases we compute and analyze the dispersion relation particularizing our calculation to specific cases of the EMSG. In Sect. 7, we analyze the stability of an exponential disk to recover the Toomre’s criteria and, then, in Sect. 8 we apply our calculations to the case of Hyper Massive Neutron Stars (HMNS). Finally, in Sect. 9 we summarize our results and conclusions.

2 Field equations of EMSG

As in any other theory of gravity, the starting point of the EMSG is the action

$$\begin{aligned} S=\frac{1}{2\gamma }\int f(R,\mathbf {T}^2)\sqrt{-g}\,d^4 x +\int {\mathcal {L}}_M \sqrt{-g}\,d^4x, \end{aligned}$$
(1)

where \(\gamma =8\pi G/c^4\), G Newton’s constant, c is the speed of light, \(\sqrt{-g}\) is the determinant of metric tensor, \({\mathcal {L}}_M\) is the matter Lagrangian density, \(\mathbf {T}^2=T_{\mu \nu }T^{\mu \nu }\), and \(T_{\mu \nu }\) is the energy–momentum tensor [41]. Notice that, we use the metric signature \((-,+,+,+)\). Working in the metric formulation of the theory, it is straightforward to vary the action with respect to the metric and find the following field equations

$$\begin{aligned} f_R R_{\mu \nu }-\frac{1}{2} g_{\mu \nu }f=\gamma T_{\mu \nu } -\Big [f_Q \theta _{\mu \nu }+(g_{\mu \nu }\square -\nabla _{\mu }\nabla _{\nu })f_R\Big ], \end{aligned}$$
(2)

where \(f=f(R,\mathbf {T}^2)\), \(f_R=\partial f/\partial R\) and \(f_Q=\partial f/\partial Q\), the \(\square \) is the usual d’Alembert operator, and for simplicity in notation we have defined \(Q\equiv \mathbf {T}^2\). On the other hand, the tensor \(\theta _{\mu \nu }\) is defined as the variation of Q with respect to the metric tensor, namely \(\theta _{\mu \nu }=\delta Q/\delta g_{\mu \nu }\). For a perfect fluid the energy–momentum tensor and \(\theta _{\mu \nu }\) are written as follows (for more details see [55, 56])

$$\begin{aligned} T_{\mu \nu }=\left( \rho +\frac{p}{c^2}\right) u_{\mu }u_{\nu }+g_{\mu \nu }p, \end{aligned}$$
(3)

and

$$\begin{aligned}&\theta _{\mu \nu }=-\left( \rho ^2 c^2+4p \rho +3 \frac{p^2}{c^2}\right) u_{\mu }u_{\nu } \end{aligned}$$
(4)
$$\begin{aligned}&Q=\rho ^2 c^4+3 p^2 \end{aligned}$$
(5)

where \(\rho \) and p are the energy density and pressure of the perfect fluid, respectively. Moreover \(u^{\mu }\) is the four velocity of the fluid. Before moving on to discuss the weak field limit of the theory, let us take the trace of field Eq. (2). The result is written as

$$\begin{aligned} f_R R-2 f=\gamma T-(f_Q \theta +3\square f_R), \end{aligned}$$
(6)

where \(\theta =g^{\mu \nu }\theta _{\mu \nu }\). Now, we have all the equations needed to perform the weak field limit of EMSG.

3 Weak field limit of EMSG

Let us compute the first order perturbations of the field equations around the Minkowski space time in order to find the governing equations for the Newtonian self-gravitating disk in the context of EMSG. To do so we write the line element in the Cartesian coordinate (ctxyz) using the perturbed metric, i.e. \(g_{\mu \nu }=\eta _{\mu \nu }+h_{\mu \nu }\) where \(|h_{\mu \nu }|\ll |g_{\mu \nu }|\), as follows

$$\begin{aligned} ds^2=-\left( 1+\frac{2\,\varPhi }{c^2}\right) c^2dt^2 +\left( 1+\frac{2\,\varPsi }{c^2}\right) (dx^2+dy^2+dz^2). \end{aligned}$$
(7)

The corresponding first order perturbations in other quantities can be written as

$$\begin{aligned}&Q=Q^0+\delta Q,\end{aligned}$$
(8)
$$\begin{aligned}&R=R^0+\delta R, \end{aligned}$$
(9)
$$\begin{aligned}&\theta _{\mu \nu }=\theta ^0_{\mu \nu }+\delta \theta _{\mu \nu }, \end{aligned}$$
(10)
$$\begin{aligned}&f=f^0+f^0_R \delta R+f^0_Q \delta Q, \end{aligned}$$
(11)
$$\begin{aligned}&f_R=f_R^0+f_{RR}^0\delta R +f_{RQ}^0 \delta Q, \end{aligned}$$
(12)
$$\begin{aligned}&f_Q=f_Q^0+f_{RQ}^0\delta R +f_{QQ}^0 \delta Q, \end{aligned}$$
(13)

where the suffix “0” indicates the background quantities, and \(f_{XY}=\partial ^2 f/\partial X\partial Y\). For the Minkowski background we have \(T_{\mu \nu }^0=0\) and consequently \(Q^0=0\). We assume that the function f(RQ) is chosen in a way that if \(T_{\mu \nu }^0=0\) in the background then the background Ricci scalar \(R^0=0\) and \(f^0=f(0,0)=0\). In this case, using the definitions of \(T_{\mu \nu }\) and \(\theta _{\mu \nu }\), it is straightforward to verify that

$$\begin{aligned}&\delta T_{\mu \nu }\simeq \rho u^0_{\mu }u^0_{\nu }, \end{aligned}$$
(14)
$$\begin{aligned}&\delta \theta _{\mu \nu }\simeq -\rho ^2 c^2 u^0_{\mu }u^0_{\nu }, \end{aligned}$$
(15)

where commonly we have assumed that in the weak field limit \(p/\rho c^2 \ll 1\). Furthermore one should note that the background velocity four-vector is given by \(u^0_{\mu }=(-c,0,0,0)\) Now, let us substitute perturbed quantities in given Eqs. (8)–(13) into Eqs.  (2) and (6). Keeping only the first order terms, Eq. (6) takes the following form

$$\begin{aligned} f^0_{RR}\square \delta R+f^0_{RQ}\square \delta Q =\frac{\gamma }{3}\delta T-\frac{f^0_Q}{3}\delta \theta +\frac{f^0_R}{3}\delta R+\frac{2f^0_{Q}}{3}\delta Q, \end{aligned}$$
(16)

where \(\delta T\) and \(\delta \theta \) are perturbations in T and \(\theta \) respectively. Hereafter, for brevity in notation and prevent confusion with temporal components of the tensors, we drop the “0” suffix. Now we use Eq. (16) to linearize the time-time component of the field Eq.  (2) as follows

$$\begin{aligned} \delta R^0_0=\frac{\gamma }{f_R}\Big (\delta T^0_0-\frac{1}{3}\delta T\Big ) -\frac{f_Q}{f_R}\Big (\delta \theta _0^0-\frac{\delta \theta }{3} +\frac{\delta Q}{6}\Big )+\frac{\delta R}{6}, \end{aligned}$$
(17)

using a standard gauge it is straightforward to show that \(\delta R^0_0=-\nabla ^2 \varPhi /c^2\). To be precise, the standard gauge is also commonly called standard gauge of post-Newtonian theory (for more details see Sec. 8.3.7 in [57]). In particular, this gauge condition allow us to simplify the perturbed field equations at first order. Using this gauge condition we obtain two results: firstly we eliminate the higher order time derivatives of metric tensor, and secondly the Poisson equations can be solved more easily.

On the other hand, in principle, we can consider the perturbed Ricci scalar as a function of \(\rho \) and p, i.e. \(\delta R=\delta R(\rho ,p)\). To see this fact more clearly, let us conveniently assume that \(f_{RQ}=0\). We will use this assumption everywhere in this paper. This means that we only deal with models that can be recast in the following form

$$\begin{aligned} f(R,Q)=f_1(R)+f_2(Q). \end{aligned}$$
(18)

In this case, keeping in mind that \(\square \delta R=\nabla ^2 \delta R\), we rewrite Eq. (16) as

$$\begin{aligned} \nabla ^2 \delta R-{\mathcal {M}}^2 \delta R= H(\rho , p), \end{aligned}$$
(19)

where the mass \({\mathcal {M}}^2\) is defined as

$$\begin{aligned} {\mathcal {M}}^2=\frac{f_R}{3 f_{RR}}, \end{aligned}$$
(20)

and the function H is

$$\begin{aligned} H(\rho ,p)=\frac{1}{3f_{RR}} \Big [\gamma \delta T-f_Q(\delta \theta -2 \delta Q)\Big ]. \end{aligned}$$
(21)

It should be noted that it is natural to expect that \(f_R\) in the background is unity. However for the sake of completeness we keep it as a free parameter in the calculations. Furthermore let us define new parameter \(\alpha \) as

$$\begin{aligned} \alpha =\frac{f_Q}{\gamma }. \end{aligned}$$
(22)

Consequently, for the general form for f(RQ) in Eq. (18) with \(f_{RR}\ne 0\), we can integrate Eq. (19) to obtain \(\delta R\) in terms of \(\rho \) and p

$$\begin{aligned} \delta R=\chi \int \frac{e^{- {\mathcal {M}} |\mathbf {r}-\mathbf {r}'|}}{|\mathbf {r} -\mathbf {r}'|}\Big [- \delta T(\mathbf {r}')+\alpha (\delta \theta (\mathbf {r}') -2\delta Q(\mathbf {r}'))\Big ]d^3 \mathbf {r}', \end{aligned}$$
(23)

where for convenience, we have defined \(\chi \equiv \frac{\gamma {\mathcal {M}}^2}{4\pi f_R}.\)

Therefore, using Eq. (17), the modified version of Poisson’s equation in EMSG can be written as

$$\begin{aligned} \nabla ^2 \varPhi =\frac{\gamma \, c^4}{2}\tilde{\rho }, \end{aligned}$$
(24)

where we have defined the density \(\tilde{\rho }\) as

$$\begin{aligned} \tilde{\rho }=-\frac{2}{f_R}\Big (\delta T^0_0-\frac{1}{3}\delta T\Big ) +\frac{2\alpha }{ f_R}\Big (\delta \theta _0^0-\frac{\delta \theta }{3} +\frac{\delta Q}{6}\Big )-\frac{\delta R}{3\gamma \,c^2}. \end{aligned}$$
(25)

On the other hand by using Eqs. (14) and (15) we have

$$\begin{aligned} \delta T^0_0\simeq -\rho c^2, ~~~~\delta T\simeq -\rho c^2,~~~\delta Q=\delta \theta =\delta \theta ^0_0\simeq \rho ^2 c^4, \end{aligned}$$
(26)

where we have applied the condition \(p\ll \rho c^2\) in the weak field limit. In GR we have \(f_Q=0\). Moreover in this case we have \(\delta R=-\gamma \delta T\). Consequently it is easy to show that \(\tilde{\rho }=\rho \), and Eq. (24) recovers the standard Poisson’s equation. For another special case, the EMSG model studied in [41] is given by \(f(R,Q)=R-\eta \mathbf {T}^2=R-\eta Q\). For this model we have \(f_{RR}=0\), \(f_R=1\), and \(f_Q=-\eta =\alpha \gamma \). Moreover from Eq. (16), one may simply verify that \(\delta R=-\gamma (\delta T+\alpha (2\delta Q-\delta \theta ))\). Therefore Eq. (25) gives

$$\begin{aligned} \tilde{\rho }=\rho (1+2\alpha \rho c^2), \end{aligned}$$
(27)

in this special case, the effects of EMSG can be included in the effective density and pressure defined as, see [43]

$$\begin{aligned}&\rho _{\mathrm{eff}} =\rho +\frac{\alpha c^2}{2}\left( 8\rho \frac{p}{c^2} +\rho ^2+3\frac{p^2}{c^4} \right) , \end{aligned}$$
(28)
$$\begin{aligned}&p_{\mathrm{eff}}=p+\frac{\alpha c^4}{2}\left( \rho ^2+3\frac{p^2}{c^4} \right) , \end{aligned}$$
(29)

More specifically, it has been shown in [43] that the governing equations of EMSG, are completely similar to GR and the only difference is that \(\rho \) and p are replaced with \(\rho _{\mathrm{eff}}\) and \(p_{\mathrm{eff}}\). In this case, in the weak field limit we can rewrite Eq. (27) as \(\tilde{\rho } =\rho _{\mathrm{eff}}+3p_{\mathrm{eff}}/c^2\). In other words the Poisson’s equation, as one may expect, takes the following form

$$\begin{aligned} \nabla ^2 \varPhi =\frac{\gamma c^4}{2} \left( \rho _{\mathrm{eff}}+3\frac{p_{\mathrm{eff}}}{c^2}\right) . \end{aligned}$$
(30)

This is similar to the corresponding equation in GR, where we take into account pressure as a source for gravity (see [58] for more details).

Now before moving on to discuss the Euler equation, let us summarize the weak field limit and write the modified Poisson’s equation for two different categories, namely EMSG models with \(f_{RR}=0\) and \(f_{RR}\ne 0\). For the first case, using Eqs. (16) and (24)–(26), we arrive at

$$\begin{aligned} \nabla ^2 \varPhi =\frac{\gamma \, c^4}{2 f_R}\Big (\rho +2\alpha \rho ^2 c^2\Big ), \end{aligned}$$
(31)

and similarly for the second case, using Eqs. (16) and (23)–(26), we find a more complicated Poisson’s equation

$$\begin{aligned} \nabla ^2 \varPhi= & {} \frac{\gamma c^4}{6 f_R}\Big [4\rho +5\alpha \rho ^2 c^2\nonumber \\&-\frac{{\mathcal {M}}^2}{4\pi }\int \frac{e^{-{\mathcal {M}}|\mathbf {r} -\mathbf {r}'|}}{|\mathbf {r}-\mathbf {r}'|}\Big (\rho (\mathbf {r}') -\alpha \rho ^2(\mathbf {r}') c^2\Big )d^3 r'\Big ].\nonumber \\ \end{aligned}$$
(32)

4 Hydrodynamics equations in weak-field limit

To find the Newtonian limit of the hydrodynamics equations, one can take the covariant derivative of the field Eq. (2) as below

$$\begin{aligned}&\left( \nabla ^{\mu } f_R\right) R_{\mu \nu } {+} f_R\left( \nabla ^{\mu }R_{\mu \nu }\right) {-}\frac{1}{2}g_{\mu \nu }\left( \nabla ^{\mu }f\right) =\gamma \nabla ^{\mu } \left( T_{\mu \nu }\right) \nonumber \\&\quad +\left( \square \nabla _{\nu }-\nabla _{\nu }\square \right) f_R-\nabla ^{\mu } \left( f_Q \theta _{\mu \nu } \right) . \end{aligned}$$
(33)

To simplify the third term in the above equation, we recall that \(f(R,Q)=f_1(R)+f_2(Q)\). Therefore one can easily verify that

$$\begin{aligned} g_{\mu \nu }\nabla ^{\mu }f=g_{\mu \nu }\left( f_R\nabla ^{\mu }R+f_Q\nabla ^{\mu }Q \right) . \end{aligned}$$
(34)

Also, the fifth term can be simplified as below [59]

$$\begin{aligned} \left( \square \nabla _{\nu }-\nabla _{\nu }\square \right) f_R=R_{\mu \nu }\nabla ^{\mu }f_R. \end{aligned}$$
(35)

Using the Bianchi identity, and after some manipulations, one can find the perturbed form of Eq. (33) as

$$\begin{aligned} \nabla ^{\mu }\left( \delta T_{\mu \nu }\right) =\alpha \left( \nabla ^{\mu } \left( \delta \theta _{\mu \nu } \right) -\frac{1}{2}\eta _{\mu \nu }\nabla ^{\mu } \left( \delta Q\right) \right) , \end{aligned}$$
(36)

where \(\delta Q=\delta T_{\mu \nu }\delta T^{\mu \nu }\). Note that, the background quantities are shown without the “0” index here. To achieve the hydrodynamics equations in the Newtonian limit, one can ignore the terms containing the pressure compared with the similar terms containing the density. In fact, the pressure plays role in the relativistic situations, which are not, of course, of interest in this study.

Let us look at the order of magnitudes. What we have assumed is: firstly, as mentioned before, the pressure can be ignored comparing with the density in our background system. Secondly, the gravitational field assumed to be weak. And finally, the velocities inside the background are slow. Using a small parameter \(\epsilon \), these assumptions can read

$$\begin{aligned} \frac{p}{\rho c^2}\simeq \frac{v^2}{c^2}\simeq \frac{\varPhi }{c^2} \propto \epsilon ^2 \end{aligned}$$
(37)

On the other hand, considering the Newtonian form of the Euler’s equation, one can see that \(\partial \mathbf {v}/\partial t\simeq (\mathbf {v}\cdot \nabla )\mathbf {v}\simeq \nabla \varPhi \) and, therefore

$$\begin{aligned} \frac{\partial }{\partial t}\simeq \mathbf {v}\cdot \nabla \propto \epsilon \end{aligned}$$
(38)

Moreover, remembering the smallness of \(\alpha \), some terms containing a multiplication of \(\alpha \) and \(\epsilon \) should be treated carefully. In fact, one can consider the same order of magnitude for these parameters, and keep the terms only up to \(O(\epsilon ^2)\). It is worth mentioning that, the parameter \(\epsilon \) is only a useful gadget to track the order of magnitudes. After finding the hydrodynamics equation, one can truly assume that \(\epsilon \rightarrow 1\).

Now, keeping in mind that \(u^{\mu }=(c,\mathbf {v})\), and by finding \(\delta T_{\mu \nu }\), \(\delta \theta _{\mu \nu }\), and \(\delta Q\), one can easily decompose the Eq. (33) to the temporal and spacial components. It is worth mentioning that, during the simplification of the components, the terms containing \(\alpha \epsilon ^2\) can be ignored. Furthermore, the terms including a temporal derivative multiplied by \(\epsilon ^2\) or \(\alpha \epsilon \) can be ignored. After some manipulations one can show that the t-component of the Eq. (33), can be written as

$$\begin{aligned}&(1+\alpha \rho c^2)\frac{\partial \rho }{\partial t} +\epsilon \left\{ (1+\alpha \rho c^2)\rho \nabla \cdot \mathbf {v}\right. \nonumber \\&\quad +\left. (1+2\alpha \rho c^2)\mathbf {v}\cdot \nabla \rho \right\} =0, \end{aligned}$$
(39)

then one may simply rewrite this equation as

$$\begin{aligned} \frac{\partial \rho }{\partial t} +\nabla \cdot (\rho \mathbf {v}) =-\frac{\alpha \rho c^2}{1+\alpha \rho c^2} \mathbf {v}\cdot \nabla \rho \end{aligned}$$
(40)

this equation, is the continuity equation in the weak-field limit of the f(RQ) gravity. It is not difficult to show that, in terms of the effective quantities defined in Eqs. (28) and (29), the continuity equation in the Newtonian limit can be written as

$$\begin{aligned} \frac{\partial \rho _{\mathrm{eff}}}{\partial t}+\nabla \cdot \biggl [\left( \rho _{\mathrm{eff}}+\frac{p_{\mathrm{eff}}}{c^2}\right) \mathbf {v}\biggr ]=0. \end{aligned}$$
(41)

Now, from the spacial components of Eq. (33) and using Eq. (40) after some manipulations one can find the Euler’s equation in the weak-field limit of this theory. To do so we obtain \(\partial \rho /\partial t\) from Eq. (40) and ignore the \(O(\epsilon ^3)\) terms. Then we substitute the result into the spatial components of Eq. (33). Therefore Euler’s equation in the weak-field limit reads

$$\begin{aligned} \epsilon ^2(\mathbf {v}\cdot \nabla )\mathbf {v}+\epsilon \left( \frac{\partial \mathbf {v}}{\partial t} +\nabla \varPhi +\frac{\nabla p}{\rho }\right) +\alpha c^4\nabla \rho =0. \end{aligned}$$
(42)

This equation can be written in terms of the effective quantities to. The result is

$$\begin{aligned} \frac{\partial \mathbf {v}}{\partial t}+(\mathbf {v}\cdot \nabla )\mathbf {v} +\nabla \varPhi +\frac{\nabla p_{\mathrm{eff}}}{\rho _{\mathrm{eff}}}=0. \end{aligned}$$
(43)

Now, we have a complete set of differential equations in the weak field limit governing the dynamics of a self-gravitating fluid. Using Eqs. (40), (42), the Poisson’s Eq. (31) (or 32), and also an equation of state, one can investigate the gravitational stability of a self gravitating fluid in the context of the f(RQ) gravity.

5 Jeans analysis in the EMSG

Let us focus on a static, infinite, homogeneous, spherically symmetric fluid in the context of the EMSG. The question is when such a system can be locally fragmented under its own gravity? To find the answer, one should find the dispersion relation by linearizing the Poisson’s Eq. (31) or (32) for the cases \(f_{RR}= 0\) and \(f_{RR}\ne 0\) respectively, and also the hydrodynamics Eqs. (40) and (42). The physical quantities are considered to be as \(\mathsf {X}=\mathsf {X}_0 +\mathsf {X}_1\), where \(\mathsf {X}_1\ll \mathsf {X}_0\) and the “0” (1) index indicates the background (perturbed) quantities. For a static background system we have \(\mathbf {v}_0=0\). Moreover, homogeneity implies that, \(\rho _0\) and \(p_0\) are constant. Also, we set the gravitational potential of the background to be constant. However these assumptions do not satisfy the background equations. Therefore, to avoid the underlying ambiguity, one may assume that the Poisson’s equation can describe only the perturbed system. This assumption is known as the Jeans swindle and is widely used even in the standard Newtonian [46], and post-Newtonian [60] cases. The resulting first order equations can be easily found as follows

$$\begin{aligned}&\frac{\partial \rho _1}{\partial t}+\rho _0\nabla \cdot \mathbf {v}_1=0, \end{aligned}$$
(44)
$$\begin{aligned}&\frac{\partial \mathbf {v}_1}{\partial t}+\nabla \varPhi _1+\frac{\nabla p_1}{\rho _0} +\alpha c^4\nabla \rho _1=0, \end{aligned}$$
(45)
$$\begin{aligned}&\nabla ^2\varPhi _1=\frac{\gamma c^4}{2f_R} \left( \rho _1+4\alpha c^2\rho _0\rho _1 \right) , \end{aligned}$$
(46)

when \(f_{RR}\ne 0\) then the last equation should be replaced by

$$\begin{aligned} \nabla ^2\varPhi _1=\frac{\gamma c^4}{6f_R} \bigg [4\rho _1+10\alpha c^2\rho _0\rho _1+{\mathcal {H}}\bigg ], \end{aligned}$$
(47)

where \({\mathcal {H}}\) is defined as

$$\begin{aligned} {\mathcal {H}}=-\frac{{\mathcal {M}}^2}{4\pi }\int \frac{e^{-{\mathcal {M}} |\mathbf {r-r'}|}}{|\mathbf {r-r'}|}\left( \rho _1(\mathbf {r'})-2\alpha c^2\rho _0\rho _1 (\mathbf {r'})\right) d^3r'. \end{aligned}$$
(48)

Now, taking the temporal derivative of Eq. (44) and the divergence of Eq. (45), and also using the Poisson’s Eq. (46), one can easily find the following result for the case of \(f_{RR}=0\)

$$\begin{aligned} \frac{1}{\rho _0}\frac{\partial ^2\rho _1}{\partial t^2} -\frac{\gamma c^4(1+4\alpha c^2\rho _0)}{2f_R}\rho _1 -\left( \frac{c_s^2}{\rho _0}+\alpha c^4 \right) \nabla ^2\rho _1=0.\nonumber \\ \end{aligned}$$
(49)

Similarly for the case of \(f_{RR}\ne 0\) one may simply find

$$\begin{aligned}&\frac{1}{\rho _0}\frac{\partial ^2\rho _1}{\partial t^2} -\frac{\gamma c^4}{6f_R}\bigg [4\rho _1+10\alpha c^2\rho _0 \rho _1+{\mathcal {H}}\bigg ]\nonumber \\&\quad -\left( \frac{c_s^2}{\rho _0}+\alpha c^4 \right) \nabla ^2\rho _1=0. \end{aligned}$$
(50)

In the spherical coordinates system \((r,\theta ,\varphi )\), using the Fourier form for the first-order perturbations as \(\rho _1=\rho _a \exp \left( i\left( \mathbf {k}\cdot \mathbf {r}-\omega t \right) \right) \), one can simplify the dispersion relations. Therefore, Eq. (50) can be straightforwardly integrated setting \(\mathbf {r-r'}=\mathbf {R}\). Without loss of generality, one can assume that the wavenumber k is along with the z axis. So, after some manipulations it can be seen

$$\begin{aligned} {\mathcal {H}}= & {} -\frac{{\mathcal {M}}^2}{4\pi }(1-2\alpha c^2\rho _0)\rho _1 \int \limits _{0}^{\infty }\int \limits _{0}^{\pi }\int \limits _{0}^{2\pi } \frac{e^{-{\mathcal {M}}R}}{R}e^{-i k R\cos \theta }\nonumber \\&\times R^2\sin \theta dR d\theta d\varphi =-\frac{(1-2\alpha c^2\rho _0)}{1+\frac{k^2}{{\mathcal {M}}^2}}\rho _1. \end{aligned}$$
(51)

Finally, the dispersion relation can be found as

$$\begin{aligned}&\omega ^2 - (c_s^2+\alpha c^4\rho _0)k^2+ \frac{\gamma c^4\rho _0}{2 f_R}\nonumber \\&\quad \times {\left\{ \begin{array}{ll} (1+4\alpha c^2\rho _0)=0,~~f_{RR}=0\\ \frac{1}{3}\left( 4 +10\alpha c^2\rho _0 -\frac{(1-2\alpha c^2\rho _0){\mathcal {M}}^2}{{\mathcal {M}}^2+k^2}\right) {=}0,~~f_{RR}\ne 0 \end{array}\right. }\nonumber \\ \end{aligned}$$
(52)

Now, let us investigate these two different cases with more detail.

5.1 The case \(f_{RR}=0\)

Let us introduce two quantities that encode the modifications of the dispersion relations:

$$\begin{aligned} {\mathcal {C}}_s^2= & {} c_s^2+\alpha c^4\rho _0, \end{aligned}$$
(53)
$$\begin{aligned} {\mathcal {G}}= & {} \frac{G}{f_R}(1+4\alpha c^2\rho _0). \end{aligned}$$
(54)

Thus, the first expression in the dispersion relation (52) can be recast as

$$\begin{aligned} \omega ^2-{\mathcal {C}}_s^2 k^2+4\pi {\mathcal {G}}\rho _0=0. \end{aligned}$$
(55)

This equation is similar to its Newtonian counterpart. In other words, turning the EMSG correction terms off, the Newtonian dispersion relation will be reproduced. From Eqs. (53)–(55) it is clear that EMSG for \(f_R=1\) and \(\alpha >0\) (\(\alpha <0\)) increases (decreases) both the sound speed and the gravitational strength effectively. However, these quantities have completely opposite impacts on the stability of the system. Therefore at the first sight it is not trivial to argue the final outcome of EMSG’s corrections on the stability of the system. Nevertheless, it is important to mention that in our perturbative analysis we assumed that \(|\alpha c^2\rho |\ll 1\). Therefore, it is clear from Eq. (54) that the EMSG effects do not change the effective gravitational strength significantly. On the other hand, the effective sound speed can be influenced substantially. One should note that we are working in a regime where the characteristic energy scale is high enough to allow EMSG effects to play role. In such a situation, since the sound speed \(c_s^2\) is not necessarily much smaller than c, the correction term in (53) is not negligible. In other words, the ratio of two terms in the right hand side of (53), i.e., \(|\alpha \rho c^2|\frac{c^2}{c_s^2}\) , should not be considered as a very small ratio. Consequently, if we consider the sound speed as the representative of pressure content of the system, then one may accordingly infer that EMSG influences the effective pressure in the system. Now it make sense to conclude from (53) that if \(\alpha >0\) (\(\alpha <0\)) then EMSG stabilizes (destabilizes) the fluid. In the following we compute a new Jeans wavenumber to clarify this issue.

By setting \(\omega =0\) in Eq. (55), we can obtain the border of stability. In this case, the new Jeans wavenumber in the context of EMSG can be obtained as

$$\begin{aligned} k_{JE}^2=\frac{4\pi {\mathcal {G}}\rho _0}{f_R {\mathcal {C}}_s^2} =k_J^2f_R^{-1}\Big (\frac{1+4\alpha c^2\rho _0}{1+\frac{\alpha c^4 \rho _0}{c_s^2}}\Big ) \end{aligned}$$
(56)

where \(k_J^2=4\pi G\rho _0/c_s^2\) is the standard Jeans wavenumber.

Now, we can recast the dispersion relation in Eq. (55) in a more useful form by introducing the following variables:

$$\begin{aligned}&\tilde{\omega }= \frac{\omega }{\sqrt{4\pi G\rho _0}}, \end{aligned}$$
(57)
$$\begin{aligned}&\tilde{k}= \frac{k}{k_J}. \end{aligned}$$
(58)

Thus, after some calculations, we obtain

$$\begin{aligned} \tilde{\omega }^2 - \left( 1+ \alpha \rho _0\frac{c^4}{c_s^2}\right) \tilde{k}^2 + \frac{1+4\alpha c^2\rho _0}{f_R}=0. \end{aligned}$$
(59)

Let us note that, once again, turning off the EMSG’s correction terms the Newtonian case is obtained.

Also, the standard Jeans wavelength, and Jeans mass can be introduced as \(\lambda _J=2\pi /k_J\), and \({\mathfrak {M}}_J =\frac{4\pi }{3}\rho _0\left( \frac{\lambda _J}{2}\right) ^3\) respectively, where \({\mathfrak {M}}_J\) is defined as the mass inside a sphere with radius \(\lambda _J/2\). It can be shown that, the modified versions of the Jeans wavelength and mass respectively are

$$\begin{aligned} \lambda _{JE}^2=\lambda _J^2\,\Big (\frac{ 1+\alpha \rho _0 \frac{c^4}{c_s^2}}{1+4\alpha c^2 \rho _0}\Big )f_R. \end{aligned}$$
(60)

and

$$\begin{aligned} {\mathfrak {M}}_{JE}={\mathfrak {M}}_{J} \left( \frac{\lambda _{JE}}{\lambda _J}\right) ^3 ={\mathfrak {M}}_{J}f_R^{3/2}\Big (\frac{ 1+\alpha \rho _0 \frac{c^4}{c_s^2}}{1+4\alpha c^2 \rho _0}\Big )^{\frac{3}{2}}. \end{aligned}$$
(61)

The fluid system can be more unstable (stable) in the context of EMSG than the Newtonian case, whenever \({\mathfrak {M}}_{JE}<{\mathfrak {M}}_{J}\) (\({\mathfrak {M}}_{JE}>{\mathfrak {M}}_{J}\)). Equation (61) shows that, deviations from the standard case directly depends on the sign and value of the free parameter \(\alpha \). It is obvious that, negative (positive) values of \(\alpha \) make the system more unstable (stable) in the context of EMSG with respect to the Newtonian case. Another less interesting point is that higher values for \(f_R\) leads to higher Jeans masses. In other words, by increasing \(f_R\) the system is stabilized. This is expected since \(f_R\) reduces the effective gravitational constant, i.e., \(G_{eff}\propto G/f_R\), and consequently weakens the destabilizing behaviour of gravitational force. However, we know that this parameter cannot deviate from unity significantly.

Before moving on to close this subsection it is interesting to mention that in the original EMSG model [41], \(\alpha =-\eta <0\). On the other hand this model leads to bouncing cosmological solutions and prevents the big bang singularity. As we showed, negative \(\alpha \) destabilizes the local perturbations and supports the local gravitational collapse. This behaviour seems completely in disagreement with the “stabilizing” behaviour of the theory in the early universe. However, one should note that here we present a non-relativistic description, while in the early universe we deal with a completely relativistic situation.

5.2 The case \(f_{RR}\ne 0\)

For sake of completeness, we also compute the Jeans scale for the more general case \(f_{RR}\ne 0\). In this case, we can recast the second relation in the dispersion relation in Eq. (52) in term of the variable in Eqs. (57) and (58):

$$\begin{aligned} \tilde{\omega }^2+ & {} \frac{1}{3f_R}\biggl [ (4+10\alpha c^2\rho _0) -\frac{1-2\alpha c^2\rho _0}{1+\tilde{k}^2k_J^2{\mathcal {M}}^{-2}}\biggr ]\nonumber \\&\quad -\biggl [1+\frac{\alpha c^4\rho _0}{c_s^2}\biggr ]\tilde{k}^2=0. \end{aligned}$$
(62)

Then, one can simply find the modified version of the Jeans wavenumber by setting \(\tilde{\omega }^2=0\) in this equation, and solve for \(\tilde{k}^2\). In this case we found two solutions. One of these solutions recovers the standard Jeans wavenumber. This solution, without any expansion with respect to \(\alpha \), reads

$$\begin{aligned} k_{JE}^2&=\frac{1}{6 {\mathcal {C}}_s^2 f_R} \biggl [2 c_s^2 k_J^2 (5 \alpha c^2 \rho _0+2)-3 {\mathcal {C}}_s^2 f_R {\mathcal {M}}^2\nonumber \\&\quad \!+\!\sqrt{(c_s^2 {\mathcal {K}}_1 \!+\! 3 f_R {\mathcal {M}}^2 \alpha c^2 \rho _0)^2 \!+\!(6 c_s k_J{\mathcal {M}} {\mathcal {C}}_s)^2 f_R {\mathcal {K}}_2}\biggr ] \end{aligned}$$
(63)

where we have defined

$$\begin{aligned}&{\mathcal {K}}_1=3 f_R {\mathcal {M}}^2-2 k_J^2 (5 \alpha c^2\rho _0+2), \end{aligned}$$
(64)
$$\begin{aligned}&{\mathcal {K}}_2=\frac{{\mathcal {G}}f_R}{G}, \end{aligned}$$
(65)

and we have chosen the solution which reproduces the standard Jeans wavenumber in the limiting case \({\mathcal {M}}\rightarrow \infty \) or equivalently \(f_{RR}\rightarrow 0\). It is worth mentioning that the limit \({\mathcal {M}}\rightarrow 0\) does not recover the Newtonian results.

As our final remark in this section, one should take this case, i.e., \(f_{RR}\ne 0\), with more care. As already mentioned, our analysis in this paper can be considered as a modification to the so-called metric f(R) gravity theory. In order to make f(R) gravity suitable for explaining the cosmic speed up, it seems necessary to include a very small scalar mass \({\mathcal {M}}\). Otherwise the extra scalar degree of freedom intrinsic in f(R) gravity is not light enough to propagate in cosmic scales. Therefore it will not be effective for explaining the late time cosmic acceleration. However it is well-known in the relevant literature that if we take small \({\mathcal {M}}\), then the theory will have serious problems in the weak field limit and cannot recover the Newtonian gravity , for a review see [32]. This is exactly what we see in our calculations for the Jeans wavenumber. In other words, we see that at the limit \({\mathcal {M}}\rightarrow 0\) Eq. (63) does not recover the Newtonian Jeans wavenumber.

To address the above mentioned problem, it is necessary to take into account screening behaviour of f(R) gravity theory [61]. Investigating stability issues in the presence of screening effects, can be considered as a separate study. In this paper we continue our analysis without including the screening effects in the calculations for \(f_{RR}\ne 0\), and put emphasis on the \(f_{RR}=0\) case which does not suffer from the above mentioned problem. Therefore, hereafter we only discuss the \(f_{RR}=0\).

6 Toomre’s criterion in EMSG

So far we have studied the stability of an infinite homogeneous medium without rotation. Nevertheless, the stability of rotating systems in EMSG is interesting in the sense that in high energy systems like HMNS, where we expect that EMSG contributions to be significant, the differential rotation is one of the main ingredients of the system. For simplicity, we restrict ourselves to rotating thin disks. On the other hand, for such a system there is already a well-known stability criterion in the standard Newtonian description known as Toomre’s stability criterion [51]. In this case, one may simply compare the stability criteria based on EMSG with the Newtonian one.

To study the gravitational stability of a thin self-gravitating fluid disk, one should find the dispersion relation of propagating perturbations. This task can be addressed by particularizing the hydrodynamics equations of EMSG, given in Sect. 4, for a thin disk. One may conveniently assume that \(\rho =\varSigma \delta (z)\), where \(\varSigma \) is the surface density and \(\delta \) is the Dirac’s Delta function. Then the continuity Eq. (40) takes the following form

$$\begin{aligned} \frac{\partial \varSigma }{\partial t} +\nabla \cdot (\varSigma \mathbf {v}) =- \frac{\alpha \varSigma \delta (z) c^2}{1+\alpha \varSigma \delta (z) c^2} \mathbf {v}\cdot \nabla \varSigma \end{aligned}$$
(66)

It is clear that the right-hand-side diverges on the plane of the disk. This is not the case in Newtonian gravity. In fact, at \(z=0\) we have \(\alpha \varSigma \delta (z) c^2\rightarrow \infty \). This means that our linearized analysis based on the main assumption that \(\alpha \rho c^2\ll 1\) is obviously violated. To skip this complexity and keep the analysis self-consistent, it is useful to assume a finite thickness for the disk. To do so, we simply assume that the density does not change in the vertical direction and is given by \(\rho (r,\varphi )=\varSigma (r,\varphi )/l\) [51], where l is a small thickness of the disk and appears as a constant in our calculations. This method leads to a powerful estimation for the effect of the thickness on the stability of the disk [62]. For \(|z|> l/2\) the matter density vanishes \(\rho =0\), and the only constraint on the thickness is \(l\gg \alpha \varSigma c^2\) everywhere on the disk. This condition guarantees the validity of our perturbative analysis. For a more careful way to include the thickness of the disk, we refer the reader to [63] and [64]. Note that we are interested to the effects of EMSG and not the full analysis of the thickness. Therefore, it turns out that following the Toomre’s method [51] is helpful here. Finally for completeness, we will shortly discuss the generalization of the other method, i.e., [63], in EMSG.

Now for our disk with finite thickness, the continuity equation and the Euler’s equation read

$$\begin{aligned}&\frac{\partial \varSigma }{\partial t} +\nabla \cdot (\varSigma \mathbf {v})\simeq -\frac{\alpha c^2 \varSigma }{l} \mathbf {v}\cdot \nabla \varSigma , \end{aligned}$$
(67)
$$\begin{aligned}&\frac{\partial \mathbf {v}}{\partial t}+(\mathbf {v}\cdot \nabla )\mathbf {v}+\nabla \varPhi +\frac{\nabla p}{\varSigma }+\frac{\alpha c^4}{l}\nabla \varSigma =0. \end{aligned}$$
(68)

Note that in the Euler’s equation the p is a pressure defined as force per unit length. On the other hand the Poisson’s equation for the case \(f_{RR}=0\) is given by the following equation

$$\begin{aligned} \nabla ^2\varPhi =\frac{\gamma {c^4 \varSigma }}{2 l\,f_R} \left( 1+\frac{2\alpha {c^2}}{l}\varSigma \right) . \end{aligned}$$
(69)

In order to have a closed set of differential equations, we assumed the equation of state (EOS) to be barotropic,i.e., \(p=p(\varSigma )\). To study the stability of a disk we have to achieve the modified version of Toomre’s criterion in the context of EMSG. To do so, one should linearize the hydrodynamics Eqs. (67)–(69) in the cylindrical coordinate system \((r,\varphi ,z)\). Equations (67) and (68) in the cylindrical coordinate are

$$\begin{aligned}&\frac{\partial \varSigma }{\partial t}+\frac{1}{r}\frac{\partial }{\partial r} \left( \varSigma r v_{r}\right) +\frac{1}{r}\frac{\partial }{\partial \varphi } \left( \varSigma v_{\varphi }\right) +\frac{\partial }{\partial z}\left( \varSigma v_{z}\right) \end{aligned}$$
(70)
$$\begin{aligned}&\qquad \quad \qquad \qquad \qquad \quad \quad \quad = - \frac{\alpha c^2 \varSigma }{l} \Big (v_r \frac{\partial \varSigma }{\partial r} + \frac{v_\varphi }{r} \frac{\partial \varSigma }{\partial \varphi }\Big ),\nonumber \\ \end{aligned}$$
(71)
$$\begin{aligned}&\frac{\partial v_{r}}{\partial t}+v_{r} \frac{\partial v_{r}}{\partial r} +\frac{v_{\varphi }}{r} \frac{\partial v_{r}}{\partial \varphi } -\frac{v_{\varphi }^{2}}{r}+v_z \frac{\partial v_r}{\partial z} \end{aligned}$$
(72)
$$\begin{aligned}&\qquad \qquad \qquad \quad \quad \quad = - \frac{\partial }{\partial r}\left( \varPhi +h\right) -\frac{\alpha c^4}{l} \frac{\partial \varSigma }{\partial r}, \end{aligned}$$
(73)
$$\begin{aligned}&\frac{\partial v_{\varphi }}{\partial t}+v_{r} \frac{\partial v_{\varphi }}{\partial r} +\frac{v_{\varphi }}{r} \frac{\partial v_{\varphi }}{\partial \varphi } +\frac{v_{\varphi } v_{r}}{r}+v_z \frac{\partial v_{\varphi }}{\partial z} \end{aligned}$$
(74)
$$\begin{aligned}&\qquad \qquad \qquad \quad \quad \quad = - \frac{1}{r}\frac{\partial }{\partial \varphi }\left( \varPhi +h\right) -\frac{\alpha c^4}{l r}\frac{\partial \varSigma }{\partial \varphi }, \end{aligned}$$
(75)
$$\begin{aligned}&\frac{\partial v_{z}}{\partial t}+v_{r} \frac{\partial v_{z}}{\partial r} +\frac{v_{\varphi }}{r} \frac{\partial v_{z}}{\partial \varphi } +v_z \frac{\partial v_{z}}{\partial z} \end{aligned}$$
(76)
$$\begin{aligned}&\qquad \qquad \qquad \quad \quad \quad = -\frac{\partial }{\partial z}\left( \varPhi +h\right) -\frac{\alpha c^4}{l }\frac{\partial \varSigma }{\partial z}, \end{aligned}$$
(77)

for more details we refer the reader to [46].

Let us find the perturbed form of Eqs. (70)–(75). We recall that, the physical quantities are considered to be as \(\mathsf {X}=\mathsf {X}_0+\mathsf {X}_1\), where \(\mathsf {X}_1\ll \mathsf {X}_0\) and the “0” (1) index represents the background (perturbed) quantity. Moreover, the background is assumed to be static and axis-symmetric, and the initial radial and vertical velocities vanish everywhere throughout the disk i.e., \(v_{r0}=v_{z0}=0\). The only non-zero velocity component is \(v_{\varphi 0}\). We assume that \(v_{\varphi 0}\) is a function of radius and does not change in the z direction. Of course, one can show that such barotropic equilibrium state with the density \(\rho (r,\varphi )=\varSigma (r,\varphi )/l\) does not exists. Therefore, for the background system we are using a generalized version of the so-called Jeans-swindle. Consequently although we do not care about the validity of the background system, the linear perturbations should satisfy all the linearized equations. The linearized versions of the continuity equation and also the radial, azimuthal and the vertical components of the Euler’s equation read (the linearized Poisson’s equation is discussed in the next subsection)

$$\begin{aligned}&\frac{\partial \varSigma _{1}}{\partial t}\!+\!\frac{1}{r}\frac{\partial }{\partial r} \left( \varSigma _{0} r v_{r1}\right) \!+\!\varOmega \frac{\partial \varSigma _{1}}{\partial \varphi } \!+\!\frac{\varSigma _{0}}{r} \frac{\partial v_{\varphi 1}}{\partial \varphi } \end{aligned}$$
(78)
$$\begin{aligned}&\qquad \quad \quad +\varSigma _0\frac{\partial v_{z1}}{\partial z} \!=\!-\frac{\alpha c^2\varSigma _0}{l}\Big (v_{r1} \frac{\partial \varSigma _0}{\partial r} \!+\!\varOmega \frac{\partial \varSigma _1}{\partial \varphi }\Big ), \end{aligned}$$
(79)
$$\begin{aligned}&\frac{\partial v_{r1}}{\partial t}+\varOmega \frac{\partial v_{r1}}{\partial \varphi } -2\varOmega v_{\varphi 1}=-\frac{\partial }{\partial r}\left( \varPhi _{1}+h_{1}\right) -\frac{\alpha c^4}{l}\frac{\partial \varSigma _1}{\partial r}, \end{aligned}$$
(80)
$$\begin{aligned}&\frac{\partial v_{\varphi 1}}{\partial t}+\varOmega \frac{\partial v_{\varphi 1}}{\partial \varphi } +\frac{\kappa ^{2}}{2\varOmega } v_{r1}=-\frac{1}{r} \frac{\partial }{\partial \varphi }\left( \varPhi _{1}+h_{1}\right) -\frac{\alpha c^4}{l\,r}\frac{\partial \varSigma _1}{\partial \varphi }, \end{aligned}$$
(81)
$$\begin{aligned}&\frac{\partial v_{z 1}}{\partial t}+\varOmega \frac{\partial v_{z 1}}{\partial \varphi } =-\frac{\partial }{\partial z}\left( \varPhi _{1}+h_{1}\right) -\frac{\alpha c^4}{l} \frac{\partial \varSigma _1}{\partial z}, \end{aligned}$$
(82)

where \(\varOmega =v_{\varphi 0}/r\) is the rotational frequency and \(v_r\) and \(v_{\varphi }\) are the radial and azimuthal components of velocity respectively. Moreover, \(\kappa =\sqrt{4\varOmega ^2 +2r\varOmega \varOmega '}\) is the epicyclic frequency and also \(h_1\) is defined as \(h_1=c_s^2 \varSigma _1/\varSigma _0\). It should be noted that the prime stands for derivative with respect to r. It is worth mentioning that, by turning the correction terms (containing \(\alpha \)) off, the Newtonian hydrodynamics equations will be reproduced. Hereafter we define a new parameter \(\alpha ^*\) as

$$\begin{aligned} \alpha ^*=\frac{\alpha }{l}. \end{aligned}$$
(83)

Now, by applying the WKB approximation one can study the local stability of the disk and benefit an elegant simplification as well. The general form of the perturbed quantities can be written as \(\mathsf {X}_1=\mathsf {X}_a \exp \left( i\left( \mathbf {k}\cdot \mathbf {r} +m\varphi -\omega t \right) \right) \), where \(\omega \) is the oscillation frequency, k is the wavenumber, and m is a positive integer which determines the symmetry of the disturbances. For a tightly wound density wave, one can show that \(|k r/m|\gg 1\), therefore, using the WKB approximation the terms containing 1/r can be neglected comparing with the analogous terms containing k. This directly means that our description works only for local perturbations and one cannot use it for global stability of the disk. It is also necessary to mention that in order to recover the standard Toomre’s criterion and regarding to the small thickness of the disk, we study perturbations propagating the \(x-y\) plane. This means that the vertical component of \(\mathbf {k}\) is zero. Finally, using Eq. (78), the continuity equation can be written as follows

$$\begin{aligned} {[}\omega -m \varOmega (1+\alpha ^* c^2 \varSigma _0)]\varSigma _a-k \varSigma _0 v_{ra} +i \varSigma _0 \frac{\partial v_{az}}{\partial z}=0. \end{aligned}$$
(84)

Moreover, considering Eqs. (80) and (81), the solutions for coefficients of Fourier expansions of perturbed velocity components can be found as below

$$\begin{aligned} v_{r a}= & {} \frac{(m\varOmega -\omega )k}{\varDelta } \left( \varPhi _a+h_a+\alpha ^*{c^4}\varSigma _a \right) , \end{aligned}$$
(85)
$$\begin{aligned} v_{\varphi a}= & {} -\frac{2B i k}{\varDelta }\left( \varPhi _a+h_a +\alpha ^*{c^4}\varSigma _a \right) , \end{aligned}$$
(86)
$$\begin{aligned} v_{za}= & {} \frac{i}{m\varOmega -\omega }\frac{\partial \varPhi _a}{\partial z}, \end{aligned}$$
(87)

where rotation Oort’s constant B and \(\varDelta \) are defined as below

$$\begin{aligned}&B(r)=-\frac{1}{2}\left( \varOmega +\frac{d(\varOmega r)}{dr} \right) , \end{aligned}$$
(88)
$$\begin{aligned}&\varDelta =\kappa ^2-(m\varOmega -\omega )^2. \end{aligned}$$
(89)

In order to find the dispersion relation, the next step is to find the potential of a WKB spiral pattern \(\varSigma _1=\varSigma _a\exp \left( i(\mathbf {k}\cdot \mathbf {r}+m\varphi -\omega t) \right) \) to determine \(\varPhi _a\) in terms of \(\varSigma _a\).

6.1 The gravitational potential of a thick WKB density wave in Newtonian gravity and EMSG

Now let us calculate the gravitational potential of a WKB density wave in the context of EMSG. The linearized version of the modified Poisson’s Eq. (69) is

$$\begin{aligned} \nabla ^2 \varPhi _1=\frac{\gamma c^4}{2 l\, f_R}(1+4 \alpha ^* c^2 \varSigma _0)\varSigma _1 =\frac{4\pi {\mathcal {G}}}{l}\varSigma _1 \end{aligned}$$
(90)

where we have used the definition \({\mathcal {G}}=G/f_R (1+4\alpha ^* c^2 \varSigma _0)\) introduced in (54) and replaced \(\alpha \) with \(\alpha ^*\). Let us note that all the calculations in this subsection also hold in Newtonian gravity, we just need to set \({\mathcal {G}}\rightarrow G\). We know that the density wave \(\varSigma _1\) in the WKB approximation at arbitrary location \(\mathbf {r}\) on the disk, can be considered a plane wave propagating in the radial direction. Therefore, without loosing of generality we take \(\mathbf {k}\) along \(\hat{\mathbf {x}}\). So finding the potential of a WKB wave reduces to finding the potential of a plane density wave in EMSG. Consequently, one may write \(\varSigma _1=\varSigma _a \exp i( k x- \omega t)\). For this plane wave we guess the potential has the following functional form

$$\begin{aligned} \varPhi _1(x,z,t)= \varPhi _a e^{ i( k x- \omega t) }f(z). \end{aligned}$$
(91)

Substituting the above solution into Eq. (90), we get the following differential equation for the function f(z)

$$\begin{aligned} \frac{d^2 f}{d z^2}-k^2 f(z)=\frac{4\pi {\mathcal {G}}}{l} \frac{\varSigma _a}{\varPhi _a}. \end{aligned}$$
(92)

This equation holds for \(|z|<l/2\). On the other hand we know that for \(|z|>l/2\) the potential is given by

$$\begin{aligned} \varPhi _{out}=\varPhi _a e^{i (k x -\omega t)} e^{- k |z|}. \end{aligned}$$
(93)

Note that hereafter we restrict ourselves to trailing density waves with \(k>0\). Eq. (92) can be simply integrated to obtain f(z), and the integration constants will be fixed using the following matching conditions

$$\begin{aligned} \varPhi _1(z=\pm l/2)=\varPhi _{out}(z=\pm l/2). \end{aligned}$$
(94)

Thus, we obtain

$$\begin{aligned} \varPhi _a f(z)= & {} \frac{ \text {cosh}(k z) \text {sech}(k l/2)-1}{k^2} \frac{4\pi {\mathcal {G}}}{l}\varSigma _a \end{aligned}$$
(95)
$$\begin{aligned}&+(1-\text {tanh}(k l/2))\text {cosh}(k z)\varPhi _a. \end{aligned}$$
(96)

On the other hand, we expect that at the plane of the disk (\(z=0\)) and in the limit \(l\rightarrow 0\), the standard thin disk potential should be recovered. Therefore let us choose \(\varPhi _a\) as follows

$$\begin{aligned} \varPhi _a=-\frac{\text {sinh}\frac{kl}{2}}{\frac{kl}{2}} \frac{2 \pi {\mathcal {G}}}{k}\varSigma _a, \end{aligned}$$
(97)

it is clear that at the limit \(l\rightarrow 0\) the standard thin disk density wave, \(\varPhi _a=-\frac{2 \pi {\mathcal {G}}}{k}\varSigma _a\) [46], is recovered. Combining Eqs. (96) and (97) we find the final form of the potential as

$$\begin{aligned} \varPhi _1=-\frac{ e^{-\frac{1}{2} k (l+2 z)} \left( 2 e^{\frac{k l}{2}+k z}-e^{2 k z}-1\right) }{k l} \frac{2 \pi {\mathcal {G}} }{k}\varSigma _1, \end{aligned}$$
(98)

and at \(z=0\) we have

$$\begin{aligned} \varPhi _1=-\frac{1-e^{-\frac{k l}{2}}}{k l/2} \frac{2 \pi {\mathcal {G}}}{k}\varSigma _1=-{\mathcal {F}}(kl) \frac{2 \pi {\mathcal {G}}}{k}\varSigma _1. \end{aligned}$$
(99)

The reduction factor \({\mathcal {F}}(kl)=\frac{ 1-\exp {(-kl/2})}{k l/2}\) is exactly the coefficient derived in [51]. This reduction coefficient can be interpreted as a decrease in the surface density. Consequently, it is well-established in the literature that thickness of the disk has stabilizing effects.

6.2 The dispersion relation and the Toomre’s parameter

It is straightforward to show that the vertical average value of \(\frac{\partial \varPhi _1}{\partial z}\), namely its integration over \((-l/2,l/2)\), vanishes. Therefore, we use the approximation \(\frac{\partial \varPhi _1}{\partial z}\simeq 0\). Consequently, we have \(v_{za}\simeq 0\). This assumption is another reason for taking Toomre’s method as an estimation and not a precise calculation. Now, substituting Eqs. (99) and (85) into Eq. (84), and confining ourselves to the plane \(z=0\), we find the following dispersion relation for axisymmetric (\(m=0\)) density waves

$$\begin{aligned} \omega ^2=\kappa ^2+{\mathcal {C}}_s^2 k^2-{\mathcal {F}}(k l)2 \pi {\mathcal {G}}\varSigma _0 k, \end{aligned}$$
(100)

where the effective sound speed \({\mathcal {C}}_s^2=c_s^2+\alpha ^* c^4\varSigma _0\) is defined in (53), again replacing \(\alpha \) with \(\alpha ^*\). At the limit \(\alpha ^*\rightarrow 0\), we have \({\mathcal {C}}_s=c_s\) and \({\mathcal {G}}=G\). Therefore, as expected, the dispersion relation (100) reduces to the standard one in Newtonian gravity [62]. As we already mentioned, for \(\alpha ^*>0\), EMSG effects can be interpreted as an increase in the sound speed and in the gravitational constant as well. Increase in the sound speed, stabilizes the system while increase in the gravitational strength promotes the instability. Therefore, a careful analysis is required to discriminate between these opposite features. To do so, let us find the generalized version of the Toomre’s criterion in EMSG.

Using the dispersion relation (100), the stability condition \(\omega ^2>0\) takes the following form

$$\begin{aligned} {\mathcal {Q}}(X)^2>-\frac{4 X^2 \left( \beta +e^{-\frac{\beta }{X}}-1\right) }{\beta }, \end{aligned}$$
(101)

where the dimensionless wavelength X is defined as \(X=k_{\mathrm{crit}}/k\) and \(k_{\mathrm{crit}}=\kappa ^2/(2 \pi {\mathcal {G}}\varSigma _0)\). Furthermore the \(\beta \) parameter as the representative of the thickness of the disk is defined as \(\beta =k_{\mathrm{crit}} l/2\). Before discussing the effects of EMSG, let us briefly review the impact of thickness on the stability of the disk. Our discussion here is true in both Newtonian gravity and EMSG. It is straightforward to verify that for \(\beta \ge 1\), the right hand side of (101) gets negative. This means that all the perturbations would be stable. Note that for large \(\beta \), the characteristic length of the system in the vertical direction increases. Therefore, one may expect the ordinary Jeans’s criterion accounts for the stability of the system. At this limit Eq. (100) is written as

$$\begin{aligned} \omega ^2\simeq {\mathcal {C}}_s^2 k^2-4 \pi {\mathcal {G}}\rho +\kappa ^2, \end{aligned}$$
(102)

if we ignore the angular momentum in the system, i.e., \(\kappa =0\), then the well-known dispersion relation already derived in Jeans analysis of an infinite medium is recovered, see Eq. (55). If the combination of the last two terms in the right hand side of (102) gets positive, or equivalently if \(\beta \ge 1\), the all the wavelengths will be stable.

The other more interesting case is \(\beta \le 1\). In this case we directly use the dispersion relation (101) to find the stability criterion for each wavelength X. The boundary of stability, namely the minimum value required for \({\mathcal {Q}}(X)\) to stabilize the wavelength X is shown in Fig. 1. In this figure, darker colors show larger values of \(\beta \). Moreover, the inner surface of each curve supposed to be the unstable area. Therefore, it is clear that, the larger values of \(\beta \) decrease the instability area. In other words, this figure directly shows that increasing the \(\beta \) parameter, increases the stability of the disk. As we mentioned, from this perspective, both Newtonian and EMSG behave in a similar way.

So far we considered the stabilizing effects of the disk thickness. Now, let us investigate our main purpose in this section: impact of EMSG on the stability of self-gravitating disks. To do so, one should note that \({\mathcal {Q}}\) in (101) has been written in terms of effective parameters \({\mathcal {C}}_s\) and \({\mathcal {G}}\). Furthermore the epicycle frequency \(\kappa \) is different from Newtonian case in the sense that it includes EMSG corrections. Therefore for comparison with Newtonian gravity, it is helpful to rewrite (101) in terms of the Newtonian Toomre’s parameter \(Q=\kappa _N c_s/\pi G\varSigma _0\). Where \(\kappa _N\) is the epicycle frequency obtained using the Newtonian gravitational force. For a given matter density \(\varSigma _0\), let us express the epicycle frequency as

$$\begin{aligned} \kappa ^2=\kappa _N^2+\delta \kappa ^2, \end{aligned}$$
(103)

where \(\delta \kappa ^2\) is the corrections to \(\kappa _N^2\) induced by EMSG. We expect this correction be proportional to \(\alpha ^*\). Accordingly we have \(X=X_N+\delta X\) and \(\beta =\beta _N +\delta \beta \), where \(\delta X\) and \(\delta \beta \) are also proportional to \(\alpha ^*\). Now, we rewrite Eq. (101) as

$$\begin{aligned} Q(X_N)^2>-\frac{4 X_N^2 \left( \beta _N f_R +e^{-\frac{\beta _N }{X_N}}-1\right) }{\beta _N f_R}+\alpha ^* \varDelta +{\mathcal {O}}(\alpha ^{*2}), \end{aligned}$$
(104)

where \(\varDelta \) is a complicated function of \(\delta \beta \), \(\delta \kappa \), \(\delta X\), \(\beta _N\), \(X_N\) and \(\kappa _N\). So we avoid to write it here. This term includes all the corrections introduced by EMSG. It is difficult to specify the sign of \(\alpha ^* \varDelta \). We know that if \(\alpha ^* \varDelta <0\) (\(>0\)) then EMSG stabilizes (destabilizes) the disk. Furthermore we need a known surface density \(\varSigma _0\) to calculate all the functions in Eq. (104) and quantify the differences of EMSG and Newtonian gravity. In the next section we study an exponential toy model in order to describe the EMSG impact on the stability of disks.

Fig. 1
figure 1

Curves from up to down belongs to \(\beta = 0.01\), 0.3, 0.6, 0.8, 0.96. This shows that thickness of the disk seriously stabilizes the disk. Of course, one should note that we have used an estimative way in our analysis

Before closing this section it is worthy to mention that we used an estimative method to find the Toomre’s criterion. Interestingly we found that the reduction factor \({\mathcal {F}}\) in the dispersion relation is the same as in Newtonian gravity. On the other hand, we know that a more precise method to include the thickness of the disk in Newtonian gravity leads to a different reduction factor as \({\mathcal {F}}=(1+ kl/2)^{-1}\) [63]. Based on what happened in out estimative method, one may expect that the above mentioned reduction factor appears in EMSG as well. In this case the dispersion relation may be written as

$$\begin{aligned} \omega ^2=\kappa ^2+{\mathcal {C}}_s^2 k^2- 2\pi {\mathcal {G}} \varSigma _0\frac{k}{1+ k l/2}. \end{aligned}$$
(105)

However, we continue working with the estimative method explained comprehensively in this section.

7 Exponential fluid disk in the context of the EMSG

Here, we are going to achieve a modified version of the Toomre’s criterion using a common toy model. In fact, this model could help us to compare the results in the context of EMSG and Newtonian gravity. The EMSG effects appear in the frequency parameter \(\kappa \) and effective parameters \({\mathcal {C}}_s\) and \({\mathcal {G}}\). Consequently, as mentioned earlier, it seems that, it is not straightforward to compare the new criterion with the standard one. However, by specifying the surface density profile \(\varSigma _0\), one can compare both theories.

The exponential surface density profile is widely used to model wide variety of astrophysical systems. In this section we take the following exponential model as a toy model to clarify some differences between EMSG and GR in the weak field limit

$$\begin{aligned} \varSigma _0=\varSigma _p e^{-2y}. \end{aligned}$$
(106)

Here, \(y=r/2R_d\) is a dimensionless radius and \(\varSigma _p\) and \(R_d\) are the central density and a characteristic length scale respectively. Taking such a density profile, by solving the Poisson’s equation in the Newtonian regime, one can show that the gravitational potential of a razor thin disk reads

$$\begin{aligned} \varPhi _0(y,z=0)=-2\pi G \varSigma _p R_d y\left[ I_0(y)K_1(y)-I_1(y)K_0(y) \right] \end{aligned}$$
(107)

where \(I_n\) and \(K_n\) for \(n=0,1\), are modified Bessel functions of the first and second kinds, respectively (For more details see [46]). We need to find the gravitational potential of a disk with a small thickness in the context of EMSG. It is clear that, Eq. (69) can be written as

$$\begin{aligned} \nabla ^2\varPhi _0=\frac{4\pi G}{l f_R}\left( \varSigma _0+2\alpha ^* c^2\varSigma _0^2 \right) . \end{aligned}$$
(108)

Note that, we will assume \(f_R=1\) in the following. However, it can be recovered in the results by replacing G by \(G/f_R\). It is clear from Eq. (108) that in order to find the potential in EMSG, we can simply add the Newtonian potentials of two separate disks with small thicknesses l, and mass densities \(\varSigma _0/l\) and \(2\alpha ^* c^2\varSigma _0^2/l\). Therefore all we need is to find the gravitational potential of a thick disk with exponential functionality in the radial direction in Newtonian gravity. To do so, we find the gravitational potential of a thin disk along the z axis and then integrate over thin disks to find the gravitational potential of a thick disk. Let us begin with finding the Newtonian gravitational potential of a razor thin and exponential disk, with the density profile given by Eq. (106). The gravitational potential for the field points that situated on the z axis could be found as

$$\begin{aligned} \varPhi _{\mathrm{out}}(r=0,z)= & {} -G\int \frac{\varSigma _0 dA}{|\mathbf {r}-\mathbf {r}'|}\nonumber \\= & {} -G\varSigma _p\int _{0}^{\infty }\int _{0}^{2\pi } \frac{e^{-r/R_d}r'dr'd\varphi }{\sqrt{r^{'2}+z^2}}, \end{aligned}$$
(109)

the integral can be simply solved, see [65], to give

$$\begin{aligned} \varPhi _{\mathrm{out}}(y=0,z)=\pi ^2 G \varSigma _p \left| z\right| \left( \mathbf {H}_{-1}\left( \frac{\left| z\right| }{R_d}\right) +Y_1\left( \frac{\left| z\right| }{R_d}\right) \right) , \end{aligned}$$
(110)

where, \(\mathbf {H}\) and Y denotes the Struve function and the Bessel function of the second kind respectively. Note that, both of these functions are well-behaved. On the other hand, the gravitational potential on the surface of the disk plane is given by Eq. (107), i.e., \(\varPhi _{\mathrm{in}}(y,z=0)=\varPhi _0(y, z=0)\). It is worth mentioning that, since the potential \(\varPhi \) is a continues function, it is easy to show that

$$\begin{aligned} \lim \limits _{y\rightarrow 0}\varPhi _{\mathrm{in}}(y,z=0) =\lim \limits _{z\rightarrow 0}\varPhi _{\mathrm{out}}(y=0,z)=-2\pi G\varSigma _p R_d. \end{aligned}$$
(111)

Now, the gravitational potential of a razor thin disk all over the space could be found combining Eqs. (110) and (111) obtaining

$$\begin{aligned} \varPhi _{\mathrm{out}}(y=0,z)&=-2\pi G \varSigma _p R_d\nonumber \\&\quad \times \left( -\frac{\pi \left| z\right| }{2R_d} \right) \left( \mathbf {H}_{-1}\left( \frac{\left| z\right| }{R_d}\right) +Y_1\left( \frac{\left| z\right| }{R_d}\right) \right) , \end{aligned}$$
(112)

where, the first term at the right hand side, denotes \(\varPhi _{\mathrm{in}}\) at \(y\rightarrow 0\) limit. Therefore, regarding Eqs. (111) and (112) and keeping in mind that the gravitational potential can be separated in terms of vertical and radial coordinates as \(\varPhi (y,z)=f_1(z)f_2(y)\) , one can find the gravitational potential of a thin disk over the whole space as follows

$$\begin{aligned} \varPhi (y,z)= & {} -2\pi G\varSigma _p R_d\left( y\left[ I_0(y)K_1(y) -I_1(y)K_0(y) \right] \right) \nonumber \\&\times \left( -\frac{\pi \left| z\right| }{2R_d} \right) \left( \mathbf {H}_{-1}\left( \frac{\left| z\right| }{R_d}\right) +Y_1\left( \frac{\left| z\right| }{R_d}\right) \right) .\nonumber \\ \end{aligned}$$
(113)

It can be shown that for the \(z\rightarrow 0\) limit, the potential of Eq. (107) will be reproduced.

In the next step, we are going to describe the calculation of the gravitational potential of a thick disk. Here, the density profile of the thick disk assumed to be

$$\begin{aligned} \rho _0(y,z)=\varSigma _0(y)\zeta (z), \end{aligned}$$
(114)

where \(\varSigma _0(r)\) is an exponential function as in (106) and \(\zeta (z)=1/l\). We have chosen this special form for \(\zeta (z)\) to be completely self-consistent with our calculations in the previous section. Of course one can use more realistic functions like \(\zeta (z)\propto e^{-\mu z}\).

As already mentioned, a thick disk can be considered as a set of many infinitesimal thin layers with thicknesses \(dz'\). The potential of each layer that situated at the vertical distance \(dz'\) from the disk, at the field point (yz) reads

$$\begin{aligned} d\tilde{\varPhi }_0(y,z)=dz' \varPhi (y,z-z')\zeta (z'). \end{aligned}$$
(115)

By adding the contributions of all layers, the result will be

$$\begin{aligned} \tilde{\varPhi }_0(y,z)=\int _{-\infty }^{\infty }dz' \varPhi (y,z-z')\zeta (z'). \end{aligned}$$
(116)

It should be noted that, since we are interested in the gravitational effects inside the disk, let us restrict ourselves to the equatorial plane \(z=0\). So, using Eqs. (113) and (116), one can see

$$\begin{aligned} \tilde{\varPhi }_0(y,z=0)&= \int _{-l/2}^{+l/2} \frac{dz'}{l}(-2\pi G\varSigma _p R_d)\nonumber \\&\quad \times \left( y\left[ I_0(y)K_1(y) -I_1(y)K_0(y) \right] \right) \nonumber \\&\quad \times \left( -\frac{\pi \left| z'\right| }{2R_d} \right) \left( \mathbf {H}_{-1}\left( \frac{\left| z'\right| }{R_d}\right) +Y_1\left( \frac{\left| z'\right| }{R_d}\right) \right) . \end{aligned}$$
(117)

Finally, after some manipulations, the integral can be analytically solved. The result reads

$$\begin{aligned} \tilde{\varPhi }_0(y)= & {} \frac{\pi c^2 \eta y (I_1(y) K_0(y)-I_0(y) K_1(y))}{2 \xi }\nonumber \\&\times \left( 8 \pi G_{1,3}^{2,0} \left( \frac{\xi ^2}{16}\bigg |\begin{array}{c} 1 \\ \frac{1}{2},\frac{3}{2},0 \\ \end{array}\right) \right. \nonumber \\&\left. -\xi ^2 \,_2F_3\left( 1,1;\frac{1}{2}, \frac{3}{2},2;-\frac{\xi ^2}{16}\right) \right) , \end{aligned}$$
(118)

where, the functions G and F that appears in this equation are the MeijerG and the generalized hypergeometric functions respectively. Also, the dimensionless constants \(\eta \) and \(\xi \) are defined as

$$\begin{aligned} \eta =G R_d\varSigma _p/c^2,~~~ \xi =l/R_d. \end{aligned}$$
(119)

It should be noted that, since the density falls off much faster along the z axis than in the radial direction within the plane, the galactic disks assumed to be thin [46]. Although we have not restricted our analysis to galactic disks, we will keep the thin disk approximation in the subsequent sections. So, the Eq. (118) could be expanded for the small values of \(\xi \). By keeping only the linear terms of \(\xi \), the result reads

$$\begin{aligned} \tilde{\varPhi }_0(y) \simeq -\frac{1}{2} \pi c^2 (\xi -4) \eta y (I_1(y) K_0(y)-I_0(y) K_1(y)).\nonumber \\ \end{aligned}$$
(120)

Again, as expected, it is clear that, turning off the contribution of thickness, the potential of a razor thin disk will be reproduced (see the Eq. 107).

Now, to find the complete form of the gravitational potential in the context of EMSG, one can apply the following replacements to Eq. (120),

$$\begin{aligned} y\rightarrow 2y,~~ R_d\rightarrow \frac{R_d}{2},~~ \varSigma _p\rightarrow 2\alpha ^* c^2\varSigma _p^2. \end{aligned}$$
(121)

The overall result is as follows

$$\begin{aligned} \tilde{\varPhi }_0(y)&=\frac{ \pi c^2 \eta y}{2} \big (-4 {\mathcal {A}} (\xi -2) (I_1(2 y) K_0(2 y)\nonumber \\&\quad -I_0(2 y) K_1(2 y))-(\xi -4)(I_1(y) K_0(y)\nonumber \\&\quad -I_0(y) K_1(y))\big ) \end{aligned}$$
(122)

where \({\mathcal {A}}\) is a dimensionless parameter defined as

$$\begin{aligned} {\mathcal {A}}=\alpha ^* \varSigma _p c^2. \end{aligned}$$
(123)

Hereafter we remove the tilde sign over the potential \(\varPhi _0\).

Now, as the next step, one can find the epicycle frequency \(\kappa \). Using the radial component of Eq. (68), i.e., Eq. (73), it can be shown that, the rotational frequency reads

$$\begin{aligned} \varOmega ^2=\frac{1}{4R_d^2 y}\left( \frac{\partial \varPhi _0}{\partial y} +\frac{1}{\varSigma _0}\frac{\partial p_0}{\partial y}+\alpha ^* c^4 \frac{\partial \varSigma _0}{\partial y} \right) . \end{aligned}$$
(124)

Using this equation one can find an analytic expression for the epicycle frequency in terms of radius. Now, following notation of the previous section, the epicycle frequency can be written as

$$\begin{aligned} \kappa ^2=\kappa _N^2+\delta \kappa ^2 \end{aligned}$$
(125)

where \(\kappa _{N}^2\) is the Newtonian part of the epicycle frequency and \(\delta \kappa ^2\) is defined to parameterize the EMSG corrections. These functions are given by the following expressions

$$\begin{aligned} \kappa _N^2&=\frac{c^2}{2 R_d^2} \bigg (\frac{\varGamma \mu e^{-2 (\varGamma -1) y} (2 (\varGamma -1) y-3)}{y}+\pi \eta (\xi -4) \nonumber \\&\quad \times ((y I_0(y){+}I_1(y)) K_1(y){-}(2 I_0(y){+}y I_1(y)) K_0(y))\bigg ) \end{aligned}$$
(126)
$$\begin{aligned} \delta \kappa ^2&= \frac{c^2{\mathcal {A}}}{2 R_d^2}\bigg ( 8 \pi \eta (\xi -2) ((2 y I_0(2 y)+I_1(2 y)) K_1(2 y) \nonumber \\&\quad -2 (I_0(2 y)+y I_1(2 y)) K_0(2 y))+\frac{e^{-2 y} (2 y-3)}{y}\bigg ) \end{aligned}$$
(127)

The new dimensionless parameter \(\mu \) is a representative of the sound speed and defined as follows

$$\begin{aligned} \mu =\frac{K\varSigma _p^{\varGamma -1}}{c^2}. \end{aligned}$$
(128)

In fact, the nature of this definition could be explored by picking an EOS. Here we have used the polytropic EOS

$$\begin{aligned} p=K\varSigma _0^{\varGamma }, \end{aligned}$$
(129)

where \(\varGamma \) is the polytropic index. Then, it is straightforward to show that the sound speed could be written as

$$\begin{aligned} c_s^2= K\varSigma _p^{\varGamma -1}\varGamma e^{-2y(\varGamma -1)}=c^2\mu \varGamma e^{-2y(\varGamma -1)}. \end{aligned}$$
(130)

Now we are in a position to define the modified version of the Toomre’s criterion.

7.1 Toomre’s criterion for an exponential disk in EMSG

First, we want to emphasis that the main aim of this paper is to study the role of EMSG. On the other hand, the effect of thickness has been studied in Sect. 6.2 in order to overcome some technical difficulties. So, assuming \(l=\xi R_d\), where \(\xi \ll 1\), one can easily expand the dispersion relation of Eq. (100) up to \(O(\xi )\). Note that, although this simplification may not include all the real physical properties, it can provide a way to track the footprints of the EMSG effects. Considering this point, the Eq. (100 could be written as following

$$\begin{aligned} \omega ^2=\kappa ^2+\left( {\mathcal {C}}_s^2 +\frac{\pi {\mathcal {G}}R_d\xi \varSigma _0}{2} \right) k^2-2\pi {\mathcal {G}}\varSigma _0 k+O(\xi ^2). \end{aligned}$$
(131)

One can see that, replacing the effective quantities \({\mathcal {C}}_s\), and \({\mathcal {G}}\), with the Newtonian values, and also ignoring the thickness l (or equivalently setting \(\xi =0\)), the Newtonian dispersion relation will be reproduced. As we already mentioned, disk thickness has stabilizing effects. This fact is clearly seen in the dispersion relation (131where the thickness parameter \(\xi \) appears with a positive sign on the right hand side. Of course one should take this description more carefully in the sense that thickness also appears in \(\kappa \). On the right hand side (hereafter RHS) of Eq. (131) all the quantities are real, therefore \(\omega ^2\) will be real too. So, if \(\omega ^2>0\), then \(\omega \) is real and the disk is stable. On the other hand, there will be some growing modes, if \(\omega ^2<0\) and therefore, the disk will be unstable. Since the RHS of Eq. (131) is a quadratic function of k, and also, the coefficient of \(k^2\) is positive (note that the correction terms assumed to be smaller than the main terms), one can seek for a k at which the RHS is minimum. This wavenumber reads

$$\begin{aligned} k_{\mathrm{min}}=\frac{\pi {\mathcal {G}} \varSigma _0}{{\mathcal {C}}_s^2} -\frac{\pi ^2 {\mathcal {G}}^2 \xi R_d \varSigma _0^2}{2 {\mathcal {C}}_s^4}. \end{aligned}$$
(132)

If the RHS of the Eq. (131) is positive for \(k_{\mathrm{min}}\), it will be positive for all other wavenumbers. Now we substitute \(k_{\mathrm{min}}\) from Eq. (132) into the RHS of Eq. (131), and expand the result for small values of \(\xi \). Furthermore, considering the definition of the modified Toomre’s parameter \({\mathcal {Q}}={\mathcal {C}}_s \kappa /\pi G\varSigma _0\), one can write the stability condition \(\omega ^2>0\) as follows

$$\begin{aligned} 1-\frac{1}{{\mathcal {Q}}^2}+\frac{\xi R_d\kappa ^2}{2\pi {\mathcal {G}} \varSigma _0{\mathcal {Q}}^4}>0. \end{aligned}$$
(133)

This inequality can be rewritten as a condition on the magnitude of \({\mathcal {Q}}\)

$$\begin{aligned} {\mathcal {Q}}(={\mathcal {Q}}_N+{\mathcal {A}}{\mathcal {Q}}_c)>1 -\frac{\kappa ^2 \xi R_d}{4\pi {\mathcal {G}}\varSigma _0}, \end{aligned}$$
(134)

where, \({\mathcal {Q}}_N\) (\({\mathcal {Q}}_c\)) is the Newtonian (correction) part of the modified Toomre’s parameter.Footnote 1 It is clear that, removing the small corrections appeared as the coefficients of \(\xi \) and \({\mathcal {A}}\), leads to the standard Toomre’s criterion. Now, regarding to the definitions \({\mathcal {C}}_s^2=c_s^2+\alpha ^* c^4\varSigma _0\), and \({\mathcal {G}}=G (1+4\alpha ^* c^2 \varSigma _0)\) , and using Eq. (125), one can express the new Toomre’s criterion (134) in terms of the Newtonian quantities. As mentioned before, holding this inequality guaranties the stability of the fluid disk against all unstable modes. In fact, the expanded form of this relation is rather complicated to be written here. However it has been written in the appendix A. Finally, the Toomre’s criterion in the context of EMSG reads

$$\begin{aligned} {\mathcal {Q}}_N>1+{\mathcal {C}}\left( \varGamma ,{\mathcal {A}},\eta ,\mu ,y,\xi \right) . \end{aligned}$$
(135)

Furthermore, \({\mathcal {C}}\left( \varGamma ,{\mathcal {A}},\eta ,\mu ,y,\xi \right) \) includes all the corrections introduced by EMSG to the local stability criterion of the rotating gaseous disks. This is one of the most important results of this paper.

Using Eq. (135), one can compare the stability of the fluid disk in the context of Newtonian gravity and EMSG. In fact, whenever \({\mathcal {C}}<0\) (\({\mathcal {C}}>0\)), the fluid disk will be more stable (unstable) in the context of EMSG. To see this fact more clearly we refer to the Fig. 2. It should be noted that, assuming \(1+{\mathcal {C}}\left( \varGamma ,{\mathcal {A}},\eta , \mu ,y,\xi \right) =RHS\), the stability condition (135) can be written as \({\mathcal {Q}}_N-RHS>0\). Therefore, in summary, to compare the Newtonian gravity and the EMSG, we will study the signature of \({\mathcal {C}}\). Moreover, the pure effects of the parameters in the context of EMSG should be tracked using the stability condition \({\mathcal {Q}}_N-RHS>0\). These situations are plotted in Fig. 2. The left panel at the first row shows the effect of thickness parameter \(\xi \) on the correction term \({\mathcal {C}}\). In fact, for a given \(\mu \) and \(\eta \) parameters, increasing \(\xi \) makes this term bigger (but with a minus sign), and so the RHS of Eq. (135) will be smaller. Therefore, one can see that, increasing \(\xi \) leads to stability in system. This is in harmony with the mentioned role for \(\beta \) in previous section. Moreover, the left panel at the second row confirms this role as well. In this panel we have shown \({\mathcal {Q}}_N-RHS\) at radius \(y=2\) for fixed values of stability parameters \(\eta \) and \(\mu \). This panel shows that increasing \(\xi \), supports the stability of the system.

Fig. 2
figure 2

The correction term \({\mathcal {C}}\) and the Toomre’s criterion in EMSG for various values of the model parameters. Note that, the Toomre’s criterion of Eq. (135) is written here as \({\mathcal {Q}}_N-RHS\). Whenever this expression is positive, the stability condition will be held. Furthermore, we take \(\varGamma =2\)

Fig. 3
figure 3

Growth rate of small perturbations in the fluid disk in the context of EMSG. It is worth mentioning that, in all panels of growth rate figures, the solid (dashed) curves show \(y=1\) (\(y=3\)) and we have assumed, \(\varGamma =2\)

To study the effect of \({\mathcal {A}}\), or equivalently the free parameter \(\alpha \) of EMSG, one can see the top right panel in Fig. 2. This panel shows that for a positive (negative) value of this parameter, increasing the magnitude of \({\mathcal {A}}\) makes the disk more stable (unstable). It is clear that, in this situation, the RHS of Eq. (135) will be smaller (bigger) and the Toomre’s criterion will be supported (opposed). This behavior is completely consistent with that explained in Sect. 5.1 for an infinite medium. As we saw in the Jeans analysis, a positive (negative) \(\alpha \) stabilizes (destabilizes) an infinite non-rotating fluid medium. This fact also can be seen in the left panel at the second row of Fig. 2 where we investigated the stability at a fixed radius \(y=2\).

It could be also interesting to investigate the role of each parameters \(\eta \) and \(\mu \) here. As mentioned before, considering Eq. (130), it is clear that, \(\mu \) shows the strength of the pressure in fluid disk. Therefore, in general, it is expected for this parameter to induce stabilizing effects. For another parameter \(\eta \), the situation is more complicated. Let us begin our discussion with looking at Eq. (126). In some astrophysical systems, in Newtonian regime, where the sound speed is much smaller than the angular velocity \(v_{\varphi }=\sqrt{r d\varPhi _0/dr}\), the first term in this equation can be ignored. However, it should be noted that, even in the Newtonian viewpoint there are some astrophysical systems like advection-dominated accretion flows [66], where can have \(c_s\simeq v_{\varphi }\). So, the epicycle frequency in this case has a coefficient of \(\eta \), and the Toomre’s parameter is proportional to \(\sqrt{\mu /\eta }\). Therefore, since \(\mu \) has stabilizing effect, it seems that \(\eta \), with a destabilizing effect, can be considered to be a representative for the gravity in system. Therefore, these parameters are useful dimensionless quantities to interpret the results and simplify the stability analysis.

The right panel at the second row of Fig. 2 devoted to studying the role of \(\eta \) and \(\mu \). This panel shows that, increasing \(\eta \) makes the disk more unstable. It may be expected, because regarding Eq. (122), increasing this parameter supports the gravitational strength. Moreover, this panel shows that, increasing \(\mu \) makes the disk more stable. Again, considering Eq. (130), it seems that \(\mu \) is a parameter that measures the strength of the pressure in system. Therefore, in general, it is natural to expect such role here. It should be noted that, regarding Eq. (126), for \(y<3/2(\varGamma -1)\), the first term could be negative. Therefore, at these radii, a high value of \(\mu \) can decrease the magnitude of \(\kappa _{N}\). As a result, increasing \(\mu \), may destabilize the inner radii. This unexpected behaviour of the pressure can be also seen in [60] and [67]. Therefore, although these two parameters have been introduced to characterize the role of gravity and pressure, there may be some exceptions that should be treated carefully.

For the sake of completeness, let us study the growth rate of the axisymmetric unstable modes. It is not difficult to rewrite Eq. (135) as

$$\begin{aligned} {\mathcal {S}}^2=-q^2(1+{\mathcal {H}})+\frac{2}{{\mathcal {Q}}}q-1, \end{aligned}$$
(136)

where \({\mathcal {S}}=i\omega /\kappa \), \(q={\mathcal {C}}_s k/\kappa \), and \({\mathcal {H}}=\pi {\mathcal {G}}\varSigma _0 \xi R_d/2{\mathcal {C}}_s^2\). The role of each parameter in the growth rate of the unstable modes for both theories EMSG and Newtonian gravity can be seen in Figs. 3, 4 and 5. To have a complete study, all cases are plotted at two different radii. The role of the theory’s free parameter \({\mathcal {A}}\) can be found in two top panels of Fig. 3. It is clear that, increasing \({\mathcal {A}}\) decrease the growth rate and consequently stabilizes the disk. This behavior is in harmony with those obtained from Fig. 2. Also, the bottom panels of this figure show that, for both positive and negative values of \({\mathcal {A}}\), the disk will be stable with the thickness parameter \(\xi \). Moreover, the role of \(\eta \) and \(\mu \) in the context of EMSG (Newtonian gravity) have been illustrated in the Fig. 4 (Fig. 5). These figures show that, increasing \(\eta \) increases the growth rates. In other words, as one may expected, this parameter makes the disk more unstable in both EMSG and Newtonian gravity. Also, these figures show that, for both theories, increasing \(\mu \), decreases the growth rates and therefore makes the disk more stable. Note that, regarding the mentioned stabilizing role of \(\mu \), this behavior was expected. Finally, it seems that, the growth rates are higher at the small radii for both theories. In some senses, it is expected. In fact at central parts we expect high surface density and consequently stronger gravity. Naturally, stronger gravity leads to higher growth rate.

Fig. 4
figure 4

Growth rate of small perturbations in the fluid disk in the context of EMSG

Fig. 5
figure 5

Growth rate of small perturbations in the fluid disk in the context of Newtonian gravity

8 Applying the results to an astrophysical system

As a realistic system, the new Toomre’s criterion can be studied in HMNSs. An HMNS is a resulting object in the merging of a neutron star binary. Because of including the strong gravitational field, and also fast movements, this system seems to be a good candidate to track the footprints of the relativistic effects in local fragmentation. However, the physical properties of HMNSs are not well known yet because of their complicated evolution. On the other hand, there are many attempts to study the main properties of HMNSs using numerical simulations in GR and also approximative methods (for example see [68] and [69]).

In this section, using a toy model recently introduced in [60] and [54], we try to roughly estimate the possibility of local fragmentation in an HMNS. In these studies, the stability parameters \(\eta \) and \(\mu \) are found for an HMNS using five models of the EOS studied in [68]. In fact, using Table 1 of [60] the gravitational local stability has been studied. Here, we focus on the model GNH3-M125 only. However, the behavior of the others are more or less similar. First, let us study the case \({\mathcal {A}}>0\). The Toomre’s criterion in the context of EMSG for this case is shown in the top side of Fig. 6. It is clear that, for a positive value of \({\mathcal {A}}\), the system will be stable at outer radii. Furthermore, increasing the value of \({\mathcal {A}}\) in this case leads to a more stable system. In this case, the edge of disk could be the most stable part of the system for a given \({\mathcal {A}}\). On the other hand, as one can see in the bottom side of Fig. 6, for negative values of \({\mathcal {A}}\), it is possible for (almost all) the system to be unstable. It should be noted that, the Toomre’s criterion by itself is not enough to conclude about the occurrence of the instability. In fact, comparing the dynamical time scales of the system can shed some light on the stability problem. To ensure about the possibility of the occurrence of the instability in a fluid system, one can compare the time scale for the perturbation growth, i.e., \(t\propto 1/|\omega |\), with the dynamical time scale of the system (for HMNS it is typically around a millisecond). It is not difficult to show that, for the unstable area in the bottom side of Fig. 6, the perturbation growth timescale for the most unstable wave-number (see Eq. 132) and the values \(R_d=4.105\) (for GNH3-M125), \(\xi =0.1\), is \(\sim 10^{-4} - 10^{-5}\) s. Therefore, the perturbation growth timescale is smaller than the dynamical timescale of an HMNS (\(\sim 10^{-3}\) s). It means that, it may be possible to occur the local instability in an HMNS system in the context of EMSG.

Fig. 6
figure 6

The Toomre’s criterion (TC), \({\mathcal {Q}}-RHS\) versus the dimensionless radius y and the model parameter \({\mathcal {A}}\) for GNH3-M125 (for more details see the Table 1 of [60]). Here \(\varGamma =2\)

Although this is a straightforward and simple outcome of our stability analysis, it can put a serious constraint on the viability of EMSG. To the best of our knowledge, no observational evidence has been reported on the existence of local fragmentation in HMNS. On the other hand EMSG with negative \(\alpha \) predicts gravitational instability in this system. Naturally, this means that EMSG with positive \(\alpha \) is more acceptable from physical point of view.

9 Discussion and conclusion

In this paper we studied the local gravitational stability of an infinite fluid (Jeans analysis) and also a differentially rotating fluid disk (Toomre’s criterion) in the context of EMSG. Firstly, by introducing the field equations of the EMSG and finding the weak field limit of this theory, we derived the modified version of the Poisson’s equation. Although, two different cases \(f_{RR}\ne 0\) and \(f_{RR}=0\) can be studied, we only focused on the latter case for practical aims. In fact, the case \(f_{RR}\ne 0\) is totally reminiscent of the weak field limit of the f(R) gravity. As we already mentioned, in this case because of an inherent non-linearity and consequently screening effects, one cannot simply linearize the field equations, for more details see [61]. We left the Jeans analysis of this case as a subject for another separate study.

By deriving the hydrodynamics equations and assuming a polytropic EOS we studied the local gravitational stability. An infinite homogeneous self gravitating fluid is the first system which is studied. In this case, by linearizing the hydrodynamics equations as well as the Poisson’s equation we found the dispersion relation, and by setting \(\omega ^2=0\) the Jeans wavenumber could be derived. By achieving the Jeans mass, we showed that the EMSG could have a stabilizing (destabilizing) effect for a positive (negative) value of the model parameter \(\alpha \). Moreover, increasing \(\alpha \) makes the system more stable for both cases \(\alpha >0\) and \(\alpha <0\).

Afterwards, we considered a fluid disk. Also, to skip some complexities and keep the analysis self-consistent, we have assumed a finite thickness for the disk. Then, by achieving the potential of a WKB density wave, and also using a perturbative method, we derived the dispersion relation. Finally, defining the modified versions of the sound speed, the gravitational constant, and the epicyclic frequency, the so-called Toomre’s criterion is achieved in the context of EMSG. Then, considering an exponential surface density profile, and also dimensionless parameters \(\eta \) (related to strength of gravity), \(\mu \) (related to the pressure in system), and \({\mathcal {A}}\) (dimensionless model parameter), the modified version of Toomre’s criterion can be rewritten in terms of the standard case. It is interesting that, the general form of this criterion could be written as \({\mathcal {Q}}_N>1+{\mathcal {C}}\left( \varGamma ,{\mathcal {A}},\eta ,\mu ,y,\xi \right) \), where an additional correction term is included here. Again, the EMSG may stabilize or destabilize the disk depending on the sign of the model parameter. However, in both cases, increasing \(\mathcal A\) will support the stability of the system. To conduct a more detailed stability analysis, we studied the rate of growing unstable modes in the disk. We showed that, for both cases \({\mathcal {A}}>0\) and \({\mathcal {A}}<0\), increasing \({\mathcal {A}}\) will makes the disk more stable. Moreover, the growth rate decreases with radius in both EMSG and Newtonian gravity.

In the last part, using a toy model which has been introduced in [54] and [60], we applied our results to an HMNS. We showed that, for a negative value of \({\mathcal {A}}\) the local fragmentation could be possible in an HMNS. However, a positive \({\mathcal {A}}\), in agreement with the observations and numerical simulations, could exclude (some parts of) the system to be locally fragmented.