# A Time Dependent Spacetime in $f(R,T)$ Gravity: Gravitational Collapse

Prabir Rudra

###### Abstract

In this note a time dependent spacetime is explored in the background of
$f(R,T)$ gravity via the gravitational collapse of a massive star. The star is
modelled by the Vaidya spacetime which is time dependent in nature. The
coupling of matter with curvature is the key feature of $f(R,T)$ theory and
here we have investigated its effects on a collapsing scenario. Two different
types of models, one involving minimal and the other involving non-minimal
coupling between matter and curvature are considered for our study. Power law
and exponential functionalities are considered as examples to check the
outcome of the gravitational collapse. Our prime objective is to explore the
nature of singularities (black hole or naked singularity) that form as an end
state of the collapse. Existence of outgoing radial null geodesics from the
central singularity was probed and such existence implied the formation of
naked singularities thus defying the cosmic censorship hypothesis. The absence
of such outgoing null geodesics would imply the formation of an event horizon
and the singularity formed becomes a black hole. Conditions under which such
possibilities occur are derived for all the models and sub-models.
Gravitational strength of the singularity is also investigated and the
conditions under which we can get a strong or a weak singularity is derived.
The results obtained are very interesting and may be attributed to the
coupling between curvature and matter. It is seen that for non-minimal
coupling there is a possibility of a globally naked singularity, whereas for a
minimal coupling scenario local nakedness is the only option. It is also found
that the singularity formed can be sufficiently weak in nature, which is
cosmologically desirable.

##  1 Introduction

For the last two decades we have been aware of the fact that our universe has
entered into a phase of accelerated expansion [1, 2]. Although this came as a
total surprise to the scientific community, extensive research have been able
to put some meaning to this observed phenomenon over the years. It is a widely
known fact that this cosmic acceleration can be explained via two different
theoretical frameworks. One is the theory of dark energy (DE) that recasts the
matter content of the universe to some exotic substance possessing negative
pressure. The other way is to modify Einstein’s theory of gravity leading to
modified gravity theories. The reader may refer to Refs.[3, 4, 5] for
extensive reviews on modified gravity theories and to the Ref.[6] for a
detailed review on DE.

One of the most popular way to modify Einstein’s gravity is by replacing the
Ricci scalar $R$ in the gravity Lagrangian of the Einstein-Hilbert action of
general relativity (GR) by an analytic function of $R$, i.e. $f(R)$, which
gives rise to $f(R)$ gravity theory. Extensive reviews on $f(R)$ gravity can
be found in Refs.[7, 8]. In Ref.[9] the authors proposed an even more generic
class of models by considering the gravitational lagrangian as an analytic
function of Ricci scalar $R$ and matter Lagrangian $L_{m}$, paving the path
for $f(R,L_{m})$ theories. Further developments in $f(R,L_{m})$ theories can
be found in Refs.[10, 11, 12]. In Ref.[13] Harko et al. proposed the $f(R,T)$
theory, where the matter Lagrangian is given by the trace $T$ of the energy-
momentum tensor $T_{\mu\nu}$. So here the gravitational Lagrangian is an
analytic function of two scalar invariants, namely the Ricci scalar $R$ and
the trace of the energy-momentum tensor $T$. Here the contributions of $T$
will come from the matter content of the universe. As a result it is found
that the field equations of $f(R,T)$ theory depends on a source term, which is
given by the variation of the energy-momentum tensor with respect to the
metric. This will in turn depend on the matter Lagrangian or the nature of
matter content of the universe. So it is obvious that for different types of
matter, such as scalar fields, perfect fluid, electromagnetic field, etc. we
will get different set of field equations. From the form of the function, it
is obvious that this theory involves coupling between matter and geometry. So
by studying this theory one can probe such coupling effects and their
consequences on various astrophysical and cosmological phenomenon. It is seen
that the covariant divergence of the energy-momentum tensor is non-zero for
this model, which leads to non-geodesic motion of the massive test particles.
This is because the coupling effects between matter and geometry induces an
extra acceleration on the particles. Thermodynamics in $f(R,T)$ gravity was
studied by Sharif and Zubair in [14]. Cosmological Evolution in $f(R,T)$
theory with collisional matter was studied in Ref.[15]. In Ref.[16] cosmic
coincidence problem was studied in the background of $f(R,T)$ gravity.
Cosmological models in $f(R,T)$ theories as phase space was explored in [17].
Dynamics of scalar perturbations in $f(R,T)$ gravity was studied by Alvarenga
et al. in [18]. Gravastars in $f(R,T)$ gravity was studied in Ref.[19]. Dark
matter from $f(R,T)$ gravity was investigated by the authors in [20].
Propagation of polar gravitational waves in $f(R,T)$ scenario was explored in
[21]. Dynamical behavior of the Tolman metrics in $f(R,T)$ gravity was studied
by Hansraj and Banerjee in [22].

Gravitational collapse is a key astrophysical phenomenon that helps us to
understand various aspects of the universe such as structure formation,
properties of stars, formation of black holes, white dwarfs, neutron stars,
etc. A star undergoes a gravitational collapse due to its own mass at the end
of its life cycle, when it has exhausted all its nuclear fuel. During its
collapse journey there are various stages at which the collapse may stop,
depending on the initial mass of the collapsing star. If the star is massive
i.e. mass$>20M_{\odot}$ ($M_{\odot}$ represents solar mass), then the collapse
does not come to a halt at any of the intermediate stages (such as white dwarf
or neutron star), but directly proceeds to form a singularity such as a black
hole (BH). The study of gravitational collapse started with Oppenheimer and
Snyder [23] in $1939$ when they explored the gravitational collapse of a dust
cloud modelled by a static Schwarzschild exterior and Friedmann interior.
Following this, Tolman [24] and Bondi [25] studied the collapse of spherically
symmetric inhomogeneous distribution of dust. Subsequently a lot of interest
was generated in this subject and numerous work related to this can be found
in literature. Some reviews in gravitational collapse can be found in
Refs.[26, 27]. Roger Penrose in 1969 proposed cosmic censorship hypothesis
(CCH) [28], where he stated that any cosmological singularity will always be
covered by an event horizon, thus censoring the singularity from an external
observer. Such a singularity (popularly called a black hole) is associated
with permanent loss of physical information allowing multiple physical states
to devolve into a single state. This is known as the black hole information
loss paradox [29, 30, 31, 32]. Over the years, in the absence of a formidable
proof of CCH, scientists started questioning its validity. As a consequence, a
search was initiated that will culminate in the discovery of a singularity
that will be free from any event horizon. This type od singularity will not
only disprove CCH but also in the absence of information loss it will enhance
our knowledge about gravity. Such a singularity is named as a naked
singularity (NS) [33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44] which is
considered to be a crucial tool in the formulation of an effective theory of
quantum gravity.

The first effective relativistic line element representing the spacetime of a
realistic star was given by P. C. Vaidya [45] in 1951. It represented the
radiation for a non-static mass, thus generalizing the static solution of
Schwarzschild. Schwarzschild’s solution basically represented the spacetime
around a spherically symmetric cold dark body with a constant mass. So it is
obvious that it could never model the spacetime outside a star. This is the
problem that was addressed by Vaidya in his phenomenal paper [45] of 1951. The
solution proposed by Vaidya was termed as Vaidya spacetime and is often
referred to as the shining or radiating Schwarzschild metric. It should be
noted that the basic difference between the two metrics is that the constant
mass parameter in the Schwarzschild metric is replaced by a time dependent
mass parameter in the Vaidya metric, which consequently becomes a time
dependent spacetime. Notable studies in Vaidya metric can be found in the
Refs.[46, 47, 48, 49, 50, 51, 52, 53, 54].

Here we are interested in exploring the gravitational collapse of a massive
star modelled by the Vaidya metric in the background of $f(R,T)$ gravity.
Collapsing scenario in the presence of coupling between matter and curvature
is expected to be an interesting proposition. Moreover the behaviour of Vaidya
spacetime has never been explored in the background of $f(R,T)$ gravity. So
there is more than enough motivation for attempting this work. We will
basically focus on the nature of the singularity formed (BH or NS) as the end
state of the collapse. We will report the conditions under which these
singularities can form in a comparative manner. We hope to obtain interesting
and new results in our collapsing scheme in the background of curvature-matter
coupling. In the next section we will report the basic equations of Vaidya
spacetime in $f(R,T)$ gravity and find solutions for the system. In section
III we will explore the collapsing scenario of a massive star. Section IV will
deal with the strength of the singularity formed and finally the paper will
end with a detailed discussion and conclusion in section V.

##  2 Vaidya spacetime in $f(R,T)$ gravity

The Einstein-Hilbert action for general relativity is given by,

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

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

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

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

$S_{f(R)}=\frac{1}{2\kappa}\int
f(R)\sqrt{-g}d^{4}x+\int\mathcal{L}_{m}\sqrt{-g}d^{4}x$ (2.3)

where $\mathcal{L}_{m}$ is the matter Lagrangian and the second integral on
the R.H.S is $S_{M}$ representing the matter fields. To obtain the action for
$f(R,T)$ gravity we further modify the action for $f(R)$ gravity by
introducing the trace of the energy-momentum tensor $T_{\mu\nu}$ in the
gravity Lagrangian as follows [13],

$S_{f(R,T)}=\frac{1}{2\kappa}\int
f(R,T)\sqrt{-g}d^{4}x+\int\mathcal{L}_{m}\sqrt{-g}d^{4}x$ (2.4)

Here $f(R,T)$ is an arbitrary function of the Ricci scalar $R$ and the trace
$T$ of the energy-momentum tensor $T_{\mu\nu}$. The energy-momentum tensor is
defined as [55],

$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}\mathcal{L}_{m})}{\delta
g^{\mu\nu}}$ (2.5)

The trace of this tensor can be given as $T=g^{\mu\nu}T_{\mu\nu}$. Taking
variation with respect to the metric we get the field equations for $f(R,T)$
gravity as,

$f_{R}(R,T)R_{\mu\nu}-\frac{1}{2}f(R,T)g_{\mu\nu}+\left(g_{\mu\nu}\Box-\nabla_{\mu}\nabla_{\nu}\right)f_{R}(R,T)=\kappa
T_{\mu\nu}-f_{T}(R,T)T_{\mu\nu}-f_{T}(R,T)\Theta_{\mu\nu}$ (2.6)

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

$\Theta_{\mu\nu}\equiv g^{\alpha\beta}\frac{\delta T_{\alpha\beta}}{\delta
g^{\mu\nu}}$ (2.7)

In the field equations $\nabla_{\mu}$ denotes covariant derivative associated
with the Levi-Civita connection of the metric and
$\Box\equiv\nabla^{\mu}\nabla_{\mu}$ is the D’Alembertian operator. Moreover
we have denoted $f_{R}(R,T)=\partial f(R,T)/\partial R$ and
$f_{T}(R,T)=\partial f(R,T)/\partial T$. The tensor $\Theta_{\mu\nu}$ can be
calculated as,

$\Theta_{\mu\nu}=-2T_{\mu\nu}+g_{\mu\nu}\mathcal{L}_{m}-2g^{\alpha\beta}\frac{\partial^{2}\mathcal{L}_{m}}{\partial
g^{\mu\nu}\partial g^{\alpha\beta}}$ (2.8)

It is seen that the above tensor depends on the matter lagrangian. For perfect
fluid the above tensor becomes,

$\Theta_{\mu\nu}=-2T_{\mu\nu}+pg_{\mu\nu}$ (2.9)

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

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

where $f(t,r)=-\left(1-\frac{m(t,r)}{r}\right)$ and using the units $G=c=1$.
The total energy momentum tensor of the field equation (2.6) is given by the
following sum,

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

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

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

and

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

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

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

and they satisfy the relations

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

Therefore, the non-vanishing components of the total energy-momentum tensor
will be as follows

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

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

$p=\omega\rho$ (2.17)

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

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

$\displaystyle
R_{00}=\frac{\left(m-r\right)m^{\prime\prime}+2\dot{m}}{2r^{2}},$
$\displaystyle~{}~{}~{}~{}R_{01}=R_{10}=\frac{m^{\prime\prime}}{2r}$
$\displaystyle R_{22}=m^{\prime},$
$\displaystyle~{}~{}R_{33}=m^{\prime}~{}\sin^{2}\theta$ (2.18)

where . and ′ represents the derivatives with respect to time coordinate
${}^{\prime}t^{\prime}$ and radial coordinate ${}^{\prime}r^{\prime}$
respectively. For this system the Ricci scalar becomes,

$R=\frac{2m^{\prime}+rm^{\prime\prime}}{r^{2}}$ (2.19)

The trace of the energy momentum tensor is calculated as,

$T=g^{\mu\nu}T_{\mu\nu}=2\left(\omega-1\right)\rho$ (2.20)

The relation between density and mass is considered as [56],

$\rho=n\times m(t,r)$ (2.21)

where $n>0$ is the particle number density.

###  2.1 Field equations

Now we consider some particular classes of $f(R,T)$ modified gravity models,
which are obtained by some explicit functional forms of $f(R,T)$. Since the
field equations depend on the nature of matter through the tensor
$\Theta_{\mu\nu}$, here we will consider the field equations for a perfect
fluid source, which will be our field of interest in this study, as discussed
in the previous section. On a broad sense we are going to discuss two types of
models.

####  2.1.1 Model-1: $f(R,T)=f_{1}(R)+f_{2}(T)$

Here we consider models of the form $f(R,T)=f_{1}(R)+f_{2}(T)$, where
$f_{1}(R)$ and $f_{2}(T)$ are arbitrary functions of $R$ and $T$ respectively.
It is straightforward to see that for $f_{1}(R)=R$ and $f_{2}(T)=0$, we can
retrieve GR from this model. Using Eq.(2.6), the gravitational field equations
for this model is given by,

$f_{1}^{\prime}(R)R_{\mu\nu}-\frac{1}{2}f_{1}(R)g_{\mu\nu}+\left(g_{\mu\nu}\Box-\nabla_{\mu}\nabla_{\nu}\right)f_{1}^{\prime}(R)=\kappa
T_{\mu\nu}+f_{2}^{\prime}(T)T_{\mu\nu}+\left[f_{2}^{\prime}(T)p+\frac{1}{2}f_{2}(T)\right]g_{\mu\nu}$
(2.22)

where ′ represents derivative with respect to the argument. Now using the
Vaidya metric given in Eq.(2.10) and using Eqs.(2.16), (2.17), (2.18) and
(2.19) in Eq.(2.25), we compute all the components of the Einstein’s field
equations for this model (taking $\kappa=1$). Here we report the $(01)$,
$(22)$ and $(33)$ components of the field equations which will be used in our
analysis. The rest of the components are reported in the appendix section of
the paper.

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

$-r\left\\{f_{1}(R)+f_{2}(T)-2\rho\left(1+f_{2}^{\prime}(T)-\omega
f_{2}^{\prime}(T)\right)\right\\}+f_{1}^{\prime}(R)m^{\prime\prime}=0$ (2.23)

2\. The $(22)$ and $(33)$ components of field equations are given by,

$r^{2}\left[f_{1}(R)+f_{2}(T)+2\omega\rho+4f_{2}^{\prime}(T)\omega\rho\right]-2f_{1}^{\prime}(R)m^{\prime}=0$
(2.24)

The above equations along with the ones reported in the appendix are the
Einstein’s field equations for $f(R,T)$ gravity in the time dependent Vaidya
spacetime for the first model.

####  2.1.2 Model-2: $f(R,T)=f_{1}(R)+f_{2}(R)f_{3}(T)$

Now we consider a second model given by $f(R,T)=f_{1}(R)+f_{2}(R)f_{3}(T)$,
where $f_{i}(R)$, $i=1,2$ are arbitrary functions of $R$ and $f_{3}(T)$ is an
arbitrary function of $T$. Here the scalar invariants $R$ and $T$ are non-
minimally coupled to each other via the second term. In order to realize GR
from this model, we should have $f_{1}(R)=R$ and either or both of $f_{2}(R)$
and $f_{3}(T)$ equal to zero. We may also take $f_{1}(R)=0,f_{2}(R)=R$ and
$f_{3}(T)=1$ to get GR from this model. Using Eq.(2.6), the gravitational
field equations for this model is given by,

$R_{\mu\nu}\left[f_{1}^{\prime}(R)+f_{2}^{\prime}(R)f_{3}(T)\right]-\frac{1}{2}f_{1}(R)g_{\mu\nu}+\left(g_{\mu\nu}\Box-\nabla_{\mu}\nabla_{\nu}\right)\left[f_{1}^{\prime}(R)+f_{2}^{\prime}(R)f_{3}(T)\right]$
$=\kappa
T_{\mu\nu}+f_{2}(R)f_{3}^{\prime}(T)T_{\mu\nu}+f_{2}(R)\left[f_{3}^{\prime}(T)p+\frac{1}{2}f_{3}(T)\right]g_{\mu\nu}$
(2.25)

Like the previous model here also we report the necessary components of the
field equations for the second model (taking $\kappa=1$):

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

$-r\left[f_{1}(R)-2\rho+f_{2}(R)\left\\{f_{3}(T)+2f_{3}^{\prime}(T)\rho\left(\omega-1\right)\right\\}\right]+\left(f_{1}^{\prime}(R)+f_{2}^{\prime}(R)f_{3}(T)\right)m^{\prime\prime}=0$
(2.26)

2\. The $(22)$ and $(33)$ components of field equations are given by,

$r^{2}\left(f_{1}(R)+f_{2}(R)f_{3}(T)+2\omega\rho+4f_{2}(R)f_{3}^{\prime}(T)\omega\rho\right)-2\left(f_{1}^{\prime}(R)+f_{2}^{\prime}(R)f_{3}(T)\right)m^{\prime}=0$
(2.27)

The rest of the components are reported in the appendix section. The above
equations along with the ones reported in the appendix are the Einstein’s
field equations for $f(R,T)$ gravity in the time dependent Vaidya spacetime
for the second model.

###  2.2 Solution of the system

In this section we will proceed to find solutions of the above systems. In
order to do that, we will have to consider special forms for the arbitrary
functions of $R$ and $T$ as examples. We will do this for both the models
studied in the previous section.

####  2.2.1 Model-1

Here we have consider four different toy models as examples in order to solve
the field equations. The model involves arbitrary functions of $R$ and $T$
coupled in a minimal way. The functional forms used in the toy models are
basically power law and exponential forms, which are observationally the most
favoured models with constraints imposed on their free parameters.

Case-1: $f_{1}(R)=g_{1}R^{\beta_{1}},~{}~{}f_{2}(T)=g_{2}T^{\beta_{2}}$, where
$g_{1},\beta_{1},g_{2},\beta_{2}$ are constants.

Here we have chosen power law forms for both $f_{1}(R)$ and $f_{2}(T)$. For
our convenience we call this the double-power (DP) model. For $g_{2}=0$ and
$g_{1}=\beta_{1}=1$, we get back GR from this model. For this case, the $22$
or $33$ component of the field equations give the following differential
equation,

$\frac{1}{r^{2}}\left[g_{1}\left\\{2\left(\beta_{1}-1\right)m^{\prime}-rm^{\prime\prime}\right\\}\left\\{\frac{2m^{\prime}+rm^{\prime\prime}}{r^{2}}\right\\}^{\beta_{1}-1}\right]+\frac{1}{1-\omega}\left[2^{\beta^{2}}g_{2}\left(2\beta_{2}\omega+\omega-1\right)\left\\{n\left(\omega-1\right)m\right\\}^{\beta_{2}}\right]$
$-2n\omega m=0$ (2.28)

Solving the above differential equation for $\beta_{1}=\beta_{2}=1$ we get,

$m(t,r)=h_{1}(t)~{}AiryAi\left[2^{1/3}r\left\\{\frac{n\left(g_{2}-\omega-3g_{2}\omega\right)}{g_{1}}\right\\}^{1/3}\right]+h_{2}(t)~{}AiryBi\left[2^{1/3}r\left\\{\frac{n\left(g_{2}-\omega-3g_{2}\omega\right)}{g_{1}}\right\\}^{1/3}\right]$
(2.29)

where $AiryAi$ and $AiryBi$ are the two Airy functions (see appendix) and
$h_{1}(t)$, $h_{2}(t)$ are arbitrary functions of time which arises from
integration. We would like to mention here that the imposed conditions
$\beta_{1}=\beta_{2}=1$ are necessary to get a solution of this system by the
known mathematical methods.

Now the $01$ component of the field equations gives the differential equation,

$2nr^{2}m-2^{\beta_{2}}g_{2}r^{2}\left(1+\beta_{2}\right)\left\\{n\left(\omega-1\right)m\right\\}^{\beta_{2}}+g_{1}\left(\frac{2m^{\prime}+rm^{\prime\prime}}{r^{2}}\right)^{\beta_{1}-1}\left\\{r\left(\beta_{1}-1\right)m^{\prime\prime}-2m^{\prime}\right\\}=0$
(2.30)

Solving this equation for $\beta_{1}=\beta_{2}=1$ we get,

$m(t,r)=h_{3}(t)e^{\frac{nr^{3}\left(1-2g_{2}\left(\omega-1\right)\right)}{3g_{1}}}$
(2.31)

where $h_{3}(t)$ is an arbitrary function of time. For some values of the
arbitrary functions it is expected that the solutions given by Eqs.(2.29) and
(2.31) will match. In fact it has been checked that both these solutions give
rise to similar scenarios in the gravitational collapse scheme which we are
going to introduce in the next section. So we are going adopt one of these
solutions depending upon the nature of genericity of the solution. Since the
solution given by Eq.(2.29) has two arbitrary functions we will use this for
our collapse study simply because it is more general in nature and can easily
generate the other solution for some well chosen initial conditions. From here
on, we will only concentrate on the differential equation and the solution
obtained from the $22$ or $33$ component of the field equations for the reason
discussed above. Now that we have obtained the mass parameter, using it in the
Vaidya metric given in Eq.(2.10), we can easily get the Vaidya spacetime for
the corresponding model in $f(R,T)$ gravity.

Case-2: $f_{1}(R)=g_{1}e^{\beta_{1}R},~{}~{}f_{2}(T)=g_{2}e^{\beta_{2}T}$,
where $g_{1},\beta_{1},g_{2},\beta_{2}$ are constants.

Here we have chosen exponential forms for both $f_{1}(R)$ and $f_{2}(T)$. For
our convenience we call this the double-exponential (DE) model. Realizing GR
from this model is difficult. Nevertheless an approximation will help us
realize the scenario. Expanding $e^{\beta_{1}R}$ in Taylor’s series and
keeping the linear terms in $R$ only, will help us realize GR for
$g_{2}=0,g_{1}=1$ and $\beta_{1}=\frac{R-1}{R}$. For this case, the $22$ or
$33$ component of the field equations gives us the differential equation,

$2nr^{2}\omega\left\\{1+2g_{2}\beta_{2}e^{2n\beta_{2}\left(\omega-1\right)m}\right\\}m+e^{\frac{\beta_{1}\left(2m^{\prime}+rm^{\prime\prime}\right)}{r^{2}}}g_{1}\left(r^{2}-2\beta_{1}m^{\prime}\right)+g_{2}r^{2}e^{2n\beta_{2}\left(\omega-1\right)m}=0$
(2.32)

This equation has got the unknown function $m$ and its derivatives in
exponential form. It is not possible to find a general solution of this
equation. So we search for approximate solutions. We expand the exponentials
in the first and the third term in Taylor series and take the linear terms
only to get the following solution for $\beta_{1}=0$,

$m(t,r)=\frac{g_{2}n\beta_{2}\left(1-3\omega\right)-n\omega+\sqrt{n^{2}\left[\omega^{2}+g_{2}^{2}\beta_{2}^{2}\left(1+\omega\right)^{2}+2g_{2}\beta_{2}\omega\left\\{3\omega-4g_{1}\beta_{2}\left(\omega-1\right)-1\right\\}\right]}}{8g_{2}n^{2}\beta_{2}^{2}\omega\left(\omega-1\right)}$
(2.33)

We see that this a constant solution for the mass parameter. Using this in
Eq.(2.10) we will get the Vaidya spacetime in $f(R,T)$ gravity for this case.

Case-3: $f_{1}(R)=g_{1}R^{\beta_{1}},~{}~{}f_{2}(T)=g_{2}e^{\beta_{2}T}$,
where $g_{1},\beta_{1},g_{2},\beta_{2}$ are constants.

Here we have chosen power law for $f_{1}(R)$ and exponential form for
$f_{2}(T)$. For our convenience we call this the power-exponential (PE) model.
For $g_{2}=0$ and $g_{1}=\beta_{1}=1$, we get back GR from this model. The
$22$ or $33$ component of the field equations gives us the differential
equation,

$2nr^{2}\left(1+2g_{2}\beta_{2}e^{2n\beta_{2}\left(\omega-1\right)m}\right)\omega
m+e^{2n\beta_{2}\left(\omega-1\right)m}g_{2}r^{2}-g_{1}\left\\{2\left(\beta_{1}-1\right)m^{\prime}-rm^{\prime\prime}\right\\}\left(\frac{2m^{\prime}+rm^{\prime\prime}}{r^{2}}\right)^{\beta_{1}-1}=0$
(2.34)

For $\beta_{1}=1$ and $\beta_{2}=0$ we get the following solution for the
above differential equation,

$m(t,r)=\frac{1}{2n\omega}\left[g_{2}\pi
AiryAi^{\prime}\left[2^{1/3}r\left(-\frac{n\omega}{g_{1}}\right)^{1/3}\right]AiryBi\left[-\frac{2^{1/3}nr\omega}{g_{1}\left(-\frac{n\omega}{g_{1}}\right)^{2/3}}\right]-g_{2}\pi
AiryAi\left[-\frac{2^{1/3}nr\omega}{g_{1}\left(-\frac{n\omega}{g_{1}}\right)^{2/3}}\right]\times\right.$

$\left.AiryBi^{\prime}\left[2^{1/3}r\left(-\frac{n\omega}{g_{1}}\right)^{1/3}\right]\right]+h_{4}(t)AiryAi\left[-\frac{2^{1/3}nr\omega}{g_{1}\left(-\frac{n\omega}{g_{1}}\right)^{2/3}}\right]+h_{5}(t)AiryBi\left[-\frac{2^{1/3}nr\omega}{g_{1}\left(-\frac{n\omega}{g_{1}}\right)^{2/3}}\right]$
(2.35)

where $AiryAi^{\prime}$, $AiryBi^{\prime}$ are derivatives of the Airy
functions with respect to the argument and $h_{4}(t)$, $h_{5}(t)$ are
arbitrary functions of time.

Case-4: $f_{1}(R)=g_{1}e^{\beta_{1}R},~{}~{}f_{2}(T)=g_{2}T^{\beta_{2}}$,
where $g_{1},\beta_{1},g_{2},\beta_{2}$ are constants.

Here we have chosen exponential form for both $f_{1}(R)$ and power law for
$f_{2}(T)$. For our convenience we call this the exponential-power (EP) model.
A similar scenario as discussed in case-2, will help us realize GR from this
model. The $22$ or $33$ component of the field equations gives us the
differential equation,

$r^{2}\left[-2n\omega
m-\frac{2^{\beta_{2}}g_{2}\left(\omega+2\beta_{2}\omega-1\right)\left\\{n\left(\omega-1\right)m\right\\}^{\beta_{2}}}{\omega-1}\right]-g_{1}\left(r^{2}-2\beta_{1}m^{\prime}\right)e^{\frac{\beta_{1}\left(2m^{\prime}+rm^{\prime\prime}\right)}{r^{2}}}=0$
(2.36)

For $\beta_{1}=0$ and $\beta_{2}=1$ we get the following solution of the above
equation,

$m(t,r)=-\frac{g_{1}}{2n\left(\omega+3g_{2}\omega-g_{2}\right)}$ (2.37)

####  2.2.2 Model-2

Now again we consider some special models as sub-cases in order to solve the
field equations. The basic difference between this model with the previous one
is that here the functions of $R$ and $T$ will be minimally coupled to each
other which is observationally the more favoured model.

Case-1:
$f_{1}(R)=g_{1}R^{\beta_{1}},~{}f_{2}(R)=g_{2}R^{\beta_{2}},~{}f_{3}(T)=g_{3}T^{\beta_{3}}$
($g_{1},\beta_{1},g_{2},\beta_{2},g_{3},\beta_{3}$ are constants)

Here we have considered power law forms for all the three functions. We call
this the triple-power (TP) model. For $g_{1}=\beta_{1}=1$ and $g_{2}=0$ or
$g_{3}=0$, we can realize GR from this model. We may also realize GR for
$g_{1}=\beta_{3}=0$ and $g_{2}=g_{3}=\beta_{2}=1$. The $22$ or $33$ component
of the field equations gives us the differential equation,

$2n\omega
m\left(2m^{\prime}+rm^{\prime\prime}\right)+g_{1}\left(\frac{2m^{\prime}+rm^{\prime\prime}}{r^{2}}\right)^{\beta_{1}}\left\\{2\left(1-\beta_{1}\right)m^{\prime}+rm^{\prime\prime}\right\\}$
$+\frac{2^{\beta_{3}}g_{2}g_{3}\left\\{n\left(\omega-1\right)m\right\\}^{\beta_{3}}\left(\frac{2m^{\prime}+rm^{\prime\prime}}{r^{2}}\right)^{\beta_{2}}\left\\{2\left(\beta_{2}+\omega-\beta_{2}\omega+2\beta_{3}\omega-1\right)m^{\prime}+r\left(\omega+2\beta_{3}\omega-1\right)m^{\prime\prime}\right\\}}{\omega-1}=0$
(2.38)

A solution for the above differential equation can be obtained for
$\beta_{1}=1,\beta_{3}=2$ and $\omega=1$ which is given below,

$m(t,r)=h_{6}(t)AiryAi\left[-\frac{2^{1/3}nr}{g_{1}\left(-n/g_{1}\right)^{2/3}}\right]+h_{7}(t)AiryBi\left[-\frac{2^{1/3}nr}{g_{1}\left(-n/g_{1}\right)^{2/3}}\right]$
(2.39)

where $h_{6}(t)$ and $h_{7}(t)$ are arbitrary functions of time. This solution
corresponds to early universe ($\omega=1$) representing stiff perfect fluid.

Case-2:
$f_{1}(R)=g_{1}e^{\beta_{1}R},~{}~{}f_{2}(R)=g_{2}e^{\beta_{2}R},~{}~{}f_{3}(T)=g_{3}T^{\beta_{3}}$
($g_{1},\beta_{1},g_{2},\beta_{2},g_{3},\beta_{3}$ are constants)

This is the double-exponential-power (DEP) model. For this model the $22$ or
$33$ component of the field equations yields the following differential
equation,

$2nr^{2}\left(\omega-1\right)\omega
m+g_{1}\left(\omega-1\right)\left(r^{2}e^{\frac{\beta_{1}\left(2m^{\prime}+rm^{\prime\prime}\right)}{r^{2}}}-2\beta_{1}e^{\beta_{1}R}m^{\prime}\right)+2^{\beta_{3}}g_{2}g_{3}\left\\{n\left(\omega-1\right)m\right\\}^{\beta_{3}}\times$
$\left\\{r^{2}\left(\omega+3\beta_{3}\omega-1\right)e^{\frac{\beta_{2}\left(2m^{\prime}+rm^{\prime\prime}\right)}{r^{2}}}-2\beta_{2}\left(\omega-1\right)e^{\beta_{2}R}m^{\prime}\right\\}=0$
(2.40)

For $\beta_{1}=0,\beta_{3}=2$ and $\omega=1$ we get the following constant
solution of the above equation,

$m(t,r)=-\frac{g_{1}}{2n}$ (2.41)

We see that this is a constant solution. Moreover this solution is valid in
the early universe for a stiff perfect fluid ($\omega=1$).

Case-3:
$f_{1}(R)=g_{1}R^{\beta_{1}},~{}~{}f_{2}(R)=g_{2}R^{\beta_{2}},~{}~{}f_{3}(T)=g_{3}e^{\beta_{3}T}$
($g_{1},\beta_{1},g_{2},\beta_{2},g_{3},\beta_{3}$ are constants)

From the choice of the functions we can see that this is a double-power-
exponential (DPE) model. For $g_{2}=0$ or $g_{3}=0$ and $g_{1}=\beta_{1}=1$ we
get back GR from this model. In this case the $(22)$ or $(33)$ components of
the field equations gives the differential equation,

$r^{2}\left[2n\omega
m+g_{1}\left(\frac{2m^{\prime}+rm^{\prime\prime}}{r^{2}}\right)^{\beta_{1}}+g_{2}g_{3}\left(1+4n\beta_{3}\omega
m\right)e^{2n\beta_{3}\left(\omega-1\right)m}\left(\frac{2m^{\prime}+rm^{\prime\prime}}{r^{2}}\right)^{\beta_{2}}\right]$
$-2m^{\prime}\left[g_{1}\beta_{1}\left(\frac{2m^{\prime}+rm^{\prime\prime}}{r^{2}}\right)^{\beta_{1}-1}+g_{2}g_{3}\beta_{2}e^{2n\beta_{3}\left(\omega-1\right)m}\left(\frac{2m^{\prime}+rm^{\prime\prime}}{r^{2}}\right)^{\beta_{2}-1}\right]=0$
(2.42)

The following solution for the above equation is obtained for
$\beta_{1}=\beta_{2}=1$ and $\omega=0$,

$m(t,r)=h_{8}(t)+h_{9}(t)r~{}~{}~{}~{}~{}~{}~{}~{}OR~{}~{}~{}~{}~{}~{}m(t,r)=\frac{\log\left(-\frac{g_{2}g_{3}}{g_{1}}\right)}{2n\beta_{3}}$
(2.43)

where $h_{8}(t)$ and $h_{9}(t)$ are arbitrary functions of time. We will use
the first expression for the mass parameter for further study because it is
evolving with $r$ and $t$ and hence is more informative for our study.
Moreover the second expression being a constant can always be realized from
the first expression using suitable initial conditions. In this sense the
first expression is more generalized and so we intend to use it in our
analysis.

Case-4:
$f_{1}(R)=g_{1}e^{\beta_{1}R},~{}~{}f_{2}(R)=g_{2}e^{\beta_{2}R},~{}~{}f_{3}(T)=g_{3}e^{\beta_{3}T}$
($g_{1},\beta_{1},g_{2},\beta_{2},g_{3},\beta_{3}$ are constants)

This is the triple-exponential (TE) model formed by three exponential
functions. For this model the $(22)$ or $(33)$ components of the field
equations give the differential equation,

$g_{1}r^{2}e^{\frac{\beta_{1}\left(2m^{\prime}+rm^{\prime\prime}\right)}{r^{2}}}+2nr^{2}\omega
m-2g_{1}\beta_{1}e^{\frac{\beta_{1}\left(2m^{\prime}+rm^{\prime\prime}\right)}{r^{2}}}m^{\prime}+g_{2}g_{3}\left(r^{2}+4nr^{2}\beta_{3}\omega
m-2\beta_{2}m^{\prime}\right)e^{2n\beta_{3}\left(\omega-1\right)m+\frac{\beta_{2}\left(2m^{\prime}+rm^{\prime\prime}\right)}{r^{2}}}=0$
(2.44)

A solution to the above equation is obtained for $\beta_{1}=\beta_{2}=1$ and
$\omega=0$ which is given below,

$m(t,r)=\frac{r^{3}}{6}+h_{10}(t)~{}~{}~{}~{}~{}~{}~{}OR~{}~{}~{}~{}~{}~{}~{}m(t,r)=\frac{\log\left(-\frac{g_{2}g_{3}}{g_{1}}\right)}{2n\beta_{3}}$
(2.45)

where $h_{10}(t)$ is an arbitrary function of time. Just like the previous
model, here also we will use the first expression for the mass parameter, for
reasons similar to the ones discussed in the previous model. It should be
noted that this solution corresponds to dust ($\omega=0$) as far as the matter
content of the universe is concerned and cosmologically this corresponds to
early universe.

##  3 Gravitational Collapse

In this section, we will devise a mechanism in order to study the
gravitational collapse of a massive star in this system. As mentioned before,
we will consider that the parent star is a massive one so that collapse
smoothly continues until a singularity (BH or NS) is formed. Our idea is to
develop a set-up, via which the nature of the singularity (BH or NS) can be
comprehensively identified. At least our aim is to derive a condition that
will govern the nature of singularity (BH or NS) resulting out of the
gravitational collapse.

Let us consider a spherical collapsing system, where the physical radius of
the $r$-th shell of the star at time $t$ is $R(t,r)$. A suitable initial
condition would be that in the epoch $t=0$, we have $R(0,r)=r$. It is obvious
that if the collapse is inhomogeneous, then different collapsing shells may
become singular at different times. We are concerned with the light photons
emerging from the singularity and travelling along the geodesics and reaching
an external observer. An event horizon will be an obstruction for these
photons and will resist them from reaching the observer. So here we will probe
the existence of such outgoing non-spacelike geodesics. Theoretically if such
geodesics possess well defined tangent at the singularity, the quantity
$dR/dr$ will definitely tend towards a finite limit with the geodesics
approaching the singularity in the past following the trajectories. When these
trajectories reach the points $(t_{0},r)=(t_{0},0)$, there is a complete
breakdown of mathematical and physical concepts and a singularity occurs at
$R(t_{0},0)=0$. At these points ideally the collapsing matter shells are
crushed to zero radius, which results in the formation of the central
singularity. This is a highly compact object, since a huge amount of mass is
packed inside an almost negligible volume. Now if we follow back the path of
the outgoing non-spacelike geodesics that are emerging from the central
singularity, it is highly probable that they will terminate in the past at the
singularity $(r=0,t=t_{0})$ where $R(t_{0},0)=0$. Therefore from our set-up,
mathematically we should have $R\rightarrow 0$ as $r\rightarrow 0$ [57].

We obtain the equation for outgoing radial null geodesics from the Vaidya
metric (2.10) by putting $ds^{2}=0$ and
$d\Omega_{2}^{2}=d\theta^{2}+\sin^{2}\theta d\phi^{2}=0$ as furnished below

$\frac{dt}{dr}=\frac{2}{\left(1-\frac{m(t,r)}{r}\right)}.$ (3.1)

The above differential equation has a singularity at $r=0,~{}t=0$.
Mathematically this means that any solution to the above equation is not
analytic at the point $r=0,~{}t=0$. Since there is a mathematical breakdown at
the singularity we are forced to study the limiting behaviour as one
approaches the singularity. To facilitate this, we consider a parameter
$X=t/r$. The idea is to study the limiting behaviour of the function $X$ as we
approach the singularity at $r=0,~{}t=0$ following the radial null geodesic.
If we denote the limiting value of $X$ by $X_{0}$ then using L’Hospital’s rule
we have

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

This will actually generate an algebraic equation in terms of $X_{0}$. The
roots of this equation will be our prime concern because they actually
represent the slopes (direction) of the tangents to the geodesics. Here we are
only interested in the real roots because we are dealing with a realistic
collapsing scenario with no connection to the complex domain. For our set-up,
any positive real root of this algebraic equation will give the direction of
the tangent to an outgoing null geodesic at the singularity. Therefore the
existence of positive real roots of this equation corresponds to a necessary
and sufficient condition for the singularity to be naked in nature. Now as we
have discussed earlier, if a single null geodesic in the $(t,r)$ plane escapes
the singularity, it would mean that a single wavefront emitted from the
singularity reaches the external observer. In such a scenario the singularity
would be visible only instantaneously to a distant observer and become a
locally naked singularity. Physically this will correspond to a situation
where the event horizon was eliminated from the picture, but only temporarily.
But this might not be enough for a complete exchange of information between
the singularity and the observer. So for a formidable exchange of information,
the singularity is to be seen for a finite period of time. This requires a
family of null geodesics escaping from the singularity thus making it globally
naked. In our mathematical set-up this can be investigated very easily from
the number of real positive roots obtained from the above algebraic equation.
The above explained comprehensive mathematical set-up for identifying the
nature of singularity formed as the end state of a gravitational collapse was
first used by Joshi, Singh and Dwivedi in several of their papers [38, 39, 57,
37]. With the mathematical tools ready, we proceed to study the models one by
one.

###  3.1 Model-1

####  3.1.1 Case-1

Using Eq.(2.29) in Eq.(3.17) we get,

$\frac{2}{X_{0}}=\begin{array}[]{l}lim\\\ \tiny t\rightarrow 0\\\ \tiny
r\rightarrow
0\end{array}\left[1-\frac{h_{1}(t)}{r}~{}AiryAi\left[2^{1/3}r\left\\{\frac{n\left(g_{2}-\omega-3g_{2}\omega\right)}{g_{1}}\right\\}^{1/3}\right]\right.$
$\displaystyle\left.-\frac{h_{2}(t)}{r}~{}AiryBi\left[2^{1/3}r\left\\{\frac{n\left(g_{2}-\omega-3g_{2}\omega\right)}{g_{1}}\right\\}^{1/3}\right]\right]$
(3.18)

Here we will consider self-similar collapsing scenario. So we consider the
following self-similar expressions for the arbitrary functions $h_{i}(t)$,
$i=1,2$

$h_{1}(t)=\xi_{1}t,~{}~{}~{}~{}~{}~{}~{}~{}h_{2}(t)=\xi_{2}t$ where $\xi_{1}$
and $\xi_{2}$ are arbitrary constants.

Using the above chosen functions in Eq.(3.18) we get the following algebraic
equation in $X_{0}$,

$\frac{1}{\Gamma(2/3)}\left(\frac{\xi_{1}}{3^{2/3}}+\frac{\xi_{2}}{3^{1/6}}\right)X_{0}^{2}-X_{0}+2=0$
(3.19)

To evaluate the above limit, we have used the values of the Airy functions
given by Eq.(6.8) in the appendix. Solving the above equation we get two
values of $X_{0}$ which are,

$X_{0_{1,2}}^{case1}=\frac{\sqrt{\Gamma(2/3)}}{2\sqrt{3}\left(\frac{\xi_{1}}{3^{2/3}}+\frac{\xi_{2}}{3^{1/6}}\right)}\left[\sqrt{3\Gamma(2/3)}\pm\sqrt{3\Gamma(2/3)-8\times
3^{1/3}\xi_{1}-8\times 3^{5/6}\xi_{2}}\right]$ (3.20)

Here we have considered and henceforth we will consider positive sign for
root1 and negative sign for root2. Since we are dealing with a realistic
situation we should have
$3^{1/3}\xi_{1}+3^{5/6}\xi_{2}\leq\frac{3\Gamma(2/3)}{8}$. Now in order to get
a NS we should have $X_{0}>0$. We list below the respective conditions in
detail.

Conditions for a local NS: $X_{0_{1}}^{case1}>0$ $\&$ $X_{0_{2}}^{case1}<0$ OR
$X_{0_{1}}^{case1}<0$ $\&$ $X_{0_{2}}^{case1}>0$

Conditions for global NS: $X_{0_{1}}^{case1}>0$ $\&$ $X_{0_{2}}^{case1}>0$

Condition for BH: $X_{0_{1}}^{case1}<0$ $\&$ $X_{0_{2}}^{case1}<0$

We see that the above conditions put constraints of $\xi_{1}$ and $\xi_{2}$.
So by clubbing this theory with observations of collapsing massive stars, we
can get bounds on the model parameters. The roots $X_{0_{1,2}}^{case1}$ have
been plotted against the parameters $\xi_{1}$ and $\xi_{2}$ in Figs.(1) and
(2).

Fig.1 Fig.2

Figs.1 and 2 show the variation of the collapse parameter $X_{0}$ for
different values of $\xi_{1}$ and $\xi_{2}$ for Case-1 of Model-1. Fig.1 shows
the variation for the first root $X_{0_{1}}^{case1}$, whereas Fig.2 shows the
variation for the second root $X_{0_{2}}^{case1}$.

####  3.1.2 Case-2

Using Eq.(2.33) in Eq.(3.17) we get,

$\frac{2}{X_{0}}=\begin{array}[]{l}lim\\\ \tiny t\rightarrow 0\\\ \tiny
r\rightarrow
0\end{array}\left[1-\frac{g_{2}n\beta_{2}\left(1-3\omega\right)-n\omega+\sqrt{n^{2}\left[\omega^{2}+g_{2}^{2}\beta_{2}^{2}\left(1+\omega\right)^{2}+2g_{2}\beta_{2}\omega\left\\{3\omega-4g_{1}\beta_{2}\left(\omega-1\right)-1\right\\}\right]}}{8g_{2}n^{2}\beta_{2}^{2}\omega\left(\omega-1\right)r}\right]$
. (3.21)

Evaluating the above limit we get $\frac{2}{X_{0}}\rightarrow\infty$, which
implies $X_{0}\rightarrow 0$. Since the mass function in this case is not a
function of $t$ and $r$, we do not get a realistic collapsing scenario for
this particular model according to our scheme of study.

####  3.1.3 Case-3

Using Eq.(2.35) in Eq.(3.17) we get,

$\frac{2}{X_{0}}=\begin{array}[]{l}lim\\\ \tiny t\rightarrow 0\\\ \tiny
r\rightarrow 0\end{array}\left[1-\frac{1}{2n\omega r}\left\\{g_{2}\pi
AiryAi^{\prime}\left[2^{1/3}r\left(-\frac{n\omega}{g_{1}}\right)^{1/3}\right]\times
AiryBi\left[-\frac{2^{1/3}nr\omega}{g_{1}\left(-\frac{n\omega}{g_{1}}\right)^{2/3}}\right]\right.\right.$

$\left.\left.-g_{2}\pi
AiryAi\left[-\frac{2^{1/3}nr\omega}{g_{1}\left(-\frac{n\omega}{g_{1}}\right)^{2/3}}\right]\times
AiryBi^{\prime}\left[2^{1/3}r\left(-\frac{n\omega}{g_{1}}\right)^{1/3}\right]\right\\}-\frac{h_{4}(t)}{r}AiryAi\left[-\frac{2^{1/3}nr\omega}{g_{1}\left(-\frac{n\omega}{g_{1}}\right)^{2/3}}\right]\right.$

$\displaystyle\left.-\frac{h_{5}(t)}{r}AiryBi\left[-\frac{2^{1/3}nr\omega}{g_{1}\left(-\frac{n\omega}{g_{1}}\right)^{2/3}}\right]\right]$
(3.22)

We consider the following functions:
$h_{4}(t)=\xi_{4}t,~{}~{}h_{5}(t)=\xi_{5}t$, where $\xi_{4}$ and $\xi_{5}$ are
arbitrary constants. Using the above chosen functions in Eq.(3.22) we get the
following algebraic equation in $X_{0}$,

$\frac{1}{\Gamma(2/3)}\left(\frac{\xi_{4}}{3^{2/3}}+\frac{\xi_{5}}{3^{1/6}}\right)X_{0}^{2}+\left(\xi_{3}-1\right)X_{0}+2=0$
(3.23)

where $\xi_{3}$ is a constant that arises as a limiting value of the second
term of the expression in Eq.(3.21). Solving the above quadratic we get,

$X_{0_{1,2}}^{case3}=\frac{\Gamma(2/3)\left[1-\xi_{3}\pm\sqrt{\left(\xi_{3}-1\right)^{2}-\frac{8}{\Gamma(2/3)}\left(\frac{\xi_{4}}{3^{2/3}}+\frac{\xi_{5}}{3^{1/6}}\right)}\right]}{2\left(\frac{\xi_{4}}{3^{2/3}}+\frac{\xi_{5}}{3^{1/6}}\right)}$
(3.24)

where we should have
$\left(\xi_{3}-1\right)^{2}\geq\frac{8}{\Gamma(2/3)}\left(\frac{\xi_{4}}{3^{2/3}}+\frac{\xi_{5}}{3^{1/6}}\right)$.
The conditions for NS or BH will be similar as discussed in Case-1.

Conditions for a local NS: $X_{0_{1}}^{case3}>0$ $\&$ $X_{0_{2}}^{case3}<0$ OR
$X_{0_{1}}^{case3}<0$ $\&$ $X_{0_{2}}^{case3}>0$

Conditions for global NS: $X_{0_{1}}^{case3}>0$ $\&$ $X_{0_{2}}^{case3}>0$

Condition for BH: $X_{0_{1}}^{case3}<0$ $\&$ $X_{0_{2}}^{case3}<0$

We see that the above conditions put numerical bounds on $\xi_{3}$, $\xi_{4}$
and $\xi_{5}$ from the perspective of a collapsing scenario. The roots
$X_{0_{1,2}}^{case3}$ have been plotted against the parameters $\xi_{4}$ and
$\xi_{5}$ in Figs.(3) and (4).

Fig.3 Fig.4

Figs.3 and 4 show the variation of the collapse parameter $X_{0}$ for
different values of $\xi_{4}$ and $\xi_{5}$ for Case-3 of Model-1. Fig.3 shows
the variation for the first root $X_{0_{1}}^{case3}$, whereas Fig.4 shows the
variation for the second root $X_{0_{2}}^{case3}$. Here we have taken
$\xi_{3}=0.5$.

####  3.1.4 Case-4

Using Eq.(2.37) in Eq.(3.17) we get,

$\displaystyle\frac{2}{X_{0}}=\begin{array}[]{l}lim\\\ \tiny t\rightarrow 0\\\
\tiny r\rightarrow
0\end{array}\left[1+\frac{g_{1}}{2nr\left(\omega+3g_{2}\omega-
g_{2}\right)}\right]$ (3.28)

Just like Case-2, here also we get $\frac{2}{X_{0}}\rightarrow\infty$, which
implies $X_{0}\rightarrow 0$. The mass function being independent of $t$ and
$r$ does not generate a realistic collapsing scenario for this particular case
according to our scheme of study.

###  3.2 Model-2

####  3.2.1 Case-1

In this model we will replace $X_{0}$ by $Y_{0}$, just to differentiate the
results from those obtained for model-1. Moreover this is just a
representational issue. The definition remains same as given in Eq.(3.17).
Using Eq.(2.39) in Eq.(3.17) we get,

$\displaystyle\frac{2}{Y_{0}}=\begin{array}[]{l}lim\\\ \tiny t\rightarrow 0\\\
\tiny r\rightarrow
0\end{array}\left[1-\frac{h_{6}(t)}{r}AiryAi\left[-\frac{2^{1/3}nr}{g_{1}\left(-n/g_{1}\right)^{2/3}}\right]-\frac{h_{7}(t)}{r}AiryBi\left[-\frac{2^{1/3}nr}{g_{1}\left(-n/g_{1}\right)^{2/3}}\right]\right]$
(3.32)

Here we consider: $h_{6}(t)=\xi_{6}t$, $h_{7}(t)=\xi_{7}t$, where $\xi_{6}$
and $\xi_{7}$ arbitrary constants. Using these functional forms in Eq.(3.32)
we get the following algebraic equation for this case,

$\frac{1}{\Gamma(2/3)}\left(\frac{\xi_{6}}{3^{2/3}}+\frac{\xi_{7}}{3^{1/6}}\right)Y_{0}^{2}-Y_{0}+2=0$
(3.33)

It should be noted that this equation is similar to the one obtained for
Case-1 in model-1. This is due to the fact that, though the mass functions
have different forms in the two cases, yet their limiting values coincide with
other and hence generate similar collapsing scenario. The solution for the
above equation is obtained as,

$Y_{0_{1,2}}^{case1}=\frac{\sqrt{\Gamma(2/3)}}{2\sqrt{3}\left(\frac{\xi_{6}}{3^{2/3}}+\frac{\xi_{7}}{3^{1/6}}\right)}\left[\sqrt{3\Gamma(2/3)}\pm\sqrt{3\Gamma(2/3)-8\times
3^{1/3}\xi_{6}-8\times 3^{5/6}\xi_{7}}\right]$ (3.34)

where $3\Gamma(2/3)\geq 8\times 3^{1/3}\xi_{6}+8\times 3^{5/6}\xi_{7}$. The
collapsing outcomes may be discussed as below,

Conditions for a local NS: $Y_{0_{1}}^{case1}>0$ $\&$ $Y_{0_{2}}^{case1}<0$ OR
$Y_{0_{1}}^{case1}<0$ $\&$ $Y_{0_{2}}^{case1}>0$

Conditions for global NS: $Y_{0_{1}}^{case1}>0$ $\&$ $Y_{0_{2}}^{case1}>0$

Condition for BH: $Y_{0_{1}}^{case1}<0$ $\&$ $Y_{0_{2}}^{case1}<0$

The above conditions put numerical bounds on $\xi_{6}$ and $\xi_{7}$ from the
perspective of a collapsing scenario. The roots $Y_{0_{1,2}}^{case1}$ have
been plotted against the parameters $\xi_{6}$ and $\xi_{7}$ in Figs.(5) and
(6).

Fig.5 Fig.6

Figs.5 and 6 show the variation of the collapse parameter $Y_{0}$ for
different values of $\xi_{6}$ and $\xi_{7}$ for Case-1 of Model-2. Fig.5 shows
the variation for the first root $Y_{0_{1}}^{case1}$, whereas Fig.6 shows the
variation for the second root $Y_{0_{2}}^{case1}$.

####  3.2.2 Case-2

Using Eq.(2.41) in Eq.(3.17) we get,

$\displaystyle\frac{2}{Y_{0}}=\begin{array}[]{l}lim\\\ \tiny t\rightarrow 0\\\
\tiny r\rightarrow 0\end{array}\left[1+\frac{g_{1}}{2nr}\right]$ (3.38)

Here we get $\frac{2}{Y_{0}}\rightarrow\infty$, which implies
$Y_{0}\rightarrow 0$. The mass function being independent of $t$ and $r$ does
not generate a realistic collapsing scenario for this particular case
according to our scheme of study. This is equivalent to the scenarios in
Case-2 and Case-4 in model-1.

####  3.2.3 Case-3

Using Eq.(2.43) in Eq.(3.17) we get,

$\displaystyle\frac{2}{Y_{0}}=\begin{array}[]{l}lim\\\ \tiny t\rightarrow 0\\\
\tiny r\rightarrow 0\end{array}\left[1-\frac{h_{8}(t)}{r}-h_{9}(t)\right]$
(3.42)

Here we consider $h_{8}(t)=\xi_{8}t$ (where $\xi_{8}$ is an arbitrary
constant). We do not need to consider any particular functional form for
$h_{9}(t)$. This is because irrespective of the form of $h_{9}(t)$, it will
always yield a constant value in the limit $t\rightarrow 0$. This gives an
additional degree of freedom to the collapsing system. The above equation
yields,

$\xi_{8}Y_{0}^{2}+\left(\xi_{9}-1\right)Y_{0}+2=0$ (3.43)

where $\xi_{9}$ is the limiting value of $h_{9}(t)$ as $t\rightarrow 0$. The
above equations yields the solution,

$Y_{0_{1,2}}^{case3}=\frac{1-\xi_{9}\pm\sqrt{\left(\xi_{9}-1\right)^{2}-8\xi_{8}}}{2\xi_{8}}$
(3.44)

where $\left(\xi_{9}-1\right)^{2}\geq 8\xi_{8}$. The conditions for different
collapse outcomes are given below,

Conditions for a local NS: $Y_{0_{1}}^{case3}>0$ $\&$ $Y_{0_{2}}^{case3}<0$ OR
$Y_{0_{1}}^{case3}<0$ $\&$ $Y_{0_{2}}^{case3}>0$

Conditions for global NS: $Y_{0_{1}}^{case3}>0$ $\&$ $Y_{0_{2}}^{case3}>0$

Condition for BH: $Y_{0_{1}}^{case3}<0$ $\&$ $Y_{0_{2}}^{case3}<0$

The above conditions put numerical bounds on $\xi_{8}$ and $\xi_{9}$ from the
perspective of gravitational collapse of a massive star. The roots
$Y_{0_{1,2}}^{case3}$ have been plotted against the parameters $\xi_{8}$ and
$\xi_{9}$ in Figs.(7) and (8).

Fig.7 Fig.8

Figs.7 and 8 show the variation of the collapse parameter $Y_{0}$ for
different values of $\xi_{8}$ and $\xi_{9}$ for Case-3 of Model-2. Fig.7 shows
the variation for the first root $Y_{0_{1}}^{case3}$, whereas Fig.8 shows the
variation for the second root $Y_{0_{2}}^{case3}$.

####  3.2.4 Case-4

Using Eq.(2.45) in Eq.(3.17) we get,

$\displaystyle\frac{2}{Y_{0}}=\begin{array}[]{l}lim\\\ \tiny t\rightarrow 0\\\
\tiny r\rightarrow
0\end{array}\left[1-\frac{r^{2}}{6}-\frac{h_{10}(t)}{r}\right]$ (3.48)

Here we consider the functional form $h_{10}(t)=\xi_{10}t$, where $\xi_{10}$
is an arbitrary constant. From the above equation we get the algebraic
equation,

$\xi_{10}Y_{0}^{2}-Y_{0}+2=0$ (3.49)

Solving the above equation we get,

$Y_{0_{1,2}}^{case4}=\frac{1\pm\sqrt{1-8\xi_{10}}}{2\xi_{10}}$ (3.50)

where $\xi_{10}\leq 1/8$. Here the conditions for NS and BH can be discussed
as below,

Conditions for a local NS: $Y_{0_{1}}^{case4}>0$ $\&$ $Y_{0_{2}}^{case4}<0$ OR
$Y_{0_{1}}^{case4}<0$ $\&$ $Y_{0_{2}}^{case4}>0$

Conditions for global NS: $Y_{0_{1}}^{case4}>0$ $\&$ $Y_{0_{2}}^{case4}>0$

Condition for BH: $Y_{0_{1}}^{case4}<0$ $\&$ $Y_{0_{2}}^{case4}<0$

The above conditions put numerical bounds on $\xi_{10}$ from the perspective
of gravitational collapse of a massive star. The roots $Y_{0_{1,2}}^{case4}$
have been plotted against the parameter $\xi_{10}$ in Fig.(9).

Fig.9

Fig.9 shows the variation of the collapse parameter $Y_{0}$ for different
values of $\xi_{10}$ for Case-4 of Model-2.

###  3.3 Numerical Analysis

In order get greater insights about the nature of the singularity formed as an
end state of the gravitational collapse for our models, we have generated
plots for the collapsing parameter ($X_{0}$ or $Y_{0}$) against the other free
parameters. According to our scheme of study, we are interested only in the
signature of the collapse parameter ($X_{0}$ or $Y_{0}$) and not in the actual
value of the parameter. In figs.(1) and (2) we have generated plots for the
Case-1 (DP model) of Model-1. We see from the figures that the first root lies
in the negative region, whereas the second root lies in the positive region.
So it can predicted that for this case, we will have a local NS. In figs.(3)
and (4) plots have been obtained for Case-3 (PE model) of Model-1. We see that
for the considered initial conditions the first root is again negative and the
second root is positive. So here also the collapse will end in a local NS.
Figures (5) and (6) show the corresponding plots for Case-1 (TP model) of
Model-2. In fig.(5) although the major portion of the surface lies in the
negative region, yet there is an array of points represented by a straight
line lying in the positive region around $\xi_{7}=0$. In fig.(6) the entire
surface lies in the positive region. So in this case, there is a possibility
to get more than a local NS. By properly adjusting the initial conditions, it
is quite possible that more than one null geodesic originating in the
singularity reach a distant observer. In case of such an event, the NS will
become global in nature. In figs.(7) and (8) we have obtained plots for Case-3
(DPE model) of Model-2. Here we see that by properly adjusting the parameters
$\xi_{8}$ and $\xi_{9}$ both positive and negative values can be realized for
both the roots. So in this case we can have BH, local NS and global NS
depending on the chosen initial conditions. Finally in fig.(9) we have
obtained plots for Case-4 (TE model) of Model-2. Here we see that irrespective
of the initial conditions one root is always positive and the other is always
negative. So here the singularity is destined to be a local NS.

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

The gravitational strength of a singularity is defined as the estimate of its
destructive capacity. We know that most theories of gravity till date have
been plagued by the existence of singularities. Though theoretical methods of
removal of such singularities have been proposed in literature, yet they are
highly exotic in nature and far from being comprehensible. It is known that
singularities are holes in the fabric of the otherwise continuous and smooth
spacetime. Now for a weak singularity the hole is shallow and an extension of
space-time is possible through the singularity. This is equivalent to a
removable discontinuity mathematically and can be a cure for the discontinuity
of spacetime at a singularity. From the above discussion it is quite clear
that one should be highly interested in finding out whether a singularity is
strong or weak in nature. According to Tipler [58] a curvature singularity is
said to be strong if any object hitting it is squeezed to zero volume. In
Ref.[58] the condition for a strong singularity has been given as,

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

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

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

where

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

and

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

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

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

where

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

and $\dot{m_{0}}$ is given by the eqn.(4.18).

Studies by Dwivedi and Joshi in Refs.[60, 38] showed that any classical
singularity in Vaidya spacetime in Einstein gravity is supposed to be a strong
curvature singularity in a very strong sense. Additionally they have also
shown that the conjecture [61] that the strong curvature singularities are
never naked is not always true. It is speculated that in the background of
$f(R,T)$ gravity the strength of the singularity may weaken due to the exotic
component arising from the modified gravity. Moreover the structure of such a
NS was studied in detail in Ref.[62] and it was shown that the singularity
admits a directional behaviour in terms of curvature growth along the
geodesics terminating in the singularity. On the contrary it was found that in
a quantum regime the singularity formed is supposed to be gravitationally
weak, thus allowing a continuous extension of the spacetime beyond the
singularity [63]. Below we study the strength of the singularities for the
different models.

###  4.1 Model-1

####  4.1.1 Case-1

Using Eqs.(2.29), (4.10) and (4.18) we have,

$S=\lim_{\tau\rightarrow
0}\tau^{2}\psi=\frac{X_{0}^{2}}{2\Gamma(2/3)}\left(\frac{\xi_{1}}{3^{2/3}}+\frac{\xi_{2}}{3^{1/6}}\right)$
(4.24)

It is obvious that the signature of the above expression is independent of the
collapsing parameter $X_{0}$ since $X_{0}^{2}>0$. So the strength of the
singularity ultimately depends on the values of the parameters $\xi_{1}$ and
$\xi_{2}$. The condition for a strong singularity is
$\frac{1}{\Gamma(2/3)}\left(\frac{\xi_{1}}{3^{2/3}}+\frac{\xi_{2}}{3^{1/6}}\right)>0$
and that for a weak singularity is
$\frac{1}{\Gamma(2/3)}\left(\frac{\xi_{1}}{3^{2/3}}+\frac{\xi_{2}}{3^{1/6}}\right)\leq
0$.

####  4.1.2 Case-2

Using Eqs.(2.33), (4.10) and (4.18) we have,

$S=\lim_{\tau\rightarrow 0}\tau^{2}\psi=0$ (4.25)

The above value shows that the singularity formed is weak in nature. In the
previous section it was seen that for this model, the constancy of mass
parameter did not assist in studying the nature of the singularity. But
whatever be the nature of the singularity formed, it should always be
gravitationally weak in nature.

####  4.1.3 Case-3

Using Eqs.(2.35), (4.10) and (4.18) we have,

$S=\lim_{\tau\rightarrow
0}\tau^{2}\psi=\frac{X_{0}^{2}}{2\Gamma(2/3)}\left(\frac{\xi_{4}}{3^{2/3}}+\frac{\xi_{5}}{3^{1/6}}\right)$
(4.26)

Similar to case-1, here the strength of singularity is independent of $X_{0}$.
The condition for a strong singularity is
$\frac{1}{\Gamma(2/3)}\left(\frac{\xi_{4}}{3^{2/3}}+\frac{\xi_{5}}{3^{1/6}}\right)>0$
and that for a weak singularity is
$\frac{1}{\Gamma(2/3)}\left(\frac{\xi_{4}}{3^{2/3}}+\frac{\xi_{5}}{3^{1/6}}\right)\leq
0$.

####  4.1.4 Case-4

Using Eqs.(2.37), (4.10) and (4.18) we have,

$S=\lim_{\tau\rightarrow 0}\tau^{2}\psi=0$ (4.27)

This is a situation similar to case-2 where the singularity is always
gravitationally weak.

###  4.2 Model-2

####  4.2.1 Case-1

Using Eqs.(2.39), (4.10) and (4.18) we have,

$S=\lim_{\tau\rightarrow
0}\tau^{2}\psi=\frac{Y_{0}^{2}}{2\Gamma(2/3)}\left(\frac{\xi_{6}}{3^{2/3}}+\frac{\xi_{7}}{3^{1/6}}\right)$
(4.28)

Here the signature of the above expression and hence the strength of the
singularity depends on the values of the parameters $\xi_{6}$ and $\xi_{7}$.
We get a strong singularity if
$\frac{1}{\Gamma(2/3)}\left(\frac{\xi_{6}}{3^{2/3}}+\frac{\xi_{7}}{3^{1/6}}\right)>0$,
and a weak singularity if
$\frac{1}{\Gamma(2/3)}\left(\frac{\xi_{6}}{3^{2/3}}+\frac{\xi_{7}}{3^{1/6}}\right)\leq
0$.

####  4.2.2 Case-2

Using Eqs.(2.41), (4.10) and (4.18) we have,

$S=\lim_{\tau\rightarrow 0}\tau^{2}\psi=0$ (4.29)

Hence the singularity is gravitationally weak in nature.

####  4.2.3 Case-3

Using Eqs.(2.43), (4.10) and (4.18) we have,

$S=\lim_{\tau\rightarrow 0}\tau^{2}\psi=\frac{1}{2}Y_{0}^{2}\xi_{8}$ (4.30)

Here we have considered no special form for the function $h_{9}(t)$ as was
done in the previous section. It is clear from the above expression that the
strength of the singularity basically depends on the signature of $\xi_{8}$.
If $\xi_{8}>0$, then the singularity is strong and if $\xi_{8}\leq 0$, then
the singularity is weak in nature. However if we do consider a special form
for the function $h_{9}(t)$, we can have a different result. If we consider
$h_{9}(t)=\gamma_{9}\log(t)$, then we have from the Eqs.(2.43), (4.10) and
(4.18),

$S=\lim_{\tau\rightarrow
0}\tau^{2}\psi=\frac{1}{2}Y_{0}^{2}\left(\xi_{8}+\frac{\gamma_{9}}{Y_{0}}\right)$
(4.31)

Now using Eqs.(2.43), (4.18), (4.19) and (4.23) we get a relation from where
the values of $X_{0}$ may be extracted. Using these values of $X_{0}$ in the
above equation we may have a different scenario for the strength of the
singularity.

####  4.2.4 Case-4

Using Eqs.(2.45), (4.10) and (4.18) we have,

$S=\lim_{\tau\rightarrow 0}\tau^{2}\psi=\frac{1}{2}Y_{0}^{2}\xi_{10}$ (4.32)

Here the strength of the singularity depends on the signature of $\xi_{10}$.
If $\xi_{10}>0$, then the singularity is strong, and if $\xi_{10}\leq 0$, the
singularity is weak.

##  5 Conclusion and Discussion

In this work, we have explored a gravitational collapse mechanism of a massive
star in $f(R,T)$ gravity. A time dependent Vaidya spacetime is used to model
the collapsing phenomenon. The Einstien’s field equations for $f(R,T)$ gravity
in the Vaidya spacetime are calculated and the corresponding solutions for the
mass parameter $m(t,r)$ are obtained. We have considered two different
category of $f(R,T)$ models, each consisting of four sub-models. The two
models are considered on the basis of the nature of coupling between the
scalar invariants $R$ and $T$. The sub-models for each model basically involve
various combinations of power and exponential functional forms. Here we
considered the collapse of a massive star $(>20M_{\odot})$, which will
invariably continue its collapse until the formation of a singularity. The
huge mass of the parent star will always keep the collapsing mass beyond the
Chandrasekhar limit $(1.3M_{\odot})$, and hence neither the electron nor the
neutron degeneracy pressure will be able to counterbalance the inward the
collapsing force. Hence the collapse will not terminate in any middle stage
like a white dwarf or a neutron star, but will continue all the way to a
singularity (BH or NS). The scheme followed for the gravitational collapse
study involved the quest for outgoing radial null geodesics from the central
singularity formed as an end state of the collapse. If such outgoing geodesics
exists then the singularity becomes a naked singularity and the formation of
the event horizon is hindered. Such a situation will definitely defy the
cosmic censorship hypothesis. Moreover depending on the number of such
escaping geodesics, we can have a locally or globally naked singularity. More
number of escaping geodesics will mean greater exposure time of the
singularity to an external observer, and hence result in a globally naked
singularity. However if no such geodesic escape from the singularity, the
collapse is destined to end in a black hole and thus favour the censorship
hypothesis. Our study predicts that in almost all the cases of model-1 we get
a locally naked singularity. Model-2 seems to be a mixed bag, predicting the
formation of black holes, local and global naked singularities depending on
the initial conditions. However in the Case-4 of Model-2 (TE model), the
collapse always results in a local naked singularity. So here it should be
noted that the nature of coupling between the scalar invariants $R$ and $T$
does play a very important role in the nature of singularity formed as an end
state of the collapse. For minimal coupling (Model-1), we see that the
collapse generally ends in a local naked singularity. But for non-minimal
coupling (Model-2), all the options (BH, local and global NS) are possible
except the TE model (Case-4). Hence these models resulting from the minimal
and non-minimal coupling between curvature and matter, can be considered as
significant counterexamples of the cosmic censorship hypothesis. But as we
know that non-minimal coupling is observationally the favoured model, the
result derived for this model-2 will be cosmologically more relevant. Moreover
we see that for minimal coupling we generally do not get the global nature of
the naked singularity, but in case of non-minimal coupling this can be a
reality. One thing which may be worrying for the reader is that the final
limiting forms ($X_{0}$ or $Y_{0}$) in the collapsing scheme does not involve
the model parameters $g_{i}$ or $\beta_{i}$, $i=1,2,3$. So how does one
differentiate the collapse outcomes between the models? We see that here the
solutions are in terms of Airy functions which are relatively complicated
mathematical forms. In the limiting scenario the argument of these functions
vanish giving constant values, which is reason we do not see any model
parameters in the limiting forms. However it should be mentioned here that the
functional forms of $X_{0}$ or $Y_{0}$ are different for different models,
which is testimony of the fact they arise from different functional forms.
Moreover the imprints of such functional forms are carried by the functions
$h_{i}(t)$, $i=1,2,3...10$, and thus the parameters $\xi_{i}$, $i=1,2,3...10$
which are present in the analysis. One more thing that the reader needs to
note is that some of the solutions derived for the special cases in model-2
are valid for early universe. So in such cases we are actually studying the
collapsing scenarios of primordial black holes that existed at the beginning
of the universe. However it should also be kept in mind that these solutions
are just specific examples to get greater insights into the bigger picture and
in no sense represent the entire story.

To complement the collapsing scheme we have studied the strength of the
singularity formed for all our models. We see that for the DP (case-1) and PE
(case-3) models of model-1, the singularity can be both gravitationally weak
or strong depending on the model parameters. However for the DE (case-2) and
EP (case-4) models, the singularity formed is always gravitationally weak. A
weak singularity will obviously be pathologically favoured because the
spacetime can be extended beyond such a singularity and we get a sense of
continuity. For model-2, we see that for the TP (case1), DPE (case-3) and TE
(case-4) models the strength of the singularity depends on the initial
conditions but for the DEP (case-2) model the singularity is always
gravitationally weak. So it is understandable that for all models, by suitably
adjusting the initial conditions, we can have a sufficiently weak singularity,
which will be cosmologically desirable, since an extension of the spacetime
beyond the singularity becomes a possibility. In principle we can create a
scenario where the singularity may be completely avoided. This is a direct
consequence of the coupling of matter with geometry and hence an intrinsic
property of $f(R,T)$ models and their exotic nature.

## Acknowledgments

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

##  6 Appendix

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

Model-1

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

$r^{4}\left\\{f_{1}(R)+f_{2}(T)-2\left(f_{2}^{\prime}(T)+1\right)\left(\rho+\sigma\right)+2f_{2}^{\prime}(T)\omega\rho\right\\}-8f_{1}^{\prime\prime\prime}(R)\dot{m}^{\prime
2}-r\left\\{8f_{1}^{\prime\prime\prime}(R)\dot{m}^{\prime}\dot{m}^{\prime\prime}\right.$
$\left.+r\left(f_{1}^{\prime}(R)rm^{\prime\prime}+m\left(r\left(f_{1}(R)+f_{2}(T)-2\rho\left(1+f_{2}^{\prime}(T)-\omega
f_{2}^{\prime}(T)\right)\right)-f_{1}^{\prime}(R)m^{\prime\prime}\right)-2f_{1}^{\prime}(R)\dot{m}\right.\right.$
$\left.\left.+2f_{1}^{\prime\prime\prime}(R)\dot{m}^{\prime\prime
2}+4f_{1}^{\prime\prime}(R)\ddot{m}^{\prime}+2f_{1}^{\prime\prime}(R)r\ddot{m}^{\prime\prime}\right)\right\\}=0$
(6.1)

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

$f_{1}^{\prime\prime\prime}(R)\left[-4m^{\prime}+r\left(m^{\prime\prime}+rm^{(3)}\right)\right]^{2}+r^{2}f_{1}^{\prime\prime}(R)\left(12m^{\prime}-6rm^{\prime\prime}+r^{3}m^{(4)}\right)=0$
(6.2)

where $(3)$ and $(4)$ in the power represents the third and fourth order
derivative with respect to $r$ respectively.

Model-2

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

$2r^{4}\sigma+2r^{3}\rho\left(r-m\right)-f_{2}(R)r^{3}\left(f_{3}(T)+2f_{3}^{\prime}(T)\omega\rho\right)\left(r-m\right)-f_{1}(R)r^{3}\left(r-m\right)+2f_{2}(R)f_{3}^{\prime}(T)r^{3}\left\\{r\left(\rho+\sigma\right)-\rho
m\right\\}$
$+\left(f_{1}^{\prime}(R)+f_{2}^{\prime}(R)f_{3}(T)\right)r^{2}\left\\{m^{\prime\prime}\left(r-m\right)-2\dot{m}\right\\}+2f_{2}^{\prime}(R)f_{3}^{\prime\prime}(T)r^{4}\lambda^{2}\dot{\rho}^{2}+4f_{2}^{\prime\prime}(R)f_{3}^{\prime}(T)r^{2}\lambda\dot{\rho}\left(2\dot{m}^{\prime}+r\dot{m}^{\prime\prime}\right)$

$+2f_{1}^{\prime\prime\prime}(R)\left(2\dot{m}^{\prime}+r\dot{m}^{\prime\prime}\right)^{2}+2f_{2}^{\prime\prime\prime}(R)f_{3}(T)\left(2\dot{m}^{\prime}+r\dot{m}^{\prime\prime}\right)^{2}+2f_{2}^{\prime}(R)f_{3}^{\prime}(T)r^{4}\lambda\ddot{\rho}+4f_{1}^{\prime\prime}(R)r^{2}\ddot{m}^{\prime}+4f_{2}^{\prime\prime}(R)f_{3}(T)r^{2}\ddot{m}^{\prime}$
$+2\left(f_{1}^{\prime\prime}(R)+f_{2}^{\prime\prime}(R)f_{3}(T)\right)r^{3}\ddot{m}^{\prime\prime}=0$
(6.3)

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

$r^{6}\left(-f_{2}^{\prime}(R)f_{3}^{\prime\prime}(T)\lambda^{2}(\rho^{\prime})^{2}-f_{2}^{\prime}(R)f_{3}^{\prime}(T)\lambda\rho^{\prime\prime}\right)-r^{3}\left[2f_{2}^{\prime\prime}(R)f_{3}^{\prime}(T)\lambda\rho^{\prime}\left(-4m^{\prime}+r\left(m^{\prime\prime}+rm^{\prime\prime\prime}\right)\right)\right]$
$-f_{1}^{\prime\prime\prime}(R)\left[-4m^{\prime}+r\left(m^{\prime\prime}+rm^{\prime\prime\prime}\right)\right]^{2}-f_{2}^{\prime\prime\prime}(R)f_{3}(T)\left[-4m^{\prime}+r\left(m^{\prime\prime}+rm^{\prime\prime\prime}\right)\right]^{2}-r^{2}f_{1}^{\prime\prime}(R)\left(12m^{\prime}-6rm^{\prime\prime}+r^{3}m^{iv}\right)$
$-r^{2}f_{2}^{\prime\prime}(R)f_{3}(T)\left(12m^{\prime}-6rm^{\prime\prime}+r^{3}m^{iv}\right)=0$
(6.4)

Airy Function[64]

Here we would like to present a short description of Airy functions for the
reader’s convenience. Airy function is a special function named after the
British astronomer George Biddell Airy (1801-1892). There are in fact, two
Airy functions $Ai(x)$ (Airy function of the first kind) and $Bi(x)$ (Airy
function of the second kind), which are linearly independent solutions of the
Airy differential equation given by,

$\frac{d^{2}y}{dx^{2}}-xy=0$ (6.5)

For real values of $x$ the Airy function of the first kind is defined by the
improper integral,

$Ai(x)=\frac{1}{\pi}\int_{0}^{\infty}\cos\left(\frac{t^{3}}{3}+xt\right)dt\equiv\frac{1}{\pi}\lim_{b\rightarrow\infty}\int_{0}^{b}\cos\left(\frac{t^{3}}{3}+xt\right)dt$
(6.6)

which is convergent. This solution is subject to the condition $y\rightarrow
0$ as $x\rightarrow\infty$. The Airy function of the second kind is defined
as,

$Bi(x)=\frac{1}{\pi}\int_{0}^{\infty}\left[exp\left(-\frac{t^{3}}{3}+xt\right)+\sin\left(\frac{t^{3}}{3}+xt\right)\right]dt$
(6.7)

This solution has the same amplitude of oscillation as Ai(x) as
$x\rightarrow-\infty$ differing in phase by $\pi/2$. The values of Airy
function $(Ai(x),Bi(x))$ and its derivatives $(Ai^{\prime}(x),Bi^{\prime}(x))$
at $x=0$ are given by,

$Ai(0)=\frac{1}{3^{2/3}\Gamma(2/3)},~{}~{}~{}~{}Bi(0)=\frac{1}{3^{1/6}\Gamma(2/3)},~{}~{}~{}Ai^{\prime}(0)=-\frac{1}{3^{1/3}\Gamma(1/3)},~{}~{}~{}Bi^{\prime}(0)=\frac{3^{1/6}}{\Gamma(1/3)}$
(6.8)

## References

  * [1] S. Perlmutter et. al. :- Astrophys. J. 517 565 (1999).
  * [2] A. G. Riess et al. :- Astron. J. 116 1009 (1998).
  * [3] S. Nojiri, S. D. Odintsov, V. K. Oikonomou :- Phys. Rep. 692 1 (2017).
  * [4] S. Nojiri and S. D. Odintsov :- Int. J. Geom. Methods Mod. Phys. 04 115 (2007).
  * [5] S. Capozziello, R. D Agostino, O. Luongo:- Int. J. Mod. Phys. D 28 1930016 (2019).
  * [6] P. Brax :- Rep. Prog. Phys. 81 016902 (2018)
  * [7] T. P. Sotiriou, V. Faraoni :- Rev. Mod. Phys. 82 451 (2010).
  * [8] A. De Felice, S. Tsujikawa :- Living Rev. Relativity 13 3 (2010).
  * [9] T. Harko, F. S. N. Lobo :- Eur. Phys. J. C. 70 373 (2010).
  * [10] R. Ribeiro, J. P ramos :- Phys. Rev. D 90 124065 (2014).
  * [11] R. P. L. Azevedo, J. P ramos :- Phys. Rev. D 94 064036 (2016).
  * [12] B. Pourhassan, P. Rudra :- Phys. Rev. D 101 084057 (2020).
  * [13] T. Harko, F. S. N. Lobo, S. Nojiri, S. D. Odintsov :- Phys. Rev. D. 84 024020 (2011).
  * [14] M. Sharif, M. Zubair :- JCAP 03 028 (2012).
  * [15] E. H. Baffou, M. J. S. Houndjo, M. E. Rodrigues, A. V. Kpadonou, J. Tossa :- Phys. Rev. D 92 8, 084043 (2015).
  * [16] P. Rudra :- Eur. Phys. J. Plus 130 4, 66 (2015).
  * [17] H. Shabani, M. Farhoudi :- P͡hys. Rev. D 88 044048 (2013).
  * [18] F. G. Alvarenga, A. de la Cruz-Dombriz, M. J. S. Houndjo, M. E. Rodrigues, D. S ez-G mez :- Phys. Rev. D 87 10, 103526 (2013).
  * [19] A. Das, S. Ghosh, B. K. Guha, S. Das, F. Rahaman :- Phys. Rev. D 95 12, 124011 (2017).
  * [20] R. Zaregonbadi, M. Farhoudi, N. Riazi :- Phys. Rev. D 94 084052 (2016).
  * [21] M. Sharif, A. Siddiqa :- Gen. Rel. Grav. 51 6, 74 (2019).
  * [22] S. Hansraj, A. Banerjee :- Phys. Rev. D 97 10, 104020 (2018).
  * [23] J. R. Oppenhiemer, H. Snyder :- Phys. Rev. 56 455 (1939).
  * [24] R. C. Tolman :- Proc. Natl. Acad. Sci. USA 20 169 (1934);
  * [25] H. Bondi :- Mon. Not. Astron. Soc. 107 410 (1947).
  * [26] P. S. Joshi, D. Malafarina :- Int. J. Mod. Phys. D 20 14, 2641 (2011).
  * [27] D. Malafarina :- Universe 3 48 (2017).
  * [28] R. Penrose :- Riv. Nuovo Cimento 1 252 (1969).
  * [29] P. Chen, Y. C. Ong, D-h. Yeom :- Phys. Rept. 603 1 (2015).
  * [30] B. Zhang, Q-y. Cai, M-s. Zhan, L. You :- Phys. Rev. D 87 4, 044006 (2013).
  * [31] L. Smolin :- Phys. Rev. D 90 2, 024074 (2014).
  * [32] H. Nikolic :- Phys. Lett. B 678 218 (2009).
  * [33] D. M. Eardley, L. Smarr :- Phys. Rev. D 19 2239 (1979).
  * [34] D. Christodoulou :- Commun. Math. Phys. 93 171 (1984).
  * [35] R. P. A. C. Newman :- Class. Quantum Grav. 3 527 (1986).
  * [36] I. H. Dwivedi, P. S. Joshi :- Class. Quantum Grav. 9 L39 (1992).
  * [37] P. S. Joshi, I. H. Dwivedi :- Phys. Rev. D 47 5357 (1993).
  * [38] P. S. Joshi, I.H. Dwivedi :- Commun. Math. Phys. 146 333 (1992).
  * [39] P. S. Joshi, T.P. Singh :- Phys. Rev. D 51 6778 (1995).
  * [40] B. Waugh, K. Lake :- Phys. Rev. D 34 2978 (1986).
  * [41] A. Ori, T. Piran :- Phys. Rev. D 42 1068 (1990).
  * [42] K. Lake :- Phys. Rev. Lett. 68 3129 (1992).
  * [43] P. Szekeres, V. Iyer :- Phys. Rev. D 47 4362 (1993).
  * [44] S. G. Ghosh, S. D. Maharaj :- Phys. Rev. D 85 124064 (2012).
  * [45] P. C. Vaidya :- Proc. Indian Acad. Sci. Sect. A 33 264 (1951).
  * [46] P. Rudra, M. Faizal, A. F. Ali :- Nucl. Phys. B 909 725 (2016).
  * [47] Y. Heydarzade, P. Rudra, F. Darabi, A. F. Ali, M. Faizal :- Phys. Lett. B 774 46 (2017).
  * [48] P. Rudra, S. Maity :- Eur. Phys. J. C 78 828 (2018).
  * [49] Y. Heydarzade, P. Rudra, B. Pourhassan, M. Faizal, A. F. Ali :- J. Cosmol. Astropart. Phys. 06 038 (2018).
  * [50] P. Rudra, U. Debnath :- Can. J. Phys. 92(11) 1474 (2014).
  * [51] P. Rudra, R. Biswas, U. Debnath :- Astrophys. Space Sci. 354 2 (2014).
  * [52] U. Debnath, P. Rudra, R. Biswas :- Astrophys. Space Sci. 339 135 (2012).
  * [53] P. Rudra, R. Biswas, U. Debnath :- Astrophys.Space Sci. 335 505 (2011).
  * [54] P. Rudra :- Nucl. Phys. B. 956 115014 (2020).
  * [55] L. D. Landau, E. M. Lifshitz :- The Classical Theory of Fields (Butterworth-Heinemann, Oxford, 1998).
  * [56] J. de Boer, J. Hartong, N. A. Obers, W. Sybesma, S. Vandoren :- SciPost Phys. 5 003 (2018)
  * [57] T. P. Singh, P.S. Joshi :- Class. Quant. Gravity 13 559 (1996).
  * [58] F. J. Tipler :- Phys. Lett. A. 64, 8 (1977)
  * [59] M. D. Mkenyeleye, R. Goswami, and S. D. Maharaj :- Phys. Rev. D 90 064034 (2014).
  * [60] I. H. Dwivedi, P. S. Joshi :- Class. Quantum Grav. 6 1599 (1989).
  * [61] F. J. Tipler, C. J. S. Clarke, G. F. R. Ellis :- General Relativity and Gravitation vol 2 p97 ed A Held (NewYork Plenum) (1980).
  * [62] I. H. Dwivedi, P. S. Joshi :- Class. Quantum Grav. 8 1339 (1991).
  * [63] A. Bonanno, B. Koch, A. Platania :- Foundations of Phys. 48 1393 (2018).
  * [64] $https://en.wikipedia.org/wiki/Airy_{f}unction$