A non-static quantum inspired spacetime in f(R) gravity: Gravity's rainbow

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.


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, L GR = R is replaced by an analytic function of R i.e. 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 F(E) and G(E) as given below, 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 energymomentum dispersion relations of classical GR from the Eqn. (1). This means that the rainbow functions are required to satisfy the relations, In the above relations the limit E S /E P → 0 corresponds to a spacetime with low energy where GR becomes dominant. It is expected that as E s → 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 κ-Minkowski noncommutative spacetime these functions have been proposed as [34,35], where η is a constant. The modified dispersion from constant velocity of light motivates the following rainbow functions [33] F(E s /E P ) = G(E s /E P ) = 1 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], 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 where e 0 = F −1 (E s /E P )ẽ 0 and e i = G −1 (E s /E P )ẽ 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.

Rainbow deformed Vaidya spacetime in f (R) gravity
The Einstein-Hilbert action of GR is given by, where κ ≡ 8π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, Beginning from the action (8) and adding a matter term S M , the total action for f (R) gravity takes the form, where ψ collectively denotes the matter fields. Taking variation with respect to the metric we get the field equations as, where T μν is given by, In the above equations a prime denotes differentiation with respect to the argument, i.e. R. ∇ μ denotes covariant derivative associated with the Levi-Civita connection of the metric and ≡ ∇ μ ∇ μ is the D'Alembertian operator. The field equations given by equation (10) can also be written in the following form, The Vaidya metric in the advanced time coordinate system is given by, where f (t, r) = − 1 − m(t,r) r . It is obvious that here we have used the units G = c = 1. Introducing rainbow deformations in the above metric we get [22], The total energy momentum tensor of the field equation (12) is given by the following sum, where T (n) μν and T (m) μν are the contributions from the Vaidya null radiation and perfect fluid respectively defined as, and T (m) μν = (ρ + p)(l μ η ν + l ν η μ ) + pg μν (17) where ρ and p are the energy density and pressure for the perfect fluid and σ is the energy density corresponding to Vaidya null radiation. In the co-moving co-ordinates (t, r, θ 1 , θ 2 , ..., θ n ), the two eigen vectors of energy-momentum tensor namely l μ and η μ are linearly independent future pointing null vectors having components and they satisfy the relations Now we impose rainbow deformations on the linearly independent future pointing null vectors l μ and η μ such that we get, satisfying the following conditions (19). Therefore, the non-vanishing components of the total energy-momentum tensor will be as follows Here we consider matter in the form of perfect barotropic fluid given by the equation of state p = ωρ (22) where 'ω' is the barotropic parameter. The non-vanishing components of the Einstein tensors are given by, where dot and dash represents derivatives with respect to time t and radial coordinate r respectively. For this system the Ricci scalar becomes,

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,
where m iv denotes the fourth partial derivative of m with respect to r.

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.

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] where a ≥ 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', We can see that for a = 0, we get an identity which corresponds to the case for GR. From the above equation we get, It should be noted that the parameter 'a' does not appear in the above differential equation. Solving the above differential equation we get, 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,

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, where α and n > 0 are constants. Using eqns. (24), (25) and (31) we get the following differential equation in terms of the mass 'm', 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: For eqn. (33) the solution is obtained as, 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, 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), 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, Hence for the power law model, the rainbow deformed Vaidya spacetime in f (R) gravity is given by,

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 nonspacelike geodesics possess well defined tangent at the singularity, 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 The equation for outgoing radial null geodesics can be obtained from equation (14) by putting 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 = 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 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.

Starobinsky's model
Using equations (29) and (42), we have 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: where γ 1 , δ 1 , ξ 1 and 1 are arbitrary constants. It should be clear that the above choices have been made depending on the definition of X 0 in equation (42) 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. (43) exists.
Using the above chosen functions in eqn. (43) we get the following algebraic equation in X 0 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.

Power law model
Using equations (39) and (42), we have Just like the previous section we choose the following: where γ 2 and δ 2 are arbitrary constants. Using the above functions in eqn. (45) we get the following algebraic equation, 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, where Here P may be positive or negative satisfying 4 (3δ So the condition for a naked singularity becomes and that for a black hole is, In the limit F(E) → 1, G(E) → 1 we can retrieve the corresponding conditions in general relativity.

Case-2 (n = 2)
Using equations (39) and (42), we have We choose the following: where γ 3 , δ 3 , γ 4 , δ 4 , ξ 3 and 3 are arbitrary constants. Using the above functions in eqn. (51) we get the following algebraic equation, where τ = Aδ 3 + Bδ 4 and ζ = Aγ 3 + Bγ 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.

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. (44). 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 δ 1 for different values of γ 1 and ξ 1 respectively. From Fig. 1 we see that for δ 1 < 0, X 0 remains in the negative level for different  values of γ 1 , thus favouring the formation of a BH. But for δ 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 γ 1 and ξ 1 are almost of the identical nature. In Figs. 3 and 4 plots are generated for X 0 against ξ 1 for different values of δ 1 and γ 1 respectively. In Fig. 3 we see that almost all the trajectories lie in the positive region throughout the entire domain of ξ 1 , thus favouring the formation of NS. Only around ξ 1 ≥ 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 δ 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 γ 1 . Moreover the transition from NS to BH starts at an earlier stage around ξ 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. (52). 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 τ for different values of ζ and ξ 3 respectively. In Fig. 6, for τ < 0, any singularity that forms will be clothed by an event horizon, thus forming a BH. On the contrary, for τ > 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 τ > 0. Similar results are visible in Fig. 7, where trajectories are obtained for different values of ξ 3 . In both the plots the tendency of NS increases with the increase in the values of parameters ζ and ξ 3 .
In Figs. 8 and 9, plots have been obtained for X 0 against the parameter ζ for different values of τ and ξ 3 respectively. In Fig. 8, we see that throughout the domain of ζ trajectories exist both in the positive and negative level. So tendency of formation of both BH and NS exists for any value of ζ and eventually depends on the initial conditions for the final outcome. But with the increase in the value of τ the tendency of formation of NS decreases. From Fig. 9 it is evident that for ξ 3 < 0, collapse results in the formation of NS, whereas for ξ 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 ξ 1 < 0 region we have three positive roots of X 0 for γ 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.

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, where R μν is the Ricci tensor, ψ is a scalar given by ψ = R μν K μ K ν , where K μ = dx μ dτ is the tangent to the non spacelike geodesics at the singularity and τ is the affine parameter. In the paper [ where In ref. [79] it has also been shown that the relation between X 0 and the limiting values of mass is given by, where and ṁ 0 is given by the eqn. (56). It was shown by Dwivedi and Joshi in Ref. [53,80] 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.

Starobinsky's model
Using eqn. (29) in the above relation (54) we get Using eqns. (29), (56) and (58) in eqn. (57) we get an equation for X 0 The above equation can be investigated for roots and using these values of X 0 in the eqn. (59) we get the conditions for which S = lim τ 2 ψ > 0, i.e. the conditions under which we get a strong singularity. For a strong singularity we have from eqn. (59) where X 0 can be obtained as solutions of eqn. (60). Since eqn. (60) 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.
Example. Here we consider the following set of numerical values for the parameters involved, Using eqn. (66) we get the condition for a strong singularity for this model as, where X 0 is obtained from eqn. (67).

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.

Starobinsky's model
The event horizon (r h ) can be obtained from the relation f (t, r) = 0, i.e., 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, For the Starobinsky's model the expression for thermalization temperature becomes, The entropy of the system is given by, where we consider πG = 1. Here G being the Newton's universal gravitation constant. Total energy can be obtained from the relation Using the equations (71), (72) and (73) we get the total energy for this model as, 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, Using the eqns. (71), (72) and (74) we get the expression Helmholtz free energy for this model as, Specific heat at constant volume is given by, Using relations (71), (72), (73) and (77) we get the specific heat at constant volume for the Starobinsky's model as, 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 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 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 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 F(E) = 1. From Fig. 11 it is seen that the range of thermalization temperature T increases with the increase in the value of δ 1 . On the contrary in Fig. 12 it is seen that the range of Internal Energy U decreases with increase in the value of δ 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 δ 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 δ 1 . Finally Fig. 15 shows the effect of gravity's rainbow on the thermodynamic parameter T . 1−a 1 (E s /E P ) . In Fig. 11 the initial conditions are taken as γ 1 = 0.1, ξ 1 = 1, 1 = 5, t = 5, a 1 = 1, E s = 1, E P = 5. In Fig. 12 the initial conditions are γ 1 = 0.1, ξ 1 = 5, 1 = 5, t = 5, a 1 = 1, E s = 1, E P = 5. In both the Figs. 13 and 14 the initial conditions are taken as γ 1 = 0.1, ξ 1 = 0.3, 1 = 3, t = 5, a 1 = 1, E s = 1, E P = 5.

Power law model
The event horizon (r h ) can be obtained from the relation f (t, r) = 0, i.e., The real positive root of the above equation gives the radius of the event horizon. Solving the above equation we get, Obviously for a realistic event horizon radius we should have: For this model the expression for thermalization temperature becomes, Using the equations (81), (72) and (73) we get the total energy for this model as, Using the eqns. (81), (72) and (82) we get the expression Helmholtz free energy for this model as, Finally using relations (81), (72), (82) and (77) we get the specific heat at constant volume for this model as, 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 ≥ 4)). In Fig. 16 we see that the thermalization temperature (T ) decreases with the increase of r. Moreover with the increase of δ 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 δ 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 δ 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 δ 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.
The real positive root of the above equation gives the radius of the event horizon. Here, the expression for thermalization temperature becomes, Using the equations (86), (72) and (73) we get the total energy for this model as, Using the eqns. (86), (72) and (87) we get the expression Helmholtz free energy for this model as, Finally using relations (86), (72), (87) and (77) we get the specific heat at constant volume for this model as, 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 δ 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 δ 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 δ 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 δ 3 , there is a corresponding increase in C. Finally in Fig. 25 we see the effect of gravity's rainbow on the thermalization temperature.

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 ≈ 10 19 GeV) and the length scale tends to the Planck length (L P ≈ 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 ρ of the matter content of the system. This is solely because ρ 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. (42), 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 → 0, t → 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.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.