1 Introduction

The general theory of relativity (GR) is an impeccably robust relativistic theory of gravity. It forms the basis for our understanding of gravitational phenomenon at small as well as large scales [1, 2]. However, it is well known that GR cannot be the ultimate theory of gravitation since it has a well defined regime of validity; for example, understanding past and future spacetime singularities are beyond the reach of GR. Suitable modification(s) of GR is(are) essential to understand singularities, or to resolve them. It is believed that a quantum theory of gravity may lead to solution to the problems affecting GR [2,3,4,5,6]. Naturally, in absence of any consensus on the theory of quantum gravity, modified theories of gravity with quantum corrections are also of interest. The Einstein- Hilbert action may be thought of as only a low energy contribution and higher curvature terms consistent with the diffeomorphism invariance may become relevant as one goes to higher energies. Higher curvature corrections should leave imprints at low energy scales which become important for low energy physics too [5,6,7]. Out of these alternate theories, we shall study the f(R) model since it has been found interesting in the cosmological studies as well. These f(R) gravity models are thought of as alternate to the dark energy models [8,9,10]. The standard way to construct these f(R) theories is to replace the Einstein- Hilbert Lagrangian by a well defined function of Ricci scalar, f(R) (for general relativity \(f(R)=R\)) [11, 12]. For a detailed review of the motivation, validity of various functional forms f(R), applications as well as shortcomings of these gravity theories have been extensively analyzed [13,14,15,16,17,18,19,20,21].

The purpose of the present work is to construct some examples of spacetimes in f(R) gravity which admit gravitationally weak spacetime singularities [22]. The physical situation which we consider here is the following: the initial mass of the collapsing star or matter configuration is so high that the force of gravity overwhelms the thermal or quantum mechanical pressures. As a result, no stable configuration such as a neutron star or a white dwarf exists during the collapse process. Generically, such collapsing matter configurations admit diverging density and curvatures at the central singularity (we must emphasize that spherically symmetry is assumed throughout the collapse). However, we shall show below that it is possible, without violating energy conditions, that during the gravitational collapse matter is radiated away at such a rate that the matter boundary never reaches its horizon. In this case, no horizon is formed at the boundary since the star radiates off all its mass before reaching the singularity. The central singularity is naked but is gravitationally weak. During this process, all the physical quantities energy density, radial and tangential pressure, pressure anisotropy, heat flux, remain regular and positive throughout the collapse. The luminosity and adiabatic index are also regular and positive, and admits maximum value when the matter approaches the singularity. Thus, for an observer at infinity observing the collapse, the configuration shall become extremely bright, reaching its maximum luminosity before turning off, indicating that it has radiated off all its mass. Solutions of such kind are not unknown and is possible in the Newtonian gravity as well. Consider a star in the Newtonian gravity which is extremely heavy to be supported by the Pauli exclusion principle alone. So, when the gravitation contraction takes place, thermal pressure may balance to some extent. But since there is no event horizon, the star shall continue to radiate all the gravitational energy to infinity. Hence all the matter contained in the star shall be converted to thermal radiation and radiated off. We show here that configurations of similar nature are also possible in f(R) gravity.

From a theoretical point of view, gravitational collapse of matter is important and recently, there has been an increasing interest to understand whether the nature of collapse is altered in modified gravity [19, 23,24,25,26]. Although it remains to study in detail if compact objects formed in theories like the f(R) gravity can be experimentally detected, several investigations have already been carried out in this direction. The idea is to use the multi-wavelength and well as the gravitational wave data to probe matter under extreme gravity (like the composition of the inner core of a neutron star) which remain unknown, as neither they can they be probed directly with astrophysical observations nor are they they within the reach of the present-day theory experiments. Alternatively, the direct gravitational wave data have also made it possible to investigate if these objects are possible in alternate gravity theories. In particular, the recent detection of GW190814, have fueled the speculation that it might be the heaviest neutron star till date [27]. It has been argued that some models of f(R) gravity can indeed explain the origin of such a large mass of (2.5–2.7)\(M_{\odot }\), even using the well- known equations of state [28, 29].

In this paper, we shall not concern ourselves with compact objects like the white dwarf or a neutron star. We shall assume that the matter is extremely massive, and such objects continue to collapse under its own gravity. The gravitational collapse phenomenon in GR show that the collapse outcome depends upon, among other quantities, the choices of mass profiles and velocity profiles of the collapsing matter. In the context of inhomogeneous LTB models in GR, these issues have been considered in great detail for various matter models including dust and viscous fluids [30, 31]. The stellar collapse of stellar collapse in f(R) gravity using different forms of f(R) function and different matter distributions may be found in  [32,33,34,35,36,37,38,39,40,41,42,43,44,45,46]. Although different type of f(R) models may be considered, only the ones which are in agreement with the standard cosmological observations should be of interest. Here, we consider three models of f(R) gravity given respectively by: \((a)~f(R)=(R+\lambda R^{2})\) [47]  ,  \((b) ~f(R)\sim R^{1+\epsilon }\), with \(\epsilon<<1\) [37] and \((c)~f(R)=R+\lambda \,\left[ \exp {(-\sigma R)}-1\right] \), with \(\lambda \sigma <1\),   [48]. All three theories play a major role in the study of the early and the late universe. For example, in the \(R+\lambda R^2\) model, the \(R^{2}\) term is weak in the weak gravity regime. But it also contributes significantly during the early universe and in strong gravity regime. Here, we study gravitational collapse of an extremely massive object (so massive that its gravitational attraction predominates thermal or quantum mechanical pressures), and so naturally, since we are in a strong gravitational field, this model may assume significance. In all these theories, we shall show that an extremely massive spherical matter cloud, admitting radial and tangential pressures, and outgoing heat flux, can collapse in a manner that no horizon is formed at the boundary (since the star radiates off all its mass before reaching the singularity) making the central singularity naked but weakly singular. To ensure that the gravitational collapse proceeds continuously without forming any stable object, we assume the radius of the 2-sphere cross-sections decrease linearly with time. The interior collapsing spacetime shall be smoothly matched with the exterior Vaidya spacetime  [49] over a timelike surface \(\Sigma \) [50].

We approach the problem as follows. In Sect. 2, we give the field equations of f(R) gravity and the junction conditions for smooth matching of the interior and the exterior spacetimes across the timelike hypersurface \(\Sigma \). This section also includes the solutions of the f(R) field equations along with the explicit expressions for physical quantities. For the solution, we use the Karmarkar condition [51,52,53]. These conditions determine gravitational potentials for static and non-static systems [54,55,56,57]. We must mention here that similar studies on f(R) gravity have been carried out in [41]. However, the solutions obtained there are restrictive in the sense that one of the metric function have been kept constant to derive the values of other metric function. On the other hand we shall show that such conditions are overly restrictive. Note that since Karmarkar condition expresses relationship between metric functions, the forms of metric functions are arbitrary, and dependent on one’s choice. However we argue that this arbitrariness may be removed if metric functions are related to matter variables. For example, if we assume a specific form of pressure anisotropy (difference of the radial and tangential pressures, denoted by \(\Delta =p_{t}-p_{r}\)), this gives rise to unambiguous set of gravitational potentials, and a class of metric representing collapsing spacetimes. We inspect the physical relevance of these exact solutions by verifying the energy conditions. The stability criteria and discussion about the luminosity and adiabatic index, radial and transverse velocity are carried out in Sect. 2.1. It is shown that a faraway observer will see a source whose luminosity is exponentially increasing until a time when it shuts off quickly. This is due to the fact that the total mass of the star radiates linearly and, as the star reaches its maximum luminosity there is no mass left to radiate. The evolution of the temperature profiles during stellar collapse is also studied in Sect. 2.2 since they play an important role in the study of transport processes in radiative gravitational collapse [58,59,60,61,62,63,64,65,66]. In Sect. 3 contains discussion of the results accompanied with concluding remarks.

2 Field equations and matching conditions

The action for the f(R) gravity is obtained by replacing the standard Einstein–Hilbert Lagrangian by a well defined function of Ricci scalar [16]

$$\begin{aligned} {\mathcal {S}}= & {} \frac{1}{2}\,\int \sqrt{-g}\,\left[ \,f(R)+2\, {\mathcal {L}}_{M}(g_{\mu \nu },\Psi _{m})\right] d^{4}x\, , \end{aligned}$$
(1)

where \(\Psi _{m}\) refers collectively to all matter fields, \( {\mathcal {L}}_{M}\) is the Lagrangian density of the matter fields \(\Psi _{m}\), g is the determinant of the metric tensor \(g_{\mu \nu }\), R is the Ricci scalar curvature and f(R) is the generic function of Ricci scalar defining the theory under consideration and (using units with \(c=1=8\pi G)\). Varying the action (1) with respect to the metric tensor \(g_{\mu \nu }\) yields the following field equations:

$$\begin{aligned} F(R)\,R_{\mu \nu }{-}\frac{1}{2}f(R)\,g_{\mu \nu }{-}\left( \nabla _{\mu }\,\nabla _{\nu }{-}g_{\mu \nu }\,\Box \right) F(R) = T^{M}_{\mu \nu },\nonumber \\ \end{aligned}$$
(2)

where \(F(R)={d\,f(R)}/{dR}\),  and \(\Box \equiv \nabla _{\mu } \nabla ^{\mu }\). This equation may also be rewritten as

$$\begin{aligned} R_{\mu \nu }-(1/2)\, g_{\mu \nu }\,R=F(R)^{-1}\,\left( T^{M}_{\mu \nu }+T^{D}_{\mu \nu }\right) \, , \end{aligned}$$
(3)

where the left side of the Eq. (3) is the usual Einstein tensor, \(T^{M}_{\mu \nu }\) and \(T^{D}_{\mu \nu }\) are the energy momentum tensor and effective energy momentum tensor having the form as:

$$\begin{aligned} T^{M}_{\mu \nu }= & {} (p_{t}+\rho )u_\mu u_\nu +p_{t} g_{\mu \nu }+(p_{r}-p_{t})X_{\mu }X_{\nu }\nonumber \\&+q_{\mu }u_{\nu }+q_{\nu }u_{\mu }\,, \end{aligned}$$
(4)
$$\begin{aligned} T^{D}_{\mu \nu }= & {} (1/2)\left[ f(R)-R\,F(R)\right] \, g_{\mu \nu }\nonumber \\&+\left( \nabla _{\mu }\nabla _{\nu }-g_{\mu \nu }\,\Box \right) F(R). \end{aligned}$$
(5)

Here, \(\rho \), \(p_{r}\) and \(p_{t}\) are the energy density, radial pressure and the tangential pressure respectively. Also, \(q^{\mu }\) \(u^{\mu }\), \(X^{\mu }\) represents the radial heat flow vector, 4-velocity vector and spacelike 4-vector respectively, which satisfy \(u_\mu u^\mu =-X_\mu X^\mu =-1\) and \(u_\mu X^\mu =u_\mu q^\mu =0\).

We now consider a general non- static shear free spherically symmetric spacetime metric given by the following form

$$\begin{aligned} ds^{2}= & {} -a(r)^2 dt^2 + b(r)^2 s(t)^2 \nonumber \\&\quad \times \left( dr^2 + r^2 d \theta ^2+ r^{2} \sin ^2{\theta } d\phi ^2 \right) . \end{aligned}$$
(6)

The forms of \(u^{\mu }\), \(X^{\mu }\) and \(q^{\mu }\) in terms of the metric (6) are

$$\begin{aligned} u^\mu = a^{-1}\,\delta ^{\mu }_{0}; X^{\mu }=(b\,s)^{-1}\,\delta ^{\mu }_{1}; q^{\mu } = (b\,s)^{-1}\,X^{\mu } , \end{aligned}$$
(7)

The magnitude of the expansion scalar \(\Theta \) and Ricci scalar for the metric (6) have the form

$$\begin{aligned} \Theta= & {} \bigtriangledown _{\mu }u^{\nu }=\frac{3\,{\dot{s}}}{a\,s}, \end{aligned}$$
(8)
$$\begin{aligned} R= & {} 6\,\frac{s\ddot{s}+{\dot{s}}^{2}}{a^2\,s^2}-\frac{2}{b^2\,s^2}\nonumber \\&\times \left[ \frac{a^{\prime \prime }}{a}- \frac{b^{\prime 2}}{b^{2}}+\frac{a^{\prime }b^{\prime }}{a\,b} +2\frac{b^{\prime \prime }}{b}+\frac{2}{r}\left( \frac{a^{\prime }}{a}+2\frac{b^{\prime }}{b} \right) \right] . \end{aligned}$$
(9)

The field equations in f(R) gravity for the metric (6), energy momentum tensor (4), (5) and (7) are

$$\begin{aligned} \rho= & {} \frac{F}{s^2}\left[ \frac{3\,{\dot{s}}^{2}}{a^2}-\frac{1}{b ^{2}}\left( \frac{2\,b^{\prime \prime }}{b} -\frac{b^{\prime }\,^{2}}{b\,^{2}}+\frac{4}{r}\frac{b^{\prime }}{b} \right) \right] \nonumber \\&+\left( \frac{f-R\,F}{2}\right) +\frac{3{\dot{s}}{\dot{F}}}{s\,a^{2}}\nonumber \\&-\frac{1}{b^{2}\,s^{2}}\left[ F^{\prime \prime } +F^{\prime }\left( \frac{b^{\prime }}{b}+\frac{2}{r}\right) \right] , \end{aligned}$$
(10)
$$\begin{aligned} {p}_{r}= & {} \frac{F}{s^2}\left[ -\frac{1}{a^{2}}\left( 2s\,\ddot{s}+{\dot{s}}^{2}\right) \right. \nonumber \\&\left. +\frac{1}{b^{2}}\left( \frac{2\,a^{\prime }\,b^{\prime }}{a\,b}+\frac{2}{r}\left( \frac{a^{\prime }}{a} +\frac{b^{\prime }}{b}\right) +\frac{b^{\prime }\,^{2}}{b\,^{2}} \right) \right] \nonumber \\&-\left( \frac{f-R\,F}{2}\right) \nonumber \\&-\frac{{\dot{F}}}{a^2}\left( \frac{{\ddot{F}}}{{\dot{F}}}+\frac{2{\dot{s}}}{s}\right) +\frac{F^{\prime }}{b^{2}\,s^{2}}\left( \frac{a^{\prime }}{a}+\frac{2}{r}+\frac{2b^{\prime }}{b}\right) , \end{aligned}$$
(11)
$$\begin{aligned} {p}_{t}= & {} \frac{F}{s^2}\left[ -\frac{1}{a^{2}}\left( 2s\,\ddot{s}+{\dot{s}}^{2}\right) \right. \nonumber \\&\left. +\frac{1}{b^{2}}\left( \frac{a^{\prime \prime }}{a}+\frac{b^{\prime \prime }}{b}-\frac{b^{\prime }\,^{2}}{b\,^{2}} +\frac{1}{r}\left( \frac{a^{\prime }}{a}+\frac{b^{\prime }}{b}\right) \right) \right] \nonumber \\&-\left( \frac{f-R\,F}{2}\right) \nonumber \\&-\frac{{\dot{F}}}{a^2}\left( \frac{{\ddot{F}}}{{\dot{F}}} +\frac{2{\dot{s}}}{s}\right) +\frac{1}{b^{2}\,s^{2}}\left( F^{\prime \prime } +F^{\prime }\left( \frac{a^{\prime }}{a}+\frac{1}{r}\right) \right) ,\nonumber \\ \end{aligned}$$
(12)
$$\begin{aligned} q= & {} -\frac{2\,a^{\prime }\,{\dot{s}}\,{\dot{F}}}{a^{2}\,b^{2}\,s^{3}} +\frac{1}{a^{2}\,b^{2}\,s^{2}} \left( {\dot{F}}\,^{\prime }-\frac{{\dot{F}} a^{\prime }}{a} -\frac{{\dot{s}} \, F^{\prime }}{s} \right) , \end{aligned}$$
(13)

where prime and dot are the derivatives with respect to r and t respectively.

Let us consider the junction conditions for the smooth matching of the interior manifold \(\mathcal {M^-}\) (Eq. (6) considered above) with the exterior manifold \(\mathcal {M^+}\) across timelike hypersurface \(\Sigma \), at \(r=r_b\). As described in [67, 68], the junction conditions for the f(R) gravity requires the matching of several geometric quantities other than the induced metric (\(h_{ij}\)) and the extrinsic curvature (\(K_{ij})\). In fact, it has been established that in f(R) gravity, the following variables must be matched at the boundary:

$$\begin{aligned} {[}h_{ij}]_{-}^{+}= & {} 0, \end{aligned}$$
(14)
$$\begin{aligned} F(R)\,[K_{ij}-(1/3)K\, h_{ij}]_{-}^{+}= & {} 0, \end{aligned}$$
(15)
$$\begin{aligned} ~[K]_{-}^{+}= & {} 0, \end{aligned}$$
(16)
$$\begin{aligned} (\partial F(R)/ \partial R)\,[\partial _{\tau }R]_{-}^{+}= & {} 0, \end{aligned}$$
(17)
$$\begin{aligned} ~[R]_{-}^{+}= & {} 0, \end{aligned}$$
(18)

where \(\tau \) represents the proper time of the timelike hypersurface, and K is the trace of the extrinsic curvature. Out of these five conditions, for the set of f(R) theories under consideration, it is sufficient to match the metric, the extrinsic curvature, the Ricci scalar, and the derivative of the Ricci as given above, determined from either sides.

The Vaidya spacetime in the outgoing coordinate is taken to be our exterior spacetime \(\mathcal {M^+}\) [49]

$$\begin{aligned} ds^{2}_{+}= & {} -\left[ 1-\frac{2\,M(v)}{{\mathbf{r}}}\right] \, dv^{2}-2dvd{\mathbf{r}}\nonumber \\&+ {\mathbf{r}}^2 \left( d \theta ^2 +\sin ^2\theta d\phi ^2\right) , \end{aligned}$$
(19)

For our later convenience, let us define the proper time, \(d\tau =a(r)_{_{\Sigma }}\, dt \). The junction condition as given (14), implies the following conditions

$$\begin{aligned} {\mathbf{r}}_{_{\Sigma }}(v)= & {} (r\,b\,s)_{_{\Sigma }} , \end{aligned}$$
(20)
$$\begin{aligned} \left( \frac{dv}{d\tau }\right) _{\Sigma }^{-2}= & {} \left( 1-\frac{2M}{\mathbf{r}}+2\frac{d\mathbf{r}}{dv}\right) _{\Sigma }, \end{aligned}$$
(21)

where \(\tau \) represents the proper time defined on the hypersurface \(\Sigma \). The normal vector fields to \(\Sigma \) are given by

$$\begin{aligned} n^{-}_{l}= & {} \left[ 0,(b\,s)_{_{\Sigma }},0,0\right] ,\nonumber \\ n^{+}_{l}= & {} \left[ 1-\frac{2M}{{\mathbf{r}}}+2\frac{d{\mathbf{r}}}{dv}\right] ^{-\frac{1}{2}}_\Sigma \left[ -\frac{d{\mathbf{r}}}{dv}\delta ^0_l+\delta ^1_l\right] _{\Sigma } . \end{aligned}$$
(22)

The extrinsic curvatures for metrics (6) and (19) are given by

$$\begin{aligned} K^{-}_{\tau \tau }= & {} -\left[ \frac{a^{\prime }}{a\,b\,s}\right] _{\Sigma }, K^{-}_{\theta \theta }=\left[ r\,b\,s\left( 1+\frac{r\,b^{\prime }}{b} \right) \right] _\Sigma , \end{aligned}$$
(23)
$$\begin{aligned} K^{+}_{\tau \tau }= & {} \left[ \frac{d^2v}{d\tau ^2}\left( \frac{dv}{d\tau }\right) ^{-1}-\left( \frac{dv}{d\tau }\right) \frac{M}{{\mathbf{r}}^2}\right] _{\Sigma },\nonumber \\ K^{+}_{\theta \theta }= & {} \left[ \left( \frac{dv}{d\tau }\right) \left( 1-\frac{2M}{{\mathbf{r}}}\right) {\mathbf{r}}-{\mathbf{r}}\frac{d{\mathbf{r}}}{d\tau }\right] _\Sigma , \end{aligned}$$
(24)
$$\begin{aligned} K^{-}_{\phi \phi }= & {} \sin ^2{\theta } K^{-}_{\theta \theta }\,\, , \,\, K^{+}_{\phi \phi }=\sin ^2{\theta } K^{+}_{\theta \theta } . \end{aligned}$$
(25)

Now, from the junction condition on \(K_{ij}\) (because of the conditions that K must satisfy on the hypersurface, matching \(K_{ij}\) is enough), we get the following. From the equality for the \(\theta \theta \) components at hypersurface \(\Sigma \), and the Eqs. (20) and (21), we obtain

$$\begin{aligned}&\left[ r\,b\,s\left( 1+\frac{r\,b^{\prime }}{b} \right) \right] _\Sigma \nonumber \\&\quad =\left[ \left( \frac{dv}{d\tau }\right) \left( 1-\frac{2M}{{\mathbf{r}}}\right) {\mathbf{r}}-{\mathbf{r}}\frac{d{\mathbf{r}}}{d\tau }\right] _\Sigma , \end{aligned}$$
(26)

and the total energy inside the boundary hypersurface \(\Sigma \), given by the Misner–Sharp mass, denoted by 2m (such that \(m=M\) on the matching hypersurface) [69, 70], where

$$\begin{aligned} m_{\Sigma }= & {} \left[ \frac{r^{3}\,{\dot{s}}^{2}\,b^{3}s}{2\,a^{2}}-\frac{r^{3}\,s\,b^{\prime }\,^{2}}{2\,b} -r^2\,s\,b^{\prime }\right] _{\Sigma }. \end{aligned}$$
(27)

Now, again from the matching of the \(K^{+}_{\tau \tau }=K^{-}_{\tau \tau }\) component we have the following equation:

$$\begin{aligned} -\left[ \frac{a^{\prime }}{a\,b\,s}\right] _{\Sigma }= & {} \left[ \frac{d^2v}{d\tau ^2}\left( \frac{dv}{d\tau }\right) ^{-1}-\left( \frac{dv}{d\tau }\right) \frac{M}{{\mathbf{r}}^2}\right] _\Sigma , \end{aligned}$$
(28)

and, substituting the relation between proper and coordinate time along with the Eqs. (20) and (27) into the Eq. (26) we have

$$\begin{aligned} \left( \frac{dv}{d\tau }\right) _{\Sigma }= & {} \left( 1+\frac{r\,b^{\prime }}{b}+\frac{r\,b\,{\dot{s}}}{a}\right) _{\Sigma }^{-1}. \end{aligned}$$
(29)

Now, differentiating (29) with respect to the \(\tau \) and using Eqs. (27) and (29), we can rewrite (28). Further, comparing with equations (11) and (13) we have the following useful form

$$\begin{aligned} \left( p_{r}+T^{D}_{rr}+b\,sT^{D}_{tr}\right) _{_{\Sigma }}= & {} (q\,b\,s)_{_{\Sigma }}. \end{aligned}$$
(30)

where,

$$\begin{aligned} T^{D}_{rr}= & {} \left( \frac{f-R\,F}{2}\right) +\frac{{\dot{F}}}{a^2}\left( \frac{{\ddot{F}}}{{\dot{F}}}+\frac{2{\dot{s}}}{s}\right) \nonumber \\&-\frac{F^{\prime }}{b^{2}\,s^{2}}\left( \frac{a^{\prime }}{a}+\frac{2}{r}+\frac{2b^{\prime }}{b} \right) , \end{aligned}$$
(31)
$$\begin{aligned} T^{D}_{tr}= & {} \frac{1}{a^{2}\,b^{2}\,s^{2}} \left( {\dot{F}}\,^{\prime }-\frac{{\dot{F}}a^{\prime }}{a}-\frac{{\dot{s}}\,F^{\prime }}{s} \right) , \end{aligned}$$
(32)

are the dark source terms. From Eq. (30), it is found that just like for general relativity, the radial pressure does not vanish at the boundary but, instead is proportional to the dissipative as well as radiative dark source terms. The extra terms \(T^{D}_{rr}\) and \(T^{D}_{tr}\) on the LHS of Eq. (30) are the dark source term and may appear due to the higher order curvature geometry of the collapsing sphere [41].

Let us now move to match K as given in (14). The expressions for the trace of extrinsic curvatures on the either sides lead to the following matching condition on the hypersurface:

$$\begin{aligned} \left[ p_{r}+T^{D}_{rr}+b\,sT^{D}_{tr} -q\,b\,s\right] _{_{\Sigma }} = 2\,[M - m]_{_{\Sigma }}, \end{aligned}$$
(33)

and naturally, this condition is identically satisfied due to the above mentioned equations. The matching of the Ricci and its proper time derivative gives the following conditions which are to be satisfied for the metric of the internal manifold (at the hypersurface \(\Sigma \)):

$$\begin{aligned} s\ddot{s}+{\dot{s}}^{2}= & {} \frac{a^2}{3b^2}\left[ \frac{a^{\prime \prime }}{a}+2\frac{b^{\prime \prime }}{b} +\frac{a^{\prime }b^{\prime }}{a\,b}\right. \nonumber \\&\left. +\frac{2}{r}\left( \frac{a^{\prime }}{a}+2\frac{b^{\prime }}{b} \right) - \frac{b^{\prime 2}}{b^{2}}\right] ,\nonumber \\ 3{\dot{s}}\ddot{s}+s\dddot{s}= & {} 0. \end{aligned}$$
(34)

The metric of the internal manifold must be chosen so as to satisfy the two conditions in (34). To determine metric functions according to all these junction conditions, it is necessary to use some auxiliary conditions. We shall see below that these equations are consistent with a collapsing time dependent internal metric. In fact, one may argue that junction conditions indeed force such a possibility. Additionally, we must also ascertain the physical viability of the spacetime metric. From Eqs. (11), (12) and (6), the pressure anisotropy factor \(\Delta =p_{t}-p_{r}\) has the form

$$\begin{aligned} \Delta= & {} \frac{F}{b^2\,s^2}\left[ \frac{b^{\prime \prime }}{b}-\frac{2\,b^{\prime \,^{2}}}{b^2}+\frac{a^{\prime \prime }}{a} -\frac{2\,a^{\prime }\,b^{\prime }}{a\,b}-\frac{1}{r}\left( \frac{a^{\prime }}{a}+\frac{b^{\prime }}{b}\right) \right] \nonumber \\&+\frac{F^{\prime \prime }}{b^2\,s^2}-\frac{F^{\prime }}{b^2\,s^2}\left( \frac{2b^{\prime }}{b}+\frac{1}{r} \right) . \end{aligned}$$
(35)

The general expression for the shear free spacetime as given in (35) is has the complicated form. To find the solution of the metric functions and mathematical simplicity, we take an adhoc form of the pressure anisotropy \(\Delta \) to be:

$$\begin{aligned} \Delta= & {} \frac{F^{\prime \prime }}{b^2\,s^2}-\frac{F^{\prime }}{b^2\,s^2}\left( \frac{2b^{\prime }}{b} +\frac{1}{r}\right) \nonumber \\&-\frac{F}{b^2\,s^2}\left[ \frac{2\,a^{\prime }\,b^{\prime }}{a\,b} -\frac{a^{\prime \prime }}{a}+\frac{a^{\prime }}{r\,a}\right] \end{aligned}$$
(36)

Although, we have chosen this form of the anisotropy in pressure \(\Delta \) for the mathematical simplicity, later we will see that they represents the physically viable solutions of the potentials. Also, this choice of \(\Delta \) is physically significant, such that \(\Delta \) is regular throughout the collapse. It must be noted that this choice of the anisotropy (36) reduces the total pressure anisotropic equation (35) as differential equation of only one function, given by

$$\begin{aligned} 0= & {} \frac{1}{s^2\,b^2}\left( \frac{b^{\prime \prime }}{b}-\frac{2\,b^{\prime \,^{2}}}{b^2}-\frac{b^{\prime }}{r\,b} \right) , \end{aligned}$$
(37)

The form of the function b(r) is

$$\begin{aligned} b(r)= & {} -2[C_3\,r^2+2\,C_4]^{-1}, \end{aligned}$$
(38)

where \(C_3\) and \(C_4\) are constant of integration.

Let us now use the fact that under certain conditions, a \((n+1)\)-dimensional space can be embedded into a pseudo Euclidean space of dimension \((n+2)\) [51]. Thus the necessary and sufficient condition for any Riemannian space to be an embedding class \(\mathbf {I}\) is the Karmarkar condition [52, 53],

$$\begin{aligned} R_{rtrt}\,R_{\theta \phi \theta \phi }=R_{r\theta r\theta }\,R_{t\phi t\phi } -R_{\theta rt\theta }\,R_{\phi rt\phi }. \end{aligned}$$
(39)

The non vanishing components of the Riemann tensor for the metric (6) are

$$\begin{aligned} R_{rtrt}= & {} a^{2}\left( \frac{a^{\prime \prime }}{a}-\frac{b^{2}s}{a^{2}}\ddot{s}-\frac{a^{\prime }}{a}\frac{b^{\prime }}{b} \right) \,,\nonumber \\ R_{\theta \phi \theta \phi }= & {} {r^{4}b^{2}s^{2}}\left( \frac{b^{2}}{a^{2}}{\dot{s}}^{2}-\frac{2b^{\prime }}{rb} -\frac{b^{\prime \,^{2}}}{b^{2}} \right) \sin ^2\theta , \end{aligned}$$
(40)
$$\begin{aligned} R_{r\theta r\theta }= & {} {r^2b^{2}s^{2}}\left( \frac{b^{2}}{a^{2}}{\dot{s}}^{2}-\frac{b^{\prime }}{rb} -\frac{b^{\prime \prime }}{b}+\frac{b^{\prime \,^{2}}}{b^{2}} \right) \,,\nonumber \\ R_{\theta rt\theta }= & {} \frac{r^{2}b^{2}s}{a}a^{\prime }{\dot{s}}, \end{aligned}$$
(41)
$$\begin{aligned} R_{t\phi t \phi }= & {} {r^2a^2b}\left( \frac{a^{\prime }}{ra}-\frac{b^2s}{a^2}\ddot{s} +\frac{a^{\prime }}{a}\frac{b^{\prime }}{b} \right) \sin ^2\theta \,,\nonumber \\ R_{\phi rt\phi }= & {} \sin ^2\theta R_{\theta rt\theta } . \end{aligned}$$
(42)

Using the expressions for Riemannian tensors from Eqs. (40)–(42) into the Eq. (39) we have

$$\begin{aligned} 0= & {} b^{2}{{\dot{s}}}^2b^{3}\left( \frac{a^{\prime \prime }}{a}-2\frac{a^{\prime }}{a}\frac{b\prime }{b} +\frac{a^{\prime \,^2}}{a^{2}}-\frac{a^{\prime }}{ra}\right) \nonumber \\&-r^{2}b^{3}s\ddot{s}\left( \frac{b^{\prime \prime }}{b}-2\frac{b^{\prime \,^{2}}}{b^{2}} -\frac{b^{\prime }}{rb}\right) \nonumber \\&+r^{2}aa^{\prime }b^{\prime \prime }\left( \frac{b^{\prime }}{b}+\frac{1}{r} \right) -r^{2}aba^{\prime \prime }\left( \frac{b^{\prime \,^{2}}}{b^2}+2\frac{b^{\prime }}{rb} \right) \nonumber \\&+raba^{\prime }\left( \frac{b^{\prime }}{rb}+2\frac{b^{\prime \,^{2}}}{b^2}\right) . \end{aligned}$$
(43)

For a given form of metric function b(r) (38), the class I condition in Eq. (43) is nonlinear. A physical relevant collapsing model must satisfy (30) and (43) simultaneously. It must be noted that simplest choice of solutions of (30) is a linear solution   [71]

$$\begin{aligned} s(t)= & {} -C_{Z} \,t , \end{aligned}$$
(44)

\(C_Z>0\). The form of the other metric function a(r) is obtained by using Eq. (38) and (44) into the class 1 condition (43)

$$\begin{aligned} a(r)= & {} \frac{1}{2\sqrt{2 C_{3}C_{4}}}[C_4^{2}\left( C_1b(r)+4C_2C_3\right) ^2-4C_Z^2]^{1/2} \end{aligned}$$
(45)

where \(C_1\) and \(C_2\) are integration constants. Surprisingly, the quantity in the numerator inside the square root, arises naturally from the matching of the Ricci scalar (and it’s derivative), given in (34) These forms of the solutions of the gravitational potentials are same as obtained in  [57] for shear free spacetime. In [57], it has been shown that for the static case, Karmarkar condition together with the pressure isotropy yields the Schwarzschild  [72] like form of the metric functions. Also, it has been shown that these set of gravitational potentials are the special class of those found in  [73]. Thus, although we have assume this particular form of \(\Delta \) (36) for the mathematical simplicity, represents the physically viable solutions.

It is now instructive to rewrite the physical quantities of the matter cloud in terms of the metric variables for a better understanding of the dynamics of spacetime during the collapse process. These expressions have been written in detail in the Appendix. The boundary condition (30) in the view of these equations in the Appendix, (63)–(65) becomes

$$\begin{aligned} 2s\,\ddot{s}+{\dot{s}}^2-2x{\dot{s}}+b\,s\,T^{D}_{rt}= & {} y-T^{D}_{rr} , \end{aligned}$$
(46)

where \(T^{D}_{rr}\) and \(T^{D}_{rt}\) are given by equations (31) and (32) respectively and the quantities x and y are

$$\begin{aligned} x= & {} \left( \frac{a^{\prime }}{b} \right) _{\Sigma },\nonumber \\ y= & {} \left( \frac{a^2}{b^2}\left[ \frac{b^{\prime \,^{2}}}{b^2}+\frac{2}{r}\frac{b^{\prime }}{b} +\frac{2a^{\prime }b^{\prime }}{ab}+\frac{2}{r}\frac{a^{\prime }}{a}\right] \right) _{\Sigma }. \end{aligned}$$
(47)

The metric functions a(r) and b(r) should not vanish during the collapsing phenomena, since otherwise the metric shall become degenerate. This also implies that their signatures remain unchanged. For second metric potential to be positive i.e. \(a(r)>0\) we must have, from (45), that

$$\begin{aligned} C_{Z}^2<C_{4}^{2}\left[ \frac{C_{1}}{C_{3}r^{2}+2C_{4}}-2C_{2}C_{3}\right] ^2. \end{aligned}$$
(48)

This equation also implies that at the center of the cloud, \(r=0\), we must have \(C_{Z}<C_{1}-2C_{2}C_{3}C_{4}\).

The graphical representations of the physical quantities (62)–(65) shows that they are well defined throughout the stellar collapse for both the f(R) models. Figures 1a–c, 2a–c, 3a–c and 4a–c shows that the density, radial pressure, tangential pressure and pressure anisotropy are positive and regular throughout the collapse for all the f(R) models and the parameters considered here. Also, as seen from the Fig. 5a–c, the heat flux increase as the collapse starts and remains positive throughout the collapse for these cases (Fig. 6).

Fig. 1
figure 1

a–c Shows the plots of the density \(\rho \) (62) w.r.t. time t and radial r coordinates for \(f(R)=R+\lambda R^{2}\) with \(\lambda =5\), \(f(R)=R^{1+\epsilon }\) with \(\epsilon =0.01\) and \(f(R)=R+\lambda \left[ \exp \left( -\gamma R\right) -1\right] \) with \(\lambda =0.1\) & \(\gamma =0.0002\) respectively. For each of these f(R) models, it remains regular as well as positive throughout the collapse

Full size image
Fig. 2
figure 2

a–c Shows the plots of the radial pressure \(p_{r}\) (63) w.r.t. time t and radial r coordinates for \(f(R)=R+\lambda R^{2}\) with \(\lambda =5\), \(f(R)=R^{1+\epsilon }\) with \(\epsilon =0.01\) and \(f(R)=R+\lambda \left[ \exp \left( -\gamma R\right) -1\right] \) with \(\lambda =0.1\) & \(\gamma =0.0002\) respectively. For each of these f(R) models, it remains regular as well as positive throughout the collapse

Full size image
Fig. 3
figure 3

a–c Shows the plots of the tangential pressure \(p_{t}\) (64) w.r.t. time t and radial r coordinates for \(f(R)=R+\lambda R^{2}\) with \(\lambda =5\), \(f(R)=R^{1+\epsilon }\) with \(\epsilon =0.01\) and \(f(R)=R+\lambda \left[ \exp \left( -\gamma R\right) -1\right] \) with \(\lambda =0.1\) & \(\gamma =0.0002\) respectively. For each of these f(R) models, it remains regular as well as positive throughout the collapse

Full size image
Fig. 4
figure 4

a–c Shows the plots of the pressure anisotropy \(\Delta \) (36) w.r.t. time t and radial r coordinates for \(f(R)=R+\lambda R^{2}\) with \(\lambda =5\), \(f(R)=R^{1+\epsilon }\) with \(\epsilon =0.01\) and \(f(R)=R+\lambda \left[ \exp \left( -\gamma R\right) -1\right] \) with \(\lambda =0.1\) & \(\gamma =0.0002\) respectively. For each of these f(R) models, it remains regular as well as positive throughout the collapse

Full size image
Fig. 5
figure 5

a–c Shows the plots of the heat flux q (65) w.r.t. time t and radial r coordinates for \(f(R)=R+\lambda R^{2}\) with \(\lambda =5\), \(f(R)=R^{1+\epsilon }\) with \(\epsilon =0.01\) and \(f(R)=R+\lambda \left[ \exp \left( -\gamma R\right) -1\right] \) with \(\lambda =0.1\) & \(\gamma =0.0002\) respectively. For each of these f(R) models, it remains positive throughout the collapse

Full size image
Fig. 6
figure 6

a Plot of the expansion scalar \(\Theta \) (66) w.r.t. time t and radial r coordinates. At the beginning of collapse \(\Theta \) has zero value and it starts decreasing and remains negative throughout the collapse. b Plot of the mass of the collapsing star (67) w.r.t. time t, and it shows that mass radiates linearly

Full size image

A related quantity of importance in this study is the total luminosity visible to an observer at infinity, which may be defined in terms of the mass loss from the boundary surface:

$$\begin{aligned} L_{\infty }= & {} -\left( \frac{dm}{dv}\right) _{\Sigma }\nonumber \\= & {} \left[ \frac{r^2\,s^{2}\,b^{2}\,p_{r}}{2} \left( 1+\frac{r\,b^{\prime }}{b}+\frac{r\,b\,{\dot{s}}}{a}\right) ^2\right] _{\Sigma }, \end{aligned}$$
(49)

where we have used the Eqs. (11), (27) and (28). Now, as soon as the black hole is formed, by definition, the luminosity of the surface is zero. From the above equation, this implies that sufficient condition for the formation of a black hole is

$$\begin{aligned} \left[ 1+\frac{r\,b^{\prime }}{b}+\frac{r\,b\,{\dot{s}}}{a}\right] _{\Sigma }=0. \end{aligned}$$
(50)

Naturally, for any static observer at asymptotic infinity, the redshift diverges at the time of formation of the black hole.

Fig. 7
figure 7

a–c Shows the plots of the luminosity (49) at \(r=r_{\Sigma }=1\), w.r.t. time t for \(f(R)=R+\lambda R^{2}\) with \(\lambda =5\), \(f(R)=R^{1+\epsilon }\) with \(\epsilon =0.01\) and \(f(R)=R+\lambda \left[ \exp \left( -\gamma R\right) -1\right] \) with \(\lambda =0.1\) & \(\gamma =0.0002\) respectively

Full size image

To show that these spacetime solutions are physically viable, we show that they satisfy the energy conditions as well. Indeed, all the energy conditions namely weak (W), null (N), dominant (D) and strong (S) hold good for the collapsing star. In the following we list these conditions [72, 74, 75]

E1 ::

\(\left( \rho +p_{r}\right) ^2-4q^{2}\,\ge \,0\)                                     (D/S/W)

E2 ::

\(\rho -p_{r}\,\ge \,0\,\,\,\, \)                         (D)

E3 ::

\(\rho -p_{r}-2p_{t}+\sqrt{\left( \rho +p_{r}\right) ^2-4q^{2}}\,\ge \,0\)                   (D)

E4 ::

\(\rho -p_{r}+\sqrt{\left( \rho +p_{r}\right) ^2-4q^{2}}\,\ge \,0\,\,\,\, \)                   (W/D)

E5 ::

\(\rho -p_{r}+2p_{t}+\sqrt{\left( \rho +p_{r}\right) ^2-4q^{2}}\,\ge \,0\)       (D/W/S)

E6 ::

\(2p_{t}+\sqrt{\left( \rho +p_{r}\right) ^2-4q^{2}}\,\ge \,0\)                               (S)

The star should also satisfy

E7 ::

\(\rho >0\),  \(p_{r}>0\),  \(p_{t}>0\), and \(\rho ^{\prime }<0\),  \(p_r^{\prime }<0\),  \(p_{t}^{\prime }<0\).

It is clear from the above conditions that E1,  E2, E3 & E7 are enough to validate the physical conditions existing inside the star. For the radiating- collapsing stellar models in f(R) gravity, Fig. 8a–c show that the energy conditions are positive and regular throughout the interior of the star.

Fig. 8
figure 8

a–c Shows the plots of energy conditions E1, E2 and E3 w.r.t. time t and radial r coordinates for \(f(R)=R+\lambda R^{2}\) with \(\lambda =5\) respectively

Full size image
Fig. 9
figure 9

a–c Shows the plots of energy conditions E1, E2 and E3 w.r.t. time t and radial r coordinates for \(f(R)=R^{1+\epsilon }\) with \(\epsilon =0.01\) respectively

Full size image
Fig. 10
figure 10

a–c Shows the plots of energy conditions E1, E2 and E3 w.r.t. time t and radial r coordinates for \(f(R)=R+\lambda \left[ \exp \left( -\gamma R\right) -1\right] \) with \(\lambda =0.1\) & \(\gamma =0.0002\) respectively

Full size image

2.1 Stability criteria

The study of dynamical instability (stability) of spherical stellar system shows that for adiabatic index \(\Gamma <4/3\) \(\left( \Gamma >4/3\right) \) the stellar system becomes unstable (stable) as the weight of the stellar system increase much faster (remains less than) than that of its pressure  [76]. Also, the causality condition imposes certain constraints on the dynamics of the stellar system such that inside the star, the radial \(V_{r}\) and the transverse \(V_{t}\) components of the speed of sound should be less than the speed of the light (\(c=1\)), so that \(0\le V_{r}\le 1\) and \(0\le V_{t}\le 1\)  [77]. Thus, to check the stability/instability of the collapsing stellar system, we need to study the behavior of the important physical quantities, adiabatic index, sound of the speed which are defined as [50, 77, 78]

$$\begin{aligned} \Gamma _{eff}= & {} \left[ \frac{\partial ( \ln p_{r})}{\partial ( \ln \rho )}\right] _{\Sigma } =\left[ \left( \frac{{\dot{p}}_{r}}{p_{r}}\right) \left( \frac{{\dot{\rho }}}{\rho }\right) \right] _{\Sigma } \end{aligned}$$
(51)
$$\begin{aligned} Vr= & {} \frac{d\,p_{r}}{d\rho }\,\, ,Vt=\frac{d\,p_{t}}{d\rho } \end{aligned}$$
(52)
Fig. 11
figure 11

a–c Shows the plots of adiabatic index at \(r=r_{\Sigma }=1\), w.r.t. time t for \(f(R)=R+\lambda R^{2}\) with \(\lambda =5\), \(f(R)=R^{1+\epsilon }\) with \(\epsilon =0.01\) and \(f(R)=R+\lambda \left[ \exp \left( -\gamma R\right) -1\right] \) with \(\lambda =0.1\) & \(\gamma =0.0002\) respectively

Full size image

Although stability may be understood from the behavior of the pressure and density variables, the quantities in (51) and (52) are considered to be better to establish stability. For \(f(R)=R+\lambda R^{2}\) model with \(\lambda =5\), Figs. 7a and 11a shows that the total luminosity and the adiabatic index are positive and increasing. Note that the adiabatic index attains a maximum value where the luminosity is maximum. This behavior of the luminosity and adiabatic index can be interpreted as follows. Any static observer at asymptotic infinity will see an exponentially radiating radial source until a time when luminosity reaches its maximum value after which it instantaneously turn off. This is due to the fact that the total mass of the star radiates linearly as seen from the Fig. 6b and when the star reaches its maximum luminosity, there is no mass left to radiates and hence the observer at rest at infinity will see sudden turn off of the light source. The similar kind of behavior were obtained in  [78]. Figure 11a shows that the effective adiabatic index is positive and less than 4/3 which implies that the considered stellar system is unstable and representing the collapsing scenario  [76]. For \(f(R)=R^{1+\epsilon }\), and \(f(R)=R+\lambda \,\left[ \exp {(-\sigma R)}-1\right] \) models, similar behavior of luminosity is obtained as that of for the first model. Figure 11b shows that the effective adiabatic index is constant function of time, and is positive and less than 4/3, which implies it represents the collapsing scenario. As we have shown graphically that the star radiates all its mass before reaching at the singularity. So, there are no trapped surfaces formed during the collapse. Which implies that neither the black hole nor naked singularity are the end state of the collapse.

2.2 Thermal properties

Earlier studies have shown that relaxation effects are important to understand dissipative gravitational collapse [59,60,61, 64, 79,80,81]. To study the temperature profiles, we consider the transport equation for the metric (6) given by  [58, 59, 82]

$$\begin{aligned} \tau h_{\mu }^{\nu }{\dot{q}}_{\nu }+q_{\mu }= & {} -k\left( h_{\mu }^{\nu } \nabla _{\nu } T+T{\dot{u}}_{\mu } \right) \end{aligned}$$
(53)
$$\begin{aligned} \tau \left( qbs\right) _{,t}+q\,a\,b\,s= & {} -\frac{k\left( aT \right) _{,r}}{bs}, \end{aligned}$$
(54)

where, \(\alpha >0\), \(\beta >0\),  \(\gamma >0\) and \(\sigma >0\), & \(h^{\mu \nu }=g^{\mu \nu }+u^{\mu }u^{\nu }\). Also,

$$\begin{aligned} \tau _{c} = \left( {\alpha /\gamma }\right) \,(T)^{-\sigma },k = \gamma \,\tau _{c}\, T^{3},\tau = \tau _{c}\,(\beta \,\gamma )/\alpha \end{aligned}$$
(55)

where \(\tau _{c}\) is the mean collision time, k is thermal conductivity and \(\tau \) represents the relaxation time respectively [59, 66]. The quantity \(\tau \) measures the strength of relaxational effects and is called the causality index. The values \(\tau =0\) or \(\beta =0\) represents the noncausal temperature profile.

Fig. 12
figure 12

a–c Shows the plots of the temperature profiles of the collapsing stellar system w.r.t. radial coordinate r for \(\sigma =0\) and for \(f(R)=R+\lambda R^{2}\) with \(\lambda =5\), \(f(R)=R^{1+\epsilon }\) with \(\epsilon =0.01\) and \(f(R)=R+\lambda \left[ \exp \left( -\gamma R\right) -1\right] \) with \(\lambda =0.1\) & \(\gamma =0.0002\) respectively

Full size image

Using conditions in Eq. (55), the the causal heat transport equation (54) becomes

$$\begin{aligned} \beta T^{-\sigma } \left( qbs\right) _{,t}+q\,a\,b\,s= & {} -\frac{\alpha \,\left( aT \right) _{,r}}{bs}\,\,T^{3-\sigma } . \end{aligned}$$
(56)

The noncausal solution of the heat transport equation (56), with \(\beta =0\) i.e. \(\tau =0\) are  [66]

$$\begin{aligned} \left( a\,T\right) ^{4}= & {} -\frac{4}{\alpha }\int a^{4}\,q\,b^2\,s^2\,dr+G(t), \qquad \quad \sigma = 0 \end{aligned}$$
(57)
$$\begin{aligned} \ln \left( a\,T\right)= & {} -\frac{1}{\alpha }\int q\,b^2\,s^2\,dr+G(t).\qquad \qquad \quad \sigma =4 \end{aligned}$$
(58)

The causal solution of the above heat transport equation (56) are  [66]

$$\begin{aligned} \left( a\,T\right) ^{4}= & {} -\frac{4}{\alpha }\left[ \beta \int a^3\,b\,s(q\,b\,s)_{,t}\,dr +\int a^{4}\,q\,b^2\,s^2\,dr\right] \nonumber \\&+G(t), \quad \sigma =0 \end{aligned}$$
(59)
$$\begin{aligned} \left( a\,T\right) ^{4}= & {} -\frac{4\beta }{\alpha }\exp \left( -\int \frac{4\,q\,b^{2}\,s^{2}}{\alpha }\,dr\right) \nonumber \\&\int a^3\,b\,s(q\,b\,s)_{,t}\,dr\,\exp \left( \int \frac{4\,q\,b^{2}\,s^{2}}{\alpha }\,dr\right) \nonumber \\&+G(t)\exp \left( -\int \frac{4\,q\,b^{2}\,s^{2}}{\alpha }\,dr\right) , \qquad \sigma =4 \end{aligned}$$
(60)

where G(t) appears as a function of integration and is determined by following boundary condition

$$\begin{aligned} \left( T^{4} \right) _{\Sigma }= & {} \left( \frac{L_{\infty }}{4\pi \delta r^{2}b^{2}s^{2}} \right) _{\Sigma }. \end{aligned}$$
(61)

where \(L_{\infty }\) is the total luminosity for an observer at infinity given by (49) and \(\delta >0\) is constant.

The Fig. 12a–c shows that for all these three f(R) models under consideration, the stellar system deviates from thermodynamical equilibrium due to the relaxation effects. Also, the causal and noncausal temperature profiles differ inside the interior of the star and causal temperature remains greater than that of the non-causal temperature.

3 Discussion of the results

In this paper, we investigated the dynamics of a collapsing matter configuration in f(R) gravity. The matter is assumed to be so massive that no stable compact objects like a white dwarf or a neutron star forms. The interior spacetime is smoothly matched to the outgoing radiation Vaidya metric across a timelike hypersurface. Incidentally, as has been noted earlier too, the matching conditions for the f(R) gravity is highly restrictive, since the geometric variables which are to matched here not only includes induced metric and the extrinsic curvatures, but also the trace of the extrinsic curvatures, and the Ricci scalar along with it’s time derivative. However, we have shown that all these matching conditions can be carried out consistently, leading to a spacetime solution which admits a collapsing scenario in which the matter cloud radiates heat flux, in such a manner that the entire matter is radiated out without forming a black hole. Although similar solutions have been reported earlier for GR, our solution incorporates these features into the collapsing models of three particular f(R) gravity models, while maintaining all the energy conditions. In particular, for all the three f(R) theories we have analyzed physical quantities like energy density, in Eq. (62), radial pressure (63) and tangential pressure in Eq. (64), pressure anisotropy, in Eq. (36) and it can be seen from the Figs. 1a–c, 2a–c, 3a–c and 4a–c that they are regular and positive throughout the collapse. From Fig. 5a–c it is also clear that the radial heat flux (65) is finite and positive throughout collapse. In particular, for \(f(R)=R+\lambda R^{2}\), Figs. 7a and 11a show that both total luminosity and the effective adiabatic index are positive and increasing and have maximum value where luminosity is maximum. This behavior of the luminosity and adiabatic index can be interpreted as follows: an observer at rest at infinity will see a exponential radiating radial source until it reaches time when luminosity reaches its maximum value and then instantaneous turn off of the radial source. This happens since the total mass of the star radiates linearly as seen from the Fig. 6b, and when the star reaches its maximum luminosity, there is no mass left to radiate and hence the observer at rest at infinity will see sudden turn off of the light source. The similar kind of behavior were obtained in GR too [78]. The Fig. 11a also shows that the effective adiabatic index is positive and less than 4/3 which implies that the considered stellar system is unstable and represents a collapsing phenomena. Also note that the Fig. 8a–c show that the energy conditions are positive and regular throughout the interior of the star.

For \(f(R)=R^{1+\epsilon }\) and \(f(R)=R+\lambda \,\left[ \exp {(-\sigma R)}-1\right] \) similar behavior is obtained as well. For example, Fig. 11b and c shows that the effective adiabatic index is constant function of time, and is positive and less than 4/3, and hence represents a collapsing scenario. We have shown graphically that the star radiates all its mass before reaching the singularity. So, there are no trapped surfaces formed during the collapse. This implies that neither a black hole nor a naked singularity exist at the end state of collapse. The Figs. 9a–c and 10a–c show that these model under consideration are physically viable. Also, the results obtained here reduces to those for GR for \(f(R)=R\) [57] (Figs. 11, 12).

Let us now comment on the nature of the central singularity. First, we note that the Ricci scalar (9) together with (44) imply that it diverges at \(t=0\), when all the matter has been radiated away. So, naturally the question arises regarding the gravitational strength of the central curvature singularity (a strong singularity would imply a naked singularity). However, in the present case, the curvature scalars go as \(t^{-2}\) which is precisely the sufficient condition for the singularity to be weak. So, our solutions represent a physically viable model where the spherically symmetric collapsing matter cloud undergoes gravitational collapse, which, during the collapse, also radiates away mass in the form of heat flux. The flux is radiated at such a rate that no horizon is ever formed and the central singularity is naked but gravitationally weak in nature.