# A Non-static Quantum inspired Spacetime in $f(R)$ Gravity: Gravity’s Rainbow

Prabir Rudra
_Department of Mathematics, Asutosh College, Kolkata-700 026, India._
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this note we explore a non-static spacetime in quantum regime in the
background of $f(R)$ gravity. The time dependent Vaidya metric which
represents the spacetime of a radiating body like star is studied in an energy
dependent gravity’s rainbow, which is a UV completion of General Relativity.
In our quest we have used gravitational collapse as the main tool. The focus
is to probe the nature of singularity (black hole or naked singularity) formed
out of the collapsing procedure. This is achieved via a geodesic study. For
our investigation we have considered two different models of $f(R)$ gravity,
namely the inflationary Starobinsky’s model and the power law model. Our study
reveals the fact that naked singularity is as good a possibility as black hole
as far as the central singularity is concerned. Via a proper fine tuning of
the initial data, we may realize both black hole or naked singularity as the
end state of the collapse. Thus this study is extremely important and relevant
in the light of the Cosmic Censorship hypothesis. The most important result
derived from the study is that gravity’s rainbow increases the tendency of
formation of naked singularities. We have also deduced the conditions under
which the singularity will be a strong or weak curvature singularity. Finally
in our quest to know more about the model we have performed a thermodynamical
study. Throughout the study we have obtained results which involve deviation
from the classical set-up. Such deviations are expected in a quantum evolution
and can be attributed to the quantum fluctuations that our model suffers from.
It is expected that this study will enhance our knowledge about quantization
of gravity and subsequently about the illusive theory of quantum gravity.

The author dedicates this work to the memory of the people who died due to
COVID-19

##  1 Introduction

More than a hundred years back, Einstein proposed his theory of general
relativity (GR) which revolutionized our idea of gravity. It provided us a
theory that could act as a powerful tool in our quest to know the universe.
With the passage of time and with extensive research glitches in the theory
started showing up. The major blow came at the turn of the last century when
the discovery of the accelerated expansion of the universe [1, 2] left GR
inconsistent at cosmological distances. Since then we have resorted to the
alternative techniques of modified gravity and dark energy to incorporate this
accelerated expansion in our theory of gravity. While the former deals with
the geometry of spacetime, the latter is concerned with the matter content of
the universe. Extensive reviews in modified gravity can be found in the
Refs.[3, 4, 5].

Many of such theories aim at modifying the linear function of scalar curvature
$R$ from its special form in GR to a more generic form. $f(R)$ gravity is one
such attempt where the gravitational lagrangian of GR, $\mathcal{L}_{GR}=R$ is
replaced by an analytic function of $R$ i.e. $\mathcal{L}_{f(R)}=f(R)$.
Choosing a suitable function for $f(R)$, one can explore the non-linear
effects of the scalar curvature on the evolution of the universe. Extensive
reviews in $f(R)$ gravity can be found in the Refs.[6, 7]. Viability of $f(R)$
dark energy models have been studied in ref.[8], where the $f(R)$ models with
a power law of $R$ has been ruled out. Author of Ref.[9] studied the interplay
between $f(R)$ theories and scalar-tensor theories via the Palatini formalism.
Formation of large scale structure in $f(R)$ gravity was studied in Ref.[10].
A reconstruction scheme for $f(R)$ theories was explored in Ref.[11]. Various
other studies related to $f(R)$ gravity can be found in Refs.[12, 13, 14, 15].

For a long time we have been searching for a theory of gravity that will be
consistent at all length scales and at all energy levels. Such a theory is
termed as the Theory of Everything (TOE). It is understood that this would
eventually boil down to a theory of quantum gravity where GR will be
reconciled to the theory of quantum mechanics (QM). Till now there have been a
few proposals for such a theory namely, Loop quantum theory [16, 17], String
theory [18, 19], Horava-Lifshitz (HL) gravity [20, 21], etc. The UV completion
of GR in the limit that GR is recovered in IR limit has eventually led to the
development of the Horava-Lifshitz gravity. An alternative mechanism of UV
completion of GR was proposed by Magueijo and Smolin in the Ref.[22] where the
geometry of spacetime depends on the energy of the test particle. This theory
is termed as Gravity’s Rainbow. Although the conceptual basis of HL gravity
and gravity’s rainbow are quite different yet they aim to pursue similar
ideas. Both the theories resorts to the modification of the usual energy-
momentum dispersion relations of the special theory of relativity in the UV
limit, such that in the IR limit the usual relations are retrieved. Following
this, the authors in Ref.[23] tried to bridge the two theories from the
conceptual background.

We know that in GR the usual energy-momentum relations are governed by the
Lorentz symmetry and so any modifications to these in the UV limit will
directly imply the violation of the symmetry. In fact different quantum
gravity approaches have shown that at high energy scales (UV limit) Lorentz
symmetry breaks down [24, 25, 26, 27]. This breakdown is expected to occur in
models like string field theory [28], discrete spacetime [29], spacetime foam
[30], non-commutative geometry [31], spin network in Loop quantum gravity
[32], etc. Magueijo and Smolin in the Ref.[33] proposed the theory of Doubly
Special Relativity (DSR), where Einstein’s Special Theory of Relativity (STR)
is generalized for spacetimes with high energy, i.e. energies compared to
Planck energy $(E_{P}=10^{19}GeV)$. In this theory Planck energy joins the
speed of light as an invariant quantity. Nevertheless this formulation was
achieved at the cost of the violation of the Lorentz symmetry. The
modifications to the field theory also suggested corresponding modification to
the equivalence principle which determines how DSR can be embedded in GR.
Introducing non-zero curvature in DSR, we get what can be called the Doubly
General Relativity (DGR) or Gravity’s Rainbow. In this theory the
modifications to the usual energy-momentum dispersion relations
$(E^{2}-p^{2}=m^{2})$ as discussed above is introduced via energy dependent
rainbow functions $\mathcal{F}(E)$ and $\mathcal{G}(E)$ as given below,

$E^{2}\mathcal{F}^{2}(E)-p^{2}\mathcal{G}^{2}(E)=m^{2}$ (1)

Here $E=E_{s}/E_{P}$, where $E_{s}$ is the maximum energy that a probe in the
spacetime can support, and $E_{P}$ is the Planck energy. From the conceptual
background of the theory it is obvious that $E_{s}$ cannot exceed $E_{P}$. The
choice of the rainbow functions should be such that they respect the
correspondence principle, i.e. at the IR limit we should be able to recover
the usual energy-momentum dispersion relations of classical GR from the
Eqn.(1). This means that the rainbow functions are required to satisfy the
relations,

$\lim\limits_{E_{s}/E_{P}\to
0}\mathcal{F}(E)=1,\qquad\lim\limits_{E_{s}/E_{P}\to 0}\mathcal{G}(E)=1.$ (2)

In the above relations the limit $E_{S}/E_{P}\to 0$ corresponds to a spacetime
with low energy where GR becomes dominant. It is expected that as $E_{s}\to
E_{P}$, effects of GR gradually fades away and quantum gravity effects become
more and more dominant. The rainbow functions are motivated from various
theoretical and phenomenological considerations in literature. Using the
results from loop quantum gravity and $\kappa$-Minkowski noncommutative
spacetime these functions have been proposed as [34, 35],

$\mathcal{F}(E_{s}/E_{P})=1,~{}~{}~{}~{}\mathcal{G}(E_{s}/E_{P})=\sqrt{1-\eta\frac{E_{s}}{E_{P}}}$
(3)

where $\eta$ is a constant. The modified dispersion from constant velocity of
light motivates the following rainbow functions [33]

$\mathcal{F}(E_{s}/E_{P})=\mathcal{G}(E_{s}/E_{P})=\frac{1}{1-a_{1}\frac{E_{s}}{E_{P}}}$
(4)

where $a_{1}$ is a constant. Moreover the observations of the hard spectra
from gamma ray bursters have been used motivate the rainbow functions [30],

$\mathcal{F}(E_{s}/E_{P})=\frac{e^{a_{2}E_{s}/E_{P}}-1}{a_{2}E_{s}/E_{P}},~{}~{}~{}~{}~{}~{}~{}\mathcal{G}(E_{s}/E_{P})=1$
(5)

where $a_{2}$ is a constant. In this study it makes sense to use the rainbow
functions given by eqn.(3) because the choice is motivated from a theory of
quantum gravity (loop quantum gravity).

The metric in gravity’s rainbow is written as

$g^{\mu\nu}(E)=\eta^{ab}e^{\mu}_{a}(E)e^{\nu}_{b}(E).$ (6)

where $e_{0}=\mathcal{F}^{-1}(E_{s}/E_{P})\tilde{e_{0}}$ and
$e_{i}=\mathcal{G}^{-1}(E_{s}/E_{P})\tilde{e_{i}}$. Here the tilde quantities
refer to the frame fields which are independent of energy thus corresponding
to the geometry explored by a low energy quanta. Eqn.(6) actually represents a
family of energy dependent metrics thus forming a rainbow. In such a spacetime
each test particle with different energy will probe a different geometry thus
following different geodesics. Due to its quantum gravity background gravity’s
rainbow have been studied extensively in recent times [36, 37, 38, 39, 40, 41,
42, 43].

Gravitational collapse is an astrophysical phenomenon that plays a central
role in the structure formation process of the universe. For this reason
gravitational collapse has been a field of interest for astrophysicists over
the years. It all started with Oppenheimer and Snyder who studied the collapse
of a dust cloud with a static Schwarzschild exterior and Friedmann like
interior [44]. Subsequently the collapse of spherically symmetric
inhomogeneous distribution of dust was studied by Tolman [45] and Bondi [46].
In 1969, Roger Penrose in a phenomenal paper [47] argued that any types of
cosmological singularity is bound to be shrouded by an event horizon thus
making it a Black Hole (BH). This proposal is popularly known in literature as
the Cosmic Censorship Hypothesis (CCH). But this idea has subsequently been
questioned in the absence of any form of rigorous proof or observational
evidence. This naturally led people to search for collapsing models that can
yield singularities which are uncensored, commonly known as Naked
Singularities (NS) [48, 49, 50, 51, 52, 53, 54]. Such an uncensored
singularity will be really interesting to study because it will leak
information that is generally hidden behind the event horizon of a BH. Not
only the solution of information paradox but a proper knowledge of NS will
eventually help us understand the gravity quanta thus allowing us to formulate
the illusive satisfactory theory of quantum gravity.

In 1951, P.C. Vaidya formulated a relativistic line element representing the
field of radiation for a non-static mass [55]. This is a generalization of the
Schwarzschild solution for non-static mass. Schwarzschild’s external solution
can represent the gravitational field of a cold dark body with a constant
mass. So it is obvious that the application of this solution to describe the
sun’s gravitational field should only be considered as approximate. Vaidya’s
metric precisely solved this problem by successfully describing the spacetime
of Sun. In fact it can represent the spacetime of any radiating mass, such as
a star. This is why it is sometimes called the radiating or shining
Schwarzschild metric. It should be mentioned here that this metric, if
expressed in radiation coordinates, differs from the Schwarzschild metric only
in that the constant mass parameter $m$ is replaced by a function of retarded
time. Notable studies using gravity’s rainbow in Vaidya’s metric can be found
in [56, 57, 58]. Other important studies on Vaidya spacetime can be found in
[59, 60, 61, 62, 63].

From the above discussion we feel the need to explore the non-static radiating
Vaidya spacetime in a quantum regime. This can be achieved by introducing
rainbow deformations in the Vaidya spacetime in the background of a gravity
theory. $f(R)$ gravity being the simplest and most obvious theory of modified
gravity at least at the mathematical level, we are inclined to consider it in
the background of our model. It is expected that the radiating star
represented by Vaidya spacetime will yield very interesting results in the
quantum limit. The model will be investigated via a study of gravitational
collapse, which is a very important astrophysical phenomenon. The study will
be eventually complimented by a thermodynamical study in the said model. The
search for the illusive theory of quantum gravity is a motivation for the
present work and we hope that our investigation will be a step towards the
better understanding of the nature of gravity quanta. The paper is organized
as follows. In Sec.2, the field equations for the rainbow deformed Vaidya
spacetime in $f(R)$ gravity are derived and a solution is obtained. Sec.3 is
dedicated to the study of gravitational collapse in the system. In Sec.4 we
have studied the thermodynamical properties of the system in detail. Finally
the paper is concluded with a conclusion and discussion in Sec.5.

##  2 Rainbow deformed Vaidya spacetime in $f(R)$ gravity

The Einstein-Hilbert action of GR is given by,

$S_{EH}=\frac{1}{2\kappa}\int d^{4}x\sqrt{-g}R$ (7)

where $\kappa\equiv 8\pi G$, $G$ is the gravitational constant, $g$ is the
determinant of the metric and $R$ is the Ricci scalar (we have considered
$c=1$). We replace the Ricci scalar, $R$ in the above action by a generalized
function of $R$ to get the action for $f(R)$ gravity,

$S=\frac{1}{2\kappa}\int d^{4}x\sqrt{-g}f(R)$ (8)

Beginning from the action (8) and adding a matter term $S_{M}$, the total
action for $f(R)$ gravity takes the form,

$S=\frac{1}{2\kappa}\int d^{4}x\sqrt{-g}f(R)+S_{M}(g_{\mu\nu},\psi)$ (9)

where $\psi$ collectively denotes the matter fields. Taking variation with
respect to the metric we get the field equations as,

$f^{\prime}(R)R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}f(R)+\left(g_{\mu\nu}\Box-\nabla_{\mu}\nabla_{\nu}\right)f^{\prime}(R)=\kappa
T_{\mu\nu}$ (10)

where $T_{\mu\nu}$ is given by,

$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\frac{\delta S_{M}}{\delta g^{\mu\nu}}$ (11)

In the above equations a prime denotes differentiation with respect to the
argument, i.e. $R$. $\nabla_{\mu}$ denotes covariant derivative associated
with the Levi-Civita connection of the metric and
$\Box\equiv\nabla^{\mu}\nabla_{\mu}$ is the D’Alembertian operator. The field
equations given by equation (10) can also be written in the following form,

$G_{\mu\nu}=R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R=\frac{\kappa
T_{\mu\nu}}{f^{\prime}(R)}+g_{\mu\nu}\frac{\left[f(R)-Rf^{\prime}(R)\right]}{2f^{\prime}(R)}+\frac{\left[\nabla_{\mu}\nabla_{\nu}f^{\prime}(R)-g_{\mu\nu}\Box
f^{\prime}(R)\right]}{f^{\prime}(R)}$ (12)

The Vaidya metric in the advanced time coordinate system is given by,

$ds^{2}=f(t,r)dt^{2}+2dtdr+r^{2}\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right)$ (13)

where $f(t,r)=-\left(1-\frac{m(t,r)}{r}\right)$. It is obvious that here we
have used the units $G=c=1$.

Introducing rainbow deformations in the above metric we get[22],

$ds^{2}=-\frac{1}{\mathcal{F}^{2}(E)}\left(1-\frac{m(t,r)}{r}\right)dt^{2}+\frac{2}{\mathcal{F}(E)\mathcal{G}(E)}dtdr+\frac{r^{2}}{\mathcal{G}^{2}(E)}\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right)$ (14)

The total energy momentum tensor of the field equation (12) is given by the
following sum,

$T_{\mu\nu}=T_{\mu\nu}^{(n)}+T_{\mu\nu}^{(m)}$ (15)

where $T_{\mu\nu}^{(n)}$ and $T_{\mu\nu}^{(m)}$ are the contributions from the
Vaidya null radiation and perfect fluid respectively defined as,

$T_{\mu\nu}^{(n)}=\sigma l_{\mu}l_{\nu}$ (16)

and

$T_{\mu\nu}^{(m)}=(\rho+p)(l_{\mu}\eta_{\nu}+l_{\nu}\eta_{\mu})+pg_{\mu\nu}$
(17)

where $\rho$ and $p$ are the energy density and pressure for the perfect fluid
and $\sigma$ is the energy density corresponding to Vaidya null radiation. In
the co-moving co-ordinates ($t,r,\theta_{1},\theta_{2},...,\theta_{n}$), the
two eigen vectors of energy-momentum tensor namely $l_{\mu}$ and $\eta_{\mu}$
are linearly independent future pointing null vectors having components

$l_{\mu}=(1,0,0,0)~{}~{}~{}~{}and~{}~{}~{}~{}\eta_{\mu}=\left(\frac{1}{2}\left(1-\frac{m}{r}\right),-1,0,0\right)$
(18)

and they satisfy the relations

$l_{\lambda}l^{\lambda}=\eta_{\lambda}\eta^{\lambda}=0,~{}l_{\lambda}\eta^{\lambda}=-1$
(19)

Now we impose rainbow deformations on the linearly independent future pointing
null vectors $l_{\mu}$ and $\eta_{\mu}$ such that we get,

$l_{\mu}=\left(\frac{1}{\mathcal{F}(E)},0,0,0\right)~{}~{}~{}~{}~{}\&~{}~{}~{}\eta_{\mu}=\left(\frac{1}{2\mathcal{F}(E)}\left(1-\frac{m(t,r)}{r}\right),-\frac{1}{\mathcal{G}(E)},0,0\right)$
(20)

satisfying the following conditions (19). Therefore, the non-vanishing
components of the total energy-momentum tensor will be as follows

$\displaystyle
T_{00}=\frac{\sigma}{\mathcal{F}^{2}(E)}+\frac{\rho}{\mathcal{F}^{2}(E)}\left(1-\frac{m(t,r)}{r}\right),$
$\displaystyle~{}~{}T_{01}=-\frac{\rho}{\mathcal{F}(E)\mathcal{G}(E)},$
$\displaystyle T_{22}=\frac{pr^{2}}{\mathcal{G}^{2}(E)},$
$\displaystyle~{}~{}T_{33}=\frac{pr^{2}\sin^{2}\theta}{\mathcal{G}^{2}(E)}$
(21)

Here we consider matter in the form of perfect barotropic fluid given by the
equation of state

$p=\omega\rho$ (22)

where ’$\omega$’ is the barotropic parameter.

The non-vanishing components of the Einstein tensors are given by,

$\displaystyle
G_{00}=\frac{\mathcal{G}(E)}{r^{3}\mathcal{F}^{2}(E)}\left\\{\mathcal{G}(E)\left(r-m\right)m^{\prime}+\mathcal{F}(E)r\dot{m}\right\\},$
$\displaystyle~{}~{}~{}~{}G_{01}=G_{10}=-\frac{\mathcal{G}(E)m^{\prime}}{r^{2}\mathcal{F}(E)}$
$\displaystyle G_{22}=-\frac{1}{2}rm^{\prime\prime},$
$\displaystyle~{}~{}G_{33}=-\frac{1}{2}rm^{\prime\prime}~{}\sin^{2}\theta$
(23)

where dot and dash represents derivatives with respect to time $t$ and radial
coordinate $r$ respectively. For this system the Ricci scalar becomes,

$R=\frac{\mathcal{G}^{2}(E)}{r^{2}}\left(2m^{\prime}+rm^{\prime\prime}\right)$
(24)

###  2.1 Field Equations

Here we report the computed field equation for $f(R)$ gravity’s rainbow in the
time dependent Vaidya spacetime. In the study that follows we have made use of
the $(11)$ component of the field equations because of its simplicity and so
we report only that component in this section. The other components of the
field equations are reported in the appendix section.

The (11)-component of field equations is given by,

$f^{\prime\prime}(R)r^{2}\left(12m^{\prime}-6m^{\prime\prime}r+m^{iv}r^{3}\right)+f^{\prime\prime\prime}(R)\mathcal{G}^{2}(E)\left\\{-4m^{\prime}+r\left(m^{\prime\prime}+m^{\prime\prime\prime}r\right)\right\\}^{2}=0$
(25)

where $m^{iv}$ denotes the fourth partial derivative of $m$ with respect to
$r$.

###  2.2 Solution of the system

To find a solution of the field equations we have to consider specific models
of $f(R)$ gravity. Here we consider two popular models, namely the
inflationary Starobinsky’s model and the power law model, and subsequently
find the solution of the system. The emergence of the cosmic structure from
the homogeneous and isotropic universe cannot be clearly explained by the
standard inflationary models, because the mechanism involved preserves the
homogeneity and isotropy at all times. A solution to this problem has been
proposed by Sudarsky et al. in Refs.[64, 65, 66, 67, 68, 69, 70, 71, 72],
where they have introduced the concept of self induced collapse hypothesis.
Here the collapse of the wave function of the inflaton mode is restricted to
occur during the inflationary period. So there are reasons to believe that
collapsing scenario in inflationary models can help us understand the quantum
evolution of universe. Moreover quantum corrections induced by rainbow
functions will be highly compatible with such a model. This motivates us to
choose the Starobinsky’s model for our study. The motivation for the power law
model is obvious as it is the most generic model of $f(R)$ gravity capable of
representing all the epochs of the universe by fine tuning the initial data.

####  2.2.1 Starobinsky’s Model

In this section we consider the popular Starobinsky’s inflationary model of
$f(R)$ gravity. The model is given as [73, 74]

$f(R)=R+aR^{2}$ (26)

where $a\geq 0$ is the only free parameter of the theory and has the
dimensions of $[mass]^{-2}$. So the parameter $a$ can also be written in the
form $a=1/M^{2}$, where mass becomes the free parameter of the gravity theory.
It can easily be seen that GR can be retrieved from the theory for $a=0$.
Using eqns. (24),(25) and (26) we get the following differential equation in
terms of the mass ’$m$’,

$2ar^{2}\left(r^{3}m^{iv}-6rm^{\prime\prime}+12m^{\prime}\right)=0$ (27)

We can see that for $a=0$, we get an identity which corresponds to the case
for GR. From the above equation we get,

$r^{3}m^{iv}-6rm^{\prime\prime}+12m^{\prime}=0,~{}~{}~{}a\neq 0$ (28)

It should be noted that the parameter ${}^{\prime}a^{\prime}$ does not appear
in the above differential equation. Solving the above differential equation we
get,

$m(t,r)=f_{1}(t)-\frac{f_{2}(t)}{r}+f_{3}(t)r^{3}+f_{4}(t)r^{4}$ (29)

where $f_{1}(t)$, $f_{2}(t)$, $f_{3}(t)$ and $f_{4}(t)$ are arbitrary
functions of time $t$. So for the Starobinsky’s model, the rainbow deformed
Vaidya spacetime in $f(R)$ gravity is given by,

$ds^{2}=-\frac{1}{\mathcal{F}^{2}(E)}\left[1-\frac{1}{r}\left(f_{1}(t)-\frac{f_{2}(t)}{r}+f_{3}(t)r^{3}+f_{4}(t)r^{4}\right)\right]dt^{2}+\frac{2}{\mathcal{F}(E)\mathcal{G}(E)}dtdr+\frac{r^{2}}{\mathcal{G}^{2}(E)}\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right)$ (30)

####  2.2.2 Power Law Model

Regarding power law model we would like to state that probably it is the most
generic model for any theory and has a wide range of applicability. Hence it
is not surprising that it will reduce to Starobinsky’s model for some physical
restrictions. But it should be mentioned here that one should be cautious in
using a generic form of power law model because for all values of the exponent
the model may not remain cosmologically viable or for all values the model may
not satisfy the standard observational tests. Amendola et al. [8] points out
such a discrepancy and rules out some basic models of f(R) theory. There are
other examples as well. In the light of these results any unhindered and free
use of the power law model will be worrying. However here our concern is to
keep our study as generic as possible and not bother about observational
compliances. Since our basic aim is to study gravitational collapse which is
an astrophysical phenomenon we have restricted our study to just astrophysical
implications.

We consider the following power law model [75, 76] of $f(R)$ gravity,

$f(R)=\alpha R^{n}~{},~{}~{}~{}~{}~{}n>0$ (31)

where $\alpha$ and $n>0$ are constants. Using eqns.(24),(25) and (31) we get
the following differential equation in terms of the mass ’$m$’,

$\displaystyle\mathcal{G}^{2}(E)\alpha~{}n\left(n-1\right)\left[\frac{\mathcal{G}^{2}(E)}{r^{2}}\left(rm^{\prime\prime}+2m^{\prime}\right)\right]^{n-3}\left[\left(n-2\right)\left(r^{2}m^{\prime\prime\prime}+rm^{\prime\prime}-4m^{\prime}\right)^{2}\right.$
$\left.+\left(rm^{\prime\prime}+2m^{\prime}\right)\left(r^{3}m^{iv}-6rm^{\prime\prime}+12m^{\prime}\right)\right]=0$
(32)

It is almost impossible to get a general solution for the above equation by
the known mathematical methods. So we may discuss some particular cases. It
can easily be seen that the equation becomes an identity for $n=0$ and $n=1$
and does not possess a viable solution. Moreover these are trivial cases for
the power law and are of very little interest. For $n=2$ this equation reduces
to the corresponding equation for the Starobinsky model discussed in the
previous section.

The equation can be separated into two component equations as given below:

$rm^{\prime\prime}+2m^{\prime}=0,~{}~{}~{}~{}~{}for~{}~{}n\geq 4$ (33)

OR

$\left(n-2\right)\left(r^{2}m^{\prime\prime\prime}+rm^{\prime\prime}-4m^{\prime}\right)^{2}+\left(rm^{\prime\prime}+2m^{\prime}\right)\left(r^{3}m^{iv}-6rm^{\prime\prime}+12m^{\prime}\right)=0$
(34)

For eqn.(33) the solution is obtained as,

$m(t,r)=f_{5}(t)-\frac{f_{6}(t)}{r}$ (35)

where $f_{5}(t)$ and $f_{6}(t)$ are arbitrary functions of time, $t$. Eqn.(34)
gives two alternative solutions only for $n=2$. The solutions are given by,

$m(t,r)=f_{7}(t)-\frac{f_{8}(t)}{r}$ (36)

OR

$m(t,r)=f_{9}(t)-\frac{f_{10}(t)}{r}+f_{11}(t)r^{3}+f_{12}(t)r^{4}$ (37)

where $f_{i}(t),~{}i=7,8,9,10,11,12$ are arbitrary functions of time.
Combining eqns.(36) and (37) we get (for $n=2$),

$m(t,r)=A\left(f_{7}(t)-\frac{f_{8}(t)}{r}\right)+B\left(f_{9}(t)-\frac{f_{10}(t)}{r}+f_{11}(t)r^{3}+f_{12}(t)r^{4}\right)$
(38)

where $A$ and $B$ are arbitrary constants. Since solution (35) and (38) are
valid for different values of $n$, they cannot be combined to get a single
expression for $m(t,r)$. So the expression for $m(t,r)$ becomes,

$m(t,r)=\left\\{\begin{array}[]{c}f_{5}(t)-\frac{f_{6}(t)}{r},~{}~{}~{}~{}when~{}n\geq
4\\\ \\\
A\left(f_{7}(t)-\frac{f_{8}(t)}{r}\right)+B\left(f_{9}(t)-\frac{f_{10}(t)}{r}+f_{11}(t)r^{3}+f_{12}(t)r^{4}\right),~{}~{}~{}~{}when~{}n=2\end{array}\right.$
(39)

Hence for the power law model, the rainbow deformed Vaidya spacetime in $f(R)$
gravity is given by,

$ds^{2}=\left\\{\begin{array}[]{c}-\frac{1}{\mathcal{F}^{2}(E)}\left[1-\frac{1}{r}\left(f_{5}(t)-\frac{f_{6}(t)}{r}\right)\right]dt^{2}+\frac{2}{\mathcal{F}(E)\mathcal{G}(E)}dtdr+\frac{r^{2}}{\mathcal{G}^{2}(E)}\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right),\\\ \\\ ~{}~{}~{}~{}when~{}n\geq 4\\\ \\\
-\frac{1}{\mathcal{F}^{2}(E)}\left[1-\frac{1}{r}\left\\{A\left(f_{7}(t)-\frac{f_{8}(t)}{r}\right)+B\left(f_{9}(t)-\frac{f_{10}(t)}{r}+f_{11}(t)r^{3}+f_{12}(t)r^{4}\right)\right\\}\right]dt^{2}+\frac{2}{\mathcal{F}(E)\mathcal{G}(E)}dtdr\\\
\\\ +\frac{r^{2}}{\mathcal{G}^{2}(E)}\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right),~{}~{}~{}~{}when~{}n=2\end{array}\right.$ (40)

##  3 Gravitational Collapse

In this section we aim to study the phenomenon of gravitational collapse of a
star modelled by the deformed Vaidya metric given in eqn.(14). We will use the
solutions obtained in eqns.(30) and (40) for the two different models
separately. The methodology to be employed to study the collapsing procedure
is the geodesic study. We will be interested in probing the existence of non-
spacelike radial geodesics emanating from the central singularity, that were
terminated in the past at the singularity $r=0$. A geodesic coming out from
the central singularity and interacting with the outside would imply the non
existence of the event horizon around the singularity. Such a singularity is
quite aptly termed as a naked singularity. Due to the absence of any form of
horizon an outside observer can receive information from such a singularity
and vice-versa. If only a single null geodesic escapes from the singularity,
it indicates the emission of a single wavefront and hence the singularity
would be visible (naked) only momentarily to a distant observer. This
singularity is locally naked. On the other hand if the NS is to be visible for
a finitely prolonged time period, a family of geodesics must leave the central
singularity, which will make the singularity globally naked. If our search for
the emanating geodesics yield negative result, then the singularity is bound
to be a BH and we will get yet another reason to upheld the cosmic censorship
hypothesis, which does not have a rigorous proof till date. On the contrary
there are quite a few works in literature [48, 53, 54] that supports the
formation of NS thus providing significant counterexamples for the CCH as
discussed earlier. In spite of these counterexamples, generically in classical
background, CCH holds good and a singularity is always censored but in this
work it is expected that the quantum nature of gravity will play its role and
support the formation of NS.

Considering the collapse to be spherical, we assume the physical radius of the
$r$-th shell of the star at time $t$ be $R(t,r)$. In the epoch $t=0$, we have
$R(0,r)=r$. If the collapse is inhomogeneous, then different shells may become
singular at different times. If the outgoing non-spacelike geodesics possess
well defined tangent at the singularity, $\frac{dR}{dr}$ will tend to a finite
limit as the geodesic approaches the singularity in the past along the
trajectories. As the trajectories reach the points $(t_{0},r)=(t_{0},0)$, the
singularity occurs, which is given by $R(t_{0},0)=0$. Physically this
corresponds to the matter shells being crushed to zero radius, resulting in
the formation of the central singularity. If we trace back the trajectories of
the outgoing non-spacelike geodesics from this central singularity, it is
likely that they will terminate in the past at the singularity $(r=0,t=t_{0})$
where $R(t_{0},0)=0$. Therefore we should have $R\rightarrow 0$ as
$r\rightarrow 0$ [77].

The equation for outgoing radial null geodesics can be obtained from equation
(14) by putting $ds^{2}=0$ and $d\Omega_{2}^{2}=d\theta^{2}+\sin^{2}\theta
d\phi^{2}=0$ as

$\frac{dt}{dr}=\frac{2\mathcal{F}(E)}{\mathcal{G}(E)\left(1-\frac{m(t,r)}{r}\right)}.$
(41)

From the above expression it is quite clear that at $r=0,~{}t=0$ there is a
singularity of the above differential equation. Suppose we consider a
parameter $X=\frac{t}{r}$. Using this parameter we can study the limiting
behaviour of the function $X$ as we approach the singularity at $r=0,~{}t=0$
along the radial null geodesic. If we denote the limiting value by $X_{0}$
then using L’Hospital’s rule we have

$\displaystyle\begin{array}[]{c}X_{0}\\\ \\\
{}\hfil\end{array}\begin{array}[]{c}=\lim X\\\ \tiny t\rightarrow 0\\\ \tiny
r\rightarrow 0\end{array}\begin{array}[]{c}=\lim\frac{t}{r}\\\ \tiny
t\rightarrow 0\\\ \tiny r\rightarrow
0\end{array}\begin{array}[]{c}=\lim\frac{dt}{dr}\\\ \tiny t\rightarrow 0\\\
\tiny r\rightarrow
0\end{array}\begin{array}[]{c}=\lim\frac{2\mathcal{F}(E)}{\mathcal{G}(E)\left(1-\frac{m(t,r)}{r}\right)}\\\
\tiny t\rightarrow 0{}{}{}{}{}{}{}{}{}{}{}{}\\\ \tiny r\rightarrow
0{}{}{}{}{}{}{}{}{}{}{}{}{}\end{array}$ (57)

This would indeed give an algebraic equation in terms of $X_{0}$. Now any
positive real root of this equation will give the direction of the tangent to
an outgoing null geodesic at the singularity. So the existence of positive
real roots to this equation is a necessary and sufficient condition for the
singularity to be naked. Now as discussed earlier, if one single null geodesic
in the $(t,r)$ plane escapes the singularity, it would mean that a single
wavefront is emitted from the singularity and hence the singularity would
appear to be naked only instantaneously to a distant observer. This would
result in a locally naked singularity. If the singularity is to be seen for a
finite period of time a family of null geodesics must escape from the
singularity thus making it globally naked. These can be investigated from the
number of real positive roots obtained from the above equation. Now we will
consider the results for the two models separately.

###  3.1 Starobinsky’s Model

Using equations (29) and (57), we have

$\displaystyle\frac{2}{X_{0}}=\begin{array}[]{l}lim\\\ \tiny t\rightarrow 0\\\
\tiny r\rightarrow
0\end{array}\frac{\mathcal{G}(E)}{\mathcal{F}(E)}\left[1-\frac{f_{1}(t)}{r}+\frac{f_{2}(t)}{r^{2}}-f_{3}(t)r^{2}-f_{4}(t)r^{3}\right]$
(61)

The functions $f_{i}(t)$, $t=1,2,3,4$ are in fact constants of integration
that we get on solving (integrating) the field equations. So in principle they
should be obtained from suitable initial or boundary conditions and physically
should represent the mass of the system. Obviously these initial and boundary
conditions will correspond to some particular physically viable cases of the
generic mathematical solutions. Therefore we should consider those particular
cases which are compatible with the present study of gravitational collapse.
Our choice for these functions should be guided by the fact that for some
suitable transformations we should be able to retrieve those chosen functions
as initial or boundary conditions for our system. We further mention that we
will only consider self similar expressions for these functions to suit our
analysis. Non self similar expressions can also be considered with
consequences. We choose the following:

$f_{1}(t)=\gamma_{1}t$, $f_{2}(t)=\delta_{1}t^{2}$, $f_{3}(t)=\xi_{1}t^{-2}$,
$f_{4}(t)=\epsilon_{1}t^{-3}$

where $\gamma_{1}$, $\delta_{1}$, $\xi_{1}$ and $\epsilon_{1}$ are arbitrary
constants. It should be clear that the above choices have been made depending
on the definition of $X_{0}$ in equation (57) such that the ratio $t/r$ can be
formed and correspondingly its limit can be evaluated. So we have actually
considered such values of the arbitrary functions of time so that the limit
given by eqn.(61) exists.

Using the above chosen functions in eqn.(61) we get the following algebraic
equation in $X_{0}$

$\delta_{1}X_{0}^{5}-\gamma_{1}X_{0}^{4}+X_{0}^{3}-\frac{2\mathcal{F}(E)}{\mathcal{G}(E)}X_{0}^{2}-\xi_{1}X_{0}-\epsilon_{1}=0$
(62)

The above equation being a five degree equation and hence it is highly
unlikely to find the solution and get the roots by the known mathematical
methods. Hence we try to plot the $X_{0}$ against the other parameters and get
an idea about its trend. Such plots have been generated in Figs.1, 2, 3 and 4.
Also in fig.5 we have shown the effect of gravity’s rainbow on the collapsing
scenario and compared the results with those of general relativity.

Fig.1 Fig.2

Fig.3 Fig.4

Figs.1 and 2 show the variation of $X_{0}$ with $\delta_{1}$ for different
values of $\gamma_{1}$ and $\xi_{1}$ respectively for the Starobinsky’s model.
The Rainbow functions are taken as $\mathcal{F}(E)=1$,
$\mathcal{G}(E)=\sqrt{1-\eta\left(\frac{E_{s}}{E_{P}}\right)^{n}}$. In Fig.1
the initial conditions are taken as $\xi_{1}=2$, $\epsilon_{1}=0.1$, $\eta=1$,
$n=2$, $E_{s}=1$ and $E_{P}=5$. In Fig.2 the initial conditions are
$\gamma_{1}=0.00001$, $\epsilon_{1}=0.1$, $\eta=1$, $n=2$, $E_{s}=1$ and
$E_{P}=5$.

Figs.3 and 4 show the variation of $X_{0}$ with $\xi_{1}$ for different values
of $\delta_{1}$ and $\gamma_{1}$ respectively for the Starobinsky’s model. In
Fig.3 the initial conditions are taken as $\gamma_{1}=0.01$, $\epsilon_{1}=5$,
$\eta=1$, $n=2$, $E_{s}=1$ and $E_{P}=5$. In Fig.4 the initial conditions are
$\delta_{1}=1$, $\epsilon_{1}=0.8$, $\eta=1$, $n=2$, $E_{s}=1$ and $E_{P}=5$.
The rainbow functions are same as that of the previous figures.

Fig.5
Fig.5 shows the effect of the rainbow functions on the outcome of collapse for
the Starobinsky’s model.

###  3.2 Power Law Model

####  3.2.1 Case-1 $(n\geq 4)$

Using equations (39) and (57), we have

$\displaystyle\frac{2}{X_{0}}=\begin{array}[]{l}lim\\\ \tiny t\rightarrow 0\\\
\tiny r\rightarrow
0\end{array}\frac{\mathcal{G}(E)}{\mathcal{F}(E)}\left[1-\frac{f_{5}(t)}{r}+\frac{f_{6}(t)}{r^{2}}\right]$
(66)

Just like the previous section we choose the following:

$f_{5}(t)=\gamma_{2}t$, $f_{6}(t)=\delta_{2}t^{2}$

where $\gamma_{2}$ and $\delta_{2}$ are arbitrary constants. Using the above
functions in eqn.(66) we get the following algebraic equation,

$\delta_{2}X_{0}^{3}-\gamma_{2}X_{0}^{2}+X_{0}-\frac{2\mathcal{F}(E)}{\mathcal{G}(E)}=0$
(67)

Solving the above equation we get only one real root for $X_{0}$. The other
two roots are complex conjugates which are of no interest in this study. The
real root is given by,

$X_{0}=\frac{1}{3\delta_{2}}\left[\gamma_{2}-\frac{2^{1/3}\left(3\delta_{2}-\gamma_{2}^{2}\right)}{P}+\frac{P}{2^{1/3}}\right],~{}~{}~{}~{}~{}~{}\delta_{2}\neq
0,~{}~{}P\neq 0$ (68)

where

$P=\left[2\gamma_{2}^{3}-9\gamma_{2}\delta_{2}+\frac{54\mathcal{F}(E)\delta_{2}^{2}}{\mathcal{G}(E)}+\sqrt{4\left(3\delta_{2}-\gamma_{2}\right)^{3}+\left(2\gamma_{2}^{3}-9\gamma_{2}\delta_{2}+\frac{54\mathcal{F}(E)\delta_{2}^{2}}{\mathcal{G}(E)}\right)^{2}}\right]^{1/3}$
(69)

Here $P$ may be positive or negative satisfying
$4\left(3\delta_{2}-\gamma_{2}\right)^{3}+\left(2\gamma_{2}^{3}-9\gamma_{2}\delta_{2}+\frac{54\mathcal{F}(E)\delta_{2}^{2}}{\mathcal{G}(E)}\right)^{2}\geq
0$.
So the condition for a naked singularity becomes

$\gamma_{2}-\frac{2^{1/3}\left(3\delta_{2}-\gamma_{2}^{2}\right)}{P}+\frac{P}{2^{1/3}}>0,~{}~{}~{}for~{}~{}\delta_{2}>0$
(70)

and that for a black hole is,

$\gamma_{2}-\frac{2^{1/3}\left(3\delta_{2}-\gamma_{2}^{2}\right)}{P}+\frac{P}{2^{1/3}}<0,~{}~{}~{}for~{}~{}\delta_{2}<0$
(71)

In the limit $\mathcal{F}(E)\rightarrow 1,\mathcal{G}(E)\rightarrow 1$ we can
retrieve the corresponding conditions in general relativity.

####  3.2.2 Case-2 $(n=2)$

Using equations (39) and (57), we have

$\displaystyle\frac{2}{X_{0}}=\begin{array}[]{l}lim\\\ \tiny t\rightarrow 0\\\
\tiny r\rightarrow
0\end{array}\frac{\mathcal{G}(E)}{\mathcal{F}(E)}\left[1-A\left(\frac{f_{7}(t)}{r}-\frac{f_{8}(t)}{r^{2}}\right)-B\left(\frac{f_{9}(t)}{r}-\frac{f_{10}(t)}{r^{2}}+f_{11}(t)r^{2}+f_{12}(t)r^{3}\right)\right]$
(75)

We choose the following:

$f_{7}(t)=\gamma_{3}t$, $f_{8}(t)=\delta_{3}t^{2}$, $f_{9}(t)=\gamma_{4}t$,
$f_{10}(t)=\delta_{4}t^{2}$, $f_{11}(t)=\xi_{3}t^{-2}$,
$f_{12}(t)=\epsilon_{3}t^{-3}$

where $\gamma_{3}$, $\delta_{3}$, $\gamma_{4}$, $\delta_{4}$, $\xi_{3}$ and
$\epsilon_{3}$ are arbitrary constants. Using the above functions in eqn.(75)
we get the following algebraic equation,

$\tau X_{0}^{5}-\zeta
X_{0}^{4}+X_{0}^{3}-\frac{2\mathcal{F}(E)}{\mathcal{G}(E)}X_{0}^{2}-B\xi_{3}X_{0}-B\epsilon_{3}=0$
(76)

where $\tau=A\delta_{3}+B\delta_{4}$ and $\zeta=A\gamma_{3}+B\gamma_{4}$.
Since it is highly unlikely to get a feasible solution of the above fifth
degree equation, we generate plots of $X_{0}$ against various parameters to
get an idea about the nature of $X_{0}$. The plots have been generated in
Figs. 6, 7, 8 and 9.

Fig.6 Fig.7

Fig.8 Fig.9

Figs.6 and 7 show the variation of $X_{0}$ with $\tau$ for different values of
$\zeta$ and $\xi_{3}$ respectively for the Power law model (Case-2). In Fig.6
the initial conditions are taken as $\xi_{3}=20$, $\epsilon_{3}=1$, $B=0.05$,
$\eta=1$, $n=2$, $E_{s}=1$ and $E_{P}=5$. In Fig.7 the initial conditions are
$\zeta=10$, $\epsilon_{3}=1$, $B=5$, $\eta=1$, $n=2$, $E_{s}=1$ and $E_{P}=5$.

Figs.8 and 9 show the variation of $X_{0}$ with $\zeta$ for different values
of $\tau$ and $\xi_{3}$ respectively for the Power law model (Case-2). In
Fig.8 the initial conditions are taken as $\xi_{3}=5$, $\epsilon_{3}=0.01$,
$B=5$, $\eta=1$, $n=2$, $E_{s}=1$ and $E_{P}=5$. In Fig.9 the initial
conditions are $\tau=0.1$, $\epsilon_{3}=0.01$, $B=5$, $\eta=1$, $n=2$,
$E_{s}=1$ and $E_{P}=5$.

Fig.10

Fig.10 shows the effect of the rainbow functions on the outcome of collapse
for the power law model (case-2).

###  3.3 Numerical Analysis of the Results

Here we will analyze the numerical results obtained in the previous section.
In case of the Starobinsky’s model we obtained a five degree equation in
$X_{0}$ given by eqn.(62). Unable to get a general solution by the known
algebraic methods, we have resorted to numerical techniques to get an idea of
the nature of $X_{0}$, which in turn determines the nature of singularity
formed. In Figs.1 and 2 we have generated contour plots for $X_{0}$ against
$\delta_{1}$ for different values of $\gamma_{1}$ and $\xi_{1}$ respectively.
From Fig.1 we see that for $\delta_{1}<0$, $X_{0}$ remains in the negative
level for different values of $\gamma_{1}$, thus favouring the formation of a
BH. But for $\delta_{1}>0$, we see trajectories of $X_{0}$ both in the
negative and positive levels. This indicates that there is a realistic chance
of formation of NS for this range. In Fig.2, we have almost a similar result
as in Fig.1. In the both the figures we see that the tendency of formation of
NS increases with the increase in the value of the parameter. So the
dependencies of collapsing procedure on $\gamma_{1}$ and $\xi_{1}$ are almost
of the identical nature.

In Figs.3 and 4 plots are generated for $X_{0}$ against $\xi_{1}$ for
different values of $\delta_{1}$ and $\gamma_{1}$ respectively. In Fig.3 we
see that almost all the trajectories lie in the positive region throughout the
entire domain of $\xi_{1}$, thus favouring the formation of NS. Only around
$\xi_{1}\geq 8$, we see that there is a transition and the possibility of the
formation of an event horizon brightens. It should also be mentioned here that
an increase in the value of $\delta_{1}$ decreases the chance of formation of
NS. A near similar trend is obtained in Fig.4 where we see that the tendency
of formation of NS increases with the increase in $\gamma_{1}$. Moreover the
transition from NS to BH starts at an earlier stage around $\xi_{1}=3$ as
compared to Fig.3. The study would lose its significance if we cannot
understand the role played by gravity’s rainbow in the collapsing scenario.
This is shown in Fig.5, where we clearly see the trajectories for different
rainbow functions and can also compare with that of general relativity. It can
be seen that gravity’s rainbow has a tendency to push the trajectories towards
the positive level thus indicating greater affinity towards the formation of
NS.

In the case-2 $(n=2)$ of the power law model we encounter another five degree
algebraic equation in $X_{0}$ given in eqn.(76). Figs.6, 7, 8 and 9 are
dedicated to the study of this equation. In Figs.6 and 7, we have generated
plots for $X_{0}$ against $\tau$ for different values of $\zeta$ and $\xi_{3}$
respectively. In Fig.6, for $\tau<0$, any singularity that forms will be
clothed by an event horizon, thus forming a BH. On the contrary, for $\tau>0$
we see that there is a fair chance of formation of NS. Eventually what type of
singularity forms depends on the initial conditions for $\tau>0$. Similar
results are visible in Fig.7, where trajectories are obtained for different
values of $\xi_{3}$. In both the plots the tendency of NS increases with the
increase in the values of parameters $\zeta$ and $\xi_{3}$.

In Figs.8 and 9, plots have been obtained for $X_{0}$ against the parameter
$\zeta$ for different values of $\tau$ and $\xi_{3}$ respectively. In Fig.8,
we see that throughout the domain of $\zeta$ trajectories exist both in the
positive and negative level. So tendency of formation of both BH and NS exists
for any value of $\zeta$ and eventually depends on the initial conditions for
the final outcome. But with the increase in the value of $\tau$ the tendency
of formation of NS decreases. From Fig.9 it is evident that for $\xi_{3}<0$,
collapse results in the formation of NS, whereas for $\xi_{3}>0$, the generic
tendency is the formation of BH. In Fig.10, we can see the effect of gravity’s
rainbow on the collapsing process, which is same as the previous case.
Gravity’s rainbow increases the tendency of formation of NS.

In almost all the figures we see that only one trajectory appears at the
positive level, thus indicating that only one null geodesic escapes the
singularity thus making it locally naked in nature. In Fig.4 we find a
striking difference from the above fact. We see that in the $\xi_{1}<0$ region
we have three positive roots of $X_{0}$ for $\gamma_{1}=3$. This would mean
that the singularity formed in this case is visible to a distant observer for
a sufficiently long amount of time to make it globally naked in nature. So the
Starobinsky’s model supports a globally naked singularity. Understanding the
role of initial data, it must be mentioned that any asymptotic nature near the
vertical axis (corresponds to a bundle of null geodesics escaping from the
singularity) would eventually result in global nakedness which may be the case
in Figs.1, 2 and 6.

###  3.4 Strength of the singularity (Curvature growth near the singularity)

The strength of singularity is defined as the measure of its destructive
capacity. The prime concern is that whether extension of space-time is
possible through the singularity or not under any situation. Following Tipler
[78] a curvature singularity is said to be strong if any object hitting it is
crushed to zero volume. In [78] the condition for a strong singularity is
given by,

$\displaystyle\begin{array}[]{c}S=\lim\tau^{2}\psi\\\ \tiny\tau\rightarrow
0\\\ \end{array}\begin{array}[]{c}=\lim\tau^{2}R_{\mu\nu}K^{\mu}K^{\nu}>0\\\
\tiny\tau\rightarrow 0\\\ \end{array}$ (81)

where $R_{\mu\nu}$ is the Ricci tensor, $\psi$ is a scalar given by
$\psi=R_{\mu\nu}K^{\mu}K^{\nu}$, where $K^{\mu}=\frac{dx^{\mu}}{d\tau}$ is the
tangent to the non spacelike geodesics at the singularity and $\tau$ is the
affine parameter. In the paper [79] Mkenyeleye et al. have shown that,

$\displaystyle\begin{array}[]{c}S=\lim\tau^{2}\psi\\\ \tiny\tau\rightarrow
0\\\
\end{array}\begin{array}[]{c}=\frac{1}{4}X_{0}^{2}\left(2\dot{m_{0}}\right)\\\
{}\tiny\hfil\\\ \end{array}$ (86)

where

$\displaystyle\begin{array}[]{c}m_{0}=\lim~{}m(t,r)\\\ \tiny t\rightarrow 0\\\
\tiny r\rightarrow 0\end{array}$ (90)

and

$\displaystyle\begin{array}[]{c}\dot{m_{0}}=\lim\frac{\partial}{\partial~{}t}\left(m(t,r)\right)\\\
\tiny t\rightarrow 0\\\ \tiny r\rightarrow 0\end{array}$ (94)

In ref. [79] it has also been shown that the relation between $X_{0}$ and the
limiting values of mass is given by,

$X_{0}=\frac{2}{1-2m_{0}^{\prime}-2\dot{m_{0}}X_{0}}$ (95)

where

$\displaystyle\begin{array}[]{c}m_{0}^{\prime}=\lim\frac{\partial}{\partial~{}r}\left(m(t,r)\right)\\\
\tiny t\rightarrow 0\\\ \tiny r\rightarrow 0\end{array}$ (99)

and $\dot{m_{0}}$ is given by the eqn.(94). It was shown by Dwivedi and Joshi
in Ref.[80, 53] that a classical singularity in Vaidya spacetime is supposed
to be a strong curvature singularity in a very strong sense. It was also shown
that the conjecture [81] that the strong curvature singularities are never
naked is not true. Moreover the structure of such NS were studied in detail in
Ref.[82] and it was shown that the singularity presents a directional
behaviour in terms of curvature growth along the singular geodesics. In a
quantum regime the singularity formed is supposed to become gravitationally
weak, thus allowing a continuous extension of the spacetime beyond the
singularity [83]. Below we study the strength of the singularities for the
different models.

####  3.4.1 Starobinsky’s Model

Using eqn.(29) in the above relation (86) we get

$\displaystyle\begin{array}[]{c}S=\lim\tau^{2}\psi\\\ \tiny\tau\rightarrow
0\\\
\end{array}\begin{array}[]{c}=\frac{1}{2}X_{0}^{2}\left[\gamma_{1}-2\delta_{1}X_{0}-\frac{2\xi_{1}}{X_{0}^{3}}-\frac{3\epsilon_{1}}{X_{0}^{4}}\right]\\\
{}\tiny\hfil\\\ \end{array}$ (104)

Using eqns.(29), (94) and (99) in eqn.(95) we get an equation for $X_{0}$

$2\delta_{1}X_{0}^{5}-2\gamma_{1}X_{0}^{4}+X_{0}^{3}-2X_{0}^{2}-2\xi_{1}X_{0}-2\epsilon_{1}=0$
(105)

The above equation can be investigated for roots and using these values of
$X_{0}$ in the eqn.(104) we get the conditions for which
$S=\lim\tau^{2}\psi>0$, i.e. the conditions under which we get a strong
singularity. For a strong singularity we have from eqn.(104)

$2\delta_{1}X_{0}^{5}-\gamma_{1}X_{0}^{4}+2\xi_{1}X_{0}+3\epsilon_{1}<0$ (106)

where $X_{0}$ can be obtained as solutions of eqn.(105). Since eqn.(105) is a
fifth degree algebraic equation in $X_{0}$, it is difficult to get a general
solution for $X_{0}$ by the known mathematical methods. However we can get
various solutions numerically. Below we give a particular example of this.

Example

Here we find out the roots of eqn.(105) for a particular set of values of the
parameters involved. Considering $\gamma_{1}=1$, $\delta_{1}=0.1$,
$\xi_{1}=0.5$, $\epsilon_{1}=10$, the roots of the eqn.(105) are found to be
$X_{0}=-1.0636-1.23735~{}i,~{}~{}-1.0636+1.23735~{}i,~{}~{}1.26098-1.52332~{}i,~{}~{}1.26098+1.52332~{}i,~{}~{}9.60523$.
Taking the real value of $X_{0}=9.60523$ and the above values of the
parameters we get $S=-42.7027<0$ from eqn.(104). Hence for this particular set
of values the singularity is weak. Similarly for a different set of initial
conditions we may have different results.

####  3.4.2 Power Law Model

Case-1 $(n\geq 4)$

Using eqn.(39) in the above relation (86) we get,

$\displaystyle\begin{array}[]{c}S=\lim\tau^{2}\psi\\\ \tiny\tau\rightarrow
0\\\
\end{array}\begin{array}[]{c}=\frac{1}{2}X_{0}^{2}\left[\gamma_{2}-2\delta_{2}X_{0}\right]\\\
{}\tiny\hfil\\\ \end{array}$ (111)

Using eqns.(39), (94) and (99) in eqn.(95) we get an equation for $X_{0}$,

$2\delta_{2}X_{0}^{3}-2\gamma_{2}X_{0}^{2}+X_{0}-2=0$ (112)

Solving the above equation we get only one real root given by,

$X_{0}=\frac{\gamma_{2}}{3\delta_{2}}-\frac{6\delta_{2}-4\gamma_{2}^{2}}{6\times
2^{1/3}\delta_{2}Q^{1/3}}+\frac{Q^{1/3}}{3\times 2^{2/3}\delta_{2}}$ (113)

where
$Q=4\gamma_{2}^{3}-9\gamma_{2}\delta_{2}+54\delta_{2}^{2}+3\sqrt{3\left(108\delta_{2}^{4}-36\gamma_{2}\delta_{2}^{3}+2\delta_{2}^{3}+16\gamma_{2}^{3}\delta_{2}^{2}-\gamma_{2}^{2}\delta_{2}^{2}\right)}$.

Using the value of $X_{0}$ from eqn.(113) in eqn.(111) we get the condition
for a strong singularity as,

$\frac{\left(-2\times 2^{2/3}\gamma_{2}^{2}+3\times
2^{2/3}\delta_{2}+\gamma_{2}Q^{1/3}-2^{1/3}Q^{2/3}\right)\left(2\times
2^{2/3}\gamma_{2}^{2}-3\times
2^{2/3}\delta_{2}+2\gamma_{2}Q^{1/3}+2^{1/3}Q^{1/3}\right)^{2}}{216\delta_{2}^{2}Q}>0$
(114)

Case-2 $(n=2)$

Using eqn.(39) in the above relation (86) we get,

$\displaystyle\begin{array}[]{c}S=\lim\tau^{2}\psi\\\ \tiny\tau\rightarrow
0\\\
\end{array}\begin{array}[]{c}=\frac{1}{2}X_{0}^{2}\left[A\left(\gamma_{3}-2\delta_{3}X_{0}\right)+B\left(\gamma_{4}-2\delta_{4}X_{0}-\frac{2\xi_{3}}{X_{0}^{3}}-\frac{3\epsilon_{3}}{X_{0}^{4}}\right)\right]\\\
{}\tiny\hfil\\\ \end{array}$ (119)

Using eqns.(39), (94) and (99) in eqn.(95) we get an equation for $X_{0}$,

$2\left(A\delta_{3}+B\delta_{4}\right)X_{0}^{5}-2\left(A\gamma_{3}+B\gamma_{4}\right)X_{0}^{4}+X_{0}^{3}-2X_{0}^{2}-2B\xi_{3}X_{0}-2B\epsilon_{3}=0$
(120)

Just like the Starobinsky’s model, it is difficult to obtain a general
solution for $X_{0}$ from the above equation. So we will resort to a numerical
solution and give a particular example of it.

Example

Here we consider the following set of numerical values for the parameters
involved,
$\gamma_{3}=1,~{}~{}\gamma_{4}=0.1,~{}~{}\delta_{3}=0.1,~{}~{}\delta_{4}=0.5,~{}~{}\xi_{3}=0.5,~{}~{}\epsilon_{3}=10,~{}~{}A=2,~{}~{}B=5$
Using these values in eqn.(120) we get the roots as,
$X_{0}=-1.26336-1.03348~{}i,~{}~{}-1.26336+1.03348~{}i,~{}~{}0.6967-1.70004~{}i,~{}~{}0.6967+1.70004~{}i,~{}~{}2.05924$.
Taking the real value of $X_{0}=2.05924$ and the above values of the
parameters we get $S=-37.177<0$ from eqn.(119). Hence for this particular set
of values the singularity is weak. But by changing the initial conditions we
may have different results.

Using eqn.(119) we get the condition for a strong singularity for this model
as,

$A\left(\gamma_{3}-2\delta_{3}X_{0}\right)+B\left(\gamma_{4}-2\delta_{4}X_{0}-\frac{2\xi_{3}}{X_{0}^{3}}-\frac{3\epsilon_{3}}{X_{0}^{4}}\right)>0$
(121)

where $X_{0}$ is obtained from eqn.(120).

##  4 A Thermodynamical Analysis of the System

Here we will study the thermodynamical properties of the models.
Thermodynamics is the heart of any physical process which deals with exchange
of heat to and from the system. So this study is crucial in our analysis as it
involves a collapsing mechanism and eventually the formation of a singularity
which can be either a BH or NS. It is a well known fact that BH thermodynamics
is a very important topic of astrophysics. In a fundamental theory of quantum
gravity, it is expected that the thermodynamic properties of a system should
emerge from a microscopic statistical description. Deriving motivation from
these facts we proceed to study the thermodynamical aspects of the BHs formed
as a result of collapse for the various models.

###  4.1 Starobinsky’s Model

The event horizon ($r_{h}$) can be obtained from the relation $f(t,r)=0$,
i.e.,

$-\frac{1}{\mathcal{F}^{2}(E)}\left[1-\frac{1}{r}\left(f_{1}(t)-\frac{f_{2}(t)}{r}+f_{3}(t)r^{3}+f_{4}(t)r^{4}\right)\right]=0$
(122)

The real positive root of the above equation gives the radius of the event
horizon. The thermalization temperature of a system is defined as the
temperature at which the system attains thermal equilibrium. The relation for
thermalization temperature is given by,

$T=\frac{1}{4\pi}\frac{d}{dr}f(t,r)|_{r=r_{h}}$ (123)

For the Starobinsky’s model the expression for thermalization temperature
becomes,

$T=\frac{1}{4\mathcal{F}^{2}(E)\pi
r^{3}}\left[-rf_{1}(t)+2f_{2}(t)+r^{4}\left(2f_{3}(t)+3rf_{4}(t)\right)\right]$
(124)

The entropy of the system is given by,

$S=\pi^{2}r_{h}^{2}$ (125)

where we consider $\pi G=1$. Here $G$ being the Newton’s universal gravitation
constant. Total energy can be obtained from the relation

$U=\int T~{}dS$ (126)

Using the equations (124), (125) and (126) we get the total energy for this
model as,

$U=\frac{\pi}{24\mathcal{F}^{2}(E)r}\left[-12rf_{1}(t)\log(r)-24f_{2}(t)+8r^{4}f_{3}(t)+9r^{5}f_{4}(t)\right]$
(127)

Helmholtz free energy is a thermodynamic potential which is the measure of the
useful work obtainable from a closed system at constant temperature and
volume. It is given by the relation,

$F_{1}=U-TS$ (128)

Using the eqns.(124), (125) and (127) we get the expression Helmholtz free
energy for this model as,

$F_{1}=-\frac{\pi}{24\mathcal{F}^{2}(E)r}\left[6rf_{1}(t)\left\\{2\log(r)-1\right\\}+36f_{2}(t)+r^{4}\left\\{4f_{3}(t)+9rf_{4}(t)\right\\}\right]$
(129)

Specific heat at constant volume is given by,

$C=\left(\frac{dU}{dT}\right)_{V}$ (130)

Using relations (124), (125), (126) and (130) we get the specific heat at
constant volume for the Starobinsky’s model as,

$C=\frac{\pi^{2}r^{2}\left[-rf_{1}(t)+2f_{2}(t)+r^{4}\left\\{2f_{3}(t)+3rf_{4}(t)\right\\}\right]}{rf_{1}(t)-3f_{2}(t)+r^{4}\left\\{f_{3}(t)+3rf_{4}(t)\right\\}}$
(131)

In Figs.11, 12, 13 and 14 the thermodynamical parameters have been plotted
against the radial coordinate $r$ for the Starobinsky’s model. In this section
we have used the rainbow functions given by eqn.(4). In the thermodynamic
parameters it is seen that the rainbow function present is $\mathcal{F}(E)$.
So if we use the rainbow functions given by eqn.(3) as used in the collapse
study, the results would correspond to that of GR (since $\mathcal{F}(E)=1$).
This is the reason why we choose a different set of values for the rainbow
functions for the thermodynamic study. Moreover by using different values of
$\mathcal{F}(E)$ we can understand the effect of gravity’s rainbow on the
thermodynamics of the system and also compare with the results of GR by taking
$\mathcal{F}(E)=1$. From Fig.11 it is seen that the range of thermalization
temperature $T$ increases with the increase in the value of $\delta_{1}$. On
the contrary in Fig.12 it is seen that the range of Internal Energy $U$
decreases with increase in the value of $\delta_{1}$. In Fig.13 the trend for
Helmholtz free energy $F_{1}$ is obtained against $r$. From the figure it is
evident that with the increase in $\delta_{1}$, there is a decrease in the
value of $F_{1}$. In Fig.14 a plot for the specific heat at constant volume
$C$ is obtained against $r$. From the figure we see that the curves are
discontinuous and each of them have two branches. We see that $C$ increases
with the increase in $\delta_{1}$. Finally Fig.15 shows the effect of
gravity’s rainbow on the thermodynamic parameter $T$.

Fig.11 Fig.12

Fig.13 Fig.14

Figs.11, 12, 13 and 14 show the variation of Thermalization temperature $(T)$,
Internal Energy $(U)$, Helmholtz free energy $(F_{1})$ and Specific Heat at
constant volume $(C)$ respectively with $r$ for different values of
$\delta_{1}$ for the Starobinsky’s model. The rainbow function is considered
as $\mathcal{F}(E)=\frac{1}{1-a_{1}(E_{s}/E_{P})}$. In Fig.11 the initial
conditions are taken as $\gamma_{1}=0.1$, $\xi_{1}=1$, $\epsilon_{1}=5$,
$t=5$, $a_{1}=1$, $E_{s}=1$, $E_{P}=5$. In Fig.12 the initial conditions are
$\gamma_{1}=0.1$, $\xi_{1}=5$, $\epsilon_{1}=5$, $t=5$, $a_{1}=1$, $E_{s}=1$,
$E_{P}=5$. In both the figures 13 and 14 the initial conditions are taken as
$\gamma_{1}=0.1$, $\xi_{1}=0.3$, $\epsilon_{1}=3$, $t=5$, $a_{1}=1$,
$E_{s}=1$, $E_{P}=5$.

Fig.15

Fig.15 shows the effect of the rainbow functions on Thermalization temperature
for the Starobinsky’s model.

###  4.2 Power Law Model

####  4.2.1 Case-1 $(n\geq 4)$

The event horizon ($r_{h}$) can be obtained from the relation $f(t,r)=0$,
i.e.,

$-\frac{1}{\mathcal{F}^{2}(E)}\left[1-\frac{1}{r}\left(f_{5}(t)-\frac{f_{6}(t)}{r}\right)\right]=0$
(132)

The real positive root of the above equation gives the radius of the event
horizon. Solving the above equation we get,

$r_{h}=\frac{1}{2}\left[f_{5}(t)\pm\sqrt{\left\\{f_{5}(t)\right\\}^{2}-4f_{6}(t)}\right]$
(133)

Obviously for a realistic event horizon radius we should have:

1\. $\left\\{f_{5}(t)\right\\}^{2}\geq 4f_{6}(t)$

2\. $f_{5}(t)\pm\sqrt{\left\\{f_{5}(t)\right\\}^{2}-4f_{6}(t)}>0\implies
f_{6}>0$.

For this model the expression for thermalization temperature becomes,

$T=\frac{1}{4\mathcal{F}^{2}(E)\pi r^{3}}\left[2f_{6}(t)-rf_{5}(t)\right]$
(134)

Using the equations (134), (125) and (126) we get the total energy for this
model as,

$U=-\frac{\pi}{2\mathcal{F}^{2}(E)}\left[\frac{2f_{6}(t)}{r}+f_{5}(t)\log(r)\right]$
(135)

Using the eqns.(134), (125) and (135) we get the expression Helmholtz free
energy for this model as,

$F_{1}=\frac{\pi}{4\mathcal{F}^{2}(E)r}\left[f_{5}(t)\left\\{r-2r\log(r)\right\\}-6f_{6}(t)\right]$
(136)

Finally using relations (134), (125), (135) and (130) we get the specific heat
at constant volume for this model as,

$C=\frac{\pi^{2}r^{2}\left[-rf_{5}(t)+2f_{6}(t)\right]}{rf_{5}(t)-3f_{6}(t)}$
(137)

Fig.16 Fig.17

Fig.18 Fig.19

Figs.16, 17, 18 and 19 show the variation of Thermalization temperature $T$
and Internal Energy $U$, Helmholtz free energy $F_{1}$ and Specific Heat at
constant volume $C$ respectively with $r$ for different values of $\delta_{2}$
for the Power law model (Case-1). The rainbow function used is
$\mathcal{F}(E)=\frac{1}{1-a_{1}(E_{s}/E_{P})}$. In both the figures 16 and 17
the initial conditions are taken as $\gamma_{2}=0.1$, $t=5$, $a_{1}=1$,
$E_{s}=1$, $E_{P}=5$. In Fig.18 the initial conditions are $\gamma_{2}=0.1$,
$t=5$, $a_{1}=1$, $E_{s}=1$, $E_{P}=5$. In Fig.19 the initial conditions are
$\gamma_{2}=0.5$, $t=5$, $a_{1}=1$, $E_{s}=1$, $E_{P}=5$.

Fig.20

Fig.20 shows the effect of the rainbow functions on Thermalization temperature
for the Power law model (Case-1).

In Figs.16, 17, 18 and 19, plots for the thermodynamical parameters have been
generated against $r$ for the power law model (Case-1 $(n\geq 4)$). In Fig.16
we see that the thermalization temperature $(T)$ decreases with the increase
of $r$. Moreover with the increase of $\delta_{2}$, the range of $T$
increases. From Fig.17 we see that the internal energy $(U)$ increases with
the increase of $r$. On the contrary, with the increase of $\delta_{2}$, there
is decrease in the range of $U$. From Fig.18 it is evident that Helmholtz free
energy $(F_{1})$ increases with $r$. Moreover with an increase in the value of
$\delta_{2}$, there is a considerable decrease in the value of $F_{1}$. In
Fig.19 the trend of specific heat at constant volume $(C)$ is exhibited
against $r$. It is seen that with an increase in the value of $\delta_{2}$
there is a decrease in $C$. Finally in Fig.20, we see the effect of the
rainbow function in the thermodynamics of the system.

####  4.2.2 Case-2 $(n=2)$

The event horizon ($r_{h}$) can be obtained from the relation $f(t,r)=0$,
i.e.,

$-\frac{1}{\mathcal{F}^{2}(E)}\left[1-\frac{1}{r}\left\\{A\left(f_{7}(t)-\frac{f_{8}(t)}{r}\right)+B\left(f_{9}(t)-\frac{f_{10}(t)}{r}+f_{11}(t)r^{3}+f_{12}(t)r^{4}\right)\right\\}\right]=0$
(138)

The real positive root of the above equation gives the radius of the event
horizon. Here, the expression for thermalization temperature becomes,

$T=\frac{1}{4\mathcal{F}^{2}(E)\pi
r^{3}}\left[-Arf_{7}(t)+2Af_{8}(t)-Brf_{9}(t)+2Bf_{10}(t)+2Br^{4}f_{11}(t)+3Br^{5}f_{12}(t)\right]$
(139)

Using the equations (139), (125) and (126) we get the total energy for this
model as,

$U=\frac{\pi}{2\mathcal{F}^{2}(E)}\left[-\log(r)\left\\{Af_{7}(t)+Bf_{9}(t)\right\\}-\frac{2}{r}\left\\{Af_{8}(t)+Bf_{10}(t)\right\\}+\frac{2}{3}Br^{3}f_{11}(t)+\frac{3}{4}Br^{4}f_{12}(t)\right]$
(140)

Using the eqns.(139), (125) and (140) we get the expression Helmholtz free
energy for this model as,

$\displaystyle
F_{1}=-\frac{\pi}{24\mathcal{F}^{2}(E)r}\left[12r\log(r)\left\\{Af_{7}(t)+Bf_{9}(t)\right\\}-6Arf_{7}(t)+36Af_{8}(t)-6Brf_{9}(t)\right.$
$\left.+36Bf_{10}(t)+4Br^{4}f_{11}(t)+9Br^{5}f_{12}(t)\right]$ (141)

Finally using relations (139), (125), (140) and (130) we get the specific heat
at constant volume for this model as,

$C=\frac{\pi^{2}r^{2}\left[-Arf_{7}(t)+2Af_{8}(t)-Brf_{9}(t)+2Bf_{10}(t)+2Br^{4}f_{11}(t)+3Br^{5}f_{12}(t)\right]}{Arf_{7}(t)-3Af_{8}(t)+Brf_{9}(t)-3Bf_{10}(t)+Br^{4}f_{11}(t)+3Br^{5}f_{12}(t)}$
(142)

In the Figs.21, 22, 23 and 24 the thermodynamical parameters for the power law
model (Case-2 $(n=2)$) are plotted against $r$. From Fig.21 it is evident that
the thermalization temperature $(T)$ decreases with the increase in $r$. With
the increase in $\delta_{3}$, there is an increase in the value of $T$. From
Fig.22 it is clear that there is an increase in the internal energy $(U)$ of
the system with an increase in size, i.e. $r$. Moreover as $\delta_{3}$
increases there is a corresponding decrease in the value of $U$. Fig.23 shows
the variation of Helmholtz free energy $(F_{1})$ against $r$. From the figure
it is seen that $F_{1}$ increases with an increase in $r$. Further it is
evident that with an increase in $\delta_{3}$ there is a decrease in $F_{1}$.
Finally in Fig.24 we see that the specific heat at constant volume $(C)$
decreases with an increase in $r$. With an increase of $\delta_{3}$, there is
a corresponding increase in $C$. Finally in Fig.25 we see the effect of
gravity’s rainbow on the thermalization temperature.

Fig.21 Fig.22

Fig.23 Fig.24

Figs.21, 22, 23 and 24 show the variation of Thermalization temperature $T$
and Internal Energy $U$, Helmholtz free energy $F_{1}$ and Specific Heat at
constant volume $C$ respectively with $r$ for different values of $\delta_{3}$
for the Power law model (Case-2). The rainbow function is taken as
$\mathcal{F}(E)=\frac{1}{1-a_{1}(E_{s}/E_{P})}$. In both the figures 21 and 22
the initial conditions are taken as $\gamma_{3}=0.1$, $\gamma_{4}=0.5$,
$\xi_{3}=0.7$, $\epsilon_{3}=2$, $\delta_{4}=10$, $A=0.5$, $B=0.8$, $t=5$,
$a_{1}=1$, $E_{s}=1$, $E_{P}=5$. In Fig.23 the initial conditions are
$\gamma_{3}=0.1$, $\gamma_{4}=0.5$, $\xi_{3}=0.7$, $\epsilon_{3}=2$,
$\delta_{4}=10$, $A=0.5$, $B=0.8$, $t=5$, $a_{1}=1$, $E_{s}=1$, $E_{P}=5$. In
Fig.24 the initial conditions are $\gamma_{3}=0.1$, $\gamma_{4}=0.5$,
$\xi_{3}=0.7$, $\epsilon_{3}=2$, $\delta_{4}=10$, $A=5$, $B=0.8$, $t=5$,
$a_{1}=1$, $E_{s}=1$, $E_{P}=5$.

Fig.25

Fig.25 shows the effect of the rainbow functions on Thermalization temperature
for the Power law model (Case-2).

##  5 Conclusion and Discussion

In this work we have investigated rainbow modified $f(R)$ gravity via a
gravitational collapse and thermodynamical study. We have considered a non-
static (time dependent) geometry, namely the Vaidya spacetime for our study.
The collapsing procedure being a time dependent phenomenon justifies our
choice. Rainbow modifications to Vaidya geometry was considered, thus making
the spacetime energy dependent in nature. As discussed earlier the rainbow
modifications to a metric is significant only in the quantum regime, i.e. when
the energy of the particle probed is comparable to the Planck energy
$(E_{P}\approx 10^{19}GeV)$ and the length scale tends to the Planck length
$(L_{P}\approx 10^{-33}cm)$. Under such a scenario the quantum fluctuations
will be at play, and will have far reaching effect on the dynamics of the
system. This was the motivation of our model and we wanted to probe such
effects via a theoretical framework.

The field equations for such a system was derived for the general $f(R)$
expression. In order to get a solution of the system, we needed to consider
various models of $f(R)$ gravity. In this work we considered the famous
inflationary Starobinsky’s model and the Power law model of $f(R)$ gravity
with suitable motivations. Solutions for each of these models were obtained
separately which gave the mass of the collapsing system as a function of both
the radial and temporal coordinates. Using them in the metric we obtained the
geometries for each of the models. Here it should be mentioned that in the
process of solving the system we have used the $(11)$ component of the field
equations given by eqn.(25). This has been done purely due to the simplicity
of the equation compared to the others with no other bias. This choice helped
us to reduce our computational efforts considerably, but in doing so we have
missed out on a very important parameter, namely the density parameter $\rho$
of the matter content of the system. This is solely because $\rho$ did not
feature in the $(11)$ component of the field equations and hence remained
absent in the subsequent analysis. So in our analysis dependence of the
collapsing procedure on the matter content of the system could not be studied,
which in turn restricted us from mapping our results to various cosmological
eras.

A geodesic study was performed in order to study the collapsing procedure of
the star in both the models. The study was aimed at enquiring the nature of
singularity resulting from the collapse. A singularity censored by an event
horizon will be a BH, but any singularity devoid of such horizon will be a
naked one and easily accessible from the outside. So the loss of information
and other problems posed by a BH will cease to exist for such type of
singularities. This motivated us to initiate such a study. However the
complexity of the systems provoked us to use numerical analysis in our study.
In this context we mention here that our choice for the functions $f_{i}(t)$,
$i=1,2,3...12$ seems to be self similar in nature. This has been done in
accordance with the definition of $X_{0}$ in eqn.(57), such that the gradient
$t/r$ can be formed without much complications. However, naturally the
question arises that what would happen if we consider non self-similar terms?
From the mathematical framework we see that under the limit $r\rightarrow
0,t\rightarrow 0$ non self-similar terms will result in the deletion of
crucial terms or in the creation of mathematically undefined terms, both of
which are undesirable. Deletion of important terms will result in loss of
information about the system while undefined terms will make further
mathematical computations next to impossible. From the study it is seen that
NS is quite a possibility along with BH for the collapsing models. Fine tuning
the initial conditions we may get NS or BH for different ranges of the
parameters involved. It should be mentioned that this is a significant
counter-example of the cosmic censorship hypothesis. Moreover it is seen that
the naked singularity formed can be globally naked for certain scenarios.

An important result that we have derived from this study is that under the
effect of gravity’s rainbow there is greater tendency of formation of NS. The
deviations produced by the presence of rainbow functions favours the formation
of NS. This is a very interesting finding given the renewed interest of
scientists in NS. It is admitted that the result obtained depends on the
initial data, which is expected from previous works like [54, 84]. Over the
years it has been seen that the role of initial data is central to any study
of gravitational collapse. We know that we have not been able to detect a NS
till date. In the light of such a comment it seems that our work loses some
significance. But on the contrary we would like to emphasize on the fact that
our work should generate widespread motivations for the search of NS. This is
more so because for quite sometime it is strongly believed that with the
failure of cosmic censorship we would come face to face with the laws of
quantum gravity, whenever gravitational collapse of a distant star results in
a NS. Thus a complete Theory Of Everything (TOE) needs a proper knowledge of
NS, which can act as a laboratory for quantum gravity.

We have also investigated the conditions for a strong singularity for
different models. The conditions have been derived and supplemented with
suitable examples where required. Finally we have concluded our analysis with
a thermodynamical study of the system. Expressions for various thermodynamical
parameters like the thermalization temperature, Internal energy, etc. have
been deduced for the models and their nature have been analyzed via various
plots. In a classical set-up it is expected that the thermalization
temperature $(T)$ should decrease with distance from the core of the star or
with an increase in the size of the star. From the definition of the Internal
energy $(U)$ it is clear that its value should increase with the size of the
star for a given $T$. Similarly the Helmholtz free energy $F_{1}$ should also
increase with an increase in $r$, given its definition and relative trends of
$T$ and $U$. Our results show deviations from the above expected results. Such
deviations can be attributed to the quantum fluctuations of our model. Our
study reveals information that characterizes a quantum evolution of the
universe. So we are hopeful on the fact that these results will be useful in
any future attempts towards the formulation of a successful theory of quantum
gravity.

## Acknowledgments

The author acknowledges the Inter University Centre for Astronomy and
Astrophysics (IUCAA), Pune, India for granting visiting associateship.

##  6 Appendix

Here we report the other components of the field equations for this model
which have not been used in our analysis.

1\. The (00)-component of field equations is given by,

$\displaystyle
f^{\prime\prime}(R)\mathcal{G}^{2}(E)\left\\{-2\mathcal{F}^{2}(E)r^{4}\left(2\ddot{m}^{\prime}+\ddot{m}^{\prime\prime}r\right)+\mathcal{G}^{2}(E)\left(m-r\right)\left(m-m^{\prime}r\right)\left(-4m^{\prime}+r\left(m^{\prime\prime}+m^{\prime\prime\prime}r\right)\right)\right.$
$\displaystyle\left.-\mathcal{F}(E)\mathcal{G}(E)r^{2}\left(m\dot{r}\left(m^{\prime\prime}+m^{\prime\prime\prime}r\right)-m\left(2\dot{m}^{\prime}+\dot{m}^{\prime\prime}r\right)+m^{\prime}\left(-4\dot{m}+r\left(2\dot{m}^{\prime}+\dot{m}^{\prime\prime}r\right)\right)\right)\right\\}$
$\displaystyle-r^{2}\left\\{-2\mathcal{F}(E)\mathcal{G}(E)f^{\prime}(R)\dot{m}r^{2}+2\mathcal{F}^{2}(E)\mathcal{G}^{4}(E)f^{\prime\prime\prime}(R)\left(2\dot{m}^{\prime}+\dot{m}^{\prime\prime}r\right)^{2}+r\left(r-m\right)\left(f^{\prime}(R)\left(-2\mathcal{G}^{2}(E)m^{\prime}+r^{2}R\right)\right.\right.$
$\left.\left.+r^{2}\left(2~{}\Box
f^{\prime}(R)-f(R)+2\kappa\rho\right)\right)+2r^{4}\kappa\sigma\right\\}=0$
(143)

where $\Box f^{\prime}(R)$ is given by,

$\displaystyle\Box
f^{\prime}(R)=\frac{1}{r^{7}}\left[f^{\prime\prime\prime}(R)\mathcal{G}^{5}(E)\left\\{4m^{\prime}-r\left(m^{\prime\prime}+rm^{\prime\prime\prime}\right)\right\\}\left\\{-2\mathcal{F}(E)r^{2}\left(2\dot{m}^{\prime}+r\dot{m}^{\prime\prime}\right)+\mathcal{G}(E)\left(m-r\right)\right.\right.$
$\displaystyle\left.\left.\left(-4m^{\prime}+r\left(m^{\prime\prime}+rm^{\prime\prime\prime}\right)\right)\right\\}+f^{\prime\prime}(R)\mathcal{G}^{3}(E)r^{2}\left\\{2r^{2}\mathcal{F}(E)\left(-2\dot{m}^{\prime}+r\left(2\dot{m}^{\prime\prime}+\dot{m}^{\prime\prime\prime}r\right)\right)\right.\right.$
$\displaystyle\left.\left.+\mathcal{G}(E)\left(-m\left(8m^{\prime}+r\left(-5m^{\prime\prime}+r^{2}m^{\prime\prime\prime\prime}\right)\right)+r\left(4m^{\prime
2}-m^{\prime}\left(-4+r\left(m^{\prime\prime}+rm^{\prime\prime\prime}\right)\right)\right.\right.\right.\right.$
$\left.\left.\left.\left.+r\left(-4m^{\prime\prime}+r\left(m^{\prime\prime\prime}+rm^{\prime\prime\prime\prime}\right)+rm^{\prime\prime\prime}\right)\right)\right)\right\\}\right]$
(144)

2\. The (22)-component of field equations is given by,

$\displaystyle-2f^{\prime\prime\prime}(R)\mathcal{G}^{5}(E)\left\\{4m^{\prime}-r\left(m^{\prime\prime}+m^{\prime\prime\prime}r\right)\right\\}\left[-2\mathcal{F}(E)r^{2}\left(2\dot{m}^{\prime}+\dot{m}^{\prime\prime}r\right)+\mathcal{G}(E)\left(m-r\right)\left\\{-4m^{\prime}+r\left(m^{\prime\prime}+m^{\prime\prime\prime}r\right)\right\\}\right]$
$\displaystyle+r^{2}\left[2f^{\prime\prime}(R)\mathcal{G}^{3}(E)\left\\{\mathcal{F}(E)r^{2}\left(6\dot{m}^{\prime}-r\left(3\dot{m}^{\prime\prime}+2\dot{m}^{\prime\prime\prime}r\right)\right)+\mathcal{G}\left(r\left(-4m^{\prime
2}+m^{\prime}\left(-8+r\left(m^{\prime\prime}+m^{\prime\prime\prime}r\right)\right)-r\left(-5m^{\prime\prime}\right.\right.\right.\right.\right.$
$\left.\left.\left.\left.\left.+m^{iv}r^{2}+rm^{\prime\prime\prime}\right)\right)+m\left(12m^{\prime}+r\left(-6m^{\prime\prime}-m^{\prime\prime\prime}r+m^{iv}r^{2}+rm^{\prime\prime\prime}\right)\right)\right)\right\\}+r^{5}\left(f(R)-f^{\prime}(R)R+2\kappa\omega\rho\right)\right]=0$
(145)

3\. The (33)-component of field equations is given by,

$\displaystyle-2f^{\prime\prime\prime}(R)\mathcal{G}^{5}(E)\left\\{4m^{\prime}-r\left(m^{\prime\prime}+m^{\prime\prime\prime}r\right)\right\\}\left[-2\mathcal{F}(E)r^{2}\left(2\dot{m}^{\prime}+\dot{m}^{\prime\prime}r\right)+\mathcal{G}\left(m-r\right)\left\\{-4m^{\prime}+r\left(m^{\prime\prime}+m^{\prime\prime\prime}r\right)\right\\}\right]$
$\displaystyle+r^{2}\left[2f^{\prime\prime}(R)\mathcal{G}^{3}(E)\left\\{12\mathcal{G}(E)mm^{\prime}-2\mathcal{G}(E)\left(3mm^{\prime\prime}+2m^{\prime}\left(2+m^{\prime}\right)\right)r+\left(\mathcal{G}m^{\prime\prime}\left(m^{\prime}+5\right)+6\mathcal{F}(E)\dot{m}^{\prime}\right)r^{2}+\left(\mathcal{G}(E)\right.\right.\right.$
$\left.\left.\left.\left(mm^{iv}+m^{\prime\prime\prime}\left(m^{\prime}-1\right)\right)-3\mathcal{F}(E)\dot{m}^{\prime\prime}\right)r^{3}-\left(\mathcal{G}(E)m^{iv}+2\mathcal{F}(E)\dot{m}^{\prime\prime\prime}\right)r^{4}\right\\}+r^{5}\left(f(R)-f^{\prime}(R)R+2\kappa\omega\rho\right)\right]=0$
(146)

4\. The $(01)$ or $(10)$-component of field equations is given by,

$\displaystyle
2f^{\prime\prime\prime}(R)\mathcal{G}^{5}(E)\left\\{4m^{\prime}-r\left(m^{\prime\prime}+m^{\prime\prime\prime}r\right)\right\\}\left[\mathcal{F}(E)r^{2}\left(2\dot{m}^{\prime}+\dot{m}^{\prime\prime}r\right)-\mathcal{G}(E)\left(m-r\right)\left\\{-4m^{\prime}+r\left(m^{\prime\prime}+m^{\prime\prime\prime}r\right)\right\\}\right]$
$\displaystyle+r^{2}\left[f^{\prime\prime}(R)\mathcal{G}^{3}(E)\left\\{-2\mathcal{F}(E)r^{3}\left(3\dot{m}^{\prime\prime}+\dot{m}^{\prime\prime\prime}r\right)+\mathcal{G}(E)\left(r\left(-4m^{\prime
2}+m^{\prime}\left(-8+r\left(m^{\prime\prime}+m^{\prime\prime\prime}r\right)\right)-2r\left(-4m^{\prime\prime}\right.\right.\right.\right.\right.$
$\displaystyle\left.\left.\left.\left.\left.+r\left(m^{\prime\prime\prime}+m^{iv}r\right)+rm^{\prime\prime\prime}\right)\right)+m\left(12m^{\prime}+r\left(-9m^{\prime\prime}+r\left(m^{\prime\prime\prime}+2m^{iv}r\right)+2rm^{\prime\prime\prime}\right)\right)\right)\right\\}+r^{5}\left(f(R)\right.\right.$
$\left.\left.-f^{\prime}(R)R-2\kappa\rho\right)\right]=0$ (147)

## References

  * [1] S. Perlmutter et. al. :- Astrophys. J. 517 565 (1999).
  * [2] A. G. Riess et al. :- Astron. J. 116 1009 (1998).
  * [3] S. Nojiri, S. D. Odintsov, V. K. Oikonomou :- Phys. Rep. 692 1 (2017).
  * [4] S. Nojiri and S. D. Odintsov :- Int. J. Geom. Methods Mod. Phys. 04 115 (2007).
  * [5] S. Capozziello, R. D Agostino, O. Luongo:- Int. J. Mod. Phys. D 28 1930016 (2019).
  * [6] T. P. Sotiriou, V. Faraoni :- Rev. Mod. Phys. 82 451 (2010).
  * [7] A. De Felice, S. Tsujikawa :- Living Rev. Relativity 13 3 (2010).
  * [8] L. Amendola, D. Polarski, S. Tsujikawa :- Phys. Rev. Lett. 98 131302 (2007).
  * [9] T. P. Sotiriou :- Classical Quantum Gravity 23 5117 (2006).
  * [10] Y.-S. Song,W. Hu, I. Sawicki :- Phys. Rev. D 75 044004 (2007).
  * [11] S. Nojiri, S. D. Odintsov :- Phys. Rev. D 74 086005 (2006).
  * [12] P. Shah, G. C. Samanta :- Eur. Phys. J. C. 79 414 (2019).
  * [13] S. D. Odintsov, V. K. Oikonomou :- Phys. Rev. D 96 104049 (2017).
  * [14] K. Bamba, S. Nojiri, S. D. Odintsov, D. S ez-G mez :- Phys. Rev. D 90 124061 (2014).
  * [15] S. Capozziello, M. De Laurentis, V. Faraoni :- Open Astron. J. 3 49 (2010).
  * [16] C. Rovelli :- Class. Quantum Grav. 8 1613 (1991).
  * [17] C. Rovelli :- Living Rev.Rel. 1 1 (1998).
  * [18] M. Kaku, K. Kikkawa :- Phys.Rev. D 10 1110 (1974).
  * [19] M. Kaku, K. Kikkawa :- Phys.Rev. D 10 1823 (1974).
  * [20] P. Horava :- Phys. Rev. Lett. 102 161301 (2009).
  * [21] P. Horava :- Phys. Rev. D 79 084008 (2009).
  * [22] J. Magueijo, L. Smolin :- Class. Quant. Grav. 21 1725 (2004).
  * [23] R. Garattini, E.N. Saridakis :- Eur. Phys. J. C 75 343 (2015).
  * [24] R. Iengo, J.G. Russo, M. Serone :- JHEP 0911 020 (2009).
  * [25] B. M. Gripaios :- JHEP 0410 069 (2004).
  * [26] J. Alfaro, P. Gonzalez, R. Avila :- Phys. Rev. D 91 105007 (2015).
  * [27] H. Belich, K. Bakke :- Phys. Rev. D 90 025026 (2014).
  * [28] V.A. Kostelecky, S. Samuel :- Phys. Rev. D 39 683 (1989).
  * [29] G. Hooft :- Class. Quant. Grav. 13 1023 (1996).
  * [30] G. Amelino-Camelia, J.R. Ellis, N.Mavromatos, D.V.Nanopoulos, S. Sarkar :- Nature 393 763 (1998).
  * [31] S.M. Carroll, J.A. Harvey, V.A. Kostelecky, C.D. Lane, T. Okamoto :- Phys. Rev. Lett. 87 141601 (2001).
  * [32] R. Gambini, J. Pullin :- Phys. Rev. D 59 124021 (1999).
  * [33] J. Magueijo, L. Smolin :- Phys. Rev. Lett. 88 19 (2002).
  * [34] G. Amelino-Camelia, J.R. Ellis, N.E. Mavromatos, D.V. Nanopoulos :- Int. J. Mod. Phys. A 12 607 (1997).
  * [35] G. Amelino-Camelia :- Living Rev. Relativ. 16 5 (2013).
  * [36] M. Dehghani :- Phys. Lett. B 777 351 (2018).
  * [37] S. H. Hendi, A.Dehghani, M. Faizal :- Nucl. Phys. B 914 117 (2017).
  * [38] S. H. Hendi, S. Panahiyan, B.E. Panah, M. Faizal, M. Momennia :- Phys. Rev. D 94 024028 (2016).
  * [39] S. H. Hendi, M. Momennia, B.E. Panah, M. Faizal :- Astrophys. J. 827 153 (2016).
  * [40] A. Ashour, M. Faizal, A.F. Ali, F. Hammad :- Eur. Phys. J. C. 76 264 (2016).
  * [41] S.H. Hendi, M. Faizal, B.E. Panah, S. Panahiyan :- Eur. Phys. J. C. 76 296 (2016).
  * [42] S.H. Hendi, M. Faizal :- Phys. Rev. D 92 044027 (2015).
  * [43] A.F. Ali, M. Faizal, B. Majumder :- Europhys. Lett. 109 20001 (2015).
  * [44] J. R. Oppenhiemer, H. Snyder :- Phys. Rev. 56 455 (1939).
  * [45] R. C. Tolman :- Proc. Natl. Acad. Sci. USA 20 169 (1934);
  * [46] H. Bondi :- Mon. Not. Astron. Soc. 107 410 (1947).
  * [47] R. Penrose :- Riv. Nuovo Cimento 1 252 (1969).
  * [48] D. M. Eardley, L. Smarr :- Phys. Rev. D 19 2239 (1979).
  * [49] D. Christodoulou :- Commun. Math. Phys. 93 171 (1984).
  * [50] R. P. A. C. Newman :- Class. Quantum Grav. 3 527 (1986).
  * [51] I. H. Dwivedi, P. S. Joshi :- Class. Quantum Grav. 9 L39 (1992).
  * [52] P. S. Joshi, I. H. Dwivedi :- Phys. Rev. D 47 5357 (1993).
  * [53] P. S. Joshi, I.H. Dwivedi :- Commun. Math. Phys. 146 333 (1992).
  * [54] P. S. Joshi, T.P. Singh :- it Phys. Rev. D 51 6778 (1995).
  * [55] P. C. Vaidya :- Proc. Indian Acad. Sci. Sect. A 33 264 (1951).
  * [56] P. Rudra, M. Faizal, A. F. Ali :- Nucl. Phys. B 909 725 (2016).
  * [57] Y. Heydarzade, P. Rudra, F. Darabi, A. F. Ali, M. Faizal :- Phys. Lett. B 774 46 (2017).
  * [58] P. Rudra, S. Maity :- Eur. Phys. J. C 78 828 (2018).
  * [59] Y. Heydarzade, P. Rudra, B. Pourhassan, M. Faizal, A. F. Ali :- J. Cosmol. Astropart. Phys. 06 038 (2018).
  * [60] P. Rudra, U. Debnath :- Can. J. Phys. 92(11) 1474 (2014).
  * [61] P. Rudra, R. Biswas, U. Debnath :- Astrophys. Space Sci. 354 2 (2014).
  * [62] U. Debnath, P. Rudra, R. Biswas :- Astrophys. Space Sci. 339 135 (2012).
  * [63] P. Rudra, R. Biswas, U. Debnath :- Astrophys.Space Sci. 335 505 (2011).
  * [64] A. Perez, H. Sahlmann, D. Sudarsky :- Classical Quantum Gravity 23 2317 (2006).
  * [65] D. Sudarsky :- Int. J. Mod. Phys. D 20 509 (2011).
  * [66] S. Landau, G. Le n, D. Sudarsky :- Phys. Rev. D 88 023526 (2013).
  * [67] A. de Un nue, D. Sudarsky :- Phys. Rev. D 78 043510 (2008).
  * [68] G. Le n, D. Sudarsky :- Classical Quantum Gravity 27 225017 (2010).
  * [69] G. Le n, A. De Un nue, D. Sudarsky :- Classical Quantum Gravity 28 155010 (2011).
  * [70] A. Diez-Tejedor, D. Sudarsky :- J. Cosmol. Astropart. Phys. 07 (2012) 045.
  * [71] S. J. Landau, C. G. Sc ccola, D. Sudarsky :- Phys. Rev. D 85 123001 (2012).
  * [72] P. Ca ate, P. Pearle, D. Sudarsky :- Phys. Rev. D 87 104024 (2013).
  * [73] A. A. Starobinsky :- Phys. Lett. B 91 99 (1980).
  * [74] G. Panotopoulos, I. Lopes :- Phys. Rev. D 97 no.2 024025 (2018).
  * [75] A. De Felice, S. Tsujikawa :- Living Rev. Relativ.13 3 (2010).
  * [76] T. P. Sotiriou, V. Faraoni :- Rev. Mod. Phys. 82 (2010).
  * [77] T. P. Singh, P.S. Joshi :- Class. Quant. Gravity 13 559 (1996).
  * [78] F. J. Tipler :- Phys. Lett. A. 64, 8 (1977)
  * [79] M. D. Mkenyeleye, R. Goswami, and S. D. Maharaj :- Phys. Rev. D 90 064034 (2014).
  * [80] I. H. Dwivedi, P. S. Joshi :- Class. Quantum Grav. 6 1599 (1989).
  * [81] F. J. Tipler, C. J. S. Clarke, G. F. R. Ellis :- General Relativity and Gravitation vol 2 p97 ed A Held (NewYork Plenum) (1980).
  * [82] I. H. Dwivedi, P. S. Joshi :- Class. Quantum Grav. 8 1339 (1991).
  * [83] A. Bonanno, B. Koch, A. Platania :- Foundations of Phys. 48 1393 (2018).
  * [84] S. Jhingan, P. S. Joshi, T. P. Singh :- it Class. Quantum Grav.13 3057 (1996).