Collapse in $f(R)$ gravity and the method of $R$ matching

Collapsing solutions in $f(R)$ gravity are restricted due to junction conditions that demand continuity of the Ricci scalar and its normal derivative across the time-like collapsing hypersurface. These are obtained via the method of $R$-matching, which is ubiquitous in $f(R)$ collapse scenarios. In this paper, we study spherically symmetric collapse with the modification term $\alpha R^2$, and use $R$-matching to exemplify a class of new solutions. After discussing some mathematical preliminaries by which we obtain an algebraic relation between the shear and the anisotropy in these theories, we consider two metric ansatzes. In the first, the collapsing metric is considered to be a separable function of the co-moving radius and time, and the collapse is shear-free, and in the second, a non-separable interior solution is considered, that represents gravitational collapse with non-zero shear viscosity. We arrive at novel solutions that indicate the formation of black holes or locally naked singularities, while obeying all the necessary energy conditions. The separable case allows for a simple analytic expression of the energy-momentum tensor, that indicates the positivity of the pressures throughout collapse, and is further used to study the heat flux evolution of the collapsing matter, whose analytic solutions are presented under certain approximations. These clearly highlight the role of modified gravity in the examples that we consider.

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.
This paper is organised as follows. In the next section, 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 section 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 section 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 section 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. Finally, section 6 ends this paper with a summary of the main results and some discussions.

II. MATHEMATICAL PRELIMINARIES AND SET UP
For a generic collapse scenario, in co-moving coordinates (t, r, θ, φ) the metric inside the spherically symmetric collapsing cloud is written as ds 2 − = −e 2ν(r,t) dt 2 + e 2ψ(r,t) dr 2 + Q 2 (r, t)dΩ 2 , where dΩ 2 = dθ 2 + sin 2 θdφ 2 . The metric outside the collapsing matter is usually represented by the Vaidya solution in terms of the retarded time u as 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 ds 2 + = −H(r)dt 2 + H(r) −1 dr 2 +r 2 dΩ 2 , H(r) = 1 − 2m r , 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) + L matter are given by: where We need to solve the modified Einstein equations with the energy momentum tensor 2 where we define the quantities where (, ) denote a symmetrization, and a semicolon denotes a covariant derivative. Here, ρ is the energy density, p r and p θ are the radial and tangential pressures, respectively, q µ = qn µ is the radial heat flow vector where n µ is a unit 4-vector along the radial direction, and u µ is the 4-velocity of the fluid. These satisfy n µ n µ = 1, u µ u µ = −1, u µ q µ = 0, u µ n µ = 0. Also, Θ = u ν ;ν is the expansion parameter and h µν = g µν + u µ u ν is the projection tensor. In this paper, we will be dealing with two situations, to be elaborated in sections 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 [20]) where we have defined the acceleration vector a µ = u µ ;ν u ν and σ αβ σ αβ = 2 3 σ 2 . This equation is valid only for f (R) models with d 2 F/dR 2 = 0, which is the case under consideration here. 3 Now, we will use the identity given by [21] and the well known relation between Riemann and Weyl tensors given by For f (R) gravity (recall that f (R) = αR 2 with d 2 F/dR 2 = 0), eq.(10) can be evaluated by using the results derived in [20] and after some algebra, we obtain (with R ;µ;ν ≡ R ;µν ), This generalises a corresponding result obtained for GR in [21]. Here, E µν is the electric part of the Weyl tensor, defined as E µν = C µνρλ u ρ u λ , with the magnetic part of the Weyl tensor vanishing identically due to spherical symmetry. Then, eliminating ρ from eq.(11) using eq.(8), we obtain where we have defined Equating eq.(12) and eq.(9) we get Finally, contracting with n γ n ν and denotingP = (Π + 2ησ), Expanding the left hand side of eq.(15), the evolution of the shear is given by the equation with a = n µ a µ . Eq.(16) is the most general evolution equation for the shear tensor in f (R) = αR 2 scenarios, with d 2 F/dR 2 = 0. The GR case corresponds here to α = 0 and has been analysed in [21]. We can make a few comments here. Now note that σ (being computed entirely from the metric) does not depend on the f (R) parameter α. This means that the the term in square brackets in eq.(16) has to be independent of α. For f (R) = αR 2 theories, this can be seen to imply that Eq. (17) gives an algebraic relation between the shear and the anisotropy in the f (R) theories that we consider. 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 [6] (see also [5]) by generalising a method developed in [22] 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 where we have defined In addition, the weak, dominant and strong energy conditions (WEC, DEC and SEC) are to be satisfied, and these are given respectively as 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, Finally, all the conditions above will need to be supplemented by the junction conditions in f (R) models [17], [18].
Recall that in GR, the standard Darmois-Israel junction conditions [16] are valid, which amount to matching of the first and second fundamental forms at the time-like hypersurface Σ : r = r 0 . These are defined, with a, b denoting the indices on the hypersurface, as, where N µ is the unit normal across the matching hypersurface, Here, e α a = ∂x α ∂y a are tangents to the matching hypersurface, and 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 [17], [18], so that the full set of matching conditions across the collapsing time-like hypersurface separating ds 2 − and ds 2 + are where [a] denotes the difference in the quantity a across the hypersurface Σ. 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 the these imply that the Misner-Sharp mass function [23], [24] given by equals the Schwarzschild mass when evaluated at the boundary Σ. From a fairly straightforward analysis, it is known that these also imply, from the metrics of eqs. (1) and (2) that Equivalently, the junction conditions imply that [19] 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).

III. 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 [18]. As another concrete example, suppose we assume that the interior metric is of the Lemaitre-Tolman-Bondi (LTB) form [25], [26], [27] given by where X(r, t) and Q(r, t) are functions of the co-moving radial coordinate and time, and dΩ 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 with k being a suitable constant. Also, the Einstein equations of GR can be shown to imply, for the general metric of eq.(28), with A(r) being an arbitrary function of the co-moving radial coordinate. 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 Since we are in co-moving coordinates we choose u µ = (1, 0, 0, 0) and n µ = (0, h(r)/a(t), 0, 0), so that the heat flux is along the radial direction, i.e q µ = qn µ . With this metric, for the model described the the Lagrangian density we can write down the energy momentum components: The Ricci scalar of the interior metric is calculated to be In order that the Ricci scalar matches smoothly to the collapsing co-moving boundary at all co-moving times, we will therefore require thatȧ 2 + 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 µ = 0, h(r)/a(t), 0, 0 , the condition In order to satisfy this for all times, one thus requiresȧ 2 + 2aä = 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 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 [28] to the study of 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. 5 Here and otherwise, a prime will refer to a derivative with respect to the radial coordinate.
Using eqs. (4) and (6), the relevant physical quantities are obtained for the metric of eq.(35) as Now, from the metric of eq.(35), it can be checked via eq.(37) that the pressure anisotropy is given by Importantly, if we demand that the pressure anisotropy vanishes identically, then it necessarily implies that R = 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 ≥ 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 choosė a(t) 2 +a(t)ä(t) = 0. To summarise, our ansatz for a solution of the metric of eq.(35) is (with b and n being constants), We will henceforth choose for simplicity, the constants b = 1/2 and n = 2 so that 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 it can be immediately seen that the condition N µ [T µν ] = 0 is satisfied at all times. Now upon using eqs. (36) and (39), we finally obtain the very simple expressions, 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 Θ 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 Θ > 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.
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 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 6 These diverge at the time of formation of the singularity, as expected 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 [5], [6] 1 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 experssions for these conditions become cumbersome. In this analysis,we choose α = 10 −3 . In figs.(1) to (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 These will satisfy the conditions u µ u µ = −1 and n µ n µ = 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 Also, using eq.(44), it can be checked that for our metric of eq.(35), the condition of eq.(15) is satisfied with With these inputs, it can be checked that eq.(16) is indeed satisfied in this case, with σ = 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 α = 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. [38], [30]). 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 [31]. 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.

IV. 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 A(r) = (1 − h(r)) −1 , and thus we have our ansatz for the interior metric 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 Since we want to match Ricci scalar and it's derivative across a junction, as the simplest choice, we puṫ Then, the Ricci scalar takes the simple form The solution of eq.(49) is 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 With this choice, from eq.(50), the Ricci scalar reads, (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 ∼ 1/r 2 . We will however concentrate on the singularity that forms due to the collapse process, in which case R ∼ 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 and the co-moving time for the formation of the apparent horizon t ah (r) is given by 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.
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 showing the central divergence at t s (r = 0) = 1. Also, we record the expression for the shear, In fig.(9), we show the behaviour of σ 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 σ ∼ r −1 at the origin for t → 1, as can be seen from eq.(57). In figs.(10), (11) and (12), we show the density ρ, the anisotropy Π = p θ − p r and the heat flux q as a 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 θ = π/2 here, so that σ 22 = σ 33 and hence we have to  (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)).
Morover, the anistropy 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 α and η. However, we will not go into the details here, as these are entirely similar to the situation that we have considered.

V. 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 [32]. 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 apriori no simple EOS that our co-moving observer will see, even in the simple case of the separable solutions presented in section 3. We will concentrate only on this class of solutions in this section, since the solution is section 4 does not correspond to realistic collaps of a dense star, as already mentioned.
First, we note that if we set the f (R) parameter α to zero, we have here Hence, at t = 0, the matter follows a barotropic equation of state with p r = p θ = γρ (remember that there is no pressure anisotropy with α → 0 as we have commented on at the end of section 3), with γ = 3/5. As the collapse proceeds, the barotropic index reduces in this case, and approaches 1/3 for t → 1. Hence, at the end of the collapse, with α → 0, the matter reduces to pure radiation.
In the general case, the situation is more complicated. Here, we have, from eq.(41), It is therefore seen that for small values of α (we have used α = 10 −3 in section 3), at the beginning of the collapse, the system is close to a barotropric fluid with γ = 0.6, but the effect of α 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 [33] (see also [34] for applications in the cosmological context). It can easily 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 [35]).
It is also of interest to consider the heat transport equation in our non-equilibrium collapsing scenario, following the pioneering work of [36]. 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 [36] (see also [37], [38]) Here, T is the local equilibrium temperature, κ is the thermal conductivity, and τ the relaxation timescale, and all these quantities must be positive, from physical conditions. Further, a µ is the acceleration vector defined after eq.(8).
In order to solve eq.(60), a number of assumptions is necessary, since κ and τ 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 [39], [40] (see also [41]) and assume that where β, γ, α 1 and σ 1 are non-negative constants, with the case β = 0 being the non-causal case (see, e.g [37] for an excellent exposition). For simplicity, we will restrict ourselves to cases with σ 1 ≤ 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 section 3 allows us to obtain analytic solutions at least in some approximations.
First of all, let us consider the non-causal case, and set β = 0. Then, we obtain the formal solution of eq.(60) as where F (t) is an apriori undetermined function of the co-moving time. The special cases σ 1 = 1, 3, 4 need to be solved separately. The results are where we have generically denoted an arbitrary function of the co-moving time by F (t). Eqs. (62) and (63) are the full set of solutions for the non-causal case, and the role of the f (R) parameter α can be easily 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 α is to decrease the temperature (compared to the α = 0 case) for σ 1 < 1 and σ 1 > 4, while it increases the temperature for 1 < σ 1 < 4. Also, from the first two equations of eq.(63), it is clear that close to the center, the effect of α vanishes for σ 1 = 1 and dominates for σ 1 = 3 with the term independent of α 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 → 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 σ 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 α is positive). In order to avoid this, we require 3 < σ 1 < 4 and the other values of σ 1 are ruled out in the class of models that we consider. Note also that the solutions for σ 1 = 3 and σ 1 = 4 do not suffer from this pathology, and hence our final set of admissible values of σ 1 is 3 ≤ σ 1 ≤ 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 → 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 ≤ σ 1 ≤ 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 [39], [40], [41]), if we ignore the last term on the right hand side of eq.(60), then with eq.(61), we obtain as a solution for σ 1 = 0, 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 section 3. For the nonseparable solutions of section 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.

VI. 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) = α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