Abstract

To see the evolution of the physical variables at the surface of a star and to establish the range for an unstable system, we investigate the dynamical instability of a spherically symmetric collapsing star in the framework of Rn gravity. For this purpose, we have reformulated the dynamical equations in f(R) theory and perturb them in a specific manner. Under this perturbation scheme, quantities like energy density, pressure and heat flux are obtained and then used in the collapse equation which yield instability range for the ratio of specific heats Γ. Shear-free condition is applied and results are reduced to R2, |$\frac{1}{R}$| and R gravity by choosing particular values of n in the f(R) model Rn.

1 INTRODUCTION

Stars that lived for millions of years and stretched to millions of kilometres in size finally collapse catastrophically within few seconds. The issue of gravitational collapse is playing a major role in both the high-energy astrophysics and modern cosmology. For instance, our recent observational profile about the presence of dark energy and its effects on the expanding cosmos is closely associated with the detected Type Ia supernovae in our Universe. These supernovae are the outcome of white dwarfs and provide information about the universe that how it is expanding with acceleration and the rate by which this acceleration is taking place.

Many theoretical models prompted one to explain the dark energy and cosmic expansion in the framework of Einstein's general theory of relativity. However, some researchers consider another possibility that gravity may deviate from general relativity at large scales. They suppose that the Einstein–Hilbert Lagrangian may be improved by some modifications that may be used to explain the late eras in the history of the universe. The theories based on such proposals are known as modified gravity theories. One of the most popular theories is the f(R) theory whose gravitational action is given by (Capozziello, Carloni & Troisi 2003)
\begin{equation} S_{f(R)}=\frac{1}{2\kappa }\int \mathrm{d}^{4}x\sqrt{-{\boldsymbol g}}f(R), \end{equation}
(1)
where f(R) is a general function of its argument Ricci scalar and g is the metric tensor. This approach of modifying action may also include some corrections from the side of quantum gravity. In the beginning, it was foregone to improve renormalizability of GR (Stelle 1977; Buchbinder, Odintsov & Shapiro 1992), and later on, it is used to understand the inflationary models of the early universe (Starobinsky 1980). It also has some motivation from M theory (Nojiri & Odintsov 2003b), in addition to the required phenomenological characteristics in the recent cosmology.

Various choices of f(R) functions have been formulated by researchers in the literature. However, only some of them satisfy both the theoretical viability criteria and experimental constraints (De Felice & Tsujikawa 2010; Sotiriou & Faraoni 2010). One of such choices is the model f(R) = Rn, where n is not restricted to be an integer (Carloni et al. 2005; Capozziello, Cardone & Troisi 2006; Clifton 2006; Leach, Carloni & Dunsby 2006; Bisabr 2010; Leon & Saridakis 2011). This model, like any f(R) model, is subject to experimental constraints; however, from the mathematical point of view, it is absolutely satisfactory to explore this theory as a toy model. One can also investigate this model in order to find some analytical or qualitative insight on exact solutions. It is mentioned here that the exponent n in Rn gravity is not an entirely free parameter if the theory is meant to constitute a realistic alternative to dark energy.

GR has played an important role in understanding the compact objects, strong gravity fields and radiations with high-energy densities and masses. Gamma-ray bursts emit energy equal to the energy produced by the Sun during its entire life, just in a few seconds. Several models have been proposed in terms of gravitational collapsing mechanism of a single massive star in order to explain such a high burst. The phenomena of gravitational collapse of large massive stars and matter clouds bring up the formation of new stars, galaxies, clusters and other structures in the Universe.

Several investigations on the topic of gravitational collapse have been made by using both the analytical and numerical approaches. However, still a lot of questions regarding complete understanding of the collapsing process remain unanswered from the theoretical and observational view, e.g. the mechanism behind supernova explosions, end state of the gravitational collapse, etc. To explore these aspects more precisely, researchers intend to use modified theories. The motivating idea behind this generalization is that the solutions of this theory include the solutions of GR with many new possible solutions. Just take an example of the Birkhoff theorem stating that the Schwarzschild solution is only the static spherically symmetric vacuum exterior solution and is no longer true in modified theories (Capozziello & Faraoni 2011).

In modified f(R) gravities, the study of gravitational collapse, black hole formation and structure formation is under progress. It is usually carried out in the context of specified f(R) models. In this regard, we have also investigated some issues of gravitational collapse with or without using f(R) models (Sharif & Kausar 2010, 2011a,b,c,d,e,f, 2012; Kausar 2013, 2014). The aim of this paper is to see how the modified Rn gravity affects the dynamics of a collapsing stars. We analyse the instability conditions of gravitational collapse for general function Rn and then assign some particular values to n in the obtained results which represent solutions of certain specific cosmological eras.

The paper is structured as follows. In Section 2, we write modified Einstein field equations and dynamical equations achieved by contracting Bianchi identities. In Section 3, we discuss the f(R) = Rn model in detail and apply perturbations on all the quantities involved in the dynamical equations. Here, we also apply shear-free condition and reduce our results to specific gravities, e.g. R2, 1/R and R gravities. In Section 4, summary is given.

2 GENERAL FORMULATION OF COLLAPSING EQUATION IN f(R) GRAVITY

Let us start with the variation of action (equation 1.1) with respect to the contravariant metric tensor gα β, which gives the following set of field equations:
\begin{eqnarray} &&{f_{,R}R_{\alpha \beta }-\frac{1}{2}fg_{\alpha \beta }-\nabla _{\alpha } \nabla _{\beta }f_{,R}+ g_{\alpha \beta } \Box f_{,R}=\kappa T^m_{\alpha \beta },}\nonumber\\ &&{(\alpha ,\beta =0,1,2,3).} \end{eqnarray}
(2)
Here, comma denotes differentiation with respect to R, i.e. f, R ≡ df(R)/dR, |$\Box =\nabla ^{\alpha }\nabla _{\alpha }$|⁠, with ∇α (covariant derivative) and |$T^m_{\alpha \beta }$| representing the matter part of the field equations. Fourth derivative in the third and fourth terms on the left-hand side introduces the name fourth-order gravity to f(R) theories. Writing above field equation (2.2) in the form of Einstein tensor, we have
\begin{equation} G_{\alpha \beta }= R_{\alpha \beta } - \frac{1}{2}R g_{\alpha \beta } =\frac{\kappa }{f_{,R}} \left[T_{\alpha \beta }^{m}+T_{\alpha \beta }^{({\rm eff})}\right], \end{equation}
(3)
where
\begin{equation} T_{\alpha \beta }^{({\rm eff})}=\frac{1}{\kappa }\left[\frac{f-Rf_{,R}}{2}g_{\alpha \beta }+\nabla _{\alpha } \nabla _{\beta }f_{,R} -g_{\alpha \beta } \Box f_{,R}\right]. \end{equation}
(4)
The above equation |$T_{\alpha \beta }^{({\rm eff})}$| contains all the geometric terms and is known as an effective stress–energy tensor. This way of adopting field equation implies that |$T^{({\rm eff})}_{\alpha \,\beta }$| will not satisfy any energy condition and consequently not a positive definite. Generally, in the comoving coordinates t and r, an interior metric depicting a spherically symmetric collapsing star or cloud depends upon the three unknown functions A(r, t), B(r, t) and C(r, t) written in the form as
\begin{equation} {\rm d}s^2_-=A^2(t,r){\rm d}t^{2}-B^2(t,r){\rm d}r^{2}-C^2(t,r)({\rm d}\theta ^{2}+\sin ^2\theta {\rm d}\phi ^{2}). \end{equation}
(5)
To discuss the properties of collapsing process at the boundary of a star, it is worthwhile to match interior metric describing collapsing star to some known metric in the exterior region. If the Misner and Sharp mass (Misner & Sharp 1964) remains conserved during the collapse of a dust cloud, then it is appropriate to match it with a Schwarzschild space–time in the exterior. If one considers an anisotropic or perfect fluid model, in that case the Misner and Sharp mass usually will not be conserved during collapse. Hence, in such a situation it is useful to match interior space–time with the generalized Vaidya space–time, where mass is not conserved and described in the form of outgoing radiations. Line element for such a metric can be written as (Chao-Guang 1995)
\begin{equation} {\rm d}s^2_+=\left(1-\frac{2M(\nu )}{r}\right){\rm d}\nu ^2+2{\rm d}r{\rm d}\nu -r^2({\rm d}\theta ^2+\sin ^2\theta {\rm d}\phi ^2), \end{equation}
(6)
where M(ν) denotes the total mass inside the considered system. The process of gravitational collapse is extremely dissipative process. Therefore, the heat flux q(t, r) must be present in the interior of a collapsing star. The energy–momentum tensor representing dissipative isotropic fluid can be written as (Herrera, Santos & Le Denmat 1989)
\begin{equation} T^m_{\alpha \beta }=(\rho +p)u_{\alpha }u_{\beta }-pg_{\alpha \beta }+ q_\alpha u_\beta +q_{\beta }u_\alpha . \end{equation}
(7)
Here, ρ denotes the energy density and p denotes the pressure as measured in the rest frame. If a fluid is isotropic in some frame, it leads to a metric which is isotropic in some other frame. The two frames will coincide and fluid will be at rest in comoving coordinates. Then, the four-velocity uα and unit four-vector χα along the radial direction are defined as
\begin{equation} u^{\alpha }=\sqrt{g^{00}}\delta ^{\alpha }_{0},\quad \chi _\alpha =\sqrt{g^{11}}\delta ^\alpha _1,\quad u^{\alpha }u_{\alpha }=1. \end{equation}
(8)
In this case, the norm of the four-velocity of the fluid is always equal to the speed of light which is taken here unity in the gravitational units. Moreover, the heat flux vector satisfies the following equations:
\begin{equation} q^{\alpha }=q\chi ^\alpha ,\quad u^{\alpha }q_{\alpha }=0. \end{equation}
(9)
The shear tensor defines the state of stress at some point of the deformed configuration of the considered system. It is denoted by σαβ and is defined as
\begin{equation} \sigma _{\alpha \beta }=u_{(\alpha ;\beta )}-\bar{a}_{(\alpha }u_{\beta )}-\frac{1}{3}\Theta h_{\alpha \beta }, \end{equation}
(10)
where |$\bar{a}_{\alpha }$| and Θ, respectively, represent four-acceleration and expansion scalar of the fluid. These are given by
\begin{equation} \bar{a}_{\alpha }=u_{\alpha ;\beta }u^{\beta },\quad \Theta =u^{\alpha }_{;\alpha },\ \ h_{\alpha \beta }=g_{\alpha \beta }-u_{\alpha }u_{\beta }. \end{equation}
(11)
Here, semicolon denotes the covariant derivative. From equation (2.10), the non-vanishing components of shear tensor are obtained as follows:
\begin{equation} \sigma _{11}=-\frac{2}{3}{B^{2}}\sigma ,\quad \sigma _{22}=\frac{1}{3}{C^{2}}\sigma ,\quad \sigma _{33}=\sigma _{22}\sin ^2\theta . \end{equation}
(12)
Contracting shear tensor by |$\sigma ^{\alpha \beta }\sigma _{\alpha \beta }=\frac{2}{3}\sigma ^2$|⁠, we have shear scalar σ as
\begin{equation} \sigma =\frac{1}{A}\left(\frac{\dot{B}}{B}-\frac{\dot{C}}{C}\right). \end{equation}
(13)
In the above equation and throughout the paper, dot and prime will denote derivative with respect to t and r, respectively. For general spherically symmetric line element (equation 2.5), the set of field equations is given in Appendix A. In order to explore all the possible aspects of gravitational collapse of a massive star, the basic strategy is to formulate a given set of dynamical equations. These equations should be equipped with some reasonable physical and regularity conditions. For example, positivity of energy density is the basic consistency requirement, while regularity of the initial data also matters by which collapsing process of stars proceeds under their own gravity. Usually initial data for collapse are supplied in the form of initial density and pressure profiles of the matter and the velocities of the concentric collapsing shells inside the stars. Mathematically, such dynamical equations can be formulated by contracting Bianchi identities in the following ways:
\begin{eqnarray} \left(\stackrel{(m)}{T^{\alpha \beta }}+\stackrel{{({\rm eff})}}{T^{\alpha \beta }}\right)_{;\beta }u_{\alpha }=0,\quad \left(\stackrel{(m)}{T^{\alpha \beta }}+\stackrel{{({\rm eff})}}{T^{\alpha \beta }}\right)_{;\beta } \chi _{\alpha }=0. \end{eqnarray}
(14)
For the considered spherically symmetric metric, we obtain the following two independent components of the above equations:
\begin{eqnarray} &&\!\!\!\frac{\dot{\rho }}{A}+\frac{q^{\prime }}{B}+\frac{q}{B}\left(\frac{A^{\prime }}{A} +\frac{2C^{\prime }}{C}\right)+\frac{(\rho +p)}{A} \left(\frac{\dot{B}}{B}+\frac{2\dot{C}}{C}\right)\nonumber\\ &&+\,\frac{A}{\kappa }F_1=0, \end{eqnarray}
(15)
\begin{eqnarray} \frac{p^{\prime }}{B}+\frac{\dot{q}}{A}+\frac{2q}{A}\left(\frac{\dot{B}}{B} +\frac{\dot{C}}{C}\right)+\frac{A^{\prime }}{AB}(\rho +p)-\frac{B}{\kappa }F_2=0, \end{eqnarray}
(16)
where F1 and F2 are lengthy expressions containing terms of effective energy–momentum tensor, and can be found in Kausar (2014). It is mentioned here that above equations are the basic focus of this manuscript which would help to discuss the proposed topic.

3 Rn GRAVITY AND PERTURBATION SCHEME

Here, we have considered the following model in modified f(R) theory of gravity (Nojiri & Odintsov 2003a),
\begin{equation} f(R)=R^n. \end{equation}
(17)
This model has been proposed to explain the recent cosmic acceleration due to the presence of dark energy. Especially, it is very helpful to represent the exponential expansion, i.e. power-law inflation and ordinary inflation with a minimally coupled scalar field. In the literature, this model is mostly pursued for n = −1 and 3/2 (Capozziello et al. 2003; Carloni et al. 2005), where it is conformally equivalent to Liouville field theory (Capozziello et al. 2003). For the values n ≠ 0, 1/2, 1, this theory gives power-law inflation in the form given below (Muller, Schmidt & Starobinsky 1990; Capozziello, Occhionero & Amendola 1993):
\begin{equation} \hat{a}(t)\rightarrow t^{\hat{\alpha }},\ \ {\rm with}\ \ \hat{\alpha }=\frac{-2n^2+3n-1}{n-2}, \end{equation}
(18)
where |$\hat{a}(t)$| is the cosmic scale factor.
Practically, a modified f(R) theory is considered to be meaningful if it admits the existence of de Sitter solution. In this respect, any arbitrary function f(R) should fulfil this requirement by satisfying the following condition (Barrow & Ottewill 1983):
\begin{equation} R_{\rm c}f_{,R_{\rm c}}=2f_{\rm c}, \end{equation}
(19)
in the background of constant scalar curvature (Hc, Rc) (Faraoni & Nadeau 2005). Here, subscript ‘c’ represents that above quantities are evaluated at the constant curvature |$R_{\rm c}=12H_{\rm c}^2$| with Hc being the Hubble constant. For the case of Rn gravity, the above condition becomes
$$ nR^n_{\rm c} = 2R^n_{\rm c}, $$
which clearly satisfied only for n = 2 or Rc = 0. Hence, this modified model does not accept a de Sitter-type solution for n ≠ 2 until a term with n = 0, corresponding to cosmological constant, would be added to the modified Lagrangian (Barrow & Ottewill 1983; Faraoni 2004). However, for any positive value of n, the Minkowski space would trivially be a solution without any such addition.
The other important condition for any f(R) model to be viable is to satisfy the following stability condition in the background of de Sitter solution (Muller, Schmidt & Starobinsky 1988):
\begin{equation} \frac{f_{,R_{\rm c}}^2-2f_{\rm c}f_{,R_{\rm c\,}R_{\rm c}}}{f_{,R_{\rm c}}f_{,R_{\rm c\,}R_{\rm c}}}\ge 0. \end{equation}
(20)
For our chosen model, this stability condition becomes
\begin{equation} R_{\rm c}\frac{2-n}{n(n-1)}\ge 0. \end{equation}
(21)
It is noticed that for 1 < n ≤ 2 and n < 0, this condition may be fulfilled, while for n = 2, it is also identically satisfied without inflicting any constraint on the Hubble constant. Moreover, the Minkowski space–time is stable for all n > 0.
The strategy of perturbing a system helps to study the dynamical process of gravitational collapse. In this scheme, all the quantities whether material or metric are supposed to have only radial dependence initially. Afterwards, they start depending on time with the inclusion of 0 < ϵ ≪ 1 in their perturbation. We can introduce this scheme as
\begin{eqnarray} A(t,r)&=&A_0(r)+\epsilon T(t)\bar{a}(r), \end{eqnarray}
(22)
\begin{eqnarray} B(t,r)&=&B_0(r)+\epsilon T(t)\bar{b}(r), \end{eqnarray}
(23)
\begin{eqnarray} C(t,r)&=&C_0(r)+\epsilon T(t)\bar{c}(r), \end{eqnarray}
(24)
\begin{eqnarray} \rho (t,r)&=&\rho _0(r)+\epsilon \bar{\rho }(t,r), \end{eqnarray}
(25)
\begin{eqnarray} p(t,r)&=&p_0(r)+\epsilon \bar{p}(t,r), \end{eqnarray}
(26)
\begin{eqnarray} q(t,r)&=&\epsilon \bar{q}(t,r), \end{eqnarray}
(27)
\begin{eqnarray} \sigma (t,r)&=&\epsilon \bar{\sigma }(t,r), \end{eqnarray}
(28)
\begin{eqnarray} R(t,r)&=&R_0(r)+\epsilon T(t)\bar{e}(r). \end{eqnarray}
(29)
Substituting equation (3.29) in the chosen f(R) model (equation 3.17), we have
\begin{eqnarray} f&=&R_0^n+\epsilon T\bar{e}nR_0^{n-1}, \end{eqnarray}
(30)
\begin{eqnarray} f_{,R}&=&nR_0^{n-1}+\epsilon T \bar{e}n(n-1)R_0^{n-2}, \end{eqnarray}
(31)
\begin{eqnarray} &&{f_{,RR}=n(n-1)R_0^{n-2}+\epsilon T\bar{e}n(n-1)(n-2)R_0^{n-3},}\nonumber \\ &&{f_{,RRR}=n(n-1)(n-2)R_0^{n-3}} \end{eqnarray}
(32)
\begin{eqnarray} +\epsilon T\bar{e}n(n-1)(n-2)(n-3)R_0^{n-4}, \end{eqnarray}
(33)
where |$\frac{{\rm d}f(R)}{{\rm d}R}=f_{,R}$|⁠, |$\frac{{\rm d}^2f(R)}{{\rm d}R^2}=f_{,RR}$| and |$\frac{{\rm d}^3f(R)}{{\rm d}R^3}=f_{,RRR}$|⁠. Assuming C0 = r and perturbing Ricci scalar curvature, we have the following differential equation:
\begin{equation} \ddot{T}-\omega (r) T=0, \end{equation}
(34)
where ω(r) is given by
\begin{eqnarray} \omega (r)&=&\left[-\frac{1}{A_0B_0}\left(\bar{a}^{\prime \prime }-\frac{\bar{a}A_0^{\prime \prime }}{A_0}-\frac{\bar{c}B_0^{\prime }}{r} -\frac{3\bar{b}B_0^{\prime }}{B_0}\frac{2\bar{b}A_0^{\prime \prime }}{B_0}\right)\right.\nonumber \\ &&- -\frac{2}{rB_0^3}\left(\bar{b}^{\prime }+\bar{c}^{\prime }B_0^{\prime }+\frac{\bar{e}}{2}+\frac{\bar{c}}{r^3}\right.\nonumber \\ &&-\frac{1}{A_0B_0^3}\left(\bar{a}^{\prime }B_0^{\prime }+\bar{b}A_0^{\prime } -\frac{\bar{a}A_0^{\prime }B_0^{\prime }}{A_0} -\frac{3\bar{b}A_0^{\prime }B_0^{\prime }}{B_0}\right)\nonumber \\ &&+\frac{1}{B_0^2r^2}\left(\bar{c}^{\prime }-\frac{\bar{b}}{B_0} -\frac{\bar{c}}{r}\right)+\left.\frac{2\bar{c}^{\prime \prime }}{rB_0^2}\right]\frac{A_0^2}{\left(\frac{\bar{b}}{B_0} -\frac{2\bar{c}}{r}\right)}. \end{eqnarray}
(35)
By keeping ω(r) to be positive, the solution of the above differential equation is given by
\begin{equation} T(t)=-{\rm e}^{\sqrt{\omega }t}. \end{equation}
(36)
When we use the static part of equations (3.22)–(3.33) in the dynamical equation (2.15), it is satisfied identically, whereas other dynamical equation (2.16) becomes
\begin{eqnarray} \frac{p_0^{\prime }}{B_0}+(\rho _0+p_0)\frac{A_0^{\prime }}{A_0B_0}+\frac{B_0}{\kappa } F_{2{\rm s}}=0, \end{eqnarray}
(37)
where F2s is the static part of equation F2 and is given by
\begin{eqnarray} F_{2{\rm s}}&=&(1-n)\frac{R_0^n}{2B_0^2}+\frac{n(n-1)}{B_0^2}R_0^{n-2}R_0^{\prime } +\frac{n}{B_0^4}\frac{A_0^{\prime }}{A_0}(n-1)\nonumber \\ &&\times \left[(n-2)R_0^{n-3}R_0^{\prime 2} +(n-2)R_0^{n-3}R_0^{\prime \prime }\right.\nonumber \\ &&\left.-\,R_0^{n-2}R_0^{\prime }\left(\frac{A_0^{\prime }}{A_0} -\frac{B_0^{\prime }}{B_0}\right)\right]\nonumber\\ &&+\, \frac{2B_0^{\prime }}{B_0^3}\left[(1-n)\frac{R_0^n}{2A_0^2} +\frac{n(n-1)}{B_0^2}R_0^{\prime } R_0^{n-2}\left(\frac{A_0^{\prime }}{A_0}+\frac{B_0^{\prime }}{B_0}\right)\right]. \nonumber\\ \end{eqnarray}
(38)
Applying full perturbations on equation (2.15), we obtain
\begin{eqnarray} &&\!\!\!\!\frac{\dot{\bar{\rho }}}{A_0}+\frac{{\bar{q}}^{\prime }}{B_0}+{\bar{q}} \left(\frac{A_0^{\prime }}{A_0}+\frac{2}{r}\right)+\frac{(\rho _0+p_0)}{A_0}\left(\frac{\bar{b}}{B_0} +\frac{2\bar{c}}{r}\right)\dot{T}\nonumber\\ &&+\,\frac{A_0}{\kappa }F_{1{\rm p}}\dot{T}=0, \end{eqnarray}
(39)
where F1p is given in Appendix A. Applying equations (3.22)–(3.33) and using equation (3.37) in the second dynamical equation, we obtain
\begin{eqnarray} &&\!\!\!\!\frac{{\bar{p}}^{\prime }}{B_0}+\frac{\dot{\bar{q}}}{A_0}+({{\bar{\rho }}}+{\bar{p}}) \frac{A_0^{\prime }}{A_0B_0}+T\left[\frac{\bar{b}}{B_0}p_0^{\prime }+\frac{1}{B_0}(\rho _0+p_0) \left(\frac{\bar{a}}{A_0}\right)^{\prime }\right.\nonumber \\ &&+\left.\frac{\bar{b}}{B_0^2}\frac{A_0^{\prime }}{A_0}(\rho _0+p_0)\right] +\frac{B_0}{\kappa }F_{2{\rm p}}=0, \end{eqnarray}
(40)
where F2p carries the perturbed part of the expression F2 and is given in Appendix A. Equation (3.40) can also be written as
\begin{eqnarray} {\bar{p}}^{\prime }+{\bar{p}}\frac{A_0^{\prime }}{A_0}&=&- \frac{{{\bar{\rho }}}A_0^{\prime }}{A_0}-\frac{\dot{\bar{q}}B_0}{A_0}+T\left[\bar{b}p_0^{\prime }+(\rho _0+p_0) \left(\frac{\bar{a}}{A_0}\right)^{\prime }\right.\nonumber \\ &&+\left.\frac{\bar{b}}{B_0}\frac{A_0^{\prime }}{A_0}(\rho _0+p_0)\right] +\frac{B_0^2}{\kappa }F_{2{\rm p}}. \end{eqnarray}
(41)
Second law of thermodynamics helps to build a connection between |$\bar{\rho }$| and |$\bar{p}$| in the form of specific heats. Harrison and Wheeler introduce an equation of state in this regard which is given by (Harrison et al. 1965; Chan et al. 1989)
\begin{equation} \bar{p}=\Gamma \frac{p_0}{\rho _0+p_0}\bar{\rho }. \end{equation}
(42)
Here, Γ is known as the adiabatic index or ratio of the heat capacity at constant pressure to the heat capacity at constant volume (Herrera et al. 1989).
Taking it constant and putting value of |$\bar{p}$| from the above equation in equation (3.41) and rearranging, it becomes
\begin{eqnarray} \Gamma &=&-\left[\left(\frac{p_0{\bar{\rho }}}{\rho _0+p_0}\right)^{\prime } +\left(\frac{p_0{{\bar{\rho }}}}{\rho _0+p_0}\right)\frac{A_0^{\prime }}{A_0}\right]^{-1}\left[ \frac{{{\bar{\rho }}}A_0^{\prime }}{A_0}+\frac{\dot{\bar{q}}B_0}{A_0}\right.\nonumber \\ &&\left.+\,T\left\lbrace \bar{b}p_0^{\prime }+(\rho _0+p_0) \left(\frac{\bar{a}}{A_0}\right)^{\prime }+\frac{\bar{b}}{B_0}\frac{A_0^{\prime }}{A_0}(\rho _0+p_0)\right\rbrace\right. \nonumber\\ &&\left.+\,\frac{B_0^2}{\kappa }F_{2{\rm p}}\right]. \end{eqnarray}
(43)
The above equation is the focus of our concern, so first we have to find quantities involved in it. Applying perturbations on equation (A2) and solving for |$\bar{q}$|⁠, we have
\begin{eqnarray} \bar{q}&=&\frac{nR_0^{n-1}\dot{T}}{\kappa A_0B_0}\left[\frac{-2}{r}\left(\bar{c}^{\prime }-\bar{c}\frac{A_0^{\prime }}{A_0}-\frac{\bar{b}}{B_0}\right)\right.\nonumber \\ &&\left.+\,\frac{(n-1)}{R_0^{n-1}}\left\lbrace (n-2)R_0^{n-3}\bar{e}R_0^{\prime }+R_0^{n-2}\bar{e}^{\prime }-\frac{A_0^{\prime }}{A_0}(n-2)R_0^{n-3}\right.\right.\nonumber\\ &&\left.\left.-\,\frac{\bar{b}}{B_0}R_0^{n-2}R_0^{\prime }\right\rbrace \right]. \end{eqnarray}
(44)
Differentiating the above equation with respect to ‘r’ and substituting it in equation (3.39), and then integrating the resulting equation with respect to ‘t’, we obtain
\begin{eqnarray*} {\bar{\rho }}&=&-\frac{nA_0}{\kappa B_0}T\left[\frac{R_0^{n-1}}{A_0B_0}\left\lbrace \frac{-2}{r}\left(\bar{c^{\prime }}-\bar{c}\frac{A_0^{\prime }}{A_0}-\frac{\bar{b}}{B_0}\right) +\frac{n-1}{R_0^{n-1}}\right.\right.\nonumber \\ &&\times \left\lbrace (n-2)R_0^{n-3}\bar{e}R_0^{\prime }+R_0^{n-2}\bar{e}^{\prime }-\frac{A_0^{\prime }}{A_0}(n-2)R_0^{n-3} \right.\nonumber \\ &&-\left.\left.\left.\frac{\bar{b}}{B_0}R_0^{n-2}R_0^{\prime }\right\rbrace \right\rbrace \right]_{,r}-\frac{nR_0^{n-1}T}{\kappa B_0^2}\left[\frac{-2}{r}\left(\bar{c^{\prime }}-\bar{c}\frac{A_0^{\prime }}{A_0} -\frac{\bar{b}}{B_0}\right)\right.\nonumber\\ &&+\,\frac{n-1}{R_0^{n-1}} \left\lbrace (n-2)\bar{e}R_0^{n-3}R_0^{\prime }+R_0^{n-2}\bar{e}^{\prime }\right.\nonumber \\ &&-\left.\left. \frac{A_0^{\prime }}{A_0}(n-2)R_0^{n-3}-\frac{\bar{b}}{B_0}R_0^{n-2}R_0^{\prime }\right\rbrace \right] \nonumber \\ &&-\,(\rho _0+p_0)\left(\frac{\bar{b}}{B_0}+\frac{2\bar{c}}{r}\right)T-\frac{A_0^2}{\kappa }F_{1{\rm p}}T. \end{eqnarray*}
As equation (3.43) is the required evolution equation, so we further analyse it under some physical conditions in the subsections below. It is worthwhile to notice that how the quantities in the above equations depend upon the chosen model and hence influence the evolution of the physical variables at the surface of a star.

3.1 Shear-free condition in Newtonian regime

Shear tensor plays a significant role in the evolution of self-gravitating bodies and in the violation of cosmic censorship. Its vanishing in the interior of collapsing stars has many interesting consequences which are brought out by many authors (Glass 1979; Collins & Wainwright 1983; Chan 1998; Joshi, Dadhich & Maartens 2002; Herrera & Santos 2003). Especially, it influences the space–time singularity during collapse and hence causes the appearance of a naked singularity which is visible to a far-away observers in the universe (Joshi et al. 2002). Moreover, non-dissipative shear-free flow leads to a well-known homologous evolution (Herrera, Santos & Wang 2008) which have a close relevance in astrophysics (Schwarzschild 1958; Kippenhahn & Weigert 1990; Hansen & Kawaler 1994). Applying shear-free condition, i.e. σ = 0 in equation (2.13), we obtain the following relation after perturbing
\begin{equation} \frac{\bar{b}}{B_0}=\frac{\bar{c}}{r}. \end{equation}
(45)
Considering Newtonian regime, we approximate metric components A0 = 1, B0 = 1, and suppose ρ0 ≫ p0. Hence, putting |$A_0^{\prime }=0,\ B_0^{\prime }=0$| and value of T in equation (3.43), we obtain
\begin{eqnarray} \Gamma <-\frac{1}{{\bar{\rho }^{\prime }}(N)} \left[ \dot{\bar{q}}(N)-{\rm e}^{\sqrt{\omega }t}\left(\frac{\bar{c}p_0^{\prime }}{r}+\rho _0\bar{a}^{\prime }\right)-\frac{F_{2{\rm p}}}{\kappa }\right]. \end{eqnarray}
(46)
Here, we have imposed a physical condition on the pressure for the collapsing system, i.e. |$p_0^{\prime }<0$|⁠. This equation shows that the adiabatic index depends upon the pressure, heat flux, density and higher order curvature terms in the Newtonian regime. Furthermore, if this equation would be satisfied, a self-gravitating star/cloud may keep on collapsing, otherwise it would be stable and there is no more collapse. The quantities |$\bar{q}$|⁠, |$\bar{\rho }$| and F2p(N) involved in the above equation corresponding to the Newtonian regime are provided in Appendix A. Now, we recall and compare our results by the numerical results of Chandrasekhar, who obtained the following three possibilities for the ratio of specific heats.
  • If Γ > 4/3, the pressure increases strongly with the contraction than the weight of the star. The resulting force is directed outwards and the star will move back towards equilibrium. Then, the star becomes dynamically stable.

  • If Γ < 4/3, the weight increases stronger than the pressure, and the star would collapse after the initial compression. This is dynamical instability.

  • If Γ = 4/3, the compression again leads to hydrostatic equilibrium.

The critical value 4/3 depends strongly on spherical symmetry and Newtonian gravitation. The number 4 in the numerator comes from the fact that the weight of the envelope in Newtonian mechanics varies as ∼r−2, and has to be distributed over the surface of the sphere, giving another r−2. The number 3 in the denominator comes from r3 in the formula for the volume of the sphere. Relativistic effects of general relativity may change the critical value of Γ and make the models less stable. As we have used an analytical approach to investigate the problem of dynamical instability of supermassive spherically symmetric collapsing stars in f(R) gravity, our result appears in the form of matter variables like pressure and density and curvature terms due to modifications in gravity. Hence, if inequality (3.46) is satisfied, then weight of the star increases very fast than the pressure and the star collapses resulting in a dynamical instability, whereas in the opposite case, the pressure in a star is strong enough than the weight of the outer layers which makes it less stable to exist.

In the following subsections, we reduce our obtained results for some specific values of n. These chosen values represent specific eras in the history of the universe.

3.2 R2 gravity

Putting power n = 2 in the proposed model, we have f(R) = R2 gravity. Presence of the R2 term in the action leads to the accelerated expansion of the universe. This was the initiative model proposed by Starobinsky in 1980 to explain inflation (Starobinsky 1982). In that work, it was shown that the f(R) model with R2 term is well consistent with the temperature anisotropies found in the cosmic microwave background data and hence can be used as a viable alternative to the scalar field models of inflation. For this model, the expression for adiabatic index will remain the same, whereas the involved quantities in equations (A8)–(A11) become those in equations (A12)–(A15) in Appendix.

3.3 |$\frac{1}{R}$| gravity

This regime is obtained by substituting n = −1 in equation (3.17). Inclusion of this term in the Einstein–Hilbert action proposed to explain the late-time acceleration of the universe due to dark energy. In this scenario, equations (A8)–(A11) will reduce to equations (A16)–(A19). However, Γ will remain the same.

3.4 Einstein gravity

In this case, we have recovered the Einstein gravity by substituting n = 1 in equation (3.17). Hence, we may consider that the following result of the instability equation represents GR result due to vanishing of all the effective curvature terms arising due to modified gravity:
\begin{eqnarray} \Gamma &<&\frac{\left[ \frac{2\omega }{r\kappa }\left(\frac{\bar{c}}{r}-\bar{c}^{\prime }\right) +\frac{\bar{c}p_0^{\prime }}{r}+\rho _0\bar{a}^{\prime }\right]}{\left[\frac{2}{r\kappa }\left(\frac{\bar{c}}{r}-\bar{c^{\prime }}\right)+\frac{1}{\kappa } \left\lbrace \frac{2}{r}\left(\frac{\bar{c}}{r}-\bar{c^{\prime }}\right)\right\rbrace _{,r} +\rho _0\frac{3\bar{c}}{r}\right]_{,r}}. \end{eqnarray}
(47)
If this inequality holds, the system becomes unstable and may start collapsing, provided that |$\frac{\bar{c}}{r}\ge \bar{c^{\prime }}$|⁠.

4 SUMMARY

In conclusion, the theories of f(R) gravity have improved our understanding of the peculiarities of GR in a broader spectrum. Their simple generalizations help us to construct viable models as an alternative to dark energy which are used as toy models to explain the accelerating expanding universe. Although presently there is no definite prediction that sets them apart from other dark energy models, still they have taught us about important aspects of relativistic theories of gravity.

In this paper, we have considered the Rn gravity and explore the issue of gravitational collapse. For this purpose, we have reformulated the dynamical equations governing the dynamical process of a collapsing star in this modified theory. The obtained equations are carrying lengthy expressions as compared to GR due to higher order curvature terms. We have applied a perturbation scheme on Rn gravity as well as on all the quantities used in the dynamical equations. After lengthy calculations, we found dynamical equation (3.43) in the form of adiabatic index Γ. This equation is the required collapse equation yielding the range for an unstable system and shows that how system depends on the heat flux, density, pressure and effective curvature terms during collapsing process. It is mentioned here that, in order to keep the system unstable and carry out the collapse, all the terms involved in the adiabatic index should be positive. So, there must be many constraints on the physical variables at the surface of a star.

We have also applied Newtonian limit approximation to make our results comparable with the Chandrasekhar results in the same regime. The results are further reduced to R2, 1/R and R gravities by giving specific values to the power n of the considered model. These meaningful values correspond to the certain periods in the history of the universe and used to explain the inflation and late-time acceleration in relativistic astrophysics. Hence, this paper provides dynamical instability of general spherically symmetric gravitational collapse in many backgrounds at the same time, even recovering the GR results too. It is worthwhile to explore such issues of astrophysics using modified gravity models in order to obtain more general and flexible results as compared to GR.

REFERENCES

Barrow
J. D.
Ottewill
A. C.
J. Phys. A
1983
, vol. 
16
 pg. 
2757
 
Bisabr
Y.
Gravit. Cosmol.
2010
, vol. 
16
 pg. 
239
 
Buchbinder
I. L.
Odintsov
S. D.
Shapiro
I. L.
Effective Action in Quantum Gravity
1992
Bristol
IoP Publishing
Capozziello
S.
Faraoni
V.
Beyond Einstein Gravity
2011
New York
Springer
Capozziello
S.
Occhionero Amendola
L.
Int. J. Mod. Phys. D
1993
, vol. 
1
 pg. 
615
 
Capozziello
S.
Carloni
S.
Troisi
A.
Res. Dev. Astron. Astrophys.
2003
, vol. 
1
 pg. 
625
 
Capozziello
S.
Cardone
V. F.
Troisi
A.
Phys. Rev. D
2006
, vol. 
73
 pg. 
104019
 
Carloni
S.
Dunsby
P. K. S.
Capozziello
S.
Troisi
A.
Class. Quantum Grav.
2005
, vol. 
22
 pg. 
4839
 
Chan
R.
MNRAS
1998
, vol. 
299
 pg. 
811
 
Chan
R.
Santos
N. O.
Kichenassamy
S.
Le Denmat
G.
MNRAS
1989
, vol. 
239
 pg. 
91
 
Chao-Guang
H.
Acta Phys. Sin.
1995
, vol. 
4
 pg. 
617
 
Clifton
T.
Class. Quantum Grav.
2006
, vol. 
23
 pg. 
7445
 
Collins
C. B.
Wainwright
J.
Phys. Rev. D
1983
, vol. 
27
 pg. 
1209
 
De Felice
A.
Tsujikawa
S.
Living Rev. Relativ.
2010
, vol. 
13
 pg. 
3
 
Faraoni
V.
Phys. Rev. D
2004
, vol. 
70
 pg. 
044037
 
Faraoni
V.
Nadeau
S.
Phys. Rev. D
2005
, vol. 
72
 pg. 
124005
 
Glass
E. N.
J. Math. Phys.
1979
, vol. 
20
 pg. 
1508
 
Hansen
C.
Kawaler
S.
Stellar Interiors; Physical Principles, Structure and Evolution
1994
Berlin
Springer-Verlag
Harrison
B. K.
Thorne
K. S.
Wakano
M.
Wheeler
J. A.
Gravitation Theory and Gravitational Collapse
1965
Chicago
University of Chicago Press
Herrera
L.
Santos
N. O.
MNRAS
2003
, vol. 
343
 pg. 
1207
 
Herrera
L.
Santos
N. O.
Le Denmat
G.
MNRAS
1989
, vol. 
237
 pg. 
257
 
Herrera
L.
Santos
N. O.
Wang
A.
Phys. Rev. D
2008
, vol. 
78
 pg. 
084026
 
Joshi
P. S.
Dadhich
N.
Maartens
R.
Phys. Rev. D
2002
, vol. 
65
 pg. 
101501
 
Kausar
H. R.
J. Cosmol. Astropart. Phys.
2013
, vol. 
01
 pg. 
07
 
Kausar
H. R.
MNRAS
2014
, vol. 
439
 pg. 
1536
 
Kippenhahn
R.
Weigert
A.
Stellar Structure and Evolution
1990
Berlin
Springer-Verlag
Leach
J. A.
Carloni
S.
Dunsby
P. K. S.
Class. Quantum Grav.
2006
, vol. 
23
 pg. 
4915
 
Leon
G.
Saridakis
E. N.
Class. Quantum Grav.
2011
, vol. 
28
 pg. 
065008
 
Misner
C. W.
Sharp
D.
Phys. Rev.
1964
, vol. 
136
 pg. 
571
 
Muller
V.
Schmidt
H.-J.
Starobinsky
A. A.
Phys. Lett. B
1988
, vol. 
202
 pg. 
198
 
Muller
V.
Schmidt
H.-J.
Starobinsky
A. A.
Class. Quantum Grav.
1990
, vol. 
7
 pg. 
1163
 
Nojiri
S.
Odintsov
S. D.
Phys. Rev. D
2003a
, vol. 
68
 pg. 
123512
 
Nojiri
S.
Odintsov
S. D.
Phys. Lett. B
2003b
, vol. 
576
 pg. 
11
 
Schwarzschild
M.
Structure and Evolution of the Stars
1958
New York
Dover
Sharif
M.
Kausar
H. R.
Mod. Phys. Lett. A
2010
, vol. 
25
 pg. 
3299
 
Sharif
M.
Kausar
H. R.
J. Cosmol. Astropart. Phys.
2011a
, vol. 
07
 pg. 
022
 
Sharif
M.
Kausar
H. R.
Int. J. Mod. Phys. D
2011b
, vol. 
20
 pg. 
2239
 
Sharif
M.
Kausar
H. R.
J. Phys. Soc. Japan
2011c
, vol. 
80
 pg. 
044004
 
Sharif
M.
Kausar
H. R.
Astrophys. Space Sci.
2011d
, vol. 
331
 pg. 
281
 
Sharif
M.
Kausar
H. R.
Astrophys. Space Sci.
2011e
, vol. 
332
 pg. 
463
 
Sharif
M.
Kausar
H. R.
Phys. Lett. B
2011f
, vol. 
697
 pg. 
1
 
Sharif
M.
Kausar
H. R.
Astrophys. Space Sci.
2012
, vol. 
337
 pg. 
805
 
Sotiriou
T. P.
Faraoni
V.
Rev. Mod. Phys.
2010
, vol. 
82
 pg. 
451
 
Starobinsky
A. A.
Phys. Lett. B
1980
, vol. 
91
 pg. 
99
 
Starobinsky
A. A.
Nonsingular Model of the Universe with the Quantum-Gravitational de Sitter Stage and Its Observational Consequences
1982
Moscow
INR Press
pg. 
58
 
Stelle
K. S.
Phys. Rev. D
1977
, vol. 
16
 pg. 
953
 

APPENDIX A

Field equations
\begin{eqnarray} G_{00}&=&\frac{\kappa A^{2}}{f_{,R}}\left[{\rho }+\frac{1}{\kappa }\left\lbrace \frac{f-Rf_{,R}}{2}+\frac{1}{B^2}\left(f_{,RRR}R^{\prime 2}+f_{,RR}R^{\prime \prime }\right)+\frac{f_{,RR}\dot{R}}{A^2}\left(\frac{2\dot{C}}{C}-\frac{\dot{B}}{B}\right) +\frac{f_{,RR}R^{\prime }}{B^2}\left(\frac{2C^{\prime }}{C}-\frac{B^{\prime }}{B}\right)\right\rbrace \right], \end{eqnarray}
(A1)
\begin{eqnarray} G_{01}&=&\frac{-\kappa }{f_{,R}}\left[qAB-\frac{1}{\kappa }\left(f_{,RRR}\dot{R}R^{\prime }+f_{,RR}\dot{R}^{\prime } -\frac{A^{\prime }}{A}f_{,RR}\dot{R}-\frac{\dot{B}}{B}f_{,RR}R^{\prime }\right)\right], \end{eqnarray}
(A2)
\begin{eqnarray} G_{11}&=&\frac{\kappa B^{2}}{f_{,R}}\left[p-\frac{1}{\kappa } \left\lbrace \frac{f-Rf_{,R}}{2}-\frac{1}{A^2}\left(f_{,RRR}\dot{R}^2+f_{,RR}\ddot{R}\right)\frac{f_{,RR}R^{\prime }}{B^2} \left(\frac{A^{\prime }}{A}+\frac{2C^{\prime }}{C}\right)\right\rbrace \right], \end{eqnarray}
(A3)
\begin{eqnarray} G_{22}&=&\frac{\kappa C^{2}}{f_{,R}}\left[p-\frac{1}{\kappa }\left\lbrace \frac{f-Rf_{,R}}{2}-\frac{1}{A^2}\left(f_{,RRR}\dot{R}^2+f_{,RR}\ddot{R}\right)+\frac{1}{B^2}\left(f_{,RRR}R^{\prime 2}+f_{,RR}R^{\prime \prime }\right)+ \frac{f_{,RR}\dot{R}}{A^2}\left(\frac{\dot{A}}{A}-\frac{\dot{B}}{B}+\frac{\dot{C}}{C}\right)\right.\right.\nonumber \\ &&+\left.\left. \frac{f_{,RR}R^{\prime }}{B^2}\left(\frac{A^{\prime }}{A}-\frac{B^{\prime }}{B}+\frac{C^{\prime }}{C}\right)\right\rbrace \right], \end{eqnarray}
(A4)
where
\begin{eqnarray} G_{00}&=&\left(\frac{2\dot{B}}{B}+\frac{\dot{C}}{C}\right) \frac{\dot{C}}{C}-\left(\frac{A}{B}\right)^2 \left[\frac{2C^{\prime \prime }}{C}+\left(\frac{C^{\prime }}{C}\right)^2-\frac{2B^{\prime }C^{\prime }}{BC} -\left(\frac{B}{C}\right)^2\right],\nonumber \\ G_{01}&=&-2\left(\frac{\dot{C^{\prime }}}{C}-\frac{\dot{C}A^{\prime }}{CA}-\frac{\dot{B}C^{\prime }}{BC}\right),\nonumber\\ G_{11}&=&-\left(\frac{B}{A}\right)^2\left[\frac{2\ddot{C}}{C}-\left(\frac{2\dot{A}}{A} -\frac{\dot{C}}{C}\right) \frac{\dot{C}}{C}\right] +\left(\frac{2A^{\prime }}{A}+\frac{C^{\prime }}{C}\right)\frac{C^{\prime }}{C}-\left(\frac{B}{C}\right)^2, \nonumber\\ G_{22}&=&-\left(\frac{C}{A}\right)^2\left[\frac{\ddot{B}}{B}-\frac{\ddot{C}}{C}-\frac{\dot{A}}{A} \left(\frac{\dot{B}}{B}+\frac{\dot{C}}{C}\right)+\frac{\dot{B}\dot{C}}{BC}\right] +\left(\frac{C}{B}\right)^2 \left[\frac{A^{\prime \prime }}{A}+\frac{C^{\prime \prime }}{C}-\frac{A^{\prime }B^{\prime }}{AB}+\left(\frac{A^{\prime }}{A}-\frac{B^{\prime }}{B}\right)\frac{C^{\prime }}{C}\right]. \end{eqnarray}
(A5)
\begin{eqnarray} \nonumber F_{1{\rm p}}&=&\frac{n(n-1)}{A_0^2B_0^2}\left[(n-2)R_0^{n-3}R_0^{\prime }+R_0^{n-2}\bar{e}^{\prime }\right]-\frac{1}{2A_0^2}\left\lbrace n(n-1)R_0^{n-1}-\frac{\bar{a}}{A_0}(1-n)R_0^n\right\rbrace -\frac{2}{A_0^2B_0^2}\left(\frac{\bar{a}}{A_0}+\frac{\bar{b}}{B_0}\right) \left[n(n-1)(n-2)R_0^{n-3}\right.\nonumber\\ &&\times \left. \left(R_0^{\prime 2}+R_0^{\prime \prime }\right)\right]+\frac{n(n-1)(n-2)}{A_0^2B_0^2} \left[\bar{e}(n-3)R_0^{n-4}R_0^{\prime }+2\bar{e}^{\prime }R_0^{n-3}+\bar{e}^{\prime \prime }R_0^{n-3}+\bar{e}(n-3)R_0^{n-4}R_0^{\prime \prime }\right]-\frac{n(n-1)}{B_0^2} \left[R_0^{\prime }R_0^{n-2}\left\lbrace \left(\frac{\bar{b}}{B_0}\right)^{\prime }\right.\right.\nonumber\\ &&-\left.\left.\left(\frac{2}{r}\right)^{\prime }\right\rbrace +\left(\frac{B_0^{\prime }}{B_0}+\frac{1}{r}\right)\left\lbrace \bar{e}(n-2)R_0^{n-3}R_0^{\prime }+ \left(\bar{e}^{\prime }-\frac{\bar{b}}{B_0}R_0^{\prime }\right)R_0^{n-2}\right\rbrace \right] +\frac{\bar{a}}{A_0^2}\left[\frac{(1-n)R_0^n}{2A_0^2}+\frac{n(n-1)(n-2)R_0^{n-3}}{A_0^2B_0^2}\left(R_0^{\prime 2}+R_0^{\prime \prime }\right)\right.\nonumber\\ &&-\left. \frac{1}{B_0^2}n(n-1)R_0^{n-2}R_0^{\prime }\left(\frac{B_0^{\prime }}{B_0}-\frac{2}{r}\right)\right] +\frac{\bar{b}}{A_0^2B_0^3}n(n-1)\left[(n-2)R_0^{n-3}(R_0^{\prime 2}+R_0^{\prime \prime })-R_0^{n-2}R_0^{\prime } \left(\frac{A_0^{\prime }}{A_0}+\frac{B_0^{\prime }}{B_0}\right)\right]-\frac{2R_0^{n-2}R_0^{\prime }}{rA_0^2B_0^2} \left(\frac{A_0^{\prime }}{A_0}+\frac{1}{r}\right)\nonumber\\ &&+\,\frac{n(n-1)}{A_0^2B_0^2}\left(\frac{2A_0^{\prime }}{A_0} +\frac{B_0^{\prime }}{B_0}+\frac{1}{r}\right)\left[(n-2)\bar{e}R_0^{n-3}R_0^{\prime }+\bar{e}^{\prime }R_0^{n-2}-\bar{e}\frac{A_0^{\prime }}{A_0}R_0^{n-2}-\frac{\bar{b}}{B_0}R_0^{\prime }R_0^{n-2}\right]. \end{eqnarray}
(A6)
\begin{eqnarray} \nonumber F_{2{\rm p}}&=&\ddot{T}n(n-1)\left[\frac{(n-2)R_0^{n-3}R_0^{\prime }+\bar{e}^{\prime }R_0^{n-2}}{A_0^2B_0^2} -R_0^{n-2}\left(\frac{\bar{e}A_0^{\prime }}{A_0}+R_0^{\prime }\frac{\bar{b}}{B_0}\!\right)\right]+ \left[\frac{-T}{2B_0^2}\left\lbrace \bar{e}n(n-1)R_0^{n-1}-\frac{\bar{b}}{B_0}(1-n)R_0^n\!\right\rbrace + \frac{\bar{e}n(n-1)}{A_0}\big\lbrace\! (n-2)\right.\nonumber\\ &&\left.\left.\times R_0^{n-3}\dot{T}+R_0^{n-2}\ddot{T}\right\rbrace +\frac{Tn(n-1)}{B_0^2} \left\lbrace \bar{e}(n-2)R_0^{n-3}R_0^{\prime }+R_0^{n-2}\left(\bar{e}^{\prime }-\frac{\bar{b}}{B_0}R_0^{\prime }\right) \left(\frac{A_0^{\prime }}{A_0}+\frac{2}{r}\right)+R_0^{n-2}R_0^{\prime } \left\lbrace \left(\frac{\bar{a}}{A_0}\right)^{\prime }+2\left(\frac{\bar{c}}{r}\right)^{\prime }\right\rbrace \right\rbrace \right]_{,r}\nonumber\\ &&+ \,\bar{e}n(n-1)[(n-2)R_0^{n-3}\dot{T}+R_0^{n-2}\ddot{T}]\frac{A_0^{\prime }}{A_0^3B_0^2} +T\frac{A_0^{\prime }}{A_0B_0^4}n(n-1) \left[\left(\frac{\bar{a}^{\prime }}{A_0^{\prime }}-\frac{\bar{a}}{A_0}-\frac{\bar{b}}{B_0}\right)\right. \left\lbrace (n-2)R_0^{n-3}R_0^{\prime }+(n-2)R_0^{n-3}R_0^{\prime }-R_0^{n-2}\right.\nonumber\\ &&\left.\times \, R_0^{\prime }\left(\frac{A_0^{\prime }}{A_0} -\frac{B_0^{\prime }}{B_0}\right)\right\rbrace +\bar{e}(n-2)(n-3) R_0^{n-4}R_0^{\prime }+\bar{e}^{\prime \prime }(n-2)R_0^{n-3}+\bar{e}(n-2)(n-3)R_0^{n-4}R_0^{\prime \prime }-R_0^{n-3}\left\lbrace \bar{e}(n-2) -\left(\bar{e}^{\prime }-\frac{\bar{b}}{B_0}R_0^{\prime }\right)\right\rbrace \nonumber\\ &&\left.\times \, \left(\frac{A_0^{\prime }}{A_0}-\frac{B_0^{\prime }}{B_0}\right) +R_0^{n-2}R_0^{\prime }\left\lbrace \left(\frac{\bar{a}}{A_0}\right)^{\prime }-\left(\frac{\bar{b}}{B_0}\right)^{\prime }\right\rbrace \right]-\bar{e}nR_0^{n-1} + \frac{2B_0^{\prime }}{B_0^3}\frac{\bar{e}n(n-1)}{A_0^2}\lbrace (n-2)R_0^{n-3}\dot{T}+R_0^{n-2}\ddot{T}\rbrace -\frac{\bar{a}}{A_0}(1-n)R_0^n\nonumber\\ &&+\, T\frac{2B_0^{\prime }}{B_0^3}\left[\left(\frac{\bar{b}^{\prime }}{B_0^{\prime }}-\frac{3\bar{b}}{B_0}\right) \left\lbrace \frac{1}{2A_0^2}(1-n)R_0^n-\frac{n(n-1)}{B_0^2}R_0^{\prime }R_0^{n-2} \left(\frac{A_0^{\prime }}{A_0}-\frac{B_0^{\prime }}{B_0}\right)\right\rbrace -\frac{n(n-1)}{B_0^2}\left\lbrace R_0^{\prime }R_0^{n-2} \left\lbrace \left(\frac{\bar{a}}{A_0}\right)^{\prime }+\left(\frac{\bar{b}}{B_0}\right)^{\prime }\right\rbrace \right.\right.\nonumber\\ &&\left.\left.+\,\left(\frac{A_0^{\prime }}{A_0}+\frac{B_0^{\prime }}{B_0}\right)\left\lbrace \bar{e}(n-2) R_0^{n-3}R_0^{\prime }+R_0^{n-2} \left(\bar{e}^{\prime }-\frac{\bar{b}}{B_0}R_0\right)\right\rbrace \right\rbrace \right]-\frac{2 \dot{T}\bar{c}}{A_0^2r}\frac{n(n-1)}{B_0^2}R_0^{\prime }R_0^{n-2}\left(\frac{A_0^{\prime }}{A_0}-\frac{1}{r}\right). \end{eqnarray}
(A7)
Quantities in the Newtonian regime
\begin{eqnarray} \bar{q}(N)&=&\frac{-nR_0^{n-1}\sqrt{\omega }\bar{e}^{\sqrt{\omega }t}}{\kappa } \left[\frac{2}{r}\left(\frac{\bar{c}}{r}-\bar{c}^{\prime }\right) +\frac{(n-1)}{R_0^{n-1}}\left\lbrace (n-2)R_0^{n-3}\bar{e}R_0^{\prime }+R_0^{n-2}\bar{e}^{\prime }-\frac{\bar{c}}{r}R_0^{n-2}R_0^{\prime }\right\rbrace \right], \end{eqnarray}
(A8)
\begin{eqnarray} {\bar{\rho }}(N)&=&\frac{nR_0^{n-1}\bar{e}^{\sqrt{\omega }t}}{\kappa }\left[\frac{2}{r}\left(\frac{\bar{c}}{r}-\bar{c^{\prime }}\right) +\frac{n-1}{R_0^{n-1}} \left\lbrace (n-2)\bar{e}R_0^{n-3}R_0^{\prime }-\left(\frac{\bar{c}}{r}-\bar{e}^{\prime }\right)R_0^{n-2}R_0^{\prime }\right\rbrace \right]+\frac{n\bar{e}^{\sqrt{\omega }t}}{\kappa } \left[R_0^{n-1}\left\lbrace \frac{2}{r}\left(\frac{\bar{c}}{r}-\bar{c^{\prime }}\right) +\frac{n-1}{R_0^{n-1}} \right.\right.\nonumber\\ &&\left.\left.\times \lbrace (n-2)R_0^{n-3}\bar{e}R_0^{\prime }- \left(\frac{\bar{c}}{r}-\bar{e}^{\prime }\right)R_0^{n-2}\rbrace \right\rbrace \right]_{,r} +\rho _{0}\frac{3\bar{c}}{r}\bar{e}^{\sqrt{\omega} }+\frac{1}{\kappa }{F_{1{\rm p}}e^{\sqrt{\omega }t},} \end{eqnarray}
(A9)
\begin{eqnarray} \nonumber\\ F_{2{\rm p}(N)}&=&-\omega e^{\sqrt{\omega }t}n(n-1)\left[(n-2)R_0^{n-3}R_0^{\prime }+\bar{e}^{\prime }R_0^{n-2} -R_0^{n-2}R_0^{\prime }\frac{\bar{c}}{r}\right]+ \left[\frac{e^{\sqrt{\omega }t}}{2}\left\lbrace \bar{e}n(n-1)R_0^{n-1}-\frac{\bar{c}}{r}(1-n)R_0^n\right\rbrace - {\bar{e}n(n-1)}\big\lbrace (n-2)\right.\nonumber\\ &&\left.\left.\times R_0^{n-3}\sqrt{\omega }e^{\sqrt{\omega }t}-R_0^{n-2}\omega e^{\sqrt{\omega }t}\right\rbrace +e^{\sqrt{\omega }t}n(n-1) \left\lbrace \bar{e}(n-2)R_0^{n-3}R_0^{\prime }+ R_0^{n-2}\left(\bar{e}^{\prime }-\frac{\bar{c}}{r}R_0^{\prime }\right) \frac{2}{r}+R_0^{n-2}R_0^{\prime } \left\lbrace \bar{a}^{\prime }+2\left(\frac{\bar{c}}{r}\right)^{\prime }\right\rbrace \right\rbrace \right]_{,r}\nonumber\\ &&-\frac{2 \sqrt{\omega }e^{\sqrt{\omega }t}\bar{c}}{r^2}n(n-1)R_0^{\prime }R_0^{n-2}, \end{eqnarray}
(A10)
where F1p(N) in equation (A9) is given by
\begin{eqnarray} F_{1{\rm p}(N)}&=&n(n-1)[(n-2)R_0^{n-3}R_0^{\prime }+\bar{e}^{\prime }]-\frac{1}{2}\left[n(n-1)R_0^{n-1}-\bar{a}(1-n)R_0^n\right]-2\left(\bar{a}+\frac{\bar{c}}{r}\right) \left[n(n-1)(n-2)R_0^{n-3}\left(R_0^{\prime 2}+R_0^{\prime \prime }\right)\right]\nonumber\\ &&+n(n-1)(n-2) \left[\bar{e}(n-3)R_0^{n-4}R_0^{\prime }+2\bar{e}^{\prime }R_0^{n-3}+\bar{e}^{\prime \prime }R_0^{n-3}+\bar{e}(n-3)R_0^{n-4}R_0^{\prime \prime }\right]-n(n-1)\left[R_0^{\prime }R_0^{n-2}\left\lbrace \left(\frac{\bar{c}}{r}\right)^{\prime } \right.\right.\nonumber\\ &&-\left.\left. \left(\frac{2}{r}\right)^{\prime }\right\rbrace +\frac{1}{r}\left\lbrace \bar{e}(n-2)R_0^{n-3}R_0^{\prime }+ \left(\bar{e}^{\prime }-\frac{\bar{c}}{r}R_0^{\prime }\right)R_0^{n-2}\right\rbrace \right] +\bar{a}\left[\frac{1}{2}(1-n)R_0^n+n(n-1)(n-2)R_0^{n-3}(R_0^{\prime 2}+R_0^{\prime \prime })\right.\nonumber\\ &&+\frac{2n}{r} (n-1)R_0^{n-2}R_0^{\prime }\Big]+\frac{\bar{c}}{r}n(n-1)\left[(n-2)R_0^{n-3}\left(R_0^{\prime 2}+R_0^{\prime \prime }\right)\right]- \frac{2R_0^{n-2}R_0^{\prime }}{r^2}+\left[(n-2)\bar{e}R_0^{n-3}R_0^{\prime }+\bar{e}^{\prime }R_0^{n-2}- \frac{\bar{c}}{r}R_0^{\prime }R_0^{n-2}\right]\frac{n(n-1)}{r}.\nonumber\\ \end{eqnarray}
(A11)
R2 gravity
\begin{eqnarray} {\bar{q}}(N)&=&-\frac{2R_0e^{\sqrt{\omega }t}}{\kappa }\left[\frac{2}{r}\left(\frac{\bar{c}}{r}-\bar{c}^{\prime }\right) +R_0^{-1}\left\lbrace \bar{e}-\frac{\bar{c}}{r}R_0^{\prime }\right\rbrace \right], \end{eqnarray}
(A12)
\begin{eqnarray} {{\bar{\rho }}}(N)&=&\frac{2R_0e^{\sqrt{\omega }t}}{\kappa }\left[\frac{2}{r}\left(\frac{\bar{c}}{r}-\bar{c^{\prime }}\right) +\frac{2}{R_0}\left(\bar{e}^{\prime }-\frac{\bar{c}}{r}R_0^{\prime }\right)\right] +\rho _0\frac{3\bar{c}}{r}e^{\sqrt{\omega }t}+ \frac{1}{\kappa }F_{1{\rm p}}e^{\sqrt{\omega }t}-\frac{2e^{\sqrt{\omega }t}}{\kappa }\left[\frac{2R_0}{r}\left(\frac{\bar{c}}{r}-\bar{c^{\prime }}\right) +\bar{e}^{\prime }-\frac{\bar{c}}{r}R_0^{\prime }\right]_{,r}, \end{eqnarray}
(A13)
\begin{eqnarray} F_{2{\rm p}(N)}&=&-2\omega e^{\sqrt{\omega }t}\left(\bar{e}^{\prime }-R_0^{\prime }\frac{\bar{c}}{r}\right)+ \frac{1}{2}\left[e^{\sqrt{\omega }t}\left(2\bar{e}R_0+\frac{\bar{c}}{r}R_0^2\right)-2\bar{e}\omega e^{\sqrt{\omega }t}-\frac{4e^{\sqrt{\omega }t}}{r} \left(\bar{e}^{\prime }-\frac{\bar{c}}{r}R_0^{\prime }\right)+R_0^{\prime }\left\lbrace \bar{a}^{\prime }+2\left(\frac{\bar{c}}{r}\right)^{\prime }\right\rbrace \right]_{,r} -\frac{4\sqrt{\omega }e^{\sqrt{\omega }t}}{r^2}R_0^{\prime }, \nonumber\\ \end{eqnarray}
(A14)
\begin{eqnarray} F_{1{\rm p}(N)}&=&2\bar{e}^{\prime }-R_0^{\prime }-\frac{\bar{a}}{2}R_0^{2}+2R_0^{\prime }\left\lbrace \left(\frac{\bar{c}}{r}\right)^{\prime }+\frac{2}{r^2}\right\rbrace +\frac{\bar{2}}{r}\left(\bar{e}^{\prime }-\frac{\bar{c}}{r}R_0^{\prime }\right)-\frac{\bar{a}}{2R_0^2}+\frac{4\bar{a}}{r}R_0^{\prime }+ \frac{4\bar{c}}{r^3}R_0^{\prime }+\frac{2\bar{e}^{\prime }}{r}-\frac{2\bar{c}}{r^2}R_0^{\prime }. \end{eqnarray}
(A15)

|$\!\!\frac{1}{R}$|gravity

\begin{eqnarray} {\bar{q}}(N)&=&\frac{R_0^{-2}e^{\sqrt{\omega }t}}{\kappa }\left[\frac{2}{r}\left(\frac{\bar{c}}{r}-\bar{c}^{\prime }\right) +\frac{-2}{R_0^{-2}}\left\lbrace -3R_0^{-4}\bar{e}R_0^{\prime }+R_0^{-3}\bar{e}^{\prime }-\frac{\bar{c}}{r}R_0^{-3}R_0^{\prime }\right\rbrace \right], \end{eqnarray}
(A16)
\begin{eqnarray} {{\bar{\rho }}}(N)&=&\frac{R_0^{-2}e^{\sqrt{\omega }t}}{\kappa }\left[\frac{2}{r}\left(\frac{\bar{c}}{r}-\bar{c^{\prime }}\right) +\frac{-2}{R_0^{-2}} \left\lbrace -3\bar{e}R_0^{-4}R_0^{\prime }-\left(\frac{\bar{c}}{r}-\bar{e}^{\prime }\right)R_0^{-3}\right\rbrace \right]-\frac{e^{\sqrt{\omega }t}}{\kappa } \left[R_0^{-2}\left\lbrace \frac{2}{r}\left(\frac{\bar{c}}{r}-\bar{c^{\prime }}\right) \right.\right.\nonumber\\ &&-\left.\left.\frac{2}{R_0^{-2}}\Big\lbrace\! -3R_0^{-4}\bar{e}R_0^{\prime }-\left(\frac{\bar{c}}{r}-\bar{e}^{\prime }\right)R_0^{-3}\rbrace \right\rbrace \right]_{,r}+ \rho _0\frac{3\bar{c}}{r}e^{\sqrt{\omega }t}+\frac{1}{\kappa }F_{1{\rm p}}e^{\sqrt{\omega }t}, \end{eqnarray}
(A17)
\begin{eqnarray} F_{2{\rm p}(N)}&=&-2\omega e^{\sqrt{\omega }t}\left[-3R_0^{-4}R_0^{\prime }+\bar{e}^{\prime }R_0^{-3} -R_0^{-3}R_0^{\prime }\frac{\bar{c}}{r}\right]+e^{\sqrt{\omega }t}\left[\left(\bar{e}R_0^{-2}- \frac{\bar{c}}{r}R_0^{-1}\right)+ 2\bar{e}\left(3R_0^{-4}\sqrt{\omega }-R_0^{-3}\omega \right) \right.\nonumber \\ &&\left.-2\left\lbrace -3\bar{e}R_0^{-4}R_0^{\prime }+\frac{2R_0^{-3}}{r}\left(\bar{e}^{\prime }-\frac{\bar{c}}{r}R_0^{\prime }\right) +R_0^{-3}R_0^{\prime }\left\lbrace \bar{a}^{\prime }+2\left(\frac{\bar{c}}{r}\right)^{\prime }\right\rbrace \right\rbrace \right]_{,r} -\frac{4\sqrt{\omega }e^{\sqrt{\omega }t}\bar{c}}{r^2}R_0^{\prime }R_0^{-3}, \end{eqnarray}
(A18)
\begin{eqnarray} \nonumber\\ F_{1{\rm p}(N)}&=&2(\bar{e}^{\prime }-3R_0^{-4}R_0^{\prime })-R_0^{-2}+\bar{a}R_0^{-1}+2\left(\bar{a}+\frac{\bar{c}}{r}\right) R_0^{-4}\left(R_0^{\prime 2}+R_0^{\prime \prime }\right)-6(-4\bar{e}R_0^{-5}R_0^{\prime }+2\bar{e}^{\prime }R_0^{-4} +\bar{e}^{\prime \prime }R_0^{-4}-4\bar{e}R_0^{-5}R_0^{\prime \prime }) \nonumber\\ && -\, 2\left[R_0^{\prime }R_0^{-3} \left\lbrace \left(\frac{\bar{c}}{r}\right)^{\prime }-\left(\frac{2}{r}\right)^{\prime }\right\rbrace +\frac{1}{r}\left\lbrace -3\bar{e}R_0^{-4}R_0^{\prime }+ \left(\bar{e}^{\prime }-\frac{\bar{c}}{r}R_0^{\prime }\right)R_0^{-3}\right\rbrace \right]+ \bar{a}\left[-R_0^{-1}-6R_0^{-4}\left(R_0^{\prime 2}+R_0^{\prime \prime }\right)+\frac{4}{r}R_0^{-3}R_0^{\prime }\right] \nonumber\\ &&+\,\frac{2\bar{c}}{r}\left[-3R_0^{-4}\left(R_0^{\prime 2}+ R_0^{\prime \prime }\right)-\frac{2R_0^{-3}R_0^{\prime }}{r^2}\right]+\frac{2}{r}\left[-3\bar{e}R_0^{-4}R_0^{\prime }+\bar{e}^{\prime }R_0^{-3}- \frac{\bar{c}}{r}R_0^{\prime }R_0^{-3}\right]. \end{eqnarray}
(A19)