# Gravitational Collapse without Singularity Formation in Brans-Dicke Theory

A. H. Ziaie<EMAIL_ADDRESS>Research Institute for Astronomy and
Astrophysics of Maragha (RIAAM), University of Maragheh, P. O. Box 55136-553,
Maragheh, Iran  H. Shabani<EMAIL_ADDRESS>Physics Department,
Faculty of Sciences, University of Sistan and Baluchestan, Zahedan, Iran  H.
Moradpour<EMAIL_ADDRESS>Research Institute for Astronomy and
Astrophysics of Maragha (RIAAM), University of Maragheh, P. O. Box 55136-553,
Maragheh, Iran

###### Abstract

In the present work we study collapse process of a homogeneous and isotropic
fluid in Brans-Dicke (BD) theory with non-vanishing spacetime torsion. In this
theory, torsion can be generated by the BD scalar field as well as the
intrinsic angular momentum (spin) of matter. Assuming the matter content of
the collapsing body to be a Weyssenhoff fluid, which is a generalization of
perfect fluid in general relativity (GR) in order to include the spin effects,
we find that in BD theory with torsion, the existence of spin effects could
avoid the spacetime singularity that forms in the original version of this
theory [1]. Numerical simulations of collapse model show that the spacetime
singularity is replaced by a non-singular bounce, the spacetime event at which
the collapse process halts at a minimum radius and then turns into an
expanding phase. Moreover, the model parameters can be set so that the
apparent horizon will never meet the boundary of the collapsing body so that
the bounce event can be detectable by external observers in the Universe.

##  1 Introduction

The process of gravitational collapse in GR has been studied for many years
since the pioneering work by Oppenheimer, Snyder, and Datt (OSD) [2] showed
that collapse of a spherical matter distribution under its own gravity leads
to the formation of a black hole (BH). In this simple study where the
collapsing matter is described by a homogeneous pressure-less fluid, a
spacetime singularity forms as the collapse end state which is necessarily
hidden by an event horizon and therefore invisible to far away observers.
Since then, several models have been considered in the literature with the aim
of understanding possible extensions of OSD collapse scenario. In the context
of GR, singularity theorems predict that under physical circumstances
(requirements for the matter content i.e., energy conditions) the final fate
of gravitational collapse scenario would be the formation of a spacetime
singularity, an event at which densities as well as curvatures grow
unboundedly and become infinite [3]. It is generally believed that in these
extreme spacetime events the classical framework of GR theory breaks down [4,
5]. Unlike the BH case, further investigations showed that for certain matter
distributions that satisfy physically reasonable conditions the spacetime
singularity that arises as the collapse final outcome is not hidden behind a
horizon and thus can causally connect to far away observers, see, e.g., [6]
and references therein. These spacetime events are called naked singularities
that the existence of which has been also reported in modified gravity
theories [7]. Although GR has proven to be very successful theory in
describing gravitational phenomena, the occurrence of spacetime singularities
signal that this theory is utilized beyond its domain of validity. Therefore,
it is widely believed that in the late stages of a collapse scenario where
super dense regimes of extreme gravity are present, GR must be replaced by a
suitable modified gravity theory in order to cure the singularity problem.
This is because the formation of a spacetime singularity can cause issues such
as future predictability and path incompleteness [8]. During the past decades,
much effort has gone into searching for physically reasonable collapse
settings being free of spacetime singularities. The outcomes of research
indicate that in comparison to singular spacetimes in GR, the presence of
additional/correction terms within the classical framework of modified gravity
theories can lead to singularity avoidance. Such non-singular scenarios have
been reported in $f(R)$ gravity [9] where it is shown that adding higher order
contributions of Ricci scalar could cure the curvature singularity in the
collapse process. In the context of Einstein-Cartan (EC) theory [10, 11] it is
shown that ingredients mimicking spin effects (associated with fermions) can
replace the spacetime singularity by a non-singular bounce. Also, in this
theory, nonsingular cosmological models have been studied in [12] and in [13]
it is shown that a non-minimal coupling between gravity and fermions can
influence the endstate of a classical collapse scenario long before quantum
gravity effects become important.

Nevertheless, it is a general belief that in very late stages of a collapse
process where super-ultra-dense regions of extreme gravity are present, the
Planck scale physics become dominant and it is the nontrivial quantum gravity
effects that finally remove the classical singularity [14, 15]. Much attempts
have been made so far in this direction in order to find a consistent model
which includes quantum effects in a collapse setting. In this regard, though a
full quantum theory of gravity has not yet been developed, concerted efforts
have been made to take into account the quantum aspects of gravity in the
regions near the singularity in terms of semi-classical effects, either by
quantizing the matter sector in the field equations of GR or by quantizing
certain limited degrees of freedom of the spacetime metric [5]. In [16],
quantum effects in gravitational collapse of a scalar field which classically
ends in formation of a naked singularity has been considered and it is shown
that these effects lead to the dissolution of the collapsing cloud before the
singularity can form. Singularity removal in inhomogeneous scalar field
collapse with quantum gravity corrections has been presented in [17]. The
possibility of singularity avoidance in a homogeneous collapse process with
modified matter source inspired by the quantum gravity effects has been
discussed in [18]. See also [19] for other nonsingular collapse scenarios and
[20] for nonsingular quantum cosmological models.

Among the modified gravity theories, the BD theory is one of the simplest and
well studied generalizations of GR where the gravitational coupling constant
is replaced by a varying scalar field [21]. This theory in which,
gravitational interactions are described by two fundamental non-matter fields,
i.e., the spacetime metric $g_{\mu\nu}$ and a non-minimally coupled scalar
field $\phi$, is an attempt towards incorporating Mach’s principle into GR
framework [21]. The BD scalar field, in the Jordan representation [22],
couples to gravity with an adjustable parameter, $\omega$, therefore acting as
a mediator between matter fields and spacetime geometry. Since its advent, the
BD theory has attracted attention of many researchers as it arises naturally
from supergravity models, string theories at low energies and dimensional
reduction of Kaluza-Klein theories [23]. In addition to its cosmological
implications [22],[24], in the past years, the issue of gravitational collapse
and its final outcome in BD theory has been studied comprehensively under
various models among which we can quote: collapse of an ideal gas in BD theory
[25], investigation of the waveform and amplitude of scalar-type gravitational
waves in spherically symmetric dust fluid collapse [26], perfect fluid as well
as null fluid collapse and formation of naked singularity [27, 28], collapse
of a self-gravitating electrically charged scalar field [29], dynamics of
string-inspired charged black holes [30] and critical behavior in BD collapse
[31]. Numerical simulation of Oppenheimer-Snyder collapse in BD theory has
been performed in [1] where the authors concluded that black holes that arise
as the end product of the collapse scenario are identical to those of GR in
final equilibrium, but are completely different from those of GR during
dynamical evolution.

Recently, a modified version of BD theory with non-vanishing torsion has been
investigated in [32] and it has been shown that in addition to the spin of
matter, the BD scalar field itself has contribution in generating the torsion
field. This theory can be viewed as somehow a unification of BD and EC
theories that we call it hereafter ECBD theory. Cosmological as well as
astrophysical aspects of ECBD theory has been investigated in [33, 34] and
higher dimensional extension of this theory has been studied in [35]. In [36],
the authors have studied the existence and stability of Einstein static
universe, a cosmological model according to which, our universe admits no
timelike singularity, it is ever existing and lives in an almost a static
state in the infinite past ($t\rightarrow-\infty$) and then evolves into an
inflationary era. In this model, the Universe emerges from the Einstein static
state rather than the initial big bang singularity. As it was mentioned
earlier, gravitational collapse in BD theory results in formation of spacetime
singularities either naked or covered by a horizon. It is therefore well
motivated to search for possible nonsingular collapse scenarios that could
emerge from ECBD theory in which the spin effects prevent singularity
formation. Our paper is then organized as follows: In section 2 we give a
brief review on BD theory and present a singular solution for the collapse
process in this theory. In section 3 we study the evolution of collapse
scenarios in ECBD theory and seek for nonsingular solutions along with
investigating the dynamics of apparent horizon. Sec. 4 is devoted to exterior
spacetime and we summarize our results in section 5.

##  2 Gravitational Collapse in BD theory

In the present section we give a brief review on collapse scenario in BD
theory without considering the effects of spin and torsion. In Jordan frame
[22], this theory is described by the following action

$S=\int d^{4}x\sqrt{-g}\left[\frac{\phi
R}{\kappa^{2}}-\frac{\omega}{\phi}g^{\alpha\beta}\nabla_{\alpha}\phi\nabla_{\beta}\phi\right]+S^{m},$
(1)

where $\nabla_{\alpha}$ denotes covariant derivative with respect to the
Christoffel connection $\Gamma^{\alpha}_{\,\,\beta\gamma}$, $S^{m}$ is the
action describing ordinary matter (i.e., any form of physical matter different
from the BD scalar field $\phi$), $\omega$ is a dimensionless parameter
representing the strength of coupling between the scalar field and the
spacetime metric and $\kappa^{2}=16\pi/c^{4}$. Varying the above action with
respect to $\phi$ and $g_{\mu\nu}$ gives the field equations in BD theory as

$\displaystyle
R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R=\frac{\kappa^{2}}{\phi}T_{\mu\nu}+\frac{\omega}{\phi^{2}}\left[\nabla_{\mu}\phi\nabla_{\nu}\phi-\frac{1}{2}g_{\mu\nu}\nabla_{\eta}\phi\nabla^{\eta}\phi\right]+\frac{1}{\phi}(\nabla_{\mu}\nabla_{\nu}\phi-
g_{\mu\nu}\Box\phi),$ (2)
$\displaystyle\Box\phi=\frac{\kappa^{2}}{2\omega+3}T^{\alpha}_{\alpha},$ (3)

where $T_{\mu\nu}$ is the energy momentum tensor (EMT) of matter fields which
fulfills the conservation equation $T^{\mu\nu}_{;\mu}=0$. We assume an
isotropic and homogeneous spacetime for the collapsing body whose matter
content obeys a perfect fluid with the EMT
$T_{\mu\nu}=pg_{\mu\nu}+(\rho+p)u_{\mu}u_{\nu}$, where $\rho$ and $p$ are
energy density and pressure of the fluid and $u_{\mu}$ is its four-vector
velocity. A linear equation of state (EoS) $p=w\rho$ is also assumed with
$-1\leq w\leq 1$. The interior spacetime can be parametrized by a spatially
non-flat Friedmann-Lemaître-Robertson-Walker (FLRW) line element, given by

$\displaystyle
ds^{2}=-dt^{2}+\frac{a^{2}(t)dr^{2}}{1-kr^{2}}+{\mathcal{R}}^{2}(t,r)d\Omega^{2},$
(4)

where $k$ denotes the spatial curvature, ${\mathcal{R}}(t,r)=ra(t)$ is the
physical radius of collapsing object and $d\Omega^{2}$ is the metric on a unit
two-sphere. The field equations (2) and (3) can then be written as (we use the
units so that $\kappa^{2}=1$)

$\displaystyle
3H^{2}+\frac{3k}{a^{2}}+3H\frac{\dot{\phi}}{\phi}-\frac{\omega}{2}\left(\frac{\dot{\phi}}{\phi}\right)^{2}=\frac{\rho}{\phi},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
H=\frac{d{a}/dt}{a}=\frac{\dot{a}}{a},$ (5) $\displaystyle
2\dot{H}+3H^{2}+\frac{k}{a^{2}}+2H\frac{\dot{\phi}}{\phi}+\frac{\omega}{2}\left(\frac{\dot{\phi}}{\phi}\right)^{2}+\frac{\ddot{\phi}}{\phi}=-\frac{p}{\phi},$
(6)
$\displaystyle\ddot{\phi}+3H\dot{\phi}=\frac{\rho-3p}{2\omega+3},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
H=\frac{da/dt}{a}=\frac{\dot{a}}{a}.$ (7)

Also, the conservation equation for EMT reads

$\frac{d\rho}{dt}+3H(\rho+p)=0\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \Longrightarrow\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\rho(a)=\rho_{i}a^{-3(1+w)},$ (8)

where $\rho_{i}$ is the initial energy density at the onset of collapse
process. The set of differential equations (5)-(7) are not independent from
each other, i.e., differentiating Eq. (5) with respect to time along with
substituting for $\ddot{a}$, $\dot{\rho}$ and the term $\dot{a}^{2}+k$ from
Eqs. (6), (8) and (5) leaves us with the evolution equation for BD scalar
field, Eq. (7). This means that only the two equations (6) and (7) are
independent. We therefore proceed to find an exact solution for these set of
differential equations. To this aim we change the differentiation variable
from $t$ to $a$, this gives

$\displaystyle
2aHH_{,a}+3H^{2}+\frac{k}{a^{2}}+\frac{\rho_{i}}{\phi}\left(w+\frac{1-3w}{2\omega+3}\right)a^{-3(1+w)}-aH^{2}\frac{\phi_{,a}}{\phi}+\frac{\omega}{2}a^{2}H^{2}\left(\frac{\phi_{,a}}{\phi}\right)^{2}=0,$
(9)
$\displaystyle\left(4aH^{2}+a^{2}HH_{,a}\right)\phi_{,a}+a^{2}H^{2}\phi_{,aa}-\rho_{i}\frac{1-3w}{2\omega+3}a^{-3(1+w)}=0.$
(10)

Now, taking the functionalities of $H(a)$ and $\phi(a)$ as

$H(a)=\alpha a^{\gamma},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \phi(a)=\beta a^{\delta},$ (11)

an exact solution can be obtained by setting

$\displaystyle\gamma=-1,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \delta=-(1+3w),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\alpha=-\frac{\sqrt{2k}}{\left[(3w-2)w\omega-(\omega+2)\right]^{\frac{1}{2}}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\beta=\frac{2+\omega(1+w(2-3w))}{2k(1+3w)(3+2\omega)}\rho_{i},$ (12)

provided that

$\left(-1\leq
w<-\frac{1}{3}\land\omega>\frac{2}{3w^{2}-2w-1}\right)\lor\left(-\frac{1}{3}<w<1\land\omega<\frac{2}{3w^{2}-2w-1}\right).$
(13)

The above condition guarantees that for positive spatial curvature,
$\alpha\in\mathbb{R}^{-}$. We therefore get the scale factor as

$a(t)=a_{i}+\alpha(t-t_{i}),\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ t_{s}=t_{i}-\frac{a_{i}}{\alpha}.$ (14)

The above solution indicates a spacetime singularity that forms at the finite
amount of time $t_{s}$. At this spacetime event, the curvature invariants such
as Kretschmann scalar

${\tt
K}=\frac{12}{a^{4}}\left[a^{2}\ddot{a}^{2}+\dot{a}^{4}+2k\dot{a}^{2}+k^{2}\right],$
(15)

and energy density (8) grow unboundedly and diverge. Hence it is a true
curvature singularity that cannot be removed by any coordinate transformation.
We note that this singular solution stands for a closed FLRW geometry ($k>0$)
and the other singular spacetimes with $k=0$ have been previously addressed in
the literature, see e.g. [27]. Our aim in the next section is to seek for
possible non-singular scenarios in modified BD theory in the presence of
spacetime torsion.

##  3 Non-Singular Collapse Scenario in ECBD theory

Consider a Riemann-Cartan manifold as the background spacetime in which the
general connection $\tilde{\Gamma}^{\alpha}_{\,\beta\mu}$ is defined. The
spacetime torsion is given by the antisymmetric part of this connection as

$T^{\alpha}_{\,\mu\nu}=\tilde{\Gamma}^{\alpha}_{\,\mu\nu}-\tilde{\Gamma}^{\alpha}_{\,\nu\mu}.$
(16)

Assuming that the connection is metric compatible, i.e.,
$\tilde{\nabla}_{\mu}g_{\alpha\beta}=0$, where $\tilde{\nabla}_{\mu}$ is the
covariant derivative with respect to $\tilde{\Gamma}^{\alpha}_{\,\beta\mu}$,
we get

$\tilde{\Gamma}^{\alpha}_{\,\beta\mu}=\Gamma^{\alpha}_{\,\beta\mu}+K^{\alpha}_{\,\beta\mu},$
(17)

where $\Gamma^{\alpha}_{\,\beta\mu}$ is the Christoffel connection and
$K^{\alpha}_{\,\beta\mu}$ is the contorsion tensor defined as

$K^{\alpha}_{\,\beta\mu}=\frac{1}{2}\left[T^{\alpha}_{\,\beta\mu}-T^{\,\,\alpha}_{\beta\,\,\mu}-T^{\,\,\beta}_{\alpha\,\,\mu}\right].$
(18)

The action in BD theory with torsion is given as [33, 34]

${\cal
S}=\int\sqrt{-g}d^{4}x\left[\frac{\phi\tilde{R}}{\kappa^{2}}-\frac{\omega}{\phi}g_{\alpha\beta}\tilde{\nabla}^{\alpha}\phi\tilde{\nabla}^{\beta}\phi+L_{m}(g_{\alpha\beta},K^{\alpha}_{\,\beta\gamma},\psi^{i})\right],$
(19)

where $L_{m}$ is the Lagrangian for mater field(s) $\psi^{i}$ [37]. Similar to
EC theory, the Ricci curvature scalar is constructed out of the general
connection and is given by

$\displaystyle\tilde{{R}}={R}(\Gamma)+\nabla_{\lambda}{K}^{\lambda\rho}\\!\\!\leavevmode\nobreak\
_{\rho}-\nabla_{\rho}{K}^{\lambda\rho}\\!\\!\leavevmode\nobreak\
_{\lambda}+{K}^{\sigma\mu}\\!\\!\leavevmode\nobreak\
_{\mu}{K}^{\lambda}\\!\\!\leavevmode\nobreak\
_{\sigma\lambda}-{K}^{\sigma\rho}\\!\\!\leavevmode\nobreak\
_{\nu}{K}^{\nu}\\!\\!\leavevmode\nobreak\ _{\sigma\rho}.$ (20)

We note that in addition to spin degrees of freedom of matter which is a
source of spacetime torsion, due to nonminial coupling between BD scalar field
and curvature, the BD field itself acts as a source of torsion field. Varying
the above action with respect to metric and BD field leaves us with the
modified BD field equations

$\displaystyle-$
$\displaystyle\\!\\!\\!\\!\phi\,{G}_{\alpha\beta}-\frac{\phi}{2}\bigg{[}{K}^{\gamma\,\,\delta}_{\,\,\beta}\,{K}_{\delta\gamma\alpha}+{K}^{\gamma\,\,\delta}_{\,\,\alpha}\,{K}_{\delta\gamma\beta}-{g}_{\alpha\beta}{K}^{\gamma\delta\epsilon}\,{K}_{\delta\epsilon\gamma}+{g}_{\alpha\beta}{K}^{\gamma\delta}_{\,\,\,\,\,\gamma}\,{K}_{\delta\,\,\epsilon}^{\,\,\epsilon}-{K}^{\gamma}_{\,\,\alpha\beta}\,{K}^{\delta}_{\,\,\gamma\delta}$
(21) $\displaystyle-$
$\displaystyle\\!\\!\\!\\!{K}^{\gamma}_{\,\,\beta\alpha}\,{K}^{\delta}_{\,\,\gamma\delta}\bigg{]}+\frac{1}{2}{K}^{\gamma}_{\,\,\beta\gamma}\nabla_{\alpha}\phi+\frac{1}{2}{K}^{\gamma}_{\,\,\alpha\gamma}\nabla_{\beta}\phi+\nabla_{\alpha}\nabla_{\beta}\phi+\frac{\omega}{\phi}\nabla_{\alpha}\phi\nabla_{\beta}\phi$
$\displaystyle+$
$\displaystyle\\!\\!\\!\\!\frac{1}{2}\left[{K}^{\gamma}_{\,\,\alpha\beta}-{K}^{\gamma}_{\,\,\beta\alpha}\right]\nabla_{\gamma}\phi+\frac{1}{2}{g}_{\alpha\beta}\left[{K}^{\delta\gamma}_{\,\,\,\,\gamma}-{K}^{\gamma\delta}_{\,\,\,\,\gamma}\right]\nabla_{\delta}\phi+{g}_{\alpha\beta}\Box\phi$
$\displaystyle-$
$\displaystyle\\!\\!\\!\\!\frac{\omega}{2\phi}{g}_{\alpha\beta}\nabla_{\gamma}\phi\nabla^{\gamma}\phi=\frac{\kappa^{2}}{2}{T}_{\alpha\beta},$
$\displaystyle-$
$\displaystyle{K}^{\alpha\beta}_{\,\,\,\,\,\alpha}\,{K}_{\beta\,\,\,\gamma}^{\,\,\,\gamma}+K^{\alpha\beta\gamma}\,{K}_{\beta\gamma\alpha}+{R}-\nabla_{\alpha}{K}^{\alpha\beta}_{\,\,\,\,\,\beta}+\nabla_{\beta}{K}^{\alpha\beta}_{\,\,\,\,\,\alpha}+\frac{2\omega}{\phi}\Box\phi-\frac{\omega}{\phi^{2}}\nabla_{\alpha}\phi\nabla^{\alpha}\phi=0.$
(22)

where
${T}_{\alpha\beta}=2\left(\delta{L}_{m}/\delta{g}^{\alpha\beta}\right)/\sqrt{-{g}}$
is the EMT of matter fields. Varying the action with respect to contorsion
leads the modified Cartan field equation as

$\displaystyle{K}^{\mu}_{\,\,\,\alpha\beta}=-\frac{\kappa^{2}}{4\phi}\left[{S}_{\beta}^{\,\,\,\mu}{u}_{\alpha}+{S}_{\alpha}^{\,\,\mu}{u}_{\beta}+{S}_{\alpha\beta}{u}^{\mu}\right]+\frac{1}{2\phi}\left[\delta_{\beta}^{\,\,\mu}\nabla_{\alpha}\phi-{g}_{\alpha\beta}\nabla^{\mu}\phi\right].$
(23)

The quantity $S_{\mu\nu}=-S_{\nu\mu}$ is a second-rank antisymmetric tensor
which is defined as the spin density tensor. This quantity in EC theory bears
the contribution due to intrinsic angular momentum (spin) of fermionic
particles within a perfect fluid [38]. Therefore a classical description of
spin in EC theory generalizes the perfect fluid of GR to the so-called
Weyssenhoff fluid [39]. The EMT on the right hand side of (21) can be written
as sum of the EMT of a usual perfect fluid and an intrinsic spin part, as [40]

$\displaystyle{T}_{\alpha\beta}=\left[(\rho+p){u}_{\alpha}{u}_{\beta}-p{g}_{\alpha\beta}\right]+\frac{\kappa^{2}}{8\phi}\left[2{S}_{\gamma\delta}\,{S}^{\gamma\delta}{u}_{\alpha}{u}_{\beta}-2{S}_{\alpha}^{\,\,\,\delta}{S}_{\beta\delta}{u}_{\gamma}{u}^{\gamma}\right]+\frac{\nabla^{\gamma}\phi}{8\phi}\bigg{[}{S}_{\beta\gamma}{u}_{\alpha}+{S}_{\alpha\gamma}{u}_{\beta}\bigg{]}.$
(24)

Substituting then for the above EMT into field equations (21) and (22) and
after rearranging, we arrive at final version of the field equations in ECBD
theory, as

$\displaystyle{G}_{\alpha\beta}$ $\displaystyle=$ $\displaystyle
T_{\alpha\beta}^{{\rm
sf}}+T_{\alpha\beta}^{\phi}=\frac{\kappa^{2}}{\phi}\left[\left(\rho+p-\frac{\kappa^{2}\sigma^{2}}{6\phi}\right){u}_{\alpha}{u}_{\beta}-{g}_{\alpha\beta}\left(p-\frac{\kappa^{2}\sigma^{2}}{12\phi}\right)\right]+\frac{3+2\omega}{2\phi^{2}}\nabla_{\alpha}\phi\nabla_{\beta}\phi$
(25) $\displaystyle+$
$\displaystyle\frac{1}{\phi}\nabla_{\alpha}\nabla_{\beta}\phi-{g}_{\alpha\beta}\left[\frac{\Box\phi}{\phi}+\frac{3+2\omega}{4\phi^{2}}\nabla_{\epsilon}\phi\nabla^{\epsilon}\phi\right],$
$\displaystyle\Box\phi$ $\displaystyle=$
$\displaystyle\frac{\kappa^{2}}{2(3+\omega)}\left[\rho-3p-\frac{3\kappa^{2}\sigma^{2}}{8\phi}\right],$
(26)

where $T_{\alpha\beta}^{{\rm sf}}$ and $T_{\alpha\beta}^{\phi}$ are energy
momentum tensors of Weyssenhoff spin fluid and BD scalar field, respectively.
The quantity $\sigma^{2}$ is the square of spin density which is obtained by
suitable averaging procedure on spin density tensor [40]

$\displaystyle\langle{S}_{\mu\nu}{S}^{\mu\nu}\rangle=2{\sigma}^{2},$
$\displaystyle\langle{S}_{\mu}^{\,\,\beta}{S}_{\nu\beta}\rangle=\frac{2}{3}\left[{g}_{\mu\nu}-{u}_{\mu}{u}_{\nu}\right]{\sigma}^{2}.$
(27)

It is noteworthy that in the process of taking the average of spin density
tensor, we assumed a cloud of randomly oriented spin particles so that, the
average of the spin density tensor is assumed to be zero, i.e.,
$\langle{S}^{\mu\nu}\rangle=0$. From Eqs. (25) and (26) we observe that the
presence of spin effects appear as a negative pressure within the EMT and we
expect that such negative pressure term could prevent singularity formation
during the dynamical evolution of the collapse scenario. To investigate such
possibility we begin by substituting metric (4) into the above field
equations, this gives ($\kappa^{2}=1$)

$\displaystyle
3\left[H^{2}+\frac{k}{a^{2}}\right]=\frac{1}{\phi}\left[\rho-\frac{\sigma^{2}}{12\phi}\right]+\frac{2\omega+3}{4}\frac{\dot{\phi}^{2}}{\phi^{2}}-3H\frac{\dot{\phi}}{\phi},$
(28)
$\displaystyle-2\frac{\ddot{a}}{a}-\frac{\dot{a}^{2}}{a^{2}}-\frac{k}{a^{2}}=\frac{1}{\phi}\left[p-\frac{\sigma^{2}}{12\phi}\right]+2H\frac{\dot{\phi}}{\phi}+\frac{2\omega+3}{4}\frac{\dot{\phi}^{2}}{\phi^{2}}+\frac{\ddot{\phi}}{\phi},$
(29)
$\displaystyle\ddot{\phi}+3H\dot{\phi}=\frac{1}{2\omega+6}\left[\rho-3p-\frac{3\sigma^{2}}{8\phi}\right].$
(30)

We also have conservation equation (8) for the fluid part and the following
relation

$\displaystyle\frac{d\sigma^{2}}{dt}+6H\sigma^{2}=0,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\Longrightarrow\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \sigma^{2}(a)=\sigma_{i}^{2}a^{-6},$ (31)

for the spin part, where $\sigma_{i}$ is the initial value for the square of
spin density.

###  3.1 Approximate solutions

We proceed to find a semi-exact solution in the late stages of the collapse
scenario where the contracting phase is about to end and then be converted to
an expanding one. Similar to the case of BD theory, only equations (29) and
(30) are independent. Using then the change of variable $d/dt=aHd/da$, these
two equations can be rewritten as

$\displaystyle 2aHH_{,a}\\!\\!\\!\\!\\!$ $\displaystyle+$
$\displaystyle\\!\\!\\!\\!\\!3H^{2}+\frac{k}{a^{2}}+\frac{\rho_{i}a^{-3(1+w)}}{\phi}\left[w-\frac{1-3w}{2(\omega+3)}\right]-aH^{2}\frac{\phi_{,a}}{\phi}$
(32) $\displaystyle+$
$\displaystyle\frac{2\omega+3}{4}a^{2}H^{2}\left(\frac{\phi_{,a}}{\phi}\right)^{2}-\frac{\sigma_{i}^{2}}{4a^{6}\phi^{2}}\left[\frac{1}{3}+\frac{3}{4(\omega+3)}\right]=0,$
$\displaystyle\Big{(}aH^{2}\\!\\!\\!\\!\\!$ $\displaystyle+$
$\displaystyle\\!\\!\\!\\!\\!a^{2}HH_{,a}\Big{)}\phi_{,a}+a^{2}H^{2}\phi_{,aa}+3aH^{2}\phi_{,a}-\frac{\rho_{i}(1-3w)}{2(\omega+3)\phi}a^{-3(1+w)}+\frac{3\sigma_{i}^{2}}{16(\omega+3)a^{6}\phi^{2}}=0.$
(33)

The above system of differential equations cannot be solved analytically,
however, we may employ some approximations in order to simplify it. Consider a
collapse scenario that has started its evolution at $t=t_{i}$ and
$a_{i}=a(t_{i})$ and has entered a “near-bounce” regime where the collapse
process with the rate $H_{nb}$ is undergoing a decelerated contracting phase
and then halts at $a=a_{b}$, i.e., $H_{nb}(a_{b})=0$. In such a situation
$H_{nb}^{2}(a)\ll H_{nb}(a)$ and $a_{b}<a_{nb}\ll a_{i}$. Note that $a_{nb}$
is the scale factor at near bounce regime so that $a_{nb}(t_{b})=a_{b}$ where
$t_{b}$ is the time at which the bounce occurs. We therefore neglect the term
$k/a^{2}$ and those containing square of collapse rate. Equations (32) and
(33) then take the form

$\displaystyle
2a_{nb}H_{nb}H^{\prime}_{nb}+\frac{\rho_{i}}{\phi_{nb}}\left(w+\frac{3w-1}{2(\omega+3)}\right)a_{nb}^{-3(1+w)}-\frac{\sigma_{i}^{2}(21+4\omega)}{48a_{nb}^{6}(\omega+3)\phi_{nb}^{2}}=0,$
(34) $\displaystyle
a_{nb}^{2}H_{nb}H^{\prime}_{nb}\phi^{\prime}_{nb}-\frac{1}{2(\omega+3)}\left[\rho_{i}(1-3w)a_{nb}^{-3(1+w)}-\frac{3\sigma_{i}^{2}}{8a_{nb}^{6}\phi_{nb}}\right]=0,$
(35)

where a prime denotes differentiation with respect to $a_{nb}$. Now, if we set
the BD parameter as $\omega=(1/2w)\left(1-9w\right)$, we obtain an exact
solution for the above system as

$\displaystyle\phi_{nb}$ $\displaystyle=$
$\displaystyle\left\\{C_{1}[(2+3w)a_{nb}]^{\frac{18w}{2+3w}}-\frac{16w\rho_{i}a_{nb}^{3(1-w)}}{(w+2)\sigma_{i}^{2}}\right\\}^{-1},$
(36) $\displaystyle H_{nb}$ $\displaystyle=$
$\displaystyle\pm\sqrt{2}\left[\frac{2+3w}{288(3w-1)}\left(A_{1}a_{nb}^{-6w}+A_{2}a_{nb}^{-\frac{3(3w^{2}+w+2)}{2+3w}}+A_{3}a_{nb}^{6\frac{3w-2}{3w+2}}\right)+C_{2}\right]^{\frac{1}{2}},$
(37)

where $C_{1}$ and $C_{2}$ are constants of integration and

$A_{1}=\frac{256w\rho_{i}^{2}}{\sigma_{i}^{2}(w+2)^{2}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
A_{2}=-\frac{64C_{1}w(3w+2)^{\frac{21w+2}{3w+2}}\rho_{i}}{4+(3w+5)w^{2}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
A_{3}=-\frac{C_{1}^{2}\sigma_{i}^{2}(3w+2)^{\frac{39w+2}{3w+2}}}{3w-2}.$ (38)

Next, we proceed to obtain $C_{2}$ using the condition $H_{nb}(a_{b})=0$. This
gives

$\displaystyle C_{2}$ $\displaystyle=$
$\displaystyle\frac{3w+2}{288a_{b}^{6w}}\Bigg{[}\frac{64C_{1}\rho_{i}wa_{b}^{\frac{3w^{2}+9w-2}{3w+2}}(3w+2)^{\frac{21w+2}{3w+2}}}{(w+2)(3w-1)(w(3w-1)+2)}-\frac{256w\rho_{i}^{2}}{(w+2)^{2}(3w-1)\sigma_{i}^{2}}$
(39) $\displaystyle+$
$\displaystyle\frac{C_{1}^{2}\sigma_{i}^{2}(3w+2)^{\frac{39w+2}{3w+2}}a_{b}^{6\frac{3w^{2}+9w-2}{3w+2}}}{9w(w-1)+2}\Bigg{]}.$

By substituting the above constant into Eq. (37) we get the near bounce
solution for the collapse rate so that the negative sign stands for
contracting regime and the positive one for expanding regime. This behavior
has been sketched in Fig. (1), where we observe that $H_{nb}<0$ in the
contracting regime (family of red curves) and after the bounce occurs an
expanding phase (family of blue curves) with $H_{nb}>0$ commences. In the next
section we try to solve Eqs. (28)-(30) using numerical techniques and find the
non-singular collapse solutions under physically reasonable conditions.

Figure 1: Behavior of the collapse rate against the scale factor in near
bounce regime for $w=0.1$, $C_{1}=1$ and different initial values of energy
and square of spin densities. the bounce occurs at $a_{b}=0.2$ where
$H_{nb}(a_{b})=0$.

###  3.2 Numerical solutions

In the present section we try to solve numerically the system of equations
(28)-(30) and estimate the behavior of scale factor and other related
quantities over time. Substituting for $\ddot{\phi}$ and $\dot{a}^{2}$ from
Eqs. (30) and (28) into Eq. (29) gives

$\displaystyle
a^{3w}\left[48(\omega+3)a^{5}\left(6\phi^{2}\ddot{a}-6\phi\dot{a}\dot{\phi}+(2\omega+3)a\dot{\phi}^{2}\right)-\sigma_{i}^{2}(16\omega+75)\right]$
(40) $\displaystyle+$ $\displaystyle
24\rho_{i}a^{3}\phi(w(6\omega+9)+2\omega+9)=0.$

The above equation along with Eq. (30) have to be solved numerically with the
initial condition on collapse velocity, $\dot{a}(t_{i})=\dot{a}_{i}$ as
provided by Eq. (28). Figure (2) shows the time evolution of scale factor for
dust collapse assuming different initial values of energy and square of spin
densities. We observe that spin effects provide a completely different
evolution for the collapse setting in comparison to the case when these
effects are absent. As the family of black curves show, the collapsing body
continues its contracting phase until the spin effects become dominant and
consequently the collapse process halts at a bounce time where the scale
factor reaches a minimum value. After the bounce, the contracting regime is
replaced by an expanding one with an increasing scale factor. This non-
singular collapse scenario can be compared to the case of singular one, i.e.,
when the spin effects are not taken into account. In this case the collapsing
object continues its contraction until reaching a spacetime singularity, where
the scale factor vanishes, see the red curve.

Figure 2: Behavior of scale factor for $k=1$, $w=0$, $\omega=9$, $a_{i}=1$,
$\phi_{i}=1.5$, $\dot{\phi}_{i}=1.9$ and different values of initial energy
density in the left panel with $\sigma_{i}^{2}=12$ and square of spin density
in the right panel with $\rho_{i}=7.5$.

Comparing the two diagrams in Fig. (2) along with the left panel of Fig. (3),
one can deduce that the competition between initial energy and square of spin
densities can affect the time at which the bounce occurs and its corresponding
scale factor. In the left panel of Fig. (2), for a fixed value of
$\sigma_{i}^{2}$, increasing $\rho_{i}$ leads to smaller values of $a_{b}$ and
$t_{b}$. This scenario also occurs for a fixed value of $\rho_{i}$ but
decreasing $\sigma_{i}^{2}$, see the right panel. In other words, if we
interpret energy density and square of spin density as agents for attractive
nature of gravity and repulsive effects, then initially the more the
gravitational attraction or correspondingly the less the repulsive effects,
the more the collapsing object needs to be contracted until the spin effects
become dominant to finally prevent the singularity formation. The opposite
case for this scenario also holds for smaller values of $\rho_{i}$ and larger
ones for $\sigma_{i}^{2}$. Likewise, in the left panel of Fig. (3) we observe
that even if $\sigma_{i}^{2}<\rho_{i}$, in later stages of the collapse
process the repulsive effects due to spin contribution overcome the
gravitational pull of gravity to finally avoid the singularity. In the right
panel, we have sketched the evolution of scale factor against time for
different values of BD parameter. This parameter represents the strength of
scalar to tensor coupling to the matter so that the smaller the value of BD
parameter, the larger the contribution of the BD scalar field to gravitational
interaction. We can therefore deduce that for fixed initial values of energy
and spin densities, the BD parameter could affect the softness of the change
from a contracting regime to an expanding one. In this regard, we may conclude
that the stronger the contribution of the BD scalar field (which here also
acts as a source of dynamical torsion) to the gravitational interaction (i.e.,
the smaller values of $\omega$), the softer the bounce occurs. The left panel
in Fig. (4) shows the time evolution of BD scalar field where we observe that
in comparison to the case of vanishing spin effects (red curve), the BD field
is finite and assumes a bounded maximum value at the bounce. This diagram can
be compared to the left panel of Fig. (3) in such a way that, the lower the
initial value for square of spin density, the higher this maximum for BD field
would be and the more the collapsing cloud needs to be contracted until the
gravity succumbs to the repulsive spin effects. The right panel in Fig. (4)
shows the behavior of Kretschmann scalar over time where we observe that this
quantity behaves regularly throughout the collapse process and reaches a
maximum finite value at the bounce. It then converges to zero in the post-
bounce regime where the cloud disperses. However, in the absence of spin
effects, this quantity grows unboundedly and diverges at the spacetime
singularity, see the red curve. For the sake of physical validity of the
collapse model, it is required that the collapse process obeys the weak energy
condition (WEC), i.e., $T_{\alpha\beta}v^{\alpha}v^{\beta}\geq 0$ for all non-
spacelike vector fields. This condition means that the energy density as
measured by any local observer is non-negative. For EMT given in Eq. (25) the
WEC leads to the following inequalities

$\rho_{{\rm sf}}\geq 0,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \rho_{{\rm sf}}+p_{{\rm sf}}\geq
0,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \rho_{\phi}\geq
0,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\phi}+p_{\phi}\geq 0,$ (41)

or in terms of EMT components, using Eqs. (28) and (29)

$\frac{\rho}{\phi}-\frac{\sigma^{2}}{12\phi^{2}}\geq 0,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
(1+w)\frac{\rho}{\phi}-\frac{\sigma^{2}}{6\phi^{2}}\geq 0,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\frac{2\omega+3}{4}\frac{\dot{\phi}^{2}}{\phi^{2}}-3H\frac{\dot{\phi}}{\phi}\geq
0,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\frac{2\omega+3}{2}\frac{\dot{\phi}^{2}}{\phi^{2}}-H\frac{\dot{\phi}}{\phi}+\frac{\ddot{\phi}}{\phi}\geq
0.$ (42)

In Fig. (5) we have plotted the behavior of WEC over time where in the left
panel we observe that the energy density of spin fluid is finite (family of
black curves) and positive hence the first inequality of WEC for spin fluid in
Eq. (41) is satisfied. As is seen in the right panel, the second inequality
also holds throughout the evolution of collapsing body. It should be noted
that though the growth rate of square of spin density is three times bigger
(for the dust case) than that of matter density, it is the presence of
$\phi^{2}$ term in the denominator of square of spin density that balances the
growth rate of densities and thus preserving the WEC. Figure (6) indicates the
behavior of energy density profile for BD field and the corresponding WEC,
i.e., the last part of Eq. (41). It is seen that $\rho_{\phi}$ reaches a
maximum positive value at the bounce point and then converges to negative
values in the post bounce regime where the collapse rate changes its sign from
contracting phase to an expanding one. It then tends to zero at later times.
This behavior also holds for $\rho_{\phi}+p_{\phi}\geq 0$. Thus, the WEC for
scalar field profile is violated in the post bounce regime, however for the
effective EMT, as given by right hand side of Eqs. (28) and (29), the WEC is
fulfilled.

We also note that when spin effects are absent i.e., $\sigma_{i}^{2}=0$, the
ECBD field equations (28)-(30) will be reduced to their BD counterparts111The
geometry is still non-Riemannian as from Eq. (23) the contorsion tensor is
nonzero., i.e., Eqs. (5)-(7), if we take the BD parameter as
$\omega\rightarrow\omega-3/2$. However, the collapse setting arising from this
case is singular as we observe from the red curves in Figs. (2)-(4), hence, it
is the mere presence of repulsive spin effects that could avoid singularity
formation in BD theory. Finally, in the limit where $\omega\rightarrow\infty$
and $\phi\rightarrow{\rm const.}$, the ECBD theory reduces to EC theory in
which a nonsingular collapse model has been studied earlier [11].

###  3.3 Dynamics of apparent horizon

An important point that needs to be investigated in every collapse setting is
the behavior of apparent horizon which is the outermost boundary of the
trapped region. Indeed, the apparent horizon determines the boundary between
null trajectories that have light-cones confined in the trapped region and
thus are causally disconnected from the rest of the universe, and those
trajectories that can propagate to far away observers. The condition for the
formation of apparent horizon is given by the requirement that the surface
${\mathcal{R}}={\rm const.}$ is null, i.e.,
$g^{\mu\nu}\partial_{\mu}{\mathcal{R}}\partial_{\nu}{\mathcal{R}}=0$. Hence, a
shell labeled by the comoving radial coordinate $r$ will become trapped if
$\dot{{\mathcal{R}}}_{ah}^{2}(t,r)=1$ [14]. For the line element (4) this
condition gives the physical areal radius of apparent horizon as

${\mathcal{R}}_{ah}(t,r)=\frac{a}{\sqrt{1+\dot{a}^{2}}}.$ (43)

The above equation describes time dependence of apparent horizon curve which
indeed implies that depending on the behavior of the collapse velocity, the
apparent horizon could either form or be avoided. This depends on how the
initial radius of the collapsing cloud, i.e.,
${\mathcal{R}}(t_{i},r_{\Sigma})=r_{\Sigma}$, is chosen. Fig. (7) shows the
dynamics of apparent horizon where in the left panel, we observe that in the
singular case, the apparent horizon is a decreasing function of time and as
$t\rightarrow t_{s}$ the areal radius of apparent horizon tends to zero to
finally cover the spacetime singularity. In the case of non-vanishing spin
effects, the apparent horizon decreases for a while in the contracting regime
reaching its first local minimum (Point A) and then increases to a maximum
value at the bounce point. It then converges to the next nonzero local minimum
(Point B) in the expanding phase and monotonically grows during the post-
bounce regime. We therefore find that the dynamics of apparent horizon is
completely different when spin effects are present and depending on initial
radius of the collapsing body it can be totally avoided. More precisely, if we
choose the boundary of the cloud at onset of collapse as $r_{\Sigma}=r_{1}$,
then as the collapse proceeds, the collapse velocity and the scale factor will
achieve the values $\left\\{\dot{a}_{1},a_{1}\right\\}$ at point A and
consequently the apparent horizon equation (43) will be satisfied. This means
that a dynamical horizon would form in the contracting regime to meet the
boundary of the collapsing object. For $r_{2}<r_{\Sigma}<r_{1}$, the horizon
equation will be fulfilled two times in the expanding regime, however, the
contracting regime will be free of horizons. For $r_{\Sigma}=r_{2}$ the
horizon equation will be satisfied once, i.e., when the collapse velocity and
the scale factor reach the values $\left\\{\dot{a}_{2},a_{2}\right\\}$ at
point B and hence only one horizon will form to intersect the boundary of the
cloud in the expanding phase. Finally for $r_{\Sigma}<r_{2}$, we have no
horizon formation throughout the collapsing and expanding regimes and thus,
the bounce can be visible to an external observer. We can therefore conclude
that, unlike the singular case, it is the boundedness of the collapse velocity
that, depending on suitable values of $r_{\Sigma}$, allows the avoidance of
horizon formation during the dynamical evolution of the collapsing object. In
the right panel we have sketched the time behavior of apparent horizon curve
for different values of square of spin density. It is therefore seen that the
lower the initial value of square of spin density, the smaller the initial
radius of the collapsing body must be taken in order to avoid horizon
formation.

Figure 3: Behavior of scale factor for $k=1$, $w=0$, $\omega=9$, $a_{i}=1$,
$\phi_{i}=1.5$, $\dot{\phi}_{i}=1.9$, $\rho_{i}=14$ and different initial
values of square of spin density in the left panel and different values of BD
parameter with $\rho_{i}=10$ and $\sigma_{i}^{2}=8$ in the right panel.

Figure 4: (Left panel) Behavior of BD scalar field for $k=1$, $w=0$,
$\omega=9$, $a_{i}=1$, $\phi_{i}=1.5$, $\dot{\phi}_{i}=1.9$, $\rho_{i}=14$ and
different initial values of square of spin density. (Right panel) Time
evolution of Kretschmann invariant for the same values of model parameters as
of the left panel.

Figure 5: (Left panel): Behavior of energy density of spin fluid for $k=1$,
$w=0$, $\omega=9$, $a_{i}=1$, $\phi_{i}=1.5$, $\dot{\phi}_{i}=1.9$,
$\rho_{i}=7.5$ and different initial values of square of spin density. (Right
panel) Weak energy condition for spin fluid for the same values of model
parameters as of the left panel. The red curve shows the behavior of these
quantities when spin effects are not taken into account.

Figure 6: (Left panel): Behavior of energy density of BD scalar field for
$k=1$, $w=0$, $\omega=9$, $a_{i}=1$, $\phi_{i}=1.5$, $\dot{\phi}_{i}=1.9$,
$\rho_{i}=7.5$ and different initial values of square of spin density. (Right
panel) Weak energy condition for BD scalar field for the same values of model
parameters as of the left panel. The red curve shows the behavior of these
quantities when spin effects are not taken into account.

Figure 7: Evolution of apparent horizon curve for $k=1$, $w=0$, $\omega=9$,
$a_{i}=1$, $\phi_{i}=1.5$, $\dot{\phi}_{i}=1.9$, $\rho_{i}=7.5$,
$\sigma_{i}^{2}=12$ (left panel) and different values of square of spin
density (right panel).

##  4 Exterior spacetime

The collapse scenario we have studied so far deals with the interior spacetime
of the contracting object. In order to complete the collapse model we need to
match the interior spacetime to a suitable exterior region whose boundary
$r=r_{\Sigma}$ is the surface of the collapsing matter. Let us take $({\rm
V}^{\pm},g^{\pm})$ as four dimensional spacetimes and $\Sigma$ as a timelike
three-dimensional smooth hypersurface separating spacetime into two regions
${\rm V}^{+}$ (exterior) and ${\rm V}^{-}$ (interior) with scalar fields
$\phi^{+}$ and $\phi^{-}$, respectively. The hypersurface $\Sigma$ then
results from isometric pasting of two hypersurfaces $\Sigma^{+}$ and
$\Sigma^{-}$, which respectively, bound the four-dimensional exterior (${\rm
V}^{+}$) and interior (${\rm V}^{-}$) spacetimes. For the interior region we
take the line element (4) in the FLRW form as

$ds_{-}^{2}=-dt^{2}+\frac{a^{2}(t)}{1-kr^{2}}dr^{2}+a^{2}(t)(r^{2}d\theta^{2}+r^{2}\sin^{2}\theta
d\varphi^{2}),$ (44)

where the interior coordinates are labeled as
$\left\\{x_{-}^{\mu}\right\\}\equiv\left\\{t,r,\theta,\varphi\right\\}$.
Utilizing generalized Israel-Darmois junction conditions [41] we proceed to
match the above line element to the exterior one through the boundary surface
$r=r_{\Sigma}$. The line element for the exterior region is chosen as
generalized Vaidya metric [42] which in retarded (exploding) null coordinates
reads

$ds_{+}^{2}=-f(v,{\mathtt{R}})dv^{2}-2dvd{\mathtt{R}}+{\mathtt{R}}^{2}(d\theta^{2}+\sin^{2}\theta
d\varphi^{2}),$ (45)

where $f(v,{\mathtt{R}})=1-2{\mathtt{M}(v,\mathtt{R})}/{\mathtt{R}}$ with
${\mathtt{M}(v,\mathtt{R})}$ being the Vaidya mass and the exterior
coordinates are labeled as
$\left\\{x_{\rm+}^{\mu}\right\\}\equiv\left\\{v,{\mathtt{R}},\theta,\varphi\right\\}$.
The intrinsic line element on hypersurface $\Sigma$ is given by [43]

$ds_{\Sigma}^{2}=-d\tau^{2}+Y^{2}(\tau)(d\theta^{2}+\sin^{2}\theta
d\varphi^{2}),$ (46)

where $y^{a}=\left\\{\tau,\theta,\varphi\right\\}$, $(a=0,2,3)$ are the
intrinsic coordinates of $\Sigma$ with $\tau$ being the time coordinate
defined on it and we have assumed the angular coordinates $\theta$ and
$\varphi$ are continuous across $\Sigma$. The governing equations of
hypersurface $\Sigma$ in the coordinates $x_{\pm}^{\mu}$ are given by

$r-r_{\Sigma}=0\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \text{in}\leavevmode\nobreak\ {\rm
V}^{+},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
{\mathtt{R}}-{\mathtt{R}}_{\Sigma}=0\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\text{in}\leavevmode\nobreak\ {\rm V}^{-},$ (47)

with the help of which we get the interior and exterior induced metrics on
$\Sigma^{-}$ and $\Sigma^{+}$, respectively as

$ds_{\Sigma^{-}}^{2}=-dt^{2}+a^{2}(t)r_{\Sigma}^{2}(d\theta^{2}+\sin^{2}\theta
d\varphi^{2}),$ (48)

and

$\displaystyle
ds_{\Sigma^{+}}^{2}=-\left[f\big{(}v(t),{\mathtt{R}}(t)\big{)}\dot{v}^{2}+2\dot{{\mathtt{R}}}\dot{v}\right]dt^{2}+{\mathtt{R}}^{2}(t)(d\theta^{2}+\sin^{2}\theta
d\varphi^{2}).$ (49)

The junction conditions for the first fundamental forms i.e., induced metrics,
then give

$ds_{\Sigma}^{2}=ds_{\Sigma^{-}}^{2}=ds_{\Sigma^{+}}^{2},$ (50)

whence we get

$f\big{(}v(t),{\mathtt{R}}(t)\big{)}\dot{v}^{2}+2\dot{{\mathtt{R}}}\dot{v}=1,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
{\mathtt{R}}(t)=r_{\Sigma}a(t)=Y(\tau),$ (51)

where an overdot denotes $d/dt$. Next, we proceed to calculate the extrinsic
curvature tensors of $\Sigma^{+}$ and $\Sigma^{-}$. To this aim we need unit
vector fields normal to these hypersurfaces which are given by

$\displaystyle
n^{-}_{\mu}=\left[0,\frac{a(t)}{\sqrt{1-r^{2}}},0,0\right],\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
n^{+}_{\mu}=\frac{1}{\left[f(v,{\mathtt{R}})\dot{v}^{2}+2\dot{{\mathtt{R}}}\dot{v}\right]^{\frac{1}{2}}}\left[-\dot{\mathtt{R}},\dot{v},0,0\right].$
(52)

Let us take $x^{\mu}=x^{\mu}(y^{a})$ as the parametric equations of the
hypersurface $\Sigma$. The extrinsic curvature or second fundamental form of
this hypersurface is a three-tensor defined as [43]

${\mathcal{K}}_{ab}\equiv e_{a}^{\nu}e_{b}^{\mu}\nabla_{\mu}n_{\nu},$ (53)

where $e_{a}^{\nu}=\partial x^{\nu}/\partial y^{a}$ are the basis vectors
tangent to the hypersurafce $\Sigma$ and the covariant derivative is taken
with respect to the Christoffel $\Gamma^{\gamma}_{\,\nu\mu}$ connection. The
extrinsic curvature tensors associated with $\Sigma^{\pm}$ then read

${\mathcal{K}}^{\pm}_{ab}=-n_{\mu}^{\pm}\left[\frac{\partial^{2}x_{\pm}^{\mu}}{\partial
y^{a}\partial y^{b}}+\tilde{\Gamma}^{\mu\pm}_{\nu\sigma}\frac{\partial
x_{\pm}^{\nu}}{\partial y^{a}}\frac{\partial x_{\pm}^{\sigma}}{\partial
y^{b}}\right].$ (54)

where $x_{+}^{\mu}(y^{a})$ and $x_{-}^{\mu}(y^{a})$ are parametric relations
for the hypersurfaces $\Sigma^{+}$ and $\Sigma^{-}$ and
$\tilde{\Gamma}^{\mu\pm}_{\nu\sigma}$ are the components of connections
associated to interior and exterior line elements. It should be noted that, to
calculate the extrinsic curvature components for the interior spacetime we
should take the connection as given by Eq. (17) while for the exterior
spacetime one may use the usual Christoffel connection.

From Eqs. (25) and (26) we get the following expression for Ricci tensor

$\displaystyle
R_{\alpha\beta}=\frac{\kappa^{2}}{\phi}\left[T_{\alpha\beta}^{{\rm
sf}}-\frac{2\omega+5}{4(\omega+3)}g_{\alpha\beta}T^{{\rm
sf}}\right]+\frac{3+2\omega}{2\phi^{2}}\nabla_{\alpha}\phi\nabla_{\beta}\phi+\frac{1}{\phi}\nabla_{\alpha}\nabla_{\beta}\phi,$
(55)

with the help of which the junction conditions take the form [41]

$\displaystyle\left[{\mathcal{K}}^{a}\,\\!\\!_{b}\right]$ $\displaystyle=$
$\displaystyle\frac{1}{\phi}\left(\Pi^{a}\,\\!\\!_{b}-\frac{2\omega+5}{4(\omega+3)}\Pi\delta^{a}\,\\!\\!_{b}\right),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\left[{\mathcal{K}}\right]=-\frac{(2\omega+3)\Pi}{4(\omega+3)\phi},$ (56)
$\displaystyle\left[\phi_{,n}\right]$ $\displaystyle=$
$\displaystyle\frac{\Pi}{2(\omega+3)},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\left[\phi\right]=0,$ (57)

where the notation $\left[\Psi\right]=\Psi^{+}|_{\Sigma}-\Psi^{-}|_{\Sigma}$
stands for the jump of given field across the hypersurface $\Sigma$, $n$
labels the coordinate normal to this surface and $\Pi_{ab}={\rm
diag}(\rho_{\rm s},p_{\rm s},p_{\rm s})$ is the stress-energy tensor of matter
fields (except the BD scalar field) on the shell located at $\Sigma$. The
quantities ${\mathcal{K}}$ and $\Pi=2p_{\rm s}-\rho_{\rm s}$ are traces of
${\mathcal{K}}^{a}\,\\!\\!_{b}$ and $\Pi^{a}\,\\!\\!_{b}$ with $\rho_{\rm s}$
and $p_{\rm s}$ being the surface energy density and surface pressure,
respectively. We also note that from equation (56) we can find the surface
stress-energy as

$\Pi^{a}\,\\!\\!_{b}={\phi}\left(\left[{\mathcal{K}}^{a}\,\\!\\!_{b}\right]-\frac{2\omega+5}{2\omega+3}\left[{\mathcal{K}}\right]\delta^{a}\,\\!\\!_{b}\right).$
(58)

The non-vanishing components of Christoffel connection for the exterior
spacetime are found as

$\displaystyle\Gamma^{v+}_{vv}=-\frac{f_{,{\mathtt{R}}}}{2},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\Gamma^{{\mathtt{R}}+}_{vv}=\frac{1}{2}\left(f_{,v}+ff_{,{\mathtt{R}}}\right),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\Gamma^{{\mathtt{R}}+}_{v{\mathtt{R}}}=\frac{f_{,{\mathtt{R}}}}{2},$ (59)

where $,{\mathtt{R}}\equiv\partial/\partial{\mathtt{R}}$ and
$,v\equiv\partial/\partial v$. Substituting the above expressions into
equation (54) along with using the second part of (52) we get

$\displaystyle{\mathcal{K}}^{+t}_{t}=\frac{\dot{v}^{2}\left[ff_{,{\mathtt{R}}}\dot{v}+f_{,v}\dot{v}+3f_{,{\mathtt{R}}}\dot{\mathtt{R}}\right]+2\left(\dot{v}\ddot{{\mathtt{R}}}-\dot{{\mathtt{R}}}\ddot{v}\right)}{2\left(f\dot{v}^{2}+2\dot{\mathtt{R}}\dot{v}\right)^{\frac{5}{2}}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
{\mathcal{K}}^{+\theta}_{\theta}$ $\displaystyle=$
$\displaystyle{\mathcal{K}}^{+\phi}_{\phi}=\frac{f\dot{v}+\dot{\mathtt{R}}}{{\mathtt{R}}\sqrt{f\dot{v}^{2}+2\dot{\mathtt{R}}\dot{v}}},$
(60)

where we note that
$\partial\theta/\partial\tau=\partial\varphi/\partial\tau=0$. To compute the
components of extrinsic curvature for $\Sigma^{-}$ we begin with the
connection for the interior spacetime. Using Eqs. (17) and (23) this
connection reads

$\tilde{\Gamma}^{\mu}_{\,\alpha\beta}=\Gamma^{\mu}_{\,\alpha\beta}-\frac{\kappa^{2}}{4\phi}\left[{S}_{\beta}^{\,\,\,\mu}{u}_{\alpha}+{S}_{\alpha}^{\,\,\mu}{u}_{\beta}+{S}_{\alpha\beta}{u}^{\mu}\right]+\frac{1}{2\phi}\left[\delta_{\beta}^{\,\,\mu}\nabla_{\alpha}\phi-{g}_{\alpha\beta}\nabla^{\mu}\phi\right].$
(61)

As we mentioned earlier, for a randomly oriented spin source
$\langle{S}^{\mu\nu}\rangle=0$, hence, Eq. (54) can be rewritten as follows

$\displaystyle{\mathcal{K}}^{-}_{ab}=-n_{\mu}^{-}\Bigg{[}\frac{\partial^{2}x_{-}^{\mu}}{\partial
y^{a}\partial y^{b}}+\Gamma^{\mu-}_{\nu\sigma}\frac{\partial
x_{-}^{\nu}}{\partial y^{a}}\frac{\partial x_{-}^{\sigma}}{\partial
y^{b}}+\frac{1}{2\phi}\left(\delta_{\sigma}^{\,\,\mu}\nabla_{\nu}\phi-{g}_{\nu\sigma}\nabla^{\mu}\phi\right)\frac{\partial
x_{-}^{\nu}}{\partial y^{a}}\frac{\partial x_{-}^{\sigma}}{\partial
y^{b}}\Bigg{]}.$ (62)

A straightforward calculation reveals that

$\displaystyle{\mathcal{K}}^{-t}\,\\!\\!_{t}$ $\displaystyle=$ $\displaystyle
0,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
{\mathcal{K}}^{-\theta}\,\\!\\!_{\theta}={\mathcal{K}}^{-\varphi}\,\\!\\!_{\varphi}=\frac{\sqrt{1-r_{\Sigma}^{2}}}{r_{\Sigma}a(t)}.$
(63)

It is not difficult to show that the above relations are the only non-
vanishing components of extrinsic curvature of $\Sigma^{-}$. For example,
bearing in mind that $\phi=\phi(t)$ along with considering nonzero components
of $\Gamma^{\mu-}_{\nu\sigma}$ we get

${\mathcal{K}}^{-}_{t\theta}={\mathcal{K}}^{-}_{t\varphi}={\mathcal{K}}^{-}_{\theta\varphi}=0.$
(64)

The jump for the normal derivative of the BD scalar field across $\Sigma$ is
found as

$\displaystyle\left[\phi_{,n}\right]$ $\displaystyle=$ $\displaystyle
n^{\\!\\!+a}\phi^{+}_{,a}-n^{\\!\\!-a}\phi^{-}_{,a}=n^{\\!\\!+v}\phi^{+}_{,v}+n^{\\!\\!+{\mathtt{R}}}\phi^{+}_{,{\mathtt{R}}}-n^{\\!\\!-r}\phi^{-}_{,r}$
(65) $\displaystyle=$
$\displaystyle\frac{1}{\left(f\dot{v}^{2}+2\dot{{\mathtt{R}}}\dot{v}\right)^{\frac{1}{2}}}\left\\{(\dot{{\mathtt{R}}}+f\dot{v})\phi^{+}_{,{\mathtt{R}}}-\dot{v}\phi^{+}_{,v}\right\\},$

where use has been made of the contravariant components of the normal vector
fields

$n^{\\!\\!-\mu}=[0,\frac{\sqrt{1-r^{2}}}{a(t)},0,0],\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
n^{\\!\\!+\mu}=\frac{1}{\left(f\dot{v}^{2}+2\dot{{\mathtt{R}}}\dot{v}\right)^{\frac{1}{2}}}[-\dot{v},\dot{{\mathtt{R}}}+f\dot{v},0,0].$
(66)

Thus, from Eqs. (56) and (57), the jump in the normal derivative of BD scalar
field and its continuity across $\Sigma$ gives

$\frac{1}{\left(f\dot{v}^{2}+2\dot{{\mathtt{R}}}\dot{v}\right)^{\frac{1}{2}}}\left\\{(\dot{{\mathtt{R}}}+f\dot{v})\phi^{+}_{,{\mathtt{R}}}-\dot{v}\phi^{+}_{,v}\right\\}=-\frac{\phi_{\Sigma}}{2\omega+3}\left\\{{\cal
K}^{+t}\,\\!\\!_{t}+2\left({\cal
K}^{+\theta}\,\\!\\!_{\theta}-K^{-\theta}\,\\!\\!_{\theta}\right)\right\\},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\phi^{+}|_{\Sigma}=\phi^{-}|_{\Sigma}=\phi_{\Sigma}.$ (67)

Once the matter distribution on boundary surface $\Sigma$ is specified, the
behavior of BD scalar field in the exterior region along with the exterior
metric function will be determined through Eqs. (51), (56) and Eq. (67). If we
assume the boundary surface to be free of energy density and pressure, i.e.,
$\rho_{\rm s}=0$ and $p_{\rm s}=0$, we find from Eq. (56) ${\cal
K}^{+\theta}\,\\!\\!_{\theta}={\cal K}^{-\theta}\,\\!\\!_{\theta}$ and ${\cal
K}^{+t}\,\\!\\!_{t}=0$. Therefore, the continuity of $(\theta,\theta)$ and
$(t,t)$ components of extrinsic curvature across $\Sigma$ leads to the
following relations

$\displaystyle f\dot{v}+\dot{{\tt R}}$ $\displaystyle=$
$\displaystyle\sqrt{1-r_{\Sigma}^{2}},$ (68)
$\displaystyle\dot{v}^{2}\left[(ff_{,{\tt R}}+f_{,v})\dot{v}+3f_{,{\tt
R}}\dot{{\tt R}}\right]$ $\displaystyle+$ $\displaystyle
2\left(\dot{v}\ddot{{\tt R}}-\dot{{\tt R}}\ddot{v}\right)=0,$ (69)

where we have used the relations given in Eq. (51). Taking derivatives of Eq.
(68) and the first part of Eq. (51) we obtain

$2\dot{{\tt R}}\ddot{v}-\dot{f}\dot{v}^{2}=0,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ 2\dot{{\tt R}}\ddot{{\tt
R}}+\dot{f}\dot{v}\left(2\dot{{\tt R}}+f\dot{v}\right)=0,$ (70)

whereby we get

$\ddot{{\tt R}}\dot{v}=-\ddot{v}\left(2\dot{{\tt R}}+f\dot{v}\right).$ (71)

Substituting for $\ddot{{\tt R}}\dot{v}$ and $\dot{{\tt R}}\ddot{v}$ from the
above relations into Eq. (69) and simplifying the result we finally get

${\cal K}^{+t}\,\\!\\!_{t}=-\frac{f_{,v}\dot{v}^{2}}{2\dot{{\tt R}}}=0,$ (72)

which implies that $f(v,{\tt R})$ must be a function of ${\tt R}$, only. Now,
solving Eq. (68) and the first part of Eq. (51) we find the four-velocity of
the boundary as

${\rm V}^{\alpha}=\left(\dot{v},\dot{\tt
R},0,0\right)=\left[\frac{\sqrt{1-r_{\Sigma}^{2}}\mp\sqrt{1-f-r_{\Sigma}^{2}}}{f},\pm\sqrt{1-f-r_{\Sigma}^{2}},0,0\right],$
(73)

where negative (positive) sign denotes contracting (expanding) regime. The
second part of Eq. (67) implies that the BD scalar field must be homogeneous
in the exterior region thus, $\phi^{+}_{,{\tt R}}=0$. Using this result in the
first part of Eq. (67) we find $\phi^{+}_{,v}=0$ implying that the BD scalar
field must be constant in the exterior region. It is therefore seen that the
exterior region is a static spacetime with dynamical boundary. From the second
part of Eq. (51) and radial component of the four-vector velocity we get
$r_{\Sigma}^{2}\dot{a}^{2}=1-f-r_{\Sigma}^{2}$. From the exterior viewpoint
the apparent horizon is located at $f=0$ and it would meet the boundary
surface if the collapse velocity satisfies

$r_{\Sigma}=\frac{1}{\sqrt{1+\dot{a}^{2}}}.$ (74)

The above expression is nothing but the condition of apparent horizon
formation in the interior spacetime, i.e., Eq. (43) taking ${\cal R}_{\rm
ah}=r_{\Sigma}a$. We therefore conclude that in the singular case where the
collapse velocity diverges, the initial radius of the boundary surface must be
set to zero in order to avoid horizon formation. This situation is obviously
nonphysical as the collapsing body has nonzero initial radius and consequently
there is no minimum value for $r_{\Sigma}$ so that the apparent horizon can be
avoided. On the other hand, if the velocity of collapse scenario is bounded,
as it happens in the nonsingular case, the boundary of the collapsing object
can be chosen so that the formation of apparent horizon is avoided, see
section (3.3). Thus, the bounce event can be visible to the observers in
exterior spacetime.

##  5 Concluding Remarks

The longstanding issue of gravitational collapse of compact objects in the
Universe has been an important subject of astronomical study. Over the past
decades, much efforts have been devoted to understanding the dynamics and
final fate of a collapse process. The models that have been proposed so far
describe both singular and nonsingular collapse scenarios, where in the former
the collapse process necessarily ends in formation of a spacetime singularity,
either naked or covered by black hole horizon, and, in the latter the
singularity event is prevented. Several research works in line of singularity
avoidance have been carried out in both classical and quantum collapse
scenarios. The results show that modifications to GR in classical regime or
quantum corrections can alter the whole collapse dynamics and finally replace
the spacetime singularity that occurs in GR by a nonsingular bounce [9, 10,
11, 14, 15, 16, 17, 18, 19, 44]. In the present work we tried to deal with
this issue in a modified version of BD theory where the spin effects of
collapsing matter are introduced as additional degree of freedom. In the
original version of this theory it is shown that the collapse process leads to
singularity formation covered by an event horizon [1]. However, when spin
effects are present, the dynamical behavior of the collapse setting would be
completely different. The study of modified field equations showed that the
spin effects come into play as a repulsive pressure that oppose against the
attractive pull of gravity and finally replace the spacetime singularity by a
nonsingular event. Moreover, the spin effects can change horizon dynamics so
that in the singular case where these effects are absent the radius of
apparent horizon decreases to finally cover the singular region. While, in the
nonsingular case there is a minimum value for the apparent horizon curve so
that a suitable choose of the initial radius of the collapsing body could
prevent horizon formation. The exterior solution was studied through a smooth
matching of interior to exterior spacetimes using Israel-Darmois junction
conditions. It was found the exterior spacetime is described by a
Schwarzschild metric in retarded null coordinates which admits a dynamical
boundary. We therefore concluded that the apparent horizon does not meet the
boundary through which the interior and exterior regions are joined, hence,
the bounce event can be detected by external observers.

Data Availability Statement Data sets generated during the current study are
available from the corresponding author on reasonable request.

## References

  * [1] M. A. Scheel, S. L. Shapiro and S. A. Teukolsky, Phys. Rev. D, 51 4208 (1995); Phys. Rev. D 51, 4236 (1995).
  * [2] J. R. Oppenheimer and H. Snyder, Phys. Rev. 56, 455 (1939); B. Datt, Zeitschrift für Physik 108, 314 (1938).
  * [3] S. W. Hawking and R. Penrose, Proc. R. Soc. A 314, 529 (1970);
S. W. Hawking, G.F.R. Ellis, “The Large Scale Structure of Space-Time”
Cambridge University Press, Cambridge, (1973);
J. M. M. Senovilla, D. Garfinkle, Class. Quantum Gravity 32, 124008 (2015).

  * [4] F. J. Tipler, C. J. S. Clarke and G. F. R. Ellis, “In General Relativity and Gravitation: One Hundred Years After the Birth of Albert Einstein,” ed. by A. Held (Plenum Press, New York, 1980);
J. M. M. Senovilla, Gen. Relativ. Gravit. 29, 701 (1997).

  * [5] D. Malafarina, Universe 3, 48 (2017).
  * [6] P. S. Joshi and D. Malafarina, Int. J. Mod. Phys. D 20, 2641 (2011).
  * [7] K. Mosani, D. Dey, P. S. Joshi, G. C. Samanta, H. Menon, V. D. Patel, rXiv:2106.01773 [gr-qc]; R. Shaikh and P. S. Joshi, Phys. Rev. D , 024033 (2018); G. Abbas, M. Tahir, Eur. Phys. J. C, 77 537 (2017); A. Chatterjee, A. Ghosh and S. Jaryal, Phys. Rev. D 106, 044049 (2022); G. Manna, EPJC 80, 813 (2020); A. H. Ziaie, H. Moradpour, S. Ghaffari, Phys. Lett. B 793 276 (2019); N. Dadhich, S. G. Ghosh, S. Jhingan, Phys. Rev. D 88, 084024 (2013); M. D. Mkenyeleye, R. Goswami, and S. D. Maharaj, Phys. Rev. D 90, 064034 (2014).
  * [8] T. Harada, Pramana 63, 741 (2004); J. B. Manchak, Phil. Sci. 81 1066 (2014); E. Curiel, Phil. Sci. Supp. Proc. 66 S119 (1999); P. S. Joshi, “Spacetime singularities, in Springer Handbook of Spacetime,” ed. by A. Ashtekar, V. Petkov. Springer Handbooks. (Springer, Berlin, Heidelberg, 2014).
  * [9] K. Bamba, S. Nojiri, S. D. Odintsov, Phys. Lett. B 698, 451 (2011).
  * [10] A. H. Ziaie, P. V. Moniz, A. Ranjbar, H. R. Sepangi, Eur. Phys. J. C 74, 3154 (2014).
  * [11] M. Hashemi, S. Jalalzadeha and A. H. Ziaie, Eur. Phys. J. C 75, 53 (2015).
  * [12] B. Kuchowicz, Gen. Relativ. Grav., 9, 511 (1978); M. Gasperini, Phys. Rev. Lett. 56, 2873 (1986); M. Gasperini Gen. Rel. Grav. 30 1703 (1998); K. Atazadeh, JCAP 06, 020 (2014); Q. Huang, P. Wu and H. Yu, Phys. Rev. D 91, 103502 (2015).
  * [13] C. Bambi, D. Malafarina, A. Marciano and L. Modesto, Phys. Lett. B 734, 27 (2014).
  * [14] P. S. Joshi, “Gravitational Collapse and Spacetime Singularities,” Cambridge University Press (2012).
  * [15] V. Husain, Mathematical Physics, 354-361 (2007); R. Torres and F. Fayos, Phys. Lett. B 747, 245 (2015); L. Modesto, J. W. Moffat and P. Nicolini, Phys. Lett. B 695, 397 (2011); L. Modesto, Phys. Rev. D 86, 044005 (2012); L. Modesto, Astron. Rev. 8, 4 (2013).
  * [16] R. Goswami, P. S. Joshi and P. Singh, Phys. Rev. Lett. 96 031302 (2006).
  * [17] V. Husain, Adv. Sci. Lett. 2, 214 (2009).
  * [18] C. Bambi, D. Malafarina and L. Modesto, Phys. Rev. D 88 044009 (2013); B. K. Tippett and V. Husain, Phys. Rev. D 84, 104031 (2011); C. Bambi, D. Malafarina, L. Modesto, Eur. Phys. J. C 74, 2767 (2014).
  * [19] R. Casadio, S. D. H. Hsu and B. Mirza, Phys. Lett. B 695 317 (2011); J. Marto, Y. Tavakoli and P. V. Moniz, Int. J. Mod. Phys. D 24, 1550025 (2015); M. Bojowald, R. Goswami, R. Maartens and P. Singh, Phys. Rev. Lett. 95, 091302 (2005); Y. Tavakoli, J. Marto and A. Dapor, Int. J. Mod. Phys. D 23 1450061 (2014); C. Bambi, D. Malafarina and L. Modesto, EPJC 74, 2767 (2014); D. Malafarina, Universe 3, 48 (2017).
  * [20] P. Singh and A. Toporensky, Phys. Rev. D 69 104008 (2004); M. Bojowald, Living Rev. Relativ. 11, 4 (2008); J.-Q. Guo and P. S. Joshi, Phys. Rev. D 94, 044063 (2016); M. Fathi, S., Jalalzadeh, P. V. Moniz, Eur. Phys. J. C 76, 527 (2016); J. M. Hernandez, M. Bellini, C. Moreno, Phys. Dark Univ. 23 100251 (2019); C. Kiefer, N. Kwidzinski, D. Piontek, Eur. Phys. J. C 79, 686 (2019); S. Jalalzadeh, Phys. Lett. B 833, 137285 (2022); S. Jalalzadeh, P. V. Moniz, “Challenging Routes In Quantum Cosmology ,” World Scientific (2022); K. P. Y. Thebault, Class. Quantum Grav. 40 055007 (2023).
  * [21] P. Jordan, Z. Phys., 157, 112 (1959);
C. Brans and R. H. Dicke, Phys. Rev., 124, 925 (1961); C. H. Brans, Phys.
Rev., 125, 2194 (1962).

  * [22] Y. Fujii and K.-Ichi Maeda, “The Scalar-Tensor Theory of Gravitation,” Cambridge University Press (2003);
V. Faraoni, “Cosmology in Scalar-Tensor Gravity,” Springer Science & Business
Media, (2004).

  * [23] P. G. O. Freund, Nucl. Phys. B 209, 146 (1982);
T. Appelquist, A. Chodos and P. G. O. Freund, “Modern Kaluza-Klein theories,”
(1987 AddisonWesley, Redwood City);
D. Z. Freedman and A. Van Proeyen, “Supergravity,” Cambridge University Press,
Cambridge, (2012);
E. S. Fradkin and A. A. Tseytlin, Phys. Lett. B 158, 316 (1985);
E. S. Fradkin and A. A. Tseytlin, Nucl. Phys. B 261, 1 (1985);
C. G. Callan , E. J. Martinec, M. J. Perry and D. Friedan, Nucl. Phys. B 262,
593 (1985);
C. G. Callan, I. R. Klebanov and M. J. Perry, Nucl. Phys. B 278, 78 (1986);
M. B. Green, J. H. Schwarz and E. Witten, “Superstring theory” (1987
Cambridge, Uk: Univ. Pr., Cambridge Monographs On Mathematical Physics);
D. B. Blaschke and M. P. Dabrowski, Entropy 2012, 14(10), 1978.

  * [24] A. Avilez and C. Skordis, Phys. Rev. Lett. 113, 011101 (2014);
C. Santos, R. Gregory, Annals of Physics 258, 111 (1997);
A. S. Majumdar, D. Gangopadhyay, L. P. Singh, MNRAS, 385, 1467 (2008);
A. Sheykhi, Phys. Rev. D 81, 023525 (2010);
J. Sola, A. G.-Valent, J. de Cruz Perez, C. M.-Pulido, Class. Quantum Grav.
37, 245003 (2020);
C. R. Almeida, O. Galkina, J. C. Fabris, Universe 2021, 7, 286.

  * [25] T. Matsuda and H. Nariai, Prog. Theor. Phys. 49, 1195 (1973).
  * [26] M. Shibata, K. Nakao and T. Nakamura Phys. Rev. D 50, 7304 (1994).
  * [27] A. H. Ziaie, K. Atazadeh, Y. Tavakoli, Class. Quant. Grav. 27 075016 (2010).
  * [28] P. Rudra, R. Biswas and U. Debnath, Astrophys. and Space Sci., 354 2101 (2014); Astrophys. and Space Sci., 339 135 (2012).
  * [29] A. Nakonieczna and D.-han Yeom, JHEP 05 (2016) 155; JHEP 02 (2016) 049.
  * [30] J. Hansen and D.-han Yeom, JCAP 09 (2015) 019.
  * [31] T. Chiba and J. Soda, Prog. Theor. Phys. 96 567 (1996).
  * [32] S.-W. Kim, Phys. Rev. D 34, 1011 (1986).
  * [33] S.-W. Kim, Phys. Lett A 124, 243 (1987).
  * [34] L. Yuanjie, L. Shijun and M. Weichan, Int. J. Theor. Phys. 32, 2173 (1993);
S.-W. Kim, IL Nuovo Cimento, 112, 363 (1997);
Y.-H. Wu and C.-H. Wang, Phys. Rev. D 86, 123519 (2012).

  * [35] S.-W. Kim and B. H. Cho, Phys. Rev. D 36, 2314 (1987).
  * [36] H. Shabani, A. H. Ziaie, Eur. Phys. J. C 79, 270 (2019).
  * [37] V. De Sabbata, M. Gasperini, “Introduction to Gravitation,” World Scientific, (1985);
V. De Sabbata and C. Sivaram, “Spin and Torsion in Gravitation,” World
Scientific, (1994).

  * [38] F. W. Hehl, P. von der Heyde, G. D. Kerlick and J. M. Nester, Rev. Mod. Phys., 48, 393 (1976).
  * [39] J. Weyssenhoff and A. Raabe, Acta Phys. Pol. 9, 7 (1947);
J. Weyssenhoff, in Max-Planck-Festschrift-1958, ed. by Kockel et al.
(Deutscher Verlap Wissensch., Berlin, 1958) p. 155;
F. Halbwachs, Theorie Relativiste des Fluides a Spin (GauthierVillars, Paris,
1960);
G. A. Maugin, Sur les fluides relativistes a spin. Ann. Inst. Henri Poincare
20, 41 (1974);
J. R. Ray, L. L. Smalley, Phys. Rev. D 27, 1383 (1983);
Y. N. Obukhov and V. A. Korotky, Class. Quantum Grav. 4, 1633 (1987);
K. Nomura, T. Shirafuji and K. Hayashi, Prog. Theor. Phys. 86, 1239 (1991);
L. L. Smalley and J. P. Krisch, Class. Quanrum Grav., 11, 2375 (1994).

  * [40] I. S. Nurgaliev, W. N. Ponomariev, Phys. Lett. B 130, 378 (1983); M. Gasperini, Phys. Rev. Lett., 56, 2873 (1986);
G. de Berredo-Peixoto, E. A. de Freitas, Class. Quantum Grav., 26, 175015
(2009).

  * [41] W. Israel, Nuovo Cimento B 44, 1 (1996); W. Israel, Nuovo Cimento B 48, 463 (1967); K. G. Suffern, J. Phys. A: Math. Gen. 15, 1599 (1982); C. Barrabes and G. F. Bressange, Class. Quantum Grav. 14, 805 (1997); F. Dahia and C. Romero, Phys. Rev. D 60 104019 (1999). J. D. Barrow and P. Parsons, Phys. Rev. D 55, 1906 (1997); L. Aviles, H. Maeda, C. Martinez, Class. Quantum Grav. 37 075022 (2020).
  * [42] N. O. Santos, Mon. Not. Roy. Astron. Soc. 216, 403 (1985); A. Wang and Y. Wu, Gen. Rel. Grav., 31 107 (1999); S. D. Maharaj, G. Govender and M. Govender, Gen. Relativ. Gravit. 44 1089 (2012).
  * [43] E. Poisson, “A Relativist’s Toolkit: The Mathematics of Black-Hole Mechanics”, Cambridge University Press (Cambridge, 2007).
  * [44] N. J. Poplawski, Phys. Rev. D 85, 107502 (2012); N. J. Poplawski, Astrophys. J. 832, 96 (2016); S. B. Medina, M. Nowakowski and D. Batic, Annals Phys. 400 64 (2019).