1 Introduction

Einstein’s general relativity (GR) is the most successful theory of gravity till date, although modifications to GR continue to attract much attention. One of the primary reasons for attempting such modifications has to do with explaining the late time acceleration of the universe. It is known that this phenomenon is compatible with GR in the presence of a dark energy component in the stress tensor. However, much work has been done over the last two decades in trying to explain cosmic acceleration of the universe without invoking dark forms of matter and energy. One such candidate theory is f(R) gravityFootnote 1 (for a sampling of the literature, see the excellent reviews [1,2,3,4,5,6,7]) obtained by modifying the Einstein–Hilbert action to one which includes a regular nonlinear function f(R) of the Ricci scalar R, i.e one in which the Lagrangian density is \(R+f(R)\), apart from the matter part. In this paper, we will deal with the specific model \(f(R) = \alpha R^2\), with \(\alpha \) being a positive constant, a model proposed in [8].

While phenomenological studies of f(R) gravity abound in the literature, there has been relatively lesser focus on collapse scenarios, where matter collapses under its own gravitational force, with the underlying theory being f(R) gravity. We briefly mention a few relevant papers to highlight the progress made thus far. In [9], the collapse process of a star was considered in modified gravity, and it was shown that a class of f(R) theories can result in the prevention of a central singularity in such a process. A generic study of collapse processes of self gravitating dust in f(R) gravity was initiated in [10]. In [11], this process was studied for the case of null dust. In the context of cosmology, collapse in modified gravity was studied in [12], while an extensive numerical analysis for black hole formation in these theories was carried out in [13]. A more recent analysis on collapsing stars in modified gravity was done in [14] (with a generalisation to conformally flat stars appearing in [15]) while results on the collapse of a perfect fluid in f(R) gravity was reported in [16].

As is well known by now, collapse situations in f(R) gravity are greatly restricted compared to their GR counterparts, due to stringent boundary conditions. In GR, such boundary conditions, known as the Darmois–Israel conditions [17, 18] require the first and second fundamental forms to match on the collapsing hypersurface, which is a time-like junction between an internal and an external region of space-time. This guarantees smooth matching of the two regions of space-time, i.e without a stress tensor at the junction. On the other hand, due to higher order nature of f(R) gravity additional conditions have to be imposed [19, 20] (see also [21]) over and above the Darmois–Israel conditions. These often require the Ricci scalar and its (normal) derivative to vanish at the boundary, for smooth matching of the collapsing region with an external Schwarzschild space-time.

This fact was exploited fairly recently in [14] to provide some realistic models of gravitational collapse in f(R) theories in which the coefficient of viscosity is turned off. The starting point of the analysis is the assumption of a specific form of a time dependent spherically symmetric metric, that depends on an arbitrary function of the radial coordinate. The modified Einstein equations in f(R) gravity are then used to constrain these functions in such a way that the extra junction conditions are satisfied, and specific choices give concrete examples of collapse scenarios in f(R) models. Importantly, as pointed out in [14], the additional junction conditions mentioned in [19, 20] render a straightforward generalisation of collapse processes in GR, to scenarios involving modified gravity, difficult. We should emphasise here that in addition to the junction conditions, the collapsing fluid must satisfy various energy conditions that we will elaborate upon in sequel. In totality, all this amounts to the fact that analysing collapsing scenarios in f(R) gravity might be a substantially complicated task.

In this paper, we present new solutions for collapse in f(R) gravity, by assuming some simple ansatzes for the metric, which is then solved by the extra junction conditions, namely the matching of the Ricci scalar and its derivative across a time-like boundary. This the R-matching method commonly used in f(R) collapse scenarios. This is elaborated upon for two cases, first when the metric consists of separable functions of the radial and the time coordinate, and second when it is not. Importantly, the second condition admits shear, and we study this in the presence of a non-zero coefficient of shear viscosity. The R-matching method gives us the full solution of the modified Einstein equations, and we are able to provide a class of realistic collapse models in f(R) gravity, consistent with all energy conditions. For separable solutions to the metric, we are able to provide simple analytic expressions for the components of the energy momentum tensor. These are then used to construct analytic solutions of the heat flux evolution equation. Further, the junction conditions in our model, in the corresponding Einstein scalar theory are analysed in the Jordan frame.

This paper is organised as follows. In the next Sect. 2, after a brief review of the necessary formalism, we write down the general evolution equation of the shear in f(R) theories. The general equation for the evolution of shear is written down and we obtain an algebraic relation between the shear and the anisotropy in f(R) collapse models, via this formula. After this, the necessary energy conditions and the junction conditions of the collapsing fluid are reviewed. With these ingredients, in Sect. 3, we construct a separable solution of the metric using the R-matching method, and show that the end state of collapse is necessarily a black hole. In this case, the collapse is shear-free. Then in Sect. 4, we extend this to non-separable solutions and construct collapsing solutions that obey all energy conditions with the end state being a (locally) naked singularity. The role of shear is commented upon, in this example. In Sect. 5, we study some physical properties of the collapsing fluid, for the separable case. The nature of the equation of state is commented upon, and the heat flux evolution equation is solved under some assumptions to clearly highlight the role of the f(R) parameter. In Sect. 6, we comment on the description of our model in the Einstein scalar theory, and revisit the junction conditions in the Jordan frame. Finally, Sect. 7 ends this paper with a summary of the main results and some discussions.

2 Mathematical preliminaries and set up

For a generic collapse scenario, in co-moving coordinates \((t,r,\theta ,\phi )\) the metric inside the spherically symmetric collapsing cloud is written as

$$\begin{aligned} ds^{2}_{-}=-e^{2\nu (r,t)}dt^{2}+e^{2\psi (r,t)}dr^{2}+Q^{2}(r,t)d\Omega ^{2}, \end{aligned}$$
(1)

where \(d\Omega ^2 = d\theta ^2 + \sin ^2\theta d\phi ^2\). The metric outside the collapsing matter is usually represented by the Vaidya solution in terms of the retarded time u as

$$\begin{aligned} ds_+^2 = -\left( 1-\frac{2m(u)}{{\tilde{r}}}\right) du^2 - 2dud{{\tilde{r}}} + {{\tilde{r}}}^2d\Omega ^2 \end{aligned}$$
(2)

In this paper, we will be interested in an exterior vacuum solution (i.e without any radiation) and hence with m(u) being a constant, the metric out side the collapsing matter can be taken to be the Schwarzschild metric, given by

$$\begin{aligned} ds^{2}_{+}= & {} -H(\tilde{r})d\tilde{t}^{2}+H(\tilde{r})^{-1}d\tilde{r}^{2}+\tilde{r}^{2}d\Omega ^{2},\nonumber \\ H(\tilde{r})= & {} 1-\frac{2m}{\tilde{r}} , \end{aligned}$$
(3)

where m is the (constant) Schwarzschild mass, so that the heat flux obtained from Eq. (1) is zero at the matching hypersurface.

The modified Einstein’s equations for a Lagrangian density \(R+f(R) + {\mathcal L}_{matter}\) are given by

$$\begin{aligned} G_{\mu \nu }= & {} \frac{1}{1+F}\left( T_{\mu \nu }+D_{\mu \nu }F(R)+\frac{1}{2}g_{\mu \nu }(f-RF)\right) ,\nonumber \\ T_{\mu \nu }= & {} -\frac{1}{\sqrt{-g}}\frac{\delta {\mathcal L}_{matter}}{\delta g^{\mu \nu }} \end{aligned}$$
(4)

where the Einstein tensor \(G_{\mu \nu } = R_{\mu \nu } - \frac{1}{2}g_{\mu \nu }R\), g is the determinant of the metric, and

$$\begin{aligned} F(R)= & {} \frac{df(R)}{dR},~D_{\mu \nu }=\nabla _{\mu }\nabla _{\nu }-g_{\mu \nu }\nabla _{\alpha }\nabla ^{\alpha }. \end{aligned}$$
(5)

We need to solve the modified Einstein equations with the energy momentum tensorFootnote 2

$$\begin{aligned} T_{\mu \nu }=\rho u_{\mu }u_{\nu } + P h_{\mu \nu } -\Pi _{\mu \nu } + 2qu_{(\mu }n_{\nu )} - 2\eta \sigma _{\mu \nu }, \end{aligned}$$
(6)

where we define the quantities

$$\begin{aligned} \sigma _{\mu \nu }= & {} u_{(\mu ;\nu )} + a_{(\mu }u_{\nu )} - \frac{1}{3}\Theta \left( g_{\mu \nu } + u_{\mu }u_{\nu }\right) ,\nonumber \\ \Pi _{\mu \nu }= & {} \Pi \left( n_{\mu }n_{\nu } - \frac{1}{3}h_{\mu \nu }\right) ,\nonumber \\ \Pi= & {} p_{\theta }-p_r ,~~ P = \frac{p_r + 2p_{\theta }}{3}, \end{aligned}$$
(7)

where (, ) denote a symmetrization, and a semicolon denotes a covariant derivative. Here, \(\rho \) is the energy density, \(p_r\) and \(p_{\theta }\) are the radial and tangential pressures, respectively, \(q^{\mu } = qn^{\mu }\) is the radial heat flow vector where \(n^{\mu }\) is a unit 4-vector along the radial direction, and \(u^{\mu }\) is the 4-velocity of the fluid. The quantity, \(a^{\mu } = u^{\mu }_{;\nu }u^{\nu }\) designates the 4-acceleration of the fluid. These satisfy \(n^{\mu }n_{\mu } = 1\), \(u^{\mu }u_{\mu } = -1\), \(u^{\mu }q_{\mu }=0\), \(u^{\mu }n_{\mu }=0\). Also, \(\Theta = u^{\nu }_{;\nu }\) is the expansion parameter and \(h_{\mu \nu } = g_{\mu \nu } + u_{\mu }u_{\nu }\) is the projection tensor.

In this paper, we will be dealing with two situations, to be elaborated in Sects. 3 and 4. In the former, we will consider shear-free collapse, with the fluid being non-geodesic. In the latter, we will consider a geodesic fluid, but with non-zero shear. It will therefore be useful for us to record the relations that connect these quantities, in the f(R) model that we consider. As we will see, we are led to some useful insights here.

To begin with, we record the Raychaudhuri equation, which reads (see, e.g [22])

$$\begin{aligned}&u^{\alpha }\Theta _{;\alpha }+\frac{1}{3}\Theta ^{2}+\frac{2}{3}\sigma ^{2}-a^{\mu }_{;\mu }\nonumber \\&\quad +\frac{1}{1+F} \left[ -\frac{1}{2}\left( R+f(R)\right) \right. \nonumber \\&\quad \left. +\left( \frac{dF}{dR}\right) h^{\mu \nu }R_{;\mu \nu }+T_{\mu \nu }u^{\mu }u^{\nu } \right] =0, \end{aligned}$$
(8)

where \(\sigma _{\alpha \beta }\sigma ^{\alpha \beta }=\frac{2}{3} \sigma ^{2}\). This equation is valid only for f(R) models with \(d^{2}F/dR^{2}=0\), which is the case under consideration here.Footnote 3 Now, we will use the identity given by [25]

$$\begin{aligned}&u^{\beta }u_{\rho }R^{\rho }_{\alpha \beta \mu }h^{\alpha }_{\gamma }h^{\mu }_{\nu }=h^{\alpha }_{\gamma }h^{\mu }_{\nu } (a_{\alpha ;\mu } - u^{\beta }\sigma _{\alpha \mu ;\beta })\nonumber \\&\quad + a_{\gamma }a_{\nu }-\frac{1}{3}u^{\beta }\Theta _{;\beta }h_{\gamma \nu }-\frac{1}{9}\Theta ^{2}h_{\gamma \nu }-\frac{2}{3}\Theta \sigma _{\mu \nu }\nonumber \\&\quad -\frac{\sigma ^{2}}{3} \left( n_{\gamma }n_{\nu }+\frac{1}{3}h_{\gamma \nu } \right) , \end{aligned}$$
(9)

and the well known relation between Riemann and Weyl tensors given by

$$\begin{aligned} R^{\mu }_{\nu \rho \sigma }= & {} C^{\mu }_{\nu \rho \sigma } + \frac{1}{2}\left( R^{\mu }_{\rho }g_{\nu \sigma } - R^{\mu }_{\sigma }g_{\nu \rho } - R_{\nu \rho }\delta ^{\mu }_{\sigma } + R_{\nu \sigma }\delta ^{\mu }_{\rho }\right) \nonumber \\&+\frac{R}{6}\left( \delta ^{\mu }_{\sigma }g_{\nu \rho } - \delta ^{\mu }_{\rho }g_{\nu \sigma }\right) . \end{aligned}$$
(10)

For f(R) gravity (recall that \(f(R) = \alpha R^2\) with \(d^{2}F/dR^{2}=0\)), Eq. (10) can be evaluated by using the results derived in [22] and after some algebra, we obtain (with \(R_{;\mu ;\nu } \equiv R_{;\mu \nu }\)),

$$\begin{aligned} u^{\beta }u_{\rho }R^{\rho }_{\alpha \beta \mu }h^{\alpha }_{\gamma }h^{\mu }_{\nu }= & {} E_{\gamma \nu }+\frac{1}{2(1+F)} \left[ -\frac{\left( R+f(R)\right) }{3}h_{\gamma \nu }\right. \nonumber \\&+\frac{dF}{dR} \left( h_{\gamma \nu }h^{\alpha \beta }R_{;\alpha \beta }-h^{\alpha }_{\gamma }h^{\mu }_{\nu }R_{;\alpha \mu } \right) \nonumber \\&+ \left. \left( \frac{2}{3}\rho h_{\gamma \nu } +\Pi _{\gamma \nu }+2\eta \sigma _{\gamma \nu } \right) \right] . \end{aligned}$$
(11)

This generalises a corresponding result obtained for GR in [25]. Here, \(E_{\mu \nu }\) is the electric part of the Weyl tensor, defined as \(E_{\mu \nu } = C_{\mu \nu \rho \lambda }u^{\rho }u^{\lambda }\), with the magnetic part of the Weyl tensor vanishing identically due to spherical symmetry. Then, eliminating \(\rho \) from Eq. (11) using Eq. (8), we obtain

$$\begin{aligned}&u^{\beta }u_{\rho }R^{\rho }_{\alpha \beta \mu }h^{\alpha }_{\gamma }h^{\mu }_{\nu }\nonumber \\&\quad = E_{\gamma \nu }+\frac{1}{2(1+F)} \left[ \frac{dF}{dR} \left( \frac{1}{3}h_{\gamma \nu }h^{\alpha \beta } R_{;\alpha \beta }\right. \right. \nonumber \\&\qquad \left. \left. -h^{\alpha }_{\gamma }h^{\mu }_{\nu }R_{;\alpha \mu }\right) + {{\hat{P}}}_{\mu \nu } \right] \nonumber \\&\qquad -\frac{1}{3}h_{\gamma \nu } \left( u^{\alpha }\Theta _{;\alpha }+\frac{1}{3}\Theta ^{2}-a^{\mu }_{;\mu }+\frac{2}{3}\sigma ^{2} \right) , \end{aligned}$$
(12)

where we have defined

$$\begin{aligned} E_{\mu \nu } = {\mathcal E}\left( n_{\mu }n_{\nu } - \frac{1}{3}h_{\mu \nu }\right) ,~~{{\hat{P}}}_{\mu \nu } = (\Pi _{\gamma \nu }+2\eta \sigma _{\gamma \nu }) \end{aligned}$$
(13)

Equating Eqs. (12) and (9) we get

$$\begin{aligned}&h^{\alpha }_{\gamma }h^{\mu }_{\nu }\left( a_{\alpha ;\mu } - u^{\beta }\sigma _{\alpha \mu ;\beta }\right) +a_{\gamma }a_{\nu }-\frac{1}{3}\sigma _{\gamma \nu }(2\Theta +\sigma )\nonumber \\&\quad = E_{\gamma \nu }+\frac{1}{2(1+F)} \left[ \frac{dF}{dR} \left( \frac{1}{3}h_{\gamma \nu }h^{\alpha \beta } R_{;\alpha \beta }\right. \right. \nonumber \\&\qquad \left. \left. -h^{\alpha }_{\gamma }h^{\mu }_{\nu }R_{;\alpha \mu } \right) +{{\hat{P}}}_{\mu \nu } \right] +\frac{1}{3}h_{\gamma \nu }a^{\mu }_{;\mu }. \end{aligned}$$
(14)

Finally, contracting with \(n^{\gamma }n^{\nu }\) and denoting \({{\hat{P}}} = \left( \Pi + 2\eta \sigma \right) \),

$$\begin{aligned}&n^{\alpha }n^{\mu }\left( a_{\alpha ;\mu }-u^{\beta }\sigma _{\alpha \mu ;\beta }+a_{\alpha }a_{\mu }\right) -\frac{2}{9}\sigma (2\Theta +\sigma )\nonumber \\&\quad =\frac{2}{3}\mathcal {E} +\frac{1}{2(1+F)} \left[ \frac{dF}{dR} \left( \frac{1}{3}h^{\alpha \mu }-n^{\alpha }n^{\mu } \right) R_{;\alpha \mu }+\frac{2}{3}{\hat{P}} \right] \nonumber \\&\qquad +\frac{1}{3}a^{\mu }_{;\mu }. \end{aligned}$$
(15)

Expanding the left hand side of Eq. (15), the evolution of the shear is given by the equation

$$\begin{aligned}&e^{-\psi }\frac{da}{dr} - \frac{2}{3}e^{-\nu }\frac{d\sigma }{dt} + a^2 -\frac{2}{9}\sigma (2\Theta +\sigma ) \nonumber \\&\quad =\frac{2}{3}\mathcal {E} +\frac{1}{2(1+F)} \left[ \frac{dF}{dR} \left( \frac{1}{3}h^{\alpha \mu }-n^{\alpha }n^{\mu } \right) R_{;\alpha \mu }+\frac{2}{3}{\hat{P}} \right] \nonumber \\&\qquad +\frac{1}{3}a^{\mu }_{;\mu }, \end{aligned}$$
(16)

with \(a=n^{\mu }a_{\mu }\). Equation (16) is the most general evolution equation for the shear tensor in \(f(R)=\alpha R^2\) scenarios, with \(d^2F/dR^2=0\). The GR case corresponds here to \(\alpha = 0\) and has been analysed in [25]. We can make a few comments here. Now note that \(\sigma \) (being computed entirely from the metric) does not depend on the f(R) parameter \(\alpha \). This means that the term in square brackets in Eq. (16) has to be independent of \(\alpha \). For \(f(R)=\alpha R^2\) theories, this can be seen to imply that

$$\begin{aligned} \sigma= & {} \frac{3}{4 \eta R}\left( \frac{1}{3}h^{\alpha \mu }-n^{\alpha }n^{\mu }\right) R_{;\alpha \mu } - \frac{\Pi }{2\eta }\nonumber \\&+ \,\frac{\left( 1+ 2\alpha R\right) }{4\eta R}\left( \frac{\partial \Pi }{\partial \alpha }\right) \end{aligned}$$
(17)

Eq.(17) gives an algebraic relation between the shear and the anisotropy in the f(R) theories that we consider.Footnote 4 To the best of our knowledge, Eqs. (16) and (17) have not appeared in the literature before, and provide useful insights into the dynamics of f(R) collapse. These equations will be identically satisfied in the explicit solutions that we will construct in sequel.

The next ingredient in our analysis will be the relevant energy conditions of the collapsing fluid. In this context, we begin from the energy momentum tensor of Eq. (6), that describes the motion of a fluid with shear, with heat flow in the radial direction. The energy conditions for such a fluid including the effects of anisotropy was obtained in [24] (see also [23]) by generalising a method developed in [26] for isotropic cases. This essentially relies on the fact that the eigenvalues of the energy momentum tensor should be real, and the resulting conditions on the fluid are given by

$$\begin{aligned} \mathrm{(i)}&~&|\rho + p_r - 2\eta \sigma _{11}| -2|q| \ge 0, \nonumber \\ \mathrm{(ii)}&~&\rho - p_r + 2p_{\theta } + \Delta + 2\eta \left( \sigma _{11} - 2\sigma _{22}\right) \ge 0\nonumber \\ \mathrm{(iii)}&~&\rho - p_r + 2p_{\theta } + \Delta + 2\eta \left( \sigma _{11} - 2\sigma _{33}\right) \ge 0 \end{aligned}$$
(18)

where we have defined

$$\begin{aligned} q= -\frac{T_{01}}{\sqrt{-g_{tt}~g_{rr}}},~~\Delta =\sqrt{(\rho +p_{r}-2\eta \sigma _{11})^{2}-4q^{2}}. \end{aligned}$$
(19)

In addition, the weak, dominant and strong energy conditions (WEC, DEC and SEC) are to be satisfied, and these are given respectively as

$$\begin{aligned} \mathrm{(iv)}&~&\rho - p_r + \Delta +2\eta \sigma _{11} \ge 0 ~~(\mathrm{WEC}) \nonumber \\ \mathrm{(v)}&~&\rho - p_r + 2\eta \sigma _{11} \ge 0 ~~(\mathrm{DEC1}) \nonumber \\ \mathrm{(vi)}&~&\rho - p_r -2p_{\theta } + \Delta +2\eta \left( \sigma _{11} + 2\sigma _{22}\right) \ge 0~~ (\mathrm{DEC2}) \nonumber \\ \mathrm{(vii)}&~&\rho - p_r -2p_{\theta } + \Delta +2\eta \left( \sigma _{11} + 2\sigma _{33}\right) \ge 0 ~~(\mathrm{DEC3}) \nonumber \\ \mathrm{(viii)}&~&2p_{\theta } + \Delta - 2\eta \left( \sigma _{22} + \sigma _{33}\right) \ge 0 ~~(\mathrm{SEC}) \end{aligned}$$
(20)

where the DEC consists of three separate conditions labeled DEC1, DEC2 and DEC3. For convenience, we record the above conditions in the case of vanishing shear, and they read,

$$\begin{aligned}&\hbox {I}.~|\rho +p_{r}|-2|q|\ge 0,~ \hbox {II}.~\rho -p_{r}+2p_{\theta }+\Delta \ge 0. \end{aligned}$$
(21)
$$\begin{aligned}&\hbox {III}.~\rho -p_{r}+\Delta \ge 0,~ \hbox {IV}\mathrm{A}.~\rho -p_{r}\ge 0,\nonumber \\&\hbox {IV}\mathrm{B}.~\rho -p_{r}-2p_{\theta }+\Delta \ge 0, \hbox {V}.~~2p_{\theta }+\Delta \ge 0. \end{aligned}$$
(22)

Finally, all the conditions above will need to be supplemented by the junction conditions in f(R) models [19, 20]. Recall that in GR, the standard Darmois–Israel junction conditions [18] are valid, which amount to matching of the first and second fundamental forms at the time-like hypersurface \(\Sigma : r = r_0\). These are defined, with a, b denoting the indices on the hypersurface, as,

$$\begin{aligned} g_{ab}= & {} g_{\alpha \beta }e^{\alpha }_{a}e^{\beta }_{b},\nonumber \\ K_{ab}= & {} \frac{1}{2}\mathcal {L}_{N}g_{ab}=\frac{1}{2}\left( g_{ab,c}N^{c}+g_{cb}N^{c}_{,a}+g_{ac}N^{c}_{,b}\right) , \end{aligned}$$
(23)

where \(N^{\mu }\) is the unit normal across the matching hypersurface. Here, \(e^{\alpha }_{a}=\frac{\partial x^{\alpha }}{\partial y^{a}}\) are tangents to the matching hypersurface, and \(\mathcal {L}_{N}g_{ab}\) is the Lie derivative of the induced metric with respect to the normal vector to the hypersurface. For f(R) collapse, the additional requirements are the continuity of the Ricci scalar and its first derivative across this hypersurface [19, 20], so that the full set of matching conditions across the collapsing time-like hypersurface separating \(ds^{2}_{-}\) and \(ds^{2}_{+}\) are

$$\begin{aligned} \left[ g_{ab}\right] =0,~\left[ K_{ab}\right] =0,~[R]=0,~N^{\mu }[\partial _{\mu } R] = 0, \end{aligned}$$
(24)

where [a] denotes the difference in the quantity a across the hypersurface \(\Sigma \). Therefore, in studying any model of f(R) collapse, we will need to impose the junction conditions of Eq. (24), in addition to the energy conditions spelt out in Eqs. (18) and (20).

The first two relations of Eq. (24) are fairly straightforward to deal with. Since the analysis is standard, we will not go into the details here, but simply state that these imply that the Misner–Sharp mass function [27, 28] given by

$$\begin{aligned} M(r,t)=\frac{Q}{2}\left[ 1-e^{-2\psi }\left( \frac{dQ}{dr}\right) ^{2}+e^{-2\nu }\left( \frac{dQ}{dt}\right) ^{2}\right] , \end{aligned}$$
(25)

equals the Schwarzschild mass when evaluated at the boundary \(\Sigma \). From a fairly straightforward analysis, it is known that these also imply, from the metrics of Eqs. (1) and (2) that

$$\begin{aligned}&\frac{Q}{2}e^{-\left( \nu +\psi \right) }\left[ 2\frac{{\dot{Q}}'}{Q} -2 \frac{{\dot{Q}}}{Q}\frac{{{{\dot{\psi }}}}}{\psi } -2\frac{\nu '}{\nu }\frac{{\dot{Q}}}{Q}\right. \nonumber \\&\quad + e^{\left( \psi - \nu \right) }\left( 2\frac{{\ddot{Q}}}{Q} - 2\frac{{\dot{Q}}}{Q}\frac{{{{\dot{\nu }}}}}{\nu } + \frac{e^{2\nu }}{Q^2} + \frac{{\dot{Q}}^2}{Q^2} \right. \nonumber \\&\quad \left. \left. - e^{2\left( \nu - \psi \right) }\left( \frac{Q'^2}{Q^2} - 2\frac{\nu '}{\nu }\frac{Q'}{Q}\right) \right) \right] \biggr |_{\Sigma }=0 \end{aligned}$$
(26)

Equivalently, the junction conditions imply that [21]

$$\begin{aligned} N^{\mu }\left[ T_{\mu \nu }\right] = 0, \end{aligned}$$
(27)

which is a familiar condition in GR. The other two relations of Eq. (24) are the essential new ingredients in this analysis. In summary, our task is to study collapse in f(R) gravity, that are restricted by eight conditions mentioned in Eqs. (18) and (20) in addition to the four junction conditions spelt out in Eq. (24). Indeed, this seems to be a formidable task, especially in cases with shear, but as we elaborate upon below, some simple solutions can nonetheless be found by utilising the constraints of Eq. (24).

3 Separable interior solutions

The extra junction conditions in f(R) gravity are in fact quite strong, and can potentially exclude several well known collapse solutions in GR. For example, the Oppenheimer–Snyder solution is not an admissible collapsing solution in modified gravity scenarios [20]. As another concrete example, suppose we assume that the interior metric is of the Lemaitre–Tolman–Bondi (LTB) form [29,30,31] given by

$$\begin{aligned} ds_-^2 = -dt^2 + X^2\left( r,t\right) dr^2 + Q^2\left( r,t\right) d\Omega ^2 \end{aligned}$$
(28)

where X(r, t) and Q(r, t) are functions of the co-moving radial coordinate and time, and \(d\Omega ^2\) is the metric on the unit two-sphere. The special case of the homogeneous Friedmann–Robertson–Walker metric is obtained from Eq. (28) by writing

$$\begin{aligned} X(r,t) = \frac{a(t)}{\sqrt{1-kr^2}},~Q(r,t) = a(t) r \end{aligned}$$
(29)

with k being a suitable constant. Also, the Einstein equations of GR can be shown to imply, for the general metric of Eq. (28),

$$\begin{aligned} X(r,t) = {\mathcal A}(r)Q'(r,t), \end{aligned}$$
(30)

with \({\mathcal A}(r)\) being an arbitrary function of the co-moving radial coordinate.Footnote 5 We will consider these two cases separately.

We first consider a separable solution for the interior metric, of the form given in Eq. (29), and assume that in co-moving coordinates, this is

$$\begin{aligned} ds_-^{2}= & {} -dt^2 + \frac{a(t)^{2}}{h(r)}dr^2 + a(t)^{2}r^{2}d\theta ^2\nonumber \\&+\, a(t)^{2}r^{2}\sin ^2\theta d\phi ^2, \end{aligned}$$
(31)

Since we are in co-moving coordinates we choose \(u^{\mu }=(1,0,0,0)\) and \(n^{\mu }=(0,\sqrt{h(r)}/a(t),0,0)\), so that the heat flux is along the radial direction, i.e \(q^{\mu }=qn^{\mu }\). With this metric, for the model described the Lagrangian density \(R+f(R)=R+\alpha R^{2}\) (where \(\alpha \) is a constant) we can write down the energy momentum components:

$$\begin{aligned} \rho= & {} (1+F)G_{00}-\frac{\alpha R^{2}}{2}\nonumber \\&\quad -\frac{1}{2ra^{2}} \left( 2rhF^{\prime \prime }+(4h+rh^{\prime }F^{\prime }) -6ra\dot{a}\dot{F} \right) , \nonumber \\ \frac{a^{2}}{h}p_{r}= & {} (1+F)G_{11}+\frac{a^{2}}{h}\frac{\alpha R^{2}}{2}+\frac{2F^{\prime }}{r}-\frac{a}{h} \left( a\dot{a} \dot{F}+a\ddot{F} \right) , \nonumber \\ a^{2}r^{2}p_{\theta }= & {} (1+F)G_{22}+\frac{\alpha R^{2}}{2}a^{2}r^{2}\nonumber \\&-\,\frac{r}{2} \left( 4a\dot{a}\dot{F}+2ra^{2}\ddot{F}-(2h+rh^{\prime })F^{\prime }-\,2rhF^{\prime \prime } \right) ,\nonumber \\ \frac{a}{\sqrt{h}}q= & {} \dot{F}^{\prime }-\frac{\dot{a}}{a}F^{\prime }. \end{aligned}$$
(32)

The Ricci scalar of the interior metric is calculated to be

$$\begin{aligned} R(r,t)=\frac{2}{a^{2}}\left( \frac{1-h-rh^{\prime }}{r^{2}}\right) +6\frac{\left( \dot{a}^{2}+a\ddot{a}\right) }{a^{2}}. \end{aligned}$$
(33)

In order that the Ricci scalar matches smoothly to the collapsing co-moving boundary at all co-moving times, we will therefore require that \(\dot{a}^{2}+a\ddot{a} = 0\) (since the second term on the right hand side of Eq. (33) is a function of time only), in which case the first term of Eq. (33) can be appropriately solved in order to fulfil the requirement that R is continuous across the matching hypersurface. However, to satisfy Eq. (27), one finds after a straightforward calculation, using the unit normal vector \(N^{\mu } = \left( 0, \sqrt{h(r)}/a(t),0,0\right) \), the condition

$$\begin{aligned} h(r) -1-r^2 (\dot{a}^{2}+2a\ddot{a})=0. \end{aligned}$$
(34)

In order to satisfy this for all times, one thus requires \(\dot{a}^{2}+2a\ddot{a} = 0\) which naturally implies that this cannot be satisfied in conjunction with the criterion for a continuous Ricci scalar across the boundary, at all co-moving times. In conclusion, what we have here is a no go scenario, namely that a simple separable form of the metric given in Eq. (31) is unsuitable for describing collapse in f(R) gravity.

The assumption of a separable solution of the form in Eq. (31) is possibly an over-simplification. We will next consider another separable form of the interior metric given by

$$\begin{aligned} ds_-^2= & {} -A(r)^2dt^2 + 2a(t)^2\left( \partial _r A(r)\right) ^2dr^2\nonumber \\&+ a(t)^2A(r)^2d\Omega ^2, \end{aligned}$$
(35)

with the energy momentum tensor having the same form as in Eq. (6). We will match this with an external Schwarzschild solution. This metric was originally considered in [32] to the study the collapse of a shear-free radiating spherically symmetric star in GR. As we elaborate below, this ansatz offers considerable simplifications in the study of collapsing stars in f(R) gravity.

To this end, we first note that the Ricci scalar here is given by

$$\begin{aligned} R = -\frac{1-6\left( {\dot{a}(t)}^2 + a(t){\ddot{a}(t)}\right) }{a(t)^2A(r)^2} \end{aligned}$$
(36)

Using Eqs. (4) and (6), the relevant physical quantities are obtained for the metric of Eq. (35) as

$$\begin{aligned} A(r)^{2} \rho= & {} (1+F) \left( \frac{1+6\dot{a}^{2}}{2a^{2}} \right) +\frac{\alpha A^{2} R^{2}}{2}\nonumber \\&- \left( \frac{A\left( (2A^{\prime 2}-AA^{\prime \prime })F^{\prime }+AA^{\prime }F^{\prime \prime } \right) -6a\dot{a}A^{\prime 3}\dot{F}}{2a^{2}A^{\prime 3}} \right) ,~\nonumber \\ 2a^{2}A^{\prime 2} p_{r}= & {} (1+F) \left( \frac{A^{\prime 2}(1-2\dot{a}^{2}-4a\ddot{a})}{A^{2}} \right) \nonumber \\&+\alpha a^{2}A^{\prime 2}R^{2}- \left( \frac{A^{\prime } \left( 2aA^{\prime }(2\dot{a}\dot{F}+a\ddot{F})-3AF^{\prime } \right) }{A^{2}} \right) ,~\nonumber \\ a^{2}A^{2}p_{\theta }= & {} (1+F) \left( \frac{1-2\dot{a}^{2}-4a\ddot{a}}{2} \right) +\frac{\alpha R^{2}a^{2}A^{2}}{2}-a(2\dot{a}\dot{F}\nonumber \\&+ a\ddot{F})- \left( \frac{A \left( (AA^{\prime \prime }-2A^{\prime 2})F^{\prime }-AA^{\prime }F^{\prime \prime } \right) }{2A^{\prime 3}} \right) , \nonumber \\ -\sqrt{2}aAAq= & {} (1+F)\frac{2\dot{a}A^\prime {}}{aA}+\frac{A^{\prime }\dot{F}}{A}+\frac{\dot{a}F^{\prime }}{a}-\dot{F}^{\prime }. \end{aligned}$$
(37)

Now, from the metric of Eq. (35), it can be checked via Eq. (37) that the pressure anisotropy is given by

$$\begin{aligned} p_{\theta } - p_r = -\frac{8\alpha \left[ 1-6\left( {\dot{a}(t)}^2 + a(t){\ddot{a}(t)}\right) \right] }{a(t)^4A(r)^4}~ \equiv -\frac{8\alpha }{a^2A^2}R \nonumber \\ \end{aligned}$$
(38)

Importantly, if we demand that the pressure anisotropy vanishes identically, then it necessarily implies that \(\alpha =0\), in which case the solution reduces to one in GR. We are therefore naturally constrained to consider situations with pressure anisotropy in f(R) scenarios.

It is then seen that in order for the Ricci scalar and its derivative to be continuous across the matching hypersurface (which we choose without loss of generality to be \(r = 1\)), it is enough for us to choose \(A(r) = (1-r)^{-n}\) with \(n\ge 1\), so that continuity of the Ricci scalar and its derivative is guaranteed at the boundary. The function a(t) is unspecified at this stage. In order to simplify the computations, we will have to make a choice, and to this end we will choose \({\dot{a}(t)}^2 + a(t){\ddot{a}(t)}=0\). To summarise, our ansatz for a solution of the metric of Eq. (35) is (with b and n being constants),

$$\begin{aligned} a(t) = \sqrt{1-2bt},~~A(r) = \frac{1}{(1 - r)^n},~~n \ge 1 ,~~b>0. \end{aligned}$$
(39)

We will henceforth choose for simplicity, the constants \(b = 1/2\) and \(n=2\) so that \(R=-\left( 1-r\right) ^4/\left( 1-t\right) \) and satisfies both the conditions on the Ricci scalar mentioned in Eq. (24) at the boundary, arbitrarily close to the time of collapse. In this notation, the collapse starts at \(t=0\) and a singularity forms at \(t=1\) where the scale factor a(t) goes to zero and the Ricci scalar diverges although all co-moving observers see an apparent horizon at \(t=1/2\), as we elaborate in a while. As usual, our interior solution is matched to an external Schwarzschild metric at \(r=1\).

Now, using the fact that the hypersurface normal is given by the vector

$$\begin{aligned} N^{\mu } = \left( 0,\frac{1}{\sqrt{2} a(t)A'(r)},0,0\right) , \end{aligned}$$
(40)

it can be immediately seen that the condition \(N^{\mu }\left[ T_{\mu \nu }\right] =0\) is satisfied at all times. Now upon using Eqs. (36) and (39), we finally obtain the very simple expressions,

$$\begin{aligned} \rho= & {} \frac{\left( 1-r\right) ^4 \left( 5-7t+2t^2 +2 \alpha (1-r)^4 (4-t)\right) }{4 \left( 1-t\right) ^3},\nonumber \\ p_r= & {} \frac{(1-r)^4 \left( 3-5t + 2t^2 + 2 \alpha (1-r)^4 (14-11t)\right) }{4 \left( 1-t\right) ^3},~\nonumber \\ p_{\theta }= & {} \frac{(1-r)^4 \left( 3-5t + 2t^2 - 2 \alpha (1-r)^4 (2-5t)\right) }{4 \left( 1-t\right) ^3},\nonumber \\ q= & {} -\frac{(1-r)^4 \left( t-1-6 \alpha (1-r)^4\right) }{\sqrt{2} (1-t)^{5/2}}. \end{aligned}$$
(41)

It is clearly seen from Eq. (41) that all the components of the stress tensor vanish at the boundary \(r=1\), and that the pressure and density are positive for all values of the co-moving radius, at all co-moving times.Footnote 6 This situation thus corresponds to the realistic collapse of a dense star in f(R) gravity.

This last statement requires some clarification. From our discussion above, it follows that the collapse reaches a singularity in co-moving time \(t=1\), when \(a(t)=0\) and the Ricci scalar diverges at all co-moving radii. This is a shell focusing singularity, which happens simultaneously for all co-moving observers. In order to determine whether the singularity is naked or not, we have to investigate the formation of trapped surfaces during the collapse process. These are the compact two-dimensional space-like surfaces such that both families of ingoing and outgoing null geodesics orthogonal to them necessarily converge. Mathematically one can find out such locations from the expansion parameter \(\Theta \) of the outgoing future-directed null geodesics. We consider a congruence of outgoing radial null geodesics having the tangent vector \((u^t, u^r, 0, 0)\). If such geodesics terminate at the singularity in the past with a definite tangent vector, then at singularity we have \(\Theta > 0\). When such curves do not exist it means that an event horizon has formed earlier than singularity, thus forming a blackhole as the end stage of the collapse process.

Fig. 1
figure 1

Condition I

Full size image
Fig. 2
figure 2

Condition II

Full size image
Fig. 3
figure 3

Condition III

Full size image
Fig. 4
figure 4

Condition IVA

Full size image
Fig. 5
figure 5

Condition IVB

Full size image
Fig. 6
figure 6

Condition V

Full size image

Now recall that for a spherically symmetric metric such as the one we are considering, the co-moving time of formation of an apparent horizon is given from the equation

$$\begin{aligned} g^{\mu \nu }\partial _{\mu }Q(r,t)\partial _{\nu }Q(r,t) =0 \end{aligned}$$
(42)

Using Eq. (35) and the ansatz of Eq. (39), it is then checked that the co-moving time formation of the apparent horizon for all co-moving observers is \(t=1/2\). Hence, the end state of the collapse process is a black hole in this case. This is also obtained by computing the boundary redshift for an observer at infinity, which diverges at the formation time of the black hole. This is obtained by writing the external Schwarzschild solution of Eq. (3) in terms of retarded time and computing the junction conditions, and the co-moving time for our collapsing scenario at which the redshift at infinity diverges is [23, 24]

$$\begin{aligned} \frac{1}{\sqrt{2}} - \frac{1}{2\sqrt{1-t}} = 0,~~ \end{aligned}$$
(43)

which yields the same result \(t=1/2\).

It remains to check the validity of the energy conditions listed in Eqs. (21) and (22). This is most conveniently done numerically, since analytical expressions for these conditions become cumbersome. In this analysis,we choose \(\alpha = 10^{-3}\). In Figs. 1, 2, 3, 4, 5, 6, we show that all the energy conditions are indeed satisfied. In all these figures, the solid red, dotted blue and dashed black curves indicate the co-moving observer at \(r=0.1\), 0.5 and 0.9, respectively. We have shown the validity of the energy conditions from \(t=0\) the \(t=1\), although it is to be noted that as we have discussed, the apparent horizon forms at \(t=1/2\) for this model.

Here, the four-velocity, and the unit vector in the radial direction are

$$\begin{aligned} u^{\mu } = \left( (1-r)^2,0,0,0\right) ,~~n^{\mu }=\left( 0,\frac{(1-r)^3}{2\sqrt{2}\sqrt{1-t}},0,0\right) \nonumber \\ \end{aligned}$$
(44)

These will satisfy the conditions \(u^{\mu }u_{\mu } = -1\) and \(n^{\mu }n_{\mu }=1\), along with those mentioned after Eq. (6). The shear tensor is identically zero in this case, as is generally true for separable solutions of the form that we consider here. It is also straightforward to check that the expansion scalar for a time-like congruence is given by

$$\begin{aligned} \Theta = -\frac{3\left( 1-r\right) ^2}{2\left( 1-t\right) }, \end{aligned}$$
(45)

Also, using Eq. (44), it can be checked that for our metric of Eq. (35), the condition of Eq. (15) is satisfied with

$$\begin{aligned} {\mathcal E} = -\frac{\left( 1-r\right) ^4}{\left( 1-t\right) },~~ \Pi = -\frac{8\alpha \left( 1-r\right) ^8}{\left( 1-t\right) ^2}. \end{aligned}$$
(46)

With these inputs, it can be checked that Eq. (16) is indeed satisfied in this case, with \(\sigma = 0\), and so is Eq. (17).

Note that here the pressure anisotropy goes to zero at the matching hypersurface as it should, but does not vanish at the origin (\(r=0\)). Interestingly, this is an artefact of f(R) gravity, as the anisotropy vanishes identically with \(\alpha = 0\), as follows from Eq. (38) or (46). In this context, we note that anisotropy in static situations (for example in compact stars) have been studied extensively (see, e.g. [34, 42]). It is well known that in such static situations, the pressure anisotropy must vanish at the center, and that a non-zero central anisotropy implies that the density at the center vanishes [35]. These conditions need not be satisfied in non-equilibrium situations that we are considering here. In this context, observe from Eq. (17) that since the shear is identically zero in this case, the anisotropy at the center is forced to be non-zero, since none of the terms in that equation vanish identically at \(r=0\). This seems to be a generic feature of f(R) collapse.

4 Non-separable interior solutions

We will now consider matching of Ricci scalar and its derivatives with a non-separable spherically symmetric metric of the form given in Eq. (30). For convenience, we write \({\mathcal A}(r) = (1-h(r))^{-1}\), and thus we have our ansatz for the interior metric

$$\begin{aligned} ds^{2}=-dt^{2}+\frac{Q^{\prime 2}}{1-h(r)}dr^2+Q^{2}(r,t)d\Omega ^{2}, \end{aligned}$$
(47)

where Q(r, t) is the co-moving radius of the collapsing matter, and h(r) is function of r only. This metric has the form of a general LTB solution. We can calculate the Ricci scalar as

$$\begin{aligned} R(r,t)=\frac{1}{Q^{2}Q^{\prime }}\frac{d}{dr}\left[ 2Q\left( \dot{Q}^{2}+Q\ddot{Q}+h(r)\right) \right] . \end{aligned}$$
(48)

Since we want to match Ricci scalar and it’s derivative across a junction, as the simplest choice, we put

$$\begin{aligned} \dot{Q}^{2}+Q\ddot{Q}=0. \end{aligned}$$
(49)

Then, the Ricci scalar takes the simple form

$$\begin{aligned} R(r,t)=\frac{1}{Q^{2}Q^{\prime }}\frac{d}{dr}\left( 2Qh\right) . \end{aligned}$$
(50)

The solution of Eq. (49) is

$$\begin{aligned} Q(r, t)= r\sqrt{g(r)-2bf(r)\left( t-t_{0}\right) }, \end{aligned}$$
(51)

where \(b >0\) is a constant and f(r) and g(r) are two (positive) function of r, which we have to choose. Without loss of generality, we will henceforth set \(b=1/2\), along with \(t_0=0\), so that our collapse process begins at the origin of the co-moving time.

Also we need to take h(r) in such a way that both Ricci scalar and it’s derivative are continuous across the junction at \(r_{0}\). We will make a simple choice here, and set

$$\begin{aligned} h(r)= & {} (r_{0}-r)^{2},~~g(r) = \left( r_0-r\right) ^{-4},\nonumber \\ f(r)= & {} \left( r_0-r\right) ^{-2}. \end{aligned}$$
(52)

With this choice, from Eq. (50), the Ricci scalar reads,

$$\begin{aligned} R = \frac{2\left( r_0-r\right) ^7\left[ 1-\left( 2r^2 - 3rr_0 + r_0^2\right) t\right] }{r^2\left[ r+r_0 - \left( r-r_0\right) ^2\left( r + 2r_0\right) t + \left( r - r_0\right) ^4r_0t^2\right] }.\nonumber \\ \end{aligned}$$
(53)

It is then seen that continuity of the Ricci scalar and its derivative is guaranteed across the co-moving boundary, which for simplicity we will now choose as \(r_0=1\). Note that at \(t=0\), there is an initial singularity at the origin, and the Ricci scalar diverges as \(R \sim 1/r^2\). We will however concentrate on the singularity that forms due to the collapse process, in which case \(R \sim 1/r^3\) near the origin, at the time of formation of the central singularity. However, we note that the process described in this section may not correspond to the realistic collapse of a dense star, contrary to the analysis of the previous section.

To this end, note that this singularity forms along the curve \(t=t_{s}(r)\) defined by

$$\begin{aligned} Q(t_{s}(r), r)=0 \quad \mathrm{i.e.}~~ t_{s}(r)=\frac{1}{\left( 1-r\right) ^2}, \end{aligned}$$
(54)

and the co-moving time for the formation of the apparent horizon \(t_{ah}(r)\) is given by

$$\begin{aligned} t_{ah}(r)=\frac{8-5r}{4\left( 2-r\right) \left( 1-r\right) ^2} \end{aligned}$$
(55)

This implies that in the reference frame of a co-moving observer (at fixed r), the singularity formation is not simultaneous (note that it was simultaneous in the case of separable solutions), rather it is a curve in the \(t-r\) plane which starts at \((t,r)=(1,0) \). If the apparent horizon starts forming at a co-moving time that is earlier than that of singularity formation, then the event horizon can fully cover the singularity and the end stage is a black hole. On the other hand, if trapped surfaces form after the singularity, then it is possible that a non-space-like geodesic might come out of the singularity to reach an external observer, and in that case the final singularity will be visible, i.e the fate of the collapse will be at least a locally naked singularity.

Fig. 7
figure 7

\(\rho \) vs t

Full size image
Fig. 8
figure 8

Q(r, t) vs t

Full size image
Fig. 9
figure 9

\(\sigma \) vs t

Full size image

In Fig. 7, we show the apparent horizon curve of Eq. (55) as a function of time. This is shown in red, and the dotted blue line is the time of formation of the central singularity, i.e \(t_s(r=0)=1\). Clearly, all co-moving observers will see the formation of the central singularity first, and therefore conclude that the collapse results in a singularity that is locally naked. In Fig. 8, we show the logarithm of Q(r, t) as a function of time. Here, the thick red, dotted blue, dashed black and dot-dashed brown curves correspond to \(r=0.001\), 0.002, 0.01 and 0.02, respectively.

The expansion scalar for a time-like geodesic congruence is calculated to be

$$\begin{aligned} \Theta= & {} \frac{3 (r-1)^4 t-(r-1)^2 (r+3)}{2 \left( (r-1)^4 t^2-(r+2) (r-1)^2 t+r+1\right) },\nonumber \\&\mathrm{i.e}~~ \Theta \big |_{r \rightarrow 0} = -\frac{3}{2(1-t)}, \end{aligned}$$
(56)

showing the central divergence at \(t_s(r=0)=1\). Also, we record the expression for the shear,

$$\begin{aligned} \sigma =r(1-r)^2[1+r-t(1-r)^2(2+r)+ t^2(1-r)^4]^{-1}\nonumber \\ \end{aligned}$$
(57)

In Fig. 9, we show the behaviour of \(\sigma \) as a function of r for \(t=0\) (thick red), 0.1 (dotted blue), 0.5 (dashed black) and 0.9 (dot dashed brown). Clearly, as the collapse progresses in co-moving time, the shear which was initially regular at the center increases near the origin, and diverges as \(\sigma \sim r^{-1}\) at the origin for \(t \rightarrow 1\), as can be seen from Eq. (57).

Fig. 10
figure 10

\(\rho \) vs r

Full size image
Fig. 11
figure 11

\(p_{\theta } - p_r\) vs r

Full size image
Fig. 12
figure 12

q vs r

Full size image
Fig. 13
figure 13

Conditions (i), (ii), (iv) at \(t=0.001\)

Full size image
Fig. 14
figure 14

Conditions (i), (ii), (iv) at \(t=0.5\)

Full size image
Fig. 15
figure 15

Conditions (i), (ii), (iv) at \(t=0.99\)

Full size image
Fig. 16
figure 16

Conditions (v), (vi), (viii) at \(t=0.001\)

Full size image
Fig. 17
figure 17

Conditions (v), (vi), (viii) at \(t=0.5\)

Full size image
Fig. 18
figure 18

Conditions (v), (vi), (viii) at \(t=0.99\)

Full size image

In Figs. 10, 11 and 12, we show the density \(\rho \), the anisotropy \(\Pi = p_{\theta } - p_r\) and the heat flux q as a function of r, for \(t=0.001\) (thick red), 0.5 (dotted blue) and 0.99 (dashed black), respectively. These are computed from Eq. (6) and we have set \(\eta = 10\) and \(\alpha = 10^{-3}\).

It remains to check whether the conditions listed in Eqs. (18) and (20) are satisfied during the collapse process. Without loss of generality, we will make a further choice \(\theta = \pi /2\) here, so that \(\sigma _{22} = \sigma _{33}\) and hence we have to look at the conditions (i), (ii) of Eq. (18) and (iv), (v), (vi) and (viii) of Eq. (20). Snapshots of the logarithms of the relevant quantities on the left hand side of the corresponding equations at \(t=0.001\), \(t=0.05\) and \(t=0.99\) are shown in Figs. 13, 14, 15, 16, 17, 18. In Figs. 13, 14 and 15, the thick red, dotted blue and dashed black lines represent the logarithms of conditions (i), (ii) and (iv) and in Figs. 16, 17 and 18,these represent the logarithms of conditions (v), (vi) and (viii). We find that all the required conditions are indeed satisfied. It is also checked that Eq. (16) is identically satisfied in this case as well. As a remark, we note that the shear vanishes at the center (vide Eq. (57)). Moreover, the anisotropy diverges at the origin (due to the singular nature of the solution at \(t=0\)) at all times during the collapse. However, Eq. (17) is identically satisfied in this case, as can be checked.

The solution discussed above collapses into a locally naked singularity, as we have said. We mention in passing that it is also possible to generate black hole solutions from the generic class of non-separable metrics that we consider here. For example, one simply needs to tune the parameter b in Eq. (51) to \(b=2\) (instead of \(b=1/2\) used in the previous example) to see that close to the center, the apparent horizon forms earlier than the singularity (at \(t=0.25\)). Again, one can check that all the energy conditions can be satisfied by suitably tuning the parameters \(\alpha \) and \(\eta \). However, we will not go into the details here, as these are entirely similar to the situation that we have considered.

5 Nature of the collapsing fluid for separable solutions

The solutions presented in the previous sections indicate collapse in f(R) gravity to black holes or to singularities that are locally naked, while obeying all the energy conditions. A natural question in this context is the physical nature of the fluids, namely if they follow an equation of state (EOS). The lack of this analysis is a drawback in many studies of gravitational collapse in f(R) theories that appear in the literature till date. In this context, note that the EOS of collapsing stars is well studied especially in the non-relativistic limit, starting from the pioneering work of [36]. In the context of f(R) collapse, such a study is somewhat difficult to envisage, but clearly we can see from Eq. (41) that there is a priori no simple EOS that our co-moving observer will see, even in the simple case of the separable solutions presented in Sect. 3. We will concentrate only on this class of solutions in this section, since the solution is Sect. 4 does not correspond to realistic collapse of a dense star, as already mentioned.

First, we note that if we set the f(R) parameter \(\alpha \) to zero, we have here

$$\begin{aligned} \frac{p_r}{\rho }\biggr |_{\alpha \rightarrow 0} = \frac{p_{\theta }}{\rho } \biggr |_{\alpha \rightarrow 0} = 1 - \frac{2}{5-2t}. \end{aligned}$$
(58)

Hence, at \(t=0\), the matter follows a barotropic equation of state with \(p_r = p_{\theta } = \gamma \rho \) (remember that there is no pressure anisotropy with \(\alpha \rightarrow 0\) as we have commented on at the end of Sect. 3), with \(\gamma = 3/5\). As the collapse proceeds, the barotropic index reduces in this case, and approaches 1/3 for \(t \rightarrow 1\). Hence, at the end of the collapse, with \(\alpha \rightarrow 0\), the matter reduces to pure radiation.

In the general case, the situation is more complicated. Here, we have, from Eq. (41),

$$\begin{aligned}&\frac{p_r}{\rho }\biggr |_{t \rightarrow 0} = \frac{3+28\alpha \left( r-1\right) ^4}{5+8\alpha \left( r-1\right) ^2},\nonumber \\&\frac{p_{\theta }}{\rho } \biggr |_{t \rightarrow 0} = \frac{3-4\alpha \left( r-1\right) ^4}{5+8\alpha \left( r-1\right) ^2},~~ \frac{p_r}{\rho }\biggr |_{t \rightarrow 1}=\frac{p_{\theta }}{\rho } \biggr |_{t \rightarrow 1} = 1.\nonumber \\ \end{aligned}$$
(59)

It is therefore seen that for small values of \(\alpha \) (we have used \(\alpha = 10^{-3}\) in Sect. 3), at the beginning of the collapse, the system is close to a barotropic fluid with \(\gamma = 0.6\), but the effect of \(\alpha \) is to increase the barotropic index to unity at the time of formation of the singularity, and at this time the speed of sound equals the speed of light. However, the latter fact is true strictly at the singularity formation time, before which the barotropic index is always less that unity. We mention in passing that a related question is whether one can envisage a situation where the fluid in question consists of two simple fluids, each of which possibly follow an equation of state. This is usually achieved for cases without shear or heat flow, by rotating the coordinate basis of the co-moving observer. This has been a popular topic in the literature, starting from the work of [37] (see also [38] for applications in the cosmological context). It can be checked for our model that this is not possible in the presence of heat flux. The intuitive reason for this is that a dissipative effect cannot be un-done by a rotation of the coordinate basis (unless there is a specific form of an equation of state which also involves the heat flux, see e.g [39]).

It is also of interest to consider the heat transport equation in our non-equilibrium collapsing scenario, following the pioneering work of [40]. The simplicity of the solutions derived in Eq. (41) in the separable metric case, allows for explicit computations of the quantities appearing in the evolution equation of the heat flux, which reads [40] (see also [41, 42])

$$\begin{aligned} \tau n_{\mu }h^{\mu \nu }q_{\nu ;\sigma }u^{\sigma } + q= & {} -\kappa n_{\mu }h^{\mu \nu }\left( T_{,\nu } + T a_{\nu }\right) \nonumber \\&-\frac{1}{2}\kappa T^2 q\left( \frac{\tau u^{\mu }}{\kappa T^2}\right) _{; \mu }. \end{aligned}$$
(60)

Here, T is the local equilibrium temperature, \(\kappa \) is the thermal conductivity, and \(\tau \) the relaxation timescale, and all these quantities must be positive, from physical conditions. Further, \(a_{\mu }\) is the acceleration vector defined after Eq. (8). In order to solve Eq. (60), a number of assumptions is necessary, since \(\kappa \) and \(\tau \) are temperature dependent quantities. There is a vast amount of literature on the topic, and we do not go into the known details here, but will simply use the results of [43, 44] (see also [45]) and assume that

$$\begin{aligned} \kappa = \gamma T^3\tau _c,~~\tau = \left( \frac{\beta \gamma }{\alpha _1}\right) \tau _c,~~ \tau _c = \left( \frac{\alpha _1}{\gamma }\right) T^{-\sigma _1}, \end{aligned}$$
(61)

where \(\beta \), \(\gamma \), \(\alpha _1\) and \(\sigma _1\) are non-negative constants, with the case \(\beta = 0\) being the non-causal case (see, e.g [41] for an excellent exposition). For simplicity, we will restrict ourselves to cases with \(\sigma _1 \le 4\).

Although Eq. (60) is in general difficult to solve, the simplicity of the form of the energy-momentum tensor for the separable solution considered in Sect. 3 allows us to obtain analytic solutions at least in some approximations. First of all, let us consider the non-causal case, and set \(\beta = 0\). Then, we obtain the formal solution of Eq. (60) as

$$\begin{aligned} T^{4-\sigma _1}= & {} -\frac{\left( 1-r\right) ^2\left( \sigma _1 - 4\right) }{\alpha _1\left( 1-t\right) \left( \sigma _1 - 3\right) } - \alpha \frac{6\left( 1-r\right) ^6\left( \sigma _1 - 4\right) }{\alpha _1\left( 1-t\right) ^2\left( \sigma _1 - 1\right) }\nonumber \\&+ \left( 1-r\right) ^{8-2\sigma _1}F\left( t\right) ~~~~~~\left( \beta = 0\right) , \end{aligned}$$
(62)

where F(t) is an a priori undetermined function of the co-moving time. The special cases \(\sigma _1 = 1, 3, 4\) need to be solved separately. The results are

$$\begin{aligned} T^3= & {} -\frac{3\left( 1-r\right) ^2}{2\alpha _1\left( 1-t\right) } + \alpha \frac{36\left( 1-r\right) ^6\log \left( 1-r\right) }{\alpha _1\left( 1-t\right) ^2}\nonumber \\&+ \left( 1-r\right) ^6 F(t)\quad \left( \beta =0, ~\sigma _1 =1\right) , \nonumber \\ T= & {} \left( 1-r\right) ^2\frac{2\log \left( 1-r\right) }{\alpha _1\left( 1-t\right) } + \alpha \frac{3\left( 1-r\right) ^6}{\alpha _1\left( 1-t\right) ^2} + \left( 1-r\right) ^2F(t)\nonumber \\&\left( \beta =0, ~\sigma _1 = 3\right) ,\nonumber \\ T= & {} \mathrm{Exp}\left[ \frac{\left( 1-r\right) ^2\left( 1-t+2\alpha \left( 1-r\right) ^4\right) }{\alpha _1\left( 1-t\right) ^2}\right] \left( 1-r\right) ^2 F(t)\nonumber \\&\left( \beta =0, ~\sigma _1 = 4\right) , \end{aligned}$$
(63)

where we have generically denoted an arbitrary function of the co-moving time by F(t). Equations (62) and (63) are the full set of solutions for the non-causal case, and the role of the f(R) parameter \(\alpha \) can be readily identified, and the increase in the core temperature as a function of time is clearly seen. In particular, we see from Eq. (62) that the role of \(\alpha \) is to decrease the temperature (compared to the \(\alpha = 0\) case) for \(\sigma _1 < 1\) and \(\sigma _1 > 4\), while it increases the temperature for \(1< \sigma _1 < 4\). Also, from the first two equations of Eq. (63), it is clear that close to the center, the effect of \(\alpha \) vanishes for \(\sigma _1=1\) and dominates for \(\sigma _1 = 3\) with the term independent of \(\alpha \) vanishing in the latter case. No further conclusions can be reached without the knowledge of the arbitrary function F(t).

However, we can make the following observation from Eq. (62). Close to the boundary, i.e as \(r \rightarrow 1\), one can always make the last term on the right hand side of this equation arbitrarily close to zero at a given co-moving time of the collapse, for \(\sigma _1<3\). The fall-off of this term with r being faster than the first term on the right hand side of Eq. (62) indicates that in such cases, there will exist a domain of the co-moving radius where the temperature will not be real (since \(\alpha \) is positive). In order to avoid this, we require \(3<\sigma _1<4\) and the other values of \(\sigma _1\) are ruled out in the class of models that we consider. Note also that the solutions for \(\sigma _1 = 3\) and \(\sigma _1=4\) do not suffer from this pathology, and hence our final set of admissible values of \(\sigma _1\) is \(3\le \sigma _1 \le 4\).

Note that in cases where the interior solution is matched with an external Vaidya metric, the arbitrary function F(t) can be determined by relating the temperature at the boundary to the luminosity there, and then equating this with the luminosity as seen by an observer at infinity, via the red-shift factor. This is not possible here, as we have matched with an external Schwarzschild solution, for which the temperature and luminosity at the boundary are automatically zero, as is evident by taking the \(r \rightarrow 1\) limit in the solutions above. F(t) can thus be determined in principle if we specify the behaviour of the core temperature as a function of time, along with the condition \(3 \le \sigma _1\le 4\) discussed above.

Finally, we make some comments about the non-causal case. Here, the analysis becomes cumbersome, and analytic solutions to the heat flow equation of Eq. (60) seem difficult to obtain. As a somewhat crude approximation (used in [43,44,45]), if we ignore the last term on the right hand side of Eq. (60), then with Eq. (61), we obtain as a solution for \(\sigma _1 = 0\),

$$\begin{aligned} T^4= & {} -\frac{4\left( 1-r\right) ^2\left( 1-t\right) + 9\beta \left( 1-r\right) ^4}{3\alpha _1\left( 1-t\right) ^2}\nonumber \\&+ \alpha \frac{24\left( 1-r\right) ^6\left[ t-1+5\beta \left( 1-r\right) ^2\log \left( 1-r\right) \right] }{\alpha _1\left( 1-t\right) ^3}\nonumber \\&+ \left( 1-r\right) ^8F(t), \end{aligned}$$
(64)

where again the arbitrary function of time can be constrained if we assume a time profile of the core temperature.

The analysis in this section was related to the separable solutions that we have used in Sect. 3. For the non-separable solutions of Sect. 4, such analyses become quite tedious, and will not provide much physical insight, as should be evident from the comments made at the beginning of that section.

6 Description in the Einstein scalar theory

It is well known that the action of the metric f(R) gravity

$$\begin{aligned} S=\frac{1}{2}\int d^{4}x \sqrt{-g} [R+f(R)] +S^{(mat)}, \end{aligned}$$
(65)

can be mapped to that of the scalar tensor theory with the scalar field \(\phi =R\) (see, e.g., [4])

$$\begin{aligned} S=\frac{1}{2}\int d^{4}x \sqrt{-g}[\Psi (\phi )R-V(\phi )] +S^{(mat)}, \end{aligned}$$
(66)

with the identifications

$$\begin{aligned} \Psi (\phi )=1+f^{\prime }(\phi ) ,~~ V(\phi )=[\phi f^{\prime }(\phi ) -f(\phi )]. \end{aligned}$$
(67)

Here \(S^{(mat)}\) denotes the matter part of the action. Varying the action of Eq. (66), one can see that the condition \(\phi =R\) is satisfied. The field equations derived from (66) can be shown to be given by

$$\begin{aligned} G_{\mu \nu }=\frac{1}{\Psi }\left[ T_{\mu \nu }^{(mat)}-\frac{1}{2}U(\Psi )g_{\mu \nu }+D_{\mu \nu }\Psi \right] , \end{aligned}$$
(68)

and,

$$\begin{aligned} 3\Box \Psi +2U(\Psi )-\Psi \frac{dU(\Psi )}{d\Psi }=T^{(mat)}, \end{aligned}$$
(69)

where \(U(\Psi )=V[\phi (\Psi )]\) and \(T_{\mu \nu }^{(mat)}, T^{(mat)}\) denotes the matter part of the energy momentum tensor and it’s trace respectively.

Here we are interested in the special case of \(f(R)=\alpha R^{2}\). Then we have \(\Psi (\phi )=1+2\alpha \phi \) and \(V(\phi )=\alpha \phi ^{2}\). In this case the equation for the scalar field (\(\phi \)) obtained from Eq. (69) simplifies to

$$\begin{aligned} 6\alpha \Box \phi -\phi =T^{(mat)}. \end{aligned}$$
(70)

Taking trace of the field Eq. (4) of f(R) gravity, one can check that this is the same as the resulting equation for R, as expected. Now for our collapse model with the separable solution of Eq. (35) along with the ansatz of Eq. (39), we see that the above equation is satisfied with \(\phi (r,t)=R(r,t)\). Thus we have an equivalent scalar-tensor theory description in Jordan frame with the potential \(V(\phi )=\alpha \phi ^{2}\). We also note that the kinetic term of the scalar field is zero.Footnote 7

Now, we turn to the description of the junction conditions in the Jordan frame description of scalar tensor theory, with corresponding ones of f(R) gravity listed in Eq. (24). Recently in [46], the authors have extensively studied the junction conditions of a wide class of scalar-tensor theories (in both the Jordan and Einstein frames) for non-null and null hypersurfaces. We only mention the relevant results here, and refer the reader to [46] for details.

The first junction condition indicates a continuity of the induced metric on the hypersurface and the scalar field

$$\begin{aligned} {[}g_{ab}]=0, ~~~ [\phi ]=0. \end{aligned}$$
(71)

The first and third conditions in Eq. (24) imply that the above conditions are satisfied by the collapsing solution we have that constructed.

The second junction condition gives the discontinuity of the derivative of the scalar field across the hypersurface \(\Sigma \), and the value of the surface energy momentum tensor as,

$$\begin{aligned}&Z=2\Psi ^{\prime }(\phi )[K], ~\mathrm{and}\nonumber \\&S_{ab}=-2\Psi (\phi ) \left( [K_{ab}]-[K]g_{ab} \right) +2Z\Psi ^{\prime }(\phi )g_{ab}, \end{aligned}$$
(72)

in terms of a scalar field \(Z=N^{\mu }[\nabla _{\mu }\phi ]\). Since \(\Psi ^{\prime }(\phi )\ne 0\) at the hypersurface, for a for a smooth matching across \(\Sigma \) the necessary and sufficient conditions are \(Z=0\) and \(S_{ab}=0\). The requirement \(S_{ab}=0\) can be satisfied with only one the following three cases listed below (see proposition 2 of [46])Footnote 8

  1. 1.

    \([K_{ab}]=Z=0\),

  2. 2.

    \(\Psi ^{\prime }(\phi )_{\Sigma }=\Psi (\phi )_{\Sigma }=Z=0\),  or

  3. 3.

    \(6\Psi ^{\prime }(\phi )_{\Sigma }^{2} + 2\Psi (\phi )_{\Sigma }=0, ~~ 3[K_{ab}]-[K]h_{ab}=0\).

We can prove these relations by taking the trace of \(S_{ab}=0\) at \(\Sigma \). We thus arrive at two conditions: (a) \([K]=0\), or (b) \(6\Psi ^{\prime }(\phi )_{\Sigma }^{2} + 2\Psi (\phi )_{\Sigma }=0\). Now, \([K]=0\) in turn implies \(Z=0\) and this in conjunction with \(S_{ab}=0\) implies \([K_{ab}]=0\). This is the case 1 above. On the other hand when (b) is satisfied, the second junction condition gives cases (2) and (3). Note that for cases 1 and 2, the condition \(Z=0\) is satisfied whenever \(S_{ab}=0\), and hence both criteria of smooth matching is fulfilled. Though for case 3, we have to make sure Z vanishes along with \(S_{ab}\).

Let us now compare these conditions with the ones in GR and the extra conditions of f(R) theories for our collapsing solutions. Here, for the \(f(R)=\alpha R^{2}\) model, \(\Psi ^{\prime }(\phi )\) is a constant (\(= 2\alpha \)), thus \([K_{ab}]=0\) automatically implies \(Z=0\), and hence the condition of case 1 is satisfied. However by the same token, \(\Psi ^{\prime }(\phi )_{\Sigma }=\Psi (\phi )_{\Sigma }\ne 0\), and hence the conditions for cases 2 and 3 are not satisfied. Thus as mentioned above, only one of the three given conditions is satisfied. Nevertheless as one can see, the two conditions 2 and 3 are somewhat “special” in the Jordan frame description of scalar tensor theory. Case 1, along with the first junction condition is just one would expect from Eq. (24), in terms of the scalar field. But in case 2 for example, we have no condition on the extrinsic curvature \(K_{\mu \nu }\) [46]. This is a big advantage in the sense that one can obtain regular solutions of the f(R) gravity only by looking at the Ricci scalar.

Since in this paper we have essentially obtained our collapsing solutions by demanding that an interior solution of f(R) gravity can be matched with the outside Vaidya metric through a timelike hypersurface, the cases 2 and 3 can give rise to many interesting solutions. For example using the mapping in Eq. (67) we can readily see that any solution of \(R+\alpha R^{n}, ~~ \alpha >0\) gravity matched with outside Vaidya solution, cannot be made to satisfy the condition of case 3. The simplest such theory can be pure cubic (\(R^{3}\)) gravity with \(R_{\Sigma }= (N^{\mu }\nabla _{\mu }R)_{\Sigma }=0\). Construction of such new solutions is left for a future work.

7 Discussions

Gravitational collapse in metric f(R) theories of gravity are greatly restricted due to the extra junction conditions that have to be invoked, and involve the continuity of the Ricci scalar and its derivatives across a time-like hypersurface on which an internal collapsing metric is matched with an external solution. This is the R-matching method commonly used in f(R) scenarios. In fact, there are a total of twelve conditions that will generically need to be satisfied, and are given in Eqs. (24), (18) and (20). It is indeed a formidable task to compute explicit collapsing solutions while satisfying all these conditions and not many exact solutions are available in the literature.

In this paper, we have constructed novel examples of f(R) collapse, by using the extra junction conditions to explicitly solve for the collapsing metric. This was done in \(f(R) = \alpha R^2\) theories of gravity in two cases, one in which we assumed a separable form of the internal metric, and the other in which this assumption was relaxed. We showed that by suitably choosing some reasonable forms of a few arbitrary functions, new examples of collapse in modified gravity can be constructed. In the latter context, we have described a collapse situation that includes the effects of shear viscosity. The generic relation between the shear viscosity and the anisotropy in our f(R) models has been derived here. We have demonstrated by explicit examples the formation of black holes or naked singularities, while satisfying all the energy conditions.

The separable solution constructed by us allows for analytical forms of the components of the energy momentum tensor. Using these, we are also able to obtain analytical solutions to the evolution equation of the heat flux. Here, the simplicity of the expressions involved allows us to focus on the effect of the f(R) parameter \(\alpha \), with certain reasonable approximations. As mentioned in the text, this situation corresponds to a realistic collapse of a dense star, with the pressures remaining positive at all times. It would be interesting to understand the evolution of the entropy in this non-equilibrium situation.

Our analysis in this paper relies on a number of explicit choices that we have made, and these have been highlighted in Sects. 3 and 4. Indeed, these choices are arbitrary and serve as examples of more general cases than what we study here, and functions different from what we have chosen should generate more physical examples of collapse scenarios in modified gravity. Further, our analysis here is limited to models with \(f(R) = R^2\). It should be interesting to apply this to more generic situations. Finally, it will be very interesting to study the new solutions arising from the junction conditions discussed in Sect. 6. We hope to report on this in a future work.